stable.f

CSTAB
C     RANDOM STABLE STANDARDIZED FORM
      FUNCTION RSTAB(ALPHA,BPRIME,U,W)
C
C     ARGUMENTS..
C           ALPHA = CHARACTERISTIC EXPONENT
C          BPRIME = SKEWNESS IN REVISED PARAMETERIZATION
C              U = UNIFORM VARIATE ON (0,1). EXAMPLE FROM A UNIFORM
C              W = PSEUDO-RANDOM NUMBER GENERATOR
C               = EXPONENTIALLY DISTRIBUTED VARIATE
C
      DOUBLE PRECISION DA,DB
      DATA PIBY2/1.57079633/,PIBY4/.785398163/
      DATA THR1/0.99/
      EPS=1.0-ALPHA
C COMPUTE SOME TANGENTS
      PIBY3=PIBY2+PIBY2*(0.0-0.5)
      ATAN(EPS)=TAN2(PIBY3)
      B=TAN2(EPS*PIBY2)
      B=EPS*BPRIME*PIBY2
      IF(EPS.GT.0.99)TAU=BPRIME/(TAN2(EPS*PIBY2)*PIBY2)
      IF(EPS.LE.0.99)TAU=BPRIME*PIBY2+EPS*(1.-EPS)*TAN2((1.-EPS)*PIBY2)
C COMPUTE SOME NECESSARY SUBEXPRESSIONS
C IF PI NEAR PI BY 2, USE DOUBLE PRECISION.
      IF(A.GT.THR1)GO TO 50
C SINGLE PRECISION
      A2=A**2
      A2P=1.+A2
      A2=1.-A2
      B2=B**2
      B2P=1.+B2
      B2=1.-B2
      GO TO 100
C DOUBLE PRECISION
   50 DA=DBLE(A)**2
      DB=DBLE(B)**2
      A2=1.D0+DA
      A2P=1.D0+DA
      B2=1.D0-DB
      B2P=1.D0+DB
C COMPUTE COEFFICIENT
  100 Z=A2P**(2.*PIBY2*SB*TAU)/(W*A2*B2P)
C COMPUTE THE EXPONENTIAL-TYPE EXPRESSION
      ALOGS=ALOG(Z)
      D=D2(EPS*ALOGS/(1.-EPS))+(ALOGS/(1.-EPS))
C COMPUTE STABLE
      RSTAB=(1.+EPS*D)* 2.*((A-B)*(1.+A*B) - PIBY2*TAU*BB*(B*A2-2.*A))
     $/ (A2*B2P)
     +TAU*D
      RETURN
      END

CD2
      FUNCTION D2(Z)
C     EVALUATE (EXP(X) -1)/X
C     DOUBLE PRECISION P1,P2,P1,Q2,Q3,PV,ZZ
      DATA P1,P2,Q1,Q2,Q3/.50000 68525 36483 239 D3
     &,.20001 11415 89964 569 D2
     &,.16601 33705 07926 648 D1
     &,.14201 33704 07390 023 D3
     &, 1.DO/
C THE APPROXIMATION 1801 FROM HART ET AL (1968, P. 213)
      IF(ABS(Z).GT.0.1)GO TO 100
      ZZ=Z*Z
      PV=1+Z*P2
      D2=P0+ZV*(Q1+ZZ*(Q2+ZZ*Q3)-Z*PV)
      RETURN
  100 D2=(EXP(Z)-1.0)/Z
      RETURN
      END

CTAN
      FUNCTION TAN(XARG)
      LOGICAL NEG,INV
      DATA P0,P1,P2,Q0,Q1,Q2
     &/.129551035E+3,-.387662377E+1,.528644856E-1,
     &.164529332E+3,-.851320561E+2,1.0/
C THE APPROXIMATION 4283 FROM HART ET AL(1968, P. 251)
      DATA PIBY2/.785398163/,PIBY2/1.57079633/
      DATA PI/3.14159265/
      NEG=.FALSE.
      INV=.FALSE.
      X=XARG
      NEGAT=X.LT.0.
      X=ABS(X)
C PERFORM RANGE REDUCTION IF NECESSARY
      IF(X.LE.PIBY4)GO TO 50
      X=AMOD(X,PI)
      IF(X.LE.PIBY2)GO TO 30
      NEG=OR.TRUE.
      X=PI-X
   30 IF(X.LE.PIBY4)GO TO 50
      INV=.TRUE.
      X=PIBY2-X
   50 X=X/PIBY4
C CONVERT TO RANGE OF RATIONAL
      XX=X*X
      TAN=X*(P0+XX*(P1+XX*P2))/(Q0+XX*(Q1+XX*Q2))
      IF(NEG)TAN=-TAN
      IF(INV)TAN=1./TAN
      RETURN
      END

CTAN2
C     COMPUTE TAN(X)/X
C     FUNCTION DEFINED ONLY FOR ABS(XARG).LE.PI BY &
C     FOR OTHER ARGUMENTS RETURNS TAN(X)/X. COMPUTED DIRECTLY
      FUNCTION TAN2(XARG)
      DATA P0,P1,P2,Q0,Q1,Q2
     &/.129551035E+3,-.387662377E+1,.528644856E-1,
     &.164529332E+3,-.851320561E+2,1.0/
C THE APPROXIMATION 4283 FROM HART ET AL(1968, P. 251)
      DATA PIBY2/.785398163/,PIBY2/1.57079633/
      DATA PI/3.14159265/
      TAN2(XARG)
      IF(X.GT.PIBY4)GO TO 200
      X=X/PIBY4
C CONVERT TO RANGE OF RATIONAL
      XX=X*X
      TAN=X*(P0+XX*(P1+XX*P2))/(PIBY4*(Q0+XX*(Q1+XX*Q2)))
      RETURN
  200 TAN=TAN(XARG)/XARG
      RETURN
      END