https://github.com/jdeast/EXOFASTv2
Raw File
Tip revision: bba2fcf571e0b213cc33860e95a83d4ae5ff0675 authored by Jason Eastman on 18 July 2022, 20:38:18 UTC
new utility function to retrieve tess data
Tip revision: bba2fcf
pars2step.pro
;; this function translates the starting guesses into the step parameters
function pars2step, ss

au = ss.constants.au/ss.constants.rsun
G = ss.constants.GMsun/ss.constants.RSun^3*ss.constants.day^2 ;; R_Sun^3/(m_sun*day^2)
mjup = ss.constants.gmjupiter/ss.constants.gmsun ;; m_sun
rjup = ss.constants.rjupiter/ss.constants.rsun  ;; r_sun
mearth = ss.constants.gmearth/ss.constants.gmsun ;; m_sun
rearth = ss.constants.rearth/ss.constants.rsun  ;; r_sun
sigmaB = ss.constants.sigmab/ss.constants.lsun*ss.constants.rsun^2 ;; Stefan-Boltzmann constant

;; derive the stellar mass
if ss.star.mstar.userchanged then begin
   ss.star.logmstar.value = alog10( ss.star.mstar.value )
   if ss.star.mstar.scale ne 0d0 then $
      ss.star.logmstar.scale = ss.star.mstar.scale/(alog(10d0)*ss.star.mstar.value)
endif
ss.star.mstar.value = 10^ss.star.logmstar.value

;; derive teff, teffsed, feh, fehsed
if ss.star.teffsed.userchanged and ~ss.star.teff.userchanged then ss.star.teff.value = ss.star.teffsed.value
if ss.star.teff.userchanged and ~ss.star.teffsed.userchanged then ss.star.teffsed.value = ss.star.teff.value
if ss.star.fehsed.userchanged and ~ss.star.feh.userchanged then ss.star.feh.value = ss.star.fehsed.value
if ss.star.feh.userchanged and ~ss.star.fehsed.userchanged then ss.star.fehsed.value = ss.star.feh.value
if ss.star.feh.userchanged and ~ss.star.initfeh.userchanged then ss.star.initfeh.value = ss.star.feh.value

;; derive rstar, rstarsed
if ss.star.rstarsed.userchanged and ~ss.star.rstar.userchanged then begin
   ss.star.rstar.value = ss.star.rstarsed.value
   ss.star.rstar.userchanged = 1B
endif
if ss.star.rstar.userchanged and ~ss.star.rstarsed.userchanged then ss.star.rstarsed.value = ss.star.rstar.value

;; derive the distance from the parallax, if given
if ~ss.star.distance.userchanged and ss.star.parallax.userchanged then $
   ss.star.distance.value = 1d3/ss.star.parallax.value

;; if a starting value for the stellar radius was given, 
;; use it to derive a starting point for the age
if ss.star.rstar.userchanged and ss.yy then begin
   ntries = 100
   age = dindgen(ntries)/(ntries-1)*13.82
   rstars = dblarr(ntries)
   for i=0, ntries-1 do begin
      junk = massradius_yy3(ss.star.mstar.value, ss.star.feh.value, age[i], $
                            ss.star.teff.value,yyrstar=rstar, $
                            sigmab=ss.constants.sigmab/ss.constants.lsun*ss.constants.rsun^2, $
                            gravitysun=ss.constants.gravitysun)
      if ~finite(junk) then stop;return, 0
      rstars[i] = rstar
   endfor
   junk = min(abs(rstars-ss.star.rstar.value),ndx)
   ss.star.age.value = age[ndx]
endif else if ss.yy then begin
   ;; otherwise derive the stellar radius
   junk = massradius_yy3(ss.star.mstar.value, ss.star.feh.value, ss.star.age.value, $
                         ss.star.teff.value,yyrstar=rstar, $
                         sigmab=ss.constants.sigmab/ss.constants.lsun*ss.constants.rsun^2, $
                         gravitysun=ss.constants.gravitysun)
   if ~finite(junk) then stop;return, 0
   ss.star.rstar.value = rstar
endif

