Sources Esope | Cast3M

su
C SU        SOURCE    CB215821  22/11/16    21:15:05     11500          
C SUFIT     SOURCE    CB215821  16/04/21    21:18:30     8920
      SUBROUTINE SU(Xbar,Sd,RB1,B2,Gma,Dlta,Xlmbd,Xi)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 Xbar, Sd, RB1, B2, Gma, Dlta,
     $ Xi, Tol, B1, B3, W, Y, W1, WM1, Z, V, A, B, X,
     $ Xlmbd
      DATA Tol /0.01D0/
C
C Calcul de Beta_1 et Beta_2 -3
C
      B1 = RB1 * RB1
      B3 = B2 - 3.D0
C
C W premiere estimation de EXP(Dlta**(-2))
C
      W = SQRT(SQRT(2.D0 * B2 - 2.8D0 * B1 - 2.D0) - 1.D0)
      IF (ABS(RB1).GT.Tol) GOTO 10
C
C Cas distribution symetrique - resultats theoriques
C
      Y = 0.D0
      GOTO 20
C
C Debut des iterations de Johnson
C
  10  W1 = W + 1.D0
      WM1 = W - 1.D0
      Z = W1 * B3
      V = W * (6.D0 + W * (3.D0 + W))
      A = 8.D0 * (WM1 *(3.D0 + W * (7.D0 + V)) - Z)
      B = 16.D0 * (WM1 * (6.D0 + V) - B3)
      Y = (SQRT(A * A - 2.D0 * B * (WM1 * (3.D0 + W *
     $ (9.D0 + W * (10.D0 + V))) - 2.D0 * W1 * Z)) - A) / B
      Z = Y * WM1 * (4.D0 * (W + 2.D0) * Y + 3.D0
     $ * W1 * W1)** 2 / (2.D0 * (2.D0 * Y + W1) ** 3)
      V = W * W
      W = SQRT(SQRT(1.D0 - 2.D0 * (1.5D0 - B2 + (B1 *
     $ (B2 - 1.5D0 - V * (1.D0 + 0.5D0 * V))) / Z)) - 1.D0)
      IF (ABS(B1 - Z).GT.Tol) GOTO 10
C
C Fin de la boucle iterative de Johnson
C
      Y = Y / W
      Y = LOG(SQRT(Y) + SQRT(Y + 1.D0))
      IF (RB1.GT.0.D0) Y = -Y
  20  X = SQRT(1.D0 / LOG(W))
      Gma = Y * X
      Dlta = X
      Y = EXP(Y)
      Z = Y * Y
      X = Sd / SQRT(0.5D0 * (W - 1.D0) * (0.5D0 * W *
     $ (Z + 1.D0 / Z) + 1.D0 ))
      Xlmbd = X
      Xi = (0.5D0 * SQRT(W) * (Y - 1.D0 / Y)) * X + Xbar
C      write(*,*) RB1,B2,Gma,Dlta,Xlmbd,Xi
      RETURN
      END