PRO MLHOM
;;; V 1.0 2011/07/04 ;;;
;;; This routine implements the M-L relations from ;;;
;;; Graefener, Vink, de Koter, & Langer (2011)     ;;;

;;; Pre-definition of variables ;;;

logL = 5.0
M    = 10.0
XH   = 0.7
LMH  = 'M'
CH   = 1
NEXT = 'Q'

;;; Select which relation to use ;;;

READ,' What do you want: 1) M->L, 2) L->M, 3) M->Gamma, 4) Gamma->M ? ', CH

IF CH EQ 1 OR CH EQ 2 THEN BEGIN

MENU1: READ,' Mass Range I) Info, L) Low, M) Mid, H) High: ', LMH

IF LMH EQ 'i' OR LMH EQ 'I' THEN BEGIN
    PRINT, ' L) Low     HOM:   2-100 HEB: 0.6-100'
    PRINT, ' M) Mid     HOM:  12-250 HEB:   8-250'
    PRINT, ' H) High    HOM: 60-4000 HEB: 60-1000'
    GOTO, MENU1
ENDIF

IF LMH EQ 'l' OR LMH EQ 'L' THEN BEGIN
    ILMhom=2
    ILMheb=7
    IMLhom=12
    IMLheb=17
    PRINT, ' Mass Range HOM:   2-100 HEB: 0.6-100'
ENDIF

IF LMH EQ 'm' OR LMH EQ 'M' THEN BEGIN
    ILMhom=1
    ILMheb=6
    IMLhom=11
    IMLheb=16
    PRINT, ' Mass Range HOM:  12-250 HEB:   8-250'
ENDIF

IF LMH EQ 'h' OR LMH EQ 'H' THEN BEGIN
    ILMhom=3
    ILMheb=8
    IMLhom=13
    IMLheb=18
    PRINT, ' Mass Range HOM: 60-4000 HEB: 60-1000'
ENDIF

ENDIF

IF CH EQ 3 OR CH EQ 4 THEN BEGIN

MENU2: READ,' Mass Range I) Info, L) Low, H) High: ', LMH

IF LMH EQ 'i' OR LMH EQ 'I' THEN BEGIN
    PRINT, ' L) Low     HOM:    2-30 HEB: 0.3-100'
    PRINT, ' H) High    HOM: 12-4000 HEB:  12-500'
    GOTO, MENU2
ENDIF

IF LMH EQ 'l' OR LMH EQ 'L' THEN BEGIN
    IGMhom=4
    IGMheb=9
    IMGhom=14
    IMGheb=19
    PRINT, ' Mass Range HOM:    2-30 HEB: 0.3-100'
ENDIF

IF LMH EQ 'h' OR LMH EQ 'H' THEN BEGIN
    IGMhom=5
    IGMheb=10
    IMGhom=15
    IMGheb=20
    PRINT, ' Mass Range HOM: 12-4000 HEB:  12-500'
ENDIF

ENDIF

