Omc -> Om3 to match xmds code, Rabi frequencies match xmds

Whatch for factor of 2 in Rabi. NEED TO DOUBLE CHECK!

1 files changed, 35 insertions, 24 deletions

diff --git a/fortran/navy_four_levels/FourLevelPulseProp_Double.f95 b/fortran/navy_four_levels/FourLevelPulseProp_Double.f95
index 04786b2..39fdfe7 100644
--- a/fortran/navy_four_levels/FourLevelPulseProp_Double.f95
+++ b/fortran/navy_four_levels/FourLevelPulseProp_Double.f95
@@ -31,9 +31,6 @@ program FourLevelPulseProp_v3_Double
 	integer  zpts
 	parameter (zpts=10)	! Number of slices in z direction
 
-	! default values for peak amplitudes
-	real*8 Om1peak_def,Om2peak_def,Omcpeak_def  ! default values for field peak amplitudes
-	parameter (Om1peak_def=0.01,Om2peak_def=5.00,Omcpeak_def=0.0)
 	! ---------------------------------------------------------------------------
 
 	character*150 fname
@@ -47,7 +44,7 @@ program FourLevelPulseProp_v3_Double
 	parameter (NWrite=200)		!number of frames to actually write
 
 
-	real*8 Lcell,Om1peak,Om2peak,Omcpeak,pi,tmax,tmin,tp,tshift,t_end,t_start,t_elapsed
+	real*8 Lcell,Om1peak,Om2peak,Om3peak,pi,tmax,tmin,tp,tshift,t_end,t_start,t_elapsed
 
 
 	integer m,n
@@ -76,7 +73,7 @@ program FourLevelPulseProp_v3_Double
 	real*8 t,z(zptsmax)
 
 
-	complex*16 Om1(zptsmax),Om2(zptsmax),Omc(zptsmax),Om_vac(zptsmax)
+	complex*16 Om1(zptsmax),Om2(zptsmax),Om3(zptsmax),Om_vac(zptsmax)
 	complex*16 rho11(zptsmax),rho12(zptsmax),rho13(zptsmax),rho14(zptsmax),rho21(zptsmax),rho22(zptsmax),rho23(zptsmax),rho24(zptsmax)
 	complex*16 rho31(zptsmax),rho32(zptsmax),rho33(zptsmax),rho34(zptsmax),rho41(zptsmax),rho42(zptsmax),rho43(zptsmax),rho44(zptsmax)
 	complex*16 rho11_last(zptsmax),rho12_last(zptsmax),rho13_last(zptsmax),rho14_last(zptsmax)
@@ -92,8 +89,21 @@ program FourLevelPulseProp_v3_Double
 	parameter( hbar=1.054571726e-34 )	! reduced  Plank constant
 	parameter( epsil=8.85418781762e-12 )	! Permittivity of free space 
 
+
+	! --------------------------------------------------------------------------
+	! EM fields Rabi frequencies
+	! --------------------------------------------------------------------------
+	!
+	parameter(Gamma_super=2*pi*6e6)	! Characteristic decay of upper level in [1/s]
+	!
+	! default values for peak amplitudes
+	real*8 Om1peak_def,Om2peak_def,Om3peak_def  ! default values for field peak amplitudes
+	parameter( Om1peak_def=0.0025*(2*pi*1e6)/Gamma_super )
+	parameter( Om2peak_def=6*(2*pi*1e6)/Gamma_super )
+	parameter( Om3peak_def=0.0*(2*pi*1e6)/Gamma_super )
+	! --------------------------------------------------------------------------
+	! Atomic numbers (based on Rubidium 87)
 	! Atomic numbers (based on Rubidium 87)
-	Gamma_super=2*pi*6e6	! Characteristic decay of upper level in [1/s]
 	lambda=794.7e-9         ! Wavelength in [m]
 	Ndensity=1e10*(1e6)        ! Density in m^-1
 	! atomic decay parameters in units of Gamma_super
@@ -117,6 +127,7 @@ program FourLevelPulseProp_v3_Double
 	ga23=0.5*(W32+W31+W21)
 	ga24=0.5*(W21+W41+W42+W43)
 	ga34=0.5*(W31+W41+W32+W42+W43)
+
 	
 	! Atomic parameters
 	write (*,*)'New version from moved folder'
@@ -132,7 +143,7 @@ program FourLevelPulseProp_v3_Double
 	! Rabi frequency are unit less and scaled by Gamma_super i.e. Om1/Gamma_super -> Om1
 	Om1peak=Om1peak_def	! Field 1 peak scaled Rabi frequency for the pump at entrance of cell
 	Om2peak=Om2peak_def	! Field 2 peak scaled Rabi frequency for the pump at entrance of cell
