Esope source files

   1 : C SBFIT     SOURCE    CHAT      05/01/13    03:12:04     5004           
   2 :       SUBROUTINE SB(Xbar, Sigma, RTB1, B2, Gamma, Delta, Xlam, Xi,
   3 :      $ FAULT)
   4 :       IMPLICIT INTEGER(I-N)
   5 :       IMPLICIT REAL*8 (A-H,O-Z)
   6 :       REAL*8 HMU(6), Deriv(4), DD(4), Xbar, Sigma, RTB1, 
   7 :      $ Xlam, Xi, TT, Tol, RB1, B1, E, U, X, Y, W, F, D, G, S, H2, T,
   8 :      $ H2A, H2B, H3, H4, RBET, BET2,B2, Gamma, Delta
   9 :       LOGICAL NEG, FAULT
  10 :       DATA TT, Tol, LIMIT /1.0E-4, 0.01D0, 50/
  11 : C
  12 : C
  13 :       RB1 = DABS(RTB1)
  14 :       B1 = RB1 * RB1
  15 :       NEG = RTB1.LT.0.D0
  16 : C
  17 : C Calcul de D = premiere approximation de Delta
  18 : C
  19 :       E = B1 + 1.D0
  20 :       U = 1.D0 / 3.D0 
  21 :       X = 0.5D0 * B1 + 1.D0
  22 :       Y = RB1 * DSQRT(0.25D0 * B1 + 1.D0)
  23 :       W = (X + Y) ** U + (X - Y) ** U - 1.D0
  24 :       F = W * W * (3.D0 + W * (2.D0 + W)) - 3.D0
  25 :       E = (B2 - E) / (F-E)
  26 :       IF (DABS(RB1).GT.Tol) GOTO 5
  27 :       F = 2.D0
  28 :       GOTO 20
  29 :   5   D = 1.D0 / DSQRT(DLOG(W))
  30 :       IF (D.LT.0.64D0) WRITE (*,*) 'D < 0.64'
  31 :       IF (D.LT.0.64D0) GOTO 10
  32 :       F = 2.D0 - 8.5245D0 / (D * (D * (D - 2.163D0) + 11.346D0))
  33 :       GOTO 20
  34 :   10  F = 1.25D0 * D
  35 :   20  F = E * F + 1.D0
  36 :       IF (F.LT.1.8D0) GOTO 25
  37 :       D = (0.626D0 * F - 0.408D0) * (3.D0 - F) ** (-0.479D0)
  38 :       GOTO 30
  39 :   25  D = 0.8D0 * (F - 1.D0)
  40 : C
  41 : C Calcul de G = premiere approximation de Gamma
  42 : C
  43 :   30  G = 0.D0
  44 :       IF (B1.LT.TT) GOTO 70
  45 :       IF (D.GT.1.D0) GOTO 40
  46 :       G = (0.7466D0 * D ** 1.7973D0 + 0.5955D0)*B1**0.485D0
  47 :       GOTO 70
  48 :   40  IF (D.LE.2.5) GOTO 50
  49 :       U = 0.0124D0
  50 :       Y = 0.5291D0
  51 :       GOTO 60
  52 :   50  U = 0.0623D0
  53 :       Y = 0.4043D0
  54 :   60  G = B1 ** (U * D + Y) * (0.9281D0+D*(1.0614D0*D-0.7077D0))
  55 :   70  M = 0
  56 : C
  57 : C Boucle iterative principale
  58 : C 
  59 :   80  M = M + 1
  60 : C      WRITE (*,*) 'M =',M
  61 :       FAULT = M.GT.LIMIT
  62 :       IF (FAULT) WRITE (*,*) 'Depassement limite dans bcle principale'
  63 :       IF (FAULT) RETURN
  64 : C
  65 : C   Calcul des 6.D0 premiers moments associes a G et D
  66 : C
  66 :       CALL MOM(G, D, HMU, FAULT)
  67 : C      WRITE (*,*) '(MOM) G =',G
  68 : C      WRITE (*,*) '(MOM) D =',D
  69 : C      IF (FAULT) WRITE (*,*) 'ERREUR DANS MOM'
  70 :       IF (FAULT) RETURN
  71 :       S = HMU(1) * HMU(1)
  72 :       H2 = HMU(2) - S
  73 :       FAULT = H2.LE.0.D0
  74 :       IF (FAULT) RETURN
  75 :       T = DSQRT(H2)
  76 :       H2A = T * H2
  77 :       H2B = H2 * H2
  78 :       H3 = HMU(3) - HMU(1) * (3.D0 * HMU(2) - 2.D0 * S)
  79 :       RBET = H3 / H2A
  80 :       H4 = HMU(4) - HMU(1) * (4.D0 * HMU(3) - HMU(1) *
  81 :      $ (6.D0 * HMU(2)  - 3.D0 * S))
  82 :       BET2 = H4 / H2B
  83 :       W = G * D
  84 :       U = D * D
  85 : C
  86 : C Calcul des derivees
  87 : C
  88 :       DO 120 J = 1,2
  89 :         DO 110 K = 1,4
  90 :           T = K
  91 :           IF (J.EQ.1) GOTO 90
  92 :           S = ((W - T) * (HMU(K) - HMU(K+1)) + (T + 1.D0) *
  93 :      $    (HMU(K+1) - HMU(K+2))) / U
  94 :           GOTO 100
  95 :   90      S = HMU(K+1) - HMU(K)
  96 :  100      DD(K) = T * S / D
  97 :  110    CONTINUE
  98 :         T = 2.D0 * HMU(1) * DD(1)
  99 :         S = HMU(1) * DD(2)
 100 :         Y = DD(2) - T
 101 :         DERIV(J) = (DD(3) - 3.D0 * (S + HMU(2) * DD(1) - T * 
 102 :      $  HMU(1)) - 1.5D0 * H3 * Y / H2) / H2A
 103 :         DERIV (J+2) = (DD(4) - 4.D0 * (DD(3) * HMU(1) + DD(1) * 
 104 :      $  HMU(3)) + 6.D0 * (HMU(2) * T + HMU(1) * (S - T * HMU(1)))
 105 :      $  - 2.D0 * H4 * Y / H2) / H2B 
 106 :  120  CONTINUE
 107 :       T = 1.D0 / (DERIV(1) * DERIV(4) - DERIV(2) * DERIV(3))
 108 :       U = (DERIV(4) * (RBET - RB1) - DERIV(2) * (BET2 - B2)) * T
 109 :       Y = (DERIV(1) * (BET2 - B2) - DERIV(3) * (RBET - RB1)) * T
 110 : C
 111 : C Calcul des nouvelles estimations de G et D
 112 : C
 113 :       G = G - U
 114 :       IF (B1.EQ.0.D0.OR.G.LT.0.D0) G = 0.D0
 115 :       D = D - Y
 116 :       IF ( DABS(U).GT.TT.OR.DABS(Y).GT.TT) GOTO 80
 117 : C
 118 : C Fin des iterations
 119 : C
 120 : C      WRITE (*,*) 'Fin des iterations '
 121 :       Delta = D
 122 :       Xlam = Sigma / DSQRT(H2)
 123 :       IF (NEG) GOTO 130
 124 :       Gamma = G
 125 :       GOTO 140
 126 :  130  Gamma = -G
 127 :       HMU(1) = 1.D0 - HMU(1)
 128 :  140  Xi = Xbar - Xlam * HMU(1)
 129 :       RETURN
 130 :       END
 131 :  
 132 :  
 133 :  
 134 :