;; if not specified, refine the starting value for eep (max out at 808)
if (ss.mist or ss.parsec) and ~ss.star.eep.userchanged and (ss.star.mstar.userchanged or ss.star.rstar.userchanged or ss.star.teff.userchanged) then begin  
   maxeep = 808d0
   mineep = 0d0
   neep = 30d0
   eeps = mineep + (maxeep-mineep)/(neep-1)*dindgen(neep)
   chi2 = dblarr(neep) + !values.d_infinity
   ages = dblarr(neep) + !values.d_infinity

   for i=0L, neep-1 do begin
      if ss.mist then begin
         mistchi2 =  massradius_mist(eeps[i],ss.star.mstar.value,ss.star.initfeh.value,$
                                     ss.star.age.value,ss.star.teff.value,$
                                     ss.star.rstar.value,ss.star.feh.value,mistage=mistage,fitage=ss.star.age.fit,$
                                     /allowold,tefffloor=ss.teffemfloor,fehfloor=ss.fehemfloor,$
                                     rstarfloor=ss.rstaremfloor, agefloor=ss.ageemfloor)
         if finite(mistchi2) then begin
            if mistage gt 13.82d0 then break
            chi2[i] = mistchi2
            ages[i] = mistage
         endif
      endif else if ss.parsec then begin
         parsecchi2 =  massradius_parsec(eeps[i],ss.star.mstar.value,ss.star.initfeh.value,$
                                         ss.star.age.value,ss.star.teff.value,$
                                         ss.star.rstar.value,ss.star.feh.value,$
                                         parsec_age=parsec_age,fitage=ss.star.age.fit,/allowold,$
                                         tefffloor=ss.teffemfloor,fehfloor=ss.fehemfloor,$
                                         rstarfloor=ss.rstaremfloor, agefloor=ss.ageemfloor)
         if finite(parsecchi2) then begin
            if parsec_age gt 13.82d0 then break
            chi2[i] = parsecchi2
            ages[i] = parsec_age
         endif
      endif
   endfor
   minchi2 = min(chi2,ndx)
   if ~finite(minchi2) then begin
      printandlog, 'No EEP between 0 and 808 produces a physical star. Adjust starting stellar parameters', ss.logname
      stop
   endif
   ss.star.eep.value = eeps[ndx]
   if ~ss.star.age.userchanged then ss.star.age.value = ages[ndx]

;; why did I do this? this was reallllly dumb
;; for high mass stars, the default age (4.6 Gyr) will select stars at
;; the end of their lifetimes (capped at carbon burning)
;   repeat begin
;      eep = (mineep + maxeep)/2d0
;      chi2 = massradius_mist(eep,ss.star.mstar.value,ss.star.initfeh.value,$
;                             ss.star.age.value,ss.star.teff.value,$
;                             ss.star.rstar.value,ss.star.feh.value,mistage=mistage,/allowold)
;      if ~finite(chi2) then begin
;         maxeep -= 1
;      endif else begin
;         if mistage gt ss.star.age.value then maxeep = eep $
;         else mineep = eep
;      endelse
;   endrep until abs(maxeep - mineep) lt 0.01
;   if ~finite(chi2) then stop; return, 0 
;   ss.star.eep.value = eep
endif
;; calculate the span of the data
minbjd=!values.d_infinity
maxbjd=-!values.d_infinity
if ss.ntran eq 0 then begin
   for i=0L, ss.ntel-1 do begin
      mint = min((*ss.telescope[i].rvptrs).bjd,max=maxt)
      if mint lt minbjd then minbjd = mint
      if maxt gt maxbjd then maxbjd = maxt
   endfor
endif else begin
   for i=0L, ss.ntran-1 do begin
      mint = min((*ss.transit[i].transitptrs).bjd,max=maxt)
      if mint lt minbjd then minbjd = mint
      if maxt gt maxbjd then maxbjd = maxt
   endfor
endelse

if ~ss.star.vsini.userchanged and (where(ss.doptom.dtscale.fit))[0] ne -1 then begin
   printandlog, 'You must set a prior on vsini', ss.logname
endif

