Cut: Top of file. / Content. C. Math-Stat Subroutine Library "libmat.f". C FUNCTION GAMIN(A,X) C. Incomplete gamma-function. C. SUBROUTINE POLYG ( Z, GLN, PSI, PSI1, PSI2, PSI3 ) C. Gamma function and logarithmic derivatives. C SUBROUTINE GETAVU( ZM1, CM2, CM3, A, V, U, IFLAG ) C Finds gamma parameters from moments. C SUBROUTINE GFRAC ( INP, A, X, P, Q, PLN, QLN) C Finds fractiles of gamma distribution. C FUNCTION GFRACT ( A, P ) C Special function of GFRAC. C FUNCTION QLFRAC ( A, QLN ) C Special function of GFRAC. C FUNCTION GLFRAC ( A, PLN ) C Special function of GFRAC. C FUNCTION GAMLN ( X, ISIGN ) C Ln gamma x. C FUNCTION XOFGAM( GAM, IDOM ) C Inverse gamma-function. C FUNCTION MXGEX (MOD, CYCN, A, V, U, C 1 BUP, TUP, WUP, XMUP, EXUP, SIGXUP, SKWXUP, CHXUP, C 1 BDW, TDW, WDW, XMDW, EXDW, SIGXDW, SKWXDW, CHXDW) C. Fitting exponential gamma extreme value distribution to exponential C gamma distributed observations. C SUBROUTINE EXTGAM(A,H,C,ALN,X,B,V,U,XMEAN,STD,SKW) C Extreme value distribution of gamma distributed variable. C FUNCTION XM( A, H, AN ) C. Finds maximum point of extreme value distribution in EXTGAM. C. FUNCTION MRICE (I, EPS, X, F, P, Q) C. Organisation of evaluation of probability density and cumulative C. probability of the Rice distribution in different regimes. C. SUBROUTINE RICE (EPSLN, X, F, P, Q) C. Basic fourmulae for evaluation of probability density and cumulative C. probability of the Rice distribution. C. SUBROUTINE NORM ( I, X, F, P, Q, IER ) C. Normal probability, density, integral and inverse. C. FUNCTION NOREX (MOD, CYCN, B, T, W, C. 1 XIM, EXI, SXI, SKW, ALG, PSI, PSI1, PSI2, PSI3) C. Extreme value distribution parameters for normal distr. sample. C. FUNCTION LNOREX (MOD, AMU, SIGMA, CYCN, B, T, W, ZM, C. 1 ZM1, SIGZ, SKWZ, ALGB, PSIB, PSI1B, PSI2B, PSI3B) C. Extreme value parameters for the normalised log-normal distribution. C. C. FUNCTION BESJ (N, X, IER) Bessel function, first kind, Bn (x). C. FUNCTION BESY (N, X, IER) Bessel function, second kind, Yn (x). C. FUNCTION BESI (N, X, IER) Hyperbolic Bessel function, first kind, In (x). C. FUNCTION ATAN4 (Y, X) C. Calculates arcus tangent y/x by giveh y and x. Cut: File Contents. / Subroutine GAMIN. Cut:************************************************************** C FUNCTION GAMIN(A,X) C C CALCULATION OF THE INCOMPLETE GAMMA FUNCTION ( INTEGRAL FROM C 0 TO X ) BY MEANS OF A RELATED POWER SERIES.. C C A I PARAMETER OF THE GAMMA FUNCTION C X I UPPER INTEGRATION LIMIT C C PROGRAMMED BY SVERRE GRAN, AVD.24 DNV AUG.1977. C TOL = 1.0 / 10.0**8 N = 0 C = 1.0 SUM = 1.0 20 N = N + 1 AN = FLOAT( N ) C = C * X / (A+AN) SUM = SUM + C F = C - TOL*SUM IF( F ) 30, 30, 20 30 GAMIN = EXP(-X) * X**A * SUM / A RETURN END Cut:**************************************************************** SUBROUTINE POLYG ( Z, GLN, PSI, PSI1, PSI2, PSI3 ) C Poly-gamma functions of a given argument Z C Z = Input argument C GLN = ALOG of complete gamma function of Z C PSI = Di-gamma of Z C PSI1 = Tri-gamma of Z C PSI2 = Tetra-gamma of Z C PSI3 = Penta-gamma of Z NRED = 0 X = Z IF(Z. GT. 5.0 ) GO TO 10 X = 5.0 - Z NRED = X X = Z + NRED 10 X2 = X*X X3 = X2*X ALX = ALOG(X) GLN = 0.9189385332045 - X + ( X - 0.5 )*ALX + * (1./12.-(1./360.-(1./1260.-(1./1680.-1./(1188.*X2) ) * /X2) /X2) /X2) /X PSI = ALX - * (0.5+(1./12.-(1./120.-(1./252.-(1./240.-1./(132.*X2) ) * /X2) /X2) /X2) /X) /X PSI1 = (1.+(0.5+(1./6.-(1./30.-(1./42.-(1./30.-5./(66.*X2) ) * /X2) /X2) /X2) /X) /X) /X PSI2 = -(1.+(1.+(0.5-(1./6.-(1./6.-(3./10.-5./(6.*X2) ) * /X2) /X2) /X2) /X) /X) /X2 PSI3 = (2.+(3.+(2.-(1.-(4./3.-(3.-10./X2) * /X2) /X2) /X2) /X) /X) /X3 IF( NRED. EQ. 0 ) GO TO 99 DO 50 I = 1, NRED X = X - 1.0 GLN = GLN - ALOG( X ) PSI = PSI - 1.0 / X PSI1 = PSI1 + 1.0 / X**2 PSI2 = PSI2 - 2.0 / X**3 PSI3 = PSI3 + 6.0 / X**4 50 CONTINUE 99 RETURN END Cut:CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SUBROUTINE GETAVU( ZM1, CM2, CM3, A, V, U, IFLAG ) C C SUBROUTINE DESCRIPTION C C This routine evaluates distribution parameters of the C Exponential Gamma Distribution Function on the basis C of the three first statistical moments of a random C variable. C The probability density is C C f(a,v,u:y) = exp(-a*(y-u)/v - exp(-(y-u)/v ) ) / (Gam(a)*v) C C CALL PARAMETERS C C ZM1 = First order moment about zero, i.e.mean value C CM2 = Second order central moment, i.e. variance C CM3 = Third order central moment. C A,V,U = Distribution parameters a,v,u respectively. C IFLAG = 0 : Normal output C IFLAG = 1 : Small skewness, normal probability suggested C V and U are then standard dev. and mean respectively C IFLAG = 2 : Iteration did not converge for some reason C or skewness coefficient is greater than 2.0. The C method does not give results. C C Programmed by S.Gran, FDIV, DnV 1982 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Tolerance in LN(T) in iteration procedure, TOL TOL = 0.001 TMAX = 1.990 NMAX = 99 C Test for small variance IFLAG = 2 A = 0.0 V = 0.0 U = 0.0 IF(CM2. LE. 0.0) GO TO 999 C Skewness coefficient T T = CM3 / CM2**1.5 AT = ABS(T) C Test for too great skewness IF( AT. GT. TMAX ) GO TO 999 C Test for normal probability case. IFLAG = 1 U = ZM1 V = SQRT( V ) IF( AT. LT. TOL ) GO TO 999 C C Prepare to find exp.gamma-parameters IF( T. LT. 0.0 ) SIGN = -1.0 IF( T. GT. 0.0 ) SIGN = +1.0 ALT = ALOG(AT) C C Starting value for A and iteration start A = 0.5 + 1.0/T**2 N = 0 20 N = N + 1 C CALL POLYG( A, GLN, PSI, PSI1, PSI2, PSI3 ) F = ALOG(-PSI2) - 1.5*ALOG(PSI1) FDER = PSI3/PSI2 - 1.5*PSI2/PSI1 DF = F - ALT AA = A - DF/FDER C C See if iteration runs wild IF( N - NMAX ) 40, 40, 999 C C See if sufficient accuracy is obtained 40 A = AA DEV = ABS(DF) - TOL IF( DEV ) 50, 50, 20 C C Normal output 50 IFLAG = 0 V = CM2/PSI1 V = SIGN * SQRT( V ) U = ZM1 + PSI * V 999 RETURN END Cut:****************************************************** SUBROUTINE GFRAC ( INP, A, X, P, Q, PLN, QLN) C Subroutine for exaluation of Gamma-distribution. C INP = 1: X is given and P, Q and QLN=-Ln Q are calculated. C INP = 2: P is given, and X, Q and QLN=-Ln Q are calculated. C INP = 3: Q is given, and X, Q and QLN=-Ln Q are calculated. C INP = 4: PLN=-Ln P is given, and X, P and Q are calculated. * C INP = 5: QLN=-Ln Q is given, and X, P and Q are calculated. * C * interchanged 3/10-1995. C Sub-programs required are: C FUNCTION GFRACT ( A, P ) C FUNCTION QLFRAC ( A, QLN ) C FUNCTION GAMIN(A,X) C FUNCTION GAMLN ( X, ISIGN ) C FUNCTION GLFRAC ( A, PLN ) Ligger paa is og er beregnet C paa fraktiler med ekstremt smaa sannsynligheter. C $verre-87 C**************************************************************** GAMLA = GAMLN(A,ISIG) GAMA = EXP(GAMLA) GO TO (100, 200, 300, 400, 500), INP C X is input. Determine the probabilities. 100 P = GAMIN ( A, X) / GAMA Q = 1.0 - P PLN = 10.**30 QLN = 10.**30 IF(Q. GT. 0.0) QLN = -ALOG(Q) IF(P. GT. 0.0) PLN = -ALOG(P) GO TO 999 C Cumulative probability as input: 200 Q = 1.0 - P PLN = 10.**30 QLN = 10.**30 IF(P. GT. 0.0) PLN = -ALOG( P ) IF(Q. GT. 0.0) QLN = -ALOG( Q ) IF(QLN. LE. 10.0) X = GFRACT( A, P ) IF(QLN. GT. 10.0) X = QLFRAC( A, QLN) GO TO 999 C Exceedance probability Q as input: 300 P = 1.0 - Q PLN = 10.**30 QLN = 10.**30 IF(P. GT. 0.0) PLN = -ALOG( P ) IF(Q. GT. 0.0) QLN = -ALOG( Q ) IF(QLN. LE. 10.0) X = GFRACT( A, P ) IF(QLN. GT. 10.0) X = QLFRAC( A, QLN) GO TO 999 C -Log P as input: 400 IF(PLN. LT. 0.00001) GO TO 410 C Moderate probabilities P = EXP( -PLN ) Q = 1.0 - P QLN = -ALOG( Q ) IF(QLN. LE. 10.0) X = GFRACT( A, P ) IF(QLN. GT. 10.0) X = QLFRAC( A, QLN) GO TO 999 C Cumulative probability very close to 1.0. 410 Q = PLN QLN = -ALOG( Q ) P = 1.0 - Q PLN = -ALOG( P ) IF(QLN. LE. 10.0) X = GFRACT( A, P ) IF(QLN. GT. 10.0) X = QLFRAC( A, QLN) GO TO 999 C -Log Q as input: 500 Q = 0.0 IF(QLN. LT. 70.) Q = EXP(-QLN) P = 1.0 - Q IF(P. GT. 0.0) PLN = -ALOG( P ) IF(QLN. LE. 10.0) X = GFRACT( A, P ) IF(QLN. GT. 10.0) X = QLFRAC( A, QLN) 999 RETURN END Cut:************************************************************* C FUNCTION GFRACT ( A, P ) C C Routine for evaluation of the inverse incomplete gamma C function, or fractiles of the gamma distribution. C A is the parameter of the gamma function. C P is the cumulative probability for which the argument C (GFRACT) is to be calculated. C C DnV, sgra nov-83, Revidert -87 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IMAX = 100 TOL = 0.0001 C Special cases: GFRACT = 0.0 IF(P. LE. 0.0) GO TO 999 GFRACT = 10.0**20 IF(P. GE. 1.0) GO TO 999 GFRACT = -ALOG(1.0 - P) IF( A. EQ. 1.0 ) GO TO 999 C Select starting value for iteration: A1 = A - 1.0 GAMA = GAMLN(A,ISIG) GAMA = EXP(GAMA) PG = P * GAMA X = (P * GAMA * A)**(1.0/A) IF(A. GT. 1.0) X = A1 C Start iteration: DO 50 ITR = 1, IMAX PP = GAMIN( A, X ) FD = X**A1 * EXP(-X) DR = (PG - PP) / FD X = X + DR IF(X. LT. 0.0) X = ABS( X ) DR = ABS(DR/X) IF(DR. LE. TOL ) GO TO 99 50 CONTINUE 99 GFRACT = X 999 RETURN END Cut:********************************************************** FUNCTION QLFRAC ( A, QLN ) C Subprogram for evaluation of fractiles of the gamma distribution C with parameter A for particularly great argument with small C exceedance probability. C Input parameter QLN = - ALOG(Q) where Q is cumulative probability. C Sverre -87. C************************************************************* TOL1 = 0.000001 TOL2 = 0.00001 ALNA = ALOG( A ) GAMLA = GAMLN( A, ISIG ) GAMA = EXP( GAMLA ) C Choose starting value for the iteration X = QLN C Start iteration to find the percentage point X DO 300 J = 1, 100 C Calculate the series [1 + (a-1)/x + (a-1)(a-2)/x**2 + . . ] SERIES = 1.0 TERM = 1.0 TERML = 10.0**8 AN = 0.0 DO 210 I = 1, 100 AN = AN + 1.0 TERM = TERM * (A - AN) / X IF(TERM. GT. TERML) GO TO 220 C WRITE(6,*) I, TERM, SERIES SERIES = SERIES + TERM REL = ABS( TERM / SERIES ) TERML = TERM IF( REL. LT. TOL1) GO TO 220 210 CONTINUE C Calculate increment DX and test new X against desired accuracy. C T is current value for -ALOG(P) approaching PLN. C Compare T with input value pln 220 T = GAMLA - (A-1.0)*ALOG(X) + X IF(SERIES. LE. 0.0) GO TO 350 T = T - ALOG( SERIES ) TDER = 1.0 / SERIES DX = (QLN - T)/TDER C WRITE(6,*) A, T, DX, X X = X + DX REL = ABS( DX / X) IF(REL. LT. TOL2) GO TO 350 300 CONTINUE 350 QLFRAC = X RETURN END Cut:********************************************************** FUNCTION GLFRAC ( A, PLN ) C Subprogram for evaluation of fractiles of the gamma distribution C with parameter A. Improved version of FUNCTION GFRACT (A,P). C Input parameter PLN = - ALOG(P) where P is cumulative probability. C NB! RUTINA LIGGER PAA IS OG BRUKES IKKE. DELS FORDI DE ANDRE C RUTINENE DEKKER ALLE PRAKTISKE FORMAAL. DELS FORDI DEN C IKKE ER SKIKKELIG TESTET UT. C Sverre -87. C************************************************************* TOL1 = 0.0000001 TOL2 = 0.00001 ALNA = ALOG( A ) GAMLA = GAMLN( A, ISIG ) GAMA = EXP( GAMLA ) C Choose starting value for the iteration XSTAR = A*GAMA/EXP(PLN) IF(XSTAR. GT. 1.0. AND. A. GT. 1.0) XSTAR = XSTAR**(1.0/A) IF(XSTAR. LT. 1.0. AND. A. LT. 1.0) XSTAR = XSTAR**(1.0/A) X = 0.1 * XSTAR IF(PLN. LT. 0.7) X = 0.1 * A C Start iteration to find the percentage point X DO 300 J = 1, 100 C Calculate the series [1 + x/(a+1) + a**2/(a+1)(a+2) + . . ] SERIES = 1.0 TERM = 1.0 AN = 0.0 DO 210 I = 1, 100 AN = AN + 1.0 TERM = TERM * X / ( A + AN ) SERIES = SERIES + TERM REL = ABS( TERM / SERIES ) IF( REL. LT. TOL1) GO TO 220 210 CONTINUE C Calculate increment DX and test new X against desired accuracy. C T is current value for -ALOG(P) approaching PLN. C Compare T with input value pln 220 T = GAMLA + ALNA - A*ALOG(X) + X - ALOG( SERIES ) TDER = -A / (X * SERIES) DX = (PLN - T)/TDER X = X + DX REL = ABS( DX / X) IF(REL. LT. TOL2) GO TO 350 300 CONTINUE 350 GLFRAC = X RETURN END Cut:....................................................................... C FUNCTION GAMLN ( X, ISIGN ) C C THIS FUNCTION CALCULATES THE NATURAL LOGARITHM OF THE GAMMA FUNCTION C WITH A GIVEN ARGUMENT X. C IF THE GAMMA FUNCTION IS SINGULAR FOR X, THE OUTPUT VALUE IS SET C EQUAL TO A VERY LARGE NUMBER ( 10**20 ) AND THE SIGN-INDICATOR ISIGN C IS SET EQUAL TO 0. C IF THE SIGN OF THE GAMMA FUNCTION IS NEGATIVE, THE SIGN-INDICATOR C IS SET EQUAL TO -1 AND THE OUTPOT VALUE IS THE NATURAL LOG. OF THE C ABSOLUTE VALUE OF THE FUNCTION. C IF THE GAMMA FUNCTION IS POSITIVE, ISIGN=1 C C PROGRAMMED IN SEPTEMBER1977 BY SVERRE GRAN, AVD 24 DNV. C C .............................................................................. C PI = 3.141592654 PI2 = 2.0 * PI C C TEST FOR SINGULARITIES AT ZERO AND NEGATIVE INTEGERS IF( X ) 10, 20, 30 C C NAGATIVE X-VALUES 10 Y = - X IY = IFIX( Y ) AIY = FLOAT( IY ) DEC = Y - AIY IF( DEC ) 40, 20, 40 20 GAMLN = 10.0**20 ISIGN = 0 RETURN C C POSITIVEX. Y IS THE POSITIVE VALUE OF X WHICH ENTERS THE CALCULATIONS C IN THE FOLLOWING. IF Y IS SMALL, THE ACTUAL EVALUATION IS CARRIED OUT C FOR A NUMBER Z, AND THE VALUE FOR Y IS FOUND BY RECURSION 30 Y = X C C SEE IF X IS SO SMALL THAT RECURSION IS REQUIRED 40 YLIM = 10.0 DIFF = YLIM - Y IF( DIFF ) 50, 50, 60 50 NRECUR = 0 Z = Y GO TO 70 60 NRECUR = IFIX( DIFF ) RECUR = FLOAT( NRECUR ) Z = Y + RECUR 70 CONTINUE C C ACTUAL EVALUATION BY AN ASYMPTOTIC FORMULA ZZ = Z * Z GAMLN = (((( 5.0/4950.0/ZZ - 1.0/1680.0)/ZZ + 1.0/1260.0)/ZZ * - 1.0/360.0) / ZZ + 1.0/12.0)/Z * + (Z - 0.5)*ALOG( Z ) - Z + 0.5 * ALOG( PI2 ) C C SEE IF RECURSION HAS BEEN REQUIRED 80 IF( NRECUR ) 100, 100, 90 90 NRECUR = NRECUR - 1 RECUR = RECUR - 1.0 ZREC = Y + RECUR GAMLN = GAMLN - ALOG( ZREC ) GO TO 80 C C SEE IF X WAS NEGATIVE 100 IF( X ) 120, 110, 110 110 ISIGN = 1.0 RETURN C 120 PIY = PI * Y SP = SIN( PIY ) ASP = ABS( SP ) SIGN = -SP / ASP ISIGN = IFIX( SIGN ) GAMLN = ALOG( PI ) - ALOG( Y ) - ALOG( ASP ) - GAMLN RETURN END Cut:****************************************************************** FUNCTION XOFGAM( GAM, IDOM ) C Purpose: To evaluate the inverse gamma function for positive arguments. C GAM is the value of the gamma function. C IDOM = 1 for x-values ABOVE the minimum point at 1.46. C =-1 for x-values BELOW the minimum point at 1.46. C XOFGAM is returned as the value which gives the gamma C function equal to GAM. C If GAM is less than 0.8856, XOFGAM = 0.0 by return. C Coded by Sverre in April 86. C*************************************************************** TOL = 0.0001 C Test for value below the allowable minimum X = 0.0 IF(GAM. LT. 0.8856) GO TO 99 C Use asymptotic expressions by extreme values X = 1. / GAM IF(IDOM. LT. 0. AND. GAM. GT. 1000.) GO TO 99 X = ALOG( GAM ) IF(IDOM. GT. 0. AND. GAM. GT. 1000000.) GO TO 99 C Relevant starting values and iteration process in the moderate range IF(IDOM. LT. 0) X = 0.001 IF(IDOM. GT. 0) X = 20. 20 CALL POLYG ( X, GLN, PSI, PSI1, PSI2, PSI3 ) Y = EXP(GLN) ERR = ABS( Y - GAM ) / Y IF(ERR. LT. TOL) GO TO 99 X = X + (GAM - Y) / (Y * PSI) GO TO 20 99 XOFGAM = X RETURN END Cut: Function XOFGAM. / Function MXGEX: FUNCTION MXGEX (MOD, CYCN, A, V, U, 1 BUP, TUP, WUP, XMUP, EXUP, SIGXUP, SKWXUP, CHXUP, 1 BDW, TDW, WDW, XMDW, EXDW, SIGXDW, SKWXDW, CHXDW) C. PURPOSE: C. Fitting an extreme value distribution of the exponential Gamma C. kind to the extreme of exponential gamma distributed observations. C. CALL PARAMETERS: C. MOD I Distribution types, 1 for general exp. gamma. 2 for Gumbel. C. CYCN I Number of observations. C. A,V,U I Exponential gamma parameters. C. Extension *UP is for the largest maximum value. C. Extension *EW is for the smallest minimum value. C. B* O Shape parameter. C. T* O Slope parameter. C. W* O Mode parameter. C. XM* O Mode of distribution, most probable largest/smallest value. C. EX* O Expected largest/smallest. C. SIGX* O Standard deviation of extreme value distribution. C. SKWX* O Skewness coefficient. C. CHX* O Characteristic largest/smallest value. C. MXGEX O Status: 0 - No results. 1-Exponential gamma distribution. C. -1 - Left skewed exponential gamma distribution (v<0). C. +1 - Right skewed exponential gamma distribution (v>0). C. SUBPROGRAMS REQUIRED: C. SUBROUTINE POLYG Poly-gamma functions. C. FUNCTION GFRAC Inverse gamma distribution. C Coding: $verre -87. Web-version -96. MXGEX = 0 IF (CYCN. LE. 0.0) GO TO 990 IF (A.LE.0.0. OR. V.EQ.0.0) GO TO 990 TOL = 0.00001 ALNN = ALOG (CYCN) C. Characteristic largest and smallest values: PROC = 1.0 / CYCN IF (V. LT. 0.0) THEN CALL GFRAC (3, A, ZUP, P, PROC, PLN, QLN) CALL GFRAC (2, A, ZDW, PROC, Q, PLN, QLN) ELSE IF (V. GT. 0.0) THEN CALL GFRAC (2, A, ZUP, PROC, P, QLN, PLN) CALL GFRAC (3, A, ZDW, Q, PROC, QLN, PLN) END IF CHXUP = U - V * ALOG (ZUP) CHXDW = U - V * ALOG (ZDW) C. Provide the poly-gamma functions of A: CALL POLYG ( A, ALGA, PSI, PSI1, PSI2, PSI3 ) GAMA = EXP (ALGA) C. Derive the exponential gamma extreme for X directly. Section 4.3.2 (v). C. Solve for the most probable largest x-value by (4.3.67-68). C. Carry out the calculations under the condition that v>0. C. In case v<0 place the distribution as the extreme small distribution. VA = ABS (V) X0 = VA * ( ALNN - ALOG (GAMA) ) / A X = X0 Y = EXP ( - X / VA ) DO 240 I = 1, 50 XOLD = X DUMS = ABS (A - Y) X = X0 - Y*VA/A - VA*ALOG ( DUMS ) / A Y = EXP ( - X / VA ) ERR = X - XOLD ERR = ABS ( ERR ) IF (ERR. L T. TOL) GO TO 250 240 CONTINUE 250 XM = X + U YM = Y G2 = ( (YM - A)**2 + YM) / V**2 G3 = ( (YM - A)*((YM - A)**2 - YM) - YM ) / V**3 B = G2**3 / G3**2 CALL POLYG (B, BLG, PSI, PSI1, PSI2, PSI3 ) GAMB = EXP ( BLG ) C. If v<0 this describes the lower tail extreme value. IF (V. LT. 0.0) THEN XMDW = 2.0 * U - XM BDW = B TDW = -G2 / G3 WDW = XMDW + TDW * ALOG (B) ABST = ABS (TDW) EXDW = WDW - TDW * PSI SIGXDW = ABST * SQRT( PSI1 ) SKWXDW = PSI2 / PSI1**1.5 MXGEX = -1 C. If v>0 this describes the upper tail extreme value. ELSE IF (V. GT. 0.0) THEN XMUP = XM BUP = B TUP = -G2 / G3 WUP = XMUP + TUP * ALOG (B) ABST = ABS (TUP) EXUP = WUP - TUP * PSI SIGXUP = ABST * SQRT( PSI1 ) SKWXUP = -PSI2 / PSI1**1.5 MXGEX = 1 END IF C. Solve for the most probable largest x-value by (4.3.69-70). C. Carry out the calculations under the condition that v<0. C. In case v<0 place the distribution as the extreme small distribution. C. In case v>0 place the distribution as the extreme LARGE distribution. VA = ABS (V) Y0 = ALNN - ALOG (GAMA) Y = Y0 X = VA * ALOG ( Y0 ) DO 280 I = 1, 50 XOLD = X DUMS = Y**A / (Y - A) DUMS = Y0 + ALOG (DUMS) X = VA * ALOG( DUMS ) Y = EXP ( X / VA ) ERR = X - XOLD ERR = ABS ( ERR ) IF (ERR. L T. TOL) GO TO 290 280 CONTINUE 290 XM = X + U YM = Y G2 = ( (YM - A)**2 + YM) / V**2 G3 = ( (YM - A)*((YM - A)**2 - YM) - YM ) / V**3 B = G2**3 / G3**2 CALL POLYG (B, BLG, PSI, PSI1, PSI2, PSI3 ) GAMB = EXP ( BLG ) C. If v<0 this describes the upper tail extreme value. IF (V. LT. 0.0) THEN XMUP = XM BUP = B TUP = -G2 / G3 WUP = XMUP + TUP * ALOG (B) ABST = ABS (TUP) EXUP = WUP - TUP * PSI SIGXUP = ABST * SQRT( PSI1 ) SKWXUP = -PSI2 / PSI1**1.5 MXGEX = -1 C. If v>0 this describes the lower tail extreme value. ELSE IF (V. GT. 0.0) THEN XMDW = 2.0 * U - XM BDW = B TDW = -G2 / G3 WDW = XMDW + TDW * ALOG (B) ABST = ABS (TDW) EXDW = WDW - TDW * PSI SIGXDW = ABST * SQRT( PSI1 ) SKWXDW = PSI2 / PSI1**1.5 MXGEX = 1 END IF 990 RETURN END Cut: Function MXGEX. / Subroutine EXTGAM Cut:************************************************************** C SUBROUTINE EXTGAM(A,H,C,ALN,X,B,V,U,XMEAN,STD,SKW) C C Subroutine which evaluates the extreme value distribution C of a gamma-distributed variable. C The extreme value distribution is given as an exponential C general gamma distribution. C C Input variables: C A,H,C Shape-, Slope- and Scale-parameter of the gamma C distribution of the basic variable. C ALN Log base 10 of the number of samples or cycles. C C Output variables: C X Most probable largest value C B,V,U Exp.general gamma parameter of the extreme C XMEAN Expected largest value C STD Standard deviation of the extreme C SKW Skewness coefficient of the extreme value dist. C C Subprograms required: C SUBROUTINE POLYG C FUNCTION XM C C Reference: C Lectures in Ocean Engineering. Article 6.6 C C Coded by S.Gran Des.1984 C Error exit for X=0.0 done 3/10-1995. C *************************************************************** AN = 10.0**ALN X = XM( A, H, AN ) IF (X. EQ. 0.0) GO TO 95 D1 = (A*H - 1.) / X - H * X**(H-1.) D2 = -(A*H - 1.)/ X**2 - H*(H-1.) * X**(H-2.) D3 = 2.*(A*H - 1.)/ X**3 - H*(H-1.)*(H-2.)* X**(H-3.) G2 = D1**2 - D2 G3 = D1**3 + D1*D2 - D3 B = G2**3 / G3**2 V = - C * G2 / G3 U = C * X + V * ALOG( B ) CALL POLYG(B, GLN, PSI, PSI1, PSI2, PSI3 ) XMEAN = U - (V * PSI) STD = V * SQRT( PSI1 ) SKW = - PSI2 / PSI1**1.5 X = C * X 95 RETURN END Cut:******************************************************** FUNCTION XM( A, H, AN ) C Function for evaluation of the most probable C largest value in a general gamma distribution C by successive iterations C Error exit if XM = 0.0. C Start values put in 3/10-1995. XM = 0.0 IF (A.LE.0.0. OR. H.EQ.0.0. OR. AN.LE.0.0) GO TO 60 CALL POLYG ( A, GLN, PSI, PSI1, PSI2, PSI3 ) FN = AN / EXP( GLN ) FNLN = ALOG (FN) AH = A*H D = (AH - 1.0)/H IF (H.GT.0.0) XM = FNLN**(1.0/H) IF (H.LT.0.0) THEN IF (AH. GT. -1.50) GO TO 60 XM = AN**(-1.0/AH) END IF DO 50 I = 1, 10 XLAST = XM IF(H. GT. 0.0) THEN XM = FN * XM**AH / (XM**H - D) XM = ALOG(XM)**(1.0/H) ELSE IF(H. LT. 0.0) THEN C. First equation below is alternative to (4.3.48) which is applied: C XM = (FN*EXP(-XM**H) + XM**(-AH+H)) / D XM = FN * EXP(-XM**H) / (D - XM**H) XM = XM**(-1/(A*H)) END IF DIFF = ABS((XM-XLAST)/XM) IF( DIFF. LT. 0.00001 ) GO TO 60 50 CONTINUE XM = 0.0 60 RETURN END Cut: Function XM. / Function MRICE. FUNCTION MRICE (I, EPS, X, F, P, Q) C. PURPOSE: C. Organisation of evaluation of probability density and cumulative C probability of the Rice distribution in different regimes. C. CALL PARAMETERS: C. I = Mode of operation: C. = 1: X given, find F, P and Q. C. = 3: P given, find F, X and Q. C. = 4: Q given, find F, X and P. C. EPS = Skewness parameter C. F = Probability density C. P = Cumulative probability C. Q = Exceedance probability Ck MRICE = Status: 0 = no response. Otherwise MRICE = I. C Coded by $g July 1996. C. Rice distribution by given argument X: MRICE = 0 IF (I. EQ. 1) THEN CALL RICE (EPS, X, F, P, Q) MRICE = 1 GO TO 99 END IF X = 0.0 F = 0.0 PI = 3.141592654 SQ2PI = 2.506628275 C. Prepare for evaluation of percentage points: IF (EPS. LT. 0.0) EPS = 0.0 IF (EPS. GT. 1.0) EPS = 1.0 IF (I. EQ. 3) THEN MRICE = 3 IF (P. LT. 0.0) P = 0.0 IF (P. GT. 1.0) P = 1.0 Q = 1.0 - P ELSE IF (I. EQ. 4) THEN MRICE = 4 IF (Q. GT. 1.0) Q = 1.0 IF (Q. LT. 0.0) Q = 0.0 P = 1.0 - Q END IF IF (MRICE. EQ. 0) GO TO 99 C. Calculate Normal case for Epsilon close to 1.0: EPSBIG = 0.9990 CALL NORM (I, XNOR, FNOR, P, Q, IER) IF (EPS. GT. EPSBIG) THEN X = XNOR F = FNOR GO TO 99 END IF C. Calculate Rayleigh case for Epsilon close to 0.0: EPSLOW = 0.0010 PQMIN = 1E-12 IF (Q.GT.0.0) XMAX2 = 2.0 * ALOG (1.0/Q) IF (Q.LE.0.0) XMAX2 = 2.0 * ALOG (1.0/PQMIN) XMAX = SQRT (XMAX2) IF (EPS. LT. EPSLOW) THEN X = XMAX F = X * Q GO TO 99 END IF C. Identify minimum acceptable X-value, XMIN: XMIN = -12.0 C. But the following curve to fix lower values of X is not good. Check later. IF (EPS.GE.0.001. AND. EPS.LT.0.005) THEN XMIN = -500.0 * EPS**1.92 ELSE IF (EPS.GE.0.005. AND. EPS.LT.0.5) THEN XMIN = -11.0 * EPS**1.2 ELSE IF (EPS.GT.0.5. AND. EPS.LT.1.0) THEN XMIN = 1.0 - 12.0*EPS END IF C. Exit for probabilities close to 0 or 1: F = 0.0 C. Low cumulative probability, Return the minimum value: IF (P. LE. PQMIN) THEN X = XMIN GO TO 99 ELSE IF (Q. LE. PQMIN) THEN C. Low exceedance value. Return a balanced asymptotic average C. between Rayleigh and normal: XMID2 = XMAX2 - (1.0 - EPS**2) * ALOG (SQ2PI * XMAX) X = SQRT (XMID2) GO TO 99 END IF C. Iterate for intermediate Epsilon-values. AC is P(0). F0 is F(0). C. X1 is first approximation to X and has the correct sign: AC = 0.5 * (1.0 - SQRT (1.0 - EPS**2)) F0 = EPS / SQ2PI X1 = (P - AC) / F0 IF (AC. GT. 0.0) PL0 = - ALOG (AC) FL0 = -F0 / PL0 C. Return if P has the theoretical AC value of P(0): IF (P. EQ. AC) THEN X = X1 F = F0 GO TO 99 C. If P<=0, then iterate along PL=-ln P, Starting at X=0.0: ELSE IF (P. LE. 0.5) THEN X = 0.0 TOLSP = 0.00001 CALL RICE (EPS, X, F, PX, QX) J = 0 DELX = ALOG (P / PX) * PX / F DO 10 J = 1, 10 DELX = ALOG (P / PX) * PX / F X = X + DELX CALL RICE (EPS, X, F, PX, QX) S = SQRT (DELX**2 + (PX-P)**2) IF (X. LT. XMIN) THEN X = XMIN GO TO 99 END IF IF (S. LT. TOLSP) GO TO 99 10 CONTINUE C. If 0.5<P then interate along QL =-ln Q, Starting at the Rayleigh value: ELSE IF (P. GT. 0.5) THEN X = XMAX TOLSQ = 0.00001 CALL RICE (EPS, X, F, PX, QX) J = 0 DELX = ALOG (QX / Q) * QX / F DO 30 J = 1, 10 DELX = ALOG (QX / Q) * QX / F X = X + DELX CALL RICE (EPS, X, F, PX, QX) S = SQRT (DELX**2 + (PX-P)**2) IF (S. LT. TOLSQ) GO TO 99 30 CONTINUE END IF C. The following operation iterates along P starting at zero C and may be used in the medium P-range for some future purpose. C IF (X1. LE. 0.0) THEN C TOLSP = 0.00001 C X = 0.0 C PX = AC C F = F0 C WRITE (6, '('' AC, F:'',2F12.7)') AC, F C DO 10 J = 1, 20 C DELX = (P - PX) / F C X = X + DELX C CALL RICE (EPS, X, F, PX, QX) C S = SQRT (DELX**2 + (PX-P)**2) C WRITE(6,'(''J,Delx,X,F,Px,S:'',I4,5F12.7)')J,DELX,X,F,PX,S C IF (S. LT. TOLSP) GO TO 99 C 10 CONTINUE 99 RETURN END Cut: Function MRICE. / Subroutine Rice. SUBROUTINE RICE (EPSLN, X, F, P, Q) C Basic fourmulae for evaluation of probability density and cumulative C probability of the Rice distribution. C EPSLN = Skewness parameter C F = Probability density C P = Cumulative probability C Q = Exceedance probability. C For epsilon < 0.00100 the normal distribution is used. C For epsilon > 0.99990 the Rayleigh distribution is used. C $verre -88, Fixed up in -96. Cline: Standard variables and values EMAX = 12.0 TOLLOW = 0.0010 TOLHIG = 0.99990 PI = 3.141592654 SQ2PI = 2.506628275 Cline: Test formulae to be applied, depending of epsilon: EPS = EPSLN IF(EPSLN. LT. 0.0) EPS = 0.0 IF(EPSLN. GT. 1.0) EPS = 1.0 IF(EPS. LT. TOLLOW) GO TO 130 IF(EPS. GT. TOLHIG) GO TO 160 Cline: General Rice case with specific value of epsilon SQEPS = SQRT(1.0 - EPS**2) AR1 = X * SQEPS / EPS AR2 = X / EPS Cline: Normal probability integral of first argument AR1 XX = ABS (AR1) IF (XX. LE. EMAX) THEN EX = EXP (-0.5*XX**2) Z = 1.0 / ( 1.0 + 0.2316419 * XX ) PH1 = Z * ( 0.127414796 + Z * (-0.142248368 + Z * 1 ( 0.7107068705 + Z * (-0.7265760135 + Z * 1 0.5307027145 ) ) ) ) * EX ELSE IF (XX. GT. EMAX) THEN PH1 = 0.0 END IF PHI1C = PH1 PHI1 = 1.0 - PH1 IF(AR1. LT. 0.) PHI1 = PH1 IF(AR1. LT. 0.) PHI1C = 1.0 - PH1 Cline: Normal probability integral of second argument AR2 XX = ABS (AR2) IF (XX. LE. EMAX) THEN EX = EXP (-0.5*XX**2) Z = 1.0 / ( 1.0 + 0.2316419 * XX ) PH2 = Z * ( 0.127414796 + Z * (-0.142248368 + Z * 1 ( 0.7107068705 + Z * (-0.7265760135 + Z * 1 0.5307027145 ) ) ) ) * EX ELSE IF (XX. GT. EMAX) THEN PH2 = 0.0 END IF PHI2C = PH2 PHI2 = 1.0 - PH2 IF(AR2. LT. 0.) PHI2 = PH2 IF(AR2. LT. 0.) PHI2C = 1.0 - PH2 Cline: Exponentials XX1 = ABS (X) EXP1 = -0.5 * X**2 IF (XX1. LE. EMAX) THEN EXP1 = EXP (EXP1) ELSE IF (XX1. GT. EMAX) THEN EXP1 = 0.0 END IF XX2 = ABS (X/EPS) EXP2 = -0.5 * XX2**2 IF (XX2. LE. EMAX) THEN EXP2 = EXP (EXP2) ELSE IF (XX2. GT. EMAX) THEN EXP2 = 0.0 END IF TERM2 = SQEPS * PHI1 * EXP1 F = EPS * EXP2 / SQ2PI + TERM2 * X IF(X. LT. 0.0) THEN P = PHI2 - TERM2 IF (P. LT. 0.0) P = 0.0 Q = 1.0 - P END IF IF(X. GE. 0.0) THEN Q = PHI2C + TERM2 IF (Q. LT. 0.0) Q = 0.0 P = 1.0 - Q END IF IF (P. EQ. 0.0. AND. Q. EQ. 1.0) F = 0.0 IF (P. EQ. 1.0. AND. Q. EQ. 0.0) F = 0.0 IF (F. LT. 0.0) F = 0.0 GO TO 999 Cline: Epsilon close to 0.0, use the Rayleigh distribution. 130 IF (X. LE. 0.0) THEN Q = 1.0 ELSE IF (X.GT.0.0. AND. X.LT.EMAX) THEN Q = EXP (-0.5*X**2) ELSE IF (X. GT. EMAX) THEN Q = 0.0 END IF P = 1.0 - Q F = X * Q GO TO 999 Cline: Epsilon close to 1.0, use the normal distribution. 160 XX = ABS (X) IF (XX. LE. EMAX) THEN EX = EXP (-0.5*XX**2) F = EX / SQ2PI Z = 1.0 / ( 1.0 + 0.2316419 * XX ) QQ = Z * ( 0.127414796 + Z * (-0.142248368 + Z * 1 ( 0.7107068705 + Z * (-0.7265760135 + Z * 1 0.5307027145 ) ) ) ) * EX ELSE IF (XX. GT. EMAX) THEN F = 0.0 QQ = 0.0 END IF IF (X. LE. 0.0) THEN P = QQ Q = 1.0 - QQ ELSE IF (X. GT. 0.0) THEN Q = QQ P = 1.0 - QQ END IF Cline: Stop and return 999 RETURN END Cut: Subroutine RICE. / Subroutine NORM. C************************************************************** SUBROUTINE NORM ( I, X, F, P, Q, IER ) C C EVALUATION OF CORRESPONDING VALUES AND ARGUMENTS OF THE NORMAL C PROBABILITY FUNCTION BY APPROXIMATION WITH RATIONAL FUNCTIONS. C CORRESPONDING VALUES ARE GIVEN FOR ALL CALL PARAMETERS. C C I = 1 GIVEN ARGUMENT X. C = 2 GIVEN PROBABILITY DENSITY F. C = 3 GIVEN CUMULATIVE PROBABILITY P. C = 4 GIVEN EXCEEDANCE PROBABILITY Q. C C IER = 0 NORMAL EXIT. C = 1 X HAS TOO LARGE POSITIVE VALUE. C = 2 X HAS TOO SMALL NEGATIVE VALUE. C = 3 INPUT VALUE IS OUT OF RANGE OF THE LOW P SIDE. C = 4 INPUT IS OUT OF RANGE ON THE LOW Q SIDE. C = 5 INPUT DENSITY IS LARGER THAN 1/SQRT(2*PI) C C PROGRAMMED BY SVERRE GRAN, DNV/FDIV/24, NOVEMBER 1978. C ...................................................................... C IER = 0 DEFLTX = 100.0 XMAX = 13.1 XMIN = -XMAX SQ2PI = 2.506628275 SQINV = 0.398942280 C GO TO ( 100, 200, 300, 400 ), I C C GIVEN X : C ********* C 100 IF( X. GT. XMAX ) GO TO 701 IF( X. LT. XMIN ) GO TO 601 IF( X ) 105, 501, 105 105 XX = -0.5 * X * X EX = EXP( XX ) F = EX / SQ2PI XX = ABS( X ) Z = 1.0 / ( 1.0 + 0.2316419 * XX ) QQ = Z * ( 0.127414796 + Z * (-0.142248368 + Z * * ( 0.7107068705 + Z * (-0.7265760135 + Z * * 0.5307027145 ) ) ) ) * EX IF( X ) 135, 135, 140 135 P = QQ Q = 1.0 - QQ RETURN 140 Q = QQ P = 1.0 - QQ RETURN C C GIVEN F : C ********* C 200 IF( F - SQINV ) 220, 500, 210 210 IER = 5 GO TO 500 220 IF( F ) 700, 700, 230 230 XX = F * SQ2PI XX = -2.0 * XX X = SQRT( XX ) GO TO 100 C C GIVEN P : C ********* C 300 IF( P ) 600, 600, 305 305 IF( 1.0 - P ) 700, 700, 310 310 IF( P - 0.5 ) 320, 500, 330 320 QQ = P GO TO 800 330 QQ = 1.0 - P GO TO 800 C C GIVEN Q : C ********* C 400 IF( Q ) 700, 700, 405 405 IF( 1.0 - Q ) 600, 600, 410 410 IF( Q - 0.5 ) 420, 500, 430 420 QQ = Q GO TO 800 430 QQ = 1.0 - Q GO TO 800 C 500 X = 0.0 501 F = SQINV P = 0.5 Q = 0.5 RETURN 600 IER = 3 X = -DEFLTX 601 F = 0.0 P = 0.0 Q = 1.0 RETURN 700 IER = 4 X = DEFLTX 701 F = 0.0 P = 1.0 Q = 0.0 RETURN 800 T = -2.0 * ALOG( QQ ) T = SQRT(T) A = 2.515517 + T * ( 0.802853 + T * 0.010328 ) B = 1 + T*( 1.432788 + T*( 0.189269 + T*0.001308 ) ) X = T - A/B XX = -0.5 * X * X F = EXP( XX ) / SQ2PI 850 IF( I - 3 ) 870, 855, 870 855 Q = 1.0 - P IF( P - 0.5 ) 860, 500, 865 860 X = -X 865 RETURN 870 P = 1.0 - Q IF( Q - 0.5 ) 885, 500, 880 880 X = -X 885 RETURN END Cut: Subroutine NORM. / Function NOREX. FUNCTION NOREX (MOD, CYCN, B, T, W, 1 XIM, EXI, SXI, SKW, ALG, PSI, PSI1, PSI2, PSI3) C PURPOSE: C. The routine calculates extreme value distribution parameters for C. the normalised normal distribution in a sample of size CYCN. C. INPUT DESCRIPTION: C. MOD Type of extreme value distribution: C. = 1 Generalised exponential distribution, parameters B,T,W. C. = 2 Double exponential distribution, parameter B=1. C. CYCN Size of sample. C. B,T,W Shape, Slope and scale parameter. C. XIM Mode of normalised variable (most probable largest). C. EXI Expected largest. C. SXI Standard deviation of normalised extreme. C. SKW Skewness of extreme value distribution. C. ALG Natural log to gamma (B). C. PSI-PSI3 Polygamma functions of argument B. C. NOREX = 0 Error output. C. 1 Generalised exponential gamma extreme value distribution. C. 2 Double exponential extreme value distribution. NOREX = 0 IF (CYCN. LT. 0.0) CYCN = -CYCN IF (CYCN. EQ. 0.0) GO TO 990 TOL = 0.00001 ALNN = ALOG (CYCN) ALNN2 = 2.0 * ALNN ALGN = ALGN / 2.3025851 SQ2PI = 2.5066282 PI = 3.141592 C. Find most probable normalised value X and distribution parameters: C. - General exponential gamma: IF (MOD. EQ. 1) THEN XIM = 1.0 DO 40 I = 1,50 XIMOLD = XIM XIM = ALNN2 - 2.0*ALOG(SQ2PI*XIM) XIM = SQRT(XIM) ERR = XIM - XIMOLD ERR = ABS(ERR) IF(ERR. LT. TOL) GO TO 50 40 CONTINUE 50 B = (XIM**2 + 1.0)**3 / ( XIM**2 * (XIM**2 - 1.0)**2 ) CALL POLYG (B, ALG, PSI, PSI1, PSI2, PSI3) T = (XIM**2 + 1.0) / (XIM * (XIM**2 - 1.0)) W = XIM + T * ALOG(B) NOREX = 1 C. - Double exponential: ELSE IF (MOD. EQ. 2) THEN XIM = ALOG (4.0 * PI * ALNN) XIM = SQRT(2.0 * ALNN - XIM) B = 1.0 CALL POLYG (B, ALG, PSI, PSI1, PSI2, PSI3) T = 1.0 / XIM W = XIM NOREX = 2 END IF SXI = T * SQRT (PSI1) EXI = W - SXI * PSI SKW = -PSI2 / (PSI1**1.5) 990 RETURN END Cut: Function NOREX. / Function LNOREX. FUNCTION LNOREX (MOD, AMU, SIGMA, CYCN, B, T, W, ZM, 1 ZM1, SIGZ, SKWZ, ALGB, PSIB, PSI1B, PSI2B, PSI3B) C PURPOSE: C. The routine calculates extreme value distribution parameters for C. the normalised log-normal distribution in a sample of size CYCN. C. Kind of extreme value distribution is a generalised exponential C. gamma distribution. See Section 4.3.2(vi). C. INPUT DESCRIPTION: C. MOD 1: Generalised Exponential Gamma. C. 2: General Gamma. C. AMU Logarithmic mean value, - distribution parameter. C. SIGMA Logarithmic standard deviation, - dist. parameter. C. CYCN Size of sample. C. B,T,W Shape, Slope and scale parameter. C. ZM Mode of normalised variable (most probable largest). C. ZM1 Expected largest. C. SIGZ Standard deviation of normalised extreme. C. SKWZ Skewness of extreme value distribution. C. ALGB Natural log to gamma (B). C. PSIB-PSI3B Polygamma functions of argument B. C. LNOREX 0: Error output. C. 1: Generalised exponential gamma extreme value distribution. LNOREX = 0 IF (CYCN. LT. 0.0) CYCN = -CYCN IF (CYCN. EQ. 0.0) GO TO 990 TOL = 0.00001 ALNN = ALOG (CYCN) ALNN2 = 2.0 * ALNN ALGN = ALGN / 2.3025851 SQ2PI = 2.5066282 PI = 3.141592 C. Derive the exponential gamma extreme for X directly. Section 4.3.2 (vi). C. Solve for the most probable largest x-value, through the C. auxiliary variable z. See (4.3.82) to (4.3.84). IF (MOD. EQ. 1) THEN ZETA = 2.0 * ( ALNN - ALOG ( SQ2PI ) ) ZETA = ABS (ZETA) ZETA = SQRT (ZETA) DO 240 I = 1,50 ZOLD = ZETA DUMS = SQ2 PI * (ZETA + SIGMA) ZETA = 2.0 * (ALNN - ALOG (DUMS)) ZETA = SQRT (ZETA) ERR = ZETA - ZOLD ERR = ABS ( ERR ) IF (ERR. LT. TOL) GO TO 250 240 CONTINUE C. Mode found. Calculate the parameters: 250 ALNZ = AMU + SIGMA * ZETA ZM = EXP (ALNZ) DLP1 = -(1.0 - AMU/SIGMA**2) / ZM - ALNZ /(SIGMA**2 * ZM) DLP2 = (1.0 - AMU/SIGMA**2 - 1.0/SIGMA**2) / ZM**2 * + ALNZ /(SIGMA*ZM)**2 DLP3 = -(2.0 - 2.0*AMU/SIGMA**2 - 3.0/SIGMA**2) / ZM**3 * - 2.0*ALNZ /(SIGMA**2 * ZM**3) G2 = DLP1**2 - DLP2 G3 = DLP1 * (DLP1**2 + DLP2) - DLP3 B = G2**3 / G3**2 T = -G2 / G3 W = ZM + T * ALOG (B) LNOREX = 1 C. Simplified, Generalised Gumdel distribution, (4.3.23-26): ELSE IF (MOD. EQ. 2) THEN QQ = 1.0 / CYCN CALL NORM (4, XI, DENS, P, QQ, IER1) B = 1.0 ZM = EXP (AMU + SIGMA * XI) DENS = DENS / (ZM * SIGMA) T = QQ / DENS W = ZM LNOREX = 2 END IF C. Expectation value ZM1, standard dev. SIGZ CALL POLYG (B, BLG, PSI, PSI1, PSI2, PSI3 ) GAMB = EXP (BLG) ZM1 = W - T * PSI SIGZ = T * SQRT (PSI1) SKWZ = -PSI2 / PSI1**1.5 990 RETURN END Cut: Function LNOREX. / Function BESJ. FUNCTION BESJ (N, X, IER) C C FUNCTION BESJ C C PURPOSE C TO CALCULATE THE BESSEL FUNCTION OF FIRST KIND AND C INTEGER ORDER C C USAGE C BESJ (N, X, IER) C C N - ORDER OF THE BESSEL FUNCTION C X - ARGUMENT C C SUBPROGRAMS REQUIRED C NONE C C INTRINSIC ARRAYS C NONE C C METHOD C THE FUNCTION IS A SHORT HAND VERSION OF THE SUBROUTINE C BESJ DESCRIBED IN THE IBM-SYSTEM 360 SCIENTIFIC SUBROUTINE C PACKAGE (360A-CM-03X) VERSION II C C ............................................................. C C C COMPUTES THE BESSEL FUNCTION OF THE FIRST KIND FOR C GIVEN ORDER AND ARGUMENT C C THE FOLLOWING PARAMETER D IS THE REQUIRED ACCURACY C D=0.00005 C BESJ = 1.0 IF(N. EQ. 0. AND. X. EQ. 0.) RETURN BESJ = 0.0 ABSX = ABS (X) IF(N. GT. 0. AND. ABSX. LT. D) RETURN IF (N. LT. 0) THEN IER = 1 RETURN END IF C 31 IF(X-15.0) 32,32,34 32 NTEST = 20.0 + 10.0*X - X** 2/3 GO TO 36 34 NTEST = 90.0 + 0.5*X 36 IF(N-NTEST)40,38,38 38 IER = 4 RETURN C 40 IER = 0 N1 = N+1 BPREV = 0.0 C C COMPUTE STARTING VALUE OF M C IF(X-5.0) 50,60,60 50 MA = X + 6.0 GO TO 70 60 MA = 1.4*X+60.0/X 70 MB = N + IFIX(X)/4 + 2 MZERO = MAX0(MA,MB) C C SET UPPER LIMIT OF M C MMAX = NTEST 100 DO 190 M=MZERO,MMAX,3 C C SET F(M),F(M-1) C FM1 = 1.0E-28 FM = 0.0 ALPHA = 0.0 IF(M-(M/2)*2) 120,110,120 110 JT = -1 GO TO 130 120 JT = 1 130 M2 = M-2 DO 160 K=1,M2 MK = M-K BMK = 2.0*FLOAT(MK)*FM1/X - FM FM = FM1 FM1 = BMK IF(MK-N-1) 150,140,150 140 BESJ = BMK 150 JT = -JT S = 1+JT 160 ALPHA = ALPHA + BMK*S BMK = 2.0*FM1/X - FM IF(N) 180,170,180 170 BESJ = BMK 180 ALPHA = ALPHA + BMK BESJ = BESJ/ALPHA IF(ABS(BESJ-BPREV)-ABS(D*BESJ)) 200,200,190 190 BPREV = BESJ IER = 3 200 RETURN END Cut: Function BESJ / Function BESY FUNCTION BESY (N, X, IER) C PURPOSE C TO CALCULATE THE BESSEL FUNCTION OF SECOND KIND AND INTEGER ORDER C C USAGE C BESY ( N, X, IER ) C C N - ORDER OF THE BESSEL FUNCTION C X - ARGUMENT C IER - ERROR INDICATOR. C C SUBPROGRAMS REQUIRED C NONE C C INTRINSIC ARRAYS C NONE C C METHOD C THE FUNCTION IS A SHORT HAND VERSION OF THE SUBROUTINE C BESY DESCRIBED IN THE IBM-SYSTEM 360 SCIENTIFIC SUBROUTINE C PACKAGE (360A-CM-03X) VERSION II C C ............................................................. C C C COMPUTES THE BESSEL FUNCTION OF SECOND KIND FOR GIVEN C ORDER AND ARGUMENT C BESY = 0.0 C C CHECK FOR ERRORS IN N AND X C IF(N) 180,10,10 10 IER = 0 IF(X) 190,190,20 20 PI = 3.141592653 C C BRANCH IF X LESS THAN OR EQUAL TO 4 C IF(X-4.0) 40,40,30 C C COMPUTE Y0 AND Y1 FOR X GREATER THAN 4 C 30 T = 4.0/X P0 = 0.3989422793 Q0 = -0.0124669441 P1 = 0.3989422819 Q1 = 0.0374008364 A = T*T B = A P0 = P0 - 0.0017530620 * A Q0 = Q0 + 0.0004564324 * A P1 = P1 + 0.0029218256 * A Q1 = Q1 - 0.00063904 * A A = A*A P0 = P0 + 0.00017343 * A Q0 = Q0 - 0.0000869791 * A P1 = P1 - 0.000223203 * A Q1 = Q1 + 0.0001064741 * A A = A*B P0 = P0 - 0.0000487613 * A Q0 = Q0 + 0.0000342468 * A P1 = P1 + 0.0000580759 * A Q1 = Q1 - 0.0000398708 * A A = A*B P0 = P0 + 0.0000173565 * A Q0 = Q0 - 0.0000142078 * A P1 = P1 - 0.000020092 * A Q1 = Q1 + 0.00001622 * A A = A*B P0 = P0 - 0.0000037043 * A Q0 = Q0 + 0.0000032312 * A P1 = P1 + 0.0000042414 * A Q1 = Q1 - 0.0000036594 * A A = SQRT(2.0*PI) B = 4.0*A P0 = A*P0 Q0 = B*Q0/X P1 = A*P1 Q1 = B*Q1/X A = X-PI/4.0 B = SQRT(2.0/(PI*X)) Y0 = B*(P0*SIN(A) + Q0*COS(A)) Y1 = B*(-P1*COS(A)+Q1*SIN(A)) GO TO 90 C C COMPUTE Y0 AND Y1 FOR X LESS THAN OR EQUAL TO 4 C 40 XX = X/2.0 X2 = XX*XX T = ALOG(XX) + 0.5772156649 SUM = 0.0 TERM = T Y0 = T C DO 70 L=1,15 IF(L-1)50,60,50 50 SUM = SUM + 1.0/FLOAT(L-1) 60 FL = L TS = T - SUM TERM = (TERM*(-X2)/FL**2)*(1.0-1.0/(FL*TS)) 70 Y0 = Y0 + TERM TERM = XX*(T-0.5) SUM = 0.0 Y1 = TERM C DO 80 L=2,16 SUM = SUM + 1.0/FLOAT(L-1) FL = L FL1 = FL - 1.0 TS = T - SUM TERM = (TERM*(-X2)/(FL1*FL))*((TS-0.5/FL)/(TS+0.5/FL1)) 80 Y1 = Y1 + TERM PI2 = 2.0/PI Y0 = PI2*Y0 Y1 = -PI2/X+PI2*Y1 C C CHECK IF ONLY Y0 OR Y1 IS DESIRED C 90 IF(N-1)100,100,130 C C RETURN EITHER Y0 OR Y1 AS REQUIRED C 100 IF(N) 110,120,110 110 BESY = Y1 GO TO 170 120 BESY = Y0 GO TO 170 C C PERFORM RECURRENCE OPERATIONS TO FIND YN(X) C 130 YA = Y0 YB = Y1 K = 1 140 T = FLOAT(2*K)/X YC = T*YB - YA IF(ABS(YC)-1.0E35) 145,145,141 141 IER = 3 RETURN 145 K = K+1 IF(K-N)150,160,150 150 YA = YB YB = YC GO TO 140 160 BESY = YC 170 RETURN 180 IER = 1 RETURN 190 IER = 2 RETURN END Cut: Function BESY / Function BESI. FUNCTION BESI (N, X, IER) C. Function for the N'th order hyperbolic bessel function of argument X. C. N and X must be non-negative. C. For moderate order and argument, the Taylor series is applied. C. For large order and argument (12<N<X) an asymptotic series is applied. C. The routine is a modified version of the old IBM System/360 C. Scientific Subroutine Package. (360A-CM-03X) Version II. C. IER = 0 Normal exit. C. IER = 1 N is negative. C. IER = 2 X is negative. C. IER = 3 Overflow, X and N greater than 55. C. Prepared by $verre Gran April 1990. C. Standard limits and values. AN = FLOAT (N) EPSLN = 1.E-25 TOL = 1.E-6 XMAX = 55.0 PI = 3.141592653 C. Check for trivial cases and erroneous input. BESI = 1.0 IER = 0 IF (N. EQ. 0. AND. X. EQ. 0.0) GO TO 990 IF (N. LT. 0) THEN IER = 1 GO TO 990 END IF IF (X. LT. 0) THEN IER = 2 GO TO 990 END IF C. See if the asymptotic series is to be applied. IF (X. GT. 12.0. AND. X. GT. AN) GO TO 110 C. Calculate Bessel function by the series. First term and initial sum value. 40 XX = X / 2.0 50 TERM = 1.0 IF (N. LE. 0) GO TO 70 DO 60 I = 1, N FI = FLOAT (I) ATERM = ABS (TERM) IF (ATERM. LT. EPSLN) THEN BESI = 0.0 GO TO 990 END IF TERM = TERM * XX / FI 60 CONTINUE 70 BESI = TERM XX = XX * XX C. Calculate the series terms and stop when the terms are small enough. DO 90 K = 1, 1000 ATERM = ABS (TERM) BTOL = ABS (BESI * TOL) IF (ATERM. LE. BTOL) GO TO 990 FK = FLOAT (K) * FLOAT (N + K) TERM = TERM * XX / FK BESI = BESI + TERM 90 CONTINUE C. Use the asymptotic series. 110 FN = 4 * N * N IF (X. GT. XMAX) THEN IER = 3 GO TO 990 END IF XX = 1.0 / (8.0 * X) TERM = 1.0 BESI = 1.0 DO 130 K = 1, 30 ATERM = ABS (TERM) BTOL = ABS (BESI * TOL) IF (ATERM. LE. BTOL) GO TO 140 IFK = (2 * K - 1)**2 FK = FLOAT (IFK) TERM = TERM * XX * (FK - FN) / FLOAT (K) BESI = BESI + TERM 130 CONTINUE C. If no significant answer is obtained after 30 terms, try the series. GO TO 40 140 BESI = BESI * EXP(X) / SQRT (2.0 * PI * X) 990 RETURN END Cut: Function BESI. / Function ATAN4. FUNCTION ATAN4 (Y, X) C Calculates arcus tangent y/x by giveh y and x. C The routine applies the ATAN2 routine, but has some logics C for x and y equal to 0, where ATAN2 usually goes into error. C $verre, -89. PI = 3.141592654 IF (X. EQ. 0. AND. Y. EQ. 0.) THEN ATAN4 = 0.0 GO TO 99 ENDIF IF (X. EQ. 0. AND. Y. GT. 0.) THEN ATAN4 = 0.5 * PI GO TO 99 ENDIF IF (X. EQ. 0. AND. Y. LT. 0.) THEN ATAN4 = - 0.5 * PI GO TO 99 ENDIF IF (X. GT. 0. AND. Y. EQ. 0.) THEN ATAN4 = 0.0 GO TO 99 ENDIF IF (X. LT. 0. AND. Y. EQ. 0.) THEN ATAN4 = PI GO TO 99 ENDIF ATAN4 = ATAN2 (Y, X) 99 RETURN END Cut: Function ATAN4. / Bottom.