program thomas implicit double precision(a-h,o-z) parameter(NDIM=1001) dimension u(NDIM),phi(NDIM),phip(NDIM) data acc/1.0d-10/ pi = acos(-1d0) write(6,1000) ' Enter z and n : ' 1000 format(/a,\$) read(5,*) z,n * xi = sqrt(9d0*pi**2/(128d0*z)) * correction by G. Faussurier (7/4/04) xi = (9d0*pi**2/(128d0*z))**(1d0/3d0) * xdphi = -(z-n)/z 100 continue write(6,1000) ' Enter two initial guesses for R : ' read(5,*) ra,rb a = ra/xi b = rb/xi i = 0 200 continue i = i + 1 if(i.eq.1) then x0 = a else if(i.eq.2) then x0 = b else if (fa*fb.gt.0.0d0) then write(6,*) ' No Solution with these guesses !' write(6,1010) ' Ra = ',xi*a,' phi(0,a)-1 = ',fa 1010 format(a,f10.6,a,f16.8) write(6,1010) ' Rb = ',xi*b,' phi(0,b)-1 = ',fb go to 100 else x0 = 0.5*(a+b) endif endif u(1) = sqrt(x0) h = -u(1)/(NDIM-1) do 210 j= 1, NDIM u(j) = u(1) + (j - 1) * h 210 continue phi(1) = 0 phip(1) = 2*xdphi/u(1) call runge(u,phi,phip,h) write(6,1020) ' Initial value for R =',xi*x0 1020 format(/a,f16.10) write(6,1030) ' phi(0) = ',phi(NDIM) 1030 format(a,f14.10/) if(i.eq.1) fa = phi(NDIM) - 1.0d0 if(i.eq.2) fb = phi(NDIM) - 1.0d0 if(i.gt.2) then fc = phi(NDIM) - 1.0d0 if (fc*fa.lt.0.0d0) then b = x0 fb = fc else a = x0 fa = fc endif endif if(abs(phi(NDIM) - 1.0d0).gt.acc) go to 200 ***** output file (1) * * thomas.dat: r[i], phi[i], N[i], Zeff[i] * * where V(r) == -Z/r + Vel(r) = -Zeff(r)/r * ***** rmax = xi*u(1)**2 open(unit=1,file='thomas.dat',form='formatted',status='unknown') rewind(1) write(1,1090) 1090 format(7x,'r',10x,'phi',10x,'N',10x,'Zeff') do 300 i = NDIM,1,-2 r = xi*u(i)**2 ph = phi(i) en = z*(0.5d0*u(i)*phip(i)-phi(i)+1) zeff = z*ph + (z-n)*r/rmax write(1,1100) r,ph,en,zeff 1100 format(4f12.5) 300 continue rewind(1) stop end double precision function f(x,y,yp) implicit double precision(a-h,o-z) parameter(NDIM=1001) if (x.eq.0.0d0) then f = 0.0d0 else f = 4.0 * x * sqrt(abs(y*y*y)) + yp / x endif return end subroutine runge(x,y,yp,h) implicit double precision(a-h,o-z) double precision k1,k2,k3,k4 ****************************************************************************** * Fourth-Order Runge-Kutta Method * Abramowitz and Stegun * Handbook of Math Functions (1964 edition) * NBS AMS 55 * Chapter 25, Equation 25.5.20 ****************************************************************************** parameter(NDIM=1001) dimension x(NDIM),y(NDIM),yp(NDIM) data six/6.0d0/,two/2.0d0/,half/0.5d0/,eighth/0.125d0/ do 100 i = 1,NDIM-2 u = x(i) v = y(i) w = yp(i) w = yp(i) k1 = h * f(u,v,w) u = x(i) + half * h v = y(i) + h * (half * yp(i) + eighth * k1) w = yp(i) + half * k1 k2 = h * f(u,v,w) u = x(i) + half * h v = y(i) + h * (half * yp(i) + eighth * k2) w = yp(i) + half * k2 k3 = h * f(u,v,w) u = x(i) + h v = y(i) + h * ( yp(i) + half * k3 ) w = yp(i) +k3 k4 = h * f(u,v,w) y(i+1) = y(i) + h * ( yp(i) + ( k1 + k2 + k3 )/six ) yp(i+1) = yp(i) + (k1 + two * k2 + two * k3 + k4)/six 100 continue * * Special treatnent for x = 0.0 * Adams-Bashforth Fourth-Order * Chapter 25, Equation 25.5.4 * i = NDIM-1 yp(i+1) = yp(i) + & h * ( 55.0d0 * f(x(i),y(i),yp(i)) & -59.0d0 * f(x(i-1),y(i-1),yp(i-1)) & +37.0d0 * f(x(i-2),y(i-2),yp(i-2)) & - 9.0d0 * f(x(i-3),y(i-3),yp(i-3)) ) /24.0d0 y(i+1) = y(i) + & h * ( 55.0d0 * yp(i) & -59.0d0 * yp(i-1) & +37.0d0 * yp(i-2) & - 9.0d0 * yp(i-3) ) /24.0d0 return end