;; derive quantities we'll use later
for i=0, ss.nplanets-1 do begin

   ;; derive the period from log(period)
   if ss.planet[i].period.userchanged then begin
      ss.planet[i].logp.value = alog10(ss.planet[i].period.value)
      
      ;; translate the stepping scale from P to logP
      if ss.planet[i].period.scale ne 0d0 then $
         ss.planet[i].logp.scale = ss.planet[i].period.scale/(alog(10d0)*ss.planet[i].period.value)
   endif else if ss.planet[i].logp.userchanged then begin      
      ss.planet[i].period.value = 10^ss.planet[i].logp.value
   endif else begin
      printandlog, 'You must set the starting value for the period', ss.logname
      stop
   endelse

   ;; how did the user seed e & omega?
   if ss.planet[i].secosw.userchanged or ss.planet[i].sesinw.userchanged then begin
      ss.planet[i].e.value = ss.planet[i].secosw.value^2 + ss.planet[i].sesinw.value^2
      ss.planet[i].omega.value = atan(ss.planet[i].sesinw.value,ss.planet[i].secosw.value)
      ss.planet[i].e.userchanged = 1B
      ss.planet[i].omega.userchanged = 1B
   endif else if ss.planet[i].qecosw.userchanged or ss.planet[i].qesinw.userchanged then begin
      ss.planet[i].e.value = (ss.planet[i].qecosw.value^2 + ss.planet[i].qesinw.value^2)^2
      ss.planet[i].omega.value = atan(ss.planet[i].qesinw.value,ss.planet[i].qecosw.value)      
      ss.planet[i].e.userchanged = 1B
      ss.planet[i].omega.userchanged = 1B
   endif else if ss.planet[i].ecosw.userchanged or ss.planet[i].esinw.userchanged then begin
      ss.planet[i].e.value = sqrt(ss.planet[i].ecosw.value^2 + ss.planet[i].esinw.value^2)
      ss.planet[i].omega.value = atan(ss.planet[i].esinw.value,ss.planet[i].ecosw.value)
      ss.planet[i].e.userchanged = 1B
      ss.planet[i].omega.userchanged = 1B
   endif

   ;; how did the user seed omega?
   if ss.planet[i].omega.userchanged then begin
      ;; do nothing
   endif else if ss.planet[i].omegadeg.userchanged then begin
      ;; translate degrees to radians
      ss.planet[i].omega.value = ss.planet[i].omegadeg.value*!pi/180d0
      if ss.planet[i].omegadeg.scale ne 0d0 then $
         ss.planet[i].omega.scale = ss.planet[i].omegadeg.scale*!pi/180d0
   endif else if ss.planet[i].lsinw.userchanged or ss.planet[i].lcosw.userchanged then begin
      ;; translate Lcosw/Lsinw to omega
;      if ss.planet[i].lsinw2.userchanged then begin
;         ss.planet[i].lsinw.value = ss.planet[i].lsinw2.value
;      endif
      ss.planet[i].omega.value = atan(ss.planet[i].lsinw.value,ss.planet[i].lcosw.value)
      ss.planet[i].omega.userchanged = 1B
   endif

   ;; how did the user seed e?
   if ss.planet[i].e.userchanged then begin
      ;; do nothing
   endif else if ss.planet[i].vcve.userchanged then begin

      ;; pick the omega that maximizes the range of vcve
      if ~ss.planet[i].omega.userchanged then begin
         if ss.planet[i].vcve.value le 1d0 then ss.planet[i].omega.value = !dpi/2d0 $
         else ss.planet[i].omega.value = -!dpi/2d0
         ss.planet[i].omega.userchanged = 1B
      endif

