Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
merging new pro files
  • Loading branch information
jdeast01 committed May 15, 2017
1 parent ff28540 commit cb5c44c
Show file tree
Hide file tree
Showing 13 changed files with 214 additions and 121 deletions.
26 changes: 13 additions & 13 deletions derivepars.pro
Expand Up @@ -61,18 +61,18 @@ for i=0, ss.nplanets-1 do begin
ss.planet[i].period.value[positive],$
ss.star.mstar.value[positive])

ss.planet[i].mp.value = ss.planet[i].mpsun.value/mjup
ss.planet[i].msini.value = ss.planet[i].mp.value*sini
ss.planet[i].mpearth.value = ss.planet[i].mpsun.value/mearth
ss.planet[i].msiniearth.value = ss.planet[i].mpearth.value*sini
ss.planet[i].q.value = ss.planet[i].mpsun.value/ss.star.mstar.value
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].arsun.value=(G*(ss.star.mstar.value+ss.planet[i].mp.value)*ss.planet[i].period.value^2/(4d0*!dpi^2))^(1d0/3d0)
ss.planet[i].ar.value = ss.planet[i].arsun.value/ss.star.rstar.value ;; unitless
ss.planet[i].a.value = ss.planet[i].arsun.value/AU ;; AU
ss.planet[i].rpsun.value = ss.planet[i].p.value*ss.star.rstar.value
ss.planet[i].rp.value = ss.planet[i].rpsun.value/rjup
ss.planet[i].rpearth.value = ss.planet[i].rpsun.value/rearth
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].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

;; time of periastron
ss.planet[i].phase.value=exofast_getphase(ss.planet[i].e.value,ss.planet[i].omega.value,/pri)
Expand Down Expand Up @@ -151,12 +151,12 @@ for i=0, ss.nplanets-1 do begin
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].mp.value/ss.planet[i].rp.value^2*9.31686171d0) ;; cgs
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

;; depth != delta if grazing (ignore limb darkening)
ss.planet[i].delta.value = ss.planet[i].p.value^2
for j=0, ss.nsteps-1 do begin
for j=0L, ss.nsteps-1L do begin
exofast_occultquad, ss.planet[i].b.value[j], 0d0, 0d0, ss.planet[i].p.value[j],mu1
ss.planet[i].depth.value[j] = 1d0-mu1[0]
endfor
Expand Down
45 changes: 35 additions & 10 deletions exofast_chi2v2.pro
Expand Up @@ -246,10 +246,12 @@ 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
chi2 += 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')
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
chi2 += massradius_yy3(ss.star.mstar.value, ss.star.feh.value, ss.star.age.value, ss.star.teff.value,yyrstar=rstar, debug=ss.debug)
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
chi2 += yychi2
;print, 'YY penalty = ' + strtrim(yychi2,2)
if ~finite(chi2) then begin
if ss.debug then print, 'star is bad'
return, !values.d_infinity
Expand All @@ -263,6 +265,11 @@ if ss.star.errscale.value le 0 then begin
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
endif

if ss.amoeba then begin
;; need some derived parameters
;; (but others will change after we change logk or p, so
Expand Down Expand Up @@ -437,11 +444,11 @@ 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,verbose=ss.debug, f0=f, fp0=fp, ep0=ep, psname=psname+'.sed.eps')
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,verbose=ss.debug, f0=f, fp0=fp, ep0=ep)
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 ~finite(sedchi2) then begin
Expand All @@ -455,6 +462,7 @@ if file_test(ss.star.fluxfile) then begin
return, !values.d_infinity
endif
chi2 += sedchi2
; print, 'SED penalty = ' + strtrim(sedchi2,2)
endif

;; RV model (non-interacting planets)
Expand Down Expand Up @@ -502,6 +510,7 @@ for j=0, ntelescopes-1 do begin

if ~finite(rvchi2) then stop
chi2 += rvchi2
; print, 'RV penalty = ' + strtrim(rvchi2,2)
; chi2 += total(((rv.rv - modelrv)/rv.err)^2)
endfor

Expand Down Expand Up @@ -532,6 +541,8 @@ for j=0, ntransits-1 do begin
u2err = 0.05d0
chi2 += ((band.u1.value-u1claret)/u1err)^2
chi2 += ((band.u2.value-u2claret)/u2err)^2
; print, 'u1 penalty = ' + strtrim(((band.u1.value-u1claret)/u1err)^2,2)
; print, 'u2 penalty = ' + strtrim(((band.u2.value-u2claret)/u2err)^2,2)
endif

;; Kepler Long candence data; create several model points and average
Expand All @@ -547,6 +558,14 @@ for j=0, ntransits-1 do begin
modelflux = dblarr(npoints) + 1d0
endelse

;; get the motion of the star due to the planet
junk = exofast_getb2(transitbjd,inc=ss.planet.i.value,a=ss.planet.ar.value,$
tperiastron=ss.planet.tp.value,$
period=ss.planet.period.value,$
e=ss.planet.e.value,omega=ss.planet.omega.value,$
q=ss.star.mstar.value/ss.planet.mpsun.value,$
x1=x1,y1=y1,z1=z1)