;;; Coefficients from Table (A.1). Zeros are added to let indices start at 1.
FC = [[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. ],$
      [ 0., 2.875, -3.966, 2.496, 2.652, -.31, -.511, 0., 0., 0. ],$
      [ 0., 1.967, -2.943, 3.755, 1.206, -.727, -.026, 0., 0., 0. ],$
      [ 0., 3.862, -2.486, 1.527, 1.247, -.076, -.183, 0., 0., 0. ],$
      [ 0., -2.688, -7.843, 2.471, 2.758, -.233, -.747, 0., 0., 0. ],$
      [ 0., -2.416, -5.118, 1.869, -.4, .064, .05, 0., 0., 0. ],$
      [ 0., 3.017, 2.446, -.306, 0., 0., 0., 0., 0., 0. ],$
      [ 0., 2.635, 2.986, -.488, 0., 0., 0., 0., 0., 0. ],$
      [ 0., 3.826, 1.619, -.099, 0., 0., 0., 0., 0., 0. ],$
      [ 0., -2.204, 1.831, .149, 0., 0., 0., 0., 0., 0. ],$
      [ 0., -1.676, 1.075, .404, 0., 0., 0., 0., 0., 0. ],$
      [ 0., 4.026, 4.277, -1., 25.48, 36.93, -2.792, -3.226, -5.317, 1.648 ],$
      [ 0., 2.582, .829, -1., 9.375, .333, .543, -1.376, -.049, .036 ],$
      [ 0., 10.05, 8.204, -1., 151.7, 254.5, -11.46, -13.16, -31.68, 2.408 ],$
      [ 0., 5.303, 5.918, -1., 16.58, -4.292, -72.89, -7.881, -13.76, 3.206 ],$
      [ 0., -14.6, 3.125, 1., 251., 15.63, 72.24, 18.2, 12.21, .781 ],$
      [ 0., 3.997, -1., 25.83, -3.268, 0., 0., 0., 0., 0. ],$
      [ 0., 3.059, -1., 14.76, -2.049, 0., 0., 0., 0., 0. ],$
      [ 0., 8.177, -1., 105.5, -10.1, 0., 0., 0., 0., 0. ],$
      [ 0., -6.144, 1., 52.54, 6.711, 0., 0., 0., 0., 0. ],$
      [ 0., -1.33, 1., 5.919, 2.475, 0., 0., 0., 0., 0. ]]
;;; FC -> FC(row,column)
FC = TRANSPOSE(FC)


;;; Apply M-L relations
IF CH EQ 1 THEN BEGIN
 READ,'     M/Msun: ', M
 READ,'XH(surface): ', XH

 REPEAT1: logM  = ALOG10(M)

;;; Eq. (9)
 LMhom = FC(ILMhom,1)+XH*FC(ILMhom,2)+(FC(ILMhom,3)+XH*FC(ILMhom,4))*logM $
       +(FC(ILMhom,5)+XH*FC(ILMhom,6))*logM*logM
;;; Eq. (10)
 LMheb = FC(ILMheb,1)+logM*FC(ILMheb,2)+logM*logM*FC(ILMheb,3)

;;; Compute Gamma_e using Eq. (8)
 G_hom  = 10.^(-4.813 + ALOG10(1.+XH) - logM + LMhom)
 G_heb  = 10.^(-4.813 + ALOG10(1.+XH) - logM + LMheb)
 G_heb0 = 10.^(-4.813 - logM + LMheb)
 

 PRINT, ' '
 PRINT, ' Input: log(M/Msun) =',  logM, ',  XH(surface)=', XH, FORMAT='(A,F8.3,A,F8.3)'
 PRINT, 'Output: HOM: log(L) =', LMhom, ', HEB: log(L) =', LMheb, FORMAT='(A,F8.3,A,F8.3)'
 PRINT, '    --> HOM: Gamma_e=', G_hom, ', HEB: Gamma_e=', G_heb, ', (XH=0.0:', G_heb0,')', FORMAT='(A,F8.3,A,F8.3,A,F8.3,A)' 
 PRINT, ' '

 READ,  'New M (m), XH (x), or Quit (q)? ', NEXT
 IF NEXT EQ 'M' OR NEXT EQ 'm' THEN BEGIN
    READ,'     M/Msun: ', M
    GOTO, REPEAT1
 ENDIF
 IF NEXT EQ 'X' OR NEXT EQ 'x' THEN BEGIN
    READ,'XH(surface): ', XH
    GOTO, REPEAT1
 ENDIF
 
ENDIF


;;; Apply L-M relations
IF CH EQ 2 THEN BEGIN
 READ,'log(L/Lsun): ', logL
 READ,'XH(surface): ', XH

 REPEAT2:

;;; Eqs. (11)+(12)
 MLhom = (FC(IMLhom,1) + XH*FC(IMLhom,2) $
        + FC(IMLhom,3)*SQRT(FC(IMLhom,4)+XH*FC(IMLhom,5)+XH*XH*FC(IMLhom,6) $
                          +(FC(IMLhom,7)+XH*FC(IMLhom,8))*logL)) $
       / (1+FC(IMLhom,9)*XH)
