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=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 (*,*)'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=0.01
!	write (*,*)'Enter peak scaled Rabi frequency for the probe at entrance of cell'
!	read (*,*)Om2peak
	Om2peak=1
!	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.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=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
        write (*,*) 'peak center at the cell begining i.e. tshift = ', tshift
	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
	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
	b3=0.5*ci*alpha2tilde*dt
	b4=0.
	b5=0.
	b6=0.
	b7=0.
	
	c1=1.
	c2=0.5*ci*alphactilde*dt
	c3=0.5*ci*alphactilde*dt
	c4=0.
	c5=0.
	c6=0.
	
	d1=1-(ga12-ci*(del2_prop-del1_prop))*dt
	d2=0.25*ci*dt
	d3=-0.25*ci*dt
	d4=-0.25*ci*dt
	d5=0.25*ci*dt
	d6=-0.25*ci*dt
	d7=-0.25*ci*dt
	
	e1=1-(ga13+ci*del1_prop)*dt
	e2=0.25*ci*dt
	e3=-0.25*ci*dt
	e4=0.25*ci*dt
	e5=-0.25*ci*dt
	
	f1=1-(ga14-ci*(del2_prop-del1_prop-delc_prop))*dt
	f2=0.25*ci*dt
	f3=-0.25*ci*dt
	f4=0.25*ci*dt
	f5=-0.25*ci*dt
	
	g1=1-(ga23+ci*del1_prop)*dt
	g2=-0.25*ci*dt
	g3=0.25*ci*dt
	g4=0.25*ci*dt
	g5=-0.25*ci*dt
	g6=0.25*ci*dt
	g7=0.25*ci*dt
	
	h1=1-(ga24+ci*delc_prop)*dt
	h2=0.25*ci*dt
	h3=0.25*ci*dt
	h4=0.25*ci*dt
	h5=0.25*ci*dt
	
	i1=1-(ga34-ci*(del2_prop-delc_prop))*dt
	i2=-0.25*ci*dt
	i3=0.25*ci*dt
	i4=0.25*ci*dt
	i5=-0.25*ci*dt
	i6=0.25*ci*dt
	i7=0.25*ci*dt
	
	j1=1-W12*dt
	j2=W21*dt
	j3=W31*dt
	j4=0.25*ci*dt
	j5=-0.25*ci*dt
	j6=0.25*ci*dt
	j7=-0.25*ci*dt
	j8=W41*dt
	
	k1=1-(W32+W31+W34)*dt
	k2=W43*dt
	k3=-0.25*ci*dt
	k4=-0.25*ci*dt
	k5=0.25*ci*dt
	k6=0.25*ci*dt
	k7=-0.25*ci*dt
	k8=-0.25*ci*dt
	k9=0.25*ci*dt
	k10=0.25*ci*dt
	
	l1=1-(W43+W42+W41)*dt
	l2=-0.25*ci*dt
	l3=0.25*ci*dt
	l4=-0.25*ci*dt
	l5=0.25*ci*dt

	NSkip=int(NFrame/NWrite)	
							
	fname='MovieParameters4level.txt'
!	write (*,*)'Enter file name to save parameters'
!	read (*,3)fname
3	format(a150)
	open(9,file=fname)
	write (9,133)Nframe,zpts,beta,NSkip,dt
!	write (*,*)'Nframe,zpts,beta,NSkip,dt'
!	write (*,133)Nframe,zpts,beta,NSkip,dt
	
133	format(1x,i10,',',i5,',',f12.2,',',i5,',',f12.2)
	close (9)
	
	fname='Movie4level.dat'
!	write (*,*)'Enter file name to save movie'
!	read (*,3)fname

	open(9,file=fname)
	fname='Movie4level_EndPoints.dat'
!	write (*,*)'Enter file name to save endpoints'
!	read (*,3)fname

	open(10,file=fname)

	do 60 n=1,Nframe
!		write (*,*)'n = ',n
		t=float(n-1)*dt
		Om1(1)=Om1peak*exp(-(t-tshift)**2/(tp*tp))	
		Om2(1)=Om2peak
		Omc(1)=Omcpeak
		Om_vac(1)=Om1(1)
		if (int(n/10).eq.0) write(fname,130)'Movie',n
		if (int(n/10).ge.1.and.int(n/100).eq.0) write (fname,131)'Movie',n
		if (int(n/10).ge.1.and.int(n/100).gt.0) write (fname,132)'Movie',n
130	format(a5,i1)
131	format(a5,i2)
132	format(a5,i3)
!		write (*,125)fname
125	format(1x,a12)
!		open(9,file=fname)

		do 345 m=1,zpts
			rho11_last(m)=rho11(m)
			rho12_last(m)=rho12(m)
			rho13_last(m)=rho13(m)
			rho14_last(m)=rho14(m)
			rho21_last(m)=rho21(m)
			rho22_last(m)=rho22(m)
			rho23_last(m)=rho23(m)
			rho24_last(m)=rho24(m)
			rho31_last(m)=rho31(m)
			rho32_last(m)=rho32(m)
			rho33_last(m)=rho33(m)
			rho34_last(m)=rho34(m)
			rho41_last(m)=rho41(m)
			rho42_last(m)=rho42(m)
			rho43_last(m)=rho43(m)
			rho44_last(m)=rho44(m)