for i=0, nplanets-1 do begin
if ss.planet[i].fittran then begin

Expand All @@ -568,8 +587,8 @@ for j=0, ntransits-1 do begin
reflect=band.reflect.value, $
dilute=band.dilute.value,$
tc=ss.planet[i].tc.value,$
rstar=ss.star.rstar.value/AU) - 1d0)

rstar=ss.star.rstar.value/AU,$
x1=x1,y1=y1,z1=z1) - 1d0)
endif

endfor
Expand All @@ -584,6 +603,7 @@ for j=0, ntransits-1 do begin

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

Expand All @@ -593,11 +613,14 @@ for j=0, ntransits-1 do begin
forprint, transitbjd, modelflux, format='(f0.8,x,f0.6)', textout=psname + '.model.' + strtrim(j,2) + '.flux.txt', /nocomment,/silent

if ~finite(transitchi2) then stop
chi2 += transitchi2
; print, 'transit penalty: ' + strtrim(transitchi2,2)

; chi2 += total(((transit.flux - modelflux)/transit.err)^2)
;transitchi2 = strtrim(total(((transit.flux - modelflux)/transit.err)^2),2)

chi2 += transitchi2
; print, 'transit penalty: ' + strtrim(transitchi2,2) + ' ' + strtrim(ss.transit[j].variance.value,2)

;print, 'transit chi2 = ' + strtrim(total(((transit.flux - modelflux)/transit.err)^2),2)
;stop
; screen = GET_SCREEN_SIZE()
; if win_state(5) then wset, 5 $
; else window, 5, xsize=xsize, ysize=ysize, xpos=screen[0]/3d0, ypos=0
Expand Down Expand Up @@ -638,6 +661,9 @@ if ss.ttvs then begin
chi2 += ((coeffs[0]-ss.planet[i].tc.value)/sigma[0])^2
chi2 += ((coeffs[1]-ss.planet[i].period.value)/sigma[1])^2

; print, ((coeffs[0]-ss.planet[i].tc.value)/sigma[0])^2
; print, ((coeffs[0]-ss.planet[i].tc.value)/sigma[0])^2

