;******************************************************* ; 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