;;; Eq. (13)
 MLheb = FC(IMLheb,1) + FC(IMLheb,2) * SQRT(FC(IMLheb,3)+FC(IMLheb,4)*logL)

 Mhom = 10.^MLhom
 Mheb = 10.^MLheb

;;; Compute Gamma using Eq. (8)
 G_hom  = 10.^(-4.813 + ALOG10(1.+XH) - MLhom + logL)
 G_heb  = 10.^(-4.813 + ALOG10(1.+XH) - MLheb + logL)
 G_heb0 = 10.^(-4.813 - MLheb + logL)

 PRINT, ' '
 PRINT, ' Input: log(L/Lsun) =',  logL, ',  XH(surface)=', XH, FORMAT='(A,F8.3,A,F8.3)'
 PRINT, 'Output: HOM: M/Msun =',  Mhom, ', HEB: M/Msun =', Mheb, FORMAT='(A,F8.3,A,F8.3)'
 PRINT, '    --> HOM: Gamma_e=', G_hom, ', HEB: Gamma_e=', G_heb, ', (XH=0.0:', G_heb0,')', FORMAT='(A,F8.3,A,F8.3,A,F8.3,A)' 
 PRINT, ' '

 READ,  'New L (l), XH (x), or Quit (q)? ', NEXT
 IF NEXT EQ 'L' OR NEXT EQ 'l' THEN BEGIN
    READ,'log(L/Lsun): ', logL
    GOTO, REPEAT2
 ENDIF
 IF NEXT EQ 'X' OR NEXT EQ 'x' THEN BEGIN
    READ,'XH(surface): ', XH
    GOTO, REPEAT2
 ENDIF

ENDIF


;;; Apply M-Gamma relations
IF CH EQ 3 THEN BEGIN
 READ,'     M/Msun: ', M
 READ,'XH(surface): ', XH

 REPEAT3: logM  = ALOG10(M)

;;; Eq. (B.7), result is log(G_4)
 GMhom =  FC(IGMhom,1) + ALOG10(1.+XH)*FC(IGMhom,2) $ 
       + (FC(IGMhom,3) + ALOG10(1.+XH)*FC(IGMhom,4))*logM $
       + (FC(IGMhom,5) + ALOG10(1.+XH)*FC(IGMhom,6))*logM*logM
;;; Eq. (B.8), result is log(G_4)
 GMheb =  FC(IGMheb,1) + logM*FC(IGMheb,2) + logM*logM*FC(IGMheb,3)

;;; Compute Gamma from G_4 using Eqs. (B.9)-(B.12)
 G4 = 10.^GMhom
  C = 1./SQRT(12.)
  D = (C/3.*SQRT(256.*G4+27.) + 0.5) / G4^2.
  E = 3.*G4*D^(2./3.)-4.
 G_hom = 1. - C/D^(1./6.) * SQRT(E/G4) $ 
              * ( SQRT(3.*SQRT(D*G4)/(E^(3./2.)*C) -1.)-1.)
 G4 = 10.^GMheb
  C = 1./SQRT(12.)
  D = (C/3.*SQRT(256.*G4+27.) + 0.5) / G4^2.
  E = 3.*G4*D^(2./3.)-4.
 G_heb0 = 1. - C/D^(1./6.) * SQRT(E/G4) $ 
               * ( SQRT(3.*SQRT(D*G4)/(E^(3./2.)*C) -1.)-1.)
;;; NOTE: G_heb0 is Gamma_e for pure He star.
;;;       We transform this to Gamma_e for a given XH(surface)
 G_heb = G_heb0 * (1. + XH)

