summaryrefslogtreecommitdiff
path: root/code_from_navy/fortran/FourLevelPrograms/FourLevelPulseProp_v2_Double.f
diff options
context:
space:
mode:
Diffstat (limited to 'code_from_navy/fortran/FourLevelPrograms/FourLevelPulseProp_v2_Double.f')
-rw-r--r--code_from_navy/fortran/FourLevelPrograms/FourLevelPulseProp_v2_Double.f1416
1 files changed, 708 insertions, 708 deletions
diff --git a/code_from_navy/fortran/FourLevelPrograms/FourLevelPulseProp_v2_Double.f b/code_from_navy/fortran/FourLevelPrograms/FourLevelPulseProp_v2_Double.f
index 3ec1e70..9ab27bf 100644
--- a/code_from_navy/fortran/FourLevelPrograms/FourLevelPulseProp_v2_Double.f
+++ b/code_from_navy/fortran/FourLevelPrograms/FourLevelPulseProp_v2_Double.f
@@ -1,708 +1,708 @@
- 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.
-
- character*150 fname
- integer nmat,npts,Nfrac,Nframe,Nframemax,NSkip,NWrite,tpts,zpts
- parameter (nmat=3,npts=100) !matrix size, number of detuning points in dispersion curve
- !REMEMBER TO CHANGE NMAT IN LMatConstruct Routine
- parameter (tpts=10,zpts=tpts+1) !Caution: funny things happened when tpts=200 (and presumably greater)
- !tpts is the number of temporal points in the cell
- parameter (Nframemax=2000000)
- parameter ( NWrite=100) !number of frames to actually write
- integer i,j,k,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,beta,c,delmax,del1_prop,del2_prop,delc_prop
- real*8 dt,dz,eta
- real*8 W12,W21,W31,W32,W41,W42,W43,W34,ga12,ga13,ga14,ga23,ga24,ga34
- real*8 Ga2,Ga4,Om_crit
- real*8 Lcell,Om1peak,Om2peak,Omcpeak,pi,tmax,tp,tshift,t_end,t_start,t_elapsed
- real*8 tpeak,tpeak_vac
- real*8 epsil,hbar,lambda
- real*8 del(npts)
- real*8 t,z(zpts)
- complex*16 yplot(nmat,npts)
- complex*16 Imat(nmat)
-
-
- complex*16 Om1(zpts),Om2(zpts),Omc(zpts),Om_vac(zpts)
- complex*16 rho11(zpts),rho12(zpts),rho13(zpts),rho14(zpts),rho21(zpts),rho22(zpts),rho23(zpts),rho24(zpts)
- complex*16 rho31(zpts),rho32(zpts),rho33(zpts),rho34(zpts),rho41(zpts),rho42(zpts),rho43(zpts),rho44(zpts)
- complex*16 rho11_last(zpts),rho12_last(zpts),rho13_last(zpts),rho14_last(zpts)
- complex*16 rho21_last(zpts),rho22_last(zpts),rho23_last(zpts),rho24_last(zpts)
- complex*16 rho31_last(zpts),rho32_last(zpts),rho33_last(zpts),rho34_last(zpts)
- complex*16 rho41_last(zpts),rho42_last(zpts),rho43_last(zpts),rho44_last(zpts)
-
- !No Om_last because we never need the previous spatial point
- complex*16 L(nmat,nmat),Linv(nmat,nmat),Ltemp(nmat,nmat)
-
- common/para/ga12,W21
-
-
- real*8 d !used by NR Routines
- integer indx(nmat)
-
-!
-! Fundamental numbers
-!
- ci=cmplx(0.,1.)
- pi=acos(-1.0)
- c=3e8
- hbar=1.055e-34
- epsil=8.85e-12
-
-!
-! Atomic numbers (based on Rubidium 85)
-!
- beta=2*pi*3e6 !in Hz
- W41=0
- W42=0.
- 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)
- lambda=780.24e-9
-
-!
-! Atomic parameters
-!
- write (*,*)'New version from moved folder'
-! write (*,*)'Enter density in m^-1'
-! read (*,*)eta
- eta=6.9e11
- alpha1=3*eta*lambda*lambda/(2*pi)
- alpha1tilde=alpha1*c/beta
- alpha2=3*eta*lambda*lambda/(2*pi)
- alpha2tilde=alpha2*c/beta
- alphac=3*eta*lambda*lambda/(2*pi)
- alphactilde=alphac*c/beta
-
-!
-! Initialize matrices and set up Identity Matrix
-!
- do 20 i=1,nmat
- do 10 j=1,nmat
- L(i,j)=cmplx(0.,0.)
- Linv(i,j)=cmplx(0.,0.)
-10 continue !j loop
- Imat(i)=cmplx(0.,0.)
- Linv(i,i)=cmplx(1.,0.) !contains identity matrix
-20 continue
-
-
-!
-! User defined numbers
-!
-! write (*,*)'Enter peak scaled Rabi frequency for the pump at entrance of cell'
-! read (*,*)Om1peak
- Om1peak=1
-! write (*,*)'Enter peak scaled Rabi frequency for the probe at entrance of cell'
-! read (*,*)Om2peak
- Om2peak=1e-2
-! write (*,*)'Enter maximum detuning in MHz for dispersion lineshape plot'
-! read (*,*) delmax
- 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
-
-! write (*,*)'Enter peak scaled Rabi frequency for the coupling field at entrance of cell'
-! read (*,*)Omcpeak
- Omcpeak=0.
-
-!
-! First plot the dispersion lineshape
-!
-! do 40 n=1,npts
-! del(n)=-delmax+2*float(n)*delmax/npts
-! call LMatConstruct(Ompeak,del(n),L) !construct the L matrix
-! call LMatConstruct(Ompeak,del(n),Ltemp) !Need a temporary because L gets destroyed
-! call ImatConstruct(Ompeak,Imat) !Need to call in loop because it gets destroyed
-! !See also note in subroutine
-! call ludcmp(L,nmat,nmat,indx,d)
-! call lubksb(L,nmat,nmat,indx,Imat) !Imat now contains psi
-! do 35 i=1,nmat
-! yplot(i,n)=Imat(i)
-!35 continue
-! if (.false.) call MatCheck(Ltemp,Linv)
-!40 continue
-! call plotit(del,yplot,nmat,npts)
-!
-! Now that the user has an idea of the dispersion, do the full propagation problem
-! write (*,*)'Enter detuning of center frequency of the coupling pulse in MHz'
-! read (*,*)delc_prop !del_prop is the detuning used for the propagation
- delc_prop=0
- delc_prop=2*pi*1e6*delc_prop/beta !Now dimensionless
-
-! write (*,*)'Enter detuning of center frequency of the pump pulse in MHz'
-! read (*,*)del2_prop !del_prop is the detuning used for the propagation
- del2_prop=0;
- del2_prop=2*pi*1e6*del2_prop/beta !Now dimensionless
-! write (*,*)'Enter detuning of center frequency of the probe pulse in MHz'
-! read (*,*)del1_prop !del_prop is the detuning used for the propagation
- del1_prop=0.
- del1_prop=2*pi*1e6*del1_prop/beta !Now dimensionless
-
-
-! write (*,*) 'Enter pulse width in nsec'
-! read (*,*)tp
- tp=1e-6
- tp=beta*tp !Now dimensionless
-! write (*,*)'Enter length of cell in m'
-! read (*,*)Lcell
- Lcell=100;
- Lcell=beta*Lcell/c !Now dimensionless
- t_start=secnds(0.E0)
- write (*,*)'t_start = ',t_start
-! XXXXXXX
-!
-! Set up initial pulse.
-!
- tshift=2*tp
- tmax=Lcell !Length of time to pass cell (no c because we're dimensionless)
- dt=tmax/tpts
- dz=dt !(no c because we're dimensionless)
-! write (*,*)'tp = ',tp
- Nframe=zpts+int(6*tp/dt)+1 !Change the number 4 to anything you want to see longer pulse evolution
- if (Nframe.ge.Nframemax) write (*,*)'Error!!!!Nframe>Nframemax'
-! write (*,*)'Nframe,tpts = ',Nframe,tpts
-!
-! Initialize matrices
-!
- Omold=cmplx(0.,0.)
- Omold_vac=cmplx(0.,0.)
- tpeak=-1
- tpeak_vac=-1
- write (*,*)'Nframe= ', Nframe
-! pause
-! do 110 n=1,Nframe
-
- do 100 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.)
-
-100 continue
-!110 continue
-
-
-
-! 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
- b2=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)*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_v2.txt'
-! write (*,*)'Enter file name to save parameters'
-! read (*,3)fname
-3 format(a150)
- open(9,name=fname)
- write (9,133)Nframe,zpts,beta,NSkip,dt
-! write (*,*)'Nframe,zpts,beta,NSkip,dt'
-! write (*,133)Nframe,zpts,beta,NSkip,dt
-
-133 format(1x,i10,',',i5,',',f12.2,',',i5,',',f12.2)
- close (9)
-
- fname='Movie4level_v2.dat'
-! write (*,*)'Enter file name to save movie'
-! read (*,3)fname
-
- open(9,name=fname)
- fname='Movie4level_EndPoints_v2.dat'
-! write (*,*)'Enter file name to save endpoints'
-! read (*,3)fname
-
- open(10,name=fname)
-
- do 60 n=1,Nframe
-! write (*,*)'n = ',n
- t=float(n-1)*dt
- Om1(1)=Om1peak
- Om2(1)=Om2peak*exp(-(t-tshift)**2/(tp*tp))
- Omc(1)=Omcpeak
- Om_vac(1)=Om2(1)
- if (int(n/10).eq.0) write(fname,130)'Movie',n
- if (int(n/10).ge.1.and.int(n/100).eq.0) write (fname,131)'Movie',n
- if (int(n/10).ge.1.and.int(n/100).gt.0) write (fname,132)'Movie',n
-130 format(a5,i1)
-131 format(a5,i2)
-132 format(a5,i3)
-! write (*,125)fname
-125 format(1x,a12)
-! open(9,name=fname)
-
- do 345 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)
-345 continue
-
- do 50 m=zpts,2,-1
- z(m)=float(m)*dz
-
- 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)
-
-
- Om_vac(m)=a1*Om_vac(m-1)
-
- if (mod(n,Nskip).eq.0) write (9,120)z(m),Om2(m),Om_vac(m),Omc(m)
-
-50 continue
- if (cdabs(Om2(zpts)).gt.cdabs(Omold)) tpeak=t
- if (cdabs(Om_vac(zpts)).gt.cdabs(Omold_vac)) tpeak_vac=t
- write (10,139) t,cdabs(Om2(zpts)),cdabs(Om_vac(zpts)) !EndPoint File
- Omold=Om2(zpts)
- Omold_vac=Om_vac(zpts)
-
-60 continue
- close(9)
- close(10)
-139 format(1x,f12.6,',',F12.6,',',F12.6)
-! write (*,*)'Medium pulse out at ',tpeak/(beta*1e-6),' microseconds'
-! write (*,*)'Vacuum pulse out at ',tpeak_vac/(beta*1e-6),' microseconds'
- write (*,*)Omcpeak,(tpeak-tpeak_vac)/(beta*1e-9)
-120 format(1x,f12.6,',',f12.6,',',f12.6,',',f12.6,',',f12.6,',',f12.6,',',f12.6)
- t_end=secnds(0.E0)
- t_elapsed=t_end-t_start
- write(*,*)'T elapsed = ',t_elapsed
- stop
- end
-
-!
- subroutine IMatConstruct(Om,Imat)
-!
-! NOTE: This subroutine actually calculates -Imat because we need to solve LPsi=-Imat.
- implicit none
- integer i,nmat
- parameter (nmat=3)
- real*8 Om
- complex*8 Imat(nmat)
-
- Imat(1)=-cmplx(0.,-0.5*Om)
- Imat(2)=-cmplx(0.,0.5*Om)
- Imat(3)=cmplx(0.,0.)
- return
- end
-
-
-!
- subroutine LMatConstruct(Om,del,L)
-!
-
- implicit none
- integer nmat
- parameter (nmat=3)
- real*8 del,ga12,Om,W21
- complex*8 L(nmat,nmat)
- common/para/ga12,W21
-
-
- L(1,1)=cmplx (-ga12,-del)
- L(1,3)=cmplx(0.,Om)
- L(2,2)=cmplx(-ga12,del)
- L(2,3)=cmplx(0.,-Om)
- L(3,1)=cmplx(0.,0.5*Om)
- L(3,2)=cmplx(0.,-0.5*Om)
- L(3,3)=cmplx(-W21,0.)
- return
- end
-
-!
- subroutine MatCheck(L,Linv)
-!
-
- implicit none
- integer i,j,k,nmat
- parameter (nmat=3)
- complex*8 L(nmat,nmat),Linv(nmat,nmat),Res(nmat,nmat)
-
- write (*,*)'L = '
- do 10 i=1,nmat
- write (*,120)(L(i,j),j=1,nmat)
-10 continue
- write (*,*)'Linv = '
- do 20 i=1,nmat
- write (*,120)(Linv(i,j),j=1,nmat)
-20 continue
- write (*,*)'Res = '
- do 50 i=1,nmat
- do 40 j=1,nmat
- Res(i,j)=cmplx(0.,0.)
- do 30 k=1,nmat
- Res(i,j)=Res(i,j)+Linv(i,k)*L(k,j)
-30 continue
-40 continue
- write (*,120)(Res(i,j),j=1,nmat)
-50 continue
-120 format(1x,3(f8.4,'+i',f8.4,' '))
-
- return
- end
-
-!
- subroutine plotit(x,y,nmat,npts)
-!
-! See MATLAB routine that will do this plotting.
-
- implicit none
- integer i,n,nmat,npts
- real*8 x(npts)
- complex*8 y(nmat,npts)
-
- write (*,*)'For now, we are just going to write the file'
- open(9, FILE='TwoLevelPulseProp.txt')
- do 10 n=1,npts
- write (9,100)x(n),(y(i,n),i=1,3)
-! write(*,100)x(n),(y(i,n),i=1,3)
-10 continue
-100 format(1x,f9.6,',',3(f9.6,',',f9.6,','))
-
-
- return
- end
-
-
-
-
-
-!********************************************************
-! Numerical Recipes
-!********************************************************
-
-
- SUBROUTINE ludcmp(a,n,np,indx,d)
- implicit none
-
- INTEGER n,np,indx(n),NMAX
- REAL*8 d,TINY
- PARAMETER (NMAX=500,TINY=1.0e-20)
- INTEGER i,imax,j,k
- REAL aamax,dum,vv(NMAX)
-! My changed variables
- complex*8 sum,dum2
- complex*8 a(np,np)
-
-
- d=1.
- do 12 i=1,n
- aamax=0.
- do 11 j=1,n
- if (cabs(a(i,j)).gt.aamax) aamax=cabs(a(i,j))
-11 continue
- if (aamax.eq.0.) pause 'singular matrix in ludcmp'
- vv(i)=1./aamax
-12 continue
- do 19 j=1,n
- do 14 i=1,j-1
- sum=a(i,j)
- do 13 k=1,i-1
- sum=sum-a(i,k)*a(k,j)
-13 continue
- a(i,j)=sum
-14 continue
- aamax=0.
- do 16 i=j,n
- sum=a(i,j)
- do 15 k=1,j-1
- sum=sum-a(i,k)*a(k,j)
-15 continue
- a(i,j)=sum
- dum=vv(i)*cabs(sum)
- if (dum.ge.aamax) then
- imax=i
- aamax=dum
- endif
-16 continue
- if (j.ne.imax)then
- do 17 k=1,n
- dum2=a(imax,k)
- a(imax,k)=a(j,k)
- a(j,k)=dum2
-17 continue
- d=-d
- vv(imax)=vv(j)
- endif
- indx(j)=imax
- if (cabs(a(j,j)).eq.0.) a(j,j)=cmplx(TINY,TINY)
- if(j.ne.n)then
- dum2=1./a(j,j)
- do 18 i=j+1,n
- a(i,j)=a(i,j)*dum2
-18 continue
- endif
-19 continue
- return
- END
-
- SUBROUTINE lubksb(a,n,np,indx,b)
- implicit none
- INTEGER n,np,indx(n)
- INTEGER i,ii,j,ll
-
-! My changed variables
- complex*8 sum
- complex*8 b(n)
- complex*8 a(np,np)
- ii=0
- do 12 i=1,n
- ll=indx(i)
- sum=b(ll)
- b(ll)=b(i)
- if (ii.ne.0)then
- do 11 j=ii,i-1
- sum=sum-a(i,j)*b(j)
-11 continue
- else if (sum.ne.0.) then
- ii=i
- endif
- b(i)=sum
-12 continue
- do 14 i=n,1,-1
- sum=b(i)
- do 13 j=i+1,n
- sum=sum-a(i,j)*b(j)
-13 continue
- b(i)=sum/a(i,i)
-14 continue
- return
- END
-
+ 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.
+
+ character*150 fname
+ integer nmat,npts,Nfrac,Nframe,Nframemax,NSkip,NWrite,tpts,zpts
+ parameter (nmat=3,npts=100) !matrix size, number of detuning points in dispersion curve
+ !REMEMBER TO CHANGE NMAT IN LMatConstruct Routine
+ parameter (tpts=10,zpts=tpts+1) !Caution: funny things happened when tpts=200 (and presumably greater)
+ !tpts is the number of temporal points in the cell
+ parameter (Nframemax=2000000)
+ parameter ( NWrite=100) !number of frames to actually write
+ integer i,j,k,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,beta,c,delmax,del1_prop,del2_prop,delc_prop
+ real*8 dt,dz,eta
+ real*8 W12,W21,W31,W32,W41,W42,W43,W34,ga12,ga13,ga14,ga23,ga24,ga34
+ real*8 Ga2,Ga4,Om_crit
+ real*8 Lcell,Om1peak,Om2peak,Omcpeak,pi,tmax,tp,tshift,t_end,t_start,t_elapsed
+ real*8 tpeak,tpeak_vac
+ real*8 epsil,hbar,lambda
+ real*8 del(npts)
+ real*8 t,z(zpts)
+ complex*16 yplot(nmat,npts)
+ complex*16 Imat(nmat)
+
+
+ complex*16 Om1(zpts),Om2(zpts),Omc(zpts),Om_vac(zpts)
+ complex*16 rho11(zpts),rho12(zpts),rho13(zpts),rho14(zpts),rho21(zpts),rho22(zpts),rho23(zpts),rho24(zpts)
+ complex*16 rho31(zpts),rho32(zpts),rho33(zpts),rho34(zpts),rho41(zpts),rho42(zpts),rho43(zpts),rho44(zpts)
+ complex*16 rho11_last(zpts),rho12_last(zpts),rho13_last(zpts),rho14_last(zpts)
+ complex*16 rho21_last(zpts),rho22_last(zpts),rho23_last(zpts),rho24_last(zpts)
+ complex*16 rho31_last(zpts),rho32_last(zpts),rho33_last(zpts),rho34_last(zpts)
+ complex*16 rho41_last(zpts),rho42_last(zpts),rho43_last(zpts),rho44_last(zpts)
+
+ !No Om_last because we never need the previous spatial point
+ complex*16 L(nmat,nmat),Linv(nmat,nmat),Ltemp(nmat,nmat)
+
+ common/para/ga12,W21
+
+
+ real*8 d !used by NR Routines
+ integer indx(nmat)
+
+!
+! Fundamental numbers
+!
+ ci=cmplx(0.,1.)
+ pi=acos(-1.0)
+ c=3e8
+ hbar=1.055e-34
+ epsil=8.85e-12
+
+!
+! Atomic numbers (based on Rubidium 85)
+!
+ beta=2*pi*3e6 !in Hz
+ W41=0
+ W42=0.
+ 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)
+ lambda=780.24e-9
+
+!
+! Atomic parameters
+!
+ write (*,*)'New version from moved folder'
+! write (*,*)'Enter density in m^-1'
+! read (*,*)eta
+ eta=6.9e11
+ alpha1=3*eta*lambda*lambda/(2*pi)
+ alpha1tilde=alpha1*c/beta
+ alpha2=3*eta*lambda*lambda/(2*pi)
+ alpha2tilde=alpha2*c/beta
+ alphac=3*eta*lambda*lambda/(2*pi)
+ alphactilde=alphac*c/beta
+
+!
+! Initialize matrices and set up Identity Matrix
+!
+ do 20 i=1,nmat
+ do 10 j=1,nmat
+ L(i,j)=cmplx(0.,0.)
+ Linv(i,j)=cmplx(0.,0.)
+10 continue !j loop
+ Imat(i)=cmplx(0.,0.)
+ Linv(i,i)=cmplx(1.,0.) !contains identity matrix
+20 continue
+
+
+!
+! User defined numbers
+!
+! write (*,*)'Enter peak scaled Rabi frequency for the pump at entrance of cell'
+! read (*,*)Om1peak
+ Om1peak=1
+! write (*,*)'Enter peak scaled Rabi frequency for the probe at entrance of cell'
+! read (*,*)Om2peak
+ Om2peak=1e-2
+! write (*,*)'Enter maximum detuning in MHz for dispersion lineshape plot'
+! read (*,*) delmax
+ 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
+
+! write (*,*)'Enter peak scaled Rabi frequency for the coupling field at entrance of cell'
+! read (*,*)Omcpeak
+ Omcpeak=0.
+
+!
+! First plot the dispersion lineshape
+!
+! do 40 n=1,npts
+! del(n)=-delmax+2*float(n)*delmax/npts
+! call LMatConstruct(Ompeak,del(n),L) !construct the L matrix
+! call LMatConstruct(Ompeak,del(n),Ltemp) !Need a temporary because L gets destroyed
+! call ImatConstruct(Ompeak,Imat) !Need to call in loop because it gets destroyed
+! !See also note in subroutine
+! call ludcmp(L,nmat,nmat,indx,d)
+! call lubksb(L,nmat,nmat,indx,Imat) !Imat now contains psi
+! do 35 i=1,nmat
+! yplot(i,n)=Imat(i)
+!35 continue
+! if (.false.) call MatCheck(Ltemp,Linv)
+!40 continue
+! call plotit(del,yplot,nmat,npts)
+!
+! Now that the user has an idea of the dispersion, do the full propagation problem
+! write (*,*)'Enter detuning of center frequency of the coupling pulse in MHz'
+! read (*,*)delc_prop !del_prop is the detuning used for the propagation
+ delc_prop=0
+ delc_prop=2*pi*1e6*delc_prop/beta !Now dimensionless
+
+! write (*,*)'Enter detuning of center frequency of the pump pulse in MHz'
+! read (*,*)del2_prop !del_prop is the detuning used for the propagation
+ del2_prop=0;
+ del2_prop=2*pi*1e6*del2_prop/beta !Now dimensionless
+! write (*,*)'Enter detuning of center frequency of the probe pulse in MHz'
+! read (*,*)del1_prop !del_prop is the detuning used for the propagation
+ del1_prop=0.
+ del1_prop=2*pi*1e6*del1_prop/beta !Now dimensionless
+
+
+! write (*,*) 'Enter pulse width in nsec'
+! read (*,*)tp
+ tp=1e-6
+ tp=beta*tp !Now dimensionless
+! write (*,*)'Enter length of cell in m'
+! read (*,*)Lcell
+ Lcell=100;
+ Lcell=beta*Lcell/c !Now dimensionless
+ t_start=secnds(0.E0)
+ write (*,*)'t_start = ',t_start
+! XXXXXXX
+!
+! Set up initial pulse.
+!
+ tshift=2*tp
+ tmax=Lcell !Length of time to pass cell (no c because we're dimensionless)
+ dt=tmax/tpts
+ dz=dt !(no c because we're dimensionless)
+! write (*,*)'tp = ',tp
+ Nframe=zpts+int(6*tp/dt)+1 !Change the number 4 to anything you want to see longer pulse evolution
+ if (Nframe.ge.Nframemax) write (*,*)'Error!!!!Nframe>Nframemax'
+! write (*,*)'Nframe,tpts = ',Nframe,tpts
+!
+! Initialize matrices
+!
+ Omold=cmplx(0.,0.)
+ Omold_vac=cmplx(0.,0.)
+ tpeak=-1
+ tpeak_vac=-1
+ write (*,*)'Nframe= ', Nframe
+! pause
+! do 110 n=1,Nframe
+
+ do 100 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.)
+
+100 continue
+!110 continue
+
+
+
+! 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
+ b2=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)*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_v2.txt'
+! write (*,*)'Enter file name to save parameters'
+! read (*,3)fname
+3 format(a150)
+ open(9,name=fname)
+ write (9,133)Nframe,zpts,beta,NSkip,dt
+! write (*,*)'Nframe,zpts,beta,NSkip,dt'
+! write (*,133)Nframe,zpts,beta,NSkip,dt
+
+133 format(1x,i10,',',i5,',',f12.2,',',i5,',',f12.2)
+ close (9)
+
+ fname='Movie4level_v2.dat'
+! write (*,*)'Enter file name to save movie'
+! read (*,3)fname
+
+ open(9,name=fname)
+ fname='Movie4level_EndPoints_v2.dat'
+! write (*,*)'Enter file name to save endpoints'
+! read (*,3)fname
+
+ open(10,name=fname)
+
+ do 60 n=1,Nframe
+! write (*,*)'n = ',n
+ t=float(n-1)*dt
+ Om1(1)=Om1peak
+ Om2(1)=Om2peak*exp(-(t-tshift)**2/(tp*tp))
+ Omc(1)=Omcpeak
+ Om_vac(1)=Om2(1)
+ if (int(n/10).eq.0) write(fname,130)'Movie',n
+ if (int(n/10).ge.1.and.int(n/100).eq.0) write (fname,131)'Movie',n
+ if (int(n/10).ge.1.and.int(n/100).gt.0) write (fname,132)'Movie',n
+130 format(a5,i1)
+131 format(a5,i2)
+132 format(a5,i3)
+! write (*,125)fname
+125 format(1x,a12)
+! open(9,name=fname)
+
+ do 345 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)
+345 continue
+
+ do 50 m=zpts,2,-1
+ z(m)=float(m)*dz
+
+ 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)
+
+
+ Om_vac(m)=a1*Om_vac(m-1)
+
+ if (mod(n,Nskip).eq.0) write (9,120)z(m),Om2(m),Om_vac(m),Omc(m)
+
+50 continue
+ if (cdabs(Om2(zpts)).gt.cdabs(Omold)) tpeak=t
+ if (cdabs(Om_vac(zpts)).gt.cdabs(Omold_vac)) tpeak_vac=t
+ write (10,139) t,cdabs(Om2(zpts)),cdabs(Om_vac(zpts)) !EndPoint File
+ Omold=Om2(zpts)
+ Omold_vac=Om_vac(zpts)
+
+60 continue
+ close(9)
+ close(10)
+139 format(1x,f12.6,',',F12.6,',',F12.6)
+! write (*,*)'Medium pulse out at ',tpeak/(beta*1e-6),' microseconds'
+! write (*,*)'Vacuum pulse out at ',tpeak_vac/(beta*1e-6),' microseconds'
+ write (*,*)Omcpeak,(tpeak-tpeak_vac)/(beta*1e-9)
+120 format(1x,f12.6,',',f12.6,',',f12.6,',',f12.6,',',f12.6,',',f12.6,',',f12.6)
+ t_end=secnds(0.E0)
+ t_elapsed=t_end-t_start
+ write(*,*)'T elapsed = ',t_elapsed
+ stop
+ end
+
+!
+ subroutine IMatConstruct(Om,Imat)
+!
+! NOTE: This subroutine actually calculates -Imat because we need to solve LPsi=-Imat.
+ implicit none
+ integer i,nmat
+ parameter (nmat=3)
+ real*8 Om
+ complex*8 Imat(nmat)
+
+ Imat(1)=-cmplx(0.,-0.5*Om)
+ Imat(2)=-cmplx(0.,0.5*Om)
+ Imat(3)=cmplx(0.,0.)
+ return
+ end
+
+
+!
+ subroutine LMatConstruct(Om,del,L)
+!
+
+ implicit none
+ integer nmat
+ parameter (nmat=3)
+ real*8 del,ga12,Om,W21
+ complex*8 L(nmat,nmat)
+ common/para/ga12,W21
+
+
+ L(1,1)=cmplx (-ga12,-del)
+ L(1,3)=cmplx(0.,Om)
+ L(2,2)=cmplx(-ga12,del)
+ L(2,3)=cmplx(0.,-Om)
+ L(3,1)=cmplx(0.,0.5*Om)
+ L(3,2)=cmplx(0.,-0.5*Om)
+ L(3,3)=cmplx(-W21,0.)
+ return
+ end
+
+!
+ subroutine MatCheck(L,Linv)
+!
+
+ implicit none
+ integer i,j,k,nmat
+ parameter (nmat=3)
+ complex*8 L(nmat,nmat),Linv(nmat,nmat),Res(nmat,nmat)
+
+ write (*,*)'L = '
+ do 10 i=1,nmat
+ write (*,120)(L(i,j),j=1,nmat)
+10 continue
+ write (*,*)'Linv = '
+ do 20 i=1,nmat
+ write (*,120)(Linv(i,j),j=1,nmat)
+20 continue
+ write (*,*)'Res = '
+ do 50 i=1,nmat
+ do 40 j=1,nmat
+ Res(i,j)=cmplx(0.,0.)
+ do 30 k=1,nmat
+ Res(i,j)=Res(i,j)+Linv(i,k)*L(k,j)
+30 continue
+40 continue
+ write (*,120)(Res(i,j),j=1,nmat)
+50 continue
+120 format(1x,3(f8.4,'+i',f8.4,' '))
+
+ return
+ end
+
+!
+ subroutine plotit(x,y,nmat,npts)
+!
+! See MATLAB routine that will do this plotting.
+
+ implicit none
+ integer i,n,nmat,npts
+ real*8 x(npts)
+ complex*8 y(nmat,npts)
+
+ write (*,*)'For now, we are just going to write the file'
+ open(9, FILE='TwoLevelPulseProp.txt')
+ do 10 n=1,npts
+ write (9,100)x(n),(y(i,n),i=1,3)
+! write(*,100)x(n),(y(i,n),i=1,3)
+10 continue
+100 format(1x,f9.6,',',3(f9.6,',',f9.6,','))
+
+
+ return
+ end
+
+
+
+
+
+!********************************************************
+! Numerical Recipes
+!********************************************************
+
+
+ SUBROUTINE ludcmp(a,n,np,indx,d)
+ implicit none
+
+ INTEGER n,np,indx(n),NMAX
+ REAL*8 d,TINY
+ PARAMETER (NMAX=500,TINY=1.0e-20)
+ INTEGER i,imax,j,k
+ REAL aamax,dum,vv(NMAX)
+! My changed variables
+ complex*8 sum,dum2
+ complex*8 a(np,np)
+
+
+ d=1.
+ do 12 i=1,n
+ aamax=0.
+ do 11 j=1,n
+ if (cabs(a(i,j)).gt.aamax) aamax=cabs(a(i,j))
+11 continue
+ if (aamax.eq.0.) pause 'singular matrix in ludcmp'
+ vv(i)=1./aamax
+12 continue
+ do 19 j=1,n
+ do 14 i=1,j-1
+ sum=a(i,j)
+ do 13 k=1,i-1
+ sum=sum-a(i,k)*a(k,j)
+13 continue
+ a(i,j)=sum
+14 continue
+ aamax=0.
+ do 16 i=j,n
+ sum=a(i,j)
+ do 15 k=1,j-1
+ sum=sum-a(i,k)*a(k,j)
+15 continue
+ a(i,j)=sum
+ dum=vv(i)*cabs(sum)
+ if (dum.ge.aamax) then
+ imax=i
+ aamax=dum
+ endif
+16 continue
+ if (j.ne.imax)then
+ do 17 k=1,n
+ dum2=a(imax,k)
+ a(imax,k)=a(j,k)
+ a(j,k)=dum2
+17 continue
+ d=-d
+ vv(imax)=vv(j)
+ endif
+ indx(j)=imax
+ if (cabs(a(j,j)).eq.0.) a(j,j)=cmplx(TINY,TINY)
+ if(j.ne.n)then
+ dum2=1./a(j,j)
+ do 18 i=j+1,n
+ a(i,j)=a(i,j)*dum2
+18 continue
+ endif
+19 continue
+ return
+ END
+
+ SUBROUTINE lubksb(a,n,np,indx,b)
+ implicit none
+ INTEGER n,np,indx(n)
+ INTEGER i,ii,j,ll
+
+! My changed variables
+ complex*8 sum
+ complex*8 b(n)
+ complex*8 a(np,np)
+ ii=0
+ do 12 i=1,n
+ ll=indx(i)
+ sum=b(ll)
+ b(ll)=b(i)
+ if (ii.ne.0)then
+ do 11 j=ii,i-1
+ sum=sum-a(i,j)*b(j)
+11 continue
+ else if (sum.ne.0.) then
+ ii=i
+ endif
+ b(i)=sum
+12 continue
+ do 14 i=n,1,-1
+ sum=b(i)
+ do 13 j=i+1,n
+ sum=sum-a(i,j)*b(j)
+13 continue
+ b(i)=sum/a(i,i)
+14 continue
+ return
+ END
+
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-16 11:06:59 +0000