Old Kepler May 29 2012

Kepler.pro

pro readdata, max
 
readcol, 'data.txt', name, KIC, Kp, t_0, t_0_unc, P, P_unc, Rp, a, T_eq, Dur, Depth, a_to_Rstar, a_to_Rstar_unc, r_to_Rstar, r_to_Rstar_unc, b, b_unc, SNR, chi, T_eff, log_g, Rad2, f_T_eff, /silent ;read in the given file with all the information
  
BJD=2456069                     ; May 21th 2012, initial start time for the plot
t=findgen(1000000)/10000+BJD          ;ends August 29 2012
dist=500                        ;pc distance from earth to the objects
f=4.84e-6                       ;convert AU to pc
ff=2.254e-8                     ;conert solar R to pc
 
for i=0, max-1 do begin         ; # of entries from 0 to 2322th all the entries in the given table
    array = [0]
    
    while round(name[i]) eq round(name[(i+1) mod max]) do begin ;omit the numbers after decimal has the same so that I can catagorize them into different systems 
        if array[n_elements(array)-1] ne i then begin ;the first pair in each sys
            array = [array, i, i+1] ;first pair, to count how many planets are in that system
        endif else begin
            array = [array, i+1]
        endelse
        i = i+1
    endwhile
    
    if n_elements(array) gt 1 then begin ;more than one planet
        array = array[1:*]      ;crop array 
        number=n_elements(array) ;# of planets in systsm
                                ; if n_elements(array) gt 1 then begin
        gamarr=fltarr(n_elements(t)) ;define as an array with same dim as t
        thetarr=fltarr(n_elements(t))
        thetarr_l=fltarr(n_elements(t))
        for j = 0, n_elements(array)-1 do begin ;in the system
            gamma=a[array[j]]*f/dist*sin(2*!pi*(t-t_0[array[j]])/P[array[j]]) ;angle seen from Earth, array[j] gives that # entry
            gamarr = [[gamarr], [gamma]]
            theta= r_to_Rstar[array[j]]*(Rad2)[array[j]]*ff/dist ;angular size of each object
            theta= fltarr(n_elements(t))+theta ; make theta has the same dim as t
            thetarr = [[thetarr],[theta]] ;same dim as gamarr
     
    
        endfor
 
 
;make the plot for each combination in a system, e.g. 1a with 1b, 1b
;with 1c and 1a with 1c
        for k = 1, (size(gamarr))[2]-2 do begin
            for l = k+1, (size(gamarr))[2]-1 do begin
                psopen, 'plots/blah'+prettynum(name[array[k-1]]) + '_' + prettynum(name[array[l-1]]) +".ps", /inch, xsize=6, ysize=6 ;prettynum is a program
 
                plot, t, alog(abs(gamarr[*,k]-gamarr[*,l])), xtitle="date [BJD]", ytitle="log difference in angle from E [rad]", title="log Theta VS time"
                item=['The two planets being plotted are ']+prettynum(name[array[k-1]])+[' and ']+prettynum(name[array[l-1]])
                legend, item, charsize =1.1
                item2=['log theta=log delta gamma','log theta=log sum of angular size']
                legend, item2, lines=[0,5], /right, /bottom
                oplot, t, alog(abs(thetarr[*,k]-thetarr[*,l])), linestyle=5 ;the critical condition, such that below this line is the time worth observing
 
                psclose
                
                theta_l=alog(abs(gamarr[*,k]-gamarr[*,l])) ; angular separation between planets
                thetarr_l=[[thetarr_l],[theta_l]] ;angular separation between planets store as array later used for n body
           
       
 
       ; thetarr_l=thetarr_l[*,1:number] ;dim=n_elements(t)Xnumber of planets
       ; for k = 1, (size(gamarr))[2]-2 do begin
       ;     for l = k+1, (size(gamarr))[2]-1 do begi
                
                tarr=[0]
                starr=[0]
                etarr=[0]
                
                                ;pairs alignment              
                if i eq 4 then begin 
                    
                    openw, unit, 'transitime.txt', /get_lun
                    
                    for t_ind=0, 100, 1./10000 do begin
                                ;while y lt (n_elements(t)+BJD) do begin
                        if (abs(gamarr[*,k]-gamarr[*,l]))[t_ind] EQ (abs(thetarr[*,k]-thetarr[*,l]))[t_ind] then begin ;both array have dim of n_elements X 1 
                            
                            tarr=[tarr,t[t_ind]]
                        endif
                    endfor
                    tarr=tarr[1:*]
                                ;separate start time and endtime
                    starr=tarr[0:*:2]
                    etarr=tarr[1:*:2]
                    
                    namename=string(prettynum(name[array[k-1]])+'_'+prettynum(name[array[l-1]]))
                    aa=string([namename,tarr])
                    printf, unit, aa
                    close, unit
                endif else begin
                                ;to save the same thing for other pairs
                    openu, unit, 'transitime.txt', /get_lun, /append
                    for t_ind=0, 100, 1./10000 do begin
                                ;while y lt (n_elements(t)+BJD) do begin
                        if (abs(gamarr[*,k]-gamarr[*,l]))[t_ind] EQ (abs(thetarr[*,k]-thetarr[*,l]))[t_ind] then begin ;both array have dim of n_elements X 1 
                            
                            tarr=[tarr,t[t_ind]]
                        endif
                    endfor
                    tarr=tarr[1:*]
                    namename=string(prettynum(name[array[k-1]])+'_'+prettynum(name[array[l-1]]))
                    aa=string([namename,tarr])
                    printf, unit, aa
                    close, unit
                endelse
                
            endfor
        endfor
        
        thetarr_l=thetarr_l[*,1:number] ;dim=n_elements(t)Xnumber of planet
        ttarr=[0]
;find the same thing for n-bodies, e.g 1a, 1b and 1c all aligned
        for j=0, number-1 do begin ;number was defined as the total number of planets in that system
            openu, unit, 'transitime.txt', /get_lun, /append
            for t_ind=0, 100, 1./10000 do begin
                sumtheta=0
                for m=0, number-1 do begin
                    sumtheta=sumtheta+thetarr[t_ind,m]
                endfor
            endfor
            for t_ind=0, 100, 1./10000 do begin
                if max(thetarr_l[t_ind,*]) eq sumtheta then begin
                    ttarr=[ttarr,t[t_ind]]	
                endif
            endfor
            ttarr=ttarr[1:*]
            namename=0
            for m=0, number-1 do begin
                namename=namename+prettynum(name[array[m]])
                namename=string(namename)
            endfor                   
            aa=string([namename,ttarr])
            printf, unit, aa
            close, unit
        endfor
    endif
    
endfor
end