adding more pro files

exofast_plotchains.pro

@@ -0,0 +1,145 @@
+pro exofast_plotchains, ss, chainfile=chainfile
+
+if n_elements(chainfile) eq 0 then chainfile = 'chain.ps'
+
+burnndx = ss.burnndx
+nsteps = ss.nsteps/ss.nchains
+medndx = (nsteps-ss.burnndx)/2
+halfsigma = erf(1d/sqrt(2d0))/2d
+lowsigndx = round(((nsteps-ss.burnndx)/2.d0 - (nsteps-ss.burnndx)*halfsigma)*ss.nchains)
+hisigndx = round(((nsteps-ss.burnndx)/2.d0 + (nsteps-ss.burnndx)*halfsigma)*ss.nchains)
+;; prepare the postscript device
+mydevice=!d.name
+set_plot, 'PS'
+aspect_ratio=1.5
+xsize = 18
+ysize=xsize/aspect_ratio
+!p.font=0
+device, filename=chainfile
+device, set_font='Times',/tt_font
+device, xsize=xsize,ysize=ysize,/color,bits=24
+charsize=1.5
+loadct,39,/silent
+red = 254
+
+allpars = []
+parnames = []
+medianpars = []
+!p.multi=[0,2,4] ;; 8 to a page
+
+chi2 = reform((*ss.chi2),ss.nsteps/ss.nchains,ss.nchains)
+ymin = min(chi2[ss.burnndx:-1,*],max=ymax)
+xtitle='!3 Chain link number'
+if ymin lt 0 then sign = ' + ' else sign = ' - '
+ytitle='!3' + textoidl('\chi^2') + sign + strtrim(abs(ymin),2)
+plot, [0],[0], psym=10, xtitle=xtitle, ytitle=ytitle,$
+      charsize=charsize,xstyle=1, yrange=[0,ymax-ymin], xrange=[0,nsteps],font=1
+for l=0L, ss.nchains-1L do $
+   oplot, chi2[*,l] - ymin, color=l*255d0/ss.nchains ;, transparency=100d0-100d0/ss.nchains
+oplot, [ss.burnndx,ss.burnndx],[-9d99,9d99]
+print, ymin
+
+for i=0, n_tags(ss)-1 do begin
+   for j=0, n_elements(ss.(i))-1 do begin
+      for k=0, n_tags(ss.(i)[j])-1 do begin
+
+         derive = 0B
+         m=-1
+         ;; this captures the detrending variables
+         if (size(ss.(i)[j].(k)))[1] eq 10 then begin
+            if (size(ss.(i)[j].(k)))[0] ne 0 then begin
+               for l=0L, n_tags(*(ss.(i)[j].(k)))-1 do begin
+                  if (size((*(ss.(i)[j].(k))).(l)))[2] eq 8 then begin 
+                     for m=0L, n_elements((*(ss.(i)[j].(k))).(l))-1 do begin
+                        if tag_exist((*(ss.(i)[j].(k))).(l)[m],'derive') then begin
+                           if (*(ss.(i)[j].(k))).(l)[m].derive then begin
+                              pars = (reform((*(ss.(i)[j].(k))).(l)[m].value,nsteps,ss.nchains))
+;                              pars = (*(ss.(i)[j].(k))).(l)[m].value[burnndx:*,*]
+                              derive = 1B
+                              label = (*(ss.(i)[j].(k))).(l)[m].label
+                              unit = (*(ss.(i)[j].(k))).(l)[m].unit
+                              latex = (*(ss.(i)[j].(k))).(l)[m].latex
+                              best = (*(ss.(i)[j].(k))).(l)[m].best
+                           endif
+                        endif
+                     endfor
+                  endif
+               endfor
+            endif            
+         endif else if n_tags(ss.(i)[j].(k)) ne 0 then begin
+            if n_tags(ss.(i)[j].(k)) ne 0 then begin
+               if tag_exist(ss.(i)[j].(k),'derive') then begin
+                  if ss.(i)[j].(k).derive then begin
+                     pars = (reform(ss.(i)[j].(k).value,nsteps,ss.nchains))
+                     derive = 1B
+                     label = ss.(i)[j].(k).label 
+                     unit = ss.(i)[j].(k).unit
+                     latex = ss.(i)[j].(k).latex
+                     best = ss.(i)[j].(k).best
+                  endif
+               endif
+            endif
+         endif
+
+         if derive then begin
+
+            if n_elements(allpars) eq 0 then allpars = transpose(pars) $
+            else allpars = [allpars,transpose(pars)]
+
+            ;; check for bad values
+            bad = where(~finite(pars),complement=good)
+            if bad[0] ne -1 then message, $
+               "ERROR: NaNs in " + label + " distribution"
+
+            ;; if angular, center distribution about the mode
+            if unit eq 'DEGREES' then halfrange = 180d0 $
+            else if unit eq 'RADIANS' then halfrange = !dpi
+            if unit eq 'DEGREES' or unit eq 'RADIANS' then begin
+
+               ;; find the mode
+               hist = histogram(pars,nbins=100,locations=x,/nan)
+               max = max(hist,modendx)
+               mode = x[modendx]
+
+               toohigh = where(pars gt (mode + halfrange))
+               if toohigh[0] ne -1 then pars[toohigh] -= 2.d0*halfrange
+
+               toolow = where(pars lt (mode - halfrange))
+               if toolow[0] ne -1 then pars[toolow] += 2.d0*halfrange
+            endif                     
+
+            ;; 68% confidence interval
+            sorted = sort(pars)
+            medvalue = pars[sorted[medndx]]
+            upper = pars[sorted[hisigndx]] - medvalue
+            lower = medvalue - pars[sorted[lowsigndx]]
+
+            if n_elements(medianpars) eq 0 then medianpars = [medvalue,upper,lower] $
+            else medianpars = [[medianpars],[[medvalue,upper,lower]]]
+
+            xmax = (medvalue + 4*upper) < max(pars)
+            xmin = (medvalue + 4*lower) > min(pars)
+            parnames = [parnames,latex]
+
+            if xmin eq xmax then begin
+               message, 'WARNING: ' + label + ' is singularly valued.',/continue
+            endif else begin
+
+               ;; plot labels
+               xtitle='!3 Chain link number'
+               ytitle='!3' + textoidl(latex + '_{' + ss.(i)[j].label + '}')
+
+               plot, [0],[0], psym=10, xtitle=xtitle, ytitle=ytitle,$
+                     charsize=charsize,xstyle=1, yrange=[xmin,xmax], xrange=[0,nsteps],font=1
+
+               for l=0L, ss.nchains-1L do $
+                  oplot, pars[*,l], color=l*255d0/ss.nchains ;, transparency=100d0-100d0/ss.nchains
+               oplot, [ss.burnndx,ss.burnndx],[-9d99,9d99]
+
+            endelse
+         endif
+      endfor
+   endfor
+endfor
+
+end