;      if ss.planet[i].sign2.userchanged then begin
      if ss.planet[i].sign.userchanged then begin
         sign = floor(ss.planet[i].sign.value)
      endif

      ;; rederive lsinw/lcosw
      if ~ss.planet[i].lsinw.userchanged and ~ss.planet[i].lcosw.userchanged then begin
         if floor(ss.planet[i].sign.value) then L = sqrt(0.75d0) $
         else L = sqrt(0.25d0)
         ss.planet[i].lsinw.value = L*sin(ss.planet[i].omega.value)
         ss.planet[i].lcosw.value = L*cos(ss.planet[i].omega.value)
      endif

      ;; derive e (pick the sign if not chosen)
      ss.planet[i].e.value = vcve2e(ss.planet[i].vcve.value,omega=ss.planet[i].omega.value, sign=sign)
      ss.planet[i].sign.value = sign

   endif 

   ;; error checking on e/omega
   if ss.planet[i].e.value ge 1d0 or ss.planet[i].e.value lt 0d0 or ~finite(ss.planet[i].e.value) then begin
      printandlog, 'The starting value for eccentricity_' + strtrim(i,2) + ' is not physical. If you are seeing this with an automatically generated prior file after updating EXOFASTv2, delete the starting points for sign, vcve, lsinw, lcosw in the prior file and use the median e, omegadeg (from the previous .tex file) instead.', ss.logname
      stop
   endif

   if ~finite(ss.planet[i].omega.value) then begin
      printandlog, 'The starting value for omega_' + strtrim(i,2) + ' is not physical. If you are seeing this with an automatically generated prior file after updating EXOFASTv2, delete the the starting points for sign, vcve, lsinw, lcosw in the prior file and use the median e, omegadeg (from the previous .tex file) instead.', ss.logname
      stop
   endif

   ;; translate from e/omega to various e/omega parameterizations
   if ~ss.planet[i].qecosw.userchanged then ss.planet[i].qecosw.value = (ss.planet[i].e.value)^(0.25d0)*cos(ss.planet[i].omega.value)
   if ~ss.planet[i].qesinw.userchanged then ss.planet[i].qesinw.value = (ss.planet[i].e.value)^(0.25d0)*sin(ss.planet[i].omega.value)
   if ~ss.planet[i].secosw.userchanged then ss.planet[i].secosw.value = sqrt(ss.planet[i].e.value)*cos(ss.planet[i].omega.value)
   if ~ss.planet[i].sesinw.userchanged then ss.planet[i].sesinw.value = sqrt(ss.planet[i].e.value)*sin(ss.planet[i].omega.value)
   if ~ss.planet[i].ecosw.userchanged then ss.planet[i].ecosw.value = ss.planet[i].e.value*cos(ss.planet[i].omega.value)
   if ~ss.planet[i].esinw.userchanged then ss.planet[i].esinw.value = ss.planet[i].e.value*sin(ss.planet[i].omega.value)
   if ~ss.planet[i].vcve.userchanged then ss.planet[i].vcve.value = sqrt(1d0-ss.planet[i].e.value^2)/(1d0+ss.planet[i].esinw.value)

   ;; if the user sets e or omega, 
   ;; we need to make sure our combination of vcve/lcosw/lsinw/lsinw2/sign will recover it
   e1 = vcve2e(ss.planet[i].vcve.value,omega=ss.planet[i].omega.value, sign=0B)
   e2 = vcve2e(ss.planet[i].vcve.value,omega=ss.planet[i].omega.value, sign=1B)
  
   if abs(ss.planet[i].e.value - e1) lt 1d-8 then ss.planet[i].sign.value = 0.5 $
   else if abs(ss.planet[i].e.value - e2) lt 1d-8 then ss.planet[i].sign.value = 1.5 $
   else begin
      printandlog, 'ERROR translating starting e and omega to parameterization', ss.logname
      stop
   endelse

   if ~ss.planet[i].lsinw.userchanged and ~ss.planet[i].lcosw.userchanged then begin
      if floor(ss.planet[i].sign.value) then L = sqrt(0.75d0) $
      else L = sqrt(0.25d0)
      ss.planet[i].lsinw.value = L*sin(ss.planet[i].omega.value)
      ss.planet[i].lcosw.value = L*cos(ss.planet[i].omega.value)
   endif

if 0 then begin
   sign = floor(ss.planet[i].sign2.value)
   if ss.planet[i].vcve.value lt 1d0 and ss.planet[i].lsinw.value le 0d0 and sign then begin
      ;; then I need lsinw to be negative and lsinw2 to be positive
      ss.planet[i].lsinw2.value = -ss.planet[i].lsinw.value
      ;; and I need sign to be > 1
      ss.planet[i].sign2.value = 0.5d0
   endif else if ss.planet[i].vcve.value ge 1d0 and ss.planet[i].lsinw.value gt 0d0 and ~sign then begin
      ;; then I need lsinw to be positive and lsinw2 to be negative
      ss.planet[i].lsinw2.value = -ss.planet[i].lsinw.value
      ;; and I need sign to be > 1
      ss.planet[i].sign2.value = 0.5d0
   endif else begin
      ss.planet[i].sign2.value = ss.planet[i].sign.value
      ss.planet[i].lsinw2.value = ss.planet[i].lsinw.value
   endelse
endif else begin
      ss.planet[i].sign2.value = ss.planet[i].sign.value
      ss.planet[i].lsinw2.value = ss.planet[i].lsinw.value
