PROGRAM Stern_Gerlach
!---------------------------------------------------------------------------------------------------
! The Monte-Carlo Simulation of the Stern-Gerlach Experiment (Fortran 90 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
!
!---------------------------------------------------------------------------------------------------
IMPLICIT NONE
INTEGER, PARAMETER :: K = 8 ! Kind spcifier (K=4 for float, K=8 for double )
REAL(KIND=K), PARAMETER :: PI = 3.14159265359_K ! Number PI
INTEGER, PARAMETER :: N = 10000 ! Number of Silver, Ag(47), atoms
REAL(KIND=K), PARAMETER :: GS = 2.0023_K ! Gyromagnetic ratio for the electron
REAL(KIND=K), PARAMETER :: EOM = 1.759E+11_K ! e/m ratio for the electron in C/kg
REAL(KIND=K), PARAMETER :: HBAR= 1.055E-34_K ! Plank constant in J.s
REAL(KIND=K), PARAMETER :: L = 0.5E+0_K ! Length in m of the region containing magnetic filed along y axis
REAL(KIND=K), PARAMETER :: D = 0.5E+0_K ! Distance in m between the end of magnetic filed and Photograhic plate
REAL(KIND=K), PARAMETER :: MAg = 1.792E-25_K ! Mass of the silver atom kg
REAL(KIND=K), PARAMETER :: kB = 1.38E-23_K ! Bolzman constant in J/K
REAL(KIND=K), PARAMETER :: T = 2000.0_K ! Temperature in Kelvin
REAL(KIND=K), PARAMETER :: XMAX= 0.050_K ! Slit width in x direction in m
REAL(KIND=K), PARAMETER :: ZMAX= 0.005_K ! Slit width in z direction in m
REAL(KIND=K), PARAMETER :: DBZ = 100.0_K ! Field gradient, dB/dz, in T/m
REAL(KIND=K) :: sz, muz, theta
REAL(KIND=K) :: s,sm,x0,z,z0,v,a,dBdz
INTEGER :: i, field, quantum
CALL Command_Line(Field,Quantum) ! Get the magnetic field from command line
if( field==0 .AND. quantum==0) PRINT *,"spin: dB/dz = 0, Sz is not quantized."
if( field==0 .AND. quantum==1) PRINT *,"spin: dB/dz = 0, Sz is quantized."
if( field>=1 .AND. quantum==0) PRINT *,"spin: dB/dz > 0, Sz is not quantized."
if( field>=1 .AND. quantum==1) PRINT *,"spin: dB/dz > 0, Sz is quantized."
Sm = HBAR*SQRT(3.0_K)/2.0_K ! Magnitude of the spin vector
CALL RANDOM_SEED ! initilize the random number generator
!-- Loop over Silver atoms
DO I=1,N
! A random angle on the x-y plane in the range [0,2pi]
!phi = 2.0_K*PI*Rnd()
if(quantum==1) then
s = +1.0
if(Rnd() < 0.5) s = -1.0
theta = acos(1.0_K/sqrt(3.0_K)) ! angle and Sz
sz = sm*cos(theta)*s ! are quantized
else
theta = acos(2.0_K*Rnd()-1.0_K) ! angle and Sz
sz = sm*cos(theta); ! are not quantized
end if
! magnetic moment
muz = -sz*GS*EoM/2.0_K
! 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_K + Rnd()*XMAX;
z0 = -ZMAX/2.0_K + 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*abs(a)*L)/v;
! print the final results to the user screen
PRINT '(4(ES12.4,3x))',x0,z
END DO
CONTAINS
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
REAL(KIND=K) FUNCTION Maxwell_Boltzman(m,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.
!--------------------------------------------------------------
REAL(KIND=K), INTENT(IN) :: m,T
REAL(KIND=K) :: r,v,vmax,vmin,vp,Ptest,fmb,fmax,kT,C
kT = kB*T
C = PI*SQRT(2.0_K)*(m/PI/kT)**1.5_K
vp = SQRT(2.0_K*kT/m)
vmin = vp-vp/4.0_K
vmax = vp+vp/4.0_K
fmax = C * vp**2 * EXP(-1.0_K)
! Rejection algorithm
DO
r = Rnd();
v = Rnd();
Ptest = fmax*r
v = vmin + (vmax-vmin)*v
Fmb = C * v**2 * exp(-0.5_K*m*v**2/kT)
IF (Fmb > Ptest) THEN
Maxwell_Boltzman = v
RETURN
END IF
END DO
END FUNCTION Maxwell_Boltzman
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
REAL(KIND=K) FUNCTION Field_Grad(op,x)
!------------------------------------------------------------
! Returns the components of the magnetic field
! gradient (in Tesla/meter) along the z direction.
!------------------------------------------------------------
INTEGER, INTENT(IN) :: op
REAL(KIND=K), INTENT(IN) :: x
IF(op==0) THEN ! uniform magnetic field in z axis
Field_Grad = 0.0_K
ELSE ! non-uniform in z but uniform in x direction
Field_Grad = DBZ
END IF
END FUNCTION Field_Grad
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
REAL(KIND=K) FUNCTION Rnd()
!--------------------------------------------------------------
! Returns an random number uniformly distributed between [0,1].
!--------------------------------------------------------------
REAL(KIND=K) :: R
CALL RANDOM_NUMBER(R)
Rnd = R
END FUNCTION Rnd
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
REAL(KIND=K) FUNCTION sign(x)
!--------------------------------------------------------------
! Returns sign of x (-1 or +1)
!--------------------------------------------------------------
REAL(KIND=K), INTENT(IN) :: x
if(x>0.0) THEN
sign = +1.0_K;
else
sign = -1.0_K;
END IF
END FUNCTION sign
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SUBROUTINE Command_Line(f,q)
!------------------------------------------------------------
! Gets and checks command line arguments, and returs a
! field (f), and quantum (q).
!------------------------------------------------------------
INTEGER, INTENT(OUT) :: f,q
CHARACTER(LEN=5) :: Argv1, Argv2
INTEGER :: Argc
Argc = COMMAND_ARGUMENT_COUNT()
IF( Argc==0 ) THEN
f = 1
q = 1
RETURN
END IF
IF(Argc > 2 .OR. Argc == 1) THEN
CALL Usage()
ELSE
CALL GET_COMMAND_ARGUMENT(1,Argv1)
CALL GET_COMMAND_ARGUMENT(2,Argv2)
IF( TRIM(argv1) == "-f0" ) THEN
f = 0;
ELSE IF( TRIM(argv1) == "-f1") THEN
f = 1;
ELSE IF( TRIM(argv1) == "-f2") THEN
f = 2;
ELSE
CALL Usage();
END IF
IF( TRIM(argv2) == "-q0" ) THEN
q = 0;
ELSE IF( TRIM(argv2) == "-q1") THEN
q = 1;
ELSE
CALL Usage();
END IF
END IF
END SUBROUTINE Command_Line
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SUBROUTINE usage()
!--------------------------------------------------------------
! Prints the usage of the program.
!--------------------------------------------------------------
PRINT *,"spin: Fortran 90 version."
PRINT *
PRINT *,"Usage: spin [-f012] [-q01]"
PRINT *
PRINT *,"Options:"
PRINT *,"Magnetic field gradient:"
PRINT *," -f0 sets up a uniform magnetic field along +z axis with B > 0"
PRINT *," -f1 sets up a non-uniform magnetic field along +z axis with dB/dz > 0 (default)"
PRINT *," -f2 sets up a non-uniform modulated by a single Gaussian"
PRINT *," magnetic field along +z axis with dB/dz > 0"
PRINT *
PRINT *,"Quantum effect:"
PRINT *," -q0 Excludes the quantum effect. Sz is not quantized."
PRINT *," That is cos(theta) is randomly oriented"
PRINT *," -q1 Includes the quantum effect. Sz is quantized."
PRINT *," That is cos(theta) = +- 1/sqrt(3)"
PRINT *
PRINT *," See for details: http://www1.gantep.edu.tr/~bingul/seminar/spin/"
PRINT *
STOP
END SUBROUTINE usage
END PROGRAM Stern_Gerlach