Demo entry 3548685

periodograma

   

Submitted by ReddTea on Jan 17, 2016 at 18:01
Language: IDL. Code size: 13.1 kB.

Function harps_lomb_peri_func,data,planet,activity=activity,bisector=bisector,random=random,noboot=noboot,$
sampling=sampling,fap=fap,log=log,mark=mark,hipp=hipp,range=range,max_search=max_search,yaxis=yaxis,
$
show_planets=show_planets,fwhm=fwhm,cont=cont,no_plot=no_plot
;Compute the Lomb-Scargle periodogram of an input RV dataset!!
;FILE is a list of files that contain the RVs, use from the
;SYSTEMIC_FILES/ directory to get the name right!! The ACTIVITY or
;BISECTOR keywords should be set if you want to look at the activity
;periodicities (S-indices for instance) or the bisector velocities to
;see if the periods correlate at all with the period of the RVs!! You
;need to set this as they are in the other systemic directory!! The
;keyword RANDOM must be set to the number of random realisations we
;need to do to get a good FAP, generally 10000 is good and is the
;default!! The keyword NOBOOT should be set to stop the bootstrap FAP
;test from being performed!! If keyword_set SAMPLING we run the
;timestamps through the periodogram analysis!! The keyword FAP should
;be set to a value of the FAP if we know it!! Keyword LOG allows the
;x-axes in the output plots to be in log scale!! Keyword MARK can be
;used to overplot a line showing an input period given by mark!!
;If the HIPP keyword is set we are using the HIPPARCOS photometry to
;hunt for periods in the raw HIPPARCOS data!!
;The PLANET input is set to the number of planets for the plot i.e.
;if there are no planet signals subtracted then this is zero!!
;
;The RANGE keyword should be set to a two element vector with a
;minimum and maximum number of days to print out, e.g. range=[1.,5.]
;would only print out the powers between frequencies of 1-5 days!!
;
;The MAX_SEARCH keyword is a two element vector to get the max peak to
;test the FAP, which helps if, for instance, the alias is stronger
;than the actual doppler peak!!
;
;The YAXIS keyword should be set to the height of the y-axis we want
;to plot. This is not necessary but useful if you want to directly
;compare two periodograms on the same scale!!
;
;The SHOW_PLANETS keyword can be set to plot the period of the planets
;on the periodogram!!
;
;Set the FWHM keyword to tell it that the file is my HARPS output file
;and so can do the fwhm instead!!
;
;Set the CONT keyword to tell it that the file is my HARPS output file
;and so can do the contrast instead!!
;
;NO_PLOT stops the plotting!!
;
;Function returns the period with the maximum power, the power, the
;probability for this power, and the alias of this frequency.


!X.THICK=4
!Y.THICK=4
!P.THICK=4
!P.CHARSIZE=2
!P.CHARTHICK=4
!Y.TITLE=2


;readcol_old,file,format='(a)',files

;Check if the periodograms directory is there!!
if keyword_set(no_plot) then goto,plot_skip
res=file_test('periodograms',/directory)
if res EQ 0 then spawn,'mkdir periodograms'
res=file_test('periodograms/bisector',/directory)
if res EQ 0 then spawn,'mkdir periodograms/bisector'
res=file_test('periodograms/hipparcos',/directory)
if res EQ 0 then spawn,'mkdir periodograms/hipparcos'

set_plot,'ps'

plot_skip:

t=data[*,0,0]
rv=data[*,1,0]
er=data[*,2,0]

if n_elements(t) LT 2 then goto,skip

if ( n_elements(t) NE n_elements(rv) ) then begin
print,"Inputs are not the same length!!"
;return
stop
endif

if keyword_set(sampling) then rv=t

;First we want to subtract the mean, compute the variance determine
;the range of frequencies to use and then initialise Pn!!

;Frequencies between 1-10 days in steps of 0.001!!
;To increase or decrease the period we increase or
;decrease the denominator and also when we increase
;the denominator we must decrease the offset by the same
;fraction i.e. increase the denominator by an order of
;magnitude, then we decrease the offset by the same
;order of magnitude.
;freq=(dindgen(9001.)/1000.D) + 0.0001D
freq=(dindgen(10000.)/10000.) + 0.0001d