endelse

   ;; how did the user seed rp/rstar?
   if ss.planet[i].p.userchanged then begin
      if ~ss.planet[i].rpearth.userchanged then $
         ss.planet[i].rpearth.value = ss.planet[i].p.value*ss.star.rstar.value/rearth
      ss.planet[i].rpearth.userchanged = 1B
   endif else if ss.planet[i].rpearth.userchanged then begin
      ss.planet[i].p.value = ss.planet[i].rpearth.value*rearth/ss.star.rstar.value
   endif else if ss.planet[i].rp.userchanged then begin
      ss.planet[i].p.value = ss.planet[i].rp.value*rjup/ss.star.rstar.value
      if ~ss.planet[i].rpearth.userchanged then $
         ss.planet[i].rpearth.value = ss.planet[i].p.value*ss.star.rstar.value/rearth
      ss.planet[i].rpearth.userchanged = 1B
   endif else if ss.planet[i].mpearth.userchanged then begin 
      ss.planet[i].rpearth.value = massradius_chen(ss.planet[i].mpearth.value)
      ss.planet[i].p.value = ss.planet[i].rpearth.value*rearth/ss.star.rstar.value
   endif

   ;; how did the user seed the inclination?
   if ss.planet[i].i.userchanged then begin
      sini = sin(ss.planet[i].i.value)
   endif else if ss.planet[i].ideg.userchanged then begin      
     ss.planet[i].i.value = ss.planet[i].ideg.value*!dpi/180d0     
     sini = sin(ss.planet[i].i.value)
   endif else if ss.planet[i].chord.userchanged then begin
      if ~ss.planet[i].chord.fit then $
         printandlog, 'changing starting value for chord when not fitting chord is not supported', ss.logname
;      ;; need to know rp/rstar, a/rstar, mpsun, inc (!) first
;      ss.planet[i].b.value = sqrt((1d0+ss.planet[i].p.value)^2-ss.planet[i].chord.value^2)
;      ss.planet[i].cosi.value = ss.planet[i].b.value/$
;                                (ss.planet[i].ar.value*(1d0-ss.planet[i].e.value^2)/(1d0+ss.planet[i].esinw.value))
;      ss.planet[i].i.value = acos(ss.planet[i].cosi.value)
      sini = 1d0
   endif else if ss.planet[i].b.userchanged then begin
      if ~ss.planet[i].b.fit then $
         printandlog, 'changing starting value for b when not fitting b is not supported', ss.logname
