92 lines
3.2 KiB
Fortran
92 lines
3.2 KiB
Fortran
FUNCTION SIGMAR(ALPHA,XMDT,TEF,OMEGA,RELR,RELT,RELZ)
|
|
C =====================================================
|
|
c
|
|
C--------------------------------------------------------------------
|
|
C The following function takes as inpute various disk parameters computed
|
|
C at a certain radius in cgs units, and outputs the surface mass density,
|
|
C assuming that the opacity is electron-scattering dominated (or that kappa
|
|
C is independent of the density). The equations were derived assuming a 1-zone
|
|
C model (i.e. rho, T_g, mu are constant with height), assuming t_r,phi =
|
|
C -alpha P, and that the dissipation is constant per unit optical depth.
|
|
C See chapter 7 of Krolik for information on the notation used here
|
|
C (especially relativistic factors).
|
|
C
|
|
C--------------------------------------------------------------------
|
|
C Uses Numerical Recipes subroutine LAGUER
|
|
C--------------------------------------------------------------------
|
|
C ALPHA - viscosity parameter
|
|
C XMDT - accretion rate (in g/s)
|
|
C TEF - temperature in Kelvins
|
|
C OMEGA - Keplerian frequency (in Hz)
|
|
C RELR,RELT, RELZ - relativistic factors
|
|
C KAPPA - opacity (cm^2/g)
|
|
C MU - mean atomic mass (g) = rho/N
|
|
C--------------------------------------------------------------------
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
REAL*8 KAPPA,MU
|
|
COMPLEX*16 COEFF(11),XGUESS
|
|
C
|
|
C We should check that the physical constants used here agree with those in
|
|
C disk195g:
|
|
C
|
|
parameter (ZERO=0.D0,ONE=1.D0,TRES=3.D0,
|
|
* FOUR=4.D0,THIRD=ONE/TRES,FOURTH=0.25D0,EPS=1.D-5)
|
|
parameter (C=2.9979D10,SIGMAB=5.6703D-5)
|
|
parameter (BK=1.3807D-16)
|
|
PI=ACOS(-ONE)
|
|
C
|
|
C We'll assume fully ionized, pure hydrogen:
|
|
C
|
|
KAPPA=0.39D0
|
|
MU=0.5D0*1.6726D-24
|
|
FAC1=RELZ*(HALF*TRES*C*OMEGA/ALPHA/KAPPA/SIGMAB/TEF**4)**2
|
|
FAC2=(HALF*KAPPA)**FOURTH*BK*TEF/MU
|
|
FAC3=XMDT*OMEGA*RELT/PI
|
|
C
|
|
C Coefficients of the equation for x^4=Sigma:
|
|
C
|
|
COEFF(1)=DCMPLX(FAC1*(HALF*FAC3)**2,ZERO)
|
|
COEFF(2)=ZERO
|
|
COEFF(3)=ZERO
|
|
COEFF(4)=ZERO
|
|
COEFF(5)=DCMPLX(-(TRES*FAC3)/(8.D0*ALPHA),ZERO)
|
|
COEFF(6)=DCMPLX(-FAC1*FAC3*ALPHA*FAC2,ZERO)
|
|
COEFF(7)=ZERO
|
|
COEFF(8)=ZERO
|
|
COEFF(9)=ZERO
|
|
COEFF(10)=DCMPLX(FOURTH*FAC2,ZERO)
|
|
COEFF(11)=DCMPLX(FAC1*(ALPHA*FAC2)**2,ZERO)
|
|
C
|
|
C At small radii, we'll approximate P_rad >> P_gas
|
|
C First, compute sigma assuming radiation pressure dominates:
|
|
C
|
|
SIGRAD=FOUR*OMEGA*C*C*RELT*RELZ/ALPHA/KAPPA**2/SIGMAB/TEF**4/RELR
|
|
C
|
|
C Next, compute sigma assuming gas pressure dominates:
|
|
C
|
|
SIGGAS=((MU*XMDT*OMEGA*RELT/PI/ALPHA/BK/TEF)**4/8./KAPPA)**0.2D0
|
|
C
|
|
C Use a starting guess which has the correct value for P_gas >> P_rad
|
|
C or P_rad >> P_gas.
|
|
C
|
|
XGUESS=DCMPLX(ONE/(ONE/SIGRAD+ONE/SIGGAS)**FOURTH,ZERO)
|
|
C
|
|
C Look for root of the 10th order equation for x:
|
|
C
|
|
CALL LAGUER(COEFF,10,XGUESS,ITS)
|
|
C
|
|
C Make sure that we haven't landed a wrong root:
|
|
C
|
|
IF(ABS(DIMAG(XGUESS)).LT.EPS.AND.DBLE(XGUESS).GT.ZERO) THEN
|
|
SIGMAR=DBLE(XGUESS)**4
|
|
ELSE
|
|
SIGMAR=ONE/(ONE/SIGRAD+ONE/SIGGAS)
|
|
WRITE(6,*) 'Surface density approximated'
|
|
ENDIF
|
|
WRITE(6,2000) TEF,SIGRAD,SIGGAS,SIGMAR
|
|
RETURN
|
|
2000 FORMAT(20(2x,1pe12.5))
|
|
END
|