consistent use of constants, now defined once in mkconstants.pro

jdeast · Jul 14, 2017 · 935bcb6 · 935bcb6
1 parent 28d3b7b
commit 935bcb6
Show file tree

Hide file tree

Showing 18 changed files with 417 additions and 1,296 deletions.
diff --git a/constants.pro b/constants.pro
diff --git a/derivepars.pro b/derivepars.pro
@@ -1,25 +1,19 @@
 pro derivepars, ss
 
-;; constants
-G = 2942.71377d0 ;; R_sun^3/(m_sun*day^2), Torres 2010
-AU = 215.094177d0 ;; R_sun
-mjup = 0.000954638698d0 ;; m_sun
-rjup = 0.102792236d0    ;; r_sun
-mearth = 0.00000300245 ;; m_sun
-rearth = 0.0091705248 ;; r_sun
-
-sigmab = 5.670373d-5/3.839d33*6.9566d10^2 ;; Stefan-boltzmann Constant (L_sun/(r_sun^2*K^4))
-sigmasb = 5.670373d-5
-
 ss.star.mstar.value = 10^ss.star.logmstar.value
-ss.star.rhostar.value = ss.star.mstar.value/(ss.star.rstar.value^3)*1.41135837d0
-ss.star.logg.value = alog10(G*ss.star.mstar.value/(ss.star.rstar.value^2)*9.31686171d0)
-ss.star.lstar.value = 4d0*!dpi*ss.star.rstar.value^2*ss.star.teff.value^4*sigmaB 
+ss.star.rhostar.value = ss.star.mstar.value/(ss.star.rstar.value^3)*ss.constants.rhosun ;; rho_sun
+ss.star.logg.value = alog10(ss.star.mstar.value/(ss.star.rstar.value^2)*ss.constants.gravitysun) ;; cgs
+ss.star.lstar.value = 4d0*!dpi*ss.star.rstar.value^2*ss.star.teff.value^4*ss.constants.sigmab/ss.constants.lsun*ss.constants.rsun^2 ;; lsun
 
 ;; derive the absolute magnitude and distance
 nsteps = n_elements(ss.star.teff.value)
 ss.star.parallax.value = 1d3/ss.star.distance.value ;; mas
 
+for j=0, ss.ntel-1 do begin
+   positive = where(ss.telescope[j].jittervar.value gt 0d0)
+   ss.telescope[j].jitter.value[positive] = sqrt(ss.telescope[j].jittervar.value[positive])
+endfor
+
 for i=0, ss.nplanets-1 do begin
 
    ;; eccentricity and argument of periastron
@@ -33,6 +27,7 @@ for i=0, ss.nplanets-1 do begin
    zero = where(ss.planet[i].e.value eq 0d0,complement=nonzero)
    if zero[0] ne -1 then ss.planet[i].omega.value[zero] = !dpi/2d0 
 
+   ss.planet[i].lambdadeg.value = ss.planet[i].lambda.value*180d0/!dpi
    ss.planet[i].omegadeg.value = ss.planet[i].omega.value*180d0/!dpi
    ss.planet[i].esinw.value = ss.planet[i].e.value*sin(ss.planet[i].omega.value)
    ss.planet[i].ecosw.value = ss.planet[i].e.value*cos(ss.planet[i].omega.value)
@@ -59,20 +54,20 @@ for i=0, ss.nplanets-1 do begin
                                                   ss.planet[i].e.value[positive],$
                                                   ss.planet[i].i.value[positive],$
                                                   ss.planet[i].period.value[positive],$
-                                                  ss.star.mstar.value[positive])
+                                                  ss.star.mstar.value[positive], G=ss.constants.G/1d6) ;; convert G to SI
 
