C ****f* shootaway/main[0.1] C C NAME C ShootAway C C AUTHOR C Rudolph Magyar C C SYNOPSIS C C FUNCTION C C INPUTS C C RESULT C C NOTES C C BUGS C C None Known. C *** C-------------------------------------- C Declare and Define C-------------------------------------- implicit none real*8 pi parameter(pi=3.14159265358979323) C define max. mesh size (used to dim arrays) integer msh parameter(msh=40000) C Counters integer i,j,nit,nitfin C The real mesh C Variables for the splined density with which we work C r(msh) the uniform r grid C dr the spacing between grid points C nmax the # of grid points used in calculations C rmax the radius at nmax real*8 r(0:msh) real*8 dr,rmax,rmin integer nmax,nmin C The wavefunction and auxillory functions C u= wf ' / wf real*8 wf(0:msh),dwfdr(0:msh) real*8 u(0:msh),dudr(0:msh) C The potential term real*8 v C The potential type label integer vtyp c energy, the ration of wf ' / wf at r=0 , the ext pot strength real*8 energy,lambda,Z C used in the shooting method real*8 emin,emax C convergence criteria real*8 delta,econv C-------------------------------------- C Set Some Defaults C-------------------------------------- write(6,*) write(6,*) & 'CODE WHICH USES THE SHOOTING METHOD TO SOLVE 1 D PROBLEMS' write(6,*) 'by Rudolph Magyar 2002 All Rights Reserved' write(6,*) delta=1d-10 econv=1d-10 nitfin=100 C-------------------------------------- C Set R Mesh Up C-------------------------------------- nmin=0 rmin=0.d0 nmax=40000 rmax=20.d0 C set up uniform r grid. dr = rmax/(real(nmax)) do i =0,nmax r(i)=real(i)*dr enddo write(6,*) 'rmin = ',r(0) write(6,*) 'rmax = ',r(nmax) write(6,*) 'npts = ',nmax C-------------------------------------- C Set the Hamiltonian Up C-------------------------------------- C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C%%%%%% The Poor Man's Hookes Hamiltonian lambda = 0d0 z = 1.d0 vtyp=1 C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C-------------------------------------- C shoot C-------------------------------------- C integrate upward first. C intial conditions u(0)=lambda C Guess eigenvalue Emax= 10000.d0 Emin= - 10000.d0 C main interation loop do nit = 1, nitfin energy=1/2d0*(emin+emax) write(6,*) 'it. number :',nit write(6,*) & ' energy guess = ',energy,' minmax =',emin,',',emax C loop over r values to solve for u(r) do i=0,nmax C calculate the first derivative at i dudr(i)= - 2 * Energy + 2*v(z,r(i),vtyp) - u(i) ** 2 C The following line is a simple first order integration C method and can be used to check that things are roughly C working. Be sure to comment out the RK method if using this. u(i+1)=u(i)+dudr(i)*dr C Call RK method. I haven't implimented this yet. C Using the info at i we get u(i+1). C The routine needs to be adopted to the tasks at hand. C call RK4F(u(i),dudr(i),r(i),dr,u(i+1),energy,z,vtyp) C Get rid of nasty numbers if(u(i+1).ne.u(i+1)) u(i+1)=0d0 C if(u(i+1).gt. 1d3) u(i+1)=1d3 C if(u(i+1).lt.-1d3) u(i+1)=-1d3 C Check to see if the trial k crosses the y axis. C If so then energy is too negative; emin=energy if ( & ( (u(i+1).gt.0d0) .and. (u(i).lt.0d0) ) & .or. & ( (u(i+1).lt.0d0).and.(u(i).gt.0d0) ) & ) & then emin=energy C write(6,*) 'A' goto 10 endif C end of r loop for u(r) enddo C Second loop to get the wavefunction. C the following is an arbitrary choice C However, since we start with wf > 0 C and we demand no nodes then we must C not allow the wf to cross the y axis. wf(0)=1.d0 do j=0,nmax C calculate the first derivative at j dwfdr(j)= wf(j)*u(j) C The following line is a simple first order integration C method and can be used to check that things are roughly C working. Be sure to comment out the RK method if using this. wf(j+1)=wf(j)+dwfdr(j)*dr C Get rid of nasty numbers if(wf(j+1).ne.wf(j+1)) wf(j+1)=0d0 if((wf(j+1).lt. 1d-30).and.(wf(j+1).gt.0d0)) wf(j+1)=0d0 C Enforce vanishing WF if the energy is converged. C In theory this should not be necessary, but due to C finite numeric this must be done. if( ( (Abs(emin-emax)).lt.econv ).and.(wf(j+1).lt.0d0)) then wf(j+1)=0.d0 endif C check to make sure that wf doesn't change sign if (wf(j+1).lt.0d0) then emax=energy C write(6,*) 'B' goto 10 endif enddo C Check to see if the WF goes within the target window if( wf(nmax) .gt. delta) then emin=energy C write(6,*) 'C' goto 10 else if( wf(nmax) .lt. -delta) then emax=energy C write(6,*) 'D' goto 10 endif C escape the loop is we passed all criterea. C Oh, how I will f90 had while loops. goto 20 C the code is send here when it is ordered to try a new energy 10 continue C end of the shooting loop enddo 20 if(nit.ne.nitfin) then write(6,*) 'Converged within ' write(6,*) 'Delta E = ',emax-emin write(6,*) 'wf(rmax) = ',delta else write(6,*) 'After ',nit,' interations' write(6,*) 'we have Delta E = ',emax-emin write(6,*) 'we have wf(rmax) = ',wf(nmax) endif write(6,*) 'Energy : ',energy C-------------------------------------- C Normalize the WF C-------------------------------------- C call integrate(wf,wfi,0,nmax,dr) C-------------------------------------- C output raw results C-------------------------------------- write(6,*) '' write(6,*) 'OUTPUT' write(6,*) '' C output to screen. write(6,*) 'ENERGY =',energy C output to file. Open (unit=2,file='wf.dat') C FILE: WF.DAT do i=0,nmax write(2,9000) i,r(i),wf(i) enddo 9000 FORMAT(I8,10E18.10) C That's it. We have reached the end. end C##################################################### C Define the potential functions C##################################################### C The potential function v(k,x,n) implicit none C v(x) ext = C 1/2*k*x^2 C the potential strength and x coordinate real*8 k,x,v C the potential type label integer n C Hooke's potential if(n.eq.1) v=1/2d0*k*x**2 C The airy potential if(n.eq.2) v=k*abs(x) RETURN END c----------------------------------------------------------------------- subroutine DSPLINE(X,Y,N,YP1,YPN,Y2) c is a copy of Numerical recipes spline routine, in real*8 c It sets up coefficients for a spline, in array Y2. It need only c be called once. (Note: yp1 and ypn often set to 1d30). c Kieron Burke, 9.15.99 c in: n=no of points c x=vector of x-values c y=vector of y-values c yp1=derivative at x1 c ypn=derivative at xn c out: y2=vector of ceofficients implicit real*8(a-h,o-z) PARAMETER (NMAX=10000) DIMENSION X(N),Y(N),Y2(N),U(NMAX) c endhead--------------------------------------------------------------- c 1 2 3 4 5 6 7 c 23456789012345678901234567890123456789012345678901234567890123456789012 IF (YP1.GT..99d30) THEN Y2(1)=0d0 U(1)=0d0 ELSE Y2(1)=-0.5d0 U(1)=(3d0/(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1) ENDIF DO 11 I=2,N-1 SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1)) P=SIG*Y2(I-1)+2d0 Y2(I)=(SIG-1.)/P U(I)=(6d0*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1)) * /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P 11 CONTINUE IF (YPN.GT..99d30) THEN QN=0. UN=0. ELSE QN=0.5d0 UN=(3d0/(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1))) ENDIF Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1d0) DO 12 K=N-1,1,-1 Y2(K)=Y2(K)*Y2(K+1)+U(K) 12 CONTINUE RETURN END c---------------------------------------------------------------------- subroutine DSPLINT(XA,YA,Y2A,N,X,Y) c is a copy of NUM REC routine splint, in real*8 c It gives the value of the splined function at point x. c dSPLINE must be called beforehand. c Kieron Burke, 9.15.99 c in: n=no of points c xa=vector of x-values c ya=vector of y-values c y2a=vector of coefs c x=chosen x c out: y=splined value for x implicit real*8(a-h,o-z) DIMENSION XA(N),YA(N),Y2A(N) c endhead--------------------------------------------------------------- c 1 2 3 4 5 6 7 c 23456789012345678901234567890123456789012345678901234567890123456789012 KLO=1 KHI=N 1 IF (KHI-KLO.GT.1) THEN K=(KHI+KLO)/2d0 IF(XA(K).GT.X)THEN KHI=K ELSE KLO=K ENDIF GOTO 1 ENDIF H=XA(KHI)-XA(KLO) IF (H.EQ.0.) PAUSE 'Bad XA input.' A=(XA(KHI)-X)/H B=(X-XA(KLO))/H Y=A*YA(KLO)+B*YA(KHI)+ * ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6d0 RETURN END c---------------------------------------------------------------------- SUBROUTINE RK4F(Y,DYDX,X,DX,YOUT,energy,z,vtyp) REAL*8 DX,X,DYDX,Y,YOUT real*8 energy,z integer vtyp real*8 k1,k2,k3,k4 real*8 v external v k1=dx*DYDX DYDX = & -2 * Energy +2*v(z,x+dx/2.d0,vtyp) **2 - (y+k1/2.d0)**2 k2=dx*DYDX DYDX = & -2 * Energy +2*v(z,x+dx/2.d0,vtyp) **2 - (y+k2/2.d0)**2 k3=dx*DYDX DYDX = & -2 * Energy +2*v(z,x+dx,vtyp) **2 - (y+k2)**2 k4=dx*DYDX yout=Y+k1/6.d0+k2/3.d0+k3/3.d0+k4/6.d0 C write(6,*) k1,k2,k3,k4,y,yout RETURN END c---------------------------------------------------------------------- subroutine integrate(dfdr,fun,A,B,dr) C Integrates dfdr from grid point A to B, C and returns Fun, the integral at each point C Using Simpson's 3/8 rule parameter (msh=32000) real*8 dfdr(0:msh) real*8 fun(0:msh) real*8 dr,x integer A,B,i Fun(A)=0.d0 Fun(A+1)=1/2.d0*(dfdr(A)+dfdr(A+1))*dr Fun(A+2)=1/3.d0*(dfdr(A)+4.d0*dfdr(A+1)+dfdr(A+2))*dr do i=A+3,B x=3/8.d0*(dfdr(i-3)+3.d0*dfdr(i-2)+3.d0*dfdr(i-1)+dfdr(i))*dr Fun(i)=x+Fun(i-3) enddo return end c 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012 c 1 2 3 4 5 6 7 c---------------------------------------------------------------------- c SNTGRL simply integrates a function from 0 to 1 on n+1 uniformly c spaced points; n must be at least 8. c The formula is from Numerical Recipes, 4.1.14. c Kieron Burke, 1.29.99 function sntgrl(n,f) implicit real*8(a-h,o-z) dimension f(0:n) sum=(-31d0*(f(0)+f(n))+11d0*(f(1)+f(n-1))-5d0*(f(2)+f(n-2)) & +f(3)+f(n-3))/48d0 do i=0,n sum=sum+f(i) enddo sntgrl=sum/real(n) return end c---------------------------------------------------------------------- c---------------------------------------------------------------------- c SNTGR2 simply integrates a function from 0 to 1 on n+1 uniformly c spaced points, producing the integral at all points inbetween. c The formula is from Numerical Recipes, 4.1.14. c Kieron Burke, 1.29.99 subroutine sntgr2(n,f,fint) implicit real*8(a-h,o-z) real*8 f(0:n),fint(0:n) dy=1d0/dfloat(24*n) fint(0)=0d0 fint(1)=dy*(f(3)-5*f(2)+19*f(1)+9*f(0)) do i=2,n-1 fint(i)=fint(i-1)+dy*(13*(f(i)+f(i-1))-(f(i+1)+f(i-2))) enddo fint(n)=fint(n-1)+dy*(9*f(n)+19*f(n-1)-5*f(n-2)+f(n-3)) return end c----------------------------------------------------------------------