gettc.pro

@@ -0,0 +1,108 @@
+pro gettc, secondary=secondary
+
+npoints = 1d5
+
+G = 2942.71377d0 ;; R_sun^3/(m_sun*day^2), Torres 2010
+nplanets = 3d0
+period = [1d0,2d0,10d0]
+period = 1+randomu(seed,3)*30d0
+period = period[sort(period)]
+tc = [0d0,3d0,5d0]
+tc = randomu(seed,3)*max(period)
+inc=[!dpi/2d0,!dpi/2d0,!dpi/2d0]
+m1 = 1d0
+m2 = [0.001d0,0.01d0,0.01d0]
+q = m1/m2
+a = (period^2*G*(m1+m2)/(4d0*!dpi^2))^(1d0/3d0)
+;e = [0.5d0,0.25d0,0.9d0]
+e = dblarr(3);[0.5d0,0.25d0,0.9d0]
+omega = [!dpi/5d0,!dpi/6d0,!dpi/7d0]
+
+
+bjd = min(tc) + dindgen(npoints)/(npoints-1)*max(period)*10
+tp = dblarr(nplanets)
+for i=0, nplanets-1L do begin
+   phase = exofast_getphase(e[i],omega[i], /primary)
+   tp[i] = tc[i] - period[i]*phase
+endfor
+
+if 0 then begin
+b = exofast_getb2(bjd, inc=inc,a=a,tperiastron=tp, period=period,e=e,omega=omega,q=q,x0=x,y0=y,z0=z,x1=x1,y1=y1,z1=z1,x2=x2,y2=y2,z2=z2)
+
+for i=0, nplanets-1 do begin
+
+   ndx = lclxtrem(b[i,*])
+!p.multi=[0,2,2]
+   plot, b[i,*]
+   oplot, b[i,ndx], b[i,ndx], psym=1, color='0000ff'x
+
+endfor
+endif
+
+a2 = a*q/(1d0+q)
+a1 = a2/q  
+tol = 1d-8 ;; 1 ms
+ntransits = 10
+
+for i=0L, nplanets-1 do begin
+
+   print
+   print, 'planet ' + strtrim(i,2)
+
+   for k=0, ntransits-1 do begin
+
+      t0 = systime(/seconds) 
+      minbjd = tc[i]-period[i]/4d0+k*period[i]
+      maxbjd = tc[i]+period[i]/4d0+k*period[i]
+      repeat begin
+
+         bjd = (minbjd + maxbjd)/2d0
+
+         x1 = 0d0
+         z1 = 0d0
+         for j=0L, nplanets-1 do begin
+            meananom = (2.d0*!dpi*(1.d0 + (bjd - tp[j])/period[j])) mod (2.d0*!dpi)
+            eccanom = exofast_keplereq(meananom, e[j])
+            trueanom = 2d0*atan(sqrt((1d0 + e[j])/(1d0 - e[j]))*tan(0.5d0*eccanom))
+
+            ;; calculate the star position in the barycentric frame
+            r1 = a1[j]*(1d0-e[j]^2)/(1d0+e[j]*cos(trueanom))
+
+            ;; rotate to observer's plane of reference
+            x1 += -r1*cos(trueanom + omega[j] + !dpi)
+            z1 += r1*sin(trueanom + omega[j] + !dpi)*sin(inc[j])
+         endfor
+
+         meananom = (2.d0*!dpi*(1.d0 + (bjd - tp[i])/period[i])) mod (2.d0*!dpi)
+         eccanom = exofast_keplereq(meananom, e[i])
+         trueanom = 2d0*atan(sqrt((1d0 + e[i])/(1d0 - e[i]))*tan(0.5d0*eccanom))
+
+         r = a2[i]*(1d0-e[i]^2)/(1d0+e[i]*cos(trueanom))
+         x = -r*cos(trueanom + omega[i])-x1
+         z = r*sin(trueanom + omega[i] + !dpi)*sin(inc[i]) -z1
+
+         ;; bracket the primary or secondary
+         if keyword_set(secondary) then begin
+            if z le 0d0 and x gt 0d0 then minbjd = bjd $
+            else if z le 0d0 and x le 0d0  then maxbjd = bjd $
+            else if z gt 0d0 and x gt 0d0 then maxbjd = bjd $
+            else if z gt 0d0 and x le 0d0 then minbjd = bjd 
+         endif else begin
+            if z le 0d0 and x gt 0d0 then maxbjd = bjd $
+            else if z le 0d0 and x le 0d0  then minbjd = bjd $
+            else if z gt 0d0 and x gt 0d0 then minbjd = bjd $
+            else if z gt 0d0 and x le 0d0 then maxbjd = bjd 
+         endelse
+
+         if x gt 0 then begin
+            maxbjd = bjd
+         endif else minbjd = bjd
+
+      endrep until abs(x) le tol
+
+      print, bjd, (bjd-tc[i]-k*period[i])*86400d0, systime(/seconds)-t0
+   endfor
+endfor
+
+
+end