-   ss.planet[i].mp.value = ss.planet[i].mpsun.value/mjup               ;; m_jupiter
-   ss.planet[i].msini.value = ss.planet[i].mp.value*sini               ;; m_jupiter
-   ss.planet[i].mpearth.value = ss.planet[i].mpsun.value/mearth        ;; m_earth
-   ss.planet[i].msiniearth.value = ss.planet[i].mpearth.value*sini     ;; m_earth
-   ss.planet[i].q.value = ss.planet[i].mpsun.value/ss.star.mstar.value ;; unitless
+   ss.planet[i].mp.value = ss.planet[i].mpsun.value/ss.constants.GMJupiter*ss.constants.GMSun    ;; m_jupiter
+   ss.planet[i].msini.value = ss.planet[i].mp.value*sini                                         ;; m_jupiter
+   ss.planet[i].mpearth.value = ss.planet[i].mpsun.value/ss.constants.GMearth*ss.constants.GMSun ;; m_earth
+   ss.planet[i].msiniearth.value = ss.planet[i].mpearth.value*sini                               ;; m_earth
+   ss.planet[i].q.value = ss.planet[i].mpsun.value/ss.star.mstar.value                           ;; unitless
 
    ss.planet[i].arsun.value=(G*(ss.star.mstar.value+ss.planet[i].mpsun.value)*ss.planet[i].period.value^2/(4d0*!dpi^2))^(1d0/3d0) ;; (a1 + a2)/rsun
    ss.planet[i].ar.value = ss.planet[i].arsun.value/ss.star.rstar.value                                                        ;; (a1 + a2)/rstar
-   ss.planet[i].a.value = ss.planet[i].arsun.value/AU                                                                          ;; AU
+   ss.planet[i].a.value = ss.planet[i].arsun.value/ss.constants.au*ss.constants.Rsun                                           ;; AU
    ss.planet[i].rpsun.value = ss.planet[i].p.value*ss.star.rstar.value                                                         ;; r_sun
-   ss.planet[i].rp.value = ss.planet[i].rpsun.value/rjup                                                                       ;; r_jupiter
-   ss.planet[i].rpearth.value = ss.planet[i].rpsun.value/rearth                                                                ;; r_earth
+   ss.planet[i].rp.value = ss.planet[i].rpsun.value/ss.constants.RJupiter                                                      ;; r_jupiter
+   ss.planet[i].rpearth.value = ss.planet[i].rpsun.value/ss.constants.REarth                                                   ;; r_earth
 
    ;; time of periastron
    ss.planet[i].phase.value=exofast_getphase(ss.planet[i].e.value,ss.planet[i].omega.value,/pri)  
@@ -150,9 +145,9 @@ for i=0, ss.nplanets-1 do begin
    ss.planet[i].psg.value = (ss.star.rstar.value+ss.planet[i].rpsun.value)/ss.planet[i].arsun.value*(1d0 - ss.planet[i].esinw.value)/(1d0-ss.planet[i].e.value^2) ;; eq 9, Winn 2010
    ss.planet[i].ps.value = (ss.star.rstar.value-ss.planet[i].rpsun.value)/ss.planet[i].arsun.value*(1d0 - ss.planet[i].esinw.value)/(1d0-ss.planet[i].e.value^2)  ;; eq 9, Winn 2010
 
-   ss.planet[i].rhop.value = ss.planet[i].mpsun.value/(ss.planet[i].rpsun.value^3)*1.41135837d0
-   ss.planet[i].loggp.value = alog10(G*ss.planet[i].mpsun.value/ss.planet[i].rpsun.value^2*9.31686171d0) ;; cgs
-   ss.planet[i].safronov.value = ss.planet[i].ar.value*ss.planet[i].q.value/ss.planet[i].p.value
+   ss.planet[i].rhop.value = ss.planet[i].mpsun.value/(ss.planet[i].rpsun.value^3)*ss.constants.rhosun ;; cgs
+   ss.planet[i].loggp.value = alog10(ss.planet[i].mpsun.value/ss.planet[i].rpsun.value^2*ss.constants.gravitysun) ;; cgs
+   ss.planet[i].safronov.value = ss.planet[i].ar.value*ss.planet[i].q.value/ss.planet[i].p.value ;; unitless
 
    ;; depth != delta if grazing (ignore limb darkening)
    ss.planet[i].delta.value = ss.planet[i].p.value^2

