63 lines
1.7 KiB
Fortran
63 lines
1.7 KiB
Fortran
SUBROUTINE CUBIC(DELTA)
|
|
C =======================
|
|
C
|
|
C Solution of the cubic equation for determination of
|
|
C the true gradient DELTA
|
|
C
|
|
C Input: A,B - coefficients; transmitted by COMMON/CUBCON
|
|
C DEL - DELTA(RAD) - DELTA(ADIAB); also transm. by CUBCON
|
|
C Output: DELTA - true gradient
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
COMMON/CUBCON/A,B,DEL,GRDADB,DELMDE,RHO,FLXTOT,GRAVD
|
|
PARAMETER (THIRD = 0.333333333333333D0)
|
|
data ipri /0/
|
|
C
|
|
C first solve the cubic equation
|
|
C A*X**3 + X**2 + B*X = DEL
|
|
C where X = (DELTA - DELTA(ELEM))**(1/2)
|
|
C
|
|
AA=THIRD/A
|
|
BB=B/A
|
|
CC=-DEL/A
|
|
P=BB*THIRD-AA*AA
|
|
Q=AA**3-(BB*AA-CC)/2.D0
|
|
D=Q*Q+P*P*P
|
|
IF(D.GT.0.) THEN
|
|
D=SQRT(D)
|
|
if(d-abs(q).lt.1.e-14*d) then
|
|
SOL=(2.D0*D)**THIRD-AA
|
|
ELSE
|
|
D1=ABS(D-Q)
|
|
D2=ABS(D+Q)
|
|
SOL=D1/(D-Q)*D1**THIRD-D2/(D+Q)*D2**THIRD-AA
|
|
END IF
|
|
ELSE
|
|
COSF=-Q/SQRT(ABS(P*P*P))
|
|
TANF=SQRT(UN-COSF*COSF)/COSF
|
|
FI=ATAN(TANF)*THIRD
|
|
SOL=2.D0*SQRT(ABS(P))*COS(FI)-AA
|
|
END IF
|
|
C
|
|
C if the previous formalism gives an unphysical solution
|
|
C x > DEL, then find the physical solution in the range (0, DEL)
|
|
C by a Newton-Raphson solution of the cubic equation
|
|
C
|
|
DELDA=SOL*(B+SOL)
|
|
IF(DELDA.GT.DEL.OR.DELDA.LT.0.) THEN
|
|
X0=sol
|
|
J=0
|
|
10 DELX=(DEL-X0*(B+X0+A*X0*X0))/(3.D0*A*X0*X0+2.D0*X0+B)
|
|
X0=X0+DELX
|
|
J=J+1
|
|
IF(ABS(DELX/X0).GT.1D-6.AND.J.LT.50) GO TO 10
|
|
SOL=X0
|
|
END IF
|
|
C
|
|
C finally, the actual gradient delta
|
|
C
|
|
DELTA=GRDADB+B*SOL+SOL*SOL
|
|
RETURN
|
|
END
|