SpectraRust/tlusty/extracted/grcor.f
2026-03-19 14:05:33 +08:00

124 lines
3.4 KiB
Fortran

SUBROUTINE GRCOR(AA,RR,XMSTAR,QCOR,TCOR,ARH,BRH,CRH,DRH)
C =======================================================
C
C Procedure for computing general-relativistic correction
C factors to gravitational factor (QGRAV) and effective
C temperature (TEFF)
C Also calculates all frour quantities in the Riffer-Herlod (RH)
C notation - ARH, BRH, CRH, DRH
C
C Input:
C AA - angular momentum (0.98 maximum)
C RR - R/R_g = r/(GM/c^2)
C Outout:
C QCOR - g-correction = C/B in RH notation
C TCOR - T-correction = (D/B)^(1/4) in RH notation
C ARH - A in RH notation
C BRH - B in RH notation
C CRH - C in RH notation
C DRH - D in RH notation
C
INCLUDE 'IMPLIC.FOR'
PARAMETER (THIRD=1.D0/3.D0, PI3=1.0471976)
C
C ----------------
C Imput parameters
C ----------------
C
C AA - specific angular momentum/mass
C of the Kerr black hole
C RR - distance/mass of the Kerr black hole
C
C -----------------------------------
C Classical case - no GR corrections
C ------------------------------------
C
if(Xmstar.gt.0.) then
arh=1.
brh=1.
crh=1.
drh=1.-sqrt(1./rr)
qcor=1.
tcor=drh**0.25
return
end if
c
C ---------------------------------
C Set correcion factors A through G (see Novikov & Thorne,'73, eq.5.4.1a-g)
C ---------------------------------
C
rror=rr
rr=abs(rr)
AA2=AA*AA
RR1=1/RR
RR12=SQRT(RR1)
RR2=RR1*RR1
A2R2=AA2*RR2
A4R4=A2R2*A2R2
A2R3=AA2*RR2*RR1
AR32=SQRT(A2R3)
C
A = 1 + A2R2 + 2*A2R3
B = 1 + AR32
C = 1 - 3*RR1 + 2*AR32
D = 1 - 2*RR1 + A2R2
E = 1 + 4*A2R2 - 4*A2R3 + 3*A4R4
C
C -------------------------------
C Set correction factor for QGRAV (see Novikov & Thorne,'73, eq.5.7.2)
C -------------------------------
C
if(rror.lt.0) QCOR = B*B*D*E/(A*A*C)
c
c correction - after Riffert and Harold
c
if(rror.gt.0) QCOR = (1. - 4.*AR32 + 3.*A2R2)/C
C
C -----------------------
C Set correction factor Q (see Page & Thorne,'73, eq.35)
C -----------------------
C
C Minimum radius for last stable circular orbit per unit mass, X0
C
Z1 = 1 + (1-AA2)**THIRD * ((1+AA)**THIRD + (1-AA)**THIRD)
Z2 = SQRT(3*AA2 + Z1*Z1)
X0 = SQRT(3 + Z2 - SQRT((3-Z1)*(3+Z1+2*Z2)))
C
C Roots of x^3 - 3x + 2a = 0
C
CA3 = THIRD * ACOS(AA)
X1 = 2*COS(CA3-PI3)
X2 = 2*COS(CA3+PI3)
X3 = -2*COS(CA3)
C
C FB = '[]' term in eq. (35) of Page&Thorne '73
C
X = SQRT(RR)
C1 = 3*(X1-AA)*(X1-AA)/(X1*(X1-X2)*(X1-X3))
C2 = 3*(X2-AA)*(X2-AA)/(X2*(X2-X1)*(X2-X3))
C3 = 3*(X3-AA)*(X3-AA)/(X3*(X3-X1)*(X3-X2))
AL0 = 1.5*AA*log(X/X0)
AL1 = log((X-X1)/(X0-X1))
AL2 = log((X-X2)/(X0-X2))
AL3 = log((X-X3)/(X0-X3))
FB = (X-X0 - AL0 - C1*AL1 - C2*AL2 - C3*AL3)
Q = FB*(1+AR32)*RR12/SQRT(1-3*RR1+2*AR32)
C ------------------------------
C Set correction factor for TEFF (see Novikov & Thorne,'73, eq.5.5.14b)
C ------------------------------
C
TCOR = (Q/B/SQRT(C))**0.25
C
C ------------------------------
C RH quantities
C ------------------------------
C
ARH = D
BRH = C
CRH = 1. - 4.*AR32 + 3.*A2R2
DRH = Q/B*SQRT(C)
C
RETURN
END