diff --git a/exofast_chi2v2.pro b/exofast_chi2v2.pro
@@ -169,7 +169,7 @@ COMMON chi2_block, ss
 ;; populate the structure with the new parameters
 if n_elements(pars) ne 0 then pars2str, pars, ss
 
-AU = 215.094177d0
+au = ss.constants.au/ss.constants.rsun
 
 ;; derive all required parameters 
 ;; (this may change depending on parameterization)
@@ -245,11 +245,12 @@ endif
 ss.star.mstar.value = 10^ss.star.logmstar.value
 ;; use the YY tracks to guide the stellar parameters
 if ~ss.noyy then begin
-   if keyword_set(psname) then begin
-      yychi2 = massradius_yy3(ss.star.mstar.value, ss.star.feh.value, ss.star.age.value, ss.star.teff.value,yyrstar=rstar, debug=ss.debug, psname=psname+'.yy.eps') 
-   endif else begin
-      yychi2 = massradius_yy3(ss.star.mstar.value, ss.star.feh.value, ss.star.age.value, ss.star.teff.value,yyrstar=rstar, debug=ss.debug)
-   endelse
+   if keyword_set(psname) then epsname = psname+'.yy.eps'
+   yychi2 = massradius_yy3(ss.star.mstar.value, ss.star.feh.value, $
+                           ss.star.age.value, ss.star.teff.value,$
+                           yyrstar=rstar, debug=ss.debug, psname=epsname, $
+                           sigmab=ss.constants.sigmab/ss.constants.lsun*ss.constants.rsun^2, $
+                           gravitysun=ss.constants.gravitysun) 
    chi2 += yychi2
    ;print, 'YY penalty = ' + strtrim(yychi2,2)
    if ~finite(chi2) then begin
@@ -265,6 +266,12 @@ if ss.star.errscale.value le 0 then begin
    return, !values.d_infinity
 endif
 
+bad = where(ss.doptom.dtscale.value le 0, nbad)
+if nbad gt 0 then begin
+   if ss.debug then print, 'dtscale is bad (' + strtrim(bad,2) + ')'
+   return, !values.d_infinity
+endif
+
 if ss.star.alpha.value lt -0.3d0 or ss.star.alpha.value gt 0.7d0 then begin
    if ss.debug then print, 'alpha is bad (' + strtrim(ss.star.alpha.value,2) + ')'
    return, !values.d_infinity
@@ -279,20 +286,19 @@ if ss.amoeba then begin
       return, !values.d_infinity
    endif
 
-   G = 1.3271244d20 ;; m^3/(M_sun*second^2) Standish 1995
    for j=0, ss.nplanets-1 do begin
       if ss.planet[j].chen then begin
          if ss.planet[j].rpearth.value le 0d0 then begin
             if ss.debug then print, 'rpearth is bad'
             return, !values.d_infinity
          endif         
          if ~ss.planet[j].fitrv and ss.planet[j].fittran then begin
-            mp = massradius_chenreverse(ss.planet[j].rpearth.value)*0.00000300245d0 ;; m_sun
-            k = ((2d0*!dpi*G)/(ss.planet[j].period.value*86400d0*(ss.star[0].mstar.value + mp)^2))^(1d0/3d0)*$
-                mp*sin(acos(ss.planet[j].cosi.value))/sqrt(1d0 - ss.planet[j].e.value^2) ;; m/s
+            mp = massradius_chenreverse(ss.planet[j].rpearth.value)*ss.constants.gmsun/ss.constants.gmearth ;; m_sun
+            k = ((2d0*!dpi*ss.constants.GMsun)/(ss.planet[j].period.value*ss.constants.day*(ss.star[0].mstar.value + mp)^2))^(1d0/3d0)*$
+                mp*sin(acos(ss.planet[j].cosi.value))/sqrt(1d0 - ss.planet[j].e.value^2)/ss.constants.meter ;; m/s
             ss.planet[j].logk.value = alog10(k)
          endif else if ss.planet[j].fitrv and ~ss.planet[j].fittran then begin