;;; Compute log(L/Lsun) using Eq. (8)
 Lhom = 4.813 - ALOG10(1.+XH) + ALOG10(G_hom) + logM
 Lheb = 4.813 - ALOG10(1.+XH) + ALOG10(G_heb) + logM

 PRINT, ' '
 PRINT, ' Input: log(M/Msun) =',  logM, ',  XH(surface)=', XH, FORMAT='(A,F8.3,A,F8.3)'
 PRINT, 'Output: HOM: Gamma_e=', G_hom, ', HEB: Gamma_e=', G_heb, ', (XH=0.0:', G_heb0,')', FORMAT='(A,F8.3,A,F8.3,A,F8.3,A)' 
 PRINT, '    --> HOM: log(L) =',  Lhom, ', HEB: log(L) =', Lheb, FORMAT='(A,F8.3,A,F8.3)'
 PRINT, ' '

 READ,  'New M (m), XH (x), or Quit (q)? ', NEXT
 IF NEXT EQ 'M' OR NEXT EQ 'm' THEN BEGIN
    READ,'     M/Msun: ', M
    GOTO, REPEAT3
 ENDIF
 IF NEXT EQ 'X' OR NEXT EQ 'x' THEN BEGIN
    READ,'XH(surface): ', XH
    GOTO, REPEAT3
 ENDIF

ENDIF


;;; Apply Gamma-M relations
IF CH EQ 4 THEN BEGIN
 READ,'    Gamma_e: ', Gamma
 READ,'XH(surface): ', XH

 REPEAT4: logG4 = ALOG10(Gamma/(1.-Gamma)^4.)

;;; Eqs. (B.13) & (B.14)
 MGhom = (FC(IMGhom,1) + ALOG10(1.+XH)*FC(IMGhom,2) $
       +  FC(IMGhom,3)*SQRT(FC(IMGhom,4) + logG4*FC(IMGhom,5) $ 
                          + ALOG10(1.+XH)^2.*FC(IMGhom,6) $
                          + (FC(IMGhom,7) + logG4*FC(IMGhom,8))*ALOG10(1.+XH))) $
       / (1. + FC(IMGhom,9)*ALOG10(1.+XH))
 Mhom = 10.^MGhom

;;; Eq. (B.15)
;;; First, apply the relation directly, i.e., for XH(surface)=0.0
 MGheb0 = FC(IMGheb,1) + FC(IMGheb,2) * SQRT(FC(IMGheb,3) + FC(IMGheb,4)*logG4)
 Mheb0  = 10.^MGheb0
;;; Next, for XH(surface) .NE. 0.0, we have to transform Gamma to the pure He-case
 GHE  = Gamma / (1. + XH) 
 logG4 = ALOG10(GHE/(1.-GHE)^4.)
 MGheb = FC(IMGheb,1) + FC(IMGheb,2) * SQRT(FC(IMGheb,3) + FC(IMGheb,4)*logG4)
 Mheb = 10.^MGheb

;;; Compute log(L/Lsun) using Eq. (8)
 Lhom  = 4.813 - ALOG10(1.+XH) + ALOG10(Gamma) + MGhom
 Lheb  = 4.813 - ALOG10(1.+XH) + ALOG10(Gamma) + MGheb
 Lheb0 = 4.813 + ALOG10(GHE) + MGheb0

 PRINT, ' '
 PRINT, ' Input:     Gamma_e =',  Gamma, ',  XH(surface)=', XH, FORMAT='(A,F8.3,A,F8.3)'
 PRINT, 'Output: HOM: M/Msun =', Mhom, ', HEB: M/Msun =', Mheb, ', (XH=0.0:', Mheb0,')', FORMAT='(A,F8.3,A,F8.3,A,F8.3,A)' 
 PRINT, '    --> HOM: log(L) =', Lhom, ', HEB: log(L) =', Lheb, ', (XH=0.0:', Lheb0,')', FORMAT='(A,F8.3,A,F8.3,A,F8.3,A)' 
 PRINT, ' '

 READ,  'New Gamma (g), XH (x), or Quit (q)? ', NEXT
 IF NEXT EQ 'G' OR NEXT EQ 'g' THEN BEGIN
    READ,'    Gamma_e: ', Gamma
    GOTO, REPEAT4
 ENDIF
 IF NEXT EQ 'X' OR NEXT EQ 'x' THEN BEGIN
    READ,'XH(surface): ', XH
    GOTO, REPEAT4
 ENDIF

ENDIF

END