;MPRVFIT.PRO created by Nathan De Lee 09/16/2008 ;This procedure is designed to fit Keplerian orbits to radial velocity ;data. It does so by running a periodogram and fitting Keplerian ;models to the highest periodogram peaks roughly following the work of ;Cumming, A. 2004, MNRAS, 354, 1165. ; ; Purpose: Fitting keplerian orbits to radial velocity data ; ; Usage: MPRVFIT,INPUT.CUR,BASENM=BASENM,CSIZE=CSIZE,GRID=[FSTEP,DIV],$ ; HIFAC=HIFAC,MANUAL=PERIOD,NPERIODS=NUM,OFFSET=[OFF1,OFF2],$ ; OVER=OVER,SB2=['IN1',IN2',etc],SLOPE=SLOPE,WIND=[PMIN,PMAX],$ ; YRANGE=[MIN,MAX],YRRANGE=[MIN,MAX],/BUILTIN,/CIRC,/COLOR,/FAP,$ ; /FAST,/FORCE,/KEPLER,/LINEAR,/LOUD,/MASSPL,/MASSTI,/NOFIT,$ ; /NOOFF,/NOTITLE,/PHASE,/PLAIN,/PLINE,/PLOG,/POST,/RESID,/RVSIM,$ ; /SCREEN,/SHOWECC,/UNFOLD,/VERB ; ; Inputs: ; INPUT.CUR -- Input RVcurve file: JD Vel Err (May be an ; array of input files) ; or ; ['INPUT1.cur','INPUT2.cur',etc.] ; ; ; Outputs: ; OUTPUT.LOG -- A Log file with the fit parameters ; ; Optional: ; BASENM -- Base file name for plots etc. Default (Input name). ; BUILTIN -- Use the built in lomb-scargle periodogram ; CIRC -- Force a circular fit ; COLOR -- Make the ps plot in color ; CSIZE -- Character size, Default 1.0 ; FAP -- Create a periodogram file: period fap ; FAST -- Use a fast but less precise FAP method ; GRID -- Do a grid search in period space around ; periodogram peaks. FSTEP is the number of ; frequency steps to go away from the peak and ; DIV is how many intra FSTEP steps to make. ; HIFAC -- Set the HIFAC variable used in periodogram search. ; KEPLER -- Do the Keplerian orbit based periodogram ; LINEAR -- Include a linear component in the Periodogram and ; Keplerian orbit fit. ; LOUD -- Output the interations in non-linear fitting ; MANUAL -- Provide a period, do not search for one. ; MASSPL -- Create a plot of companion mass versus ; host star mass. ; MASSTI -- Display mass estimate in title. ; NOFIT -- Only run periodogram ; NOOFF -- Do not fit offsets. ; NOTITLE -- Removes title from plot ; NPERIODS -- Number of periods to search (Max 15, Default 2). ; OFFSET -- Apply this offset to the RVs (May be a list ; if multiple input files). ; OVER -- Over sampling compared to Nyquist (Default 5). ; PHASE -- Allows this procedure to read in a phased ; RVcurve.pha: JD Vel Err Phase ; PLAIN -- Removes fit numbers from plot ; PLINE -- Plot lines for two highest peak periods ; PLOG -- Make period axis of periodogram log ; POST -- Produce a postscript output ; RESID -- Create an O-C plot ; RVSIM -- Accept RVSIM style input files and print a ; visibility plot. Input: JD RV ERR VIS SPEC ; SB2 -- An input array similar to input, that gives ; the rvcurves for the second component. ; SCREEN -- Make a plot on the screen ; SHOWECC -- Shows ecc and argp instead of ecosw esinw ; SLOPE -- Give an array of slopes equal to nperiods to ; use as a starting point for the search. ; UNFOLD -- Create unfolded RV plots ; VERBOSE -- Print verbose output ; WINDOW -- The period searching window in days, Default ; is based on the observations. If two ; elements, then [per_low,per_hi] if 3 elements ; then [per_low,per_hi,freq_step] ; YRANGE -- Specify y range for unfolded plots, but will ; do folded plot correctly as well. ; YRRANGE -- Specify y range for residual plots. ; ; Requires: ; MPFIT.PRO URL: http://www.physics.wisc.edu/~craigm/idl/fitting.html ; ASTROLIB URL: http://idlastro.gsfc.nasa.gov/ ; FIXPS.PRO UFL: http://www.idlcoyote.com/documents/programs.php ; ; ; Example: ; Run Lomb-Scargle periodogram and fit Keplerian models to the ; two highest periodogram peaks. Output a phase folded color postscript. ; IDL> mprvfit,'input.cur',/POST,/COLOR,/PLINE,/VERBOSE ; ; Revision History: ; Version 1.0 released to public on April 4th, 2013 by Nathan De Lee ;Inputs: T transit ;Outputs: T periastron function convert_t2p,date,period,arg_p,ecc ;Using equations from Kane & von Braun (2009) ;E = 2*arctan ( sqrt(1-e)/sqrt(1+e) *tan((pi/2 - argp) /2) ;At epoch of transit arg_p + f = pi/2 eanomoly = sqrt(1-ecc)/sqrt(1+ecc) * tan((!DPI/2 - arg_p)/2) eanomoly = 2 * atan(eanomoly) manomoly = eanomoly - ecc*sin(eanomoly) correct = period * manomoly / (2 * !DPI) date = date - correct return,date end ;Inputs: JD, Epoch of Periastron, Period, Eccentricity ;Ouputs: eccentric anomaly function eanomaly,data,ep_peri,period,ecc ;The eccentric anomaly is related to time through the mean anomaly. These ;relations are M = 2*pi/period *(t-t_0) and M = E - e * sin E. ;This cannot be solved analytically, so the Newton-Raphson iteration method ;must be used. The formulation for this equation is then: ;E_(n+1) = M - e(E_n*cos E_n - sin E_n) / (1 - e*cos E_n) ;Charles & Tatum 1998 suggest that E_0 should equal pi ;Convergence criteria (same as for RVSIM) conv = 1.0d-4 ;m must go between 0 and 2pi m = (data - ep_peri)/period m = 2D*!DPI * (m - floor(m)) eold = dblarr(n_elements(m))+!DPI enew = dblarr(n_elements(m)) ;Run the Newton-Raphson method for i=0L,n_elements(m)-1 do begin n=0 while(1) do begin enew[i] = (m[i] - ecc * (eold[i] * cos(eold[i]) - sin(eold[i]))) $ /(1-ecc*cos(eold[i])) ;Added in the $eold < $conv to avoid $eold of 0 issues. if(eold[i] lt conv) then break if(abs((enew[i] - eold[i])/eold[i]) lt conv) then $ break if(n > 50) then begin print, 'Too many iterations!' break endif eold[i] = enew[i] n++ endwhile endfor return,enew end function finalfit, $ period, $ date, $ rv, $ err, $ amp, $ v0, $ slope, $ offset, $ loud=loud, $ linear=linear, $ circ=circ, $ manual=manual, $ dstar=dstar, $ nooff=nooff, $ is_numdatapoints_insufficent=is_numdatapoints_insufficent nparam = 8 + n_elements(offset) ;Set noisiness of MP_MPFIT if(keyword_set(loud)) then $ soft = 0 $ else $ soft = 1 ;Setup parameters for fit parinfo = replicate({value:0.D, fixed:0, limited:[0,0], $ limits:[0.D,0]}, nparam) wfreq = (2.0D * !DPI / period) ;T0 should be after minimum date and within a period. Tmid = min(date) + period/2.0D parinfo[0:7].value = [wfreq,amp,0.01D,0.01D,Tmid,v0,slope,0.0D] ;Fill in the values for the other offsets parinfo[8:n_elements(offset)+7].value = offset ;frequency is fixed in manual mode if(n_elements(manual) gt 0) then $ parinfo[0].fixed = 1 ;Offset of data set 1 is fixed parinfo[8].fixed = 1 ;ecosw and esinw go between -1 and 1 and K should be greater than 0 ;T0 should be after minimum date and within a period. parinfo[0].limited = [1,1] parinfo[1].limited[0] = 1 parinfo[7].limited[0] = 1 parinfo[2].limited = [1,1] parinfo[3].limited = [1,1] parinfo[4].limited = [1,1] parinfo[0].limits = [.20D*!DPI/period,20.0D * !DPI/period] parinfo[1].limits[0] = 2*amp/3 parinfo[2].limits = [-1.0D,1.0D] parinfo[3].limits = [-1.0D,1.0D] parinfo[4].limits = [min(date),min(date)+ period] ;If it was not a linear fit, then free the slope to be fit, else keep ;it fixed. if(keyword_set(linear)) then $ parinfo[6].fixed = 0 $ else $ parinfo[6].fixed = 1 ;If circular set then ecosw and esinw are both 0 and fixed if(keyword_set(circ)) then begin parinfo[2].fixed = 1 parinfo[3].fixed = 1 parinfo[2].value = 0.0 parinfo[3].value = 0.0 endif ;Set semi2 to 0 unless sb2 case sb2_elem = where(dstar eq 1, scount) if(scount gt 0) then begin parinfo[7].value = (max(rv[sb2_elem]) - min (rv[sb2_elem]))/3 endif else begin parinfo[7].fixed = 1 parinfo[7].value = 0.0 endelse ;If nooff then don't fit offsets if(keyword_set(nooff)) then $ parinfo[8:n_elements(offset)+7].fixed = replicate(1,n_elements(offset)) ;Setup structure to take x,y,and err functargs = {date:date, rv:rv, err:err} ;Setup structure to store results fit = {params:dblarr(nparam),chi:0.D,cov:dblarr(nparam,nparam),dof:0.D,period:0.D,tperi:0.D,perror:dblarr(nparam)} ;;============================================================================== ;; 1. If we have enough data points, do the fit. ;;============================================================================== if ~keyword_set(is_numdatapoints_insufficent) then begin fit.params = MP_MPFIT('mporbit', PARINFO=parinfo,FUNCTARGS=functargs,BESTNORM=chi,DOF=dof,COV=cov,PERROR=perror,MAXITER=200,QUIET=soft,ERRMSG=errmsg) endif else begin chi=!values.d_nan dof=!values.d_nan endelse ;Print errors if necessary if(keyword_set(loud) and keyword_set(errmsg)) then $ print,errmsg ;Store parameter errors if(keyword_set(perror)) then begin fit.perror = perror fit.cov = cov endif ;Store chi2 info fit.chi = chi fit.dof = dof ;Store Period fit.period = (2.0D * !DPI / fit.params[0]) ;Store T_peri ;Get eccentricity and argp ecc = sqrt(fit.params[2]^2 + fit.params[3]^2) if(ecc eq 0.0) then $ argp = 0.0 $ else $ argp = atan(fit.params[3],fit.params[2]) fit.tperi = convert_t2p(fit.params[4],fit.period,argp,ecc) ;Send fit values back return,fit end ;This function finds all the peaks in the periodogram by searching for ;local maxima. ;Inputs: Lomb-Scargle Power zls ; Number of periods to return ; ;Outputs: Array of indices for zls function peakfinder,zls,nperiods first = 1 uphill = 0 for i = 0L,n_elements(zls)-1 do begin ;Initialize things if(first) then begin oldzls = zls[i] first = 0 endif ;Are we going up hill? if(zls[i] gt oldzls) then $ uphill = 1 ;Have we creasted the hill if(zls[i] lt oldzls and uphill) then begin ;Save the index of the peak if(n_elements(peakindice) eq 0) then $ peakindice = i-1 $ else $ peakindice = [[peakindice],[i-1]] ;Going down hill uphill = 0 endif oldzls = zls[i] endfor ;Find nperiods highest peaks zls_power = zls[peakindice] best_indice = peakindice[(reverse(sort(zls_power)))[0:nperiods-1]] return,best_indice end ;This function returns the False Alarm Probability FAP based on the ;work of Baluev (2008) for Cummings (2004) style (2) periodograms. It ;works for both with and without linear trends assuming the power ;calculation is done correctly. ;Inputs: Lomb-Scargle Power zls ; Dates of observation dates ; Angular frequencies searched wfreq ; Highest searched frequency ; Degrees of Freedom in sine fit nu ; Big array, which is a representaiton of the periodogram ;Outputs: FAP or if no big provided returns 2 element big array function get_cfap,zls,date,wfreq,f_hi,nu,linear=linear,fast=fast,big=big ;Gaussian probability (don't do math in exponents) pexp = -nu/2d0 prob = (1d0 + (2d0*zls/nu))^(pexp) if(keyword_set(fast)) then begin ;Determine the 50% and 99% False Alarm Probability (FAP) ;I use M defined as 2*Highest_frequency * T(obervation period) ;Cumming 2004 num_freq = 2*f_hi*(max(date)-min(date)) ;Determine which method to use to calculate FAP ;Percentiles above 99% should use - ln(F/M) ;Percentiles below 99% should use - ln(1-(1-F)^(1/M)) ;Deal with insanely high powers if(zls gt 50) then $ best_fap = num_freq*2E-22 $ else $ best_fap = num_freq * prob if(best_fap gt 0.01 ) then $ best_fap = 1.0D - ((1.0D - prob)^num_freq) return,best_fap endif ; Next 5 lines are what i added in Oct 19th,2009 ; We should calculate A from the arrays to avoid not alias free situation ; EIGENVEC is a IDL function to calculate eigenvalues tobs2 = date date=date-min(date) meant = total(date)/n_elements(date) meant2 = total(date*date)/n_elements(date) ;Only run this section once if (n_elements(big) ne 2) then begin bigA = 0d0 for nfre = 0d0, n_elements(wfreq)-1d0,1d0 do begin ; The next are to calculate a and b for a floating point situation arrayX = MAKE_ARRAY(1,n_elements(date), /double, VALUE = 1d0) arrayW = [TRANSPOSE(cos(wfreq[nfre]*date)),TRANSPOSE(sin(wfreq[nfre]*date))] arrayXandW = [arrayX,arrayW] arrayWtheta = [TRANSPOSE(-2*!pi*date*sin(wfreq[nfre]*date)),TRANSPOSE(2*!pi*date*cos(wfreq[nfre]*date))] arrayL = IDENTITY(n_elements(date))-arrayXandW##INVERT(TRANSPOSE(arrayXandW)##arrayXandW)##TRANSPOSE(arrayXandW) arrayMD = INVERT(TRANSPOSE(arrayXandW)##arrayXandW) arrayMD = arrayMD[(n_elements(arrayMD)^0.5-2):(n_elements(arrayMD)^0.5-1),(n_elements(arrayMD)^0.5-2):(n_elements(arrayMD)^0.5-1)] arrayMD2 = arrayL##arrayWtheta##arrayMD##TRANSPOSE(arrayWtheta)##arrayL eval3 = LA_HQR(LA_ELMHES(arrayMD2,/double), /DOUBLE) evalp = eval3[where(abs(real_PART(eval3)) gt 1d-6)] if(n_elements(evalp) ne 2 ) then begin print,'stop because of eval3 = ', eval3 stop endif a = (REAL_PART(evalp[0]))^(-0.5d0) b = (REAL_PART(evalp[1]))^(-0.5d0) ; The next are to calculate ap and bp for the linear trend situation arrayX = [arrayX,TRANSPOSE(date)] arrayW = [TRANSPOSE(cos(wfreq[nfre]*date)),TRANSPOSE(sin(wfreq[nfre]*date))] arrayXandW = [arrayX,arrayW] arrayWtheta = [TRANSPOSE(-2*!pi*date*sin(wfreq[nfre]*date)),TRANSPOSE(2*!pi*date*cos(wfreq[nfre]*date))] arrayL = IDENTITY(n_elements(date))-arrayXandW##INVERT(TRANSPOSE(arrayXandW)##arrayXandW)##TRANSPOSE(arrayXandW) arrayMD = INVERT(TRANSPOSE(arrayXandW)##arrayXandW) arrayMD = arrayMD[(n_elements(arrayMD)^0.5-2):(n_elements(arrayMD)^0.5-1),(n_elements(arrayMD)^0.5-2):(n_elements(arrayMD)^0.5-1)] arrayMD2 = arrayL##arrayWtheta##arrayMD##TRANSPOSE(arrayWtheta)##arrayL eval3 = LA_HQR(LA_ELMHES(arrayMD2,/double), /DOUBLE) evalp = eval3[where(abs(real_PART(eval3)) gt 1d-6)] if(n_elements(evalp) ne 2 ) then begin print,'stop because of eval3 = ', eval3 stop endif ap = (REAL_PART(evalp[0]))^(-0.5d0) bp = (REAL_PART(evalp[1]))^(-0.5d0) pied = 4d0*(!pi*a*b+(a-b)^2d0)/(a+b)-0.5d0*a*b/(a+b)*((a-b)^2d0/((a+b)^2d0+!pi*a*b)) piedp = 4d0*(!pi*ap*bp+(ap-bp)^2d0)/(ap+bp)-0.5d0*ap*bp/(ap+bp)*((ap-bp)^2d0/((ap+bp)^2d0+!pi*ap*bp)) if (nfre eq 0) then begin bigA = 0d0 +wfreq[nfre]/(2d0*!pi)*pied/(a*b) bigAp= 0d0 +wfreq[nfre]/(2d0*!pi)*piedp/(ap*bp) endif else begin bigA = bigA +(wfreq[nfre]-wfreq[nfre-1])/(2d0*!pi)*pied/(a*b) bigAp = bigAp +(wfreq[nfre]-wfreq[nfre-1])/(2d0*!pi)*piedp/(ap*bp) endelse endfor date=tobs2 return,[bigA,bigAp] endif else begin bigA = big[0] bigAp = big[1] endelse ; print,'bigA from integration = ',bigA,a,b,pied ; print,'bigA from modified = ',bigAp,ap,bp,piedp date=tobs2 ;A is an integration ;Do the floating point case if(not keyword_set(linear)) then begin fapsingle = prob ; bigw = f_hi*(4d0*!pi*(meant2-meant^2d0))^0.5d0 ; print,'bigA in alias free situation = ',2*!pi^1.5d0*bigw ; taufromb = (2d0/nu)^0.5*gamma(nu/2d0+1d0)/gamma(nu/2d0+0.5d0)*fapsingle*(zls)^(0.5d0)*bigw ; best_fap[i]=1d0-1d0*exp(-taufromb)+fapsingle*exp(-taufromb) taufromb2 = gamma(nu/2d0+1d0)/gamma(nu/2d0+0.5d0)/(2d0*!pi)*bigA*fapsingle*(2d0*zls/!pi/nu)^0.5d fap_22=1d0-1d0*exp(-taufromb2)+fapsingle*exp(-taufromb2) fap = fap_22 endif if(keyword_set(linear)) then begin ;fap_33 is a method for the situation you have a liner trend fapsingle2 = (1d0 + (2d0*zls/nu))^(-nu/2d0) taufromb3 = gamma(nu/2d0+1d0)/gamma(nu/2d0+0.5d0)/(2d0*!pi)*bigAp*fapsingle2*(2d0*zls/!pi/nu)^0.5d fap_33=1d0-1d0*exp(-taufromb3)+fapsingle2*exp(-taufromb3) fap = fap_33 endif ;; print,'fap single : ',fapsingle,(1d0-fapsingle) ;*exp(-taufromb3) ;; print,'best period :',best_period ;; print,'fap calculated for linear trend situation: ',fap_33,taufromb3 ;; print,'fap calculated for floating point situation : ', fap_22,taufromb2 ;; print,'fap calculated from alias free approximation situation: ',best_fap[i],taufromb return,fap end ;This function returns the Periodogram Power (zls) for the input FAP ;it uses a simple search algorithm coupled with the get_cfap function ;Inputs: FAP that you want the ZLS for ; Dates of observation dates ; Angular frequencies searched wfreq ; Highest searched frequency ; Degrees of Freedom in sine fit nu ; Low z to start search ;Outputs: 3 element array ZLS FAP iter function get_czls,fap,date,wfreq,f_hi,nu,zf_low,linear=linear,fast=fast,big=big ;Search tolerance factor zls_tol_factor = 0.02 zls_tol = zls_tol_factor * fap ;Start the search box limits zf_hi = 100000 fap_low = get_cfap(zf_low,date,wfreq,f_hi,nu,linear=linear,fast=fast,big=big) fap_hi = get_cfap(zf_hi,date,wfreq,f_hi,nu,linear=linear,fast=fast,big=big) ;Search through each of the ZLS i=0 startc = systime(/seconds) while(zf_hi - zf_low gt zls_tol) do begin zf_cur =(zf_hi+zf_low) / 2 fap_cur = get_cfap(zf_cur,date,wfreq,f_hi,nu,linear=linear,fast=fast,big=big) if(fap_cur lt fap) then $ zf_hi = zf_cur $ else $ zf_low = zf_cur if(abs(fap - fap_cur) lt zls_tol) then begin return,[zf_cur,fap_cur,i] endif i = i + 1 ;Useful for speed determination ; print,i,systime(/seconds) - startc,format='("Iter: ",I3,1X,F8.3,"sec")' endwhile end ;Inputs: period_array,semi_array,ecc_arr,file_basename ;Outputs: None function massplot, $ period_arr, $ semi_arr,$ ecc_arr,$ base, $ is_numdatapoints_insufficent=is_numdatapoints_insufficent, $ notitle=notitle, $ color=color entry_device = !d.name set_plot,'PS' if(keyword_set(color)) then $ device,filename=base+'_m.ps',/landscape,/color $ else $ device,filename=base+'_m.ps',/landscape,color=0 ;Define plot symobol circles use psym 8 to make this work plotsym, 0, 1, /fill ;Make thing fancy for plotting, but remember what the old values were remember = [!P.FONT,!P.THICK,!P.CHARTHICK,!X.THICK,!Y.THICK,!P.CHARSIZE] ;Pick a fancy font !P.FONT=2 ;Darken up all the lines !P.THICK=4 !P.CHARTHICK=4 !P.CHARSIZE=1.35 !X.THICK=4 !Y.THICK=4 ;Make the actual plot ;Turn off titles if(keyword_set(notitle)) then $ mytitle = '' $ else $ mytitle = base ;Make host mass array (in solar masses) mass = indgen(30)*0.1 + 0.1 ;;============================================================================== ;; 1. Check for insufficient data points. ;;============================================================================== if ~keyword_set(is_numdatapoints_insufficent) then begin ;;============================================================================== ;; 2. Make a mass plot for each period ;;============================================================================== for i=0L,n_elements(period_arr)-1 do begin mplan = semi_arr[i] * (mass * 1.98892e30)^(.66666)*(period_arr[i]*86400)^(.3333333) * sqrt(1-ecc^2) mplan /= (2*3.14159*6.673e-11)^(.333333) ;Make it all Jupiter Masses mplan /= 1.8987e27 ;Plot in mass planet versus host mass plot,mass,mplan,xtitle="Host Mass (in Solar Masses)", $ ytitle = "M sin i (in Jupiter Masses)",title=mytitle,/nodata,/ynozero oplot,mass,mplan,psym=0 endfor endif else begin for i=0L,n_elements(period_arr)-1 do begin plot,mass,0.0*mass,xtitle="Host Mass (in Solar Masses)", $ ytitle = "M sin i (in Jupiter Masses)",title=mytitle,/nodata,/ynozero oplot,mass,0.0*mass,psym=0 endfor endelse ;Close file, set back to x11 device,landscape=0,/close_file fixps,base+'_m.ps' ;Set things back to normal set_plot,entry_device !P.FONT=remember[0] !P.THICK=remember[1] !P.CHARTHICK=remember[2] !X.THICK=remember[3] !Y.THICK=remember[4] !P.multi=0 return,1 end function mporbit,params,date=date,rv=rv,err=err common dataset, dset, dstar ;Separate out the parinfo structure period = (2.0D * !DPI / params[0]) k = params[1] ecosw = params[2] esinw = params[3] ep_peri = params[4] v0 = params[5] slope = params[6] k2 = params[7] ;Apply the offset to the rv for i=8L,n_elements(params)-1 do begin num = i-8 rv[where(dset eq num)] = rv[where(dset eq num)] + params[i] endfor ;Split data up for two stars in sb2 case ;Determine if we are in the sb2 case sb2_elem = where(dstar eq 1, scount) date1 = date[where(dstar eq 0)] rv1 = rv[where(dstar eq 0)] err1 = err[where(dstar eq 0)] if(scount gt 0) then begin date2 = date[where(dstar eq 1)] rv2 = rv[where(dstar eq 1)] err2 = err[where(dstar eq 1)] endif ;Get eccentricity and argp ecc = sqrt(ecosw^2 + esinw^2) if(ecc eq 0.0) then $ argp = 0.0 $ else $ argp = atan(esinw,ecosw) ;I need E the eccentric anomaly. This will be returned from function eanomaly enew = eanomaly(date1,ep_peri,period,ecc) ;Get the true anomaly in both sine and cosine. These are identical equations ;just used sin^2 + cos^2 = 1 to transform cosf = (cos(enew) - ecc)/(1.D - ecc*cos(enew)) sinf = sqrt(1.D - ecc^2)*sin(enew)/(1.D - ecc*cos(enew)) ;Find the radial velocity ; vel = v0 + k*(cos(argp + f) + ecc * cos(argp)) ;Since I don't have f must use the cosine sum formula ovel1 = v0 + k*(cos(argp)*cosf -sin(argp)*sinf +ecc*cos(argp)) + slope*date1 ;Determine the standard deviant deviates1 = (rv1 - ovel1) / err1 ;If sb2 case remember that the two orbits share identical values ;except for semi-amplitude and argp (which is PI off). if(scount gt 0) then begin enew = eanomaly(date2,ep_peri,period,ecc) cosf = (cos(enew) - ecc)/(1.D - ecc*cos(enew)) sinf = sqrt(1.D - ecc^2)*sin(enew)/(1.D - ecc*cos(enew)) ovel2 = v0 + k2*(cos(argp+!DPI)*cosf -sin(argp+!DPI)*sinf +ecc*cos(argp+!DPI)) + slope*date2 deviates2 = (rv2 - ovel2) / err2 return,[deviates1,deviates2] endif else $ return,deviates1 end function mporbit2,params,date=date,rv=rv,err=err common dataset, dset, dstar ;Separate out the parinfo structure period = (2.0D * !DPI / params[0]) k = params[1] ecosw = params[2] esinw = params[3] ep_transit = params[4] v0 = params[5] slope = params[6] k2 = params[7] ;Apply the offset to the rv for i=8L,n_elements(params)-1 do begin num = i-8 rv[where(dset eq num)] = rv[where(dset eq num)] + params[i] endfor ;Split data up for two stars in sb2 case ;Determine if we are in the sb2 case sb2_elem = where(dstar eq 1, scount) date1 = date[where(dstar eq 0)] rv1 = rv[where(dstar eq 0)] err1 = err[where(dstar eq 0)] if(scount gt 0) then begin date2 = date[where(dstar eq 1)] rv2 = rv[where(dstar eq 1)] err2 = err[where(dstar eq 1)] endif ;Get eccentricity and argp ecc = sqrt(ecosw^2 + esinw^2) if(ecc eq 0.0) then $ argp = 0.0 $ else $ argp = atan(esinw,ecosw) ;Get the true anomaly using Karami (2009) astron. Nachr ;Determine phase phase = (date1 - ep_transit) / period phase = (phase - floor(phase)) f = 2 *!DPI * phase + (2*ecc-ecc^3/4)*sin(2*!DPI * phase) + 5.0/4.0*ecc^2*sin(4.0*!DPI*phase) + 13.0/12.0*ecc^2*sin(6*!DPI*phase) ;Find the radial velocity ; vel = v0 + k*(cos(argp + f) + ecc * cos(argp)) ovel1 = v0 + k*(cos(argp + f) +ecc*cos(argp)) + slope*date1 ;Determine the standard deviant deviates1 = (rv1 - ovel1) / err1 ;If sb2 case remember that the two orbits share identical values ;except for semi-amplitude and argp (which is PI off). if(scount gt 0) then begin enew = eanomaly(date2,ep_peri,period,ecc) cosf = (cos(enew) - ecc)/(1.D - ecc*cos(enew)) sinf = sqrt(1.D - ecc^2)*sin(enew)/(1.D - ecc*cos(enew)) ovel2 = v0 + k2*(cos(argp+!DPI +f) +ecc*cos(argp+!DPI)) + slope*date2 deviates2 = (rv2 - ovel2) / err2 return,[deviates1,deviates2] endif else $ return,deviates1 end ;A variant of the Lomb-Scargle periodogram ;Input the dates, rv, and error from the radial velocity curve ;Also determine the period window wind= [pmin:pmax] function periodogram1, $ date, $ rv, $ err, $ period1, $ base, $ hifac=hifac, $ nperiods=nperiods, $ over=over, $ wind=wind, $ builtin=builtin, $ color=color, $ fap=fap, $ fast=fast, $ is_numdatapoints_insufficent=is_numdatapoints_insufficent, $ linear=linear, $ pline=pline, $ plog=plog, $ post=post, $ notitle=notitle, $ screen=screen, $ verb=verb ;;============================================================================== ;; 0. Integrity check ;;============================================================================== ;;============================================================================== ;; 1. Set up internal variables ;;============================================================================== ;Necessary to pass info to the sinusoid and sinusoidl functions ;Warning this could impact parallel processing common frequency,sfreq ;Set default if(n_elements(wind) lt 2) then begin ;The lowest frequency is just 1/(t_max - t_min) i.e. the period equal ;to the observation time-span. f_low = 1.0D/(max(date) - min(date)) ;The highest frequency worth searching is sevaral times the Nyquist ;frequency the standard is 4 times. Nyquist = N_obs/(2*(t_max - t_min)) ;f_hi = hifac * Nyquist f_hi = f_low * n_elements(date)*hifac/2.0D ;Change back to periods wind = [1.0D/f_hi,1.0D/f_low] endif else begin f_hi = 1.0D/wind[0] f_low = 1.0D/wind[1] endelse ;Using a sampling rate based on equations from Press 1992 eq 13.8.9. ;I am oversampling by 10 if (n_elements(wind) lt 3) then $ f_step = 1.0D /(over * (max(date) - min(date))) $ else $ f_step = wind[2] ;I want this in radians/day w = 2PI * f w_step = f_step* 2.0D * !DPI ;Make an array of points to test wfreq = dindgen(floor((1.0D/wind[0] - 1.0D/wind[1])/f_step)) wfreq = wfreq * w_step + (2.0D * !DPI/wind[1]) if(keyword_set(verb)) then begin print, 'Beginning period search...' startp = systime(/seconds) endif ;;============================================================================== ;; 1.1. Check for insufficient data points. If so, skip to section 7. ;;============================================================================== if ~keyword_set(is_numdatapoints_insufficent) then begin ;;============================================================================== ;; 2. Do the fit ;;============================================================================== ;;============================================================================== ;; 2.1. The fit including a sine fit plus a linear slope ;;============================================================================== if(keyword_set(linear)) then begin ;Calculate the least-squares fit at each frequency to the equation ;f_j = A*cos(w*t_j) + B*sin(W*t_j) + Cx + D chi_circ = dblarr(n_elements(wfreq)) coeff_a = dblarr(n_elements(wfreq)) coeff_b = dblarr(n_elements(wfreq)) coeff_c = dblarr(n_elements(wfreq)) coeff_d = dblarr(n_elements(wfreq)) for i=0L, n_elements(wfreq)-1 do begin sfreq = wfreq[i] coeff = svdfit(date,rv,4,MEASURE_ERRORS=err,FUNCTION_NAME='sinusoidl',CHISQ=chi,YFIT=f_j,/double,status=status) ;Record the coeff coeff_a[i] = coeff[0] coeff_b[i] = coeff[1] coeff_c[i] = coeff[2] coeff_d[i] = coeff[3] ;Record the chisq values chi_circ[i] = chi endfor ;Now we want to determine the power in the periodogram. ;Do a linear fit to the data lin_coeff = svdfit(date,rv,2,MEASURE_ERRORS=err,CHISQ=chi_lin,/double,status=status) ;Now get the LS periodogram power Z(w) following Cumming 2004 ;with the modification that instead of the mean, I'm going to ;use the linear fit. ;z(w) = (chi_lin^2-chi_circ^2)*nu/(2*chi_circ^2) ;Where nu is degrees of fredom N-4 (due to linear trend) nu = n_elements(rv) - 4.0D z_ls = (chi_lin-chi_circ)*nu/(2.0D*chi_circ) ;Original lomb-scargle z definition ; z_ls = (chi_mean-chi_circ)/(2.0D) ;;============================================================================== ;; 2.2. The fit including only a sine fit ;;============================================================================== endif else begin ;Calculate the least-squares fit at each frequency to the equation ;f_j = A*cos(w*t_j) + B*sin(W*t_j) + C chi_circ = dblarr(n_elements(wfreq)) coeff_a = dblarr(n_elements(wfreq)) coeff_b = dblarr(n_elements(wfreq)) coeff_c = dblarr(n_elements(wfreq)) coeff_d = dblarr(n_elements(wfreq)) for i=0L, n_elements(wfreq)-1 do begin sfreq = wfreq[i] coeff = svdfit(date,rv,3,MEASURE_ERRORS=err,FUNCTION_NAME='sinusoid',CHISQ=chi,YFIT=f_j,status=status) ;Record the coeff. I make the coefficients compatible with the linear ;formulation by inserting 0.0 for the slope, which make the mean ;equibalent to the Y offset. coeff_a[i] = coeff[0] coeff_b[i] = coeff[1] coeff_c[i] = 0.0d coeff_d[i] = coeff[2] ;Record the chisq values chi_circ[i] = chi endfor ;Determine the chi^2 of a fit to the mean ;First determine aweighted mean ; wmean = total(rv)/n_elements(rv) wmean = total(rv * 1/err^2)/total(1/err^2) ;Get the chi^2 with respect to the weighted mean chi_mean = total((rv-wmean)^2/err^2) ;Now get the LS periodogram power Z(w) following Cumming 2004 ;z(w) = (chi_mean^2-chi_circ^2)*nu/(2*chi_circ^2) ;Where nu is degrees of fredom N-3 nu = n_elements(rv) - 3.0D ;This is a modified L-S power periodgram z2 in Baluev,2007, which is ;what Cumming used as well. z_ls = (chi_mean-chi_circ)*nu/(2.0D*chi_circ) ;Original lomb-scargle z definition ;z_ls = (chi_mean-chi_circ)/(2.0D) endelse if(keyword_set(verb)) then begin endp = systime(/seconds) print, 'End period search '+string(endp - startp,format='(F8.3,"secs")') print, 'Beginning FAP search...' startf = systime(/seconds) endif ;;============================================================================== ;; 3. Calculate the FAP's ;;============================================================================== ;Calculate the FAP precentage for the nperiods best z(w) best_z = dblarr(nperiods) ;Find the peaks best_index = peakfinder(z_ls,nperiods) ;Record info for each peak best_z = z_ls[best_index] best_wfreq = wfreq[best_index] best_period = 2.0D*!DPI / best_wfreq best_a = coeff_a[best_index] best_b = coeff_b[best_index] best_c = coeff_c[best_index] best_d = coeff_d[best_index] ;Setup array to hold faps best_fap = dblarr(n_elements(best_z)) ;Prep Big keyword if fast not selected if (not keyword_set(fast)) then $ big = get_cfap(best_z[0],date,wfreq,f_hi,nu,linear=linear,fast=fast) for i=0L,n_elements(best_z)-1 do begin best_fap[i] = get_cfap(best_z[i],date,wfreq,f_hi,nu,linear=linear,fast=fast,big=big) endfor ;;============================================================================== ;; 4. Output the FAP's to file ;;============================================================================== ;Create a FAP periodogram file if(keyword_set(fap)) then begin ;Okay let's open output file openw,lun,base+'.fap', /get_lun for i=0L,n_elements(wfreq)-1 do begin printf,lun,2*!DPI/wfreq[i], z_ls[i],get_cfap(z_ls[i],date,wfreq,f_hi,nu,linear=linear,fast=fast,big=big) endfor free_lun,lun endif ;;============================================================================== ;; 5. Plot the periodogram ;;============================================================================== ;;============================================================================== ;; 5.1. Determine the 50%, 99% and 99.9% cutoff for plotting ;;============================================================================== fap_50 = dblarr(n_elements(wfreq)) fap_01 = dblarr(n_elements(wfreq)) fap_001 = dblarr(n_elements(wfreq)) trouble = 0 zfinfo = get_czls(.50,date,wfreq,f_hi,nu,0.01,linear=linear,fast=fast,big=big) fap_50 = fap_50 + zfinfo[0] if(keyword_set(verb)) then $ print,'fap_50 ',zfinfo if(zfinfo[0] eq 0) then begin bottom = 0.01 trouble=1 endif $ else bottom = zfinfo[0] zfinfo = get_czls(.01,date,wfreq,f_hi,nu,bottom,linear=linear,fast=fast,big=big) fap_01 = fap_01 + zfinfo[0] if(keyword_set(verb)) then $ print,'fap_01 ',zfinfo if(zfinfo[0] eq 0) then begin bottom = 0.01 trouble=1 endif $ else bottom = zfinfo[0] zfinfo = get_czls(.001,date,wfreq,f_hi,nu,bottom,linear=linear,fast=fast,big=big) fap_001 = fap_001 + zfinfo[0] if(keyword_set(verb)) then $ print,'fap_001 ',zfinfo if(keyword_set(verb)) then begin endf = systime(/seconds) if(trouble) then $ print, 'FAP search failed, should not be trusted!' print, 'End FAP search '+string(endf - startf,format='(F8.3,"secs")') print, 'Begin Fitting...' endif if(keyword_set(verb)) then begin for i=0L,nperiods-1 do begin fmt = '(%"Period%d\nZ(w): %f\nFreq: %f\nPeriod: %f\nFAP: %f\nSlope: %f\nInter: %15.3f\n")' print,i+1,best_z[i],best_wfreq[i],best_period[i],best_fap[i],best_c[i],best_d[i],format=fmt endfor endif ;Run the IDL internal Lomb-Scargle periodogram I should consider fitting if(keyword_set(builtin)) then begin built = lnp_test(date,rv,hifac=hifac,ofac=over,wk1=bfreq,wk2=bz,jmax=jmax) fmt = '(%"Period1\nZ(w): %f\nFreq: %f\nPeriod: %f\nFAP: %f\n")' print,bz[jmax],2*!dpi*bfreq[jmax],1/bfreq[jmax],built[1],format=fmt !p.multi=[0,1,2,0,0] endif if(keyword_set(screen)) then begin ;Ensure that I'm using decomposed colors device,get_decomposed=sdecomp device,decompose=1 ;I want to plot in periodogram plot,2*!DPI/wfreq,z_ls,xtitle="Period", $ ytitle = "Z(w)", /nodata oplot,2*!DPI/wfreq,z_ls,psym=0,color='0000FF'XL oplot,2*!DPI/wfreq,fap_50,psym=0,linestyle=1,color='0000FF'XL oplot,2*!DPI/wfreq,fap_01,psym=0,linestyle=2,color='0000FF'XL device,decompose=sdecomp endif ;;============================================================================== ;; 5.2. Plot to screen ;;============================================================================== if(keyword_set(builtin) and keyword_set(screen)) then begin ;I want to plot in periodogram for built in ;Ensure that I'm using decomposed colors device,get_decomposed=sdecomp device,decompose=1 ;Percentiles above 99% should use - ln(F/M) ;Percentiles below 99% should use - ln(1-(1-F)^(1/M)) ;Again shouldn't do math in exponents. exp1 = 1/num_freq sfap_01 = dblarr(n_elements(bfreq)) - alog(.01/num_freq) sfap_50 = dblarr(n_elements(bfreq)) - alog(1-(1-0.5)^(exp1)) plot,1/bfreq,bz,xtitle="Period", $ ytitle = "Built in Z(w)",/nodata,xrange=[min(2*!DPI/wfreq),max(2*!DPI/wfreq)] oplot,1/bfreq,bz,psym=0,color='0000FF'XL oplot,2*!DPI/wfreq,sfap_50,psym=0,linestyle=1,color='0000FF'XL oplot,2*!DPI/wfreq,sfap_01,psym=0,linestyle=2,color='0000FF'XL device,decompose=sdecomp endif if(keyword_set(screen)) then begin !p.multi = 0 print,'Press key to continue...' junk = get_kbrd(1) endif ;;============================================================================== ;; 5.3. Plot to .ps file: periodogram for this RV curve ;;============================================================================== if(keyword_set(post)) then begin entry_device = !d.name set_plot,'PS' if(keyword_set(color)) then $ device,filename=base+'.ps',/landscape,/color $ else $ device,filename=base+'.ps',/landscape,color=0 ;Load graphics colors red = [0, 255, 0, 255, 0, 255, 0, 255, 0, 255, 255, 112, 219, 127, 0, 255] grn = [0, 0, 255, 255, 255, 0, 0, 255, 0, 187, 127, 219, 112, 127, 163, 171] blu = [0, 255, 255, 0, 0, 0, 255, 255, 115, 0, 127, 147, 219, 127, 255, 127] tvlct, red, grn, blu, 0 ;Set color names names = ['Black','Magenta','Cyan','Yellow','Green','Red','Blue','White',$ 'Navy','Gold','Pink','Aquamarine', 'Orchid', 'Gray','Sky','Beige'] ;Define plot symobol circles use psym 8 to make this work plotsym, 0, 1, /fill ;Make thing fancy for plotting, but remember what the old values were remember = [!P.FONT,!P.THICK,!P.CHARTHICK,!X.THICK,!Y.THICK,!P.CHARSIZE] ;Pick a fancy font !P.FONT=2 ;Darken up all the lines !P.THICK=4 !P.CHARTHICK=4 !P.CHARSIZE=1.35 !X.THICK=4 !Y.THICK=4 ;Make the actual plot ;Turn off titles if(keyword_set(notitle)) then $ mytitle = '' $ else $ mytitle = base ;Set up arrays to hold period names pnames = ['Period1','Period2','Period3','Period4','Period5'] plines = [5,1,4,6,7] pcolor = replicate(5,nperiods) pnames = pnames[0:nperiods-1] plines = plines[0:nperiods-1] pcolor = pcolor[0:nperiods-1] ;I want to plot in periodogram plot,2*!DPI/wfreq,z_ls,xtitle="Period", $ ytitle = "Z(w)",title=mytitle,/nodata,xlog=plog,yrange=[-1,max(z_ls)*1.1],ystyle=1 oplot,2*!DPI/wfreq,z_ls,psym=0 oplot,2*!DPI/wfreq,fap_50,psym=0 if (fap_01[0] ne 0) then $ oplot,2*!DPI/wfreq,fap_01,psym=0,linestyle=2 if (fap_001[0] ne 0) then $ oplot,2*!DPI/wfreq,fap_001,psym=0,linestyle=3 if(keyword_set(pline)) then begin for i=0L,nperiods-1 do begin plots,[best_period[i],best_period[i]],!y.crange,linestyle=plines[i],color=5,thick=6 endfor idlastro_legend,['50% FAP','1% FAP','0.1% FAP',pnames],linestyle=[0,2,3,plines],color=[0,0,0,pcolor],/right endif else $ idlastro_legend,['50% FAP','1% FAP','0.1% FAP'],linestyle=[0,2,3],/right ;Plot in HJD plot,date-2450000,rv,xtitle="Julian Date - 2450000", $ ytitle = "Radial Velocity (m/s)",title=mytitle,/nodata,/ynozero oplot,date-2450000,rv,psym=8,color=5 ;Add errorbars if appropriate errplot,date-2450000,rv-err,rv+err,color=5 ;Close file, set back to x11 device,landscape=0,/close_file fixps,base+'.ps' ;Set things back to normal set_plot,entry_device !P.FONT=remember[0] !P.THICK=remember[1] !P.CHARTHICK=remember[2] !P.CHARSIZE=remember[5] !X.THICK=remember[3] !Y.THICK=remember[4] !P.multi=0 endif ;;============================================================================== ;; 6. Return values ;;============================================================================== ;Return useful information if(keyword_set(builtin)) then begin ;Note I use my best sinusoid fit for the first value, and the builtin ;for the second value also the amplitude and v0 come from my best fit. period1.wind = [wind[0],wind[1],f_step] period1.z = [best_z[0],bz[jmax]] period1.wfreq = [best_wfreq[0],2*!dpi*bfreq[jmax]] period1.period=[best_period[0],1/bfreq[jmax]] period1.fap=[best_fap[0],built[1]] ;!!!!!!!!!I Should refit at the builtin frequency, not using my frequency!!!!! period1.amp=sqrt(best_a^2 + best_b^2) period1.slope = best_c period1.v0 = best_d endif else begin period1.wind = [wind[0],wind[1],f_step] period1.z = best_z period1.wfreq = best_wfreq period1.period=best_period period1.fap=best_fap period1.amp=sqrt(best_a^2 + best_b^2) period1.slope = best_c period1.v0 = best_d endelse ;;============================================================================== ;; 7. Return values and plot for the case of insufficient data ;;============================================================================== endif else begin if(keyword_set(post)) then begin entry_device = !d.name set_plot,'PS' if(keyword_set(color)) then $ device,filename=base+'.ps',/landscape,/color $ else $ device,filename=base+'.ps',/landscape,color=0 ;Load graphics colors red = [0, 255, 0, 255, 0, 255, 0, 255, 0, 255, 255, 112, 219, 127, 0, 255] grn = [0, 0, 255, 255, 255, 0, 0, 255, 0, 187, 127, 219, 112, 127, 163, 171] blu = [0, 255, 255, 0, 0, 0, 255, 255, 115, 0, 127, 147, 219, 127, 255, 127] tvlct, red, grn, blu, 0 ;Set color names names = ['Black','Magenta','Cyan','Yellow','Green','Red','Blue','White',$ 'Navy','Gold','Pink','Aquamarine', 'Orchid', 'Gray','Sky','Beige'] ;Define plot symobol circles use psym 8 to make this work plotsym, 0, 1, /fill ;Make thing fancy for plotting, but remember what the old values were remember = [!P.FONT,!P.THICK,!P.CHARTHICK,!X.THICK,!Y.THICK,!P.CHARSIZE] ;Pick a fancy font !P.FONT=2 ;Darken up all the lines !P.THICK=4 !P.CHARTHICK=4 !P.CHARSIZE=1.35 !X.THICK=4 !Y.THICK=4 ;Make the actual plot ;Turn off titles if(keyword_set(notitle)) then $ mytitle = '' $ else $ mytitle = base plot,2*!DPI/wfreq,0.0*wfreq,xtitle="Period", $ ytitle = "Z(w)",title=mytitle,/nodata,xlog=plog,yrange=[-1,1],ystyle=1 ;Plot in HJD plot,date-2450000,rv,xtitle="Julian Date - 2450000", $ ytitle = "Radial Velocity (m/s)",title=mytitle,/nodata,/ynozero oplot,date-2450000,rv,psym=8,color=5 ;Close file, set back to x11 device,landscape=0,/close_file fixps,base+'.ps' ;Set things back to normal set_plot,entry_device !P.FONT=remember[0] !P.THICK=remember[1] !P.CHARTHICK=remember[2] !P.CHARSIZE=remember[5] !X.THICK=remember[3] !Y.THICK=remember[4] !P.multi=0 endif ;Create a FAP periodogram file if(keyword_set(fap)) then begin ;Okay let's open output file openw,lun,base+'.fap', /get_lun for i=0L,n_elements(wfreq)-1 do begin printf,lun,2*!DPI/wfreq[i], !values.d_nan, !values.d_nan endfor free_lun,lun endif period1.period = !values.d_nan period1.amp[0] = (max(rv) - min(rv)) / 2 period1.v0[0] = period1.amp[0] + min(rv) period1.slope = !values.d_nan endelse return,period1 end ;A variant of the Lomb-Scargle periodogram which uses a fit to a ;Keplerian orbit ;Input the dates, rv, and error from the radial velocity curve ;Also determine the period window wind= [pmin:pmax] function periodogram2, $ date, $ rv, $ err, $ period1, $ period2, $ is_numdatapoints_insufficent=is_numdatapoints_insufficent, $ wind=wind, $ verb=verb, $ screen=screen ;;============================================================================== ;; 0. Integrity check ;;============================================================================== ;;============================================================================== ;; 0.1. Check for insufficient data points. If so, return ;; immediately, with return value equal to the input. ;;============================================================================== if is_numdatapoints_insufficent eq 1 then return,period1 ;;============================================================================== ;; 1. Set up internal variables ;;============================================================================== ;Set over sampling rate over = 3.0D hifac = 10.0D ;Using a sampling rate based on equations from Press 1992 eq 13.8.9. ;I am oversampling by 10 f_step = 1.0D /(over * (max(date) - min(date))) ;I want this in radians/day w = 2PI * f w_step = f_step* 2.0D * !DPI ;Make an array of points to test wfreq = dindgen(floor((1.0D/wind[0] - 1.0D/wind[1])/f_step)) wfreq = wfreq * w_step + (2.0D * !DPI/wind[1]) ;Setup parameters for fit parinfo = replicate({value:0.D, fixed:0, limited:[0,0], $ limits:[0.D,0]}, 6) parinfo[*].value = [wfreq[0],period1.amp[0],0.D,0.D,0.D,period1.v0[0]] ;frequency is fixed parinfo[0].fixed = 1 ;Eccentricity is between 0 and 1 and argp is between 0 and 2pi parinfo[1].limited = [1,0] parinfo[2].limited = [1,1] ;parinfo[4].limited = [1,1] parinfo[1].limits =[2*period1.amp[0]/3] parinfo[2].limits = [0.D,1.0D] ;parinfo[4].limits = [0.D,2.0D*!DPI] ;Setup structure to take x,y,and err functargs = {date:date, rv:rv, err:err} chi_kep = dblarr(n_elements(wfreq)) ;;============================================================================== ;; 2. Do the fit ;;============================================================================== for i=0L, n_elements(wfreq)-1 do begin parinfo[0].value = wfreq[i] params = MP_MPFIT('mporbit2', PARINFO=parinfo,FUNCTARGS=functargs,BESTNORM=chi,MAXITER=200,/QUIET) ;Record the chisq values chi_kep[i] = chi endfor ;Determine the chi^2 of a fit to the mean ;First determine a weighted mean ; wmean = total(rv)/n_elements(rv) wmean = total(rv * 1/err^2)/total(1/err^2) ;Get the chi^2 with respect to the weighted mean chi_mean = total((rv-wmean)^2/err^2) ;Now get the Keplarian periodogram power Z(w) following Cumming 2004 ;z_kep(w) = (chi_mean-chi_kep)*nu/(2*chi_kep) ;Where nu is degrees of fredom N-5 nu = n_elements(rv) - 5.0D z_kep = (chi_mean-chi_kep)*nu/(4.0D*chi_kep) ;Calculate the FAP precentage for the two best z(w) best_z = dblarr(2) best_z[0] = max(z_kep,best_index) ;To find the second best, elimate 2/3 the oversampling rate of points ;from either side of the first peak. if(best_index-ceil(over*2.0/3.0) lt 0) then $ gmin = 0 $ else $ gmin = best_index-ceil(over*2.0/3.0) if(best_index+ceil(over*2/3) gt n_elements(z_kep) - 1) then $ gmax = n_elements(z_kep) - 1 $ else $ gmax = best_index+ceil(over*2/3) sec_best = [z_kep[0:gmin],z_kep[gmax:*]] best_z[1] = max(sec_best) sec_index = where(best_z[1] eq z_kep,count) if(count ne 1) then $ message, 'More than one second best z_kep(w)' best_wfreq = [wfreq[best_index],wfreq[sec_index]] best_period = 2.0D*!DPI / best_wfreq if(keyword_set(verb)) then $ print,"Keplerian Periodogram:" fmt = '(%"Period1\t\t\tPeriod2\nZ(w): %f\tZ(w): %f\nFreq: %f\tFreq: %f\nPeriod: %f\tPeriod: %f\nFAP: %f\tFAP: %f\n")' ; print,best_z,best_wfreq,best_period,best_fap,format=fmt print,best_z,best_wfreq,best_period,[0,0],format=fmt ;;============================================================================== ;; 3. Plot to screen ;;============================================================================== if(keyword_set(screen)) then begin ;Ensure that I'm using decomposed colors device,get_decomposed=sdecomp device,decompose=1 ;I want to plot in periodogram plot,2*!DPI/wfreq,z_kep,xtitle="Period", $ ytitle = "Z_Kep(w)", /nodata oplot,2*!DPI/wfreq,z_kep,psym=0,color='0000FF'XL ; oplot,wfreq,fap_50,psym=0,linestyle=1,color='0000FF'XL ; oplot,wfreq,fap_01,psym=0,linestyle=2,color='0000FF'XL device,decompose=sdecomp print,'Press key to continue...' junk = get_kbrd(1) endif ;;============================================================================== ;; 4. Return useful information ;;============================================================================== period2.wind = wind period2.z = best_z period2.wfreq = best_wfreq period2.period=best_period ; period2.fap=best_fap period2.fap=[0.,0.] ; period2.amp=sqrt(best_a^2 + best_b^2) period2.amp=period1.amp ; period2.v0 = best_c period2.v0 = period1.v0 return,period2 end ;Inputs: JD,Visibility,file_basename ;Outputs: None function visplot,date,vis,base,notitle=notitle,color=color entry_device = !d.name set_plot,'PS' if(keyword_set(color)) then $ device,filename=base+'_v.ps',/landscape,/color $ else $ device,filename=base+'_v.ps',/landscape,color=0 ;Define plot symobol circles use psym 8 to make this work plotsym, 0, 1, /fill ;Make thing fancy for plotting, but remember what the old values were remember = [!P.FONT,!P.THICK,!P.CHARTHICK,!X.THICK,!Y.THICK,!P.CHARSIZE] ;Pick a fancy font !P.FONT=2 ;Darken up all the lines !P.THICK=4 !P.CHARTHICK=4 !P.CHARSIZE=1.35 !X.THICK=4 !Y.THICK=4 ;Make the actual plot ;Turn off titles if(keyword_set(notitle)) then $ mytitle = '' $ else $ mytitle = base ;Plot in Visibility as a function of RV plot,date-2450000,vis,xtitle="Julian Date - 2450000", $ ytitle = "Visibility",title=mytitle,/nodata,/ynozero oplot,date-2450000,vis,psym=8 ;Close file, set back to x11 device,landscape=0,/close_file fixps,base+'_v.ps' ;Set things back to normal set_plot,entry_device !P.FONT=remember[0] !P.THICK=remember[1] !P.CHARTHICK=remember[2] !P.CHARSIZE=remember[5] !X.THICK=remember[3] !Y.THICK=remember[4] !P.multi=0 return,1 end ;Inputs: JD,Period,Semi-Amplitude,Eccentricity,Epoch of periastron, ; Argument of periapsis,Intrisic Stellar velocity,Linear AB, ; Sinusodial ABC ;Outputs: Radial Velocity function rvorbit,data,period,k,ecosw,esinw,ep_peri,v0,slope,vlin=vlin,vsin=vsin ;Get eccentricity and argp ecc = sqrt(ecosw^2 + esinw^2) if(ecc eq 0.0) then $ argp = 0.0 $ else $ argp = atan(esinw,ecosw) ;I need E the eccentric anomaly. This will be returned from function eanomaly enew = eanomaly(data,ep_peri,period,ecc) ;Get the true anomaly in both sine and cosine. These are identical equations ;just used sin^2 + cos^2 = 1 to transform cosf = (cos(enew) - ecc)/(1.0D - ecc*cos(enew)) sinf = sqrt(1.0D - ecc^2)*sin(enew)/(1.0D - ecc*cos(enew)) ;Find the radial velocity ; vel = v0 + k*(cos(argp + f) + ecc * cos(argp)) ;Since I don't have f must use the cosine sum formula ovel = v0 + k*(cos(argp)*cosf - sin(argp)*sinf +ecc*cos(argp)) + slope*data ;Add linear trend if(n_elements(vlin) eq 2) then $ ovel += vlin[0] * data + vlin[1] ;Add Sinusoidal trend if(n_elements(vsin) eq 3) then $ ovel += vsin[0] * sin(2*!DPI*data/vsin[1] + vsin[2]) return,ovel end ;Necessary for creating the periodogram function sinusoid,X,M common frequency return,[ [sin(sfreq*x)],[cos(sfreq*x)],[1.0]] end ;Necessary for creating the periodogram function sinusoidl,X,M common frequency return,[ [sin(sfreq*x)],[cos(sfreq*x)],[x],[1.0]] end pro mprvfit, $ input, $ basenm=basenm, $ grid=grid, $ hifac=hifac, $ manual=manual, $ nperiods=nperiods, $ offset=offset, $ over=over, $ sb2=sb2, $ slope=slope, $ window=wind, $ yrange=yrange, $ yrrange=yrrange, $ builtin=builtin, $ circ=circ, $ color=color, $ csize=csize, $ fap=fap, $ fast=fast, $ force=force, $ kepler=kepler, $ linear=linear, $ loud=loud, $ masspl=masspl, $ massti=massti, $ nofit=nofit, $ nooff=nooff, $ notitle=notitle, $ phase=phase, $ plain=plain, $ pline=pline, $ plog=plog, $ post=post, $ resid=resid, $ rvsim=rvsim, $ screen=screen, $ showecc=showecc, $ unfold=unfold, $ verbose=verb ;;============================================================================== ;; 1. Set up internal variables ;;============================================================================== ;Make things standardized compile_opt idl2 ;Check to make sure there are enough parameters. if(n_params() ne 1) then $ message, "MPRVFIT,INPUT.CUR,BASENM=BASENM,CSIZE=CSIZE,GRID=[FSTEP,DIV],HIFAC=HIFAC,MANUAL=PERIOD,NPERIODS=NUM,OFFSET=[OFF1,OFF2],OVER=OVER,SB2=['IN1',IN2',etc],SLOPE=SLOPE,WIND=[PMIN,PMAX],YRANGE=[MIN,MAX],YRRANGE=[MIN,MAX],/BUILTIN,/CIRC,/COLOR,/FAP,/FAST,/FORCE,/KEPLER,/LINEAR,/LOUD,/MASSPL,/MASSTI,/NOFIT,/NOOFF,/NOTITLE,/PHASE,/PLAIN,/PLINE,/PLOG,/POST,/RESID,/RVSIM,/SCREEN,/SHOWECC,/UNFOLD,/VERB" ;Set over sampling rate if (n_elements(over) eq 0) then $ over = 5.0D ;The standard is 4, but I'm going to up it to 10 if (n_elements(hifac) eq 0) then $ hifac = 10.0D ;Set the default number of FSTEPS and DIV if (n_elements(grid) ne 2) then $ grid = [0,1] if (grid[1] lt 1) then $ grid[1] = 1 ;Set default character size if (n_elements(csize) ne 1) then $ csize = 1.0 ;Set default for nperiods if (n_elements(nperiods) lt 1) then begin if(n_elements(manual) lt 1) then $ nperiods = 2 $ else $ nperiods = n_elements(manual) endif if (n_elements(nperiods) gt 15) then $ nperiods = 15 ;Account for slope if(n_elements(slope) eq 0) then $ slope = dblarr(nperiods) ;Flag for insufficient data points for a fit is_numdatapoints_insufficent = 0 ;;============================================================================== ;; 2. Read in RV curves ;;============================================================================== ;;============================================================================== ;; 2.1. Prepare some variables ;;============================================================================== ;I need to know the total number of radial velocity points tot_lines = 0 staronelines = 0 for i=0L,n_elements(input)-1 do begin tot_lines += file_lines(input[i]) staronelines += file_lines(input[i]) endfor ;Need to include files involving SB2 as well if(n_elements(sb2) gt 0) then begin for i=0L,n_elements(sb2)-1 do $ tot_lines += file_lines(sb2[i]) startwolines = tot_lines - staronelines comb_input = [input,sb2] endif else $ comb_input = input ;Set the base name for output file (i.e. remove extension) if(n_elements(basenm) gt 0) then $ base = basenm $ else begin pos = stregex(input[0],'\.[^\.]*$') if(pos ne -1) then $ base = strmid(input[0],0,pos) $ else $ base = input[0] if(n_elements(comb_input) gt 1 ) then $ base = base + '_cmb' endelse ;The columns: date,rv,error,phase,visibility,file,source tmp_data = dblarr(7,tot_lines) current = 0 ;;============================================================================== ;; 2.2 Use readcol on the files ;;============================================================================== ;Open radial velocity curve file(s) for i=0L,n_elements(comb_input)-1 do begin if(file_test(comb_input[i]) ne 1) then $ message,comb_input[i] +' not found!' ;Use READCOL to get info into arrays, which will later be fed into the ;data array formating errors nlines = FILE_LINES(comb_input[i]) if (keyword_set(phase)) then $ readcol,comb_input[i],v1,v2,v3,v4,format='D,D,D,D',/silent $ else if (keyword_set(rvsim)) then $ readcol,comb_input[i],v1,v2,v3,v4,v5,format='D,D,D,D,A',/silent $ else $ readcol,comb_input[i],v1,v2,v3,format='D,D,D',/silent ;Combine data into large array if (keyword_set(phase)) then begin tmp_data[0,current:current+nlines-1] = v1 tmp_data[1,current:current+nlines-1] = v2 tmp_data[2,current:current+nlines-1] = v3 tmp_data[3,current:current+nlines-1] = v4 endif else if (keyword_set(rvsim)) then begin tmp_data[0,current:current+nlines-1] = v1 tmp_data[1,current:current+nlines-1] = v2 tmp_data[2,current:current+nlines-1] = v3 tmp_data[4,current:current+nlines-1] = v4 endif else begin tmp_data[0,current:current+nlines-1] = v1 tmp_data[1,current:current+nlines-1] = v2 tmp_data[2,current:current+nlines-1] = v3 endelse ;Mark which file the data came from tmp_data[5,current:current+nlines-1] = replicate(i,nlines) ;Record how many lines have been done current += nlines endfor if(n_elements(sb2) gt 0) then begin ;Mark which type star1 or star2 a file came from tmp_data[6,staronelines:tot_lines-1] = replicate(1,tot_lines-staronelines) endif ;;============================================================================== ;; 2.3. Condition the data ;;============================================================================== ;Get rid of bad points ;Remove non-finite data. good = where(finite(tmp_data[1,*]) and finite(tmp_data[2,*]) and (tmp_data[6,*] ne 1),ngood) data = tmp_data[*,good] ;Put all data into an array good = where(finite(tmp_data[1,*]) and finite(tmp_data[2,*]),ngood) data_all = tmp_data[*,good] ;Sort data by date data = data[*,sort(data[0,*])] data_all = data_all[*,sort(data_all[0,*])] ;;============================================================================== ;; 3. Plot to _v.ps file: visibility ;;============================================================================== if(keyword_set(rvsim) and keyword_set(post)) then $ junk = visplot(data[0,*],data[4,*],base,notitle=notitle,color=color) ;;============================================================================== ;; 4. Set up more internal variables ;;============================================================================== ;Setup return structure period1 = {wind:[0.0D,0.0D,0.0D],z:dblarr(nperiods),wfreq:dblarr(nperiods),period:dblarr(nperiods),fap:dblarr(nperiods),amp:dblarr(nperiods),v0:dblarr(nperiods),slope:dblarr(nperiods)} periodf = {wind:[0.0D,0.0D,0.0D],z:dblarr(nperiods),wfreq:dblarr(nperiods),period:dblarr(nperiods),fap:dblarr(nperiods),amp:dblarr(nperiods),v0:dblarr(nperiods),slope:dblarr(nperiods)} ;;============================================================================== ;; 5. Add any designated offsets that were passed in ;;============================================================================== ;Deal with offsets in data ;Setup my arrays rvcenter = dblarr(n_elements(comb_input)) corr_rv = dblarr(n_elements(data_all[1,*])) if(n_elements(comb_input) gt 1) then begin if(n_elements(offset) gt 0) then begin if(n_elements(offset) ne n_elements(comb_input)) then $ message,'Number of offsets does not equal number of files!' ;Apply the correct offset to each data set. corr_rv = dblarr(n_elements(data_all[1,*])) for i=0L,n_elements(input)-1 do begin corr_rv[where(data_all[5,*] eq i)] = reform(data_all[1,where(data_all[5,*] eq i)]) + offset[i] endfor endif else begin if (~keyword_set(nooff)) then begin ;Find my own offsets by finding the median of each dataset offset = dblarr(n_elements(comb_input)) for i=0L,n_elements(comb_input)-1 do begin rvcenter[i] = (max(reform(data_all[1,where(data_all[5,*] eq i)])) - min(reform(data_all[1,where(data_all[5,*] eq i)])))/2 + min(reform(data_all[1,where(data_all[5,*] eq i)])) offset[i] = rvcenter[0] - rvcenter[i] corr_rv[where(data_all[5,*] eq i)] = reform(data_all[1,where(data_all[5,*] eq i)]) + offset[i] endfor endif else begin offset = dblarr(n_elements(comb_input)) corr_rv = reform(data_all[1,*]) endelse endelse endif else begin offset = dblarr(n_elements(input)) corr_rv = reform(data_all[1,*]) endelse ;;============================================================================== ;; 6. Check integrity of RV curve: are there insufficient data for the fit? ;;============================================================================== ;Check to make sure that number of points is greater than or equal to ;the number of parameters if(keyword_set(linear) and n_elements(corr_rv) lt n_elements(input) + 7) then begin if ~keyword_set(FORCE) then begin print,n_elements(corr_rv),format='(%"%d is not enough points for a fit.")' return endif else begin is_numdatapoints_insufficent=1 endelse endif else begin if(n_elements(corr_rv) lt n_elements(input) + 6) then begin if ~keyword_set(FORCE) then begin print,n_elements(corr_rv),format='(%"%d is not enough points for a fit.")' return endif else begin is_numdatapoints_insufficent=1 endelse endif endelse ;;============================================================================== ;; 7. Set up more internal variables ;;============================================================================== ;Necessary to pass info to the mporbit function so that the offset can ;be applied to each data set. Also used to determine source star in ;sb2 case. ;Warning this could impact parallel processing common dataset,dset,dstar dset = reform(data_all[5,*]) dstar = reform(data_all[6,*]) ;Show offsets if(n_elements(comb_input) gt 1 and keyword_set(verb)) then begin print,'V0: ',rvcenter print,'Offsets: ',offset endif ;;============================================================================== ;; 8. Create the periodogram ;;============================================================================== ;;============================================================================== ;; 8.1. Automatic period assignment by periodogram ;;============================================================================== if(not keyword_set(manual)) then begin period1 = periodogram1( $ reform(data[0,*]), $ corr_rv[where(dstar eq 0)], $ reform(data[2,*]), $ period1, $ base, $ over=over, $ hifac=hifac, $ nperiods=nperiods, $ wind=wind, $ verb=verb, $ builtin=builtin, $ color=color, $ fap=fap, $ fast=fast, $ is_numdatapoints_insufficent=is_numdatapoints_insufficent, $ linear=linear, $ screen=screen, $ plog=plog, $ pline=pline, $ post=post, $ notitle=notitle) if(keyword_set(kepler)) then begin ;Create the modified periodogram period2 = period1 period2 = periodogram2( $ reform(data[0,*]), $ corr_rv, $ reform(data[2,*]), $ period1, $ period2, $ is_numdatapoints_insufficent=is_numdatapoints_insufficent, $ wind=wind, $ verb=verb, $ screen=screen) periodf = period2 endif else $ periodf = period1 ;;============================================================================== ;; 8.2. Manual period assignment by input argument ;;============================================================================== endif else begin if (n_elements(manual) ne nperiods) then $ message, "Number of periods must equal nperiods!" if (n_elements(slope) ne nperiods) then $ message, "Number of slopes must equal nperiods!" periodf.period = manual periodf.amp[0] = (max(corr_rv) - min(corr_rv)) / 2 periodf.v0[0] = periodf.amp[0] + min(corr_rv) periodf.slope = slope endelse ;;============================================================================== ;; 8.2. Check log file ;;============================================================================== ;Okay let's open output file openw,lun,base+'.log', /get_lun ;;============================================================================== ;; 8.3. Set up more internal variables ;;============================================================================== npts = n_elements(data_all[1,*]) if (keyword_set(linear)) then $ mylin = 1 $ else $ mylin = 0 ;;============================================================================== ;; 8.4. Print results of the period fit ;;============================================================================== ;!p.multi=[0,1,2,0,0] if(keyword_set(kepler)) then begin fmt = '(A,1X,F8.3,1X,F8.3,1X,F8.3,1X,F8.3)' printf,lun,base,period1.period,period2.period,format=fmt endif else begin printf,lun,'#File ',format='(A,$)' for i=0L,nperiods-1 do $ printf,lun,i+1,format='(" Period",I-2,$)' printf,lun,format='(" ",$)' for i=0L,nperiods-1 do $ printf,lun,i+1,format='("FAP",I-2," ",$)' printf,lun,' WIND1 WIND2 FSTEP GRID DIV OVER HIFAC LIN Npts' printf,lun,'#Fit Period Err Semi Err eCosW Err eSinW Err Epoch Err V0 Err Inter Err Slope Err Chi DOF Med_Err Med_Vis RMS M_Jup Err' fmt = '(A-40,1X,'+string(nperiods)+'(F9.3,1X),'+string(nperiods)+'(G9.4,1X),G9.4,1X,F8.3,1X,F8.3,1X,G9.4,1X,I3,1X,I3,1X,I3,1X,I3,1X,I2,1X,I3)' printf,lun,base,period1.period,period1.fap,period1.wind[0],period1.wind[1],period1.wind[2],grid[0],grid[1],over,hifac,mylin,npts,format=fmt endelse ;;============================================================================== ;; 8.5. Exit program here if no fitting required ;;============================================================================== ;We are done if nofit is set if(keyword_set(nofit)) then begin free_lun,lun return endif ;;============================================================================== ;; 9. Do Keplerian orbit fit. ;;============================================================================== ;;============================================================================== ;; 9.1. Set up more internal variables ;;============================================================================== ;Setup structure to store results nparam = 8 + n_elements(offset) anonfit = {params:dblarr(nparam),chi:0.D,cov:dblarr(nparam,nparam),dof:0.D,period:0.D,tperi:0.D,perror:dblarr(nparam)} ;Determine number of structures we need in this array the fit array ;is a 2-dim array where fit[j,k] j represents period peak being fit ;and k is the number of grid fits to that period fbest = intarr(nperiods) numfit = long(grid[0]*grid[1]) fit = replicate(anonfit,nperiods,numfit*2 + 1) g_step = dblarr(nperiods) ;Using a sampling rate based on equations from Press 1992 eq 13.8.9. ;I am oversampling by 10 f_step = 1.0D /(over * (max(data[0,*]) - min(data[0,*]))) p_step = f_step / grid[1] ;;============================================================================== ;; 9.2. Execute MPFIT to fit the Keplerian orbit. ;;============================================================================== ;We look at all peaks. for j=0L,nperiods-1 do begin for k=-numfit,numfit do begin ;Do final fit on the data fit[j,k+numfit]=finalfit( $ 1/(1/periodf.period[j]+k*p_step), $ reform(data_all[0,*]), $ reform(data_all[1,*]), $ reform(data_all[2,*]), $ periodf.amp[0], $ periodf.v0[0], $ periodf.slope[j], $ offset, $ loud=loud, $ linear=linear, $ circ=circ, $ manual=manual, $ dstar=dstar, $ nooff=nooff, $ is_numdatapoints_insufficent=is_numdatapoints_insufficent) endfor endfor ;Find the lowest chisq for each period for j=0L,nperiods-1 do begin junk1 = min(fit[j,*].chi,junk2,/nan) fbest[j] = junk2 ;Figure out the error in period g_step[j] = abs(1/(p_step + 1/fit[j,fbest[j]].period) - fit[j,fbest[j]].period) endfor ;;============================================================================== ;; 9.3. Convert the raw MPFIT results to desired derived quantities. ;;============================================================================== ;Now print information and make models for the best fits ecc = dblarr(nperiods) ecc_err = dblarr(nperiods) argp = dblarr(nperiods) argp_err = dblarr(nperiods) mplan = dblarr(nperiods) mplan_err = dblarr(nperiods) period_err = dblarr(nperiods) v0 = dblarr(nperiods) v0_err = dblarr(nperiods) best_per = dblarr(nperiods) best_amp = dblarr(nperiods) best_ecc = dblarr(nperiods) fixed_rv = dblarr(nperiods,n_elements(corr_rv[where(dstar eq 0)])) unfixed_rv = dblarr(nperiods,n_elements(corr_rv[where(dstar eq 0)])) mod_date = dblarr(nperiods,100) mod_phase = dblarr(nperiods,100) mod_rv = dblarr(nperiods,100) rvphase = dblarr(nperiods,n_elements(data[0,*])) ;;============================================================================== ;; 9.3.1. Loop over periods. ;;============================================================================== for j=0L,nperiods-1 do begin ;Get period error period_err[j] = 2.0 *!DPI * (fit[j,fbest[j]].params[0])^(-2) * fit[j,fbest[j]].perror[0] ;Get V0 middate = (data[0,n_elements(data[0,*])-1] - data[0,0])/2 + data[0,0] v0[j] = fit[j,fbest[j]].params[6]*middate+fit[j,fbest[j]].params[5] ;V0 error v0_err[j] = sqrt((fit[j,fbest[j]].perror[6])^2*middate^2 + (fit[j,fbest[j]].perror[5])^2 + 2* middate*fit[j,fbest[j]].cov[5,6]) ;Get eccentricity and argp and errors ecc[j] = sqrt(fit[j,fbest[j]].params[2]^2 + fit[j,fbest[j]].params[3]^2) ecc_err[j] = sqrt((fit[j,fbest[j]].params[2]^2 * fit[j,fbest[j]].perror[2]^2 + $ fit[j,fbest[j]].params[3]^2 * fit[j,fbest[j]].perror[3]^2)/ $ (fit[j,fbest[j]].params[2]^2 + fit[j,fbest[j]].params[3]^2)+ $ 2 *fit[j,fbest[j]].params[2] *fit[j,fbest[j]].params[3]/ $ (fit[j,fbest[j]].params[2]^2 + fit[j,fbest[j]].params[3]^2)* $ fit[j,fbest[j]].cov[2,3]) if(ecc[j] eq 0.0) then begin ecc_err[j] = 0.0 argp[j] = 0.0 argp_err[j] = 0.0 endif else begin argp[j] = atan(fit[j,fbest[j]].params[3],fit[j,fbest[j]].params[2]) if (argp[j] lt 0) then $ argp[j] = argp[j] + 2*!DPI argp_err[j] = sqrt(((fit[j,fbest[j]].perror[3]/fit[j,fbest[j]].params[2])^2+$ (fit[j,fbest[j]].params[3]^2*(fit[j,fbest[j]].perror[2]/$ fit[j,fbest[j]].params[2]^2)^2)/(1+ $ (fit[j,fbest[j]].params[3]^2 / $ fit[j,fbest[j]].params[2]^2))^2 + $ 2 *(1/(1+(fit[j,fbest[j]].params[3]^2 / $ fit[j,fbest[j]].params[2]^2)))^2 * $ abs(1/(fit[j,fbest[j]].params[3] *fit[j,fbest[j]].params[2]))$ * fit[j,fbest[j]].cov[2,3])) ;print,argp_err[j] endelse ;Get the mass for each fit assuming 1 solar mass host mplan[j] = fit[j,fbest[j]].params[1] * (1.0 * 1.98892e30)^(.66666)*(fit[j,fbest[j]].period*86400)^(.3333333)*sqrt(1-ecc[j]^2) mplan[j] /= (2*3.14159*6.673e-11)^(.333333) ;Make it all Jupiter Masses mplan[j] /= 1.8987e27 ;I will figure out the errors later mplan_err[j] = 0.0 ;Fix offset for rvcorr for l=0L,n_elements(input)-1 do begin fixed_rv[j,where(data[5,*] eq l)] = reform(data[1,where(data[5,*] eq l)]) + fit[j,fbest[j]].params[8+l] endfor ;Leave data intact for unfolded data unfixed_rv[j,*] = fixed_rv[j,*] ;;============================================================================== ;; 9.3.2. Prepare data for overplotting model curve. ;;============================================================================== ;;============================================================================== ;; 9.3.2.1. Check for insufficient data points. ;;============================================================================== if ~keyword_set(is_numdatapoints_insufficent) then begin ;Subtract off linear slope fixed_rv[j,*] = fixed_rv[j,*] - fit[j,fbest[j]].params[5] - fit[j,fbest[j]].params[6]*data[0,*] ;Create models based on these results. ;Create Phased model mod_date[j,*] = (fit[j,fbest[j]].period/100.0D * dindgen(100)) + fit[j,fbest[j]].params[4] new = (mod_date[j,*] - fit[j,fbest[j]].params[4])/fit[j,fbest[j]].period mod_phase[j,*] = new - floor(new) mod_rv[j,*] = rvorbit(mod_date[j,*],fit[j,fbest[j]].period,fit[j,fbest[j]].params[1],fit[j,fbest[j]].params[2],fit[j,fbest[j]].params[3],fit[j,fbest[j]].params[4],0,0) ;Phase the data new = (data[0,*] - fit[j,fbest[j]].params[4])/fit[j,fbest[j]].period rvphase[j,*] = new - floor(new) ;Get median error med_err = median(data[2,*]) ;If possible get median visibility if (keyword_set(rvsim)) then $ med_vis = median(data[4,*]) $ else $ med_vis = 0.0 ;Get model points for each data point ocmod_rv = rvorbit(data[0,*],fit[j,fbest[j]].period,fit[j,fbest[j]].params[1],fit[j,fbest[j]].params[2],fit[j,fbest[j]].params[3],fit[j,fbest[j]].params[4],fit[j,fbest[j]].params[5],fit[j,fbest[j]].params[6]) ocdiff = unfixed_rv[j,*] - ocmod_rv cresid = stddev(ocdiff) ;;============================================================================== ;; 9.3.2.2. If insufficient data points, manually make the quantities needed. ;;============================================================================== endif else begin ;Phase the data new = (data[0,*]) / (max(data[0,*])-min(data[0,*])) rvphase[j,*] = new - floor(new) ;Get median error med_err = median(data[2,*]) ;If possible get median visibility if (keyword_set(rvsim)) then $ med_vis = median(data[4,*]) $ else $ med_vis = 0.0 cresid = 0.0 endelse ;Print results to screen if(keyword_set(verb)) then begin print,j+1,format='(%"\nResults of Orbital Fit %d")' print,'Red Chi2: ',fit[j,fbest[j]].chi / fit[j,fbest[j]].dof print,'DOF: ',fit[j,fbest[j]].dof print,'RMS: ',cresid print,'Period: ',fit[j,fbest[j]].period,' +/- ',period_err[j] print,'Semi-Amp: ',fit[j,fbest[j]].params[1],' +/- ',fit[j,fbest[j]].perror[1] if(keyword_set(showecc)) then begin print,'ecc: ',ecc[j],' +/- ',ecc_err[j] print,'argp: ',argp[j],' +/- ',argp_err[j] endif else begin print,'ecos_w: ',fit[j,fbest[j]].params[2],' +/- ',fit[j,fbest[j]].perror[2] print,'esin_w: ',fit[j,fbest[j]].params[3],' +/- ',fit[j,fbest[j]].perror[3] endelse print,'Epoch Transit: ',fit[j,fbest[j]].params[4],' +/- ',fit[j,fbest[j]].perror[4] print,'Epoch Peri: ',fit[j,fbest[j]].tperi,' +/- ',fit[j,fbest[j]].perror[4] print,'V0: ',v0[j], ' +/- ',v0_err[j] print,'Y Inter: ',fit[j,fbest[j]].params[5],' +/- ',fit[j,fbest[j]].perror[5] print,'Slope: ',fit[j,fbest[j]].params[6],' +/- ',fit[j,fbest[j]].perror[6] print,'Mass_jup: ',mplan[j],' +/- ',mplan_err[j] if(n_elements(sb2) gt 0) then $ print,'Semi-amp2: ',fit[j,fbest[j]].params[7],' +/- ',fit[j,fbest[j]].perror[7] for k=0L,n_elements(comb_input) -1 do $ print,k,fit[j,fbest[j]].params[8+k],fit[j,fbest[j]].perror[8+k],format='(%"Offset %2d: %f +/- %f")' endif ;Fill arrays for mass plot best_per[j]=fit[j,fbest[j]].period best_amp[j]=fit[j,fbest[j]].params[1] best_ecc[j]=ecc[j] ;Print results of Keplerian fit fmt = '(A4,I1,1X,D,1X,D,1X,D,1X,D,1X,D,1X,D,1X,D,1X,D,1X,D,1X,D,1X,D,1X,D,1X,D,1X,D,1X,D,1X,D,1X,D,1X,I4,1X,F8.3,1X,F8.3,1X,G9.4,1X,F8.3,1X,F8.3)' printf,lun,'Fit',j,fit[j,fbest[j]].period,period_err[j],fit[j,fbest[j]].params[1],fit[j,fbest[j]].perror[1],fit[j,fbest[j]].params[2],fit[j,fbest[j]].perror[2],fit[j,fbest[j]].params[3],fit[j,fbest[j]].perror[3],fit[j,fbest[j]].params[4],fit[j,fbest[j]].perror[4],v0[j],v0_err[j],fit[j,fbest[j]].params[5],fit[j,fbest[j]].perror[5],fit[j,fbest[j]].params[6],fit[j,fbest[j]].perror[6],fit[j,fbest[j]].chi,fit[j,fbest[j]].dof,med_err,med_vis,cresid,mplan[j],mplan_err[j],format=fmt endfor ;;============================================================================== ;; 9.4. Plot the results. ;;============================================================================== ;Plot the data ;A list of symbols for different inputs sym = [0,4,4,5,3,0,8,4,5,3] sfill = [1,0,1,0,1,0,1,0,1,0] ;;============================================================================== ;; 9.4.1. Plot to _m.ps file: mass plot ;;============================================================================== ;Make mass plot if(keyword_set(post) and keyword_set(masspl)) then $ junk = massplot(best_per, $ best_amp, $ best_ecc, $ base, $ is_numdatapoints_insufficent=is_numdatapoints_insufficent, $ notitle=notitle, $ color=color) ;Load graphics colors red = [0, 255, 0, 255, 0, 255, 0, 255, 0, 255, 255, 112, 219, 127, 0, 255] grn = [0, 0, 255, 255, 255, 0, 0, 255, 0, 187, 127, 219, 112, 127, 163, 171] blu = [0, 255, 255, 0, 0, 0, 255, 255, 115, 0, 127, 147, 219, 127, 255, 127] tvlct, red, grn, blu, 0 ;Set color names names = ['Black','Magenta','Cyan','Yellow','Green','Red','Blue','White',$ 'Navy','Gold','Pink','Aquamarine', 'Orchid', 'Gray','Sky','Beige'] ;;============================================================================== ;; 9.4.2. Plot to screen. ;;============================================================================== ;Plot things to screen if wanted if(keyword_set(screen)) then begin ;Ensure that I'm using decomposed colors device,get_decomposed=sdecomp device,decompose=1 j=0 while(j lt nperiods) do begin ;I want to plot in phase plot,rvphase[j,*],fixed_rv[j,*],xtitle="Phase", $ ytitle = "Radial Velocity (m/s)",xrange=[0,2],/nodata,/ynozero for i=0L, n_elements(input)-1 do begin ;Define plot symbol circles use psym 8 to make this work plotsym, sym[i], 1, fill=sfill[i],thick=6 oplot,rvphase[j,where(data[5,*] eq i)],fixed_rv[j,where(data[5,*] eq i)],psym=8,color='0000FF'XL oplot,rvphase[j,where(data[5,*] eq i)]+1,fixed_rv[j,where(data[5,*] eq i)],psym=8,color='0000FF'XL errplot,rvphase[j,where(data[5,*] eq i)],fixed_rv[j,where(data[5,*] eq i)]-data[2,[where(data[5,*] eq i)]],fixed_rv[j,where(data[5,*] eq i)]+data[2,[where(data[5,*] eq i)]],color='0000FF'XL errplot,rvphase[j,where(data[5,*] eq i)]+1,fixed_rv[j,where(data[5,*] eq i)]-data[2,[where(data[5,*] eq i)]],fixed_rv[j,where(data[5,*] eq i)]+data[2,[where(data[5,*] eq i)]],color='0000FF'XL oplot,mod_phase[j,*],mod_rv[j,*],color='00FFFF'XL oplot,mod_phase[j,*]+1,mod_rv[j,*],color='00FFFF'XL print,'Press P for previous or any key for next: ' keys = get_kbrd(1) if (keys eq 'p') then $ j = j - 1 $ else $ j = j + 1 endfor ;Set things back to normal device,decompose=sdecomp endwhile endif ;To make it easier to figure out what is being done when, I am ;splitting the folded and unfolded plots into two separate for loops: ;;============================================================================== ;; 9.4.3. Plot to _c.ps file: phase-folded ;;============================================================================== if(keyword_set(post)) then begin entry_device = !d.name set_plot,'PS' if(keyword_set(color)) then $ device,filename=base+'_c.ps',/landscape,/color $ else $ device,filename=base+'_c.ps',/landscape,color=0 ;Make thing fancy for plotting, but remember what the old values were remember = [!P.FONT,!P.THICK,!P.CHARTHICK,!X.THICK,!Y.THICK,!P.CHARSIZE] ;Pick a fancy font !P.FONT=2 ;Darken up all the lines !P.THICK=4 !P.CHARTHICK=4 !P.CHARSIZE=csize !X.THICK=4 !Y.THICK=4 ;;============================================================================== ;; 9.4.3.1. Check for insufficient data points. ;;============================================================================== if ~keyword_set(is_numdatapoints_insufficent) then begin ;;============================================================================== ;; 9.4.3.2. Make a phased RV curve for each period ;;============================================================================== for j=0L,nperiods-1 do begin ;I want to plot in phase ;Turn off titles if(keyword_set(notitle)) then $ mytitle = '' $ else begin if (keyword_set(massti)) then $ mytitle = base + string(mplan[j],format='(" Mass: ",F7.2," M_Jup")') $ else $ mytitle = base endelse ;Change offset base on fitted semi-amp if(fit[j,fbest[j]].params[1] lt 2000) then $ poff = 200 $ else $ poff = 2000 if(keyword_set(yrange)) then begin u_cdiff1 = yrange[1] - yrange[0] u_cdiff2 = yrange[0] - min(reform(unfixed_rv[j,*])) + min(reform(fixed_rv[j,*])) my_yrange = [u_cdiff2,u_cdiff1 + u_cdiff2] my_ystyle = 1 endif else begin my_yrange = [min([reform(fixed_rv[j,*]),reform(mod_rv[j,*])]- poff/4), $ max([reform(fixed_rv[j,*]),reform(mod_rv[j,*])])+poff] my_ystyle = 0 endelse plot,rvphase[j,*],fixed_rv[j,*],xtitle="Phase", $ ytitle = "Radial Velocity (m/s)",xrange=[0,2], $ yrange=my_yrange,ystyle = my_ystyle,title=mytitle,/nodata,/ynozero ;+string(j+1,format='(" Peak: ",I)'),/nodata,/ynozero for i=0L, n_elements(input)-1 do begin ;Define plot symbol circles use psym 8 to make this work plotsym, sym[i], 1, fill=sfill[i],thick=6 oplot,rvphase[j,where(data[5,*] eq i)],fixed_rv[j,where(data[5,*] eq i)],psym=8,color=i oplot,rvphase[j,where(data[5,*] eq i)]+1,fixed_rv[j,where(data[5,*] eq i)],psym=8,color=i errplot,rvphase[j,where(data[5,*] eq i)],fixed_rv[j,where(data[5,*] eq i)]-data[2,[where(data[5,*] eq i)]],fixed_rv[j,where(data[5,*] eq i)]+data[2,[where(data[5,*] eq i)]],color=i errplot,rvphase[j,where(data[5,*] eq i)]+1,fixed_rv[j,where(data[5,*] eq i)]-data[2,[where(data[5,*] eq i)]],fixed_rv[j,where(data[5,*] eq i)]+data[2,[where(data[5,*] eq i)]],color=i endfor ;Plot two cycles (phase 0 to 1 and 1 to 2). oplot,[reform(mod_phase[j,*]),reform(mod_phase[j,*]+1)],[reform(mod_rv[j,*]),reform(mod_rv[j,*])],color=6 if (~keyword_set(plain)) then begin xyouts,.15,.9,string(fit[j,fbest[j]].chi/fit[j,fbest[j]].dof,fit[j,fbest[j]].dof,format='("Red Chi2:",9X,F9.3," DOF:",1X,I3)'),charsize=1.25,/normal xyouts,.15,.87,string(fit[j,fbest[j]].period,g_step[j],format='("Period:",12X,F10.3," +/- ",g10.3)'),charsize=1.25,/normal xyouts,.15,.84,string(fit[j,fbest[j]].params[1],fit[j,fbest[j]].perror[1],format='("Semi-Amp:",7X,F9.3," +/- ",F9.3)'),charsize=1.25,/normal if(keyword_set(showecc)) then begin xyouts,.15,.81,string(ecc[j],ecc_err[j],format='("Ecc:",19X,F9.3," +/- ",F9.3)'),charsize=1.25,/normal xyouts,.15,.78,string(argp[j],argp_err[j],format='("Argp:",18X,F9.3," +/- ",F9.3)'),charsize=1.25,/normal endif else begin xyouts,.15,.81,string(fit[j,fbest[j]].params[2],fit[j,fbest[j]].perror[2],format='("ecosw:",14X,F9.3," +/- ",F9.3)'),charsize=1.25,/normal xyouts,.15,.78,string(fit[j,fbest[j]].params[3],fit[j,fbest[j]].perror[3],format='("esinw:",15X,F9.3," +/- ",F9.3)'),charsize=1.25,/normal endelse xyouts,.55,.9,string(med_err,fit[j,fbest[j]].params[1]/med_err,format='("Med !9s!3:",11X,F9.3," K/!9s!3:",1X,F9.3)'),charsize=1.25,/normal xyouts,.55,.87,string(fit[j,fbest[j]].params[4]-2450000,fit[j,fbest[j]].perror[4],format='("Epoch Peri:",2X,F9.3," +/- ",g10.4)'),charsize=1.25,/normal xyouts,.55,.84,string(v0[j],v0_err[j],format='("V_0:",13X,F9.2," +/- ",g10.4)'),charsize=1.25,/normal xyouts,.55,.81,string(fit[j,fbest[j]].params[5],fit[j,fbest[j]].perror[5],format='("Y Inter:",9X,g9.3," +/- ",g10.4)'),charsize=1.25,/normal xyouts,.55,.78,string(fit[j,fbest[j]].params[6],fit[j,fbest[j]].perror[6],format='("Slope:",13X,F9.3," +/- ",g10.4)'),charsize=1.25,/normal endif endfor ;;============================================================================== ;; 9.4.3.3. If not enough points for a fit, fake the phased RV curve ;; (phase calculated using period = time baseline) ;;============================================================================== endif else begin for j=0L,nperiods-1 do begin plot,rvphase[j,*],fixed_rv[j,*],xtitle="Phase", $ ytitle = "Radial Velocity (m/s)",xrange=[0,2], $ yrange=my_yrange,ystyle = my_ystyle,title=mytitle,/nodata,/ynozero oplot,rvphase[j,*],fixed_rv[j,*],psym=8 endfor endelse ;Set things back to normal set_plot,entry_device !P.FONT=remember[0] !P.THICK=remember[1] !P.CHARTHICK=remember[2] !P.CHARSIZE=remember[5] !X.THICK=remember[3] !Y.THICK=remember[4] ;Close file, set back to x11 set_plot,'PS' device,landscape=0,/close_file fixps,base+'_c.ps' set_plot,entry_device endif ;;============================================================================== ;; 9.4.4. Plot to _u.ps file: unfolded ;;============================================================================== ;Make the actual plot if(keyword_set(unfold) and keyword_set(post)) then begin entry_device = !d.name set_plot,'PS' if(keyword_set(color)) then $ device,filename=base+'_u.ps',/landscape,/color $ else $ device,filename=base+'_u.ps',/landscape,color=0 ;Define plot symbol circles use psym 8 to make this work plotsym, 0, 1, /fill ;Make thing fancy for plotting, but remember what the old values were remember = [!P.FONT,!P.THICK,!P.CHARTHICK,!X.THICK,!Y.THICK,!P.CHARSIZE] ;Pick a fancy font !P.FONT=2 ;Darken up all the lines !P.THICK=4 !P.CHARTHICK=4 !P.CHARSIZE=csize !X.THICK=4 !Y.THICK=4 for j=0L,nperiods-1 do begin ;Turn off titles if(keyword_set(notitle)) then $ mytitle = '' $ else begin if (keyword_set(massti)) then $ mytitle = base + string(mplan[j],format='(" Mass: ",F7.2," M_Jup")') $ else $ mytitle = base endelse ;;============================================================================== ;; 9.4.4.1. Check for insufficient data points. ;;============================================================================== if ~keyword_set(is_numdatapoints_insufficent) then begin ;Create the models for unfolded data range = data[0,n_elements(data[0,*])-1] - data[0,0] dpoints = range/fit[j,fbest[j]].period * 100 unmod_date = (fit[j,fbest[j]].period/100 * dindgen(dpoints)) + data[0,0] unmod_rv = rvorbit(unmod_date,fit[j,fbest[j]].period,fit[j,fbest[j]].params[1],fit[j,fbest[j]].params[2],fit[j,fbest[j]].params[3],fit[j,fbest[j]].params[4],fit[j,fbest[j]].params[5],fit[j,fbest[j]].params[6]) ;Change offset base on fitted semi-amp if(fit[j,fbest[j]].params[1] lt 2000) then $ poff = 200 $ else $ poff = 2000 endif else begin range = data[0,n_elements(data[0,*])-1] - data[0,0] dpoints = range/(max(data[0,*])-min(data[0,*])) * 100 unmod_date = ((max(data[0,*])-min(data[0,*]))/100 * dindgen(dpoints)) + data[0,0] unmod_rv = 0.0*unmod_date poff = 200 endelse if(keyword_set(yrange)) then begin my_yrange = yrange my_ystyle = 1 endif else begin my_yrange = [min([reform(unfixed_rv[j,*]),unmod_rv]- poff/4),$ max([reform(unfixed_rv[j,*]),unmod_rv])+poff] my_ystyle = 0 endelse plot,data[0,*]-2450000,unfixed_rv[j,*],xtitle="JD-2450000",$ ytitle = "Radial Velocity (m/s)",$ yrange=my_yrange,ystyle=my_ystyle,$ title=mytitle,/nodata,/ynozero oplot,unmod_date-2450000,unmod_rv,color=6 for i=0L, n_elements(input)-1 do begin ;Define plot symbol circles use psym 8 to make this work plotsym, sym[i], 1, fill=sfill[i],thick=6 oplot,data[0,where(data[5,*] eq i)]-2450000,unfixed_rv[j,where(data[5,*] eq i)],psym=8,color=i errplot,data[0,where(data[5,*] eq i)]-2450000,unfixed_rv[j,where(data[5,*] eq i)]-data[2,[where(data[5,*] eq i)]],unfixed_rv[j,where(data[5,*] eq i)]+data[2,[where(data[5,*] eq i)]],color=i endfor if (~keyword_set(plain)) then begin xyouts,.15,.9,string(fit[j,fbest[j]].chi/fit[j,fbest[j]].dof,fit[j,fbest[j]].dof,format='("Red Chi2:",9X,F9.3," DOF:",1X,I3)'),charsize=1.25,/normal xyouts,.15,.87,string(fit[j,fbest[j]].period,g_step[j],format='("Period:",12X,F10.3," +/- ",g10.3)'),charsize=1.25,/normal xyouts,.15,.84,string(fit[j,fbest[j]].params[1],fit[j,fbest[j]].perror[1],format='("Semi-Amp:",7X,F9.3," +/- ",F9.3)'),charsize=1.25,/normal if(keyword_set(showecc)) then begin xyouts,.15,.81,string(ecc[j],ecc_err[j],format='("Ecc:",19X,F9.3," +/- ",F9.3)'),charsize=1.25,/normal xyouts,.15,.78,string(argp[j],argp_err[j],format='("Argp:",18X,F9.3," +/- ",F9.3)'),charsize=1.25,/normal endif else begin xyouts,.15,.81,string(fit[j,fbest[j]].params[2],fit[j,fbest[j]].perror[2],format='("ecosw:",14X,F9.3," +/- ",F9.3)'),charsize=1.25,/normal xyouts,.15,.78,string(fit[j,fbest[j]].params[3],fit[j,fbest[j]].perror[3],format='("esinw:",15X,F9.3," +/- ",F9.3)'),charsize=1.25,/normal endelse xyouts,.55,.9,string(med_err,fit[j,fbest[j]].params[1]/med_err,format='("Med !9s!3:",11X,F9.3," K/!9s!3:",1X,F9.3)'),charsize=1.25,/normal xyouts,.55,.87,string(fit[j,fbest[j]].params[4]-2450000,fit[j,fbest[j]].perror[4],format='("Epoch Peri:",2X,F9.3," +/- ",g10.4)'),charsize=1.25,/normal xyouts,.55,.84,string(v0[j],v0_err[j],format='("V_0:",13X,F9.2," +/- ",g10.4)'),charsize=1.25,/normal xyouts,.55,.81,string(fit[j,fbest[j]].params[5],fit[j,fbest[j]].perror[5],format='("Y Inter:",9X,g9.3," +/- ",g10.4)'),charsize=1.25,/normal xyouts,.55,.78,string(fit[j,fbest[j]].params[6],fit[j,fbest[j]].perror[6],format='("Slope:",13X,F9.3," +/- ",g10.4)'),charsize=1.25,/normal endif endfor ;Set things back to normal set_plot,entry_device !P.FONT=remember[0] !P.THICK=remember[1] !P.CHARTHICK=remember[2] !P.CHARSIZE=remember[5] !X.THICK=remember[3] !Y.THICK=remember[4] ;Close file, set back to x11 set_plot,'PS' device,landscape=0,/close_file fixps,base+'_u.ps' set_plot,entry_device endif ;;============================================================================== ;; 9.4.5. Plot to _o.ps file: residuals ;;============================================================================== ;Make the Residual (O-C) plot if(keyword_set(resid) and keyword_set(post)) then begin ;Okay let's open a residual output file openw,lun1,base+'.res', /get_lun printf,lun1,'#JD RV(m/s) Mod(m/s) Diff(m/s) File Fit' entry_device = !d.name set_plot,'PS' if(keyword_set(color)) then $ device,filename=base+'_o.ps',/landscape,/color $ else $ device,filename=base+'_o.ps',/landscape,color=0 ;Define plot symbol circles use psym 8 to make this work plotsym, 0, 1, /fill ;Make thing fancy for plotting, but remember what the old values were remember = [!P.FONT,!P.THICK,!P.CHARTHICK,!X.THICK,!Y.THICK,!P.CHARSIZE] ;Pick a fancy font !P.FONT=2 ;Darken up all the lines !P.THICK=4 !P.CHARTHICK=4 !P.CHARSIZE=csize !X.THICK=4 !Y.THICK=4 for j=0L,nperiods-1 do begin ;Turn off titles if(keyword_set(notitle)) then $ mytitle = '' $ else begin if (keyword_set(massti)) then $ mytitle = base + string(mplan[j],format='(" Mass: ",F7.2," M_Jup")') $ else $ mytitle = base endelse ;Get model points for each data point ocmod_rv = rvorbit(data[0,*],fit[j,fbest[j]].period,fit[j,fbest[j]].params[1],fit[j,fbest[j]].params[2],fit[j,fbest[j]].params[3],fit[j,fbest[j]].params[4],fit[j,fbest[j]].params[5],fit[j,fbest[j]].params[6]) ocdiff = unfixed_rv[j,*] - ocmod_rv cresid = stddev(ocdiff) ;Write out residual information fmt = '(F14.6,1X,G10.6,1X,G10.6,1X,G10.6,1X,I2,1X,I2)' for k=0L,n_elements(ocmod_rv)-1 do begin printf,lun1,data[0,k],unfixed_rv[j,k],ocmod_rv[k],ocdiff[k],data[5,k],j,format=fmt endfor ;Change offset base on oc difference if(max(ocdiff) lt 2000) then $ poff = 200 $ else $ poff = 2000 IF(KEYWORD_SET(yrrange)) THEN BEGIN my_yrrange = yrrange my_ystyle = 1 ENDIF ELSE BEGIN my_yrrange = [min(ocdiff) - poff/4,max(ocdiff) + poff] my_ystyle = 2 ENDELSE plot,data[0,*]-2450000,ocdiff,xtitle="JD-2450000",$ ytitle = "O - C (m/s)",$ yrange=my_yrrange,$ ystyle=my_ystyle,$ title=mytitle,/nodata,/ynozero ;Draw line at 0 plots,!x.crange,[0,0],linestyle=1 for i=0L, n_elements(input)-1 do begin ;Define plot symbol circles use psym 8 to make this work plotsym, sym[i], 1, fill=sfill[i],thick=6 oplot,data[0,where(data[5,*] eq i)]-2450000,ocdiff[where(data[5,*] eq i)],psym=8,color=i errplot,data[0,where(data[5,*] eq i)]-2450000,ocdiff[where(data[5,*] eq i)]-data[2,[where(data[5,*] eq i)]],ocdiff[where(data[5,*] eq i)]+data[2,[where(data[5,*] eq i)]],color=i endfor if (~keyword_set(plain)) then begin xyouts,.15,.9,string(fit[j,fbest[j]].chi/fit[j,fbest[j]].dof,fit[j,fbest[j]].dof,format='("Red Chi2:",9X,F9.3," DOF:",1X,I3)'),charsize=1.25,/normal xyouts,.15,.87,string(fit[j,fbest[j]].period,g_step[j],format='("Period:",12X,F10.3," +/- ",g10.3)'),charsize=1.25,/normal xyouts,.15,.84,string(fit[j,fbest[j]].params[1],fit[j,fbest[j]].perror[1],format='("Semi-Amp:",7X,F9.3," +/- ",F9.3)'),charsize=1.25,/normal if(keyword_set(showecc)) then begin xyouts,.15,.81,string(ecc[j],ecc_err[j],format='("Ecc:",19X,F9.3," +/- ",F9.3)'),charsize=1.25,/normal xyouts,.15,.78,string(argp[j],argp_err[j],format='("Argp:",18X,F9.3," +/- ",F9.3)'),charsize=1.25,/normal endif else begin xyouts,.15,.81,string(fit[j,fbest[j]].params[2],fit[j,fbest[j]].perror[2],format='("ecosw:",14X,F9.3," +/- ",F9.3)'),charsize=1.25,/normal xyouts,.15,.78,string(fit[j,fbest[j]].params[3],fit[j,fbest[j]].perror[3],format='("esinw:",15X,F9.3," +/- ",F9.3)'),charsize=1.25,/normal endelse xyouts,.15,.75,string(cresid,format='("RMS:",15X,F9.3)'),charsize=1.25,/normal xyouts,.55,.9,string(med_err,fit[j,fbest[j]].params[1]/med_err,format='("Med !9s!3:",11X,F9.3," K/!9s!3:",1X,F9.3)'),charsize=1.25,/normal xyouts,.55,.87,string(fit[j,fbest[j]].params[4]-2450000,fit[j,fbest[j]].perror[4],format='("Epoch Peri:",2X,F9.3," +/- ",g10.4)'),charsize=1.25,/normal xyouts,.55,.84,string(v0[j],v0_err[j],format='("V_0:",13X,F9.2," +/- ",g10.4)'),charsize=1.25,/normal xyouts,.55,.81,string(fit[j,fbest[j]].params[5],fit[j,fbest[j]].perror[5],format='("Y Inter:",9X,g9.3," +/- ",g10.4)'),charsize=1.25,/normal xyouts,.55,.78,string(fit[j,fbest[j]].params[6],fit[j,fbest[j]].perror[6],format='("Slope:",13X,F9.3," +/- ",g10.4)'),charsize=1.25,/normal endif endfor ;Set things back to normal set_plot,entry_device !P.FONT=remember[0] !P.THICK=remember[1] !P.CHARTHICK=remember[2] !P.CHARSIZE=remember[5] !X.THICK=remember[3] !Y.THICK=remember[4] ;Close file, set back to x11 set_plot,'PS' device,landscape=0,/close_file fixps,base+'_o.ps' set_plot,entry_device free_lun,lun1 endif ;Close logfile !p.multi=0 free_lun,lun end