! ********************************************************* FUNCTION bessi(n,x) implicit double precision (a-h,o-z) INTEGER n,IACC PARAMETER (IACC=40,BIGNO=1.0e10,BIGNI=1.0e-10) CU USES bessi0 INTEGER j,m ! if (n.lt.2) stop 'bad argument n in bessi' ! --------------------------------------------------------- ! RCF modified 31 Oct 00; G. Penn 30 Jun 01 ! if (x.eq.0.) then ! if( ABS(x) .lt. 1.e-10 ) then if( (ABS(x)/2.)**n .lt. 1.e-20 ) then ! --------------------------------------------------------- bessi=0. else tox=2.0/abs(x) bip=0.0 bi=1.0 bessi=0. m=2*((n+int(sqrt(float(IACC*n))))) do 11 j=m,1,-1 bim=bip+float(j)*tox*bi bip=bi bi=bim if (abs(bi).gt.BIGNO) then bessi=bessi*BIGNI bi=bi*BIGNI bip=bip*BIGNI endif if (j.eq.n) bessi=bip 11 continue bessi=bessi*bessi0(x)/bi if (x.lt.0..and.mod(n,2).eq.1) bessi=-bessi endif return END ! ********************************************************* FUNCTION bessi0(x) implicit double precision (a-h,o-z) common/FLAG/iflag ! RCF modified 10 MAR 01 SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,3.5156229d0,3.0899424d0, *1.2067492d0,0.2659732d0,0.360768d-1,0.45813d-2/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,0.1328592d-1, *0.225319d-2,-0.157565d-2,0.916281d-2,-0.2057706d-1,0.2635537d-1, *-0.1647633d-1,0.392377d-2/ ! --------------------------------------------------------- ! RCF modified 31 Oct 00; 10 Mar 01 if( ABS(x) .lt. 1.e-10 ) then bessi0 = 1. return else if( ABS(x) .gt. 100. ) then bessi0 = 0. iflag = -74 return end if ! --------------------------------------------------------- if (abs(x).lt.3.75) then y=(x/3.75)**2 bessi0=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))) else ax=abs(x) y=3.75/ax bessi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y* *(q7+y*(q8+y*q9)))))))) endif return END c ******************************************************** FUNCTION bessi1(x) implicit double precision (a-h,o-z) common/FLAG/iflag ! RCF modified 10 Mar 01 SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0, *0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,-0.3988024d-1, *-0.362018d-2,0.163801d-2,-0.1031555d-1,0.2282967d-1,-0.2895312d-1, *0.1787654d-1,-0.420059d-2/ ! --------------------------------------------------------- ! RCF modified 31 Oct 00; 10 Mar 01; 26 Apr 01 if( ABS(x) .lt. 1.e-10 ) then bessi1 = 0. return else if( ABS(x) .gt. 100. ) then bessi1 = 0. iflag = -74 return end if ! --------------------------------------------------------- if (abs(x).lt.3.75) then y=(x/3.75)**2 bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))) else ax=abs(x) y=3.75/ax bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y* *(q7+y*(q8+y*q9)))))))) if(x.lt.0.)bessi1=-bessi1 endif return END c ******************************************************** FUNCTION bessj0(x) implicit double precision (a-h,o-z) SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4, *s5,s6 DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4, *-.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-1, *.1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/ DATA r1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0, *651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/,s1,s2, *s3,s4,s5,s6/57568490411.d0,1029532985.d0,9494680.718d0, *59272.64853d0,267.8532712d0,1.d0/ ! if(abs(x).lt.8.)then y=x**2 bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+y* *(s4+y*(s5+y*s6))))) else ax=abs(x) z=8./ax y=z**2 xx=ax-.785398164 bessj0=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y* *p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5))))) endif return END c ******************************************************** FUNCTION bessj1(x) implicit double precision (a-h,o-z) SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4, *s5,s6 DATA r1,r2,r3,r4,r5,r6/72362614232.d0,-7895059235.d0, *242396853.1d0,-2972611.439d0,15704.48260d0,-30.16036606d0/,s1,s2, *s3,s4,s5,s6/144725228442.d0,2300535178.d0,18583304.74d0, *99447.43394d0,376.9991397d0,1.d0/ DATA p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4, *.2457520174d-5,-.240337019d-6/, q1,q2,q3,q4,q5/.04687499995d0, *-.2002690873d-3,.8449199096d-5,-.88228987d-6,.105787412d-6/ ! if(abs(x).lt.8.)then y=x**2 bessj1=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+ *y*(s4+y*(s5+y*s6))))) else ax=abs(x) z=8./ax y=z**2 xx=ax-2.356194491 bessj1=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y* *p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))*sign(1.d0,x) endif return END c ******************************************************** FUNCTION bessk(n,x) implicit double precision (a-h,o-z) INTEGER n CU USES bessk0,bessk1 INTEGER j ! if (n.lt.2) stop 'bad argument n in bessk' tox=2.0/x bkm=bessk0(x) bk=bessk1(x) do 11 j=1,n-1 bkp=bkm+j*tox*bk bkm=bk bk=bkp 11 continue bessk=bk return END c ******************************************************** FUNCTION bessk0(x) implicit double precision (a-h,o-z) common/FLAG/iflag ! RCF modified 10 MAR 01 CU USES bessi0 DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7 DATA p1,p2,p3,p4,p5,p6,p7/-0.57721566d0,0.42278420d0,0.23069756d0, *0.3488590d-1,0.262698d-2,0.10750d-3,0.74d-5/ DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,-0.7832358d-1,0.2189568d-1, *-0.1062446d-1,0.587872d-2,-0.251540d-2,0.53208d-3/ ! --------------------------------------------------------- ! RCF modified 10 Mar 01 if( x .lt. 1.e-3 ) then bessk0 = 0. iflag = -75 return else if( x .gt. 100. ) then bessk0 = 0. return end if ! --------------------------------------------------------- if (x.le.2.0) then y=x*x/4.0 bessk0=(-log(x/2.0)*bessi0(x))+(p1+y*(p2+y*(p3+y*(p4+y*(p5+y* *(p6+y*p7)))))) else y=(2.0/x) bessk0=(exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y* *q7)))))) endif return END c ******************************************************** FUNCTION bessk1(x) implicit double precision (a-h,o-z) common/FLAG/iflag ! RCF modified 10 MAR 01 CU USES bessi1 DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7 DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,0.15443144d0,-0.67278579d0, *-0.18156897d0,-0.1919402d-1,-0.110404d-2,-0.4686d-4/ DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,0.23498619d0,-0.3655620d-1, *0.1504268d-1,-0.780353d-2,0.325614d-2,-0.68245d-3/ ! --------------------------------------------------------- ! RCF modified 10 Mar 01 if( x .lt. 1.e-3 ) then bessk1 = 0. iflag = -75 return else if( x .gt. 100. ) then bessk1 = 0. return end if ! --------------------------------------------------------- if (x.le.2.0) then y=x*x/4.0 bessk1=(log(x/2.0)*bessi1(x))+(1.0/x)*(p1+y*(p2+y*(p3+y*(p4+y* *(p5+y*(p6+y*p7)))))) else y=2.0/x bessk1=(exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y* *q7)))))) endif return end c ******************************************************** FUNCTION bessy0(x) implicit double precision (a-h,o-z) CU USES bessj0 SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4, *s5,s6 DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4, *-.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-1, *.1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/ DATA r1,r2,r3,r4,r5,r6/-2957821389.d0,7062834065.d0, *-512359803.6d0,10879881.29d0,-86327.92757d0,228.4622733d0/,s1,s2, *s3,s4,s5,s6/40076544269.d0,745249964.8d0,7189466.438d0, *47447.26470d0,226.1030244d0,1.d0/ if(x.lt.8.)then y=x**2 bessy0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+y* *(s4+y*(s5+y*s6)))))+.636619772*bessj0(x)*log(x) else z=8./x y=z**2 xx=x-.785398164 bessy0=sqrt(.636619772/x)*(sin(xx)*(p1+y*(p2+y*(p3+y*(p4+y* *p5))))+z*cos(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5))))) endif return END c ******************************************************** FUNCTION bessy1(x) implicit double precision (a-h,o-z) CU USES bessj1 SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4, *s5,s6,s7 DATA p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4, *.2457520174d-5,-.240337019d-6/, q1,q2,q3,q4,q5/.04687499995d0, *-.2002690873d-3,.8449199096d-5,-.88228987d-6,.105787412d-6/ DATA r1,r2,r3,r4,r5,r6/-.4900604943d13,.1275274390d13, *-.5153438139d11,.7349264551d9,-.4237922726d7,.8511937935d4/,s1,s2, *s3,s4,s5,s6,s7/.2499580570d14,.4244419664d12,.3733650367d10, *.2245904002d8,.1020426050d6,.3549632885d3,1.d0/ if(x.lt.8.)then y=x**2 bessy1=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+ *y*(s4+y*(s5+y*(s6+y*s7))))))+.636619772*(bessj1(x)*log(x)-1./x) else z=8./x y=z**2 xx=x-2.356194491 bessy1=sqrt(.636619772/x)*(sin(xx)*(p1+y*(p2+y*(p3+y*(p4+y* *p5))))+z*cos(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5))))) endif return END ! ********************************************************* C C THREE DIMENSIONAL BSPLINE ROUTINES C c Report all suspected bugs to J. Scott Berg c or c c $Date: 2003/03/28 17:51:07 $ c $Revision: 1.3 $ c $Source: /home/jsberg/cvsroot/bsplinef/bsplinef3.f,v $ c C THESE ARE LIKE THE ONE-DIMENSIONAL ROUTINES, EXCEPT THAT THREE C SETS OF ABCISSAS ARE SUPPLIED, NX, N, AND K ARE THREE-ELEMENT C ARRAYS, AND THE VALUES TO INTERPOLATE ARE SUPPLIED IN AN ARRAY C WITH THE FIRST VARIABLE VARYING MOST RAPIDLY. C C DOUBLE PRECISION FUNCTION BS3EV(X,C,N,K,T,S) C C INPUTS: C C X(3) (DP) = POINT AT WHICH TO EVALUATE SPLINE SERIES C C(N(1)*N(2)*N(3)) (DP) = SPLINE COEFFICIENTS C N(3) (I) = NUMBERS OF SPLINE BASIS FUNCTIONS C K(3) (I) = SPLINE ORDERS C T(N(1)+K(1)+N(2)+K(2)+N(3)+K(3)) (DP) = KNOT SEQUENCES C S(K(2)*K(3)+K(1)) (DP) = SCRATCH ARRAY C ! ********************************************************* DOUBLE PRECISION FUNCTION BS3EV(X,C,N,K,T,S) C C EVALUATE INTERPOLATION OVER 3 DIMENSIONS. C X IS POINT TO EVALUATE AT C C IS COEFFICIENT ARRAY C N ARE THE NUMBER OF COEFFICIENTS IN EACH DIMENSION C K IS BSPLINE ORDER C T ARE KNOT SEQUENCES, STACKED SEQUENTIALLY C S IS SCRATCH ARRAY, SIZE K(2)*K(3)+K(1) C IMPLICIT NONE INTEGER K(3),N(3) DOUBLE PRECISION C(*),S(*),T(*),X(3) INTEGER I1,I2,I3,J1,J2,J3,KK,NS,BSIVL,IC I1 = BSIVL(X(1),T(1),N(1),K(1)) I2 = BSIVL(X(2),T(N(1)+K(1)+1),N(2),K(2)) I3 = BSIVL(X(3),T(N(1)+K(1)+N(2)+K(2)+1),N(3),K(3)) NS = K(2)*K(3) DO 100 J3=1,K(3) DO 100 J2=1,K(2) IC=N(1)*(I2+J2-K(2)-1+N(2)*(I3+J3-K(3)-1)) DO 120 J1=1,K(1) 120 S(J1+NS)=C(I1+J1-K(1)+IC) DO 130 KK=K(1)-1,1,-1 DO 130 J1=I1,I1-KK+1,-1 130 S(J1-I1+K(1)+NS) = (S(J1-I1+K(1)+NS)*(X(1)-T(J1))+ $ S(J1-I1+K(1)-1+NS)*(T(J1+KK)-X(1)))/ $ (T(J1+KK)-T(J1)) 100 S(J2+K(2)*(J3-1)) = S(K(1)+NS) NS = N(1)+K(1) DO 200 J3=1,K(3) J1=K(2)*J3-I2 DO 200 KK=K(2)-1,1,-1 DO 200 J2=I2,I2-KK+1,-1 200 S(J2+J1) = (S(J2+J1)*(X(2)-T(J2+NS))+ $ S(J2+J1-1)*(T(J2+KK+NS)-X(2))) $ /(T(J2+KK+NS)-T(J2+NS)) NS = NS+N(2)+K(2) J2 = K(3)-I3 DO 300 KK=K(3)-1,1,-1 DO 300 J3=I3,I3-KK+1,-1 300 S((J3+J2)*K(2)) = (S((J3+J2)*K(2))*(X(3)-T(J3+NS))+ $ S((J3+J2-1)*K(2))*(T(J3+KK+NS)-X(3)))/ $ (T(J3+KK+NS)-T(J3+NS)) BS3EV = S(K(2)*K(3)) END ! ********************************************************* INTEGER FUNCTION BSIVL(X,T,N,K) C C RETURN I SUCH THAT T(I) <= X < T(I+1) C IMPLICIT NONE INTEGER K,N DOUBLE PRECISION X,T(N+K) INTEGER I,J,L I = K J = N+1 100 IF (I+1.EQ.J) GOTO 300 L = (I+J)/2 IF (X.LT.T(L)) GOTO 200 I = L GOTO 100 200 J = L GOTO 100 300 BSIVL = I END c ******************************************************** FUNCTION dfridr(func,x,h,err) implicit double precision (a-h,o-z) INTEGER NTAB PARAMETER (CON=1.4,CON2=CON*CON,BIG=1.E30,NTAB=10,SAFE=2.) EXTERNAL func CU USES func INTEGER i,j REAL*8 a(NTAB,NTAB) ! if(h.eq.0.) stop 'h must be nonzero in dfridr' hh=h a(1,1)=(func(x+hh)-func(x-hh))/(2.0*hh) err=BIG do 12 i=2,NTAB hh=hh/CON a(1,i)=(func(x+hh)-func(x-hh))/(2.0*hh) fac=CON2 do 11 j=2,i a(j,i)=(a(j-1,i)*fac-a(j-1,i-1))/(fac-1.) fac=CON2*fac errt=max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1))) if (errt.le.err) then err=errt dfridr=a(j,i) endif 11 continue if(abs(a(i,i)-a(i-1,i-1)).ge.SAFE*err)return 12 continue return END c ******************************************************** FUNCTION elle(phi,ak) implicit double precision (a-h,o-z) CU USES rd,rf ! s=sin(phi) cc=cos(phi)**2 q=(1.-s*ak)*(1.+s*ak) elle=s*(rf(cc,q,1.d0)-((s*ak)**2)*rd(cc,q,1.d0)/3.) return END c ******************************************************** FUNCTION ellf(phi,ak) implicit double precision (a-h,o-z) CU USES rf ! s=sin(phi) ellf=s*rf(cos(phi)**2,(1.-s*ak)*(1.+s*ak),1.d0) return END c ******************************************************** FUNCTION ellpi(phi,en,ak) implicit double precision (a-h,o-z) CU USES rf,rj ! s=sin(phi) enss=en*s*s cc=cos(phi)**2 q=(1.-s*ak)*(1.+s*ak) ellpi=s*(rf(cc,q,1.d0)-enss*rj(cc,q,1.d0,1.+enss)/3.) return END c ******************************************************** FUNCTION factrl(n) implicit double precision (a-h,o-z) INTEGER n CU USES gammln INTEGER j,ntop REAL*8 a(33) SAVE ntop,a DATA ntop,a(1)/0,1./ ! if (n.lt.0) then stop 'negative factorial in factrl' else if (n.le.ntop) then factrl=a(n+1) else if (n.le.32) then do 11 j=ntop+1,n a(j+1)=j*a(j) 11 continue ntop=n factrl=a(n+1) else factrl=exp(gammln(n+1.d0)) endif return END c ******************************************************** FUNCTION gammln(xx) implicit double precision (a-h,o-z) INTEGER j real*8 cof(6) SAVE cof,stp DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, *-.5395239384953d-5,2.5066282746310005d0/ ! x=xx y=x tmp=x+5.5d0 tmp=(x+0.5d0)*log(tmp)-tmp ser=1.000000000190015d0 do 11 j=1,6 y=y+1.d0 ser=ser+cof(j)/y 11 continue gammln=tmp+log(stp*ser/x) return END c ******************************************************** FUNCTION gasdev(idum) implicit double precision (a-h,o-z) INTEGER idum CU USES ran1 INTEGER iset SAVE iset,gset DATA iset/0/ ! if (iset.eq.0) then 1 v1=2.*ran1(idum)-1. v2=2.*ran1(idum)-1. rsq=v1**2+v2**2 if(rsq.ge.1..or.rsq.eq.0.)goto 1 fac=sqrt(-2.*log(rsq)/rsq) gset=v1*fac gasdev=v2*fac iset=1 else gasdev=gset iset=0 endif return END c ******************************************************** SUBROUTINE hunt(xx,n,x,jlo) implicit double precision (a-h,o-z) INTEGER jlo,n REAL*8 xx(n) INTEGER inc,jhi,jm LOGICAL ascnd ! ascnd=xx(n).gt.xx(1) if(jlo.le.0.or.jlo.gt.n)then jlo=0 jhi=n+1 goto 3 endif inc=1 if(x.ge.xx(jlo).eqv.ascnd)then 1 jhi=jlo+inc if(jhi.gt.n)then jhi=n+1 else if(x.ge.xx(jhi).eqv.ascnd)then jlo=jhi inc=inc+inc goto 1 endif else jhi=jlo 2 jlo=jhi-inc if(jlo.lt.1)then jlo=0 else if(x.lt.xx(jlo).eqv.ascnd)then jhi=jlo inc=inc+inc goto 2 endif endif 3 if(jhi-jlo.eq.1)return jm=(jhi+jlo)/2 if(x.gt.xx(jm).eqv.ascnd)then jlo=jm else jhi=jm endif goto 3 END c ******************************************************** integer function IMPMID(x0,v) C C Implicit midpoint integrator. Assumes step size 1. C Takes a 6-D vector x0 which is modified, and a subroutine v C The subroutine v is the vector field, and takes two arguments: C a 6-D vector where the vector field is evaluated, and a 6-D C vector in which the vector field is returned. C Returns 0 on success, 1 if we don't converge in maxit iterations C ! Ref. J.S. Berg 22 Nov 2004 ! double precision x0(6) double precision e0(6),e1(6),v0(6),v1(6),x1(6) logical b(6),ball integer i,j,maxit,rc,v parameter (maxit=100) external v rc = v(x0,v0) if (rc.ne.0) then impmid = rc return end if do 100 i=1,6 x1(i) = x0(i)+v0(i) b(i) = .false. 100 continue rc = v(x1,v1) if (rc.ne.0) then impmid = rc return end if do 200 i=1,6 e1(i) = abs(v0(i))+abs(v1(i)) v1(i) = 0.5d0*(v0(i)+v1(i)) 200 continue j = 0 300 do 400 i=1,6 x1(i) = x0(i)+0.5d0*v1(i) v0(i) = v1(i) e0(i) = e1(i) 400 continue rc = v(x1,v1) if (rc.ne.0) then impmid = rc return end if ball = .true. do 500 i=1,6 e1(i) = abs(v1(i)-v0(i)) b(i) = b(i) .or. e1(i).eq.0d0 .or. e1(i).ge.e0(i) ball = ball .and. b(i) 500 continue j = j + 1 if (.not.ball.and.j.lt.maxit) goto 300 if (j.lt.maxit) then impmid = 0 else impmid = -1 end if do 600 i=1,6 x0(i) = x0(i) + v1(i) 600 continue end ! ******************************************************** subroutine LOCATE2(xx,n,x,j) ! ! Locate position J of variable X in the array XX(1:N) ! Separate array into regions [1) [2) [3) ... [N) ! x=xx(1) => j=1 x=xx(n) => j=n ! Return j=0 if x is out of the array range ! implicit real*8(a-h,o-z) real*8 xx(n) ! if( x.lt.xx(1) .or. x.gt.xx(n) ) then j = 0 go to 900 end if ! jl = 0 ! lower bracket jh = n + 1 ! upper bracket ! 10 continue if( jh-jl .gt. 1 ) then ! continue working jm = (jl + jh) / 2 ! bisect the bracket if( x .ge. xx(jm) ) then jl = jm else jh = jm end if go to 10 end if ! j = jl ! 900 continue end c ******************************************************** subroutine LUDCMP(a,n,np,indx,d) ! ! Make LU decomposition of matrix ! implicit double precision (a-h,o-z) INTEGER n,np,indx(n),NMAX REAL*8 a(np,np) PARAMETER (NMAX=500,TINY=1.0e-20) INTEGER i,imax,j,k REAL*8 vv(NMAX) ! d=1. do 12 i=1,n aamax=0. do 11 j=1,n if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) 11 continue if (aamax.eq.0.) then ! singular matrix in ludcmp d = 0. go to 900 end if 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)*abs(sum) if (dum.ge.aamax) then imax=i aamax=dum endif 16 continue if (j.ne.imax)then do 17 k=1,n dum=a(imax,k) a(imax,k)=a(j,k) a(j,k)=dum 17 continue d=-d vv(imax)=vv(j) endif indx(j)=imax if(a(j,j).eq.0.)a(j,j)=TINY if(j.ne.n)then dum=1./a(j,j) do 18 i=j+1,n a(i,j)=a(i,j)*dum 18 continue endif 19 continue ! 900 continue return END c ******************************************************** FUNCTION poidev(xm,idum) implicit double precision (a-h,o-z) INTEGER idum PARAMETER (PI=3.141592654) CU USES gammln,ran1 SAVE alxm,g,oldm,sq DATA oldm /-1./ ! if (xm.lt.12.)then if (xm.ne.oldm) then oldm=xm g=exp(-xm) endif em=-1 t=1. 2 em=em+1. t=t*ran1(idum) if (t.gt.g) goto 2 else if (xm.ne.oldm) then oldm=xm sq=sqrt(2.*xm) alxm=log(xm) g=xm*alxm-gammln(xm+1.) endif 1 y=tan(PI*ran1(idum)) em=sq*y+xm if (em.lt.0.) goto 1 ! --------------------------------------------------------- ! RCF modified 28 Sep 99 ! prevent large value of em (e.g. 559) from causing ! a floating point underflow in the calculation of t, below if( em .gt. 30. ) go to 1 ! --------------------------------------------------------- em=int(em) t=0.9*(1.+y**2)*exp(em*alxm-gammln(em+1.)-g) if (ran1(idum).gt.t) goto 1 endif poidev=em return END c ******************************************************** SUBROUTINE polcof(xa,ya,n,cof) implicit double precision (a-h,o-z) INTEGER n,NMAX REAL*8 cof(n),xa(n),ya(n) PARAMETER (NMAX=15) CU USES polint INTEGER i,j,k REAL*8 x(NMAX),y(NMAX) ! do 11 j=1,n x(j)=xa(j) y(j)=ya(j) 11 continue do 14 j=1,n call polint(x,y,n+1-j,0.d0,cof(j),dy) xmin=1.e38 k=0 do 12 i=1,n+1-j if (abs(x(i)).lt.xmin)then xmin=abs(x(i)) k=i endif if(x(i).ne.0.)y(i)=(y(i)-cof(j))/x(i) 12 continue do 13 i=k+1,n+1-j y(i-1)=y(i) x(i-1)=x(i) 13 continue 14 continue return END c ******************************************************** SUBROUTINE polint(xa,ya,n,x,y,dy) implicit double precision (a-h,o-z) INTEGER n,NMAX REAL*8 xa(n),ya(n) PARAMETER (NMAX=10) INTEGER i,m,ns REAL*8 c(NMAX),d(NMAX) ! ns=1 dif=abs(x-xa(1)) do 11 i=1,n dift=abs(x-xa(i)) if (dift.lt.dif) then ns=i dif=dift endif c(i)=ya(i) d(i)=ya(i) 11 continue y=ya(ns) ns=ns-1 do 13 m=1,n-1 do 12 i=1,n-m ho=xa(i)-x hp=xa(i+m)-x w=c(i+1)-d(i) den=ho-hp if(den.eq.0.)stop 'failure in polint' den=w/den d(i)=hp*den c(i)=ho*den 12 continue if (2*ns.lt.n-m)then dy=c(ns+1) else dy=d(ns) ns=ns-1 endif y=y+dy 13 continue return END ! ********************************************************* SUBROUTINE qgaus(func,a,b,ss) implicit double precision (a-h,o-z) ! REAL a,b,ss,func EXTERNAL func INTEGER j ! REAL dx,xm,xr,w(5),x(5) dimension w(5),x(5) SAVE w,x DATA w/.2955242247,.2692667193,.2190863625,.1494513491, *.0666713443/ DATA x/.1488743389,.4333953941,.6794095682,.8650633666, *.9739065285/ xm=0.5*(b+a) xr=0.5*(b-a) ss=0 do 11 j=1,5 dx=xr*x(j) ss=ss+w(j)*(func(xm+dx)+func(xm-dx)) 11 continue ss=xr*ss return END ! ********************************************************* SUBROUTINE qromb(func,a,b,ss) implicit double precision (a-h,o-z) INTEGER JMAX,JMAXP,K,KM EXTERNAL func ! --------------------------------------------------------- ! RCF modified 5 Oct 00 ! PARAMETER (EPS=1.e-6, JMAX=20, JMAXP=JMAX+1, K=5, KM=K-1) PARAMETER (EPS=5.e-6, JMAX=20, JMAXP=JMAX+1, K=5, KM=K-1) ! --------------------------------------------------------- CU USES polint,trapzd INTEGER j REAL*8 h(JMAXP),s(JMAXP) ! h(1)=1. do 11 j=1,JMAX call trapzd(func,a,b,s(j),j) if (j.ge.K) then call polint(h(j-KM),s(j-KM),K,0.d0,ss,dss) if (abs(dss).le.EPS*abs(ss)) return ! --------------------------------------------------------- ! RCF modified 5 Oct 00 ! Doesn't converge for symmetric case with area ~0 if( ABS(ss) .le. eps ) return ! --------------------------------------------------------- endif s(j+1)=s(j) h(j+1)=0.25*h(j) 11 continue stop 'too many steps in qromb' END c ******************************************************** FUNCTION ran1(idum) implicit double precision (a-h,o-z) INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV PARAMETER (IA=16807,IM=2147483647,AM=1.d0/IM,IQ=127773,IR=2836, *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=2.3d-16,RNMX=1.d0-EPS) INTEGER j,k,iv(NTAB),iy SAVE iv,iy DATA iv /NTAB*0/, iy /0/ ! if (idum.le.0.or.iy.eq.0) then idum=max(-idum,1) do 11 j=NTAB+8,1,-1 k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM if (j.le.NTAB) iv(j)=idum 11 continue iy=iv(1) endif k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM j=1+iy/NDIV iy=iv(j) iv(j)=idum ran1=min(AM*iy,RNMX) return END c ******************************************************** FUNCTION rc(x,y) implicit double precision (a-h,o-z) PARAMETER (ERRTOL=0.0012,TINY=1.69e-38,SQRTNY=1.3e-19,BIG=3.E37, *TNBG=TINY*BIG,COMP1=2.236/SQRTNY,COMP2=TNBG*TNBG/25.,THIRD=1./3., *C1=.3,C2=1./7.,C3=.375,C4=9./22.) ! if(x.lt.0..or.y.eq.0..or.(x+abs(y)).lt.TINY.or.(x+ *abs(y)).gt.BIG.or.(y.lt.-COMP1.and.x.gt.0..and.x.lt.COMP2))stop *'invalid arguments in rc' if(y.gt.0.)then xt=x yt=y w=1. else xt=x-y yt=-y w=sqrt(x)/sqrt(xt) endif 1 continue alamb=2.*sqrt(xt)*sqrt(yt)+yt xt=.25*(xt+alamb) yt=.25*(yt+alamb) ave=THIRD*(xt+yt+yt) s=(yt-ave)/ave if(abs(s).gt.ERRTOL)goto 1 rc=w*(1.+s*s*(C1+s*(C2+s*(C3+s*C4))))/sqrt(ave) return END c ******************************************************** FUNCTION rd(x,y,z) implicit double precision (a-h,o-z) PARAMETER (ERRTOL=0.0015,TINY=1.e-25,BIG=4.5E21,C1=3./14., & C2=1./6.,C3=9./22.,C4=3./26.,C5=.25*C3,C6=1.5*C4) ! if(min(x,y).lt.0..or.min(x+y,z).lt.TINY.or.max(x,y, *z).gt.BIG)stop 'invalid arguments in rd' xt=x yt=y zt=z sum=0. fac=1. 1 continue sqrtx=sqrt(xt) sqrty=sqrt(yt) sqrtz=sqrt(zt) alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz sum=sum+fac/(sqrtz*(zt+alamb)) fac=.25*fac xt=.25*(xt+alamb) yt=.25*(yt+alamb) zt=.25*(zt+alamb) ave=.2*(xt+yt+3.*zt) delx=(ave-xt)/ave dely=(ave-yt)/ave delz=(ave-zt)/ave if(max(abs(delx),abs(dely),abs(delz)).gt.ERRTOL)goto 1 ea=delx*dely eb=delz*delz ec=ea-eb ed=ea-6.*eb ee=ed+ec+ec rd=3.*sum+fac*(1.+ed*(-C1+C5*ed-C6*delz*ee)+delz*(C2*ee+delz*(-C3* *ec+delz*C4*ea)))/(ave*sqrt(ave)) return END c ******************************************************** FUNCTION rf(x,y,z) implicit double precision (a-h,o-z) PARAMETER (ERRTOL=0.0025,TINY=1.5e-38,BIG=3.E37,THIRD=1./3., *C1=1./24.,C2=.1,C3=3./44.,C4=1./14.) ! if(min(x,y,z).lt.0..or.min(x+y,x+z,y+z).lt.TINY.or.max(x,y, *z).gt.BIG)stop 'invalid arguments in rf' xt=x yt=y zt=z 1 continue sqrtx=sqrt(xt) sqrty=sqrt(yt) sqrtz=sqrt(zt) alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz xt=.25*(xt+alamb) yt=.25*(yt+alamb) zt=.25*(zt+alamb) ave=THIRD*(xt+yt+zt) delx=(ave-xt)/ave dely=(ave-yt)/ave delz=(ave-zt)/ave if(max(abs(delx),abs(dely),abs(delz)).gt.ERRTOL)goto 1 e2=delx*dely-delz**2 e3=delx*dely*delz rf=(1.+(C1*e2-C2-C3*e3)*e2+C4*e3)/sqrt(ave) return END c ******************************************************** FUNCTION rj(x,y,z,p) implicit double precision (a-h,o-z) PARAMETER (ERRTOL=0.0015,TINY=2.5e-13,BIG=9.E11,C1=3./14., & C2=1./3.,C3=3./22.,C4=3./26.,C5=.75*C3,C6=1.5*C4,C7=.5*C2, & C8=C3+C3) CU USES rc,rf ! if(min(x,y,z).lt.0..or.min(x+y,x+z,y+z,abs(p)).lt.TINY.or.max(x,y, *z,abs(p)).gt.BIG) then stop 'invalid arguments in rj' end if ! sum=0. fac=1. if(p.gt.0.)then xt=x yt=y zt=z pt=p else xt=min(x,y,z) zt=max(x,y,z) yt=x+y+z-xt-zt a=1./(yt-p) b=a*(zt-yt)*(yt-xt) pt=yt+b rho=xt*zt/yt tau=p*pt/yt rcx=rc(rho,tau) endif 1 continue sqrtx=sqrt(xt) sqrty=sqrt(yt) sqrtz=sqrt(zt) alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz alpha=(pt*(sqrtx+sqrty+sqrtz)+sqrtx*sqrty*sqrtz)**2 beta=pt*(pt+alamb)**2 sum=sum+fac*rc(alpha,beta) fac=.25*fac xt=.25*(xt+alamb) yt=.25*(yt+alamb) zt=.25*(zt+alamb) pt=.25*(pt+alamb) ave=.2*(xt+yt+zt+pt+pt) delx=(ave-xt)/ave dely=(ave-yt)/ave delz=(ave-zt)/ave delp=(ave-pt)/ave if(max(abs(delx),abs(dely),abs(delz),abs(delp)).gt.ERRTOL)goto 1 ea=delx*(dely+delz)+dely*delz eb=delx*dely*delz ec=delp**2 ed=ea-3.*ec ee=eb+2.*delp*(ea-ec) rj=3.*sum+fac*(1.+ed*(-C1+C5*ed-C6*ee)+eb*(C7+delp*(-C8+delp*C4))+ *delp*ea*(C2-delp*C3)-C2*delp*ec)/(ave*sqrt(ave)) if (p.le.0.) rj=a*(b*rj+3.*(rcx-rf(xt,yt,zt))) return END ! ********************************************************* FUNCTION rtsafe(funcd,x1,x2,xacc) implicit double precision (a-h,o-z) INTEGER MAXIT EXTERNAL funcd PARAMETER (MAXIT=100) INTEGER j ! call funcd(x1,fl,df) call funcd(x2,fh,df) if((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.))stop *'root must be bracketed in rtsafe' if(fl.eq.0.)then rtsafe=x1 return else if(fh.eq.0.)then rtsafe=x2 return else if(fl.lt.0.)then xl=x1 xh=x2 else xh=x1 xl=x2 endif rtsafe=.5*(x1+x2) dxold=abs(x2-x1) dx=dxold call funcd(rtsafe,f,df) do 11 j=1,MAXIT if(((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f).ge.0..or. abs(2.* *f).gt.abs(dxold*df) ) then dxold=dx dx=0.5*(xh-xl) rtsafe=xl+dx if(xl.eq.rtsafe)return else dxold=dx dx=f/df temp=rtsafe rtsafe=rtsafe-dx if(temp.eq.rtsafe)return endif if(abs(dx).lt.xacc) return call funcd(rtsafe,f,df) if(f.lt.0.) then xl=rtsafe else xh=rtsafe endif 11 continue stop 'rtsafe exceeding maximum iterations' return END ! ********************************************************* SUBROUTINE trapzd(func,a,b,s,n) implicit double precision (a-h,o-z) INTEGER n EXTERNAL func INTEGER it,j ! if (n.eq.1) then s=0.5*(b-a)*(func(a)+func(b)) else it=2**(n-2) tnm=it del=(b-a)/tnm x=a+0.5*del sum=0. do 11 j=1,it sum=sum+func(x) x=x+del 11 continue s=0.5*(s+(b-a)*sum/tnm) endif return END ! ********************************************************* subroutine TRILINEAR(u,v,w,jx,jy,jz,fg,mx,my,mz,f) ! ! Do tri-linear interpolation ! implicit double precision (a-h,o-z) dimension fg(mx,my,mz) ! f = (1.-u)*(1.-v)*(1.-w) * fg(jx,jy,jz) & + (1.-u)*(1.-v)*w * fg(jx,jy,jz+1) & + (1.-u)*v*w * fg(jx,jy+1,jz+1) & + (1.-u)*v*(1.-w) * fg(jx,jy+1,jz) & + u*(1.-v)*(1.-w) * fg(jx+1,jy,jz) & + u*(1.-v)*w * fg(jx+1,jy,jz+1) & + u*v*w * fg(jx+1,jy+1,jz+1) & + u*v*(1.-w) * fg(jx+1,jy+1,jz) ! return end ! ********************************************************* subroutine TRIPOLY(kx,ky,kz,xg,yg,zg,fg,mp,np,pp,kk, & x,y,z,h,f,df) ! ! Polynomial interpolation on 3-dimensional grid ! ! cf. Num. Rec., 2nd ed, p.103,118; uses POLINT ! ! kx,ky,kz: starting location in 3D grid arrays ! xg(mp),yg(np),zg(pp): full coordinate arrays ! fg(mp,np,pp): full field array ! kk: polynomial order {1,2,3} ! x,y,z: observation point ! h: curvature ! f: field value at observation point ! df: error estimate for returned field ! implicit double precision (a-h,o-z) integer pp dimension fg(mp,np,pp),xg(mp),yg(np),zg(pp),zgh(pp) dimension f1tmp(4),f2tmp(4),f3tmp(4) ! do 300 k=1,kk+1 ! loop over z grid points do 200 j=1,kk+1 ! loop over y grid points ! ! Load temporary array f1tmp with x grid fields for fixed y and z do 100 i=1,kk+1 f1tmp(i) = fg(kx+i-1,ky+j-1,kz+k-1) 100 continue ! ! f2tmp = field interpolated over x as a function of y for fixed z call POLINT(xg(kx),f1tmp,kk+1,x,f2tmp(j),df) ! 200 continue ! ! f3tmp = field interpolated over y as a function of z call POLINT(yg(ky),f2tmp,kk+1,y,f3tmp(k),df) ! 300 continue ! ! In cylidrical coordinate system we need to rescale the step sizes ! for different radial positions do i=1,pp zgh(i) = zg(i) * (1. + h*x) end do ! ! f = field interpolated over z call POLINT(zgh(kz),f3tmp,kk+1,z,f,df) ! return end