C23456789012345678901234567890123456789012345678901234567890123456789012
C
C LAST MODIFIED October, 2000
C MODIFICATIONS: (OCT. 23, 1989) RAND, RANDN, IRANDP, JACOB3, DET, LUDCMP ADDED
C WARNING ADDED TO JACOBI
C (OCT. 26, 1989) CHANGE IN FUNCTION SUBROUTINES TO BE COMPATIBLE WITH 3B15.
C (MARCH 14, 1991) RAND CHANGED TO RANDU. RAN0 ADDED
C (MARCH 15, 1991) CALLING SEQUENCE OF RANDOM NUMBER GENERATORS RANDN,
C AND IRANDP CHANGED. ERF CHANGED TO ERF0. DIMENSION Y(1) IN AMOEBA
C CHANGED TO Y(2)
C (APRIL 29, 1991) RAN0 MODIFIED TO RETURN IN IDUM THE LAST VALUE OF ISEED.
c (May 19,1992) dimensions changed in ludcmp, changed to double precision
C (July,1995) corrected typo in zbrent, t0l1 should be tol1
c (October, 1999) Subroutine zbrent: zbrent=0. added to error return (when FB*FA.GT.0.) in zbrent
c Subroutine ludcmp: implicit double precision statement put before parameter statement
c (October, 2000) Warning to JACOBI is modified
c***********************************************************************
C
C S U B R O U T I N E A M O E B A
C
C MODIFICATION OF PRESS ET AL: NUMBERICAL RECIPES, PAGES 292-293
C SIMPLEX METHOD OF NELDER AND MEAD
C
C CALLING SEQUENCE:
C P--MP BY (NDIM+1) MATRIX WHOSE COLUMNS ARE LINEARLY INDEPENDENT
C AND CLOSE TO INITAL GUESS. ON OUTPUT, THE SIMPLEX SPANNED BY
C THE COLUMNS OF P CONTAINS THE MINIMUM FOUND BY AMOEBA.
C Y--VECTOR OF LENGTH (NDIM+1) CONTAINING THE VALUES OF THE FUNCTION
C AT THE COLUMNS OF P.
C MP--ROW DIMENSION OF P AS SPECIFIED IN CALLING PROGRAM'S DIMENSION
C STATEMENT
C NDIM--NUMBER OF PARAMETERS TO BE MINIMIZED
C FTOL--INPUT: FUNCTION VALUES AT THE COLUMNS OF P SHOULD BE WITHIN
C FTOL OF MINIMUM VALUE
C OUTPUT: ESTIMATED ACHIEVED VALUE OF FTOL
C FUNK--NAME OF FUNCTION TO BE EVALUATED
C ITER--INPUT: MAXIMUM NUMBER OF ITERATIONS.
C OUTPUT: ACTUAL NUMBER OF ITERATIONS.
C WORK--WORK VECTOR OF LENGTH 3*NDIM
C
C
SUBROUTINE AMOEBA(P,Y,MP,NDIM,FTOL,FUNK,ITER,WORK)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
PARAMETER(ALPHA=1.0,BETA=0.5,GAMMA=2.0)
DIMENSION P(MP,1),Y(2),WORK(NDIM,3)
MPTS=NDIM+1
ITMAX=ITER
DO 2 I=1,MPTS
2 Y(I)=FUNK(P(1,I))
ITER=0
1 ILO=1
IF (Y(1).GT.Y(2)) THEN
IHI=1
INHI=2
ELSE
IHI=2
INHI=1
ENDIF
DO 11 I=1,MPTS
IF (Y(I).LT.Y(ILO)) ILO=I
IF (Y(I).GT.Y(IHI)) THEN
INHI=IHI
IHI=I
ELSE IF (Y(I).GT.Y(INHI))THEN
IF (I.NE.IHI) INHI=I
ENDIF
11 CONTINUE
RTOL=2.*ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))
IF ((RTOL.LT.FTOL).OR.(ITER.GE.ITMAX)) THEN
FTOL=RTOL
IF (ILO.EQ.1) RETURN
DO 3 I=1,NDIM
TEMP=P(I,1)
P(I,1)=P(I,ILO)
3 P(I,ILO)=TEMP
TEMP=Y(1)
Y(1)=Y(ILO)
Y(ILO)=TEMP
RETURN
ENDIF
C
ITER=ITER+1
DO 12 J=1,NDIM
12 WORK(J,1)=0.
DO 14 I=1,MPTS
IF (I.NE.IHI) THEN
DO 13 J=1,NDIM
13 WORK(J,1)=WORK(J,1)+P(J,I)
ENDIF
14 CONTINUE
DO 15 J=1,NDIM
WORK(J,1)=WORK(J,1)/NDIM
15 WORK(J,2)=(1.+ALPHA)*WORK(J,1)-ALPHA*P(J,IHI)
YPR=FUNK(WORK(1,2))
IF (YPR.LE.Y(ILO)) THEN
DO 16 J=1,NDIM
16 WORK(J,3)=GAMMA*WORK(J,2)+(1.-GAMMA)*WORK(J,1)
YPRR=FUNK(WORK(1,3))
IF (YPRR.LT.Y(ILO)) THEN
DO 17 J=1,NDIM
17 P(J,IHI)=WORK(J,3)
Y(IHI)=YPRR
ELSE
DO 18 J=1,NDIM
18 P(J,IHI)=WORK(J,2)
Y(IHI)=YPR
ENDIF
ELSE IF (YPR.GE.Y(INHI)) THEN
IF (YPR.LT.Y(IHI)) THEN
DO 19 J=1,NDIM
19 P(J,IHI)=WORK(J,2)
Y(IHI)=YPR
ENDIF
DO 21 J=1,NDIM
21 WORK(J,3)=BETA*P(J,IHI)+(1.-BETA)*WORK(J,1)
YPRR=FUNK(WORK(1,3))
IF (YPRR.LT.Y(IHI)) THEN
DO 22 J=1,NDIM
22 P(J,IHI)=WORK(J,3)
Y(IHI)=YPRR
ELSE
DO 24 I=1,MPTS
IF (I.NE.ILO) THEN
DO 23 J=1,NDIM
WORK(J,2)=0.5*(P(J,I)+P(J,ILO))
23 P(J,I)=WORK(J,2)
Y(I)=FUNK(WORK(1,2))
ENDIF
24 CONTINUE
ENDIF
ELSE
DO 25 J=1,NDIM
25 P(J,IHI)=WORK(J,2)
Y(IHI)=YPR
ENDIF
GO TO 1
END
C***********************************************************************
C*** FAST ASCENDING SORT ROUTINE ***
SUBROUTINE SORT(A,N)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION A(N)
M2=2*N-1
5 M2=INT(M2/2)
IF(M2.EQ.0)GOTO 10
K=N-M2
DO 20 I=1,K
DO 30 J2=1,I,M2
N2=I-(J2-1)
L2=N2+M2
IF(A(L2).GT.A(N2)) GOTO 20
B2=A(N2)
A(N2)=A(L2)
A(L2)=B2
30 CONTINUE
20 CONTINUE
GOTO 5
10 CONTINUE
RETURN
END
C***********************************************************************
SUBROUTINE SL(A,N,IPT)
C N=NUMBER OF ROWS OF A TO BE SORTED
C
C This routine sorts a vector A into ascending order.
C The contents of the vector IPT are shuffled accordingly.
C
C ... don't bother asking how it works.
C
C_______________________________________________________________________
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION A(N),IPT(N)
M2=2*N-1
5 M2=INT(M2/2)
IF(M2.EQ.0)GOTO 10
K=N-M2
DO 20 I=1,K
DO 30 J2=1,I,M2
N2=I-(J2-1)
L2=N2+M2
IF(A(L2).GT.A(N2))GOTO 20
B2=A(N2)
IB3=IPT(N2)
A(N2)=A(L2)
IPT(N2)=IPT(L2)
A(L2)=B2
IPT(L2)=IB3
30 CONTINUE
20 CONTINUE
GOTO 5
10 CONTINUE
RETURN
END
C***********************************************************************
SUBROUTINE JACOBI(A,N,NA,D,V,NV,WORK,NROT)
C Computes the eigenvalues and eigenvectors of the N by N matrix A.
C Modification of JACOBI in Press et al, Numerical Recipes, pp 346-348.
C A is stored in FULL storage mode.
C NA=row dimension of A in declarations of calling program
C D output vector of eigenvalues in ascending order
C V=output matrix whose columns are the eigenvectors of A
C NV=row dimension of V in calling program
C WORK=work vector of length at least 2*N
C NROT=output number of Jacobi revolutions
C
C **************************************************************
C * *
C * W A R N I N G -- A C H T U N G *
C * *
C * JACOBI OVERWRITES THE UPPER TRIANGULAR PORTION OF MATRIX A *
C * ITS LAST STEP IS TO RESET THE UPPER TRIANGULAR PORTION TO *
C * THE LOWER TRIANGULAR PORTION. THUS IF A IS SYMMETRIC (AS *
C * IT SHOULD BE), IT IS PRESERVED. *
C * *
C **************************************************************
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION A(NA,1),D(1),V(NV,1),WORK(1)
DO 12 IP=1,N
DO 11 IQ=1,N
11 V(IP,IQ)=0.D0
12 V(IP,IP)=1.
DO 13 IP=1,N
WORK(IP)=A(IP,IP)
D(IP)=WORK(IP)
IPN=IP+N
13 WORK(IPN)=0.D0
NROT=0
DO 24 I=1,50
SM=0.D0
DO 15 IP=1,N-1
DO 15 IQ=IP+1,N
15 SM=SM+ABS(A(IP,IQ))
IF (SM.EQ.0.) GO TO 30
IF (I.LT.4) THEN
THRESH=0.2D0*SM/N**2
ELSE
THRESH=0.D0
ENDIF
DO 22 IP=1,N-1
DO 22 IQ=IP+1,N
G=100.*ABS(A(IP,IQ))
IF ((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP))).AND.
& (ABS(D(IQ))+G.EQ.ABS(D(IQ)))) THEN
A(IP,IQ)=0.D0
ELSE IF (ABS(A(IP,IQ)).GT.THRESH) THEN
H=D(IQ)-D(IP)
IF (ABS(H)+G.EQ.ABS(H)) THEN
T=A(IP,IQ)/H
ELSE
THETA=0.5D0*H/A(IP,IQ)
T=1.D0/(ABS(THETA)+SQRT(1.D0+THETA**2))
IF (THETA.LT.0.D0) T=-T
ENDIF
C=1.D0/SQRT(1.D0+T**2)
S=T*C
TAU=S/(1.D0+C)
H=T*A(IP,IQ)
IPN=IP+N
WORK(IPN)=WORK(IPN)-H
IQN=IQ+N
WORK(IQN)=WORK(IQN)+H
D(IP)=D(IP)-H
D(IQ)=D(IQ)+H
A(IP,IQ)=0.D0
IF (IP.GT.1) THEN
DO 16 J=1,IP-1
G=A(J,IP)
H=A(J,IQ)
A(J,IP)=G-S*(H+G*TAU)
16 A(J,IQ)=H+S*(G-H*TAU)
ENDIF
IF (IQ.GT.(IP+1)) THEN
DO 17 J=IP+1,IQ-1
G=A(IP,J)
H=A(J,IQ)
A(IP,J)=G-S*(H+G*TAU)
17 A(J,IQ)=H+S*(G-H*TAU)
ENDIF
IF (IQ.LT.N) THEN
DO 18 J=IQ+1,N
G=A(IP,J)
H=A(IQ,J)
A(IP,J)=G-S*(H+G*TAU)
18 A(IQ,J)=H+S*(G-H*TAU)
ENDIF
DO 19 J=1,N
G=V(J,IP)
H=V(J,IQ)
V(J,IP)=G-S*(H+G*TAU)
19 V(J,IQ)=H+S*(G-H*TAU)
NROT=NROT+1
ENDIF
22 CONTINUE
DO 23 IP=1,N
IPN=IP+N
WORK(IP)=WORK(IP)+WORK(IPN)
D(IP)=WORK(IP)
23 WORK(IPN)=0.D0
24 CONTINUE
NROT=-NROT
C
30 DO 31 IP=1,N-1
DO 31 IQ=IP+1,N
31 A(IP,IQ)=A(IQ,IP)
DO 33 I=1,N-1
K=I
P=D(I)
DO 41 J=I+1,N
IF (D(J).LT.P) THEN
K=J
P=D(J)
ENDIF
41 CONTINUE
IF (K.NE.I) THEN
D(K)=D(I)
D(I)=P
DO 32 J=1,N
P=V(J,I)
V(J,I)=V(J,K)
32 V(J,K)=P
ENDIF
33 CONTINUE
RETURN
END
C
SUBROUTINE JACOB2(A,N,NA,D,V,NV,WORK,NROT)
C Single precision form of JACOBI
DIMENSION A(NA,1),D(1),V(NV,1),WORK(1)
DO 12 IP=1,N
DO 11 IQ=1,N
11 V(IP,IQ)=0.D0
12 V(IP,IP)=1.
DO 13 IP=1,N
WORK(IP)=A(IP,IP)
D(IP)=WORK(IP)
IPN=IP+N
13 WORK(IPN)=0.D0
NROT=0
DO 24 I=1,50
SM=0.D0
DO 15 IP=1,N-1
DO 15 IQ=IP+1,N
15 SM=SM+ABS(A(IP,IQ))
IF (SM.EQ.0.) GO TO 30
IF (I.LT.4) THEN
THRESH=0.2D0*SM/N**2
ELSE
THRESH=0.D0
ENDIF
DO 22 IP=1,N-1
DO 22 IQ=IP+1,N
G=100.*ABS(A(IP,IQ))
IF ((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP))).AND.
& (ABS(D(IQ))+G.EQ.ABS(D(IQ)))) THEN
A(IP,IQ)=0.D0
ELSE IF (ABS(A(IP,IQ)).GT.THRESH) THEN
H=D(IQ)-D(IP)
IF (ABS(H)+G.EQ.ABS(H)) THEN
T=A(IP,IQ)/H
ELSE
THETA=0.5D0*H/A(IP,IQ)
T=1.D0/(ABS(THETA)+SQRT(1.D0+THETA**2))
IF (THETA.LT.0.D0) T=-T
ENDIF
C=1.D0/SQRT(1.D0+T**2)
S=T*C
TAU=S/(1.D0+C)
H=T*A(IP,IQ)
IPN=IP+N
WORK(IPN)=WORK(IPN)-H
IQN=IQ+N
WORK(IQN)=WORK(IQN)+H
D(IP)=D(IP)-H
D(IQ)=D(IQ)+H
A(IP,IQ)=0.D0
IF (IP.GT.1) THEN
DO 16 J=1,IP-1
G=A(J,IP)
H=A(J,IQ)
A(J,IP)=G-S*(H+G*TAU)
16 A(J,IQ)=H+S*(G-H*TAU)
ENDIF
IF (IQ.GT.(IP+1)) THEN
DO 17 J=IP+1,IQ-1
G=A(IP,J)
H=A(J,IQ)
A(IP,J)=G-S*(H+G*TAU)
17 A(J,IQ)=H+S*(G-H*TAU)
ENDIF
IF (IQ.LT.N) THEN
DO 18 J=IQ+1,N
G=A(IP,J)
H=A(IQ,J)
A(IP,J)=G-S*(H+G*TAU)
18 A(IQ,J)=H+S*(G-H*TAU)
ENDIF
DO 19 J=1,N
G=V(J,IP)
H=V(J,IQ)
V(J,IP)=G-S*(H+G*TAU)
19 V(J,IQ)=H+S*(G-H*TAU)
NROT=NROT+1
ENDIF
22 CONTINUE
DO 23 IP=1,N
IPN=IP+N
WORK(IP)=WORK(IP)+WORK(IPN)
D(IP)=WORK(IP)
23 WORK(IPN)=0.D0
24 CONTINUE
NROT=-NROT
C
30 DO 31 IP=1,N-1
DO 31 IQ=IP+1,N
31 A(IP,IQ)=A(IQ,IP)
DO 33 I=1,N-1
K=I
P=D(I)
DO 41 J=I+1,N
IF (D(J).LT.P) THEN
K=J
P=D(J)
ENDIF
41 CONTINUE
IF (K.NE.I) THEN
D(K)=D(I)
D(I)=P
DO 32 J=1,N
P=V(J,I)
V(J,I)=V(J,K)
32 V(J,K)=P
ENDIF
33 CONTINUE
RETURN
END
C***********************************************************************
C
C SUBROUTINE X2TAB FOR CHI STAT LOOK-UP
C
C***********************************************************************
C
SUBROUTINE X2TAB(PC,IDF,X2VAL)
REAL TABLE(22,5),CENT(5)
DATA CENT/.5,.7,.9,.95,.99/
DATA TABLE/
+ .455,1.386,2.366,3.357,4.351,5.348,6.346,7.344,8.343
+ ,9.342,10.341,11.34,12.34,13.339,14.339,15.338,16.338,
+ 17.338,18.338,19.337,20.337,21.337,
& 1.074,2.408,3.665,4.878,6.064,7.231,8.383,
+ 9.524,10.656,11.781,12.899,14.011,15.119,16.222,
+ 17.322,18.418,19.511,20.601,21.689,22.775,23.858,24.939,
& 2.706,4.605,6.251,7.779,9.236,10.645,12.017,13.362,
+ 14.684,15.987,17.275,18.549,19.812,21.064,22.307,
+ 23.542,24.769,25.989,27.204,28.412,29.615,30.813,
& 3.841,5.991,7.815,9.488,11.07,12.592,14.067,15.507,
+ 16.919,18.307,19.675,21.026,22.362,23.685,24.996,
+ 26.296,27.587,28.869,30.144,31.41,32.671,33.924,
& 6.635,9.21,11.341,13.277,15.086,16.812,18.475,
+ 20.09,21.666,23.209,24.725,26.217,27.688,29.141,
+ 30.578,32.,33.409,34.805,36.191,37.566,38.932,40.289/
C
DO 10 J=1,5
IF(ABS(CENT(J)-PC).LT..0005)GOTO 20
10 CONTINUE
20 CONTINUE
40 X2VAL=TABLE(IDF,J)
RETURN
END
C***********************************************************************
C SPECIAL FUNCTION ROUTINES
C MODIFIED FROM PRESS ET AL: NUMERICAL RECIPES
C***********************************************************************
FUNCTION GAMMP(A,X,IER)
C INCOMPLETE GAMMA FUNCTION
IF (X.LT.0..OR.A.LE.0.) THEN
IER=100
GAMMP=0.
RETURN
ENDIF
IF (X.LT.A+1.) THEN
CALL GSER(GAMMPT,A,X,GLN,IER)
GAMMP=GAMMPT
ELSE
CALL GCF(GAMMCF,A,X,GLN,IER)
GAMMP=1.-GAMMCF
ENDIF
RETURN
END
C
FUNCTION GAMMQ(A,X,IER)
IF (X.LT.0..OR.A.LE.0.) THEN
IER=100
GAMMQ=0.
RETURN
ENDIF
IF (X.LT.A+1.) THEN
CALL GSER(GAMSER,A,X,GLN,IER)
GAMMQ=1.-GAMSER
ELSE
CALL GCF(GAMMQT,A,X,GLN,IER)
GAMMQ=GAMMQT
ENDIF
RETURN
END
C
SUBROUTINE GSER(GAMSER,A,X,GLN,IER)
C INCOMPLETE GAMMA FUNCTION USING SERIES REPRESENTATION; GLN=LN GAMMA(A)
PARAMETER (ITMAX=100,EPS=3.E-7)
GLN=GAMMLN(A)
IF (X.LE.0.) THEN
IF (X.LT.0.) IER=100
GAMSER=0.
RETURN
ENDIF
AP=A
SUM=1./A
DEL=SUM
DO 11 N=1,ITMAX
AP=AP+1.
DEL=DEL*X/AP
SUM=SUM+DEL
IF (ABS(DEL).LT.ABS(SUM)*EPS) GO TO 1
11 CONTINUE
IER=101
1 GAMSER=SUM*EXP(-X+A*LOG(X)-GLN)
RETURN
END
C
SUBROUTINE GCF(GAMMCF,A,X,GLN,IER)
C CONTINUED FRACTION CALCULATION OF INCOMPLETE GAMMA FUNCTION
PARAMETER (ITMAX=100,EPS=3.E-7)
GLN=GAMMLN(A)
GOLD=0.
A0=1.
A1=X
B0=0.
B1=1.
FAC=1.
DO 11 N=1,ITMAX
AN=FLOAT(N)
ANA=AN-A
A0=(A1+A0*ANA)*FAC
B0=(B1+B0*ANA)*FAC
ANF=AN*FAC
A1=X*A0+ANF*A1
B1=X*B0+ANF*B1
IF (A1.NE.0.) THEN
FAC=1./A1
G=B1*FAC
IF (ABS((G-GOLD)/G).LT.EPS) GO TO 1
GOLD=G
ENDIF
11 CONTINUE
IER=102
1 GAMMCF=EXP(-X+A*ALOG(X)-GLN)*G
RETURN
END
C
FUNCTION GAMMLN(XX)
DOUBLE PRECISION COF(6),STP,HALF,ONE,FPF,X,TMP,SER,PI,DGAMM,Z
DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0,
&-1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/
DATA HALF,ONE,FPF,PI/.5D0,1.D0,5.5D0,3.141592653589793/
IF (XX.EQ.1.) THEN
GAMMLN=0.
RETURN
ENDIF
X=ABS(XX-ONE)
Z=X
TMP=X+FPF
TMP=(X+HALF)*LOG(TMP)-TMP
SER=ONE
DO 11 J=1,6
X=X+ONE
SER=SER+COF(J)/X
11 CONTINUE
DGAMM=TMP+LOG(STP*SER)
IF (XX.GT.1.) THEN
GAMMLN=DGAMM
ELSE
GAMMLN=LOG(PI*Z/SIN(PI*Z))-DGAMM
ENDIF
RETURN
END
C
FUNCTION ERF0(X,IER)
IF (X.LT.0) THEN
ERF0=-GAMMP(.5,X**2,IER)
ELSE
ERF0=GAMMP(.5,X**2,IER)
ENDIF
RETURN
END
C
FUNCTION ERFC(X,IER)
IF (X.LT.0) THEN
ERFC=1.+GAMMP(.5,X**2,IER)
ELSE
ERFC=GAMMQ(.5,X**2,IER)
ENDIF
RETURN
END
C
FUNCTION ERFCC(X)
Z=ABS(X)
T=1./(1.+0.5*Z)
ERFCC=T*EXP(-Z*Z-1.26551223+T*(1.00002368+T*(.37409196+
& T*(.09678418+T*(.27886807+T*(-1.13520398+
& T*(1.48851587+T*(-.82215223+T*.17087277))))))))
IF (X.LT.0.) ERFCC=2.-ERFCC
RETURN
END
C
SUBROUTINE XDCH(PLEV,DF,XCHI,IER)
DIMENSION DF(1)
PLEV=GAMMP(DF(1)/2.,XCHI/2.,IER)
RETURN
END
C
FUNCTION BETAI(A,B,X,IER)
IF (X.LE.0.) THEN
BETAI=0.
IF (X.LT.0.) IER=105
ELSE IF (X.GE.1.) THEN
BETAI=1.
IF (X.GT.1.) IER=105
ELSE
BT=EXP(GAMMLN(A+B)-GAMMLN(A)-GAMMLN(B)+A*ALOG(X)+B*ALOG(1.-X))
IF (X.LT.((A+1.)/(A+B+2.))) THEN
BETAI=BT*BETACF(A,B,X,IER)/A
ELSE
BETAI=1.-BT*BETACF(B,A,1.-X,IER)/B
ENDIF
ENDIF
RETURN
END
C
FUNCTION BETACF(A,B,X,IER)
PARAMETER (ITMAX=100,EPS=3.E-7)
AM=1.
BM=1.
AZ=1.
QAB=A+B
QAP=A+1.
QAM=A-1.
BZ=1.-QAB*X/QAP
DO 11 M=1,ITMAX
EM=M
TEM=EM+EM
D=EM*(B-M)*X/((QAM+TEM)*(A+TEM))
AP=AZ+D*AM
BP=BZ+D*BM
D=-(A+EM)*(QAB+EM)*X/((A+TEM)*(QAP+TEM))
APP=AP+D*AZ
BPP=BP+D*BZ
AOLD=AZ
AM=AP/BPP
BM=BP/BPP
AZ=APP/BPP
BZ=1.
IF (ABS(AZ-AOLD).LT.EPS*ABS(AZ)) GO TO 1
11 CONTINUE
IER=106
1 BETACF=AZ
RETURN
END
C
SUBROUTINE XDF(PLEV,D,XF,IER)
DIMENSION D(2)
X=D(2)/(D(2)+D(1)*XF)
PLEV=1.-BETAI(D(2)/2.,D(1)/2.,X,IER)
RETURN
END
C
SUBROUTINE XDN(PLEV,DUMMY,XN,IER)
DIMENSION DUMMY(1)
X=XN/1.4142136
PLEV=.5+.5*ERF0(X,IER)
RETURN
END
C
C***********************************************************************
C
C Root finding using Brent algorithms
C Modification of Press et al, Numerical Recipes, pages 253-254
C
FUNCTION ZBRENT(FUNC,FVAL,DUMMY,X1,X2,TOL,NER)
PARAMETER (ITMAX=100,EPS=3.E-8)
DIMENSION DUMMY(1)
EXTERNAL FUNC
C
A=X1
B=X2
NER=0
IER=0
CALL FUNC(FA,DUMMY,A,IER)
FA=FA-FVAL
IF (IER.NE.0) NER=NER+1
IER=0
CALL FUNC(FB,DUMMY,B,IER)
FB=FB-FVAL
IF (IER.NE.0) NER=NER+1
IF (FB*FA.GT.0.) THEN
NER=-NER-1
zbrent=0.
RETURN
ENDIF
FC=FB
DO 11 ITER=1,ITMAX
IF (FB*FC.GT.0.) THEN
C=A
FC=FA
D=B-A
E=D
ENDIF
IF (ABS(FC).LT.ABS(FB)) THEN
A=B
B=C
C=A
FA=FB
FB=FC
FC=FA
ENDIF
TOL1=2*EPS*ABS(B)+0.5*TOL
XM=.5*(C-B)
IF (ABS(XM).LE.TOL1.OR.FB.EQ.0.) THEN
ZBRENT=B
RETURN
ENDIF
IF (ABS(E).GE.TOL1.AND.ABS(FA).GT.ABS(FB)) THEN
S=FB/FA
IF (A.EQ.C) THEN
P=2.*XM*S
Q=1.-S
ELSE
Q=FA/FC
R=FB/FC
P=S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.))
Q=(Q-1.)*(R-1.)*(S-1.)
ENDIF
IF (P.GT.0.) Q=-Q
P=ABS(P)
IF (2.*P .LT. MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN
E=D
D=P/Q
ELSE
D=XM
E=D
ENDIF
ELSE
D=XM
E=D
ENDIF
A=B
FA=FB
IF (ABS(D).GT.TOL1) THEN
B=B+D
ELSE
B=B+SIGN(TOL1,XM)
ENDIF
IER=0
CALL FUNC(FB,DUMMY,B,IER)
FB=FB-FVAL
IF (IER.NE.0) NER=NER+1
11 CONTINUE
NER=-1000*NER
ZBRENT=B
RETURN
END
C
C***********************************************************************
C INVERSE DISTRIBUTION ROUTINES
C
SUBROUTINE XIDN(PLEV,XN,IER)
DIMENSION P(14),TABLE(14),DUMMY(1)
EXTERNAL XDN
DATA TOL/.0001/
DATA P/.5,.6,.7,.8,.85,.9,.95,.975,.9875,.99,.995,.9975,.999,
& .9995/
DATA TABLE/0.,.253,.524,.842,1.036,1.282,1.645,1.960,2.240,
& 2.326,2.576,2.807,3.090,3.291/
C
IER=0
QLEV=.5+ABS(PLEV-.5)
DO 100 I=1,14
IF (ABS(QLEV-P(I)).LT..00049) THEN
XN=TABLE(I)
GO TO 200
ENDIF
IF (QLEV.LT.P(I)) THEN
XN=ZBRENT(XDN,QLEV,DUMMY,TABLE(I-1),TABLE(I),TOL,IER)
GO TO 200
ENDIF
100 CONTINUE
XN=TABLE(14)
IER=111111
200 IF (PLEV.GE.0.5) RETURN
XN=-XN
RETURN
END
C
SUBROUTINE XIDCH(PLEV,DF,XCHI,IER)
EXTERNAL XDCH
DIMENSION D(1),P(5)
DATA TOL/.0001/,P/.5,.7,.9,.95,.99/
IER=0
IF (DF.LE.0.) THEN
XCHI=0.
IER=111111
ENDIF
D(1)=DF
IDF=DF+.001
IF ((ABS(FLOAT(IDF)-DF).LT..001).AND.(IDF.LE.22)) THEN
DO 10 I=1,5
IF (ABS(PLEV-P(I)).LT..00049) THEN
CALL X2TAB(PLEV,IDF,XCHI)
RETURN
ENDIF
IF (I.NE.1) THEN
IF ((PLEV.GT.P(I-1)).AND.(PLEV.LT.P(I))) THEN
CALL X2TAB(P(I-1),IDF,A)
CALL X2TAB(P(I),IDF,B)
XCHI=ZBRENT(XDCH,PLEV,D,A,B,TOL,IER)
RETURN
ENDIF
ENDIF
10 CONTINUE
ENDIF
CALL XIDN(PLEV,B,IER)
B=DF+SQRT(2.*DF)*B
IF (IER.NE.0) B=DF
B=MAX(1.,B)
A=B
20 CALL XDCH(PA,D,A,IER)
IF (PA.GE.PLEV) THEN
A=.9*A
GO TO 20
ENDIF
30 CALL XDCH(PB,D,B,IER)
IF (PB.LE.PLEV) THEN
B=1.1*B
GO TO 30
ENDIF
XCHI=ZBRENT(XDCH,PLEV,D,A,B,TOL,IER)
RETURN
END
C
SUBROUTINE XIDF(PLEV,D1,D2,XF,IER)
EXTERNAL XDF
DIMENSION D(2)
DATA TOL/.0001/
D(1)=D1
D(2)=D2
IF ((PLEV.LE.0.).OR.(PLEV.GE.1.).OR.(D1.LE.0.).OR.(D2.LE.0.)) THEN
XF=1.
IER=111111
RETURN
ENDIF
CALL XIDCH(PLEV,D(1),XF,IER)
A=XF
10 CALL XDF(PA,D,A,IER)
IF (PA.GE.PLEV) THEN
A=A/2.
GO TO 10
ENDIF
B=XF
20 CALL XDF(PB,D,B,IER)
IF (PB.LE.PLEV) THEN
B=2.*B
GO TO 20
ENDIF
XF=ZBRENT(XDF,PLEV,D,A,B,TOL,IER)
RETURN
END
C
C***********************************************************************
C
REAL FUNCTION RANDU(IX)
C
C UNIFORM RANDOM NUMBER GENERATOR.
C REF: SCHRAGE, L. 1979 'A MORE PORTABLE FORTRAN RANDOM NUMBER
C GENERATOR'
C A. C. M. TRANSACTIONS ON MATHEMATICAL SOFTWARE V5 132-138.
C
INTEGER A,P,IX,B15,B16,XHI,XALO,LEFTLO,FHI,K
DATA A/16807/,B15/32768/,B16/65536/,P/2147483647/
XHI=IX/B16
XALO=(IX-XHI*B16)*A
LEFTLO=XALO/B16
FHI=XHI*A+LEFTLO
K=FHI/B15
IX=(((XALO-LEFTLO*B16)-P)+(FHI-K*B15)*B16)+K
IF (IX.LT.0) IX=IX+P
RANDU=FLOAT(IX)*4.656612875E-10
RETURN
END
C
C***********************************************************************
C
SUBROUTINE RANDN(N,X,ISEED,FRAND)
C
C STANDARD NORMAL RANDOM NUMBER GENERATOR
C
C FRAND IS A FUNCTION SUBROUTINE GENERATING UNIFORM RANDOM NUMBERS
C
DIMENSION X(1)
IF (N.LE.0) RETURN
NGEN=1
ISEED1=ISEED
C
5 DO 10 I=1,1000
X1=2.*FRAND(ISEED)-1.
X2=2.*FRAND(ISEED)-1.
S=X1*X1+X2*X2
IF (S.LT.1.) GO TO 20
10 CONTINUE
WRITE(6,1000) ISEED1
RETURN
20 X1=X1*SQRT(-2.*ALOG(S)/S)
X2=X2*SQRT(-2.*ALOG(S)/S)
C
X(NGEN)=X1
IF (NGEN.EQ.N) RETURN
NGEN=NGEN+1
X(NGEN)=X2
IF (NGEN.EQ.N) RETURN
NGEN=NGEN+1
GO TO 5
C
1000 FORMAT('0ERROR IN RANDN. SEED NUMBER: ',I12)
END
C
C***********************************************************************
C
INTEGER FUNCTION IRANDP(RLAM,ISEED,FRAND)
C
C POISSON RANDOM NUMBER GENERATOR.
C
C FRAND IS A FUNCTION SUBROUTINE GENERATING UNIFORM RANDOM NUMBERS
C
1000 FORMAT('0ERROR IN IRANDP. SEED NUMBER: ',I12)
C
ISEED1=ISEED
S=0.
DO 10 I=1,1000
S=S-ALOG(FRAND(ISEED))
IF (S.GT.RLAM) GO TO 20
10 CONTINUE
WRITE(6,1000) ISEED1
IRANDP=1000
RETURN
20 IRANDP=I-1
RETURN
END
C
C***********************************************************************
C
FUNCTION RAN0(IDUM)
C MODIFICATION OF RANDOM NUMBER MIXING ROUTINE FROM
C 'NUMERICAL RECIPES'
DIMENSION V(97)
DATA IFF /0/
IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
IFF=1
ISEED=ABS(IDUM)
DO 11 J=1,97
DUM=RANDU(ISEED)
11 CONTINUE
DO 12 J=1,97
V(J)=RANDU(ISEED)
12 CONTINUE
Y=RANDU(ISEED)
ENDIF
J=1+INT(97.*Y)
IF(J.GT.97.OR.J.LT.1) THEN
WRITE(*,*) 'ERROR IN SUBROUTINE RAN0'
J=17
ENDIF
Y=V(J)
RAN0=Y
V(J)=RANDU(ISEED)
IDUM=ISEED
RETURN
END
C
C***********************************************************************
C
SUBROUTINE JACOB3(A,N,NA,D,V,NV,WORK,NROT)
C Single precision form of JACOBI, eigenvalues are not sorted
C Computes the eigenvalues and eigenvectors of the N by N matrix A.
C Modification of JACOBI in Press et al, Numerical Recipes, pp 346-348.
C A is stored in FULL storage mode.
C NA=row dimension of A in declarations of calling program
C D output vector of eigenvalues in ascending order
C V=output matrix whose columns are the eigenvectors of A
C NV=row dimension of V in calling program
C WORK=work vector of length at least 2*N
C NROT=output number of Jacobi revolutions
C
DIMENSION A(NA,1),D(1),V(NV,1),WORK(1)
DO 12 IP=1,N
DO 11 IQ=1,N
11 V(IP,IQ)=0.D0
12 V(IP,IP)=1.
DO 13 IP=1,N
WORK(IP)=A(IP,IP)
D(IP)=WORK(IP)
IPN=IP+N
13 WORK(IPN)=0.D0
NROT=0
DO 24 I=1,50
SM=0.D0
DO 15 IP=1,N-1
DO 15 IQ=IP+1,N
15 SM=SM+ABS(A(IP,IQ))
IF (SM.EQ.0.) GO TO 30
IF (I.LT.4) THEN
THRESH=0.2D0*SM/N**2
ELSE
THRESH=0.D0
ENDIF
DO 22 IP=1,N-1
DO 22 IQ=IP+1,N
G=100.*ABS(A(IP,IQ))
IF ((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP))).AND.
& (ABS(D(IQ))+G.EQ.ABS(D(IQ)))) THEN
A(IP,IQ)=0.D0
ELSE IF (ABS(A(IP,IQ)).GT.THRESH) THEN
H=D(IQ)-D(IP)
IF (ABS(H)+G.EQ.ABS(H)) THEN
T=A(IP,IQ)/H
ELSE
THETA=0.5D0*H/A(IP,IQ)
T=1.D0/(ABS(THETA)+SQRT(1.D0+THETA**2))
IF (THETA.LT.0.D0) T=-T
ENDIF
C=1.D0/SQRT(1.D0+T**2)
S=T*C
TAU=S/(1.D0+C)
H=T*A(IP,IQ)
IPN=IP+N
WORK(IPN)=WORK(IPN)-H
IQN=IQ+N
WORK(IQN)=WORK(IQN)+H
D(IP)=D(IP)-H
D(IQ)=D(IQ)+H
A(IP,IQ)=0.D0
IF (IP.GT.1) THEN
DO 16 J=1,IP-1
G=A(J,IP)
H=A(J,IQ)
A(J,IP)=G-S*(H+G*TAU)
16 A(J,IQ)=H+S*(G-H*TAU)
ENDIF
IF (IQ.GT.(IP+1)) THEN
DO 17 J=IP+1,IQ-1
G=A(IP,J)
H=A(J,IQ)
A(IP,J)=G-S*(H+G*TAU)
17 A(J,IQ)=H+S*(G-H*TAU)
ENDIF
IF (IQ.LT.N) THEN
DO 18 J=IQ+1,N
G=A(IP,J)
H=A(IQ,J)
A(IP,J)=G-S*(H+G*TAU)
18 A(IQ,J)=H+S*(G-H*TAU)
ENDIF
DO 19 J=1,N
G=V(J,IP)
H=V(J,IQ)
V(J,IP)=G-S*(H+G*TAU)
19 V(J,IQ)=H+S*(G-H*TAU)
NROT=NROT+1
ENDIF
22 CONTINUE
DO 23 IP=1,N
IPN=IP+N
WORK(IP)=WORK(IP)+WORK(IPN)
D(IP)=WORK(IP)
23 WORK(IPN)=0.D0
24 CONTINUE
NROT=-NROT
30 RETURN
END
C***********************************************************************
REAL FUNCTION DET(A,N,NA,IWORK,WORK)
C DETERMINANT OF A N BY N MATRIX A, STORED WITH ROW DIMENSION NA. IWORK
C AND WORK ARE WORK VECTORS OF LENGTH N.
C MODIFICATION OF PRESS ET AL, PAGE 39.
C NOTE: A IS DESTROYED BY DET.
C IER IS AN ERROR RETURN FROM LUDCMP.
DIMENSION A(NA,1), IWORK(1), WORK(1)
CALL LUDCMP(A,N,NA,IWORK,DETT,WORK,IER)
DET=DETT
IF (IER.NE.0) THEN
DET=0.
RETURN
ENDIF
DO 11 J=1,N
DET=DET*A(J,J)
11 CONTINUE
RETURN
END
C**********************************************************************
SUBROUTINE LUDCMP(A,N,NP,INDX,D)
C MODIFICATION OF PRESS ET AL.
C LU DECOMPOSITION OF N BY N MATRIX A, STORED WITH NP ROWS.
C A IS DESTROYED BY THIS SUBROUTINE
C INDX AND VV ARE WORK VECTORS OF LENGTH N.
C IER=1 INDICATES A IS ALL ZEROES.
implicit double precision (a-h,o-z)
PARAMETER (TINY=1.0E-20, ndatmx = 100)
DIMENSION A(np,np),INDX(np),vv(ndatmx)
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.) pause 'singular matrix'
VV(I)=1./AAMAX
12 CONTINUE
DO 19 J=1,N
if (j.gt.1) then
DO 14 I=1,J-1
SUM=A(I,J)
if (i.gt.1) then
DO 13 K=1,I-1
SUM=SUM-A(I,K)*A(K,J)
13 CONTINUE
A(I,J)=SUM
endif
14 CONTINUE
endif
AAMAX=0.
DO 16 I=J,N
SUM=A(I,J)
if (j.gt.1)then
DO 15 K=1,J-1
SUM=SUM-A(I,K)*A(K,J)
15 CONTINUE
A(I,J)=SUM
endif
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
RETURN
END