if ss.debug or keyword_set(psname) then begin
if keyword_set(psname) then begin
;; astrobetter.com tip on making pretty IDL plots
Expand Down Expand Up @@ -727,7 +753,6 @@ if ss.debug then print, ss.star.rstar.value, pars, chi2, format='(' + strtrim(n_
;; if this stop is triggered, you've found a bug!!
if ~finite(chi2) then stop


return, chi2

end
Expand Down
72 changes: 59 additions & 13 deletions exofast_demc.pro
Expand Up @@ -155,7 +155,7 @@ pro exofast_demc, bestpars,chi2func,pars,chi2=chi2, tofit=tofit,$
nchains=nchains, angular=angular,$
burnndx=burnndx, removeburn=removeburn, $
gelmanrubin=gelmanrubin, tz=tz, maxgr=maxgr, mintz=mintz, $
derived=derived, resumendx=resumendx
derived=derived, resumendx=resumendx,stretch=stretch

;; default values
if n_elements(maxsteps) eq 0 then maxsteps = 1d5
Expand Down Expand Up @@ -223,8 +223,9 @@ for j=0, nchains-1 do begin
;; start 5 steps from best value. If that's not allowed
;; (infinite chi^2), slowly approach 0 steps away from the best value
;; niter should never be larger than ~5000
if j eq 0 then factor = 0d0 else factor = 1d0 ;; first chain starts at best fit
pars[0:nfit-1,0,j] = bestpars[tofit] + $
5d0/exp(niter/1000d)*scale*call_function(randomfunc,seed,nfit,/normal)
5d0/exp(niter/1000d)*scale*call_function(randomfunc,seed,nfit,/normal)*factor
newpars[tofit] = pars[0:nfit-1,0,j]
;; find the chi^2
chi2[0,j] = call_function(chi2func, newpars, determinant=det, derived=dpar)
Expand All @@ -247,6 +248,8 @@ print

if n_elements(dpar) ne 0 then oldderived = reform(derived[*,0,*],nderived,nchains)

a = 2d0 ;; a > 1

nextrecalc = 100L
npass = 1L
nstop = 0L
Expand All @@ -255,6 +258,9 @@ tz0 = 0d0
tzsteps = 0L
alreadywarned = 0L
t0 = systime(/seconds)
minburnndx = 0L

idealacceptancerate = 0.23d0 ;; automatically scale "a" to get ideal acceptance rate

;; start MCMC chain
for i=resumendx,maxsteps-1L do begin
Expand All @@ -266,18 +272,38 @@ for i=resumendx,maxsteps-1L do begin
;; automatically thin the chain (saves memory)
for k=0L, nthin-1L do begin

;; differential evolution mcmc step -- see Ter Braak 2006
repeat r1=floor(call_function(randomfunc,seed)*(nchains)) $
until r1 ne j ;; a random chain
repeat r2=floor(call_function(randomfunc,seed)*(nchains)) $
until r2 ne j and r2 ne r1 ;; another random chain
newpars[tofit] = oldpars[tofit] + $
gamma*(pars[*,i-1,r1]-pars[*,i-1,r2]) + $
(call_function(randomfunc,seed,nfit)-0.5d0)*scale/10d0

if keyword_set(stretch) then begin
;; affine invariant "stretch move" step
;; see Goodman & Weare, 2010 (G10), eq 7
;; http://msp.org/camcos/2010/5-1/camcos-v5-n1-p04-p.pdf
repeat r1=floor(call_function(randomfunc,seed)*(nchains)) $
until r1 ne j ;; a random chain

;; a random number between 1/a and a, weighted toward lower
;; values to satisfy the symmetry condition (G10, eq 8)
z = ((a-1d0)*call_function(randomfunc,seed) + 1d0)^2d0/2d0

;; use the most recent step (Goodman, pg 69)
if r1 lt j then ndx = i else ndx = i-1
;; define the new step (eq 7 in G10)
newpars[tofit] = pars[*,ndx,r1] + z*(oldpars[tofit]-pars[*,ndx,r1])
fac = z^(nfit-1) ;; GW10, last equation, pg 70
endif else begin
;; differential evolution mcmc step (Ter Braak, 2006)
repeat r1=floor(call_function(randomfunc,seed)*(nchains)) $
until r1 ne j ;; a random chain
repeat r2=floor(call_function(randomfunc,seed)*(nchains)) $
until r2 ne j and r2 ne r1 ;; another random chain
newpars[tofit] = oldpars[tofit] + $
gamma*(pars[*,i-1,r1]-pars[*,i-1,r2]) + $
(call_function(randomfunc,seed,nfit)-0.5d0)*scale/10d0
fac = 1d0
endelse


;; calculate the chi^2 of the new step
newchi2 = call_function(chi2func,newpars,determinant=det,derived=dpar)
C = olddet[j]/det*exp((oldchi2 - newchi2)/2d0)
C = fac*olddet[j]/det*exp((oldchi2 - newchi2)/2d0)

;; accept the step; update values
if call_function(randomfunc,seed) lt C then begin
Expand Down Expand Up @@ -378,6 +404,26 @@ for i=resumendx,maxsteps-1L do begin
;; print the progress, cumulative acceptance rate, and time remaining
acceptancerate = strtrim(string(naccept/double(i*nchains*nthin)*100,$
format='(f6.2)'),2)

;; Limited testing suggests this is unnecessary (counter-productive, even)
; ;; tune the stretch scale (a) to achieve the optimal acceptance rate (~23%)
; ;; if we've taken fewer than 5% of the total steps
; if double(i)/maxsteps lt 0.05 and keyword_set(stretch) then begin
; if naccept/(double(i*nchains*nthin)) lt idealacceptancerate*0.75 or $
; naccept/(double(i*nchains*nthin)) gt idealacceptancerate*1.25 then begin
;
; print
; print, 'Acceptance rate (' + acceptancerate + ') not ideal -- scaling stretch scale from ' + strtrim(a,2) + ' to ' + strtrim(a*naccept/(double(i*nchains*nthin))/idealacceptancerate,2)
; a = a*naccept/(double(i*nchains*nthin))/idealacceptancerate
;
; ;; more continous scaling
; if accepted then a += (1d0-idealacceptancerate)/i $
; else a -= idealacceptancerate/i
;
; minburnndx = i
; endif
; endif

timeleft = (systime(/seconds)-t0)*(maxsteps/(i+1d)-1d)
units = ' seconds '
if timeleft gt 60 then begin
Expand Down Expand Up @@ -418,7 +464,7 @@ medchi2 = min(median(chi2[0.1*nstop:nstop,*],dimension=1))
;; good chains must have some values below this median
goodchains = where(minchi2 lt medchi2,ngood)

burnndx = 0L
burnndx = minburnndx
for j=0L, ngood-1 do begin
tmpndx = (where(chi2[0:nstop,goodchains[j]] lt medchi2))(0)
if tmpndx gt burnndx then burnndx = tmpndx
Expand Down
8 changes: 2 additions & 6 deletions exofast_getb.pro
Expand Up @@ -53,9 +53,7 @@

function exofast_getb, bjd, i=i, a=a, tperiastron=tperiastron, Period=P, $
e=e, omega=omega, x=x, y=y, z=z, $
lonascnode=lonascnode, q=q

if n_elements(q) eq 0 then q = !values.d_infinity
lonascnode=lonascnode

;; calculate the mean anomaly corresponding to each observed time
meananom = (2.d0*!dpi*(1.d0 + (bjd - Tperiastron)/P)) mod (2.d0*!dpi)
Expand All @@ -71,10 +69,8 @@ endif else begin
e=0.d0
endelse

atot = a ;; a1 + a2

;; calculate the corresponding (x,y) coordinates of planet
r = atot*(1d0-e^2)/(1d0+e*cos(trueanom))
r = a*(1d0-e^2)/(1d0+e*cos(trueanom))

;; as seen from observer
x = -r*cos(trueanom + omega)
Expand Down

0 comments on commit cb5c44c

Please sign in to comment.