diff options
Diffstat (limited to 'fortran/navy_four_levels/FourLevelPulseProp_Double.f95')
| -rwxr-xr-x | fortran/navy_four_levels/FourLevelPulseProp_Double.f95 | 708 |
1 files changed, 708 insertions, 0 deletions
diff --git a/fortran/navy_four_levels/FourLevelPulseProp_Double.f95 b/fortran/navy_four_levels/FourLevelPulseProp_Double.f95 new file mode 100755 index 0000000..06f7353 --- /dev/null +++ b/fortran/navy_four_levels/FourLevelPulseProp_Double.f95 @@ -0,0 +1,708 @@ + program FourLevelPulseProp_v3_Double +! +! Written by: Dr. Frank A. Narducci +! Written on: May 12, 2008 +! +! This program calculates the propagation of a pulse of arbitrary strength +! through a two level medium. The equations used are the full equations +! based on the Risken-Numedal discretization technique. +! +! This program only "watches" the evolution of the pulse in the cell. This is +! due to the constraint that the cell is very small relative to the pulse lengths +! that we want to use. If we watched the pulse outside the cell and then increased the +! resolution within the cell, the increased burden outside the cell because huge. +! +! v2 Notes: This program is based on the dimensionless equatiosn derived on 5/16 + implicit none +! +! Double Notes: This program is the same as TwoLevelPulseProp_v2 but with double precision + complex ci +! +! ThreeLevel Notes: This program is the same as TwoLevelPulseProp_v2_Double but now for a +! three level system +! +! Four Level Notes: Valid to first order in dt +! +! Version 3: Make a step function in the coupling field. + + character*150 fname + integer nmat,npts,Nfrac,Nframe,Nframemax,NSkip,NWrite,tpts,zpts + parameter (nmat=3,npts=100) !matrix size, number of detuning points in dispersion curve + !REMEMBER TO CHANGE NMAT IN LMatConstruct Routine + parameter (tpts=10,zpts=tpts+1) !Caution: funny things happened when tpts=200 (and presumably greater) + !tpts is the number of temporal points in the cell + parameter (Nframemax=2000000) + parameter ( NWrite=100) !number of frames to actually write + integer i,j,k,m,n + complex*16 a1,a2,a3,a4,a5,a6 + complex*16 b1,b2,b3,b4,b5,b6,b7 + complex*16 c1,c2,c3,c4,c5,c6 + complex*16 d1,d2,d3,d4,d5,d6,d7 + complex*16 e1,e2,e3,e4,e5 + complex*16 f1,f2,f3,f4,f5 + complex*16 g1,g2,g3,g4,g5,g6,g7 + complex*16 h1,h2,h3,h4,h5 + complex*16 i1,i2,i3,i4,i5,i6,i7 + complex*16 j1,j2,j3,j4,j5,j6,j7,j8 + complex*16 k1,k2,k3,k4,k5,k6,k7,k8,k9,k10 + complex*16 l1,l2,l3,l4,l5 + + complex *16 Omold, Omold_vac + real*8 alpha1,alpha2,alpha1tilde,alpha2tilde,alphac,alphactilde,beta,c,delmax,del1_prop,del2_prop,delc_prop + real*8 dt,dz,eta + real*8 W12,W21,W31,W32,W41,W42,W43,W34,ga12,ga13,ga14,ga23,ga24,ga34 + real*8 Ga2,Ga4,Om_crit + real*8 Lcell,Om1peak,Om2peak,Omcpeak,pi,tmax,tp,tshift,t_end,t_start,t_elapsed + real*8 tpeak,tpeak_vac + real*8 epsil,hbar,lambda + real*8 del(npts) + real*8 t,z(zpts) + complex*16 yplot(nmat,npts) + complex*16 Imat(nmat) + + + complex*16 Om1(zpts),Om2(zpts),Omc(zpts),Om_vac(zpts) + complex*16 rho11(zpts),rho12(zpts),rho13(zpts),rho14(zpts),rho21(zpts),rho22(zpts),rho23(zpts),rho24(zpts) + complex*16 rho31(zpts),rho32(zpts),rho33(zpts),rho34(zpts),rho41(zpts),rho42(zpts),rho43(zpts),rho44(zpts) + complex*16 rho11_last(zpts),rho12_last(zpts),rho13_last(zpts),rho14_last(zpts) + complex*16 rho21_last(zpts),rho22_last(zpts),rho23_last(zpts),rho24_last(zpts) + complex*16 rho31_last(zpts),rho32_last(zpts),rho33_last(zpts),rho34_last(zpts) + complex*16 rho41_last(zpts),rho42_last(zpts),rho43_last(zpts),rho44_last(zpts) + + !No Om_last because we never need the previous spatial point + complex*16 L(nmat,nmat),Linv(nmat,nmat),Ltemp(nmat,nmat) + + common/para/ga12,W21 + + + real*8 d !used by NR Routines + integer indx(nmat) + +! +! Fundamental numbers +! + ci=cmplx(0.,1.) + pi=acos(-1.0) + c=3e8 + hbar=1.055e-34 + epsil=8.85e-12 + +! +! Atomic numbers (based on Rubidium 85) +! + beta=2*pi*3e6 !in Hz + W41=0 + W42=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.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 + 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) t,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,',',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 + |