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