realspin

float1 t,dt,angle,m=1,measure_angle;
float3 pos,spd,nor,I,B,F;
int stat[2],counter=0,c3=0,measure,B_change;


void init()
{
        pos=float3(100,0,0);
        angle=M_PI/2*((float)(rand()%10000))/10000.0;
        spd=float3(0.0,cos(angle)*10.0,sin(angle)*10.0);
        B=float3(0,0,1);

        counter=0;
        B_change=200000+((rand()%100)*1000+(rand()%1000));
        measure=B_change+200000+(rand()%100)*1000+(rand()%1000);
}
void stern_gerl()
{
    //ennyi szazalekban all at a spin down-rol up-ra a QM szerint, ha a ket magneses ter a fokban all egymashoz kepest
    measure_angle=60.0;//45=14.6% 60=25% 90=50% ennyi szazalekban all at a spin down-rol up-ra
   
    t=sin((measure_angle/2.0)*M_PI*2/360.0);
    t*=t;
    printf("%Lf fok %Lf \n",measure_angle,t);

    dt=0.01;
    init();


    stat[0]=0;
    stat[1]=0;



    while(1)
    {
        pos=pos+spd*dt;
        nor=normalize(pos);
        pos=nor*100.0;//egy gomb feluleten mozog a pont
        spd=spd-nor*dot(spd,nor);
        spd=normalize(spd)*10.0;
   
        F=cross(spd,B);//lorentz ero  *e most nincs, egyszerusitett modell
        spd=spd+F*dt*0.01/m;


        counter++;
        if(counter==measure)//impulzusmomentum iranya a magneses terhez kepest
        {
            I=cross(pos,spd);

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

            init();
        }
        else
        if(counter==B_change)//magneses ter elfordul
        {
            angle=measure_angle*M_PI*2/360.0;
            B=float3(sin(angle),0.0,cos(angle));
        }

        c3--;
        if(c3<0)//neha statisztika kiiras
        {
            c3=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]));
        }
    }
   

}

    int n=0;
    float1 h=6.626e-34;
    float3 Bv=float3(0,0,0);
    for(int y=0;y<900;y++)
    for(int x=0;x<3600;x++)
    {
        float1 a1=(float1)x*M_PI*2/3600.0;
        float1 a2=(float1)y*M_PI*2/3600.0;
        Bv=Bv+float3(cos(a1)*cos(a2),sin(a1)*cos(a2),sin(a2))*h;
        n++;
    }
    Bv=Bv/n;
    printf("%Lf  \n",atanl(h/Bv.z)*180.0/M_PI);
    printf("%Lf  \n",((h/Bv.z)/(M_PI/2)));
 

http://en.wikipedia.org/wiki/Bohr_magneton

int main()
{
    long double
        e=1.6e-19,
        h=6.626e-34,
        m=9.1e-31,
        c=3e8,
        alf=1.0/137.0,
        r,hbar,v,l,dT,T,B,b,v2,Bb,lc,f,amplitudec,xe=0,ye=0,time,phase,amplitudex,amplitudey;
    int i,n=1000000;
    float3 rv,vv,Bv;

    Bv=float3(0,0,0);

    hbar=h/(M_PI*2);
    r=hbar/(m*c*2);
    Bb=r*c*e;
    printf("%.11Le\n",Bb);


    Bb=e*hbar/(m*2);
    printf("%.11Le\n",Bb);





    lc=h/(m*c);
    amplitudec=lc/(M_PI);


n=1000000;
    v=c*alf;
    b=1.0/sqrtl(1.0-v*v/(c*c));
    l=h/(m*v*b);
   

    f=v/l;
    T=1.0/f;
    dT=T/n;
    B=0;

   

    for(i=-1;i<n;i++)
    {
        time=dT*(long double)i-T/2;

        phase=(time*c-time*v)*M_PI/lc;
        amplitudex=cos(phase)*0.5*amplitudec;
        amplitudey=sin(phase)*0.5*amplitudec;
        phase=(time*c+time*v)*M_PI/lc;
        amplitudex+=cos(phase)*0.5*amplitudec;
        amplitudey+=sin(phase)*0.5*amplitudec;

        if(i>=0)
        {
            rv=float3(amplitudex,amplitudey,(float1)0);
            vv=float3(amplitudex-xe,amplitudey-ye,dT*v)/dT;

            r=sqrtl(sqr(amplitudex)+sqr(amplitudey));
            v2=sqrtl(sqr(amplitudex-xe)+sqr(amplitudey-ye))/dT;

           
            float3 q=cross(vv,rv)*e;
            //if(q.x>0.0)// x komponens meres   
                Bv=Bv+q;

//            if(i%100000==0)    {        printf("%Le \n",c/v2);}
            B+=e*v2*r;
        }
        xe=amplitudex;
        ye=amplitudey;
    }

    B/=n;
    Bv=Bv/n;
    printf("%.11Le\n",B/2);
    printf("%.11Le\n",B/Bb);


    B=length(Bv);

    printf("%.11Le\n",B/2);
    printf("%.11Le\n",B/Bb);


    return 0;
}