gridpars.pro

@@ -0,0 +1,60 @@
+pro gridpars
+
+files = file_search('sed/nextgenfin/nextgenfin2/*.idl',count=nfiles)
+
+
+teff = []
+logg = []
+feh = []
+alpha = []
+
+for i=0L, nfiles-1 do begin
+
+   filename = file_basename(files[i],'.NextGen.spec.idl')
+   len = strlen(filename)
+   teff = [teff,double(strmid(filename,3,len-15))*100d0]
+   logg = [logg,double(strmid(filename,len-12,4))]
+   feh = [feh,double(strmid(filename,len-8,4))]
+   alpha = [alpha,double(strmid(filename,len-4,4))]
+
+endfor
+
+tefforig = teff
+loggorig = logg
+fehorig = feh
+alphaorig = alpha
+
+teff = teff[sort(teff)]
+teffs = teff[uniq(teff)]
+
+logg = logg[sort(logg)]
+loggs = logg[uniq(logg)]
+
+feh = feh[sort(feh)]
+fehs = feh[uniq(feh)]
+
+alpha = alpha[sort(alpha)]
+alphas = alpha[uniq(alpha)]
+
+nteffs = n_elements(teffs)
+nloggs = n_elements(loggs)
+nfehs = n_elements(fehs)
+nalphas = n_elements(alphas)
+grid = bytarr(nteffs,nloggs,nfehs,nalphas)
+for i=0,nteffs-1L do begin
+   for j=0, nloggs-1L do begin
+      for k=0, nfehs-1L do begin
+         for l=0, nalphas-1L do begin
+            match = where(tefforig eq teffs[i] and loggorig eq loggs[j] and $
+                          fehorig eq fehs[k] and alphaorig eq alphas[l])
+            if match[0] ne -1 then begin
+               grid[i,j,k,l] = 1B
+               oplot, teffs[i],loggs[j]
+            endif
+         endfor
+      endfor
+   endfor
+endfor
+
+stop
+end

mutual_eclipse.pro

@@ -0,0 +1,23 @@
+;; (x,y,z) = (0,0,0) is the center of the star
+;; x -- An NPlanets x NTimes array of X positions of each planet, in stellar radii
+;; y -- An NPlanets x NTimes array of Y positions of each planet, in stellar radii
+;; z -- An NPlanets x NTimes array of Z positions of each planet, in stellar radii
+;; p -- An NPlanets array of the radius of each planet, in stellar radii
+;; u1 -- linear limb darkening parameter
+;; u2 -- quadratic limb darkening parameter
+
+function mutual_eclipse, x, y, z, p, u1, u2
+
+sz = size(x)
+nplanets = sz[1]
+ntimes = sz[2]
+flux = dblarr(ntimes) + 1d0
+
+for i=0, nplanets-1 do begin
+
+
+endfor
+
+return, flux
+
+end