PRO cal_hb,Z,teff,logL,Mass,logmdot

;  ------------------------------------------------
;  -                                              -
;  - MASS LOSS RECIPE for HORIZONTAL BRANCH Stars -
;  -                                              -
;  ------------------------------------------------
;
;   The Recipe is presented in:
;
;   "Jorick S. Vink, & Santi Cassisi, 2002, A&A in press, [astro-ph/0207037]"
;
;
;   The assumed Terminal Wind velocity (vinf) equals the 
;   escape velocity (vesc), which implies Mdot is maximal. 
;   As discussed in: 
;   
;   Jorick S. Vink, & Alex de Koter, 2002, A&A in press, [astro-ph/0207170]"
;   
;   If one knows vinf > vesc, one may adopt the (vinf/vesc)
;   dependence from the massive OB star computations in:
;
;   "New theoretical mass-loss rates of O and B stars"
;   Jorick S. Vink, A. de Koter, and H.J.G.L.M. Lamers, 2000, A&A 362, 295
; 
; ------------------------------------------------------------

if n_params(0) eq 0 then begin
   print,'cal_hb,Z,teff,logL,Mass'
   print,''
   print,'input:'
   print,'   Z         - Metallicity Z compared to galactic value of Z=0.02
   print,'               E.g. for Z = 0.002 insert: 0.1 
   print,' '
   print,'   teff      - effective temperature in Kelvin'
   print,' '   
   print,'   logL      - log of Luminosity (in solar units)'
   print,' '
   print,'   mass      - Evolutionary Mass (in solar units)'
   print,' '
   print,'output:'
   print,' '
   print,'   mdot      - Predicted mass-loss rate'
   return
endif
;
; ------------------------------------------------------------
;
msun      = 1.989D33
lsun      = 3.827D33
pi        = 3.1415
stefb     = 5.670D-5
grav      = 6.670D-8

lum       = 10^logL
mstar     = mass*msun
lumi      = lum*lsun

if (Z le 0) then begin
 print,'Z value needs to be higher than 0!'
 goto, Mistake
endif

logZ      = alog10(Z)

     logmdot = (-11.70) $
                + (2.13 * alog10(lum/(10^1.5)))$
                - (1.09 * alog10(mass/0.5)) $ 
                + (1.07 * alog10(teff/20000.)) $
                + 0.97   *  logZ

print,'the predicted mass loss rate is: log(mdot) = ',logmdot

Mistake:
return
End