****************************************************** * 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