#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
/* ---------------------------------------------------------------------------------------------------
! The Monte-Carlo Simulation of the Stern-Gerlach Experiment (ANSI C Version)
!
! In quantum mechanics, the Stern.Gerlach experiment is a celebrated experiment (performed in 1922)
! 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 Potographic Plate |
! | Source |
! | |
! | | | S | || |
! | | |__________________| || |
! | \|/ | || |
! | --O--------> || |
! | /|\ | __________________ || |
! | | | | || |
! | T | | N | || |
! | . . . |
! | .<------- L ------>.<-- D -->. |
! | |
! +--------------------------------------------------------------+
!
! The program works on the command line, usage is:
!
! spin [-f012] [-q01] [default: -f1 and -q1]
!
! Options:
!
! Magnetic field gradient:
! -f0 sets up a uniform magnetic field along +z axis with B > 0
! -f1 sets up a non-uniform magnetic field along +z axis with dB/dz > 0 (default)
! -f2 sets up a non-uniform modulated by a single Gaussian
! magnetic field along +z axis with dB/dz > 0
!
! Quantum effect:
! -q0 Excludes the quantum effect. Sz is not quantized.
! That is cos(theta) is randomly oriented
! -q1 Includes the quantum effect. Sz is quantized. (default)
! That is cos(theta) = +- 1/sqrt(3)
!
! 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 : April 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
double Rnd(void);
double Maxwell_Boltzman(double, double);
double Field_Grad(int, double);
double sign(double);
void Command_Line(int, char **, int *, int *);
void usage(void);
int main(int argc, char **argv)
{
const int N = 10000; // Number of Silver, Ag(47), atoms
double T = 2000.0; // Temperature in Kelvin
double sz, muz, theta;
double s,sm,x0,z,z0,v,a,dBdz;
int i, field, quantum;
// Get the magnetic field type from command line
Command_Line(argc, argv, &field, &quantum);
if(!field && !quantum) printf("spin: dB/dz = 0, Sz is not quantized.\n");
if(!field && quantum) printf("spin: dB/dz = 0, Sz is quantized.\n");
if( field && !quantum) printf("spin: dB/dz > 0, Sz is not quantized.\n");
if( field && quantum) printf("spin: dB/dz > 0, Sz is quantized.\n");
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;
// print the final results to the user screen
printf("%12.4e %12.4e\n",x0,z);
}
} /* main */
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
double Maxwell_Boltzman(double m, double 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 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/4.0;
vmax = vp+vp/4.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 Field_Grad(int op, double x){
//------------------------------------------------------------
// Returns the components of the magnetic field
// gradient (in Tesla/meter) along the z direction.
//------------------------------------------------------------
double fg;
if(op==0) fg = 0.0; // uniform magnetic field in z axis
else if(op==1) fg = DBZ; // non-uniform in z but uniform in x direction
else
fg = DBZ*exp(-25*x*x/XMAX/XMAX);
// fg = DBZ*(1-4*x*x/XMAX/XMAX);
return fg;
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
double Rnd(void){
//--------------------------------------------------------------
// Returns an random number uniformly distributed between [0,1].
//--------------------------------------------------------------
double r;
r = (double) rand()/RAND_MAX;
return r;
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
double sign(double x){
//--------------------------------------------------------------
// Returns sign of x (-1 or +1)
//--------------------------------------------------------------
if(x>0.0) return +1.0;
else return -1.0;
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
void Command_Line(int argc, char **argv, int *f, int *q){
//------------------------------------------------------------
// Gets and checks command line arguments, and returs a
// field (f), and quantum (q).
//------------------------------------------------------------
if(argc == 1){
*f = 1;
*q = 1;
return;
}
if(argc > 3 || argc == 2) usage();
else{
if( strcmp(argv[1],"-f0")==0 ) *f = 0;
else if( strcmp(argv[1],"-f1")==0 ) *f = 1;
else if( strcmp(argv[1],"-f2")==0 ) *f = 2;
else usage();
if( strcmp(argv[2],"-q0")==0 ) *q = 0;
else if( strcmp(argv[2],"-q1")==0 ) *q = 1;
else usage();
}
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
void usage(void){
//--------------------------------------------------------------
// Prints the usage of the program.
//--------------------------------------------------------------
printf("spin: ANSI C version.\n\n");
printf("Usage: spin [-f012] [-q01]\n\n");
printf("Options:\n");
printf("Magnetic field gradient:\n");
printf(" -f0 sets up a uniform magnetic field along +z axis with B > 0\n");
printf(" -f1 sets up a non-uniform magnetic field along +z axis with dB/dz > 0 (default)\n");
printf(" -f2 sets up a non-uniform modulated by a single Gaussian\n");
printf(" magnetic field along +z axis with dB/dz > 0\n\n");
printf("Quantum effect:\n");
printf(" -q0 Excludes the quantum effect. Sz is not quantized.\n");
printf(" That is cos(theta) is randomly oriented\n");
printf(" -q1 Includes the quantum effect. Sz is quantized.\n");
printf(" That is cos(theta) = +- 1/sqrt(3)\n\n");
printf(" See for details: http://www1.gantep.edu.tr/~bingul/seminar/spin/\n\n");
exit(1);
}