#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
/* ---------------------------------------------------------------------------------------------------
! The Monte-Carlo Simulation of the Stern-Gerlach Experiment (ROOT Version)
!
! In quantum mechanics, the Stern.Gerlach experiment is a celebrated experiment (performed in 1921)
! on the deflection of particles, often used to illustrate basic principles of quantum mechanics.
! It can be used to demonstrate that electrons and atoms have intrinsically quantum properties, and
! how measurement in quantum mechanics affects the system being measured.
! This is a basic simulation of an "ideal" Stern-Gerlach experiment to spatially
! separate the electrons (spin-1/2) particles.
!
! Experimental set up
!
! +--------------------------------------------------------------+
! | |
! | Ag Slit Magnetic field Region Photographic Plate |
! | Source |
! | |
! | | | S | || |
! | | |__________________| || |
! | \|/ | || |
! | --O--------> || |
! | /|\ | __________________ || |
! | | | | || |
! | T | | N | || |
! | . . . |
! | .<------- L ------>.<-- D -->. |
! | |
! +--------------------------------------------------------------+
!
! The program can be executed on the root terminal as follows:
!
! .x spin.C(field, quantum) [default: field=1 and quantum=1]
!
! Examples:
!
! .x spin.C(0) field = 0 --> dB/dz = 0 (uniform magnetic field)
! .x spin.C(1) field = 1 --> dB/dZ > 0 (non-uniform magnetic field along z axis)
! .x spin.C(2) field = 2 --> dB/dZ > 0 (non-uniform magnetic field along z axis)
! .x spin.C(1,0) field = 1 & quantum=0 --> dB/dZ > 0 and do not consider quantum effect
!
! See for details: http://www1.gantep.edu.tr/~bingul/seminar/spin/
!
! All versions : http://www1.gantep.edu.tr/~bingul/seminar/spin/spin.f90 (Fortran 90)
! http://www1.gantep.edu.tr/~bingul/seminar/spin/spin.c (ANSI C)
! http://www1.gantep.edu.tr/~bingul/seminar/spin/spin.C (ROOT)
!
! Author : Dr. Ahmet Bingul <Ahmet.Bingul (at) gantep.edu.tr>, <Ahmet.Bingul (at) cern.ch>
!
! Date : Mar 2008
!
--------------------------------------------------------------------------------------------------- */
#define XMAX 0.05 // Slit width in x direction in m
#define ZMAX 0.005 // Slit width in z direction in m
#define kB 1.38E-23 // Bolzman constant in J/K
#define PI 3.14159265359 // Number PI
#define GS 2.0023 // Gyromagnetic ratio for the electron
#define EoM 1.759E+11 // e/m ratio for the electron in C/kg
#define HBAR 1.055E-34 // Plank constant in J.s
#define L 1.000E+0 // Length in m of the region containing magnetic filed along y axis
#define D 0.100E+0 // Distance in m between the end of magnetic filed and Photograhic plate
#define MAg 1.792E-25 // Mass of the silver atom kg
#define DBZ 100.0 // Field gradient, dB/dz, in T/m
void spin(Int_t field=1, Int_t quantum=1)
{
const Int_t N = 100000; // Number of Silver, Ag(47), atoms
Double_t T = 2000.0; // Temperature in Kelvin
Double_t X[N],Z[N];
Double_t sz, muz, phi, theta;
Double_t s,sm,x0,z,z0,v,a,dBdz;
Int_t i;
sm = HBAR*sqrt(3.0)/2.0; // Magnitude of the spin vector
srand(time(NULL)); // initilize the random number generator
// Loop over Silver atoms
for(i=0; i<N; i++)
{
// A random angle on the x-y plane in the range [0,2pi]
//phi = 2.0*PI*Rnd();
if(quantum){
s = +1.0;
if(Rnd() < 0.5) s = -1.0;
theta = acos(1.0/sqrt(3.0)); // angle and Sz
sz = sm*cos(theta)*s; // are quantized
}
else{
theta = acos(2.0*Rnd()-1.0); // angle and Sz
sz = sm*cos(theta); // are not quantized
}
// magnetic moment
muz = -sz*GS*EoM/2.0;
// Velocity of a Silver atom taken from a Maxwell-Boltzman distribution
v = Maxwell_Boltzman(MAg,T);
// Populate uniformly Silver atoms in x-z plane at y=0
x0 = -XMAX/2.0 + Rnd()*XMAX;
z0 = -ZMAX/2.0 + Rnd()*ZMAX;
// Get the field gradient
dBdz = Field_Grad(field,x0);
// Acceceleration and final position along z axis
a = muz*dBdz/MAg;
z = z0 + 0.5*a*(L/v)*(L/v) + sign(a)*D*sqrt(2.0*fabs(a)*L)/v;
X[i] = x0*100; // in cm
Z[i] = z*100; // in cm
// print the final results to the user screen
// printf("%12.4e %12.4e\n",x0,z);
}
TGraph *gr = new TGraph(N,X,Z);
TCanvas *c1 = new TCanvas("c1","Stern-Gerlach Experiment",10,10,400,400);
c1->SetFillColor(0);
//gr->SetMaximum(10);
//gr->SetMinimum(-10);
gr->GetXaxis()->SetTitle("X (cm)");
gr->GetYaxis()->SetTitle("Z (cm)");
gr->GetXaxis()->CenterTitle();
gr->GetYaxis()->CenterTitle();
gr->SetMarkerSize(0.15);
gr->SetMarkerStyle(21);
gr->SetMarkerColor(2);
gr->Draw("AP");
if(!field && !quantum) gr->SetTitle("SGE: dB/dz = 0, Sz is not quantized.");
if(!field && quantum) gr->SetTitle("SGE: dB/dz = 0, Sz is quantized");
if( field && !quantum) gr->SetTitle("SGE: dB/dz > 0, Sz is not quantized");
if( field && quantum) gr->SetTitle("SGE: dB/dz > 0, Sz is quantized");
} /* main */
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Double_t Maxwell_Boltzman(Double_t m, Double_t T){
//--------------------------------------------------------------
// Returns, for an atom of mass m (kg) and temperature T (K),
// a velocity in m/s which is randomly selected from
// a Maxwell-Boltzman Distribution function, Fmb, around the
// most probable velocity, vp. To do this,
// the Monte-Carlo Rejection method is used.
//--------------------------------------------------------------
Double_t r,v,vmax,vmin,vp,Ptest,fmb,fmax,kT,C;
kT = kB*T;
C = PI*sqrt(2.0) * pow(m/(PI*kT),1.5);
vp = sqrt(2.0*kT/m);
vmin = vp-vp/8.0;
vmax = vp+vp/8.0;
fmax = C * vp*vp * exp(-1.0);
// Rejection algorithm
while(1)
{
r = Rnd();
v = Rnd();
Ptest = fmax*r;
v = vmin + (vmax-vmin)*v;
fmb = C * v*v * exp(-0.5*m*v*v/kT);
if (fmb > Ptest) return v;
}
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Double_t Field_Grad(Int_t op, Double_t x){
//------------------------------------------------------------
// Returns the components of the magnetic field
// gradient (in Tesla/meter) along the z direction.
//------------------------------------------------------------
Double_t Field_Grad;
if(op==0) Field_Grad = 0.0; // uniform magnetic field in z axis
else if(op==1) Field_Grad = DBZ; // non-uniform in z but uniform in x direction
else
Field_Grad = DBZ*exp(-25*x*x/XMAX/XMAX);
// Field_Grad = DBZ*(1-4*x*x/XMAX/XMAX);
return Field_Grad;
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Double_t Rnd(void){
//--------------------------------------------------------------
// Returns an random number uniformly distributed between [0,1].
//--------------------------------------------------------------
return (Double_t) gRandom->Rndm();
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Double_t sign(Double_t x){
//--------------------------------------------------------------
// Returns sign of x (-1 or +1)
//--------------------------------------------------------------
if(x>0.0) return +1.0;
else return -1.0;
}