;*******************************************************
; Compute power spectrum PSD - Calibration one-sided PSD
;*******************************************************

pro compute_spec,ndata,interv,dataobs,frequency,power

interv = float(interv)
sampling = double((1./(ndata*interv))) 
ng = where(abs(dataobs) gt 0.0,ngood)
facfil = float(ndata)/float(ngood)
avg = total(dataobs(ng))/n_elements(dataobs(ng))
dataobs(ng) = dataobs(ng)-avg
psdm = fft(dataobs,-1)
psd = abs(psdm)
psd2 = psd*psd

n2top = ndata/2
power = dblarr(n2top+1)
power = psd2(0:n2top)*2
power(0) = psd2(0)
power(n2top) = psd2(n2top)
power = power*facfil/sampling
frequency = double(lindgen(n2top+1)*sampling)     

RETURN 
END


;*******************************************************
; Routine to fill the gaps in GOLF with BiSON data
; DS 11/16
;*******************************************************

pro merge_golf,nday=nday,day0=day0,dayshift=dayshift,fnmin=fmin,fmax=fmax,ploteps=ploteps

;*******************************************************
; Inputs: 
; - nday: length in day of sub-series analyzed (default 365)
; - day0: day when the gap starts in the sub series (default 100)
; - dayshift: length in day of zeros added in sub series (default 36)
; - fmin: minimum frequency in microHz to compute mean power (default 1000)
; - fmax: maximum frequency in microHz to compute mean power (default 2000)
; - ploteps: output eps plot (default no eps plot = 0)
;******************************************************* 

if n_elements(nday) eq 0 then nday = 365.
if n_elements(day0) eq 0 then day0 = 100.
if n_elements(dayshift) eq 0 then dayshift = 36.
if n_elements(fmin) eq 0 then fmin = 1000.
if n_elements(fmax) eq 0 then fmax = 2000.
if n_elements(ploteps) eq 0 then ploteps = 0

loadct,39
!p.multi=0

; Load GOLF PM1 and PM2 data
pm1 = readfits('tseries/GOLF_Res_bw_rw_1_960411_140305_F.fits',/silent)
pm2 = readfits('tseries/GOLF_Res_bw_rw_2_960411_140305_F.fits',/silent)
res0 = (pm1+pm2)/2d
dt = 20d
index = lindgen(n_elements(res0))*1d
time0 = index*dt/86400d

; Resample at 40 seconds to match BiSON temporal sampling
index = reform(index,2,n_elements(res0)/2)
res0 = mean(res0(index),dimension=1,/double,/nan)
time0 = mean(time0(index),dimension=1,/double,/nan)
dt = 40d

jd0 = julday(4,11,1996,0,0,0) ; GOLF starting day
time0 = time0+jd0

; Load BiSON data
res1 = readfits('tseries/BiSON_allsites-waverage-fill.fits.gz',/silent)
time1 = transpose(res1(0,*))
res1 = transpose(res1(1,*))

; Select common datasets
ngood = where(time1 ge jd0)
if total(ngood) ne -1 then begin
   res1 = res1(ngood)
   time1 = time1(ngood)
   nz = size(res1)
endif

ngood = where(time1 le max(time0) + (time0(1)-time0(0)))
if total(ngood) ne -1 then begin
   time1 = time1(ngood)
   res1 = res1(ngood)
endif

; Interpolate BiSON on GOLF time grid
res1_interp = interpol(res1,time1,time0)
res1 = res1_interp
time1 = time0

; Analyze one-year series
nseries = floor(n_elements(time0)*dt/86400./nday)
mpw0 = dblarr(nseries)
mpw0z = dblarr(nseries)
mpw0f = dblarr(nseries)
mpw1 = dblarr(nseries)
factorcal = dblarr(nseries)
t0_0 = dblarr(nseries)
for k = 0, nseries-1 do begin
   y0 = res0(k*nday*86400./dt:(k+1)*nday*86400./dt-1)
   y1 = res1(k*nday*86400./dt:(k+1)*nday*86400./dt-1)
   t0 = time0(k*nday*86400./dt:(k+1)*nday*86400./dt-1)
   t1 = time1(k*nday*86400./dt:(k+1)*nday*86400./dt-1)
   t0_0(k) = time0(k*nday*86400./dt)

   ng0=where(y0 ne 0) 
   ng1=where(y1 ne 0)
   factorcal(k) = stddev(y0(ng0))/stddev(y1(ng1)) ; First estimate of calibration factor

   ; Replace zeros by BiSON data
   y0z = y0
   y0f = y0
   day0 = 100.
   y0z(day0*86400./dt:(day0+dayshift)*86400./dt) = 0.
   y0f(day0*86400./dt:(day0+dayshift)*86400./dt) = y1(day0*86400./dt:(day0+dayshift)*86400./dt) * factorcal(k)

   ; Compute associated power spectra
   nz=size(y0)
   compute_spec,nz(1),dt,y0,f0,pw0
   compute_spec,nz(1),dt,y0z,f0,pw0z
   compute_spec,nz(1),dt,y0f,f0,pw0f
   compute_spec,nz(1),dt,y1,f1,pw1
   f0 = f0*1e6
   f1 = f1*1e6

   ; Smoothing of the power spectra
   npt_smooth = 1501
   pw0_sm = smooth(pw0,npt_smooth,/edge_truncate)
   pw0z_sm = smooth(pw0z,npt_smooth,/edge_truncate)
   pw0f_sm = smooth(pw0f,npt_smooth,/edge_truncate)
   pw1_sm = smooth(pw1,npt_smooth,/edge_truncate)

   ; Plot time series
   if ploteps eq 1 then begin
      fname = string(format='("psd_fill_",i3.3,".eps")',k+1)
      nimp,name=fname,/eps,/paper,/portrait
   endif
   !p.multi=[0,2,1]
   plot,t0-time0(0),y0,xtitle='!6Day',ytitle='Residuals',title='GOLF (black) & GOLF_FILLED (red)' 
   oplot,t0-time0(0),y0f,col=250
   oplot,[day0+t0(0)-time0(0),day0+t0(0)-time0(0)],[-20,20],lines=2
   oplot,[day0+dayshift+t0(0)-time0(0),day0+dayshift+t0(0)-time0(0)],[-20,20],lines=2

  ; Plot power spectra
   plot,f0(0:n_elements(f0)-1:100),pw0_sm(0:n_elements(f0)-1:100),/xlog,/ylog,xr=[50,8e3],yr=[20,2e4],xtitle='!6Frequency',ytitle='PSD',xs=1,ys=1
   oplot,f0(0:n_elements(f0)-1:100),pw0f_sm(0:n_elements(f0)-1:100),col=250   
  
   ng = where(f0 ge fmin and f0 le fmax)
   mpw0(k) = mean(pw0(ng))
   mpw0z(k) = mean(pw0z(ng))     
   mpw0f(k) = mean(pw0f(ng))      
   mpw1(k) = mean(pw1(ng))      

   if ploteps eq 1 then fimp

endfor

!p.multi=0
plotsym,0,2,/fill
plot,t0_0+dayshift/2.-time0(0),(mpw0-mpw0f)/mpw0f*100,xtitle='!6Day',ytitle='Relative difference in power (%)',ps=8,xr=[0,8000],yr=[-4,3]
oplot,[0,1e4],[0,0],lines=2
p = poly_fit(t0_0+dayshift-time0(0),(mpw0-mpw0f)/mpw0f*100,1,yfit=yfit)
oplot,t0_0+dayshift/2.-time0(0),yfit,col=250
stop

end