******************************************************
* CAHK V2.0 is a simple code for converting popular  *
* stellar activity indicators into estimates of      *
* stellar age.                                       *
*                                                    *
* Most routines are based on Mamajek & Hillenbrand   *
* (2008) and references therein. Relevant            *
* references are also listed in the subroutines in   *
* this code.                                         *
*                                                    *
* The author (EEM) has made this code available as   *
* a convenience for those who would like to          *
* estimate stellar ages using rotation/activity      *
* diagnostics for solar-type stars (approximately    *
* (~F7-K2 types).                                    *
* Caveat Emptor.                                     *
* Garbage In, Garbage Out.                           *
*                                                    *
* The expression for the ratio of total-to-selective *
* absorption R = A(V)/E(B-V) comes from Olson (1975, *
* PASP, 87, 349).                                    *
*                                                    *
* If you make use of this code, please *read* and    *
* cite:                                              *
* Mamajek & Hillenbrand, 2008, ApJ, 687, 1264        *
*                                                    *
* If you run into any obvious problems/errors,       *
* please email the author at:                        *
* emamajek pas rochester edu                         *
*                                                    *
* to compile and run (on most *NIX machines):        *
* > f77 -o cahk cahk.f                               *
* > ./cahk                                           *
*                                                    *
* V1.2 notes: added code for reporting upper and     *
* lower limits on age for R'HK input, for when R'HK  *
* is outside of useful range of activity-rotation    *
* relation (i.e. logR'HK > -4.35 or < -5.00). Also,  *
* included output regarding upper and lower limits   *
* on inferred R'HK and Rx when user inputs period.   *
*                                                    *
* V1.1 E Mamajek 3/3/2009  minor edits               *
* V1.2 E Mamajek 8/27/2009 minor edits               *
* V1.3 E Mamajek 7/21/2010 minor edits               *
* V2.0 E Mamajek 11/11/2011 major changes, including *
*      standarding Av<=>E(B-V) relation, updated all *
*      colors and bolometric corrections. Includes   *
*      updated parameters from table at:             *
* http://www.pas.rochester.edu/~emamajek/            *
*           EEM_dwarf_UBVIJHK_colors_Teff.dat        *
*                                                    *
******************************************************
      program cahk
      implicit double precision (a-z)
      character*1 qual,qual2,qual0
      character*4 spt
      character*11 satstat
      integer option,option2,option3
      double precision tsogyr,tsocgyr

      tsun = 5778.0d0 ! Teff for Sun [K]
* initialize
      r   = 0.0d0
      av  = 0.0d0
      ebv = 0.0d0
      bv  = 0.0d0
      bvo = 0.0d0
      bcv = 0.0d0

      qual0 = ' '

      write(*,*)" MADE SOME CHANGES - NEED TO TEST ALL OPTIONS! "
      write(*,*)" "
      write(*,*)"CAHK.F V2.0 (11/11/2011) Eric Mamajek"
      write(*,*)" "
      write(*,*)" What age indicator for F7-K2 dwarfs do you have?"
      write(*,*)" [-1] vsin(i) = projected rotation velocity (km/s)"
      write(*,*)"  [0] log(Lx/Lbol) "
      write(*,*)"  [1] ROSAT PSPC count-rate & HR1 hardness ratio"
      write(*,*)"  [2] Rotation Period (days) [gyrochronology]"
      write(*,*)"  [4] 1 Angstrom Ca K-index " 
      write(*,*)"  [5] Solar Sunspot Number" 
      write(*,*)"  [6] Mt. Wilson S-value "
      write(*,*)"  [7] log(R HK) (unprimed) " 
      write(*,*)"[>=8] log(R'HK) (primed) " 
      read(*,*)option

      if (option .eq. -1) then                ! INPUT = VSINI 
         write(*,*)' Enter Teff, logL, vsini [km/s]:'
         read(*,*)teff, logl, vsini
         lum = 10.0d0**logl
         rsun = dsqrt(lum)*(teff/tsun)**(-2.0d0)
         write(*,288)' Radius = ',rsun,' Rsun'
         per = 50.5788d0*rsun/vsini
         call bv_teff(teff,bvo)
         write(*,288)"   Assume (B-V)o = ",bvo

         logper = dlog10(per)
         call gyro_age_barnes(bvo,per,logtb,ulogtb,tbmyr,utbmyr)
         call gyro_age_mamajek(bvo,per,uper,logtm,ulogtm,tmmyr,utmmyr)
         call gyro_age_meibom(bvo,per,logtmb,ulogtmb,tmbmyr,utmbmyr)
         tbgyr  = 10.0d0**logtb/1.0d9
         tmgyr  = 10.0d0**logtm/1.0d9
         tmbgyr = 10.0d0**logtmb/1.0d9
         write(*,*)" "
         write(*,107)per
 107     format(" upper limit on rot. period = ",f6.2," days")
         write(*,*)" "

*         write(*,*)" UPPER LIMITS TO GYROCHRONOLOGY AGE "
*         write(*,*)" Gyro age (Barnes 2007 fit) "
*         write(*,298)"          log(age/yr) < ",logtb," ; Age = ",tbgyr,
*     + " [Gyr]"
*         write(*,*)" Gyro age (Meibom+ 2009 fit) "
*         write(*,298)"          log(age/yr) < ",logtmb," ; Age = ",
*     +   tmbgyr," [Gyr]" 

         write(*,*)" Gyro age (Mamajek & Hillenbrand 2008 fit) "
         write(*,298)"          log(age/yr) < ",logtm," ; Age = ",tmgyr,
     + " [Gyr] "

         write(*,*)" " 
         p100 = 8.4275d0*(bvo-0.4d0)**0.601d0

         write(*,*)" The 100 Myr gyrochrone roughly defines the minimum"
         write(*,*)" period of applicability of the gyro calibration."
         write(*,298)" Thresh. Per. @ 100Myr = ",p100," [days] "
         write(*,*)" " 

      endif

      if (option .ne. 1) then 
       write(*,*)"What color information do you have?"
       write(*,*)" [0] Intrinsic (B-V)o (Johnson system)"
       write(*,*)" [1] Observed Bt, Vt (Tycho system), with either Av or 
     + E(B-V) (Johnson)"
       write(*,*)" [2] Observed (B-V), with either Av or E(B-V) (Johnson
     +)"
       write(*,*)" [3] Intrinsic (V-Ks)o (Johnson,2MASS)"
       write(*,*)"[>4] Spectral Type  [other] "
       read(*,*)option2

* COLOR OPTIONS [all for getting (B-V)o for the various calibrations]

        if (option2 .eq. 0) then     ! OPTION ZERO
          write(*,*)"Enter intrinsic (B-V)o color [(B-V)Sun=0.64)]:"
          read(*,*)bvo

        elseif (option2 .eq. 1) then ! OPTION ONE / TYCHO INPUT 
          write(*,*)"Enter Tycho B, error(B), V, error(V) [mag]"
          read(*,*)bt,ubt,vt,uvt
          call tychophot(bt,ubt,vt,uvt,vj,uvj,bv,ubv)

          write(*,*)"Do you have [0] E(B-V) or [other] Av ?"
          read(*,*)option3

          if (option3 .eq. 0) then ! you have E(B-V)
            write(*,*)"Enter reddening E(B-V) [mag]:"
            read(*,*)ebv
            bvo = bv - ebv 
            call Olson75_R_BVO_EBV(bvo,ebv,r,av)

          else                     ! you have Av
            write(*,*)"Enter extinction Av [mag]:"
            read(*,*)av
            call Olson75_AV_BV_IterSolv_BVO_EBV(av,bv,bvo,ebv,r)

          endif
          write(*,288)"           (B-V) = ",bv
          write(*,288)"          (B-V)o = ",bvo
          write(*,288)"          E(B-V) = ",ebv
          write(*,288)" R = A(V)/E(B-V) = ",r
          write(*,288)"              Av = ",av

        elseif (option2 .eq. 2) then ! OPTION TWO 
          write(*,*)"Do you have [0] E(B-V) or [other] Av ?"
          read(*,*)option3

          if (option3 .eq. 0) then ! you have (B-V) and E(B-V)
            write(*,*)"Enter reddening E(B-V) [mag]:"
            read(*,*)ebv
            write(*,*)"Enter observed B-V:"
            read(*,*)bv
            bvo = bv - ebv ! intrinsic B-V 
            call Olson75_R_BVO_EBV(bvo,ebv,r,av)
          else                           ! you have (B-V) and Av 
            write(*,*)"Enter extinction Av [mag]"
            read(*,*)av
            write(*,*)"Enter observed B-V:"
            read(*,*)bv
            call Olson75_AV_BV_IterSolv_BVO_EBV(av,bv,bvo,ebv,r)
          endif
          write(*,288)"           (B-V) = ",bv
          write(*,288)"          (B-V)o = ",bvo
          write(*,288)"          E(B-V) = ",ebv
          write(*,288)" R = A(V)/E(B-V) = ",r
          write(*,288)"              Av = ",av

        elseif (option2 .eq. 3) then ! OPTION2 = 3 ; V-Ks
          write(*,*)" Do you have"
          write(*,*)"     [0] V-Ks color (Johnson,2MASS)"
          write(*,*)" [other] Johnson V, 2MASS Ks magnitudes"
          read(*,*)option3

          if (option3 .eq. 0) then 
             write(*,*)" Enter intrinsic V-Ks color"
             read(*,*)vk
          else 
             write(*,*)" Enter Johnson V magnitude "
             read(*,*)vmag
             write(*,*)" Enter 2MASS Ks magnitude "
             read(*,*)ksmag
             vk = vmag - ksmag
          endif 
          call vk_bv(vk,bvo,bk)
          write(*,288)"         (V-Ks)o = ",vk 
          write(*,288)"          (B-V)o = ",bvo

        else ! option2 > 4 [SPT] 
          write(*,*)"Enter Spectral Type, e.g. 'F7.5' or 'K1'"
          write(*,*)"[No luminosity class; half-types OK]"
          read(*,*)spt
          call spt_teff_bvo_bcv(spt,teff,bvo,bcv)
          write(*,287)"           Teff = ",teff
          write(*,288)"         (B-V)o = ",bvo
          write(*,288)"            BCv = ",bcv
          if (bv .eq. -9d99) then
            write(*,*)"ERROR - invalid spectral type - assume B-V=0.6"
            bvo = 0.6d0
          endif
        endif
        if (bvo .le. 0.4 .or. bvo .ge. 0.9) then
          write(*,*)" WARNING: COLOR IS OUTSIDE OF CALIBRATION RANGE"
          write(*,*)" FOR ACTIVITY-ROTATION RELATIONS..."
        endif

      else ! for Lx/Lbol equation
        write(*,*)"Do you know:"
        write(*,*)"[1] (B-V)o, Av"
        write(*,*)"[2] Teff,Av"
        write(*,*)"[3] SpT,Av"
        write(*,*)"[ELSE] SpT,(B-V)obs"
        read(*,*)integer

        if (integer .eq. 1) then 
          write(*,*)"Av, Vmag, plx[mas], (B-V)o, ct/s, HR1 "
          read(*,*)av,v,par,bvo,cts,hr1
          call Olson75_AV_BV0_IterSolv_BV_EBV(av,bvo,bv,ebv,r)
          call teff_bv(bvo,teff,logt)
          call teff2bcv(teff,bcv)
          write(*,288)"           (B-V) = ",bv
          write(*,288)"          (B-V)o = ",bvo
          write(*,288)"          E(B-V) = ",ebv
          write(*,288)" R = A(V)/E(B-V) = ",r
          write(*,288)"              Av = ",av
          write(*,287)"            Teff = ",teff
          write(*,288)"             BCv = ",bcv
        elseif(integer .eq. 2) then
          write(*,*)"Teff, Av, Vmag, plx[mas], ct/s, HR1 "
          read(*,*)teff,av,v,par,cts,hr1
          call bv_teff(teff,bvo)
          call teff2bcv(teff,bcv)
          call Olson75_AV_BV0_IterSolv_BV_EBV(av,bvo,bv,ebv,r)
          write(*,288)"           (B-V) = ",bv
          write(*,288)"          (B-V)o = ",bvo
          write(*,288)"          E(B-V) = ",ebv
          write(*,288)" R = A(V)/E(B-V) = ",r
          write(*,288)"              Av = ",av
          write(*,287)"            Teff = ",teff
          write(*,288)"             BCv = ",bcv
        elseif(integer .eq. 3) then
          write(*,*)"Av, SpT, Vmag, plx[mas], ct/s, HR1 "
          read(*,*)av,spt,v,par,cts,hr1
          call spt_teff_bvo_bcv(spt,teff,bvo,bcv)
          call Olson75_AV_BV0_IterSolv_BV_EBV(av,bvo,bv,ebv,r)
          write(*,288)"           (B-V) = ",bv
          write(*,288)"          (B-V)o = ",bvo
          write(*,288)"          E(B-V) = ",ebv
          write(*,288)" R = A(V)/E(B-V) = ",r
          write(*,288)"              Av = ",av
          write(*,287)"            Teff = ",teff
          write(*,288)"             BCv = ",bcv
        else
          write(*,*)" V, SpT, (B-V), plx[mas], ct/s, HR1"
          read(*,*)v,spt,bv,par,cts,hr1
          call spt_teff_bvo_bcv(spt,teff,bvo,bcv)
          call Olson75_R_BVO_EBV(bvo,ebv,r,av)
          write(*,288)"           (B-V) = ",bv
          write(*,288)"          (B-V)o = ",bvo
          write(*,288)"          E(B-V) = ",ebv
          write(*,288)" R = A(V)/E(B-V) = ",r
          write(*,288)"              Av = ",av
        endif 

        if (av .lt. 0.00) then
          av = 0.00
        endif
      endif

* 0 : LOG(L_X/L_BOL) = LOG(R_X) 

      if (option .eq. 0) then
         write(*,*)"Enter Log(Lx/Lbol):"
         read(*,*)logrx
         call age_rx_mamajek(logrx,logtx) ! direct age estimate
         txgyr = 10.0d0**logtx/1.00000d9
         write(*,*)" "
         write(*,*)" X-ray age (Mamajek & Hillenbrand 2008 eqn. A3) "
         write(*,298)" log(age/yr) = ",logtx," ; Age = ",txgyr,
     + " [Gyr]"

         if (logrx .ge. -7.0d0 .and. logrx .le. -4.0d0) then 
          call xray_rossby(logrx,logro) ! MH2008 eqn. 9
          ro  = 10.0d0**logro ! Rossby number 
          uro = 0.25d0
         else
          write(*,*)' LogLx/Lbol out of range (-7 to -4) for '
          write(*,*)' fit to Rossby number (MH2008 eqn. 9)'
          ro  = 0.0d0
          uro = 0.0d0
         endif 

         call lxlbol_rhk_mamajek(logrx,lgrphk)
         write(*,*)" " 
         write(*,*)"Corresponding R'HK index from MH2008 eqn. A1"
         write(*,288)" LogR'HK = ",lgrphk," [dex]" 

* 1 : ROSAT PSPC & HR1 
      elseif (option .eq. 1) then 
         call xflux(v,av,bcv,par,cts,hr1,fx,logrx,loglx,mv,logl)
         dist = 1000.0/par
         write(*,288)"                   Dist = ",dist," [pc]"
         write(*,289)" 0.1-2.4 keV X-ray flux = ",fx," [erg/s/cm^2] " 
         write(*,288)"           log(Lx/Lbol) = ",logrx," [dex]"
         write(*,288)"                log(Lx) = ",loglx," [dex]"
         write(*,288)"                     Mv = ",mv," [mag]"
         write(*,288)"            log(L/Lsun) = ",logl," [dex]"
         call age_rx_mamajek(logrx,logtx) ! direct age estimate
         txgyr = 10.0d0**logtx/1.00000d9
         write(*,*)" X-ray age (Mamajek & Hillenbrand 2008 fit) "
         write(*,298)"            log(age/yr) = ",logtx
         write(*,298)"                    Age = ",txgyr," [Gyr]"
         call lxlbol_rhk_mamajek(logrx,lgrphk)
         write(*,*)"Corresponding Ca HK index from MH2008 eqn. A1"
         write(*,288)" LogR'HK = ",lgrphk," [dex]" 

* 2 : PERIOD (Barnes 2007 and Mamajek & Hillenbrand 2008)
      elseif (option .eq. 2) then ! ROTATION PERIOD - GYRO AGE
         write(*,*)"Enter the Rotation Period & uncertainty [days]:"
         read(*,*)per,uper
         logper = dlog10(per)
         ulogper = (dlog10(per+uper)-dlog10(per-uper))/2.0d0
         call gyro_age_barnes(bvo,per,logtb,ulogtb,tbmyr,utbmyr)
         call gyro_age_mamajek(bvo,per,uper,logtm,ulogtm,tmmyr,utmmyr)
         call gyro_age_meibom(bvo,per,logtmb,ulogtmb,tmbmyr,utmbmyr)
         tbgyr  = 10.0d0**logtb/1.0d9
         tmgyr  = 10.0d0**logtm/1.0d9
         tmbgyr = 10.0d0**logtmb/1.0d9
         utbgyr = utbmyr/1000.0d0
         utmgyr = utmmyr/1000.0d0
         utmbgyr = utmbmyr/1000.0d0
         write(*,*)
*         write(*,*)" Gyro age (Barnes 2007 fit) "
*         write(*,298)"          log(age/yr) = ",logtb," ; Age = ",tbgyr,
*     + " [Gyr] (+- ",utbgyr,")"
*         write(*,*)" Gyro age (Meibom+ 2009 fit) "
*         write(*,298)"          log(age/yr) = ",logtmb," ; Age = ",
*    +     tmbgyr," [Gyr] (+- ",utmbgyr,")" 
         write(*,*)" Gyro age (Mamajek & Hillenbrand 2008 fit) "
         write(*,298)"          log(age/yr) = ",logtm," ; Age = ",tmgyr,
     + " [Gyr] (+- ",utmgyr,")"


         if (logtb .le. 8.0d0 .or. logtm .le. 8.0d0 .or. 
     +       logtmb .le. 8.0d0) then
             write(*,*)' WARNING: Age upper limit is ~100 Myr'
             write(*,*)'       if quoted gyro age is <100 Myr!'
         endif

         write(*,*)" " 
         p100 = 8.4275d0*(bvo-0.4d0)**0.601d0

         write(*,*)" The 100 Myr gyrochrone roughly defines the minimum"
         write(*,*)" period of applicability of the gyro calibration."
         write(*,298)" Thresh. Per. @ 100Myr = ",p100," [days] "
         write(*,*)" " 

         call tau_turnover_noyes84(bvo,tauc_noyes84,lgtauc_noyes84)
         logtauc = dlog10(tauc_noyes84)
         ro = per/tauc_noyes84
         logro = dlog10(ro) 
         write(*,*)" Convective turnover times from Noyes et al (1984)"
         write(*,288)" Tau(Conv. turnover) : ",tauc_noyes84," [days] "
         write(*,288)"          log10(tau) : ",logtauc," [log10 days]"
         write(*,288)"   Rossby Number Ro  : ",ro
         write(*,288)"           log10(Ro) : ",logro
         write(*,288)"          log10(Per) : ",logper," [log10 days]" 
         call saturated_pizzolato03(bvo,per,psat,satstat)
         write(*,*)" " 
         write(*,*)" What activity should this period correlate to?"
         uro = 0.0
         call rhk_rossby_mamajek(ro,uro,lgrphk,ulgrphk,rms,unc)
         call rhk_lxlbol_mamajek(lgrphk,logrx)

         lgrphkold = lgrphk

         write(*,*)' '
         if (lgrphk .gt. -4.35d0) then
          lgrphk = -4.35
          qual  = '>'
          write(*,*)" Activity is out of useful calibration range. "
          write(*,*)" Pay attention to upper and lower limits below."
          write(*,*)" Calculating threshold values for logR'HK = -4.35"  
         elseif (lgrphk .lt. -5.00d0) then
          logrphk = -5.00
          qual  = '<'
          write(*,*)" Activity is out of useful calibration range."
          write(*,*)" Pay attention to upper and lower limits below."
          write(*,*)" Calculating threshold values for logR'HK = -5.00"
         else ! in normal range of activity-rotation relation
          qual  = ' '
         endif

         write(*,299)"      log(R'HK) : ",qual,lgrphk," [dex]"
         write(*,299)"    log(Lx/Lbol) : ",qual,logrx," [dex]"

         goto 999

* 4  1-ANGSTROM K-INDEX
      elseif (option .eq. 4) then ! K-index (measured for Sun)
         write(*,*)" Enter 1-Angstron K-index "
         read(*,*)kindex
         s = 1.475d0*kindex + 0.041d0 ! Radick et al 1998 fit
         write(*,288)"Mt. Wilson S = ",s," (fit from Radick+1998)"
         x = 0.63 - bvo
         if (x .gt. 0) then
            dlogc = 0.135*x - 0.814*x**2 + 6.03*x**3
         else
            dlogc = 0.0
         endif
         lgccf  = 1.13*bvo**3 - 3.91*bvo**2 + 2.84*bvo - 0.47 !Middlekoop 1982
         lgccfp = lgccf + dlogc ! Noyes et al 1984 p. 774
         ccf    = 10.0d0**lgccfp
         rhk    = 1.340e-4*ccf*s
         call calc_rhkphot(bvo,lgrphot,rphot)
         rphk   = rhk - rphot
         lgrphk = dlog10(rphk)
         write(*,288)"     log(R'HK) = ",lgrphk


* 5 & 6: MT WILSON S-VALUE OR SUNSPOT NUMBER
      elseif (option .eq. 5 .or. option .eq. 6) then
         if (option .eq. 6) then 
           write(*,*)"Enter Mt. Wilson S-value: "
           read(*,*)s
         else ! if option = 5, input = sunspot number
           write(*,*)"Enter Solar Sunspot Number (between 0 and ~170)"
           read(*,*)sunspot
           s = 0.171d0 + 1.2d-4*sunspot ! Baliunas+95 ApJ 450,896
           write(*,*)"Converting to Mt. Wilson S (Baliunas et al. 1995)"
           write(*,288)"               S = ",s

         endif 
         x = 0.63d0 - bvo
         if (x .gt. 0) then
            dlogc = 0.135d0*x - 0.814d0*x**2 + 6.03d0*x**3
         else
            dlogc = 0.0d0
         endif
         lgccf = 1.13*bvo**3.0d0-3.91*bvo*bvo+2.84*bvo-0.47 !Middlekoop82
         lgccfp = lgccf + dlogc ! Noyes84 p. 774
         ccf    = 10.0d0**lgccfp
         rhk    = 1.340d-4*ccf*s
         call calc_rhkphot(bvo,lgrphot,rphot)
         rphk   = rhk - rphot
         lgrphk = dlog10(rphk)
         write(*,288)"log(R'HK) = ",lgrphk

         sn = (s-0.171d0)/1.20d-4 ! Baliunas+ 1995 ApJ 450, 896
         write(*,288)"Pred. Sunspot # (NSO) = ",sn

* 7 LOG(R HK) UNPRIMED 
      elseif (option .eq. 7) then
         write(*,*)"ENTER log(R HK) "
         read(*,*)lgrhk
         call calc_rhkphot(bvo,lgrphot,rphot)
         rhk    = 10.0d0**lgrhk
         rphk   = rhk - rphot
         lgrphk = dlog10(rphk)
         write(*,*)"photosphere-corrected LogR'HK = ",lgrphk

* ELSE = >8 R'HK CHROMOSPHERIC ACTIVITY INDEX
      else
         write(*,*)"Enter log(R'HK) "
         read(*,*)lgrphk
      endif 

      lgrphk0   = lgrphk ! save original lgrphk
      lgrphkold = lgrphk ! "

      write(*,*)" Ages using chromospheric activity-age relations..." 
      write(*,*)"                               log(Age)    Age "
      write(*,*)"           ACTIVITY => AGE        [yr]    [Gyr]" 

* SODERBLOM91 - converts log(R'HK) -> log(age)
      call soderblom91(lgrphk,logtso)
      tsogyr = (10.00000d0**logtso)/1.0000d9 
      write(*,901)"Soderblom+ (1991,line): ",logtso,tsogyr

* SODERBLOM91 - constant star-formation rate - conv log(R'HK) -> log(age)
      call soderblom91csfr(lgrphk,logtsoc)
      tsocgyr = (10.00000d0**logtsoc)/1.0000d9
      write(*,901)"Soderblom+ (1991,CSFR): ",logtsoc,tsocgyr

* AGEDONAHUE93 - converts log(R'HK) -> log(age)
      call agedonahue(lgrphk,logtdo)
      tdogyr = (10.0000d0**logtdo)/1.0000d9
      write(*,901)"       Donahue (1993): ",logtdo,tdogyr

* AGELACHAUME99 - converts log(R'HK) -> log(age)
      call agelachaume(lgrphk,bvo,logtla)
      tlagyr = (10.0000d0**logtla)/1.0000d9
      write(*,901)"      Lachaume+ (1999): ",logtla,tlagyr

      set ulgrphk = 0.05d9
      call agemamajek(lgrphk,ulgrphk,logtma,ulogtma)
      tmagyr = (10.0000d0**logtma)/1.0000d9
      write(*,901)"Mamajek & Hillenbrand(2008): ",logtma,tmagyr
      write(*,*)" Note: these ages should *not* be interpreted as"
      write(*,*)" independent estimates. The MH2008 calibration is"
      write(*,*)" definitely an improvement over previous"
      write(*,*)" calibrations."

* IF ACTIVITY VALUE IS OUTSIDE OF CALIBRATED RANGE, THEN GIVE
* UPPER OR LOWER LIMIT (SECTION ADDED 8/27/2009 by EEM)

      write(*,*)' '
      if (lgrphk .gt. -4.35d0) then
       lgrphk = -4.35d0
       qual  = '<'
       qual2 = '>'
       write(*,*)" Activity is out of useful calibration range. "
       write(*,*)" Pay attention to upper and lower limits below."
       write(*,*)" Calculating threshold values for logR'HK = -4.35."

      elseif (lgrphk .lt. -5.00d0) then
       lgrphk = -5.00d0
       qual  = '>'
       qual2 = '<'
       write(*,*)" Activity is out of useful calibration range."
       write(*,*)" Pay attention to upper and lower limits below."
       write(*,*)" Calculating threshold values for logR'HK = -5.00."
      else ! in normal range of activity-rotation relation
       qual  = ' '
       qual2 = ' '
      endif

* NOW INFER ROTATIONAL PROPERTIES

      if (option .le. 1 .or. option .ge. 4) then 

        call tau_turnover_noyes84(bvo,tauc_noyes84,lgtauc_noyes84)

        if (option .eq. 0 .or. option .eq. 1) then ! X-ray data
          write(*,*)'Using Rx-Ro fit to Mt. Wilson stars...'
          call xray_rossby(logrx,logro)
          ro = 10.0d0**logro
        else ! not X-ray input data 
          call rossby_rhk_mamajek(lgrphk,ulgrphk,ro,uro,rms,unc)
          call rhk_lxlbol_mamajek(lgrphkold,logrx) 
        endif 

        per = ro*tauc_noyes84
        call gyro_age_mamajek(bvo,per,uper,logtm,ulogtm,tmmyr,utmmyr)    
        tmgyr = 10.0d0**logtm/1.0d9

        write(*,*)" " 
        write(*,*)" ACTIVITY => ROTATION => AGE "
        write(*,*)" (improvement over ACTIVITY => AGE !)"
        write(*,299)"logR'HK : ",qual0,lgrphkold," [dex]" 
        write(*,299)"Rossby Number : ",qual,ro," = Per/Tau_conv."
        write(*,299)"Convect. Turnover time : ",qual0,tauc_noyes84," [da
     +y](from Noyes+84)"
        write(*,299)"Rotation Period : ",qual,per," [days]" 
        write(*,299)"log(Lx/Lbol) : ",qual0,logrx," [dex]"
        write(*,299)"Gyro age log(age/yr) : ",qual,logtm," [dex] *BEST*"
        write(*,299)"Gyro age : ",qual,tmgyr," [Gyr] *BEST*" 
      endif 

 287  format(a25,f7.1,a)
 288  format(a25,f7.3,a)
 289  format(a25,d12.6,a)
 296  format(a25,d7.3,a30) 
 297  format(a25,f7.1,a30)
 298  format(a25,f7.3,a8,f5.2,a11,f5.2,a1)
 299  format(a25,a1,f6.2,a)
 900  format(3(f6.2,1x))
 901  format(a30,1x,f8.2,1x,f8.2)
 902  format(f6.3)
 999  end

*********************************************************
* LXLBOL_RHK_MAMAJEK - convert Lx/Lbol into R'HK index  *
* using fit to Donahue-Baliunas sample of stars with    *
* accurate periods, logR'HK, and logR_X = logLx/Lbol.   *
* The sample is X_ray-unbiased as *all* of the          *
* Donahue-Baliunas stars are detected in RASS.          *
* I"ve also confirmed that the fit does a good job      *
* of fitting the high activity end. Fit is published in *
* appendix of Mamajek & Hillenbrand (2008)              * 
* INPUT:                                                *
* logrx = log(Lx/Lbol)                                  *
* OUTPUT:                                               *
* lgrphk = log(R'HK)                                    *
*                                                       *
* E. Mamajek 7/10/2008                                  *
* 4/8/2008 - final value for Mamajek & Hillenbrand 2008 *
*********************************************************
      subroutine lxlbol_rhk_mamajek(logrx,lgrphk)
      implicit double precision (a-z)
      b  = -4.5356
      xo =  4.920 
      m  =  0.2890
      lgrphk = b + m*(logrx + xo) 
      end

*********************************************************
* RHK_LXLBOL_MAMAJEK - convert R'HK into logLx/Lbol     *
* using fit to Donahue-Baliunas sample of stars with    *
* accurate periods, logR'HK, and logR_X = logLx/Lbol.   *
* The sample is X_ray-unbiased as *all* of the          *
* Donahue-Baliunas stars are detected in RASS.          *
* I"ve also confirmed that the fit does a decent job    *
* of fitting the high activity end. Fit is published in *
* appendix of Mamajek & Hillenbrand (2008)              * 
* INPUT:                                                *
* lgrphk = log(R'HK)                                    *
* OUTPUT:                                               *
* logrx = log(Lx/Lbol)                                  *
* E. Mamajek 7/10/2008                                  *
*********************************************************
      subroutine rhk_lxlbol_mamajek(lgrphk,logrx)
      implicit double precision (a-z)
      b  = -4.9005
      xo =  4.53
      m  =  3.4601
      logrx = b + m*(lgrphk + xo) 
      end

***************************************************
* AGEDONAHUE - convert R'HK to age using relation *
* of Donahue (1993; PhD thesis) as reported in    *
* Henry et al. (1996; AJ 111, 439).               *
* INPUT:                                          *
* lgrphk = log10(R'HK)                            *
* OUTPUT:                                         *
* logage = log10(age) [yr]                        *
*                                                 *
* E. Mamajek 2/17/2003                            *
***************************************************
      subroutine agedonahue(lgrphk,logage)
      implicit double precision (a-z)
      r5 = 1e5*10**(lgrphk)
      x  = r5
      a0 = 10.725
      a1 = -1.334
      a2 =  0.4085
      a3 = -0.0522
      logage = a0 + a1*x + a2*x*x + a3*x*x*x
      if (logage .gt. 10) then
         write(*,*)"WARNING: older than Galactic disk!"
      elseif (logage .lt. 6) then
         write(*,*)"WARNING: younger than T Tauri stars!"
      endif
      end

***************************************************
* SODERBLOM91 - calculate log(age) according to   *
* the third relation listed by Soderblom+ 1991    *
* (ApJ 375, 722) on p. 734 of his paper.          *
*                                                 *       
* INPUT:                                          *
* lgrphk = log10(R'HK)                            *
* OUTPUT:                                         *
* logage = log10(age) [yr]                        *
*                                                 *
* Note that Soderblom had 2 alternative relations *
* based on subsets of data:                       *
*  (1)    logts1 = -2.02*lgrphk - 0.31            *
*  (2)    logts2 = -1.34*lgrphk + 2.99            *
* The one in the subroutine is his 3rd and final. *
*                                                 *
* E. Mamajek 2/20/2003                            *
***************************************************
      subroutine soderblom91(lgrphk,logage)
      implicit double precision (a-z)
      logage = -1.50*lgrphk + 2.25
      if (logage .gt. 10) then
         write(*,*)"WARNING: older than Galactic disk!"
      elseif (logage .lt. 6) then
         write(*,*)"WARNING: younger than T Tauri stars!"
      endif
      end

***************************************************
* SODERBLOM91CSFR - calculate log(age) according  *
* to an analytic fit to the activity-age relation *
* from Soderblom et al 1991 (ApJ 375, 722) which  *
* assumed a constant star-formation history (and  *
* took into account disk heating). The analytical *
* form was not published, and was provided by DRS *
* (priv. comm.).                                  *
*                                                 *       
* INPUT:                                          *
* lgrphk = log10(R'HK)                            *
* OUTPUT:                                         *
* logage = log10(age) [yr]                        *
*                                                 *
* E. Mamajek 6/30/2008                            * 
***************************************************
      subroutine soderblom91csfr(lgrphk,logage)
      implicit double precision (a-z)
      logage = -8568.2 - 7037.6*lgrphk - 2164.6*lgrphk**2 
     + - 295.76*lgrphk**3 - 15.144*lgrphk**4
      if (logage .gt. 10) then
         write(*,*)"WARNING: older than Galactic disk!"
      elseif (logage .lt. 6) then
         write(*,*)"WARNING: younger than T Tauri stars!"
      endif
      end

***************************************************
* AGELACHAUME - calculate log(age) from Ca R'HK   *
* according to  eqn. 6 of Lachaumme+ (1999).      *
*                                                 *       
* INPUT:                                          *
* lgrphk = log10(R'HK)                            *
* bvo    = unreddened (B-V) color                 *
* OUTPUT:                                         *
* logage = log10(age) [yr]                        *
*                                                 *
* E. Mamajek 2/20/2003                            *
***************************************************
      subroutine agelachaume(lgrphk,bvo,logage)
      implicit double precision (a-z)
      alpha = 9.36d0
      beta  = 5.0d0
      gamma = -12.5d0
      x = -1*lgrphk - 4.6d0

      if (bvo .le. 0.6) then 
         write(*,*)"B-V is too blue for Lachaumme calibration"
      endif
      if (x .lt. 0.0d0) then
         logage = alpha + beta*x
      elseif (x .ge. 0.0d0 .and. x .lt. 0.2d0) then
         logage = alpha + beta*x + gamma*x*x
      else
         logage = alpha + 0.2*beta + 0.04*gamma
      endif
      if (logage .gt. 10) then
         write(*,*)"WARNING: older than Galactic disk!"
      elseif (logage .lt. 6) then
         write(*,*)"WARNING: younger than T Tauri stars!"
      endif
      end

***************************************************
* AGEMAMAJEK - convert R'HK to age using relation *
* from Mamajek & Hillenbrand (2008).              * 
* INPUT:                                          *
* lgrphk = log10(R'HK)                            *
* OUTPUT:                                         *
* logage = log10(age) [yr]                        *
* E Mamajek 7/10/2008                             * 
***************************************************
      subroutine agemamajek(lgrphk,ulgrphk,logage,ulogage)
      implicit double precision (a-z)
      a0 =  -38.053d0
      a1 =  -17.912d0
      a2 =   -1.6675d0 
      logage = a0 + a1*lgrphk + a2*lgrphk**2
      ulogage = dsqrt((a1+2.0d0*a2*lgrphk)**2.0*ulgrphk**2.0d0) 
      if (logage .gt. 10) then
         write(*,*)"WARNING: older than Galactic disk!"
      elseif (logage .lt. 6) then
         write(*,*)"WARNING: younger than T Tauri stars!"
      endif
      end

****************************************************
* AGE_RX_MAMAJEK - estimate age directly from R_X. *
* from Mamajek & Hillenbrand (2008; Appendix)      *
* INPUT:                                           *
* logrx = log(Lx/Lbol) [dex]                       *
* OUTPUT:                                          *
* logage = log10(age) [yr]                         *
* E Mamajek 7/7/2008                               *
****************************************************
      subroutine age_rx_mamajek(logrx,logage)
      implicit double precision (a-z)
      a0 =  1.20d0
      a1 = -2.307d0
      a2 = -0.1512d0
      logage = a0 + a1*logrx + a2*logrx*logrx
      end

***********************************************
* CALC_RHKPHOT calculates the contribution to *
* R_HK from the star"s photosphere.           *
* From Noyes et al. 1984 ApJ 279, 763         *
* INPUT                                       *
* bvo = intrinsic B-V color [mag]             *
* OUTPUT                                      *
* lgrphot = log(R_phot)                       *
* rphot   = R_phot                            *
*                                             *
* E Mamajek 7/10/2008                         *
***********************************************
      subroutine calc_rhkphot(bvo,lgrphot,rphot)
      implicit double precision (a-z)
      a0 = -4.898
      a1 = +1.918
      a2 = -2.893
      lgrphot = a0 + a1*bvo**2 + a2*bvo**3
      write(*,*)"Using Noyes et al. (1984) photospheric correction"
      if (bvo .lt. 0.44) then
         write(*,*)"Too blue for photospheric subtraction calibration"
      endif
      if (bvo .gt. 0.9) then
         write(*,*)"Too red for photospheric subtraction calibration"
      endif
      rphot = 10.0d0**lgrphot
      return
      end


****************************************************************
* TYCHOPHOT - convert Tycho photometry to Johnson photometry   *
* INPUT:                                                       *
* vt, uvt     = Tycho V mag(& uncertainty)                     *
* bt, ubt     = Tycho B mag(& uncertainty)                     *
* OUTPUT:                                                      *
* btvt, ubtvt = (B-V)_Tycho(& uncertainty)                     *
* bvj, ubvj   = (B-V)_Johnson(& uncertainty)                   *
* vj, uvj     = V_Johnson(& uncertainty)                       *
*                                                              *
* The polynomial fits are to the data of Bessel 2000 (PASP     *
* 112, 961; Table 2.), and are given in Mamajek+ 2002 (AJ 124, *
* 1670). The output uncertainties assume the color conversions *
* have no uncertainty (i.e. the uncertainties are solely from  *
* from the Bt,Vt errors). The B-V color conversion is an       *
* improvement on the HIP-Tyc Catalog"s linear regression fit   *
* (B-V)_Johnson = 0.90*(B-V)_Tycho. This linear fit tends to   *
* be too blue in Johnson B-V by ~0.05 mags amongst FGKs, & the *
* new polynomial fit removes this systematic effect. A new     *
* subroutine for calculating the color error was added 10/02.  *
* This follows the equations in color error given in Vol. 4 ,  *
* Sec. 9.5 of the Hipparcos & Tycho catalog. Note that the     *
* polynomial in Mamajek+ 2002 is in error (one of the terms has*
* the wrong sign), and the corrected polynomial is used here). *
* The corrected polynomial fit is also given at:               *
* http://www.aerith.net/astro/color_conversion.html            * 
* and in an erratum by the author:                             *
* Mamajek, Meyer, Liebert 2006, AJ, 131, 2360                  *
* 2006AJ....131.2360M                                          * 
*                                                              *
* written 8/14/2000 Eric Mamajek                               *
* updated 2/2/2002  "            better V-Vt fit               *
* updated 10/4/2002 "            added B-V error subroutine    *
****************************************************************
      subroutine tychophot(bt,ubt,vt,uvt,vj,uvj,bvj,ubvj)
      implicit double precision (a-z)
      character*32 dummy,filen1,filen2
      dif = 0.0033401186 ! polyfit dif at (B-V)t = 0.5
                         ! used to "glue" the fits together smoothly
      a4 =  0.0000000
      a3 = -0.01997773   ! Bessell 2000 Bt-Vt -> V-Vt polyfit(3)
      a2 =  0.05486121   ! IRAF polyfit 2/21/02 (good for -0.25 < Bt-Vt < 2)
      a1 = -0.1334168    !
      a0 =  9.669378E-4
      call colorerror(bt,ubt,vt,uvt,btvt,ubtvt)
      if (btvt .ge. 0.50) then ! Bessell 2000 Bt-Vt -> (B-V)-(Bt-Vt)
         b0 = -2.482282E-4     ! new piecemeal fit (2/21/02)
         b1 = -0.007813261     ! DOUBLE-CHECKED, this term is correct
         b2 = -0.1489144
         b3 =  0.03383945
         b4 =  0.000000
         b5 =  0.0000000
      else
         b0 = -0.006015661
         b1 = -0.1068928
         b2 = 0.1458791
         b3 = 0.0000000
         b4 = 0.0000000
         b5 = 0.0000000
      endif
      x2 = btvt*btvt
      x3 = x2*btvt
      x4 = x3*btvt
      x5 = x4*btvt
      vj  = vt + a4*x4 + a3*x3 + a2*x2 + a1*btvt + a0
      uvj = dsqrt((4*a4*x3 + 3*a3*x2 + 2*a2*btvt +
     +  a1)**2*ubtvt*ubtvt + uvt*uvt)
      bvj = btvt + b5*x5 + b4*x4 + b3*x3 + b2*x2 +
     +      b1*btvt + b0
      ubvj = (5*b5*x4 + 4*b4*x3 + 3*b3*x2 +
     +       2*b2*btvt + b1 + 1)*ubtvt
      junkv = vj - vt
      junk = bvj-btvt
      bj  = vj + bvj
      ubj = dsqrt(ubvj**2 - uvj**2)
      write(*,997)" (B_T - V_T) = ",btvt," +- ",ubtvt," [mag] (Tycho)"
      write(*,996)" (V_J - V_T) = ",junkv,"    (Johnson minus Tycho)"
      write(*,996)" (B-V)_J - (B-V)_T = ",junk,"  (Johnson minus Tycho)"
      write(*,997)" V = ",vj," +- ",uvj," [mag] (Johnson)"
      write(*,997)" B = ",bj," +- ",ubj," [mag] (Johnson)"
      write(*,997)" (B-V)_J = ",bvj," +- ",ubvj," [mag] (Johnson)"
      write(*,*)" "
 996  format(a20,f7.3,a29)
 997  format(a20,f7.3,a4,f5.3,a20)
      return
      end

***********************************************************
* COLORERROR - calculates the uncertainty in a color      *
* using the formulae from the Hipparcos & Tycho catalogs  *
* S1.3, p. 48, formula 1.3.12. Note that the input        *
* magnitude errors are understood to the B and V errors   *
* towards brighter magnitudes - the default values given  *
* in the Hipparcos and Tycho catalogs (fields T33 & T35). *
* INPUT:                                                  *
* b, ub = magnitude1 (and unc.) [mag]                     *
* v, uv = magnitude2 (and unc.) [mag]                     *
* OUTPUT:                                                 *
* bv,ubv = color (and unc.) [mag]                         *
* Eric Mamajek 10/4/2002                                  *
* updated 10/15/2002 using eqns from S1.3 of ESA Hip/Tyc  *
***********************************************************
      subroutine colorerror(b,ub,v,uv,bv,ubv)
      implicit double precision (a-z)
      n = 1.0d0  ! number of sigma, the formula is generalized for any n
      bv = b - v
      ubv = 2.5d0*dsqrt((dlog10(1.0d0+n-n*10**(0.4d0*ub)))**2.0d0 +
     +                (dlog10(1.0d0-n+n*10.0d0**(0.4d0*uv)))**2.0d0)
      return
      end

********************************************************
* GYRO_AGE_BARNES calculates the gryochronology age via
* Barnes 2007. Note that the original formula gave gyro
* age in log(age/Myr), so I add 6 to get log(age/yr).
* INPUT
* bv  = intrinsic B-V color [mag]
* per = rotation period [days]
* OUTPUT
* logt = gyro age [yr]
* E Mamajek 4/27/2007
********************************************************
      subroutine gyro_age_barnes(bv,per,logt,ulogt,tmyr,utmyr)
      implicit double precision (a-z)
* apologies for nomenclature here -- "bv" is the intrinsic stellar
* color, whereas "bvo" is the color singularity in the gyro equations.-EEM
* constants from Barnes paper
      n = 0.5189d0
      a = 0.7725d0
      b = 0.601d0 
      bvo = 0.40d0 
* equation 3 from Barnes (2007)
      logp = dlog10(per)
      loga = dlog10(a)
      x    = bv - bvo 
      logx = dlog10(x)      
      logt = 6.0d0 + (logp - loga - b*logx)/n 
      tmyr = 10.0d0**((logp - loga - b*logx)/n)
* errors
      c = 4.5 
      d = 0.5*dlog(tmyr)**2
      e = 2.0d0*per**0.6d0
      f = (0.6d0/x)**2.0
      frac = 0.02d0 * dsqrt(c + d + e + f)
      utmyr = frac*tmyr 
      logt1 = dlog10(6.0d0 + tmyr + utmyr)
      logt2 = dlog10(6.0d0 + tmyr - utmyr)
      ulogt = dabs(logt1 - logt2)/2.0d0
      return
      end

**********************************************************
* GYRO_AGE_MAMAJEK calculates the gryochronology age via *
* new gyrochronology parameters from Mamajek &           *
* Hillenbrand (2008). The formalism is the same as that  *
* of Barnes (2007).                                      *
* INPUT                                                  *
* bv  = intrinsic B-V color [mag]                        *
* per = rotation period [days]                           *
* OUTPUT                                                 *
* logt = gyro age [yr]                                   *
* E Mamajek 1/28/2008                                    *
**********************************************************
      subroutine gyro_age_mamajek(bv,per,uper,logt,ulogt,tmyr,utmyr)
      implicit double precision (a-z)
* apologies for nomenclature here -- "bv" is the intrinsic stellar
* color, whereas "bvo" is the color singularity in the gyro equations.-EEM
* parameters from Mamajek & Hillenbrand (2008)
      n   = 0.5657d0
      a   = 0.4072d0
      b   = 0.3257d0
      bvo = 0.4952d0
* calculated constants
      loga = dlog10(a)
      x    = bv - bvo
      logx = dlog10(x)
* equation 3 from Barnes (2007)
      logp   = dlog10(per)
      logphi = dlog10(per+uper)
      logplo = dlog10(per-uper)

      logt   = 6.0d0 + (logp   - loga - b*logx)/n
      logthi = 6.0d0 + (logphi - loga - b*logx)/n
      logtlo = 6.0d0 + (logplo - loga - b*logx)/n
      ulogt2 = dabs(logthi - logtlo)/2.0d0

      tmyr   = 10.0d0**((logp   - loga - b*logx)/n)
      tmyrhi = 10.0d0**((logphi - loga - b*logx)/n)
      tmyrlo = 10.0d0**((logplo - loga - b*logx)/n)
      utmyr2 = dabs(tmyrhi - tmyrlo)/2.0d0
*      write(*,*)' logt = ',tmyr,' +- ',utmyr2,' Myr '

* errors
      c = 4.5d0
      d = 0.5*dlog(tmyr)**2.0d0
      e = 2.0d0*per**0.6d0
      f = (0.6d0/x)**2.0d0
      frac = 0.02d0 * dsqrt(c + d + e + f)
      utmyr = frac*tmyr
      logt1 = dlog10(6.0d0 + tmyr + utmyr)
      logt2 = dlog10(6.0d0 + tmyr - utmyr)
      ulogt = dabs(logt1 - logt2)/2.0d0
      if (bv .le. bvo .or. bv .ge. 1.4) then
         write(*,*)" WARNING: B-V color out of range for gyro age"
      endif 
      end

********************************************************
* TAU_TURNOVER_NOYES84 - calculate convective turnover * 
* time from (B-V) color following Noyes et al (1984;   *
* ApJ 279, 763).                                       *
* INPUT:                                               *
* bvo   = unreddened B-V color [mag] (0.4-1.4)         *
* OUTPUT:                                              *
* lgtauc = log10 of convective turnover time for       *
*          alpha = 1.9 (mixing length/scale ht.)       *
* E. Mamajek 7/10/2008                                 *
********************************************************
      subroutine tau_turnover_noyes84(bvo,tauc,lgtauc)
      implicit double precision (a-z)
      a0 =  1.362d0
      a1 = -0.166d0
      a2 =  0.025d0
      a3 = -5.323d0
      b0 =  1.362d0
      b1 = -0.140d0
      x = 1.00d0 - bvo
      if (x .gt. 0) then
         lgtauc = a0 + a1*x + a2*x*x + a3*x*x*x
      else                     ! actually, if x < 0
         lgtauc = b0 + b1*x
      endif
      if (bvo .gt. 1.4d0 .or. bvo .lt. 0.4d0) then
         write(*,*)"B-V outside of range for calculation
     + of convective turnover parameter time"
      endif
      tauc = 10.0d0**lgtauc
      return
      end

***************************************************************
* SATURATED_PIZZOLATO03 tests what rotation/saturation        *
* regime the star is probably in, given their cluster         *
* study (Pizzolato et al. 2003 A&A 397, 147)                  *
* INPUT: bvo, per = (B-V)o color [mag] & rotation period [day]*
* OUTPUT: psat = rotation period for saturated Lx/Lbol [day]  *
*        satstat = "saturated" or "unsaturated"               *
* E Mamajek 4/29/2007                                         *
***************************************************************
      subroutine saturated_pizzolato03(bvo,per,psat,satstat)
      implicit double precision(a-z)
      character*11 satstat
      if (bvo .ge. 0.50 .and. bvo .lt. 0.60) then
         psat = 1.6d0
      elseif (bvo .ge. 0.60 .and. bvo .lt. 0.71) then
         psat = 1.3d0
      elseif (bvo .ge. 0.71 .and. bvo .lt. 0.84) then
         psat = 1.7d0
      elseif (bvo .ge. 0.84 .and. bvo .lt. 0.91) then
         psat = 2.4d0
      elseif (bvo .ge. 0.91 .and. bvo .lt. 1.06) then
         psat = 2.4d0
      elseif (bvo .ge. 1.06 .and. bvo .lt. 1.36) then
         psat = 3.5d0
      elseif (bvo .ge. 1.36 .and. bvo .lt. 1.54) then
         psat = 3.7d0
      else
         write(*,*)" outside Pizzolatto color range"
         write(*,*)" or poorly defined..."
         write(*,*)" ignore saturation status!"
      endif 

      if (per .lt. psat) then
         satstat = "SATURATED!"
      else
         satstat = "UNSATURATED"
      endif
      write(*,*)" Period for X-ray saturation for this color is: "
      write(*,998)" P(sat) = ",psat," [day] (Pizzolato et al. 2003)"
      write(*,*)" Rotation-activity Status: ",satstat
 998  format(a16,f6.2,a)
 999  format(a10,f5.1,a30)
      return
      end


**********************************************************
* RHK_ROSSBY_MAMAJEK estimates chromospheric             *
* activity logR'HK from the Rossby Number.               *
* Fits from Mamajek & Hillenbrand (2008)                 *
*                                                        *
* INPUT                                                  *
* ro,uro = rossby number (period/turnover time) & unc.   *
* OUTPUT                                                 *
* lgrphk = logR'HK                                      *
* ulgrphk = uncertainty in logR'HK from calibration     *
* rms = rms scatter in logR'HK as f(Ro)                  *
* unc = total uncertainty in predicted logR'HK value     *
*                                                        *
* E Mamajek 7/10/2008                                    *
**********************************************************
      subroutine rhk_rossby_mamajek(ro,uro,lgrphk,ulgrphk,rms,unc)
      implicit double precision(a-z)
      if (ro .ge. 0.3168) then ! INACTIVE SEQUENCE
       b   = -4.5220d0
       ub  =  0.0049d0
       m   = -0.3372d0
       um  =  0.0112d0
       c   = -0.8143d0
       uc  =  0.0000d0 ! assume no error in constant
       rms =  0.0374d0 ! from fit
      else !                     ACTIVE SEQUENCE
       b   = -4.2328d0
       ub  =  0.0221d0
       m   = -1.4508d0
       um  =  0.1307d0
       c   = -0.2330d0
       uc  =  0.0000d0 ! assume no error in constant
       rms =  0.16d0 ! from fit      
      endif
      x  = ro + c
      ux = dsqrt(uro*uro + uc*uc)
      lgrphk = b + m*x 
      ulgrphk = dsqrt(ub*ub + x*x*um*um + m*m*ux*ux)
      unc = dsqrt(ulgrphk**2 + rms**2)
      return
      end 

***************************************************
* XFLUX  - calculates X-ray flux, Lx/Lbol,        *
* log(Lx) and log(L/Lsun).                        * 
*                                                 *
* INPUT                                           *
* v        = reddened V magnitude [mag]           *
* av       = reddening at V [mag]                 *
* bcv      = BC [mag] (should be negative, and    *
*                  on scale where BC(Sun)=-0.07)  *
* par      = parallax [mas]                       *
* ct       = PSPC X-ray counts [ct/s]             * 
* hr1      = Hardness ratio HR1 [unitless]        *
* OUTPUT                                          *
* fx       = X-ray flux [erg/s/cm^2]              *
* fxfbol   = X-ray/bolometric flux ratio          *
* loglx    = log(X-ray luminosity) [erg/s]        * 
* logl     = log(L/Lsun) [dex]                    * 
*                                                 *
* CONSTANTS                                       *
* I follow the formalism of Bessell+ 1998 for     *
* the solar parameters: m_v=-26.76, M_v=4.81,     *
* M_bol = 4.74, m_bol => -26.83, f_sun =          *
* 1.371e-6 erg/s/cm2. These values lead to a      *
* flux for a zero-bolometic-magnitude object of   *
* 2.54d-5 erg/s/cm^2.                             *
* X-ray energy Conversion Factor (CF) equation    *
* from Fleming et al 1995 ApJS 99,701 (p. 706)    *
* Relation applies from soft coronae of M dwarfs  *
* to harder coronae of RS CVns                    *
*                                                 *
* Eric Mamajek 2/19/2000                          *
* updated 3/18/2004 - input V and Av instead of Vo*
***************************************************
      subroutine xflux(v,av,bcv,par,ct,hr1,fx,fxfbol,loglx,mv,logl)
      implicit double precision (a-z)
* CONSTANTS
      pi   = 3.14159265358d0
      pccm = 3.08567758131d18 ![cm/pc]
      fo   = 2.541192d-5      ! flux [erg/s/cm^2] for apparent m_bol = 0
* CALCULATIONS
      shell  = 4.00000d0*pi*(1000.0000d0/par)**2*pccm*pccm
      mbol   = v - av + bcv
      mv     = (v - av) + 5.0d0 + 5.0d0*dlog10(par/1000.0d0)
      fbol   = fo*10.00000d0**(-0.40000d0*mbol)
* ENERGY CONVERSION FACTOR (USE FLEMING RELATION UP TO A POINT...)
      ecf    = (8.31000d0+5.3000d0*hr1)*1.000000d-12
      if (ecf .gt. 1d-11) ecf = 1.0000d-11
      fx     = ct*ecf
      fxfbol = dlog10(fx/fbol)
      loglx  = dlog10(shell*fx)
      logl   = dlog10(shell*fbol) - 33.58451214d0 ! log(L/Lsun) 
      return
      end

*************************************************
* TEFF_BVO converts Teff to (B-V) from          *
* polynomial fits to the EEM's table of colors, *
* Teff, BCv, which were constructed for each    *
* dwarf spectral subtype (November 2011         *
* calibration).                                 *
* E Mamajek 12/20/2007                          *
* 11/12/2011: updated calibration               * 
*************************************************
      subroutine teff_bv(bvo,teff,logt)
      implicit double precision (a-z)
      x = bvo 

      if (bvo.ge.0.294d0.and.bvo.lt.1.431d0) then ! F/G/K 
       logt =   3.980289674446d+00
     +        - 4.879651194352d-01*x
     +        + 2.896561027685d-01*x*x
     +        - 9.980740833102d-02*x*x*x ! rms(logt) = 0.0021 dex
      elseif (bvo.ge.-0.32d0.and.bvo.lt.0.294d0) then ! O9/B/A
       logt = 3.982132981567d+00
     +        - 6.201546711048d-01*x
     +        + 2.363160147690d+00*x**2.0d0
     +        - 1.847467702212d+01*x**3.0d0
     +        + 4.574031907113d+01*x**4.0d0
     +        + 3.010467522731d+02*x**5.0d0
     +        - 1.218478000884d+03*x**6.0d0
     +        - 1.641608134464d+03*x**7.0d0
     +        + 7.490965463874d+03*x**8.0d0
      elseif(bvo.ge.1.431d0.and.bvo.le.2.2d0) then ! Mstars
        logt =  -1.481780822244d+02   ! rms = 0.00377 dex
     +         + 4.424977211837d+02*x
     +         - 5.100804625843d+02*x**2.0d0
     +         + 2.908361843212d+02*x**3.0d0
     +         - 8.211083043327d+01*x**4.0d0
     +         + 9.188300215523d+00*x**5.0d0
      else 
        write(*,*)"Out of color range: B-V = ",bvo
      endif 
      teff = 10.0d0**logt
      end 

****************************************************************
* TEFF2BCV - convert temperature to bolometric correction      *
* (for V-band, for dwarfs). Scale BCv(Sun)=-0.07 (within <0.01 *
* mag.                                                         *
* 2011 update now uses BCv from table in file:                 *
* http://www.pas.rochester.edu/~emamajek/                      *
*        EEM_dwarf_UBVIJHK_colors_Teff.dat                     *
* INPUT:                                                       *
* temp - Temperature, linear scale (K)                         *
* OUTPUT:                                                      *
* bcv   - bolometric correction                                *
* CONSTANTS:                                                   *
*                                                              *
* Eric Mamajek 8/6/2001                                        *
* 1/23/2008  made into subroutine                              *
* 11/11/2011 updated using current calibration                 *
****************************************************************
      subroutine teff2bcv(teff,bcv)
      implicit double precision (a-z)
      x = dlog10(teff)
      if (x .ge. 4.0d0) then ! OB stars
         bcv =  -2.917533499083d+03
     +         + 2.648732019608d+03*x
     +         - 8.961693571914d+02*x**2.0d0
     +         + 1.340624628847d+02*x**3.0d0
     +         - 7.495867840331d+00*x**4.0d0 ! rms = 0.016 mag
      elseif (x .lt. 4.0d0 .and. x .ge. 3.365) then 
         bcv =  -5.723152145608d+03
     +         + 5.771910900332d+03*x
     +         - 2.194565588806d+03*x**2.0d0
     +         + 3.729814492323d+02*x**3.0d0
     +         - 2.391602736625d+01*x**4.0d0 ! rms = 0.0252 mag
      else
         bcv = 9.99d99
         write(*,*)' Teff out of range for BCv calibration'
      endif
      return
      end

********************************************
* BV_TEFF estimates B-V from Teff via a    *
* polynomial fit to SpT vs. B-V vs. Teff   *
* table:                                   *
* http://www.pas.rochester.edu/~emamajek/  *
*       EEM_dwarf_UBVIJHK_colors_Teff.dat  *
*                                          *
* E Mamajek 12/20/2007                     *
* 11/11/2011 updated using table           *
********************************************
      subroutine bv_teff(teff,bvo)
      implicit double precision (a-z)
      x = dlog10(teff)
      if (x .ge. 3.60 .and. x .le. 4.519) then ! hotter than K7 or Teff > 4000K 
         bvo =  -1.799063006760d+04
     +          + 2.230395497671d+04*x
     +          - 1.101669663126d+04*x*x 
     +          + 2.711100828313d+03*x*x*x
     +          - 3.325161422798d+02*x*x*x*x
     +          + 1.626525488477d+01*x*x*x*x*x !rms = 0.008877 mag
      elseif (x .ge. 3.365 .and. x .lt. 3.60) then ! K7-M8
         bvo =   9.114227053887d+04
     +          - 3.798566330649d+04*x
     +          - 8.542687609443d+03*x**2.0d0
     +          + 4.760097798496d+03*x**3.0d0
     +          - 1.185194830297d+03*x**4.0d0
     +          + 2.941866441422d+02*x**5.0d0
     +          + 8.612517453309d+01*x**6.0d0
     +          - 1.793258057138d+01*x**7.0d0
     +          - 1.041643131032d+00*x**8.0d0
     +          - 2.146491294267d+00*x**9.0d0
     +          - 2.889177185003d-01*x**10.0d0
     +          + 2.164615211383d-01*x**11.0d0
     +          + 2.501242140201d-02*x**12.0d0
     +          - 8.321147979327d-03*x**13.0d0 ! fit rms = 0.0159 mag
      elseif (logt .lt. 3.365 .or. logt .gt. 4.519) then 
        write(*,*)"Out of Teff range: Teff = ",teff
      endif 
      return
      end 

********************************************************
* BCV_SPT outputs bolometric correction to V-band,     *
* B-V, and Teff as a function of spectral type.        *
* Input spectral types can be normal subtypes (e.g.    *
* "M2") or half subtypes (e.g. "M2.5").                *
* Bolometric corrections are on system where BCv       *
* for solar Teff (5778K) is BCv = -0.07 and Mbol(sun)  *
* = 4.75 (which follows IAU convention to ~0.001 mag). *
* INPUT                                                *
* spt = spectral type (4 char.; e.g."B2  " or "M2.5")  *
* OUTPUT                                               *
* teff = effective temperature [K]                     *
* bvo  = unreddened B-V color: (B-V)o [mag]            *
* bvc  = bolometric correction to V [mag]              *
*                                                      *
* E Mamajek 3/15/2005                                  *
* 11/12/2011: updated Teff, B-V, & BCv to new values   *
********************************************************
      subroutine spt_teff_bvo_bcv(spt,teff,bvo,bcv)
      implicit double precision (a-z)
      character*4 spt
      bvo  = -9.999d99
      teff = -9.999d99
      bcv  = -9.999d99
* BCv values ; revised BCvs 11/5/2011
      if(spt.eq."O9  ") bcv = -3.09
      if(spt.eq."O9.5") bcv = -3.02
      if(spt.eq."B0  ") bcv = -2.99
      if(spt.eq."B0.5") bcv = -2.79
      if(spt.eq."B1  ") bcv = -2.57
      if(spt.eq."B1.5") bcv = -2.39
      if(spt.eq."B2  ") bcv = -2.01
      if(spt.eq."B2.5") bcv = -1.72
      if(spt.eq."B3  ") bcv = -1.58
      if(spt.eq."B3.5") bcv = -1.535
      if(spt.eq."B4  ") bcv = -1.49
      if(spt.eq."B4.5") bcv = -1.41
      if(spt.eq."B5  ") bcv = -1.33
      if(spt.eq."B5.5") bcv = -1.225
      if(spt.eq."B6  ") bcv = -1.12
      if(spt.eq."B6.5") bcv = -1.07
      if(spt.eq."B7  ") bcv = -1.02
      if(spt.eq."B7.5") bcv = -0.875
      if(spt.eq."B8  ") bcv = -0.73
      if(spt.eq."B8.5") bcv = -0.595
      if(spt.eq."B9  ") bcv = -0.46
      if(spt.eq."B9.5") bcv = -0.38
      if(spt.eq."A0  ") bcv = -0.17
      if(spt.eq."A0.5") bcv = -0.14
      if(spt.eq."A1  ") bcv = -0.11
      if(spt.eq."A1.5") bcv = -0.08
      if(spt.eq."A2  ") bcv = -0.05
      if(spt.eq."A2.5") bcv = -0.035
      if(spt.eq."A3  ") bcv = -0.02
      if(spt.eq."A3.5") bcv = -0.01
      if(spt.eq."A4  ") bcv =  0.00
      if(spt.eq."A4.5") bcv =  0.005
      if(spt.eq."A5  ") bcv =  0.01
      if(spt.eq."A5.5") bcv =  0.015
      if(spt.eq."A6  ") bcv =  0.02
      if(spt.eq."A6.5") bcv =  0.03
      if(spt.eq."A7  ") bcv =  0.04
      if(spt.eq."A7.5") bcv =  0.04
      if(spt.eq."A8  ") bcv =  0.04
      if(spt.eq."A8.5") bcv =  0.04
      if(spt.eq."A9  ") bcv =  0.04
      if(spt.eq."A9.5") bcv =  0.035
      if(spt.eq."F0  ") bcv =  0.03
      if(spt.eq."F0.5") bcv =  0.03
      if(spt.eq."F1  ") bcv =  0.03
      if(spt.eq."F1.5") bcv =  0.025
      if(spt.eq."F2  ") bcv =  0.02
      if(spt.eq."F2.5") bcv =  0.015
      if(spt.eq."F3  ") bcv =  0.01
      if(spt.eq."F3.5") bcv =  0.005
      if(spt.eq."F4  ") bcv =  0.00
      if(spt.eq."F4.5") bcv =  0.00
      if(spt.eq."F5  ") bcv =  0.00
      if(spt.eq."F5.5") bcv = -0.005
      if(spt.eq."F6  ") bcv = -0.01
      if(spt.eq."F6.5") bcv = -0.015
      if(spt.eq."F7  ") bcv = -0.02
      if(spt.eq."F7.5") bcv = -0.025
      if(spt.eq."F8  ") bcv = -0.03
      if(spt.eq."F8.5") bcv = -0.035
      if(spt.eq."F9  ") bcv = -0.04
      if(spt.eq."F9.5") bcv = -0.045
      if(spt.eq."G0  ") bcv = -0.05
      if(spt.eq."G0.5") bcv = -0.055
      if(spt.eq."G1  ") bcv = -0.06
      if(spt.eq."G1.5") bcv = -0.065
      if(spt.eq."G2  ") bcv = -0.07
      if(spt.eq."G2.5") bcv = -0.075
      if(spt.eq."G3  ") bcv = -0.08
      if(spt.eq."G3.5") bcv = -0.09
      if(spt.eq."G4  ") bcv = -0.10
      if(spt.eq."G4.5") bcv = -0.10
      if(spt.eq."G5  ") bcv = -0.10
      if(spt.eq."G5.5") bcv = -0.105
      if(spt.eq."G6  ") bcv = -0.11
      if(spt.eq."G6.5") bcv = -0.115
      if(spt.eq."G7  ") bcv = -0.12
      if(spt.eq."G7.5") bcv = -0.125
      if(spt.eq."G8  ") bcv = -0.13
      if(spt.eq."G8.5") bcv = -0.15
      if(spt.eq."G9  ") bcv = -0.17
      if(spt.eq."G9.5") bcv = -0.18
      if(spt.eq."K0  ") bcv = -0.19
      if(spt.eq."K0.5") bcv = -0.205
      if(spt.eq."K1  ") bcv = -0.22
      if(spt.eq."K1.5") bcv = -0.245
      if(spt.eq."K2  ") bcv = -0.27
      if(spt.eq."K2.5") bcv = -0.325
      if(spt.eq."K3  ") bcv = -0.38
      if(spt.eq."K3.5") bcv = -0.45
      if(spt.eq."K4  ") bcv = -0.52
      if(spt.eq."K4.5") bcv = -0.555
      if(spt.eq."K5  ") bcv = -0.59
      if(spt.eq."K5.5") bcv = -0.66
      if(spt.eq."K6  ") bcv = -0.73
      if(spt.eq."K6.5") bcv = -0.865
      if(spt.eq."K7  ") bcv = -1.00
      if(spt.eq."K7.5") bcv = -1.015
      if(spt.eq."K8  ") bcv = -1.03
      if(spt.eq."K8.5") bcv = -1.085
      if(spt.eq."K9  ") bcv = -1.14
      if(spt.eq."K9.5") bcv = -1.165
      if(spt.eq."M0  ") bcv = -1.19
      if(spt.eq."M0.5") bcv = -1.30
      if(spt.eq."M1  ") bcv = -1.41
      if(spt.eq."M1.5") bcv = -1.47
      if(spt.eq."M2  ") bcv = -1.53
      if(spt.eq."M2.5") bcv = -1.76
      if(spt.eq."M3  ") bcv = -2.99
      if(spt.eq."M3.5") bcv = -2.16
      if(spt.eq."M4  ") bcv = -2.33
      if(spt.eq."M4.5") bcv = -2.62
      if(spt.eq."M5  ") bcv = -2.91
      if(spt.eq."M5.5") bcv = -3.20
      if(spt.eq."M6  ") bcv = -3.49
      if(spt.eq."M6.5") bcv = -3.88
      if(spt.eq."M7  ") bcv = -4.27
      if(spt.eq."M7.5") bcv = -4.39
      if(spt.eq."M8  ") bcv = -4.51
      if(spt.eq."M8.5") bcv = -4.745
      if(spt.eq."M9  ") bcv = -4.98

      if(spt.eq."O9  ") bvo = -0.311  
      if(spt.eq."O9.5") bvo = -0.305  
      if(spt.eq."B0  ") bvo = -0.300   
      if(spt.eq."B0.5") bvo = -0.295  
      if(spt.eq."B1  ") bvo = -0.278   
      if(spt.eq."B1.5") bvo = -0.250  
      if(spt.eq."B2  ") bvo = -0.210   
      if(spt.eq."B2.5") bvo = -0.193  
      if(spt.eq."B3  ") bvo = -0.178   
      if(spt.eq."B3.5") bvo = -0.172  
      if(spt.eq."B4  ") bvo = -0.165   
      if(spt.eq."B4.5") bvo = -0.161  
      if(spt.eq."B5  ") bvo = -0.156   
      if(spt.eq."B5.5") bvo = -0.148  
      if(spt.eq."B6  ") bvo = -0.140   
      if(spt.eq."B6.5") bvo = -0.134  
      if(spt.eq."B7  ") bvo = -0.128   
      if(spt.eq."B7.5") bvo = -0.119  
      if(spt.eq."B8  ") bvo = -0.109   
      if(spt.eq."B8.5") bvo = -0.090  
      if(spt.eq."B9  ") bvo = -0.070   
      if(spt.eq."B9.5") bvo = -0.050  
      if(spt.eq."A0  ") bvo =  0.000   
      if(spt.eq."A0.5") bvo =  0.022  
      if(spt.eq."A1  ") bvo =  0.043   
      if(spt.eq."A1.5") bvo =  0.059  
      if(spt.eq."A2  ") bvo =  0.074   
      if(spt.eq."A2.5") bvo =  0.082  
      if(spt.eq."A3  ") bvo =  0.090   
      if(spt.eq."A3.5") bvo =  0.115  
      if(spt.eq."A4  ") bvo =  0.140   
      if(spt.eq."A4.5") bvo =  0.150  
      if(spt.eq."A5  ") bvo =  0.160   
      if(spt.eq."A5.5") bvo =  0.165  
      if(spt.eq."A6  ") bvo =  0.170   
      if(spt.eq."A6.5") bvo =  0.190  
      if(spt.eq."A7  ") bvo =  0.210   
      if(spt.eq."A7.5") bvo =  0.232  
      if(spt.eq."A8  ") bvo =  0.253   
      if(spt.eq."A8.5") bvo =  0.254  
      if(spt.eq."A9  ") bvo =  0.255   
      if(spt.eq."A9.5") bvo =  0.275  
      if(spt.eq."F0  ") bvo =  0.294  
      if(spt.eq."F0.5") bvo =  0.314  
      if(spt.eq."F1  ") bvo =  0.334  
      if(spt.eq."F1.5") bvo =  0.354  
      if(spt.eq."F2  ") bvo =  0.374  
      if(spt.eq."F2.5") bvo =  0.382  
      if(spt.eq."F3  ") bvo =  0.389  
      if(spt.eq."F3.5") bvo =  0.401  
      if(spt.eq."F4  ") bvo =  0.412  
      if(spt.eq."F4.5") bvo =  0.425  
      if(spt.eq."F5  ") bvo =  0.438  
      if(spt.eq."F5.5") bvo =  0.461  
      if(spt.eq."F6  ") bvo =  0.484  
      if(spt.eq."F6.5") bvo =  0.497  
      if(spt.eq."F7  ") bvo =  0.510  
      if(spt.eq."F7.5") bvo =  0.520  
      if(spt.eq."F8  ") bvo =  0.530  
      if(spt.eq."F8.5") bvo =  0.541  
      if(spt.eq."F9  ") bvo =  0.552  
      if(spt.eq."F9.5") bvo =  0.570  
      if(spt.eq."G0  ") bvo =  0.588  
      if(spt.eq."G0.5") bvo =  0.596  
      if(spt.eq."G1  ") bvo =  0.604  
      if(spt.eq."G1.5") bvo =  0.623  
      if(spt.eq."G2  ") bvo =  0.642  
      if(spt.eq."G2.5") bvo =  0.652  
      if(spt.eq."G3  ") bvo =  0.661  
      if(spt.eq."G3.5") bvo =  0.668  
      if(spt.eq."G4  ") bvo =  0.674  
      if(spt.eq."G4.5") bvo =  0.677  
      if(spt.eq."G5  ") bvo =  0.680  
      if(spt.eq."G5.5") bvo =  0.692  
      if(spt.eq."G6  ") bvo =  0.704  
      if(spt.eq."G6.5") bvo =  0.709  
      if(spt.eq."G7  ") bvo =  0.713  
      if(spt.eq."G7.5") bvo =  0.725  
      if(spt.eq."G8  ") bvo =  0.737  
      if(spt.eq."G8.5") bvo =  0.757  
      if(spt.eq."G9  ") bvo =  0.777  
      if(spt.eq."G9.5") bvo =  0.797  
      if(spt.eq."K0  ") bvo =  0.816  
      if(spt.eq."K0.5") bvo =  0.832  
      if(spt.eq."K1  ") bvo =  0.847  
      if(spt.eq."K1.5") bvo =  0.870  
      if(spt.eq."K2  ") bvo =  0.893  
      if(spt.eq."K2.5") bvo =  0.942  
      if(spt.eq."K3  ") bvo =  0.990  
      if(spt.eq."K3.5") bvo =  1.045  
      if(spt.eq."K4  ") bvo =  1.100  
      if(spt.eq."K4.5") bvo =  1.117  
      if(spt.eq."K5  ") bvo =  1.134  
      if(spt.eq."K5.5") bvo =  1.196  
      if(spt.eq."K6  ") bvo =  1.257  
      if(spt.eq."K6.5") bvo =  1.297  
      if(spt.eq."K7  ") bvo =  1.336  
      if(spt.eq."K7.5") bvo =  1.359  
      if(spt.eq."K8  ") bvo =  1.382  
      if(spt.eq."K8.5") bvo =  1.400  
      if(spt.eq."K9  ") bvo =  1.418  
      if(spt.eq."K9.5") bvo =  1.425  
      if(spt.eq."M0  ") bvo =  1.431   
      if(spt.eq."M0.5") bvo =  1.458 
      if(spt.eq."M1  ") bvo =  1.484 
      if(spt.eq."M1.5") bvo =  1.494 
      if(spt.eq."M2  ") bvo =  1.503 
      if(spt.eq."M2.5") bvo =  1.528 
      if(spt.eq."M3  ") bvo =  1.553 
      if(spt.eq."M3.5") bvo =  1.606 
      if(spt.eq."M4  ") bvo =  1.658 
      if(spt.eq."M4.5") bvo =  1.766 
      if(spt.eq."M5  ") bvo =  1.874 
      if(spt.eq."M5.5") bvo =  1.953 
      if(spt.eq."M6  ") bvo =  2.031 
      if(spt.eq."M6.5") bvo =  2.075 
      if(spt.eq."M7  ") bvo =  2.119 
      if(spt.eq."M7.5") bvo =  2.158 
      if(spt.eq."M8  ") bvo =  2.197
      if(spt.eq."M8.5") bvo =  2.199
      if(spt.eq."M9  ") bvo =  2.200 

* Teff values
      if (spt.eq."O9  ") teff = 33000
      if (spt.eq."O9.5") teff = 31900
      if (spt.eq."B0  ") teff = 31500
      if (spt.eq."B0.5") teff = 28500
      if (spt.eq."B1  ") teff = 26000
      if (spt.eq."B1.5") teff = 24000
      if (spt.eq."B2  ") teff = 20500
      if (spt.eq."B2.5") teff = 18300
      if (spt.eq."B3  ") teff = 17300
      if (spt.eq."B3.5") teff = 17000
      if (spt.eq."B4  ") teff = 16700
      if (spt.eq."B4.5") teff = 16200
      if (spt.eq."B5  ") teff = 15700
      if (spt.eq."B5.5") teff = 15100
      if (spt.eq."B6  ") teff = 14500
      if (spt.eq."B6.5") teff = 14200
      if (spt.eq."B7  ") teff = 13900
      if (spt.eq."B7.5") teff = 13150
      if (spt.eq."B8  ") teff = 12400
      if (spt.eq."B8.5") teff = 11700
      if (spt.eq."B9  ") teff = 11000
      if (spt.eq."B9.5") teff = 10600
      if (spt.eq."A0  ") teff =  9550
      if (spt.eq."A0.5") teff =  9375
      if (spt.eq."A1  ") teff =  9200
      if (spt.eq."A1.5") teff =  8980
      if (spt.eq."A2  ") teff =  8760
      if (spt.eq."A2.5") teff =  8655
      if (spt.eq."A3  ") teff =  8550
      if (spt.eq."A3.5") teff =  8410
      if (spt.eq."A4  ") teff =  8270
      if (spt.eq."A4.5") teff =  8175
      if (spt.eq."A5  ") teff =  8080
      if (spt.eq."A5.5") teff =  8040
      if (spt.eq."A6  ") teff =  8000
      if (spt.eq."A6.5") teff =  7900
      if (spt.eq."A7  ") teff =  7800
      if (spt.eq."A7.5") teff =  7650
      if (spt.eq."A8  ") teff =  7500
      if (spt.eq."A8.5") teff =  7470
      if (spt.eq."A9  ") teff =  7440
      if (spt.eq."A9.5") teff =  7320
      if (spt.eq."F0  ") teff =  7200
      if (spt.eq."F0.5") teff =  7115
      if (spt.eq."F1  ") teff =  7030
      if (spt.eq."F1.5") teff =  6920
      if (spt.eq."F2  ") teff =  6810
      if (spt.eq."F2.5") teff =  6765
      if (spt.eq."F3  ") teff =  6720
      if (spt.eq."F3.5") teff =  6680
      if (spt.eq."F4  ") teff =  6640
      if (spt.eq."F4.5") teff =  6575
      if (spt.eq."F5  ") teff =  6510
      if (spt.eq."F5.5") teff =  6425
      if (spt.eq."F6  ") teff =  6340
      if (spt.eq."F6.5") teff =  6290
      if (spt.eq."F7  ") teff =  6240
      if (spt.eq."F7.5") teff =  6195
      if (spt.eq."F8  ") teff =  6150
      if (spt.eq."F8.5") teff =  6090
      if (spt.eq."F9  ") teff =  6040
      if (spt.eq."F9.5") teff =  5990
      if (spt.eq."G0  ") teff =  5940
      if (spt.eq."G0.5") teff =  5910
      if (spt.eq."G1  ") teff =  5880
      if (spt.eq."G1.5") teff =  5830
      if (spt.eq."G2  ") teff =  5780
      if (spt.eq."G2.5") teff =  5740
      if (spt.eq."G3  ") teff =  5700
      if (spt.eq."G3.5") teff =  5670
      if (spt.eq."G4  ") teff =  5640
      if (spt.eq."G4.5") teff =  5630
      if (spt.eq."G5  ") teff =  5620
      if (spt.eq."G5.5") teff =  5600
      if (spt.eq."G6  ") teff =  5580
      if (spt.eq."G6.5") teff =  5550
      if (spt.eq."G7  ") teff =  5520
      if (spt.eq."G7.5") teff =  5505
      if (spt.eq."G8  ") teff =  5490
      if (spt.eq."G8.5") teff =  5415
      if (spt.eq."G9  ") teff =  5340
      if (spt.eq."G9.5") teff =  5300
      if (spt.eq."K0  ") teff =  5260
      if (spt.eq."K0.5") teff =  5215
      if (spt.eq."K1  ") teff =  5170
      if (spt.eq."K1.5") teff =  5105
      if (spt.eq."K2  ") teff =  5040
      if (spt.eq."K2.5") teff =  4925
      if (spt.eq."K3  ") teff =  4810
      if (spt.eq."K3.5") teff =  4705
      if (spt.eq."K4  ") teff =  4600
      if (spt.eq."K4.5") teff =  4550
      if (spt.eq."K5  ") teff =  4500
      if (spt.eq."K5.5") teff =  4400
      if (spt.eq."K6  ") teff =  4300
      if (spt.eq."K6.5") teff =  4150
      if (spt.eq."K7  ") teff =  4000
      if (spt.eq."K7.5") teff =  3980
      if (spt.eq."K8  ") teff =  3960
      if (spt.eq."K8.5") teff =  3910
      if (spt.eq."K9  ") teff =  3860
      if (spt.eq."K9.5") teff =  3835
      if (spt.eq."M0  ") teff =  3810
      if (spt.eq."M0.5") teff =  3725
      if (spt.eq."M1  ") teff =  3640
      if (spt.eq."M1.5") teff =  3570 
      if (spt.eq."M2  ") teff =  3500
      if (spt.eq."M2.5") teff =  3390
      if (spt.eq."M3  ") teff =  3280
      if (spt.eq."M3.5") teff =  3180
      if (spt.eq."M4  ") teff =  3080
      if (spt.eq."M4.5") teff =  2975
      if (spt.eq."M5  ") teff =  2870
      if (spt.eq."M5.5") teff =  2785
      if (spt.eq."M6  ") teff =  2700
      if (spt.eq."M6.5") teff =  2600
      if (spt.eq."M7  ") teff =  2500
      if (spt.eq."M7.5") teff =  2470
      if (spt.eq."M8  ") teff =  2440
      if (spt.eq."M8.5") teff =  2380
      if (spt.eq."M9  ") teff =  2320
      if (spt.eq."M9.5") teff =  2285
      if (spt.eq."L0  ") teff =  2250
      if (spt.eq."L0.5") teff =  2225
      if (spt.eq."L1  ") teff =  2200
      if (spt.eq."L1.5") teff =  2130
      if (spt.eq."L2  ") teff =  2060
      if (spt.eq."L2.5") teff =  1990
      if (spt.eq."L3  ") teff =  1920
      if (spt.eq."L3.5") teff =  1855
      if (spt.eq."L4  ") teff =  1790
      if (spt.eq."L4.5") teff =  1770
      if (spt.eq."L5  ") teff =  1750
      if (spt.eq."L5.5") teff =  1713
      if (spt.eq."L6  ") teff =  1675
      if (spt.eq."L6.5") teff =  1638
      if (spt.eq."L7  ") teff =  1600
      if (spt.eq."L7.5") teff =  1500
      if (spt.eq."L8  ") teff =  1400
      if (spt.eq."L8.5") teff =  1388
      if (spt.eq."L9  ") teff =  1375
      if (spt.eq."L9.5") teff =  1363
      if (spt.eq."T0  ") teff =  1350
      if (spt.eq."T0.5") teff =  1338
      if (spt.eq."T1  ") teff =  1325
      if (spt.eq."T1.5") teff =  1313
      if (spt.eq."T2  ") teff =  1300
      if (spt.eq."T2.5") teff =  1263
      if (spt.eq."T3  ") teff =  1225
      if (spt.eq."T3.5") teff =  1188
      if (spt.eq."T4  ") teff =  1150 
      if (spt.eq."T4.5") teff =  1113
      if (spt.eq."T5  ") teff =  1075
      if (spt.eq."T5.5") teff =  1038
      if (spt.eq."T6  ") teff =  1000
      if (spt.eq."T6.5") teff =   975
      if (spt.eq."T7  ") teff =   950
      if (spt.eq."T7.5") teff =   850
      if (spt.eq."T8  ") teff =   750
      if (spt.eq."T8.5") teff =   700
      if (spt.eq."T9  ") teff =   650
      if (spt.eq."T9.5") teff =   500
      if (spt.eq."Y0  ") teff =   350

* Error message if invalid spectral type
      if (bvo.lt.-1d99.or.teff.lt.-1d99.or.bcv.lt.-1d99) then
         write(*,*)"ERROR - invalid spectral type"
      endif
      return
      end

***************************************************
* ROSSBY_RHK_MAMAJEK estimates Rossby Number      *
* from the chromospheric activity logR'HK.        *
* Fit made to solar-type dwarfs (0.5 < B-V < 0.9, *
* |delta(Mv)| < 1 mag) with periods and logR'HK.  *
* Published in Mamajek & Hillenbrand (2008).      *
* INPUT                                           *
* lgrphk  = logR'HK [dex]                         *
* ulgrphk = uncertainty in logR'HK [dex]          *
* OUTPUT                                          *
* ro = Rossby number predicted from logR'HK       *
* uro = uncertainty in Ro from calibration        *
* rms = rms scatter in Ro from fit                *
* unc = total uncertainty in predicted Ro value   *
*                                                 *
* E Mamajek 7/10/2008                             *
***************************************************
      subroutine rossby_rhk_mamajek(lgrphk,ulgrphk,ro,uro,rms,unc)
      implicit double precision(a-z)
      if (lgrphk .le. -4.3543d0) then ! ACTIVE REGIME
       b   =  0.8084d0
       ub  =  0.0144d0
       m   = -2.9657d0
       um  =  0.0984d0
       c   =  4.520d0
       uc  =  0.000d0 ! assume no error in constant
       rms =  0.158d0 ! from fit
      else !                    VERY ACTIVE REGIME
       b   =  0.2325d0
       ub  =  0.0153d0
       m   = -0.6894d0
       um  =  0.0631d0
       c   =  4.232d0
       uc  =  0.000d0 ! assume no error in constant
       rms =  0.10d0  ! from fit
      endif
      x  = lgrphk + c
      ux = dsqrt(ulgrphk*ulgrphk + uc*uc)
      ro = b + m*x 
      uro = dsqrt(ub*ub + x*x*um*um + m*m*ux*ux)
      unc = dsqrt(uro**2 + rms**2)
      return 
      end 

**********************************************************
* GYRO_AGE_MEIBOM calculates the gryochronology age via  *
* new gyrochronology parameters from Meibom+2009 (ApJ    *
* 695, 679). The formalism is the same as that of        *
* Barnes (2007).                                         *
* INPUT                                                  *
* bv  = intrinsic B-V color [mag]                        *
* per = rotation period [days]                           *
* OUTPUT                                                 *
* logt = gyro age [yr]                                   *
* E Mamajek 2/7/2011                                     *
**********************************************************
      subroutine gyro_age_meibom(bv,per,logt,ulogt,tmyr,utmyr)
      implicit double precision (a-z)
* apologies for nomenclature here -- "bv" is the intrinsic stellar
* color, whereas "bvo" is the color singularity in the gyro equations.-EEM
* parameters from Meibom et al. 2009
      n   = 0.52d0  ! see equation 3 of Meibom09
      a   = 0.770d0 ! "a" in Mamajek08 is called "c" in Meibom09 eqn. 3
      b   = 0.553d0 ! "b" in Mamajek08 is called "f" in Meibom09 eqn. 3
      bvo = 0.472d0 ! "c" in Mamajek08 is called "d" in Meibom09 eqn. 3
* equation 3 from Barnes (2007)
      logp = dlog10(per)
      loga = dlog10(a)
      x    = bv - bvo
      logx = dlog10(x)
      logt = 6.0d0 + (logp - loga - b*logx)/n
      tmyr = 10.0d0**((logp - loga - b*logx)/n)
* errors
      c = 4.5d0
      d = 0.5*dlog(tmyr)**2.0d0
      e = 2.0d0*per**0.6d0
      f = (0.6d0/x)**2.0d0
      frac = 0.02d0 * dsqrt(c + d + e + f)
      utmyr = frac*tmyr
      logt1 = dlog10(6.0d0 + tmyr + utmyr)
      logt2 = dlog10(6.0d0 + tmyr - utmyr)
      ulogt = dabs(logt1 - logt2)/2.0d0
      if (bv .le. bvo .or. bv .ge. 1.4) then
         write(*,*)" WARNING: B-V color out of range for gyro age"
      endif 
      end

*******************
* ROSSBY_XRAY estimates logRx = logLx/Lbol from Rossby number (using
* Noyes+1984 values as function of B-V) using a fit to an
* X-ray-unbiased sample of solar-type stars -- the 37 Mt. Wilson
* program stars with well-measured rotation periods from
* Donahue+1996. *All* of the Donahue+1996 stars have been detected
* in X-rays by some survey, and their Lx and Rx values have been
* calculated by EEM following Fleming+ 1995. EEM finds a fit
* in the unsaturated regime of:
* log Rx = -5.066(+-0.042) - 2.617(+-0.182)*(log Ro + 0.02861)
* rms(logRx) = +-0.25 dex.
* I assume that log Rx(saturated) = -3.16 following Wright+2011.
* E Mamajek 4/18/2011
**********************
      subroutine rossby_xray(logro,logrx)
      implicit double precision (a-z)
      a = -5.066d0
      b = -2.617d0
      c = -0.02861d0
      if (logro .le. -0.7569d0) then ! saturated regime
        logrx = -3.16d0
      else  ! unsaturated regime
        logrx = a + b*(logro - c)
      endif
      if (logro .gt. 0.44) then
        write(*,*)" WARNING: logRo > 0.44 ; beyond Ro-Rx fit"
      endif
      continue
      end 

***************************************************************
* XRAY_ROSSBY estimates logRo from logRx = logLx/Lbol (using
* Noyes+1984 values as function of B-V) using a fit to an
* X-ray-unbiased sample of solar-type stars -- the 37 Mt. Wilson
* program stars with well-measured rotation periods from
* Donahue+1996. *All* of the Donahue+1996 stars have been detected
* in X-rays by some survey, and their Lx and Rx values have been
* calculated by EEM following Fleming+ 1995. EEM finds a fit
* in the unsaturated regime of:
* log Ro = (-0.029+-0.016) - 0.382(+-0.025)(log Rx + 5.0655)
* I assume that log Rx(saturated) = -3.16 following Wright+2011.
* E Mamajek 4/18/2011
**********************
      subroutine xray_rossby(logrx,logro)
      implicit double precision (a-z)
      a = -0.029d0
      b = -0.3821d0
      c = -5.0655d0
      if (logrx .ge. -3.16d0) then ! saturated regime
        logro = -0.7569d0 
        write(*,*)' WARNING: logRx > -3.16 ; SATURATED'
        write(*,*)' WARNING: Rotation period & age are upper limits'
      else  ! unsaturated regime
        logro = a + b*(logrx - c)
      endif
      if (logro .gt. 0.44) then
        write(*,*)" WARNING: logRo > 0.44 ; beyond Ro-Rx fit"
      endif
      continue
      end 

**************************************
* OLSON75_R_BVO_EBV calculates R value from B-Vo and E(B-V)
* and uncertainty
* E Mamajek 10/3/2011
**************************************
      subroutine Olson75_R_BVO_EBV(bvo,ebv,r,av)
      implicit double precision(a-z)
      a   = 3.25d0
      b   = 0.25d0
      c   = 0.05d0
      r   = a + b*bvo + c*ebv
      av  = r*ebv
      return
      end


**************************************
* OLSON75_AV_BV_IterSolv_BVO_EBV 
* given Av and observed B-V, iteratively solve
* for (B-V)o, E(B-V), R.
* E Mamajek 11/12/2011
**************************************
      subroutine Olson75_AV_BV_IterSolv_BVO_EBV(av,bv,bvo,ebv,r)
      implicit double precision(a-z)
      a    = 3.25d0
      b    = 0.25d0
      c    = 0.05d0
      diff = 9.99d99  
      ebv  = av/3.25  ! initial estimate of E(B-V)
      bvo  = bv - ebv ! initial estimate of (B-V)o
      do while (diff .ge. 1.0d-5)
        bvoold = bvo 
        r      = a + b*bvo + c*ebv
        ebv    = av/r
        bvo    = bv - ebv 
        diff   = dabs(bvoold - bvo)
      enddo
      return
      end

**************************************
* OLSON75_AV_BV0_IterSolv_BV_EBV 
* given Av and (B-V)o, iteratively solve
* for (B-V), E(B-V), and R. 
* E Mamajek 11/12/2011
**************************************
      subroutine Olson75_AV_BV0_IterSolv_BV_EBV(av,bvo,bv,ebv,r)
      implicit double precision(a-z)
      a    = 3.25d0
      b    = 0.25d0
      c    = 0.05d0
      diff = 9.99d99  
      ebv  = av/3.25   ! initial estimate of E(B-V)
      bv   = bvo + ebv ! initial estimate of (B-V)
      do while (diff .ge. 1.0d-5)
        bvold  = bv  
        r      = a + b*bvo + c*ebv
        ebv    = av/r
        bv     = bvo + ebv 
        diff   = dabs(bvold - bv)
*        write(*,288)'    B-V = ',bv
*        write(*,288)' E(B-V) = ',ebv
      enddo
 288  format(a25,f7.3,a)
      return
      end

*****************************************
* V-K => B-V   
* Blue fit to HIP d<75pc(S/N>8) stars with JHKerror<0.05mag(AAA)
* [-0.08 > V-K > 5.63 ; -0.05 < B-V < 1.69] rms_B-V = 0.041
* Red fit to Reid CNS d<25pc stars [5.0 < V-K < 7.5] rms=0.077mag
* welded together fits at V-K = 5.1.
* E Mamajek 12/23/2010
* 10/28/2011: updated B-star fit with line
*****************************************
      subroutine vk_bv(vk,bv,bk)
      implicit double precision (a-z)
      vksplit = 5.1d0 ! where to weld fits...
      if (vk .lt. -0.07) then
        bv = -3.498615113853e-02 + 2.910310027313d-01*vk
        bk = vk + bv  
      elseif (vk .ge. -0.07 .and. vk .le. vksplit) then
        bv = -6.258354d-3 + 4.290261d-1*vk - 3.505073d-2*vk**2.0d0 +
     +  3.241736d-2*vk**3.0d0 - 1.062028d-2*vk**4.0d0 + 
     +  9.332255d-4*vk**5.0d0
        bk = vk + bv
      elseif (vk .ge. vksplit .and. vk .le. 8.5d0) then ! Reid/CNS3 fit
        a0 = -0.9001797d0
        a1 =  0.6831021d0
        a2 = -0.03775829d0
        bv = a0 + a1*vk + a2*vk*vk ! rms = 0.077 mag
        bk = vk + bv 
      else
        bv = -0.0d0
        bk = -0.0d0
      endif
      return
      end
