summaryrefslogtreecommitdiff
path: root/code_from_navy/fortran95
diff options
context:
space:
mode:
Diffstat (limited to 'code_from_navy/fortran95')
-rw-r--r--code_from_navy/fortran95/FourLevelPulseProp_v2_Double.f95708
-rw-r--r--code_from_navy/fortran95/FourLevelPulseProp_v3_Double.f95704
-rw-r--r--code_from_navy/fortran95/Makefile21
3 files changed, 1433 insertions, 0 deletions
diff --git a/code_from_navy/fortran95/FourLevelPulseProp_v2_Double.f95 b/code_from_navy/fortran95/FourLevelPulseProp_v2_Double.f95
new file mode 100644
index 0000000..4b2e86d
--- /dev/null
+++ b/code_from_navy/fortran95/FourLevelPulseProp_v2_Double.f95
@@ -0,0 +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,file=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,file=fname)
+ fname='Movie4level_EndPoints_v2.dat'
+! write (*,*)'Enter file name to save endpoints'
+! read (*,3)fname
+
+ open(10,file=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,file=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.) write (*,*) '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
+
diff --git a/code_from_navy/fortran95/FourLevelPulseProp_v3_Double.f95 b/code_from_navy/fortran95/FourLevelPulseProp_v3_Double.f95
new file mode 100644
index 0000000..395ab17
--- /dev/null
+++ b/code_from_navy/fortran95/FourLevelPulseProp_v3_Double.f95
@@ -0,0 +1,704 @@
+ 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=100,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
+ 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,Ompeak
+ 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=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)
+ lambda=780.24e-9
+
+!
+! Atomic parameters
+!
+! write (*,*)'Enter density in m^-1'
+! read (*,*)eta
+ eta=6.9e13
+ 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=.01
+! 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
+ Ompeak=0.1
+
+!
+! 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=1;
+ Lcell=beta*Lcell/c !Now dimensionless
+ t_start=secnds(0.E0)
+! 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(4*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
+! 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=W12*dt
+ j3=W31*dt
+ j4=0.25*ci*dt
+ j5=-0.25*ci*dt
+ j6=0.25*ci*dt
+ j7=-0.25*ci*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_v3.txt'
+! write (*,*)'Enter file name to save parameters'
+! read (*,3)fname
+3 format(a150)
+ open(9,file=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_v3.dat'
+! write (*,*)'Enter file name to save movie'
+! read (*,3)fname
+
+ open(9,file=fname)
+ fname='Movie4level_EndPoints_v3.dat'
+! write (*,*)'Enter file name to save endpoints'
+! read (*,3)fname
+
+ open(10,file=fname)
+
+ do 60 n=1,Nframe
+ t=float(n-1)*dt
+ Om1(1)=Om1peak
+ Om2(1)=Om2peak*exp(-(t-tshift)**2/(tp*tp))
+ Omc(1)=Omcpeak*exp(-(t-tshift)**2/(tp*tp))
+ 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,file=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)
+
+ 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*Omc(m)*(rho44_last(m)-rho22_last(m))
+ rho24(m)=rho24(m)+h4*conjg(Om2(m-1))*rho34_last(m)+h5*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.) write (*,*) '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
+
diff --git a/code_from_navy/fortran95/Makefile b/code_from_navy/fortran95/Makefile
new file mode 100644
index 0000000..b55915f
--- /dev/null
+++ b/code_from_navy/fortran95/Makefile
@@ -0,0 +1,21 @@
+# -*- make -*-
+# FILE: "/home/evmik/src/my_src/Nresonances/fortran/navy_four_levels/Makefile"
+# LAST MODIFICATION: "Sun, 10 Jul 2011 11:55:16 -0400 (evmik)"
+# (C) 2011 by Eugeniy Mikhailov, <evgmik@gmail.com>
+# $Id:$
+
+f95files = $(wildcard *.f95)
+runfiles = $(f95files:%.f95=%.run)
+
+executables: $(runfiles)
+
+$(runfiles): %.run : %.f95
+ gfortran -o $@ $<
+
+clean:
+ rm -f $(runfiles)
+ rm -f Movie4level_EndPoints_v*.dat
+ rm -f Movie4level_v*.dat
+ rm -f MovieParameters4level_v*.txt
+
+
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-16 11:30:30 +0000