From ba7aadf317232c69a71d40dacc80772c11cbc297 Mon Sep 17 00:00:00 2001
From: Eugeniy Mikhailov <evgmik@gmail.com>
Date: Wed, 6 Jul 2011 15:17:16 -0400
Subject: dos 2 unix encoding

---
 .../FourLevelPulseProp_v2_Double.f                 | 1416 ++++++++++----------
 .../FourLevelPulseProp_v3_Double.f                 | 1408 +++++++++----------
 .../FourLevelPulseProp_v3_Double.f.bak             | 1406 +++++++++----------
 .../compilefourlevelpulseprop_v2_double.bat        |    4 +-
 .../compilefourlevelpulseprop_v3_double.bat        |    4 +-
 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
-- 
cgit v1.2.3