;range1=1./range
;range1=range
;diff=abs(range1[0]-range1[1])
;freq=(freq1*diff)+range1[0]; - 0.0001d

;freq=1./((dindgen(10000)/(10000./diff))+range1[0])

;print, freq[0], freq[n_elements(freq)-1]
;print, freq[*]

z = rv - mean(rv)
var = stddev(rv)
n = n_elements(freq)
Pn = fltarr(n)

for i=0,n-1 do begin

w=2.D*!dpi*freq[i]
if w GT 0. then begin

twt = 2.D*w*t
tau = (atan(total(sin(twt)),total(cos(twt)))/2.D)/w
wtmt = w*(t - tau)
pn[i] = (total(z*cos(wtmt))^2.D)/total(cos(wtmt)^2.D) + $
(total(z*sin(wtmt))^2.D)/total(sin(wtmt)^2.D)

endif
if w LE 0. then pn[i] = (total(z*t)^2.D)/total(t^2.D)

endfor

;Now we just normalise by the variance to compute the probabilities!!
Pn=Pn/2.D/(var^2.D)
Prob = 1.D - ( 1.D - exp(-Pn) )^n
for i=0,n-1 do begin ;Accomodate any round-off error!!
if Prob[i] LT 0.001 then Prob[i] = n*exp(-Pn[i])
endfor

;---------------------------------
;Now we bootstrap the FAP!!
per=1.D/freq
res=reverse(sort(Pn))
Pn_1=Pn[res] & per_1=per[res] & Prob_1=Prob[res]

if keyword_set(noboot) then goto,skip_boot
if ~keyword_set(noboot) AND keyword_set(fap) then goto,skip_boot

if keyword_set(max_search) then begin
res_test=where(per_1 GE max_search[0] AND per_1 LE max_search[1])
Pn_1[0]=max(Pn_1[res_test])
endif

count=0.
if ~keyword_set(random) then random=10000
random=float(random)

out_power=dblarr(random)

for ii=0,random-1 do begin

rv1=scramble_jenks(rv)
z = rv1 - mean(rv1)
var = stddev(rv1)
Pn1 = fltarr(n)

for i=0,n-1 do begin

w=2.D*!dpi*freq[i]
if w GT 0. then begin

twt = 2.D*w*t
tau = (atan(total(sin(twt)),total(cos(twt)))/2.D)/w
wtmt = w*(t - tau)
pn1[i] = (total(z*cos(wtmt))^2.D)/total(cos(wtmt)^2.D) + $
(total(z*sin(wtmt))^2.D)/total(sin(wtmt)^2.D)

endif
if w LE 0. then pn1[i] = (total(z*t)^2.D)/total(t^2.D)

endfor

;Now we just normalise by the variance to compute the probabilities!!
Pn1=Pn1/2.D/(var^2.D)
Prob1 = 1.D - ( 1.D - exp(-Pn1) )^n
for i=0,n-1 do begin ;Accomodate any round-off error!!
if Prob1[i] LT 0.001 then Prob1[i] = n*exp(-Pn1[i])
endfor

;Now we test if the peak power is higher!!
res=reverse(sort(Pn1))
Pn1=Pn1[res]
if Pn1[0] GE Pn_1[0] then count=count + 1.
out_power[ii]=Pn1[0]

endfor

print,format='(a,e,a)','The bootstrap FAP for the strongest peak is =
',(count / random) * 100.,'%!!'

;Calculate the 1% FAP boundary!!
nums=fix(random * 0.01)
out_power1=out_power[reverse(sort(out_power))]
fap_1p=out_power1[nums-1]
faps=fltarr(n_elements(per))+fap_1p[0]
print,format='(a,d)','The 1% FAP level is = ',fap_1p

;Calculate the 0.1% FAP boundary!!
nums=fix(random * 0.001)
out_power1=out_power[reverse(sort(out_power))]
fap_01p=out_power1[nums-1]
faps1=fltarr(n_elements(per))+fap_01p[0]
print,'The 1% FAP level is =',fap_01p

skip_boot:

if keyword_set(log) then begin
per=alog10(per)
mark=alog10(mark)
endif