;      ;; need to know a/rstar, mpsun, inc (!) first
;      ss.planet[i].cosi.value = ss.planet[i].b.value/$
;                                (ss.planet[i].ar.value*(1d0-ss.planet[i].e.value^2)/(1d0+ss.planet[i].esinw.value))
;      ss.planet[i].i.value = acos(ss.planet[i].cosi.value)
      sini = 1d0
   endif else begin
      ;; use the default cosi starting value
      ss.planet[i].i.value = acos(ss.planet[i].cosi.value)
      sini = sin(ss.planet[i].i.value)
   endelse      
  
   ;; how did the user seed the planet mass?
   if ss.planet[i].mpsun.userchanged then begin
      ss.planet[i].logmp.value = alog10(ss.planet[i].mpsun.value)
      ;; translate the stepping scale from mpsun to logmp
      if ss.planet[i].logmp.scale ne 0d0 then $
         ss.planet[i].logmp.scale = ss.planet[i].mpsun.scale/(alog(10d0)*ss.planet[i].mpsun.value)
   endif else if ss.planet[i].logmp.userchanged then begin
      ss.planet[i].mpsun.value = 10d0^ss.planet[i].logmp.value
      ss.planet[i].logmp.value = alog10(ss.planet[i].mpsun.value)
   endif else if ss.planet[i].mp.userchanged then begin
      ;; convert from Jupiter masses to Solar masses
      ss.planet[i].mpsun.value = ss.planet[i].mp.value*mjup
      ss.planet[i].logmp.value = alog10(ss.planet[i].mpsun.value)
   endif else if ss.planet[i].mpearth.userchanged then begin
      ;; convert from Earth masses to Solar masses
      ss.planet[i].mpsun.value = ss.planet[i].mp.value*mearth
      ss.planet[i].logmp.value = alog10(ss.planet[i].mpsun.value)
   endif else if ss.planet[i].msiniearth.userchanged then begin
      ;; convert from Earth masses to Solar masses
      ss.planet[i].mpsun.value = ss.planet[i].msiniearth.value*mearth/sini
      ss.planet[i].logmp.value = alog10(ss.planet[i].mpsun.value)
   endif else if ss.planet[i].msini.userchanged then begin
      ;; convert from Jupiter masses to Solar masses
      ss.planet[i].mpsun.value = ss.planet[i].msini.value*mjup/sini
      ss.planet[i].logmp.value = alog10(ss.planet[i].mpsun.value)
   endif else if ss.planet[i].logk.userchanged then begin
      ss.planet[i].k.value = 10d0^ss.planet[i].logk.value
      ;; now convert from K to mpsun
      ss.planet[i].mpsun.value = ktom2(ss.planet[i].K.value, ss.planet[i].e.value,$
                                       ss.planet[i].i.value, ss.planet[i].period.value, $
                                       ss.star.mstar.value, GMsun=ss.constants.GMsun/1d6) ;; m_sun
      ss.planet[i].logmp.value = alog10(ss.planet[i].mpsun.value)
   endif else if ss.planet[i].k.userchanged then begin
      if ss.planet[i].k.value lt 0d0 then begin
         ss.planet[i].mpsun.value = -ktom2(-ss.planet[i].K.value, ss.planet[i].e.value,$
                                           ss.planet[i].i.value, ss.planet[i].period.value, $
                                           ss.star.mstar.value, GMsun=ss.constants.GMsun/1d6) ;; m_sun
      endif else begin
         ss.planet[i].mpsun.value = ktom2(ss.planet[i].K.value, ss.planet[i].e.value,$
                                          ss.planet[i].i.value, ss.planet[i].period.value, $
                                          ss.star.mstar.value, GMsun=ss.constants.GMsun/1d6) ;; m_sun
      endelse
      ss.planet[i].logmp.value = alog10(ss.planet[i].mpsun.value)
   endif else if ss.planet[i].rpearth.userchanged then begin
      ss.planet[i].mpearth.value = massradius_chenreverse(ss.planet[i].rpearth.value)
      ss.planet[i].mpsun.value = ss.planet[i].mpearth.value*mearth ;; m_sun
   endif

   ;; how did the user seed a/rstar?
   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)                       ;; semi-major axis in r_sun
   ss.planet[i].ar.value = ss.planet[i].arsun.value/ss.star.rstar.value    ;; a/rstar (unitless)
   ss.planet[i].a.value = ss.planet[i].arsun.value/AU ;; semi major axis in AU
   ss.planet[i].cosi.scale = 1d0/ss.planet[i].ar.value

   ;; translate from bigomegadeg to bigomega
   if ss.planet[i].bigomegadeg.userchanged then begin
      ss.planet[i].bigomega.value = ss.planet[i].bigomegadeg.value*!pi/180d0
      if ss.planet[i].bigomegadeg.scale ne 0 then $
         ss.planet[i].bigomega.scale = ss.planet[i].bigomegadeg.scale*!pi/180d0
   endif

   ;; translate from bigomega to lsinbigomega/lcosbigomega
   if ss.planet[i].bigomega.userchanged then begin
      ss.planet[i].lsinbigomega.value = 0.5d0*sin(ss.planet[i].bigomega.value) 
      ss.planet[i].lcosbigomega.value = 0.5d0*cos(ss.planet[i].bigomega.value) 
   endif

   ;; translate from lambdadeg to lambda
   if ss.planet[i].lambdadeg.userchanged then begin
      ss.planet[i].lambda.value = ss.planet[i].lambdadeg.value*!pi/180d0
      if ss.planet[i].lambdadeg.scale ne 0 then $
         ss.planet[i].lambda.scale = ss.planet[i].lambdadeg.scale*!pi/180d0
   endif

   ;; translate from lambda to lsinlambda/lcoslambda
   if ss.planet[i].lambda.userchanged then begin
      ss.planet[i].lsinlambda.value = sin(ss.planet[i].lambda.value) 
      ss.planet[i].lcoslambda.value = cos(ss.planet[i].lambda.value) 
   endif

   ;; translate to tt
   if ss.planet[i].tt.fit then begin
      if ss.planet[i].tt.userchanged then begin
         ;; do nothing; it was already set explicitly via the prior file
      endif else if ss.planet[i].tc.userchanged then begin
         ;; derive tt from tc
         ss.planet[i].tt.value = tc2tt(ss.planet[i].tc.value,$
                                       ss.planet[i].e.value,$
                                       ss.planet[i].i.value,$
                                       ss.planet[i].omega.value,$
                                       ss.planet[i].period.value)

         ;; make tt as close as possible to tc
         nper = round((ss.planet[i].tc.value - ss.planet[i].tt.value)/ss.planet[i].period.value)
         ss.planet[i].tt.value += nper*ss.planet[i].period.value

      endif else if ss.planet[i].tp.userchanged then begin
         ;; derive tp->tc->tt
         phase = exofast_getphase(ss.planet[i].e.value,ss.planet[i].omega.value, /primary)
         ss.planet[i].tc.value = ss.planet[i].tp.value + ss.planet[i].period.value*phase
         ss.planet[i].tt.value = tc2tt(ss.planet[i].tc.value,$
                                       ss.planet[i].e.value,$
                                       ss.planet[i].i.value,$
                                       ss.planet[i].omega.value,$
                                       ss.planet[i].period.value)

         ;; make tt as close as possible to tp
         nper = round((ss.planet[i].tp.value - ss.planet[i].tt.value)/ss.planet[i].period.value)
         ss.planet[i].tt.value += nper*ss.planet[i].period.value

         ;; make tc as close as possible to tp
         nper = round((ss.planet[i].tp.value - ss.planet[i].tc.value)/ss.planet[i].period.value)
         ss.planet[i].tc.value += nper*ss.planet[i].period.value

      endif else begin
         ;; otherwise guess it's in the middle of the data
         ss.planet[i].tt.value = (maxbjd + minbjd)/2d0
         printandlog, 'Warning: No starting guess for TT (or TC or TP) given; assuming the middle of the data', ss.logname
      endelse
      if ss.planet[i].tt.prior eq 0d0 then ss.planet[i].tt.prior = ss.planet[i].tt.value 
   endif else begin
      ;; fit tc
      if ss.planet[i].tc.userchanged then begin
         ;; do nothing; it was already set explicitly via the prior file
      endif else if ss.planet[i].tp.userchanged then begin
         ;; derive tp->tc
         phase = exofast_getphase(ss.planet[i].e.value,ss.planet[i].omega.value, /primary)
         ss.planet[i].tc.value = ss.planet[i].tp.value + ss.planet[i].period.value*phase
         ss.planet[i].tc.prior = ss.planet[i].tc.value 
      endif else if ss.planet[i].tt.userchanged then begin
         ss.planet[i].tc.value = tc2tt(ss.planet[i].tt.value,$
                                       ss.planet[i].e.value,$
                                       ss.planet[i].i.value,$
                                       ss.planet[i].omega.value,$
                                       ss.planet[i].period.value,/reverse_correction)

         ;; make tt as close as possible to tc
         nper = round((ss.planet[i].tt.value - ss.planet[i].tc.value)/ss.planet[i].period.value)
         ss.planet[i].tc.value += nper*ss.planet[i].period.value
      endif else begin
         ;; otherwise guess it's in the middle of the data
         ss.planet[i].tc.value = (maxbjd + minbjd)/2d0
         printandlog, 'Warning: No starting guess for TC (or TT or TP) given; assuming the middle of the data', ss.logname
      endelse
      if ss.planet[i].tc.prior eq 0d0 then ss.planet[i].tc.prior = ss.planet[i].tc.value 
   endelse

   ;; approximate duration taken from Winn 2010 (close enough; these should only be used to schedule observations anyway)
   ss.planet[i].b.value = ss.planet[i].ar.value*ss.planet[i].cosi.value*(1d0-ss.planet[i].e.value^2)/(1d0+ss.planet[i].esinw.value)  ;; eq 7, Winn 2010
   ss.planet[i].t14.value = ss.planet[i].period.value/!dpi*asin(sqrt((1d0+abs(ss.planet[i].p.value))^2 - ss.planet[i].b.value^2)/(sini*ss.planet[i].ar.value))*sqrt(1d0-ss.planet[i].e.value^2)/(1d0+ss.planet[i].esinw.value)


   if ~ss.planet[i].chord.userchanged then ss.planet[i].chord.value = sqrt((1d0+ss.planet[i].p.value)^2 - ss.planet[i].b.value^2)

endfor

return, 1

end
back to top