diff options
| -rwxr-xr-x[-rw-r--r--] | fortran/navy_four_levels/FourLevelPulseProp_Double.f95 (renamed from fortran/navy_four_levels/FourLevelPulseProp_v2_Double.f95) | 10 | ||||
| -rw-r--r-- | fortran/navy_four_levels/FourLevelPulseProp_v3_Double.f95 | 708 | ||||
| -rw-r--r-- | fortran/navy_four_levels/Makefile | 6 |
3 files changed, 8 insertions, 716 deletions
diff --git a/fortran/navy_four_levels/FourLevelPulseProp_v2_Double.f95 b/fortran/navy_four_levels/FourLevelPulseProp_Double.f95 index 045f4a3..06f7353 100644..100755 --- a/fortran/navy_four_levels/FourLevelPulseProp_v2_Double.f95 +++ b/fortran/navy_four_levels/FourLevelPulseProp_Double.f95 @@ -347,7 +347,7 @@ NSkip=int(NFrame/NWrite) - fname='MovieParameters4level_v2.txt' + fname='MovieParameters4level.txt' ! write (*,*)'Enter file name to save parameters' ! read (*,3)fname 3 format(a150) @@ -359,12 +359,12 @@ 133 format(1x,i10,',',i5,',',f12.2,',',i5,',',f12.2) close (9) - fname='Movie4level_v2.dat' + fname='Movie4level.dat' ! write (*,*)'Enter file name to save movie' ! read (*,3)fname open(9,file=fname) - fname='Movie4level_EndPoints_v2.dat' + fname='Movie4level_EndPoints.dat' ! write (*,*)'Enter file name to save endpoints' ! read (*,3)fname @@ -484,7 +484,7 @@ 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) + if (mod(n,Nskip).eq.0) write (9,120) t,z(m),Om2(m),Om_vac(m),Omc(m) 50 continue if (cdabs(Om2(zpts)).gt.cdabs(Omold)) tpeak=t @@ -500,7 +500,7 @@ ! 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) +120 format(1x,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 diff --git a/fortran/navy_four_levels/FourLevelPulseProp_v3_Double.f95 b/fortran/navy_four_levels/FourLevelPulseProp_v3_Double.f95 deleted file mode 100644 index 5672a42..0000000 --- a/fortran/navy_four_levels/FourLevelPulseProp_v3_Double.f95 +++ /dev/null @@ -1,708 +0,0 @@ - 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=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 - 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=1 - 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(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_v3.txt' -! write (*,*)'Enter file name to save parameters' -! read (*,3)fname -3 format(a150) - open(9,file=fname) - write (9,133)Nframe,zpts,beta,NSkip,dt -! write (*,*)'Nframe,zpts,beta,NSkip,dt' -! write (*,133)Nframe,zpts,beta,NSkip,dt - -133 format(1x,i10,',',i5,',',f12.2,',',i5,',',f12.2) - close (9) - - fname='Movie4level_v3.dat' -! write (*,*)'Enter file name to save movie' -! read (*,3)fname - - open(9,file=fname) - fname='Movie4level_EndPoints_v3.dat' -! write (*,*)'Enter file name to save endpoints' -! read (*,3)fname - - open(10,file=fname) - - do 60 n=1,Nframe -! write (*,*)'n = ',n - t=float(n-1)*dt - Om1(1)=Om1peak - Om2(1)=Om2peak*exp(-(t-tshift)**2/(tp*tp)) - Omc(1)=Omcpeak - Om_vac(1)=Om2(1) - if (int(n/10).eq.0) write(fname,130)'Movie',n - if (int(n/10).ge.1.and.int(n/100).eq.0) write (fname,131)'Movie',n - if (int(n/10).ge.1.and.int(n/100).gt.0) write (fname,132)'Movie',n -130 format(a5,i1) -131 format(a5,i2) -132 format(a5,i3) -! write (*,125)fname -125 format(1x,a12) -! open(9,file=fname) - - do 345 m=1,zpts - rho11_last(m)=rho11(m) - rho12_last(m)=rho12(m) - rho13_last(m)=rho13(m) - rho14_last(m)=rho14(m) - rho21_last(m)=rho21(m) - rho22_last(m)=rho22(m) - rho23_last(m)=rho23(m) - rho24_last(m)=rho24(m) - rho31_last(m)=rho31(m) - rho32_last(m)=rho32(m) - rho33_last(m)=rho33(m) - rho34_last(m)=rho34(m) - rho41_last(m)=rho41(m) - rho42_last(m)=rho42(m) - rho43_last(m)=rho43(m) - rho44_last(m)=rho44(m) -345 continue - - do 50 m=zpts,2,-1 - z(m)=float(m)*dz - - Om1(m)=a1*Om1(m-1)+a2*rho31_last(m)+a3*rho31_last(m-1) - - Om2(m)=b1*Om2(m-1)+b2*rho32_last(m)+b3*rho32_last(m-1) - - Omc(m)=c1*Omc(m-1)+c2*rho42_last(m)+c3*rho42_last(m-1) - - rho11(m)=j1*rho11_last(m)+j2*rho22_last(m)+j3*rho33_last(m)+j4*conjg(Om1(m))*rho31_last(m) - rho11(m)=rho11(m)+j5*Om1(m)*rho13_last(m)+j6*conjg(Om1(m-1))*rho31_last(m)+j7*Om1(m-1)*rho13_last(m) - rho11(m)=rho11(m)+j8*rho44_last(m) - - rho12(m)=d1*rho12_last(m)+d2*conjg(Om1(m))*rho32_last(m)+d3*Om2(m)*rho13_last(m) - rho12(m)=rho12(m)+d4*Omc(m)*rho14_last(m)+d5*conjg(Om1(m-1))*rho32_last(m) - rho12(m)=rho12(m)+d6*Om2(m-1)*rho13_last(m)+d7*Omc(m-1)*rho14_last(m) - - rho13(m)=e1*rho13_last(m)+e2*conjg(Om1(m))*(rho33_last(m)-rho11_last(m)) - rho13(m)=rho13(m)+e3*conjg(Om2(m))*rho12_last(m)+e4*conjg(Om1(m-1))*(rho33_last(m)-rho11_last(m)) - rho13(m)=rho13(m)+e5*conjg(Om2(m-1))*rho12_last(m) - - rho14(m)=f1*rho14_last(m)+f2*conjg(Om1(m))*rho34_last(m)+f3*conjg(Omc(m))*rho12_last(m) - rho14(m)=rho14(m)+f4*conjg(Om1(m-1))*rho34_last(m)+f5*conjg(Omc(m-1))*rho12_last(m) - - rho21(m)=conjg(rho12(m)) - -! rho22(m) needs to be calculated lower down - - rho23(m)=g1*rho23_last(m)+g2*conjg(Om1(m))*rho21_last(m)+g3*conjg(Om2(m))*(rho33_last(m)-rho22_last(m)) - rho23(m)=rho23(m)+g4*conjg(Omc(m))*rho43_last(m)+g5*conjg(Om1(m-1))*rho21_last(m) - rho23(m)=rho23(m)+g6*conjg(Om2(m-1))*(rho33_last(m)-rho22_last(m))+g7*conjg(Omc(m-1))*rho43_last(m) - - rho24(m)=h1*rho24_last(m)+h2*conjg(Om2(m))*rho34_last(m)+h3*conjg(Omc(m))*(rho44_last(m)-rho22_last(m)) - rho24(m)=rho24(m)+h4*conjg(Om2(m-1))*rho34_last(m)+h5*conjg(Omc(m-1))*(rho44_last(m)-rho22_last(m)) - - rho31(m)=conjg(rho13(m)) - - rho32(m)=conjg(rho23(m)) - - rho33(m)=k1*rho33_last(m)+k2*rho44_last(m)+k3*conjg(Om1(m))*rho31_last(m)+k4*conjg(Om2(m))*rho32_last(m) - rho33(m)=rho33(m)+k5*Om1(m)*rho13_last(m)+k6*Om2(m)*rho23_last(m)+k7*conjg(Om1(m-1))*rho31_last(m) - rho33(m)=rho33(m)+k8*conjg(Om2(m-1))*rho32_last(m)+k9*Om1(m-1)*rho13_last(m)+k10*Om2(m-1)*rho23_last(m) - - rho34(m)=i1*rho34_last(m)+i2*conjg(Omc(m))*rho32_last(m)+i3*Om1(m)*rho14_last(m)+i4*Om2(m)*rho24_last(m) - rho34(m)=rho34(m)+i5*conjg(Omc(m-1))*rho32_last(m)+i6*Om1(m-1)*rho14_last(m)+i7*Om2(m-1)*rho24_last(m) - - rho41(m)=conjg(rho14(m)) - - rho42(m)=conjg(rho24(m)) - - rho43(m)=conjg(rho34(m)) - - rho44(m)=l1*rho44_last(m)+l2*conjg(Omc(m))*rho42_last(m)+l3*Omc(m)*rho24_last(m) - rho44(m)=rho44(m)+l4*conjg(Omc(m-1))*rho42_last(m)+l5*Omc(m-1)*rho24(m) - - rho22(m)=1-rho11(m)-rho33(m)-rho44(m) - - - rho11_last(m)=rho11(m) - rho12_last(m)=rho12(m) - rho13_last(m)=rho13(m) - rho14_last(m)=rho14(m) - rho21_last(m)=rho21(m) - rho22_last(m)=rho22(m) - rho23_last(m)=rho23(m) - rho24_last(m)=rho24(m) - rho31_last(m)=rho31(m) - rho32_last(m)=rho32(m) - rho33_last(m)=rho33(m) - rho34_last(m)=rho34(m) - rho41_last(m)=rho41(m) - rho42_last(m)=rho42(m) - rho43_last(m)=rho43(m) - rho44_last(m)=rho44(m) - - - Om_vac(m)=a1*Om_vac(m-1) - - if (mod(n,Nskip).eq.0) write (9,120)z(m),Om2(m),Om_vac(m),Omc(m) - -50 continue - if (cdabs(Om2(zpts)).gt.cdabs(Omold)) tpeak=t - if (cdabs(Om_vac(zpts)).gt.cdabs(Omold_vac)) tpeak_vac=t - write (10,139) t,cdabs(Om2(zpts)),cdabs(Om_vac(zpts)) !EndPoint File - Omold=Om2(zpts) - Omold_vac=Om_vac(zpts) - -60 continue - close(9) - close(10) -139 format(1x,f12.6,',',F12.6,',',F12.6) -! write (*,*)'Medium pulse out at ',tpeak/(beta*1e-6),' microseconds' -! write (*,*)'Vacuum pulse out at ',tpeak_vac/(beta*1e-6),' microseconds' - write (*,*)Omcpeak,(tpeak-tpeak_vac)/(beta*1e-9) -120 format(1x,f12.6,',',f12.6,',',f12.6,',',f12.6,',',f12.6,',',f12.6,',',f12.6) - t_end=secnds(0.E0) - t_elapsed=t_end-t_start - write(*,*)'T elapsed = ',t_elapsed - stop - end - -! - subroutine IMatConstruct(Om,Imat) -! -! NOTE: This subroutine actually calculates -Imat because we need to solve LPsi=-Imat. - implicit none - integer i,nmat - parameter (nmat=3) - real*8 Om - complex*8 Imat(nmat) - - Imat(1)=-cmplx(0.,-0.5*Om) - Imat(2)=-cmplx(0.,0.5*Om) - Imat(3)=cmplx(0.,0.) - return - end - - -! - subroutine LMatConstruct(Om,del,L) -! - - implicit none - integer nmat - parameter (nmat=3) - real*8 del,ga12,Om,W21 - complex*8 L(nmat,nmat) - common/para/ga12,W21 - - - L(1,1)=cmplx (-ga12,-del) - L(1,3)=cmplx(0.,Om) - L(2,2)=cmplx(-ga12,del) - L(2,3)=cmplx(0.,-Om) - L(3,1)=cmplx(0.,0.5*Om) - L(3,2)=cmplx(0.,-0.5*Om) - L(3,3)=cmplx(-W21,0.) - return - end - -! - subroutine MatCheck(L,Linv) -! - - implicit none - integer i,j,k,nmat - parameter (nmat=3) - complex*8 L(nmat,nmat),Linv(nmat,nmat),Res(nmat,nmat) - - write (*,*)'L = ' - do 10 i=1,nmat - write (*,120)(L(i,j),j=1,nmat) -10 continue - write (*,*)'Linv = ' - do 20 i=1,nmat - write (*,120)(Linv(i,j),j=1,nmat) -20 continue - write (*,*)'Res = ' - do 50 i=1,nmat - do 40 j=1,nmat - Res(i,j)=cmplx(0.,0.) - do 30 k=1,nmat - Res(i,j)=Res(i,j)+Linv(i,k)*L(k,j) -30 continue -40 continue - write (*,120)(Res(i,j),j=1,nmat) -50 continue -120 format(1x,3(f8.4,'+i',f8.4,' ')) - - return - end - -! - subroutine plotit(x,y,nmat,npts) -! -! See MATLAB routine that will do this plotting. - - implicit none - integer i,n,nmat,npts - real*8 x(npts) - complex*8 y(nmat,npts) - - write (*,*)'For now, we are just going to write the file' - open(9, FILE='TwoLevelPulseProp.txt') - do 10 n=1,npts - write (9,100)x(n),(y(i,n),i=1,3) -! write(*,100)x(n),(y(i,n),i=1,3) -10 continue -100 format(1x,f9.6,',',3(f9.6,',',f9.6,',')) - - - return - end - - - - - -!******************************************************** -! Numerical Recipes -!******************************************************** - - - SUBROUTINE ludcmp(a,n,np,indx,d) - implicit none - - INTEGER n,np,indx(n),NMAX - REAL*8 d,TINY - PARAMETER (NMAX=500,TINY=1.0e-20) - INTEGER i,imax,j,k - REAL aamax,dum,vv(NMAX) -! My changed variables - complex*8 sum,dum2 - complex*8 a(np,np) - - - d=1. - do 12 i=1,n - aamax=0. - do 11 j=1,n - if (cabs(a(i,j)).gt.aamax) aamax=cabs(a(i,j)) -11 continue - if (aamax.eq.0.) write (*,*) 'singular matrix in ludcmp' - vv(i)=1./aamax -12 continue - do 19 j=1,n - do 14 i=1,j-1 - sum=a(i,j) - do 13 k=1,i-1 - sum=sum-a(i,k)*a(k,j) -13 continue - a(i,j)=sum -14 continue - aamax=0. - do 16 i=j,n - sum=a(i,j) - do 15 k=1,j-1 - sum=sum-a(i,k)*a(k,j) -15 continue - a(i,j)=sum - dum=vv(i)*cabs(sum) - if (dum.ge.aamax) then - imax=i - aamax=dum - endif -16 continue - if (j.ne.imax)then - do 17 k=1,n - dum2=a(imax,k) - a(imax,k)=a(j,k) - a(j,k)=dum2 -17 continue - d=-d - vv(imax)=vv(j) - endif - indx(j)=imax - if (cabs(a(j,j)).eq.0.) a(j,j)=cmplx(TINY,TINY) - if(j.ne.n)then - dum2=1./a(j,j) - do 18 i=j+1,n - a(i,j)=a(i,j)*dum2 -18 continue - endif -19 continue - return - END - - SUBROUTINE lubksb(a,n,np,indx,b) - implicit none - INTEGER n,np,indx(n) - INTEGER i,ii,j,ll - -! My changed variables - complex*8 sum - complex*8 b(n) - complex*8 a(np,np) - ii=0 - do 12 i=1,n - ll=indx(i) - sum=b(ll) - b(ll)=b(i) - if (ii.ne.0)then - do 11 j=ii,i-1 - sum=sum-a(i,j)*b(j) -11 continue - else if (sum.ne.0.) then - ii=i - endif - b(i)=sum -12 continue - do 14 i=n,1,-1 - sum=b(i) - do 13 j=i+1,n - sum=sum-a(i,j)*b(j) -13 continue - b(i)=sum/a(i,i) -14 continue - return - END - diff --git a/fortran/navy_four_levels/Makefile b/fortran/navy_four_levels/Makefile index b55915f..f9cf9d2 100644 --- a/fortran/navy_four_levels/Makefile +++ b/fortran/navy_four_levels/Makefile @@ -14,8 +14,8 @@ $(runfiles): %.run : %.f95 clean: rm -f $(runfiles) - rm -f Movie4level_EndPoints_v*.dat - rm -f Movie4level_v*.dat - rm -f MovieParameters4level_v*.txt + rm -f Movie4level_EndPoints*.dat + rm -f Movie4level*.dat + rm -f MovieParameters4level*.txt |