program FourLevelPulseProp_v3_Double ! ! Written by: Dr. Frank A. Narducci ! Written on: May 12, 2008 ! ! This program calculates the propagation of a pulse of arbitrary strength ! through a two level medium. The equations used are the full equations ! based on the Risken-Numedal discretization technique. ! ! This program only "watches" the evolution of the pulse in the cell. This is ! due to the constraint that the cell is very small relative to the pulse lengths ! that we want to use. If we watched the pulse outside the cell and then increased the ! resolution within the cell, the increased burden outside the cell because huge. ! ! v2 Notes: This program is based on the dimensionless equatiosn derived on 5/16 implicit none ! ! Double Notes: This program is the same as TwoLevelPulseProp_v2 but with double precision complex ci ! ! ThreeLevel Notes: This program is the same as TwoLevelPulseProp_v2_Double but now for a ! three level system ! ! Four Level Notes: Valid to first order in dt ! ! Version 3: Make a step function in the coupling field. ! --------------------------------------------------------------------------- ! SYSTEM DEFAULTS ! --------------------------------------------------------------------------- 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 integer nmat, npts, Nframe, tpts_max, NSkip, NWrite, zptsmax integer iSkipped; parameter (nmat=3,npts=100) !matrix size, number of detuning points in dispersion curve !tpts is the number of temporal points in the cell !zpts is the number of z output slices - 1 parameter (zptsmax=200) ! max Number of slices in z direction parameter (tpts_max=200) ! max Number of slices in t direction 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 integer m,n complex*16 a1,a2,a3,a4,a5,a6 complex*16 b1,b2,b3,b4,b5,b6,b7 complex*16 c1,c2,c3,c4,c5,c6 complex*16 d1,d2,d3,d4,d5,d6,d7 complex*16 e1,e2,e3,e4,e5 complex*16 f1,f2,f3,f4,f5 complex*16 g1,g2,g3,g4,g5,g6,g7 complex*16 h1,h2,h3,h4,h5 complex*16 i1,i2,i3,i4,i5,i6,i7 complex*16 j1,j2,j3,j4,j5,j6,j7,j8 complex*16 k1,k2,k3,k4,k5,k6,k7,k8,k9,k10 complex*16 l1,l2,l3,l4,l5 complex*16 Omold, Omold_vac real*8 alpha1,alpha2,alpha1tilde,alpha2tilde,alphac,alphactilde,Gamma_super,c,delmax,del1_prop,del2_prop,delc_prop real*8 dt,dz,Ndensity real*8 W12,W21,W31,W32,W41,W42,W43,W34,ga12,ga13,ga14,ga23,ga24,ga34 real*8 Ga2,Ga4,Om_crit real*8 tpeak,tpeak_vac real*8 epsil,hbar,lambda real*8 t,z(zptsmax) complex*16 Om1(zptsmax),Om2(zptsmax),Omc(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) complex*16 rho21_last(zptsmax),rho22_last(zptsmax),rho23_last(zptsmax),rho24_last(zptsmax) complex*16 rho31_last(zptsmax),rho32_last(zptsmax),rho33_last(zptsmax),rho34_last(zptsmax) complex*16 rho41_last(zptsmax),rho42_last(zptsmax),rho43_last(zptsmax),rho44_last(zptsmax) ! Fundamental numbers parameter( ci=cmplx(0.,1.) ) ! imaginary one parameter( pi=acos(-1.0) ) ! pi=3.1415.... parameter( c=3e8 ) ! speed of light parameter( hbar=1.054571726e-34 ) ! reduced Plank constant parameter( epsil=8.85418781762e-12 ) ! Permittivity of free space ! 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 W41=0 W42=1 W43=0 W32=1. W31=1. W21=.001 W12=W21 W34=0 ga12=0.5*(W21+W12) ga13=0.5*(W31+W12+W32) ga14=0.5*(W41+W42+W12) 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' alpha1=3*Ndensity*lambda*lambda/(2*pi) alpha1tilde=alpha1*c/Gamma_super alpha2=3*Ndensity*lambda*lambda/(2*pi) alpha2tilde=alpha2*c/Gamma_super alphac=3*Ndensity*lambda*lambda/(2*pi) alphactilde=alphac*c/Gamma_super ! User defined numbers ! 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 ! Maximum detuning in MHz for dispersion lineshape plot delmax=0 Ga4=(W41+W42+W43) Ga2=W21 Om_crit=sqrt(Om1peak**4+4*Om1peak*Om1peak*(Ga4*Ga4+Ga2*Ga4)) Om_crit=Om_crit-Om1peak*Om1peak-2*Ga2*Ga4 Om_crit=sqrt(Om_crit/2) write (*,*)'Om_crit = ',Om_crit ! Now that the user has an idea of the dispersion, do the full propagation problem ! Detuning of center frequency of the coupling pulse in MHz delc_prop=0 delc_prop=2*pi*1e6*delc_prop/Gamma_super !Now dimensionless ! Detuning of center frequency of the pump pulse in MHz del2_prop=0. del2_prop=2*pi*1e6*del2_prop/Gamma_super !Now dimensionless ! Detuning of center frequency of the probe pulse in MHz del1_prop=0. del1_prop=2*pi*1e6*del1_prop/Gamma_super !Now dimensionless ! Pulse width in nsec tp=1e-6 tmax=tp*5 tmin = -tmax tshift=0 tp=Gamma_super*tp !Now dimensionless tmax=Gamma_super*tmax tmin=Gamma_super*tmin tshift=Gamma_super*tshift ! Length of cell in m Lcell=0.07 Lcell=Gamma_super*Lcell/c !Now dimensionless t_start=secnds(0.E0) write (*,*)'t_start = ',t_start ! ! Set up initial pulse. ! write (*,*) 'peak center at the cell begining i.e. tshift = ', tshift write(*,*) 'tmax = ', tmax write(*,*) 'Lcell = ', Lcell dz=Lcell/(zpts-1) !(no c because we're dimensionless) ! It is crucial that dz = dt in unitless coordinates, there is built in ! cancellation of some term on grid because of it. See Simon's note. dt = dz Nframe=(tmax-tmin)/dz+1 if (Nframe.ge.tpts_max) write (*,*)'Error!!!!Nframe>tpts_max' ! ! Initialize matrices ! Omold=cmplx(0.,0.) Omold_vac=cmplx(0.,0.) tpeak=-1 tpeak_vac=-1 write (*,*)'Nframe= ', Nframe do m=1,zpts Om1(m)=cmplx(0.,0.) Om2(m)=cmplx(0.,0.) Omc(m)=cmplx(0.,0.) Om_vac(m)=cmplx(0.,0.) rho11(m)=cmplx(1.,0.) !Change this to change the initial condition rho12(m)=cmplx(0.,0.) rho13(m)=cmplx(0.,0.) rho14(m)=cmplx(0.,0.) rho21(m)=cmplx(0.,0.) rho22(m)=cmplx(0.,0.) rho23(m)=cmplx(0.,0.) rho24(m)=cmplx(0.,0.) rho31(m)=cmplx(0.,0.) rho32(m)=cmplx(0.,0.) rho33(m)=cmplx(0.,0.) rho34(m)=cmplx(0.,0.) rho41(m)=cmplx(0.,0.) rho42(m)=cmplx(0.,0.) rho43(m)=cmplx(0.,0.) rho44(m)=cmplx(0.,0.) end do ! Propagation co-efficients ! a1=1. a2=0.5*ci*alpha1tilde*dt a3=0.5*ci*alpha1tilde*dt a4=0. a5=0. a6=0. b1=1. b2=0.5*ci*alpha2tilde*dt b3=0.5*ci*alpha2tilde*dt b4=0. b5=0. b6=0. b7=0. c1=1. c2=0.5*ci*alphactilde*dt c3=0.5*ci*alphactilde*dt c4=0. c5=0. c6=0. d1=1-(ga12-ci*(del2_prop-del1_prop))*dt d2=0.25*ci*dt d3=-0.25*ci*dt d4=-0.25*ci*dt d5=0.25*ci*dt d6=-0.25*ci*dt d7=-0.25*ci*dt e1=1-(ga13+ci*del1_prop)*dt e2=0.25*ci*dt e3=-0.25*ci*dt e4=0.25*ci*dt e5=-0.25*ci*dt f1=1-(ga14-ci*(del2_prop-del1_prop-delc_prop))*dt f2=0.25*ci*dt f3=-0.25*ci*dt f4=0.25*ci*dt f5=-0.25*ci*dt g1=1-(ga23+ci*del1_prop)*dt g2=-0.25*ci*dt g3=0.25*ci*dt g4=0.25*ci*dt g5=-0.25*ci*dt g6=0.25*ci*dt g7=0.25*ci*dt h1=1-(ga24+ci*delc_prop)*dt h2=0.25*ci*dt h3=0.25*ci*dt h4=0.25*ci*dt h5=0.25*ci*dt i1=1-(ga34-ci*(del2_prop-delc_prop))*dt i2=-0.25*ci*dt i3=0.25*ci*dt i4=0.25*ci*dt i5=-0.25*ci*dt i6=0.25*ci*dt i7=0.25*ci*dt j1=1-W12*dt j2=W21*dt j3=W31*dt j4=0.25*ci*dt j5=-0.25*ci*dt j6=0.25*ci*dt j7=-0.25*ci*dt j8=W41*dt k1=1-(W32+W31+W34)*dt k2=W43*dt k3=-0.25*ci*dt k4=-0.25*ci*dt k5=0.25*ci*dt k6=0.25*ci*dt k7=-0.25*ci*dt k8=-0.25*ci*dt k9=0.25*ci*dt k10=0.25*ci*dt l1=1-(W43+W42+W41)*dt l2=-0.25*ci*dt l3=0.25*ci*dt l4=-0.25*ci*dt l5=0.25*ci*dt NSkip=int(NFrame/NWrite) fname='MovieParameters4level.txt' ! File name to save parameters open(9,file=fname) write (9,133)Nframe,zpts,Gamma_super,NSkip,dt 133 format(1x,i10,',',i5,',',f12.2,',',i5,',',f12.2) close (9) fname='Movie4level.dat' ! File name to save movie open(9,file=fname) fname='Movie4level_EndPoints.dat' ! File name to save endpoints open(10,file=fname) iSkipped=NSkip; ! since we want initial time points to be at output as well do n=1,Nframe iSkipped = iSkipped+1 t=tmin+float(n-1)*dt Om1(1)=Om1peak*exp(-(t-tshift)**2/(tp*tp)) Om2(1)=Om2peak Omc(1)=Omcpeak Om_vac(1)=Om1(1) do m=1,zpts rho11_last(m)=rho11(m) rho12_last(m)=rho12(m) rho13_last(m)=rho13(m) rho14_last(m)=rho14(m) rho21_last(m)=rho21(m) rho22_last(m)=rho22(m) rho23_last(m)=rho23(m) rho24_last(m)=rho24(m) rho31_last(m)=rho31(m) rho32_last(m)=rho32(m) rho33_last(m)=rho33(m) rho34_last(m)=rho34(m) rho41_last(m)=rho41(m) rho42_last(m)=rho42(m) rho43_last(m)=rho43(m) rho44_last(m)=rho44(m) end do do m=zpts,2,-1 z(m)=float(m-1)*dz Om_vac(m)=a1*Om_vac(m-1) Om1(m)=a1*Om1(m-1)+a2*rho31_last(m)+a3*rho31_last(m-1) 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) 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) 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) rho21(m)=conjg(rho12(m)) ! rho22(m) needs to be calculated lower down rho23(m)=g1*rho23_last(m)+g2*conjg(Om1(m))*rho21_last(m)+g3*conjg(Om2(m))*(rho33_last(m)-rho22_last(m)) rho23(m)=rho23(m)+g4*conjg(Omc(m))*rho43_last(m)+g5*conjg(Om1(m-1))*rho21_last(m) rho23(m)=rho23(m)+g6*conjg(Om2(m-1))*(rho33_last(m)-rho22_last(m))+g7*conjg(Omc(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)) rho31(m)=conjg(rho13(m)) rho32(m)=conjg(rho23(m)) rho33(m)=k1*rho33_last(m)+k2*rho44_last(m)+k3*conjg(Om1(m))*rho31_last(m)+k4*conjg(Om2(m))*rho32_last(m) 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) rho41(m)=conjg(rho14(m)) rho42(m)=conjg(rho24(m)) 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) rho22(m)=1-rho11(m)-rho33(m)-rho44(m) rho11_last(m)=rho11(m) rho12_last(m)=rho12(m) rho13_last(m)=rho13(m) rho14_last(m)=rho14(m) rho21_last(m)=rho21(m) rho22_last(m)=rho22(m) rho23_last(m)=rho23(m) rho24_last(m)=rho24(m) rho31_last(m)=rho31(m) rho32_last(m)=rho32(m) rho33_last(m)=rho33(m) rho34_last(m)=rho34(m) rho41_last(m)=rho41(m) rho42_last(m)=rho42(m) rho43_last(m)=rho43(m) rho44_last(m)=rho44(m) 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)) 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)) end if if (cdabs(Om2(zpts)).gt.cdabs(Omold)) tpeak=t if (cdabs(Om_vac(zpts)).gt.cdabs(Omold_vac)) tpeak_vac=t if ( iSkipped.ge.NSkip ) then ! Outputting every Nskip frame ! Outputting Rabi frequency at the beginning and the end of ! the cell for given time to the EndPoint File write (10,139) t/Gamma_super, cdabs(Om1(1)), cdabs(Om_vac(1)), cdabs(Om1(zpts)), cdabs(Om_vac(zpts)) iSkipped = 0 end if Omold=Om2(zpts) Omold_vac=Om_vac(zpts) end do close(9) close(10) 139 format(1x,E15.9,',',E15.9,',',E15.9,',',E15.9,',',E15.9) 120 format(E15.9,',',E15.9,',',E15.9,',',E15.9,',',E15.9,',',E15.9,',',E15.9,',',E15.9,',',E15.9) t_end=secnds(0.E0) t_elapsed=t_end-t_start write(*,*)'T elapsed = ',t_elapsed stop end