if keyword_set(range) then begin
res=where(per_1 GE range[0] AND per_1 LE range[1])
per_1=per_1[res]
Pn=Pn[res]
Prob_1=Prob_1[res]
endif


;To print out we should calculate the primary aliases!!
alias=1.d / (1.d - (1.d/per_1))



;Plot the periodogram!!
if keyword_set(no_plot) then goto,skip
if n_elements(freq) GT 1 then begin
if ~keyword_set(sampling) AND ~keyword_set(hipp) then
device,/landscape,filename='periodograms/'+sub_file+'_'+strcompress(planet,/remove_all)+'_'+strcompress(file_counter,/remove_all)+'.ps'
if keyword_set(sampling) AND ~keyword_set(hipp) then
device,/landscape,filename='periodograms/'+sub_file+'_'+strcompress(planet,/remove_all)+'_'+strcompress(file_counter,/remove_all)+'_sampling.ps'
if ~keyword_set(sampling) AND keyword_set(hipp) then
device,/landscape,filename='periodograms/hipparcos/'+sub_file+'_'+strcompress(planet,/remove_all)+'_'+strcompress(file_counter,/remove_all)+'_hipp.ps'
if keyword_set(sampling) AND keyword_set(hipp) then
device,/landscape,filename='periodograms/hipparcos/'+sub_file+'_'+strcompress(planet,/remove_all)+'_'+strcompress(file_counter,/remove_all)+'_hipp_sampling.ps'

if ~keyword_set(yaxis) then begin
if ~keyword_set(log) then
plot,per,Pn,/xlog,ytitle='Power',xtitle='Orbital Period
(Days)',xs=1,xticks=4 ;,yr=[0.,12.]
if keyword_set(log) then
plot,per,Pn,ytitle='Power',xtitle='log!d10!n
(period/days)',xs=1,xticks=4,xr=[1.,2.5],title='HD162255' ;,yr=[0.,4.]
if keyword_set(mark) then oplot,[mark,mark],[0.,1D9],linestyle=2

if ~keyword_set(noboot) AND ~keyword_set(fap) then oplot,per,faps
if ~keyword_set(noboot) AND ~keyword_set(fap) then
xyouts,[1000.],[fap_1p[0]+0.25],'1% FAP',charsize=1.5
if ~keyword_set(noboot) AND ~keyword_set(fap) then
oplot,per,faps1,linestyle=2
if ~keyword_set(noboot) AND ~keyword_set(fap) then
xyouts,[1000.],[fap_01p[0]+0.25],'0.1% FAP',charsize=1.5
if keyword_set(fap) then oplot,per,fltarr(n_elements(per))+fap[0],linestyle=2
if keyword_set(fap) then xyouts,[1000.],[fap[0]+0.25],'1% FAP',charsize=1.5
endif
if keyword_set(yaxis) then begin
if ~keyword_set(log) then
plot,per,Pn,/xlog,ytitle='Power',xtitle='Orbital Period
(Days)',xs=1,xticks=4,yr=[0.,yaxis],ys=1
if keyword_set(log) then
plot,per,Pn,ytitle='Power',xtitle='log!d10!n
(period/days)',xs=1,xticks=4,xr=[1.,2.5],title='HD162255',yr=[0.,yaxis],ys=1
if keyword_set(mark) then oplot,[mark,mark],[0.,1D9],linestyle=2

if ~keyword_set(noboot) AND ~keyword_set(fap) then oplot,per,faps,linestyle=2
if ~keyword_set(noboot) AND ~keyword_set(fap) then
xyouts,[1000.],[fap_1p[0]+0.25],'1% FAP',charsize=1.5
if keyword_set(fap) then oplot,per,fltarr(n_elements(per))+fap[0],linestyle=2
if keyword_set(fap) then xyouts,[1000.],[fap[0]+0.25],'1% FAP',charsize=1.5
endif

if keyword_set(show_planets) then begin
for jj=0,n_elements(show_planets)-1 do begin
oplot,[show_planets[jj],show_planets[jj]],[0.,1000000.],linestyle=2,thick=6
endfor
endif