-            rp = massradius_chen(ss.planet[j].mpearth.value)*0.0091705248d0 ;; r_sun
+            rp = massradius_chen(ss.planet[j].mpearth.value)*ss.constants.rearth/ss.constants.rsun ;; r_sun
             ss.planet[j].p.value = rp/ss.star.rstar.value
          endif
       endif
@@ -318,6 +324,7 @@ for i=0, n_elements(priors[0,*])-1 do begin
       return, !values.d_infinity
    endif
 
+
    chi2 += ((ss.(priors[0,i])[priors[1,i]].(priors[2,i]).value - $
              ss.(priors[0,i])[priors[1,i]].(priors[2,i]).prior)/$
             ss.(priors[0,i])[priors[1,i]].(priors[2,i]).priorwidth)^2
@@ -327,8 +334,9 @@ for i=0, n_elements(priors[0,*])-1 do begin
 ;          ss.(priors[0,i])[priors[1,i]].(priors[2,i]).priorwidth, $
 ;          ((ss.(priors[0,i])[priors[1,i]].(priors[2,i]).value - $
 ;            ss.(priors[0,i])[priors[1,i]].(priors[2,i]).prior)/$
-;           ss.(priors[0,i])[priors[1,i]].(priors[2,i]).priorwidth)^2
-
+;           ss.(priors[0,i])[priors[1,i]].(priors[2,i]).priorwidth)^2,$
+;          ss.(priors[0,i])[priors[1,i]].(priors[2,i]).lowerbound,$ 
+;          ss.(priors[0,i])[priors[1,i]].(priors[2,i]).upperbound
 
 endfor
 
@@ -441,15 +449,14 @@ endfor
 
 ;; fit the SED
 if file_test(ss.star.fluxfile) then begin
-   if keyword_set(psname) then begin
-      sedchi2 = exofast_sed(ss.star.fluxfile, ss.star.teff.value, ss.star.rstar.value,$
-                            ss.star.av.value, ss.star.distance.value, $
-                            logg=ss.star.logg.value,met=ss.star.feh.value,alpha=ss.star.alpha.value,verbose=ss.debug, f0=f, fp0=fp, ep0=ep, psname=psname+'.sed.eps')
-   endif else begin
-      sedchi2 = exofast_sed(ss.star.fluxfile, ss.star.teff.value, ss.star.rstar.value,$
-                            ss.star.av.value, ss.star.distance.value, $
-                            logg=ss.star.logg.value,met=ss.star.feh.value,alpha=ss.star.alpha.value,verbose=ss.debug, f0=f, fp0=fp, ep0=ep)
-   endelse 
+   if keyword_set(psname) then epsname = psname+'.sed.eps'
+   sedchi2 = exofast_sed(ss.star.fluxfile, ss.star.teff.value, $
+                         ss.star.rstar.value,$
+                         ss.star.av.value, ss.star.distance.value, $
+                         logg=ss.star.logg.value,met=ss.star.feh.value,$
+                         alpha=ss.star.alpha.value,verbose=ss.debug, $
+                         f0=f, fp0=fp, ep0=ep, psname=epsname, $
+                         pc=ss.constants.pc, rsun=ss.constants.rsun)
 
    if ~finite(sedchi2) then begin
       if ss.debug then print, 'sed is bad'
@@ -470,7 +477,7 @@ for j=0, ntelescopes-1 do begin
 
    rv = *(ss.telescope[j].rvptrs)
 
-   if (where(rv.err^2 + ss.telescope[j].jitter.value le 0d0))[0] ne -1 then return, !values.d_infinity
+   if (where(rv.err^2 + ss.telescope[j].jittervar.value le 0d0))[0] ne -1 then return, !values.d_infinity
 
    modelrv = dblarr(n_elements(rv.rv))
    for i=0, nplanets-1 do begin