-	Omcpeak=Omcpeak_def	! Field 3 peak scaled Rabi frequency for the pump at entrance of cell
+	Om3peak=Om3peak_def	! Field 3 peak scaled Rabi frequency for the pump at entrance of cell
 	! Maximum detuning in MHz for dispersion lineshape plot
 	delmax=0
 	Ga4=(W41+W42+W43)
@@ -196,7 +207,7 @@ program FourLevelPulseProp_v3_Double
 	do m=1,zpts
 		Om1(m)=cmplx(0.,0.)
 		Om2(m)=cmplx(0.,0.)
-		Omc(m)=cmplx(0.,0.)
+		Om3(m)=cmplx(0.,0.)
 
 		Om_vac(m)=cmplx(0.,0.)
 		rho11(m)=cmplx(1.,0.)	!Change this to change the initial condition
@@ -333,7 +344,7 @@ program FourLevelPulseProp_v3_Double
                 t=tmin+float(n-1)*dt
                 Om1(1)=Om1peak*exp(-(t-tshift)**2/(tp*tp))	
                 Om2(1)=Om2peak
-                Omc(1)=Omcpeak
+                Om3(1)=Om3peak
                 Om_vac(1)=Om1(1)
 
 		do m=1,zpts
@@ -364,33 +375,33 @@ program FourLevelPulseProp_v3_Double
 
 			Om2(m)=b1*Om2(m-1)+b2*rho32_last(m)+b3*rho32_last(m-1)
 			
-			Omc(m)=c1*Omc(m-1)+c2*rho42_last(m)+c3*rho42_last(m-1)
+			Om3(m)=c1*Om3(m-1)+c2*rho42_last(m)+c3*rho42_last(m-1)
 
 			rho11(m)=j1*rho11_last(m)+j2*rho22_last(m)+j3*rho33_last(m)+j4*conjg(Om1(m))*rho31_last(m)
 			rho11(m)=rho11(m)+j5*Om1(m)*rho13_last(m)+j6*conjg(Om1(m-1))*rho31_last(m)+j7*Om1(m-1)*rho13_last(m)
 			rho11(m)=rho11(m)+j8*rho44_last(m)
 			
 			rho12(m)=d1*rho12_last(m)+d2*conjg(Om1(m))*rho32_last(m)+d3*Om2(m)*rho13_last(m)
-			rho12(m)=rho12(m)+d4*Omc(m)*rho14_last(m)+d5*conjg(Om1(m-1))*rho32_last(m)
-			rho12(m)=rho12(m)+d6*Om2(m-1)*rho13_last(m)+d7*Omc(m-1)*rho14_last(m)
+			rho12(m)=rho12(m)+d4*Om3(m)*rho14_last(m)+d5*conjg(Om1(m-1))*rho32_last(m)
+			rho12(m)=rho12(m)+d6*Om2(m-1)*rho13_last(m)+d7*Om3(m-1)*rho14_last(m)
 
 			rho13(m)=e1*rho13_last(m)+e2*conjg(Om1(m))*(rho33_last(m)-rho11_last(m))
 			rho13(m)=rho13(m)+e3*conjg(Om2(m))*rho12_last(m)+e4*conjg(Om1(m-1))*(rho33_last(m)-rho11_last(m))
 			rho13(m)=rho13(m)+e5*conjg(Om2(m-1))*rho12_last(m)
 			
-			rho14(m)=f1*rho14_last(m)+f2*conjg(Om1(m))*rho34_last(m)+f3*conjg(Omc(m))*rho12_last(m)
-			rho14(m)=rho14(m)+f4*conjg(Om1(m-1))*rho34_last(m)+f5*conjg(Omc(m-1))*rho12_last(m)
+			rho14(m)=f1*rho14_last(m)+f2*conjg(Om1(m))*rho34_last(m)+f3*conjg(Om3(m))*rho12_last(m)
+			rho14(m)=rho14(m)+f4*conjg(Om1(m-1))*rho34_last(m)+f5*conjg(Om3(m-1))*rho12_last(m)
 			
 			rho21(m)=conjg(rho12(m))
 			
 			! rho22(m) needs to be calculated lower down
 
 			rho23(m)=gc1*rho23_last(m)+gc2*conjg(Om1(m))*rho21_last(m)+gc3*conjg(Om2(m))*(rho33_last(m)-rho22_last(m))
-			rho23(m)=rho23(m)+gc4*conjg(Omc(m))*rho43_last(m)+gc5*conjg(Om1(m-1))*rho21_last(m)
-			rho23(m)=rho23(m)+gc6*conjg(Om2(m-1))*(rho33_last(m)-rho22_last(m))+gc7*conjg(Omc(m-1))*rho43_last(m)
+			rho23(m)=rho23(m)+gc4*conjg(Om3(m))*rho43_last(m)+gc5*conjg(Om1(m-1))*rho21_last(m)
+			rho23(m)=rho23(m)+gc6*conjg(Om2(m-1))*(rho33_last(m)-rho22_last(m))+gc7*conjg(Om3(m-1))*rho43_last(m)
 			