device,/close
print,''
if ~keyword_set(sampling) AND ~keyword_set(hipp) then print,'Plot
saved to periodograms/'+sub_file+'_'+strcompress(planet,/remove_all)+'_'+strcompress(file_counter,/remove_all)+'.ps'
if keyword_set(sampling) AND ~keyword_set(hipp) then print,'Plot
saved to periodograms/'+sub_file+'_'+strcompress(planet,/remove_all)+'_'+strcompress(file_counter,/remove_all)+'_sampling.ps'
if ~keyword_set(sampling) AND keyword_set(hipp) then print,'Plot
saved to periodograms/hipparcos/'+sub_file+'_'+strcompress(planet,/remove_all)+'_'+strcompress(file_counter,/remove_all)+'_hipp.ps'
if keyword_set(sampling) AND keyword_set(hipp) then print,'Plot
saved to periodograms/hipparcos/'+sub_file+'_'+strcompress(planet,/remove_all)+'_'+strcompress(file_counter,/remove_all)+'_hipp_sampling.ps'
print,''

;Write out the data too of the top 20 Periods!!
if ~keyword_set(sampling) AND ~keyword_set(hipp) then
openw,1,'periodograms/'+sub_file+'_'+strcompress(planet,/remove_all)+'_'+strcompress(file_counter,/remove_all)+'.dat'
if keyword_set(sampling) AND ~keyword_set(hipp) then
openw,1,'periodograms/'+sub_file+'_'+strcompress(planet,/remove_all)+'_'+strcompress(file_counter,/remove_all)+'_sampling.dat'
if ~keyword_set(sampling) AND keyword_set(hipp) then
openw,1,'periodograms/hipparcos/'+sub_file+'_'+strcompress(planet,/remove_all)+'_'+strcompress(file_counter,/remove_all)+'_hipp.dat'
if keyword_set(sampling) AND keyword_set(hipp) then
openw,1,'periodograms/hipparcos/'+sub_file+'_'+strcompress(planet,/remove_all)+'_'+strcompress(file_counter,/remove_all)+'_hipp_sampling.dat'

;Write out the data!!
if ~keyword_set(sampling) AND ~keyword_set(hipp) then
openw,2,'periodograms/'+sub_file+'_'+strcompress(planet,/remove_all)+'_'+strcompress(file_counter,/remove_all)+'_raw.dat'
if keyword_set(sampling) AND ~keyword_set(hipp) then
openw,2,'periodograms/'+sub_file+'_'+strcompress(planet,/remove_all)+'_'+strcompress(file_counter,/remove_all)+'_sampling_raw.dat'
if ~keyword_set(sampling) AND keyword_set(hipp) then
openw,2,'periodograms/hipparcos/'+sub_file+'_'+strcompress(planet,/remove_all)+'_'+strcompress(file_counter,/remove_all)+'_hipp_raw.dat'
if keyword_set(sampling) AND keyword_set(hipp) then
openw,2,'periodograms/hipparcos/'+sub_file+'_'+strcompress(planet,/remove_all)+'_'+strcompress(file_counter,/remove_all)+'_hipp_sampling_raw.dat'


printf,2,'	Period [days]	Power'
for j=0,n_elements(per)-1 do begin
printf,2,format='(d,a,d)',per[j],' ',Pn[j]
endfor
close,2


;Use to set a range of frequencies to output!!
if keyword_set(range) then begin
res=where(per_1 GE range[0] AND per_1 LE range[1])
per_1=per_1[res] & Pn_1=Pn_1[res] & Prob_1=Prob_1[res]
endif

printf,1,'	Period [days]	Power	Probability Alias[days]'
for j=0,199 do begin
printf,1,format='(d,a,d,a,d,a,d)',per_1[j],' ',Pn_1[j],'
',Prob_1[j],' ',alias[j]
endfor
if ~keyword_set(noboot) AND ~keyword_set(fap) then
printf,1,format='(a,e,a)','The bootstrap FAP for the strongest peak is
= ',(count / random) * 100.,'%!!'
if ~keyword_set(noboot) AND ~keyword_set(fap) then
printf,1,format='(a,d)','The 1% FAP level is = ',fap_1p
close,1
endif

skip:

return,[per_1[0],Pn_1[0],Prob_1[0],alias[0]]

END

This snippet took 0.03 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).