nemmtomm

#include <stdio.h>
#include <stdlib.h>
#include <math.h>


#include <X11/Xlib.h>
#include <assert.h>
#include <unistd.h>


Display *dpy;
Window w;
GC gc;


#define num double

void def(num * v1,num x,num y,num z) {v1[0]=x;v1[1]=y;v1[2]=z;}

void cross(num * v1,num * v2,num * v3)
{
v1[0] = (v2[1]*v3[2])-(v3[1]*v2[2]);
v1[1] = (v2[2]*v3[0])-(v3[2]*v2[0]);
v1[2] = (v2[0]*v3[1])-(v3[0]*v2[1]);
};
num dot(num * v1,num *v2) { num t = v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];return t;}

void normalize(num *v1,num *v2)
{
num t = sqrt(dot(v2,v2));
v1[0] = v2[0]/t;
v1[1] = v2[1]/t;
v1[2] = v2[2]/t;
}




void pont(int x,int y,int szin)
{
XSetForeground(dpy,gc,szin);
XDrawPoint(dpy, w, gc, x,y);
}
void vonal(int x1,int y1,int x2,int y2,int szin)
{
XSetForeground(dpy,gc,szin);
XDrawLine(dpy, w, gc, x1,y1,x2,y2);
}
num sqr(num n)
{
return n*n;
}




num t,dt,szog,meresi_szog,m=1;
num pozicio[3], sebesseg[3], normal[3], I[3], B[3], F[3];
int stat[2],szamlalo=0,sz3=0,meres,B_valtas;


void set()
{
def(pozicio,100,0,0);
szog=M_PI/2*((float)(rand()%10000))/10000.0;
def(sebesseg,0.0,cos(szog)*10.0,sin(szog)*10.0);
def(B,0,0,1);

szamlalo=0;
B_valtas=100000+((rand()%1000)*100+(rand()%1000));
meres=B_valtas+50000+(rand()%1000)*50+(rand()%1000);
}
void stern_gerlach()
{
int a;

for(a=0;a<360;a+=10)
{
meresi_szog=a;

t=sin((meresi_szog/2.0)*M_PI*2/360.0);
t*=t;
printf("%f fok %f \n",meresi_szog,t);

pont(a,500-(int)(t*200.0),0xff0000);

dt=0.02;
set();


stat[0]=0;
stat[1]=0;
int szin=0;



while(1)
{
pozicio[0] += sebesseg[0]*dt;
pozicio[1] += sebesseg[1]*dt;
pozicio[2] += sebesseg[2]*dt;

normalize(normal,pozicio);//r=100
pozicio[0] = normal[0]*100.0;
pozicio[1] = normal[1]*100.0;
pozicio[2] = normal[2]*100.0;

num h=dot(sebesseg,normal);
sebesseg[0] -= normal[0]*h;
sebesseg[1] -= normal[1]*h;
sebesseg[2] -= normal[2]*h;

normalize(sebesseg,sebesseg);
sebesseg[0] *= 10.0;
sebesseg[1] *= 10.0;
sebesseg[2] *= 10.0;

cross(F,sebesseg , B );
h=dt*0.02/m;
sebesseg[0] += F[0]*h;
sebesseg[1] += F[1]*h;
sebesseg[2] += F[2]*h;


szamlalo++;
if(szamlalo==meres)//impulzusmomentum iranya a magneses terhez kepest
{
cross(I,pozicio,sebesseg);

t=dot(B,I);
if(t<0.0) stat[0]++;//ellentetes
else stat[1]++;//megegyezo

set();
}
else
if(szamlalo==B_valtas)//magneses ter elfordul
{
szog=meresi_szog*M_PI*2/360.0;
def(B,sin(szog),0.0,cos(szog));
}


#if 0
//beloves
pont(
500+(int)pozicio[1],
500+(int)pozicio[0],szin++);
cross(I,pozicio,sebesseg);
vonal(
500+(int)pozicio[1],
500+(int)pozicio[0],
500+(int)pozicio[1]+I[1]*0.1,
500+(int)pozicio[0]+I[0]*0.1,
0xff00ff );
#endif


sz3--;
if(sz3<0)//neha statisztika kiiras
{
sz3=1000000;

if(stat[1]) printf("%d/%d/%d %d %\n",stat[0],stat[1],stat[0]+stat[1],stat[0]*100/(stat[0]+stat[1]));
}
if((stat[0]+stat[1])>100) break;
}
num st=stat[0]*100/(stat[0]+stat[1]);
pont(a,500-(int)(st*2.0),0xffff00);
XFlush(dpy);
}



}


int main()
{

dpy = XOpenDisplay((0));
w = XCreateSimpleWindow(dpy, DefaultRootWindow(dpy), 0,0, 1000, 800, 0,0,0);

XSelectInput(dpy, w, StructureNotifyMask);
XMapWindow(dpy, w);

gc = XCreateGC(dpy, w, 0, (0));
XSetForeground(dpy,gc,0);

for(;;) { XEvent e; XNextEvent(dpy, &e); if (e.type == MapNotify)break; }


stern_gerlach();


XFlush(dpy);
getchar();

return 0;

}