345	continue

		do 50 m=zpts,2,-1
			z(m)=float(m)*dz
			
	Om1(m)=a1*Om1(m-1)+a2*rho31_last(m)+a3*rho31_last(m-1)

	Om2(m)=b1*Om2(m-1)+b2*rho32_last(m)+b3*rho32_last(m-1)
	
	Omc(m)=c1*Omc(m-1)+c2*rho42_last(m)+c3*rho42_last(m-1)

	rho11(m)=j1*rho11_last(m)+j2*rho22_last(m)+j3*rho33_last(m)+j4*conjg(Om1(m))*rho31_last(m)
	rho11(m)=rho11(m)+j5*Om1(m)*rho13_last(m)+j6*conjg(Om1(m-1))*rho31_last(m)+j7*Om1(m-1)*rho13_last(m)
	rho11(m)=rho11(m)+j8*rho44_last(m)
	
	rho12(m)=d1*rho12_last(m)+d2*conjg(Om1(m))*rho32_last(m)+d3*Om2(m)*rho13_last(m)
	rho12(m)=rho12(m)+d4*Omc(m)*rho14_last(m)+d5*conjg(Om1(m-1))*rho32_last(m)
	rho12(m)=rho12(m)+d6*Om2(m-1)*rho13_last(m)+d7*Omc(m-1)*rho14_last(m)

	rho13(m)=e1*rho13_last(m)+e2*conjg(Om1(m))*(rho33_last(m)-rho11_last(m))
	rho13(m)=rho13(m)+e3*conjg(Om2(m))*rho12_last(m)+e4*conjg(Om1(m-1))*(rho33_last(m)-rho11_last(m))
	rho13(m)=rho13(m)+e5*conjg(Om2(m-1))*rho12_last(m)
	
	rho14(m)=f1*rho14_last(m)+f2*conjg(Om1(m))*rho34_last(m)+f3*conjg(Omc(m))*rho12_last(m)
	rho14(m)=rho14(m)+f4*conjg(Om1(m-1))*rho34_last(m)+f5*conjg(Omc(m-1))*rho12_last(m)
	
	rho21(m)=conjg(rho12(m))
	
!	rho22(m) needs to be calculated lower down

	rho23(m)=g1*rho23_last(m)+g2*conjg(Om1(m))*rho21_last(m)+g3*conjg(Om2(m))*(rho33_last(m)-rho22_last(m))
	rho23(m)=rho23(m)+g4*conjg(Omc(m))*rho43_last(m)+g5*conjg(Om1(m-1))*rho21_last(m)
	rho23(m)=rho23(m)+g6*conjg(Om2(m-1))*(rho33_last(m)-rho22_last(m))+g7*conjg(Omc(m-1))*rho43_last(m)
	
	rho24(m)=h1*rho24_last(m)+h2*conjg(Om2(m))*rho34_last(m)+h3*conjg(Omc(m))*(rho44_last(m)-rho22_last(m))
	rho24(m)=rho24(m)+h4*conjg(Om2(m-1))*rho34_last(m)+h5*conjg(Omc(m-1))*(rho44_last(m)-rho22_last(m))
	
	rho31(m)=conjg(rho13(m))
	
	rho32(m)=conjg(rho23(m))
	
	rho33(m)=k1*rho33_last(m)+k2*rho44_last(m)+k3*conjg(Om1(m))*rho31_last(m)+k4*conjg(Om2(m))*rho32_last(m)
	rho33(m)=rho33(m)+k5*Om1(m)*rho13_last(m)+k6*Om2(m)*rho23_last(m)+k7*conjg(Om1(m-1))*rho31_last(m)
	rho33(m)=rho33(m)+k8*conjg(Om2(m-1))*rho32_last(m)+k9*Om1(m-1)*rho13_last(m)+k10*Om2(m-1)*rho23_last(m)
	
	rho34(m)=i1*rho34_last(m)+i2*conjg(Omc(m))*rho32_last(m)+i3*Om1(m)*rho14_last(m)+i4*Om2(m)*rho24_last(m)
	rho34(m)=rho34(m)+i5*conjg(Omc(m-1))*rho32_last(m)+i6*Om1(m-1)*rho14_last(m)+i7*Om2(m-1)*rho24_last(m)
	
	rho41(m)=conjg(rho14(m))

	rho42(m)=conjg(rho24(m))

	rho43(m)=conjg(rho34(m))
	
	rho44(m)=l1*rho44_last(m)+l2*conjg(Omc(m))*rho42_last(m)+l3*Omc(m)*rho24_last(m)
	rho44(m)=rho44(m)+l4*conjg(Omc(m-1))*rho42_last(m)+l5*Omc(m-1)*rho24(m)
	
	rho22(m)=1-rho11(m)-rho33(m)-rho44(m)
	
	
	rho11_last(m)=rho11(m)
	rho12_last(m)=rho12(m)
	rho13_last(m)=rho13(m)
	rho14_last(m)=rho14(m)
	rho21_last(m)=rho21(m)
	rho22_last(m)=rho22(m)
	rho23_last(m)=rho23(m)
	rho24_last(m)=rho24(m)
	rho31_last(m)=rho31(m)
	rho32_last(m)=rho32(m)
	rho33_last(m)=rho33(m)
	rho34_last(m)=rho34(m)
	rho41_last(m)=rho41(m)
	rho42_last(m)=rho42(m)
	rho43_last(m)=rho43(m)
	rho44_last(m)=rho44(m)


	Om_vac(m)=a1*Om_vac(m-1)

			if (mod(n,Nskip).eq.0) write (9,120) t,z(m),cdabs(Om1(m)),cdabs(Om2(m)),cdabs(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(Om1(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(f12.6,',',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