dos 2 unix encoding

5 files changed, 2119 insertions, 2119 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
+      
diff --git a/code_from_navy/fortran/FourLevelPrograms/FourLevelPulseProp_v3_Double.f b/code_from_navy/fortran/FourLevelPrograms/FourLevelPulseProp_v3_Double.f
index 0d9547a..2a08309 100644
--- a/code_from_navy/fortran/FourLevelPrograms/FourLevelPulseProp_v3_Double.f
+++ b/code_from_navy/fortran/FourLevelPrograms/FourLevelPulseProp_v3_Double.f
@@ -1,704 +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

-	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_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

-		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,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)

-

-	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.) 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=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
+	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_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
+		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,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)
+
+	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.) 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
+      
diff --git a/code_from_navy/fortran/FourLevelPrograms/FourLevelPulseProp_v3_Double.f.bak b/code_from_navy/fortran/FourLevelPrograms/FourLevelPulseProp_v3_Double.f.bak
index cda5337..905c861 100644
--- a/code_from_navy/fortran/FourLevelPrograms/FourLevelPulseProp_v3_Double.f.bak
+++ b/code_from_navy/fortran/FourLevelPrograms/FourLevelPulseProp_v3_Double.f.bak
@@ -1,703 +1,703 @@
-	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

-	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

-

-!

-!	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_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

-		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,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)

-

-	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.) 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=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
+	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
+
+!
+!	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_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
+		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,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)
+
+	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.) 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
+      
diff --git a/code_from_navy/fortran/FourLevelPrograms/compilefourlevelpulseprop_v2_double.bat b/code_from_navy/fortran/FourLevelPrograms/compilefourlevelpulseprop_v2_double.bat
index ec2e6a3..76a6b87 100644
--- a/code_from_navy/fortran/FourLevelPrograms/compilefourlevelpulseprop_v2_double.bat
+++ b/code_from_navy/fortran/FourLevelPrograms/compilefourlevelpulseprop_v2_double.bat
@@ -1,2 +1,2 @@
-call c:/g77/g77setup.bat

-C:/g77/bin/G77 C:\FortranFiles\FourLevel\FourLevelPulseProp_v2_double.F -lslatec -o C:\FortranFiles\FourLevel\FourLevelPulseProp_v2_double.exe
\ No newline at end of file
+call c:/g77/g77setup.bat
+C:/g77/bin/G77 C:\FortranFiles\FourLevel\FourLevelPulseProp_v2_double.F -lslatec -o C:\FortranFiles\FourLevel\FourLevelPulseProp_v2_double.exe
diff --git a/code_from_navy/fortran/FourLevelPrograms/compilefourlevelpulseprop_v3_double.bat b/code_from_navy/fortran/FourLevelPrograms/compilefourlevelpulseprop_v3_double.bat
index a4dd2a8..0c7c618 100644
--- a/code_from_navy/fortran/FourLevelPrograms/compilefourlevelpulseprop_v3_double.bat
+++ b/code_from_navy/fortran/FourLevelPrograms/compilefourlevelpulseprop_v3_double.bat
@@ -1,2 +1,2 @@
-call c:/g77/g77setup.bat

-C:/g77/bin/G77 C:\FortranFiles\FourLevel\FourLevelPulseProp_v3_double.F -lslatec -o C:FourLevelPulseProp_v3_double.exe
\ No newline at end of file
+call c:/g77/g77setup.bat
+C:/g77/bin/G77 C:\FortranFiles\FourLevel\FourLevelPulseProp_v3_double.F -lslatec -o C:FourLevelPulseProp_v3_double.exe