Esope source files

   1 : C JHNSON    SOURCE    CHAT      05/01/13    00:49:36     5004           
   2 :       Subroutine JNSN(Xbar,Sd,RB1,BB2,Itype,Gamma,Delta,Xlam,Xi,
   3 :      * IFault)
   4 :       IMPLICIT INTEGER(I-N)
   5 :       IMPLICIT REAL*8 (A-H,O-Z)
   6 :       REAL*8 Xbar,Sd,RB1,BB2,Gamma,Delta,Xlam,Xi,Tol,
   7 :      * B1,B2,Y,X,U,W
   8 :       LOGICAL fault
   9 :       DATA Tol /0.01D0/
  10 :       IFault = 1
  11 : C
  12 : C Ecart-type plus petit que 0 --> IMPOSSIBLE
  13 : C
  14 :       IF (Sd.LT.0.D0) RETURN
  15 : C
  16 : C Initialisation
  17 : C
  18 :       IFault = 0
  19 :       Xi = 0.D0
  20 :       Xlam = 0.D0
  21 :       Gamma = 0.D0
  22 :       Delta = 0.D0
  23 : C
  24 : C Ecart-type = 0 --> Distribution de type 5 : type St
  25 : C
  26 :       IF (Sd.GT.0.D0) GOTO 10
  27 :       IType = 5
  28 :       Xi = Xbar
  29 :       RETURN
  30 : C
  31 : C Calcul de Beta_1 et Beta_2
  32 : C
  33 :   10  B1 = RB1 * RB1
  34 :       B2 = BB2
  35 :       fault = .FALSE.
  36 : C
  37 : C Voir si une distribution lognormale est demandee
  38 : C
  39 :       IF (B2.GE.0.D0) GOTO 30
  40 :   20  IF (DABS(RB1).LE.Tol) GOTO 70
  41 :       GOTO 80
  42 : C
  43 : C Position par rapport aux frontieres du domaine
  44 : C
  45 : C  domaine superieur
  46 : C
  47 :   30  IF (B2.GT.(B1 + Tol + 1.D0)) GOTO 60
  48 : C
  49 : C  domaine impossible
  50 : C
  51 :       IF (B2.LT.B1 + 1.D0) GOTO 50
  52 : C
  53 : C Distribution St
  54 : C
  55 :   40  IType = 5
  56 :       Y = 0.5D0 + 0.5D0 * DSQRT(1.D0 - 4.D0/(B1 + 4.D0))
  57 :       IF (RB1.GT.0.D0)  Y = 1.D0 - Y
  58 :       X = Sd / DSQRT(Y*(1.D0 - Y))
  59 :       Xi = Xbar - Y * X
  60 :       Xlam = Xi + X
  61 :       Delta = Y
  62 :       RETURN
  63 :   50  IFault = 2
  64 :       Itype = 6
  65 :       RETURN
  66 :   60  IF (DABS(RB1).GT.Tol.OR.DABS(B2 - 3.D0).GT.Tol) GOTO 80
  67 : C
  68 : C Distribution Normale
  69 : C
  70 :   70  IType = 4
  71 :       Delta = 1.D0 / Sd
  72 :       Gamma = - Xbar / Sd
  73 :       RETURN
  74 : C
  75 : C Test de position par rapport a la droite log-normale
  76 : C
  77 :   80  U = 1.D0 / 3.D0
  78 :       X = 0.5D0 * B1 + 1.D0
  79 :       Y = RB1 * DSQRT(0.25D0 * B1 + 1.D0)
  80 :       W = (X + Y) ** U + (X - Y) ** U - 1.D0
  81 :       U = W * W * (3.D0 + W * (2.D0 + W)) - 3.D0
  82 :       IF (B2.LT.0.D0.OR.fault) B2 = U
  83 :       X = U - B2
  84 :       IF (DABS(X).GT.TOL) GOTO 90
  85 : C
  86 : C Distribution Log-Normale
  87 : C 
  88 :       IType = 1
  89 :       Xlam =  SIGN(1.D0,rb1)
  90 :       U = Xlam * Xbar
  91 :       X = 1.D0 / DSQRT(DLOG(W))
  92 :       Delta = X
  93 :       Y = 0.5D0 * X * DLOG(W * (W - 1.D0) / (Sd * Sd))
  94 :       Gamma = Y
  95 :       Xi = U - DEXP((0.5D0 / X - Y) /X)
  96 :       RETURN
  97 : C
  98 : C Distribution Sb ou Su
  99 : C
 100 :   90  IF (X.GT.0.D0) GOTO 100
 101 :       ITYPE = 2
 101 :       CALL SU(Xbar,Sd,RB1,B2,Gamma,Delta,Xlam,Xi)
 102 :       RETURN
 103 :  100  ITYPE = 3
 103 :       CALL SB(Xbar,Sd,RB1,B2,Gamma,Delta,Xlam,Xi,FAULT)
 104 :       IF(.NOT.FAULT) RETURN
 105 : C
 106 : C Echec. Convergence non atteinte. Perturbation des param
 107 : C   
 108 :       IFAULT = 3
 109 :       IF (B2.GT.B1 + 2.D0) GOTO 20
 110 :       GOTO 40
 111 :       END
 112 :  
 113 :  
 114 :  
 115 :  
 116 :