@@ -495,6 +502,7 @@ for j=0, ntelescopes-1 do begin
                                u1=0d0,t0=t0,deltarv=deltarv)
 
       endif
+
    endfor
    ;; add instrumental offset, slope, and quadratic term
    modelrv += ss.telescope[j].gamma.value + ss.star.slope.value*(rv.bjd-t0) + ss.star.quad.value*(rv.bjd-t0)^2
@@ -506,7 +514,8 @@ for j=0, ntelescopes-1 do begin
    if keyword_set(psname) then $
       forprint, rv.bjd, modelrv, format='(f0.8,x,f0.6)', textout=psname + '.model.' + strtrim(j,2) + '.rv.txt', /nocomment,/silent
 
-   rvchi2 = exofast_like((*ss.telescope[j].rvptrs).residuals,ss.telescope[j].jitter.value,rv.err,/chi2)
+;   print, ss.telescope[j].jittervar.value,ss.telescope[j].jitter.value, sqrt(ss.telescope[j].jittervar.value)
+   rvchi2 = exofast_like((*ss.telescope[j].rvptrs).residuals,ss.telescope[j].jittervar.value,rv.err,/chi2)
 
    if ~finite(rvchi2) then stop
    chi2 += rvchi2
@@ -523,6 +532,26 @@ if (where(ss.planet.fitrv))[0] ne -1 then begin
    endif
 endif
 
+;; Doppler Tomography Model
+for i=0, ss.ndt-1 do begin
+   if keyword_set(psname) then epsname = psname + 'dt.eps'
+   dtchi2 = dopptom_chi2(*(ss.doptom[i].dtptrs),$
+                         ss.planet[(*(ss.doptom[i].dtptrs)).planetndx].tc.value, $
+                         ss.planet[(*(ss.doptom[i].dtptrs)).planetndx].period.value, $
+                         ss.planet[(*(ss.doptom[i].dtptrs)).planetndx].e.value,$
+                         ss.planet[(*(ss.doptom[i].dtptrs)).planetndx].omega.value, $
+                         ss.planet[(*(ss.doptom[i].dtptrs)).planetndx].cosi.value, $
+                         ss.planet[(*(ss.doptom[i].dtptrs)).planetndx].p.value,$
+                         ss.planet[(*(ss.doptom[i].dtptrs)).planetndx].ar.value,$
+                         ss.planet[(*(ss.doptom[i].dtptrs)).planetndx].lambda.value, $
+                         ss.star.logg.value, ss.star.teff.value, ss.star.feh.value,$
+                         ss.star.vsini.value/1d3,ss.star.macturb.value/1d3,$
+                         ss.doptom[i].dtscale.value, debug=ss.debug,/like,psname=epsname)
+
+   if ~finite(dtchi2) then return, !values.d_infinity
+   chi2 += dtchi2
+endfor
+
 ;; Transit model
 for j=0, ntransits-1 do begin
 
@@ -588,7 +617,7 @@ for j=0, ntransits-1 do begin
                                     dilute=band.dilute.value,$
                                     tc=ss.planet[i].tc.value,$
                                     rstar=ss.star.rstar.value/AU,$
-                                    x1=x1,y1=y1,z1=z1) - 1d0)
+                                    x1=x1,y1=y1,z1=z1,au=au) - 1d0)
       endif
 
    endfor
@@ -598,12 +627,18 @@ for j=0, ntransits-1 do begin
    ;; now integrate the model points (before detrending)
    if ninterp gt 1 then modelflux = total(modelflux,2)/ninterp
 
