Sources Esope | Cast3M

jhnson
C JHNSON    SOURCE    CB215821  16/04/21    21:17:19     8920
      Subroutine JNSN(Xbar,Sd,RB1,BB2,Itype,Gamma,Delta,Xlam,Xi,
     * IFault)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 Xbar,Sd,RB1,BB2,Gamma,Delta,Xlam,Xi,Tol,
     * B1,B2,Y,X,U,W
      LOGICAL fault
      DATA Tol /0.01D0/
      IFault = 1
C
C Ecart-type plus petit que 0 --> IMPOSSIBLE
C
      IF (Sd.LT.0.D0) RETURN
C
C Initialisation
C
      IFault = 0
      Xi = 0.D0
      Xlam = 0.D0
      Gamma = 0.D0
      Delta = 0.D0
C
C Ecart-type = 0 --> Distribution de type 5 : type St
C
      IF (Sd.GT.0.D0) GOTO 10
      IType = 5
      Xi = Xbar
      RETURN
C
C Calcul de Beta_1 et Beta_2
C
  10  B1 = RB1 * RB1
      B2 = BB2
      fault = .FALSE.
C
C Voir si une distribution lognormale est demandee
C
      IF (B2.GE.0.D0) GOTO 30
  20  IF (ABS(RB1).LE.Tol) GOTO 70
      GOTO 80
C
C Position par rapport aux frontieres du domaine
C
C  domaine superieur
C
  30  IF (B2.GT.(B1 + Tol + 1.D0)) GOTO 60
C
C  domaine impossible
C
      IF (B2.LT.B1 + 1.D0) GOTO 50
C
C Distribution St
C
  40  IType = 5
      Y = 0.5D0 + 0.5D0 * SQRT(1.D0 - 4.D0/(B1 + 4.D0))
      IF (RB1.GT.0.D0)  Y = 1.D0 - Y
      X = Sd / SQRT(Y*(1.D0 - Y))
      Xi = Xbar - Y * X
      Xlam = Xi + X
      Delta = Y
      RETURN
  50  IFault = 2
      Itype = 6
      RETURN
  60  IF (ABS(RB1).GT.Tol.OR.ABS(B2 - 3.D0).GT.Tol) GOTO 80
C
C Distribution Normale
C
  70  IType = 4
      Delta = 1.D0 / Sd
      Gamma = - Xbar / Sd
      RETURN
C
C Test de position par rapport a la droite log-normale
C
  80  U = 1.D0 / 3.D0
      X = 0.5D0 * B1 + 1.D0
      Y = RB1 * SQRT(0.25D0 * B1 + 1.D0)
      W = (X + Y) ** U + (X - Y) ** U - 1.D0
      U = W * W * (3.D0 + W * (2.D0 + W)) - 3.D0
      IF (B2.LT.0.D0.OR.fault) B2 = U
      X = U - B2
      IF (ABS(X).GT.TOL) GOTO 90
C
C Distribution Log-Normale
C
      IType = 1
      Xlam =  SIGN(1.D0,rb1)
      U = Xlam * Xbar
      X = 1.D0 / SQRT(LOG(W))
      Delta = X
      Y = 0.5D0 * X * LOG(W * (W - 1.D0) / (Sd * Sd))
      Gamma = Y
      Xi = U - EXP((0.5D0 / X - Y) /X)
      RETURN
C
C Distribution Sb ou Su
C
  90  IF (X.GT.0.D0) GOTO 100
      ITYPE = 2
      CALL SU(Xbar,Sd,RB1,B2,Gamma,Delta,Xlam,Xi)
      RETURN
 100  ITYPE = 3
      CALL SB(Xbar,Sd,RB1,B2,Gamma,Delta,Xlam,Xi,FAULT)
      IF(.NOT.FAULT) RETURN
C
C Echec. Convergence non atteinte. Perturbation des param
C
      IFAULT = 3
      IF (B2.GT.B1 + 2.D0) GOTO 20
      GOTO 40
      END