Sources Esope | Cast3M

sb
C SB        SOURCE    CB215821  22/11/16    21:15:04     11500          
C SBFIT     SOURCE    CB215821  16/04/21    21:18:22     8920
      SUBROUTINE SB(Xbar, Sigma, RTB1, B2, Gamma, Delta, Xlam, Xi,
     $ FAULT)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 HMU(6), Deriv(4), DD(4), Xbar, Sigma, RTB1,
     $ Xlam, Xi, TT, Tol, RB1, B1, E, U, X, Y, W, F, D, G, S, H2, T,
     $ H2A, H2B, H3, H4, RBET, BET2,B2, Gamma, Delta
      LOGICAL NEG, FAULT
      DATA TT, Tol, LIMIT /1.0E-4, 0.01D0, 50/
C
C
      RB1 = ABS(RTB1)
      B1 = RB1 * RB1
      NEG = RTB1.LT.0.D0
C
C Calcul de D = premiere approximation de Delta
C
      E = B1 + 1.D0
      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
      F = W * W * (3.D0 + W * (2.D0 + W)) - 3.D0
      E = (B2 - E) / (F-E)
      IF (ABS(RB1).GT.Tol) GOTO 5
      F = 2.D0
      GOTO 20
  5   D = 1.D0 / SQRT(LOG(W))
      IF (D.LT.0.64D0) WRITE (*,*) 'D &lt; 0.64'
      IF (D.LT.0.64D0) GOTO 10
      F = 2.D0 - 8.5245D0 / (D * (D * (D - 2.163D0) + 11.346D0))
      GOTO 20
  10  F = 1.25D0 * D
  20  F = E * F + 1.D0
      IF (F.LT.1.8D0) GOTO 25
      D = (0.626D0 * F - 0.408D0) * (3.D0 - F) ** (-0.479D0)
      GOTO 30
  25  D = 0.8D0 * (F - 1.D0)
C
C Calcul de G = premiere approximation de Gamma
C
  30  G = 0.D0
      IF (B1.LT.TT) GOTO 70
      IF (D.GT.1.D0) GOTO 40
      G = (0.7466D0 * D ** 1.7973D0 + 0.5955D0)*B1**0.485D0
      GOTO 70
  40  IF (D.LE.2.5) GOTO 50
      U = 0.0124D0
      Y = 0.5291D0
      GOTO 60
  50  U = 0.0623D0
      Y = 0.4043D0
  60  G = B1 ** (U * D + Y) * (0.9281D0+D*(1.0614D0*D-0.7077D0))
  70  M = 0
C
C Boucle iterative principale
C
  80  M = M + 1
C      WRITE (*,*) 'M =',M
      FAULT = M.GT.LIMIT
      IF (FAULT) WRITE (*,*) 'Depassement limite dans bcle principale'
      IF (FAULT) RETURN
C
C   Calcul des 6.D0 premiers moments associes a G et D
C
      CALL MOM(G, D, HMU, FAULT)
C      WRITE (*,*) '(MOM) G =',G
C      WRITE (*,*) '(MOM) D =',D
C      IF (FAULT) WRITE (*,*) 'ERREUR DANS MOM'
      IF (FAULT) RETURN
      S = HMU(1) * HMU(1)
      H2 = HMU(2) - S
      FAULT = H2.LE.0.D0
      IF (FAULT) RETURN
      T = SQRT(H2)
      H2A = T * H2
      H2B = H2 * H2
      H3 = HMU(3) - HMU(1) * (3.D0 * HMU(2) - 2.D0 * S)
      RBET = H3 / H2A
      H4 = HMU(4) - HMU(1) * (4.D0 * HMU(3) - HMU(1) *
     $ (6.D0 * HMU(2)  - 3.D0 * S))
      BET2 = H4 / H2B
      W = G * D
      U = D * D
C
C Calcul des derivees
C
      DO 120 J = 1,2
        DO 110 K = 1,4
          T = K
          IF (J.EQ.1) GOTO 90
          S = ((W - T) * (HMU(K) - HMU(K+1)) + (T + 1.D0) *
     $    (HMU(K+1) - HMU(K+2))) / U
          GOTO 100
  90      S = HMU(K+1) - HMU(K)
 100      DD(K) = T * S / D
 110    CONTINUE
        T = 2.D0 * HMU(1) * DD(1)
        S = HMU(1) * DD(2)
        Y = DD(2) - T
        DERIV(J) = (DD(3) - 3.D0 * (S + HMU(2) * DD(1) - T *
     $  HMU(1)) - 1.5D0 * H3 * Y / H2) / H2A
        DERIV (J+2) = (DD(4) - 4.D0 * (DD(3) * HMU(1) + DD(1) *
     $  HMU(3)) + 6.D0 * (HMU(2) * T + HMU(1) * (S - T * HMU(1)))
     $  - 2.D0 * H4 * Y / H2) / H2B
 120  CONTINUE
      T = 1.D0 / (DERIV(1) * DERIV(4) - DERIV(2) * DERIV(3))
      U = (DERIV(4) * (RBET - RB1) - DERIV(2) * (BET2 - B2)) * T
      Y = (DERIV(1) * (BET2 - B2) - DERIV(3) * (RBET - RB1)) * T
C
C Calcul des nouvelles estimations de G et D
C
      G = G - U
      IF (B1.EQ.0.D0.OR.G.LT.0.D0) G = 0.D0
      D = D - Y
      IF ( ABS(U).GT.TT.OR.ABS(Y).GT.TT) GOTO 80
C
C Fin des iterations
C
C      WRITE (*,*) 'Fin des iterations '
      Delta = D
      Xlam = Sigma / SQRT(H2)
      IF (NEG) GOTO 130
      Gamma = G
      GOTO 140
 130  Gamma = -G
      HMU(1) = 1.D0 - HMU(1)
 140  Xi = Xbar - Xlam * HMU(1)
      RETURN
      END