-			rho24(m)=h1*rho24_last(m)+h2*conjg(Om2(m))*rho34_last(m)+h3*conjg(Omc(m))*(rho44_last(m)-rho22_last(m))
-			rho24(m)=rho24(m)+h4*conjg(Om2(m-1))*rho34_last(m)+h5*conjg(Omc(m-1))*(rho44_last(m)-rho22_last(m))
+			rho24(m)=h1*rho24_last(m)+h2*conjg(Om2(m))*rho34_last(m)+h3*conjg(Om3(m))*(rho44_last(m)-rho22_last(m))
+			rho24(m)=rho24(m)+h4*conjg(Om2(m-1))*rho34_last(m)+h5*conjg(Om3(m-1))*(rho44_last(m)-rho22_last(m))
 			
 			rho31(m)=conjg(rho13(m))
 			
@@ -400,8 +411,8 @@ program FourLevelPulseProp_v3_Double
 			rho33(m)=rho33(m)+k5*Om1(m)*rho13_last(m)+k6*Om2(m)*rho23_last(m)+k7*conjg(Om1(m-1))*rho31_last(m)
 			rho33(m)=rho33(m)+k8*conjg(Om2(m-1))*rho32_last(m)+k9*Om1(m-1)*rho13_last(m)+k10*Om2(m-1)*rho23_last(m)
 			
-			rho34(m)=i1*rho34_last(m)+i2*conjg(Omc(m))*rho32_last(m)+i3*Om1(m)*rho14_last(m)+i4*Om2(m)*rho24_last(m)
-			rho34(m)=rho34(m)+i5*conjg(Omc(m-1))*rho32_last(m)+i6*Om1(m-1)*rho14_last(m)+i7*Om2(m-1)*rho24_last(m)
+			rho34(m)=i1*rho34_last(m)+i2*conjg(Om3(m))*rho32_last(m)+i3*Om1(m)*rho14_last(m)+i4*Om2(m)*rho24_last(m)
+			rho34(m)=rho34(m)+i5*conjg(Om3(m-1))*rho32_last(m)+i6*Om1(m-1)*rho14_last(m)+i7*Om2(m-1)*rho24_last(m)
 			
 			rho41(m)=conjg(rho14(m))
 
@@ -409,8 +420,8 @@ program FourLevelPulseProp_v3_Double
 
 			rho43(m)=conjg(rho34(m))
 			
-			rho44(m)=l1*rho44_last(m)+l2*conjg(Omc(m))*rho42_last(m)+l3*Omc(m)*rho24_last(m)
-			rho44(m)=rho44(m)+l4*conjg(Omc(m-1))*rho42_last(m)+l5*Omc(m-1)*rho24(m)
+			rho44(m)=l1*rho44_last(m)+l2*conjg(Om3(m))*rho42_last(m)+l3*Om3(m)*rho24_last(m)
+			rho44(m)=rho44(m)+l4*conjg(Om3(m-1))*rho42_last(m)+l5*Om3(m-1)*rho24(m)
 			
 			rho22(m)=1-rho11(m)-rho33(m)-rho44(m)
 			
@@ -434,13 +445,13 @@ program FourLevelPulseProp_v3_Double
 
 
 			if ( iSkipped.ge.NSkip ) then ! Outputting every Nskip frame
-				write (9,120) t/Gamma_super,z(m)/Gamma_super*c,cdabs(Om1(m)),cdabs(Om2(m)),cdabs(Omc(m)),cdabs(Om_vac(m))
+				write (9,120) t/Gamma_super,z(m)/Gamma_super*c,cdabs(Om1(m)),cdabs(Om2(m)),cdabs(Om3(m)),cdabs(Om_vac(m))
 			end if
 
                 end do
 		! the very first z point is skipped in above cycle we compencete for it here
 		if ( iSkipped.ge.NSkip ) then ! Outputting every Nskip frame
-			write (9,120) t/Gamma_super,z(1)/Gamma_super*c,cdabs(Om1(1)),cdabs(Om2(1)),cdabs(Omc(1)),cdabs(Om_vac(1))
+			write (9,120) t/Gamma_super,z(1)/Gamma_super*c,cdabs(Om1(1)),cdabs(Om2(1)),cdabs(Om3(1)),cdabs(Om_vac(1))
 		end if
 
 		if (cdabs(Om2(zpts)).gt.cdabs(Omold))	 tpeak=t