+;   if band.thermal.fit then begin
+;      plot, transitbjd, transit.flux, /ynoz
+;      oplot, transitbjd, modelflux,color='0000ff'x
+;   endif
+
    ;; detrending
    modelflux += total(transit.detrendadd*replicate(1d0,n_elements(transit.bjd))##transit.detrendpars.value,1)  
 
    ;; chi^2
    transitchi2 = exofast_like(transit.flux - modelflux,ss.transit[j].variance.value,transit.err,/chi2)
 
+
    (*ss.transit[j].transitptrs).residuals = transit.flux - modelflux
    (*ss.transit[j].transitptrs).model = modelflux
 

diff --git a/exofast_getmcmcscale.pro b/exofast_getmcmcscale.pro
@@ -240,7 +240,9 @@ if bad[0] ne -1 then mcmcscale[bad,0] = mcmcscale[bad,1]
 bad = where(~finite(mcmcscale[*,1]) or mcmcscale[*,1] eq 0d0)
 if bad[0] ne -1 then mcmcscale[bad,1] = mcmcscale[bad,0]
 bad = where(~finite(mcmcscale) or mcmcscale eq 0d0,nbad)
-if bad[0] ne -1 then return, -1
+
+
+if bad[0] ne -1 then stop;return, -1
 
 ;; return the average of high and low
 return, total(mcmcscale,2)/2.d0

diff --git a/exofast_tran.pro b/exofast_tran.pro
@@ -1,7 +1,11 @@
 function exofast_tran, time, inc, ar, tp, period, e, omega, p, u1, u2, f0, $
                        rstar=rstar, thermal=thermal, reflect=reflect, $
-                       dilute=dilute, tc=tc, q=q,x1=x1,y1=y1,z1=z1
-AU = 215.094177d0
+                       dilute=dilute, tc=tc, q=q,x1=x1,y1=y1,z1=z1, au=au
+
+if n_elements(thermal) eq 0 then thermal = 0
+if n_elements(reflect) eq 0 then reflect = 0
+if n_elements(dilute) eq 0 then dilute = 0
+if n_elements(AU) eq 0 then AU = 215.094177d0
 
 ;; if we have the stellar radius, we can convert time to the
 ;; target's barycentric frame
@@ -12,12 +16,7 @@ endif else transitbjd = time
 
 ;; the impact parameter for each BJD
 z = exofast_getb2(transitbjd, i=inc, a=ar, tperiastron=tp, period=period,$
-                  e=e,omega=omega,z2=depth,x2=x,y2=y,q=q);,x0=x0,y0=y0,z0=z0)
-
-;z = exofast_getb(transitbjd, i=inc, a=ar, tperiastron=tp, period=period,$
-;                  e=e,omega=omega,z=depth,x=x,y=y)
-
-;stop
+                  e=e,omega=omega,z2=depth,x2=x,y2=y,q=q)
 
 if arg_present(z1) then depth += z1
 if arg_present(x1) then x += x1
@@ -32,7 +31,7 @@ if primary[0] ne - 1 then begin
 endif
 
 ;; calculate the fraction of the planet that is visible for each time
-if arg_present(thermal) or arg_present(reflect) then begin
+if thermal ne 0d0 or reflect ne 0d0 then begin
    planetvisible = dblarr(n_elements(time)) + 1d0
    if secondary[0] ne - 1 then begin
       exofast_occultquad, z[secondary]/p, 0, 0, 1d0/p, mu1
@@ -41,17 +40,15 @@ if arg_present(thermal) or arg_present(reflect) then begin
 endif
 
 ;; thermal emission from planet (isotropic)
-if arg_present(thermal) then modelflux += 1d-6*thermal*planetvisible
+if thermal ne 0d0 then modelflux += 1d-6*thermal*planetvisible
 
 ;; phase-dependent reflection off planet
-if arg_present(reflect) then $
+if reflect ne 0d0 then $
    modelflux+=1d-6*reflect*cos(2d0*!dpi*(transitbjd-tc)/period)*planetvisible
 
-;; dilution due to neighboring star
-if arg_present(dilute) then modelflux *= ((1d0-dilute)+dilute)
-
-;; normalization
-modelflux *= f0
+;; normalization and dilution due to neighboring star
+if dilute ne 0d0 then modelflux = f0*(modelflux*(1d0-dilute)+dilute) $
+else modelflux *= f0
 
 return, modelflux