34 lines
843 B
Fortran
34 lines
843 B
Fortran
FUNCTION BETAH(R)
|
|
C =================
|
|
C
|
|
C Determination of the total pressure scale height
|
|
C Solution of the transcendental equation by the Newton-Raphson method
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
PARAMETER (UN=1.D0,
|
|
* PISQ=1.77245385090551D0)
|
|
IF(R.LT.0.88) THEN
|
|
BET0=PISQ/2.D0/R
|
|
ELSE
|
|
BET0=UN+UN/3.D0/R/R
|
|
END IF
|
|
C
|
|
ITER=0
|
|
BETA=BET0
|
|
10 ITER=ITER+1
|
|
B1=BETA-UN
|
|
RB1=R*B1
|
|
BSQ=SQRT(BETA*B1)
|
|
ERF1=ERFCX(R*BSQ)
|
|
ERF2=ERFCX(RB1)
|
|
RHS=BSQ/B1*(UN-ERF1)+EXP(-R*RB1)*ERF2
|
|
DP=R/PISQ*(2.D0-EXP(-R*BETA*RB1))+(UN-ERF1)/2.D0/B1/BSQ+
|
|
* R*R*EXP(-R*RB1)*ERF2
|
|
DBETA=(RHS-2.D0/PISQ*BETA*R)/DP
|
|
DEL=DBETA/BETA
|
|
BETA=BETA+DBETA
|
|
IF(ABS(DEL).GT.1.D-5.AND.ITER.LE.10) GO TO 10
|
|
BETAH=BETA
|
|
RETURN
|
|
END
|