C ****f* dataking/main[1.0] C C NAME C DATAKING C C AUTHOR C Rudolph Magyar C C SYNOPSIS C Process Datas. C Returns the Integral and Derive. C C FUNCTION C Integrate and Derive Data from C An Arbitrary Mesh C C INPUTS C file: IN.DAT C r(i), f(r) C C RESULT C screen: diagnostics, integral value of all data C file: OUT.DAT C r(i), f'(r), f(r), int_r(1)^r(i) f(r) dr(i) C C NOTES C Absolutely No Guaranties C C BUGS C C None Known. C *** C-------------------------------------- C Declare and Define C-------------------------------------- implicit none C define max. no. of lam. steps (used to dim arrays) integer msh parameter (msh=10000) C Counters integer i,j C number of points to read in and spline. integer nmax,halfmax C useful for the uniform mesh real*8 rmin,rmax real*8 range,step C Arrays C the real radial value of the point in question real*8 rin(msh),runif(msh) C Functions real*8 Fin(msh) real*8 f(msh),dfdr(msh),fint(msh) C used by the spline routine real*8 coeffs(msh) real*8 temp C-------------------------------------- C Set some defaults C-------------------------------------- C-------------------------------------- C Get Input Data C-------------------------------------- C Note we read in the R and F(R)'s Open (unit=1,file='in.dat') do i =1,msh read(1,*,err=10,end=10) rin(i),fin(i) nmax=i enddo 10 continue rmin=rin(1) rmax=rin(nmax) C Switch order in the file if not from C smallest r to largest if (rmin.gt.rmax) then C note the clever use of integer division halfmax=nmax/2 do i=1,halfmax temp=rin(i) rin(i)=rin(nmax+1-i) rin(nmax+1-i)=temp temp=fin(i) fin(i)=fin(nmax+1-i) fin(nmax+1-i)=temp enddo temp=rmin rmin=rmax rmax=temp endif write(6,*) 'rmin = ',rmin write(6,*) 'rmax = ',rmax,' at ',nmax do i=1,nmax C rin(i)=Log(rin(i)) fin(i)=Log(fin(i)) enddo C-------------------------------------- C Set up the unif mesh C-------------------------------------- range=rmax-rmin step=range/real(nmax-1) do i=1,nmax runif(i)=rmin+(i-1)*step enddo C do i=1,nmax C runif(i)=Log(runif(i)) C enddo C-------------------------------------- C spline input onto new uniform grid C-------------------------------------- C Notice we loop over each input function seperately C get spline coeffs call dspline(rin,fin,nmax,1d30,1d30,coeffs) c do spline do i=1,nmax call dsplint(rin,fin,coeffs,nmax,runif(i),f(i)) enddo C undo the logging do i=1,nmax f(i)=Exp(f(i)) enddo C undo logging of r coordinate C do i=1,nmax C rout(i)=Exp(rout(i)) C enddo C-------------------------------------- C get deriv and integral of the data C-------------------------------------- C note the n-1 to agree with the subroutines indexing convetion call sntgr2(nmax-1,f,fint) do i=1,nmax fint(i)=fint(i)*range enddo C note the n-1 to agree with the subroutines indexing convetion call stddff(nmax-1,f,dfdr) do i=1,nmax dfdr(i)=dfdr(i)/range enddo C-------------------------------------- C output results C-------------------------------------- C output to screen. write(6,*) 'Integral = ', fint(nmax)-fint(1) C output to file. Open (unit=2,file='out.dat') C FILE: FILEOUT.DAT C Data: do i=1,nmax write(2,9000) runif(i),dfdr(i),f(i),fint(i) enddo 9000 FORMAT(10E18.10) C That's it. We have reached the end. write(6,*) 'DONE.' 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)) 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)+ 1 ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6d0 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 double precision(a-h,o-z) dimension f(0:n),fint(0:n) dy=1.d0/dfloat(24*n) fint(0)=0.d0 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---------------------------------------------------------------------- c---------------------------------------------------------------------- c STDDFF simply differentiates a function from 0 to 1 on n+1 uniformly c spaced points, using 6-point formulas c Heiko Appel, 5.18.99 c---------------------------------------------------------------------- subroutine stddff(n,f,df) implicit real*8(a-h,o-z) dimension f(0:n),df(0:n),iwt0(0:5),iwt1(0:5),iwt2(0:5) data iwt0/-274,600,-600,400,-150,24/ data iwt1/-24,-130,240,-120,40,-6/ data iwt2/6,-60,-40,120,-30,4/ con=real(n)/120d0 df(0)=0d0 do j=0,5 df(0)=df(0)+f(j)*iwt0(j)*con enddo c df(0)=(-3*f(0)+4*f(1)-f(2))*real(n)/2d0 df(1)=0d0 do j=0,5 df(1)=df(1)+f(j)*iwt1(j)*con enddo c df(1)=(f(2)-f(0))*real(n)/2d0 do i=2,n-3 df(i)=0d0 do j=0,5 df(i)=df(i)+f(i-2+j)*iwt2(j)*con enddo enddo df(n-2)=0d0 do j=0,5 df(n-2)=df(n-2)-f(n-5+j)*iwt2(5-j)*con enddo df(n-1)=0d0 do j=0,5 df(n-1)=df(n-1)-f(n-5+j)*iwt1(5-j)*con enddo c df(n-1)=(f(n)-f(n-2))*real(n)/2d0 df(n)=0d0 do j=0,5 df(n)=df(n)-f(n-5+j)*iwt0(5-j)*con enddo c df(n)=(3*f(n)-4*f(n-1)+f(n-2))*real(n)/2d0 return end c-----------------------------------------------------------------------