;calculate periodogramm
;(power spectral density estimation)
;input: ft = FFT of some function
;       T = observation period
;output: spectral density function (periodogram)
;       frequency = frequency of periodogramm
function powspec,ft,T=T,frequency=freq
  
                                ;calculate periodogramm
  nel=n_elements(ft)
  n2=nel/2
  p=fltarr(n2+1)
  p(0)=1./nel^2*abs(ft(0))^2
  ik=indgen(n2-1)+1
  p(ik)=1./nel^2*(abs(ft(ik))^2+abs(ft(nel-ik))^2)
  p(n2)=1./nel^2*abs(ft(n2))^2
  
                                ;calculate frequency axis
  freq=(findgen(n2+2)/nel)*T/nel
  
  return,p
end

;procedure to print data to file
pro prnt,file
  
  fpng=file+'.png'
  xyouts,/normal,color=!p.color,0.02,0.02,fpng
  img=tvrd(true=1)
  write_png,fpng,img
  print,'Wrote image to: '+fpng
  stop
end

;solution of FT-exercise
pro solution
  !p.charsize=1.5
                                ;read data
  restore,'gll_data.sav'
  help,/st,data
  
                                ;plot magnetic field versus radial distance
  plot,data.r,data.b,/yst,xtitle='radial distance R_J',ytitle='mag. field [nT]'
  prnt,'orig_data'
  
  nel=n_elements(data.b)        ;# of elements of observation
  T=data.t(nel-1)               ;period of observation
  
                                ;try dumb fft and look at the power
                                ;spectrum
  ff=fft(data.b)
  spec1=powspec(ff,T=T,frequency=f)
  plot,f,spec1,/xtype,xrange=[.001,1],xtitle='frequency',ytitle='power'  
  prnt,'orig_spectrum'
  
                                ;try to detrend data first
  b_lin=data.b-(1-findgen(nel)/(nel))*(data.b(0)-data.b(nel-1))-data.b(nel-1)
  
                                ;plot detrended data
  plot,data.r,b_lin,/yst,xtitle='radial distance R_J',ytitle='mag. field [nT]'
  prnt,'lindetrend_data'
  
  ff=fft(b_lin)                 ;fft of linear detrended data
  spec1=powspec(ff,T=T,frequency=f)
  plot,f,spec1,/xtype,xrange=[.001,1],xtitle='frequency',ytitle='power'
  prnt,'lindetrend_spec'
  
                                ;try detrending with median filter
  b_med=data.b-median(data.b,100)
  plot,data.r,b_med,/yst,xtitle='radial distance R_J',ytitle='mag. field [nT]'
  prnt,'meddetrend_data'
  
  ff=fft(b_med)                 ;fft of median detrended data
  spec1=powspec(ff,T=T,frequency=f)
  plot,f,spec1,/xtype,xrange=[.001,1],xtitle='frequency',ytitle='power'
  prnt,'meddetrend_spec'
  
                                ;use first 500 data points only
  nu=500
  iu=indgen(nu)
  b_lin500=data.b(iu)- $
    (1-findgen(nu)/(nu))*(data.b(iu(0))-data.b(iu(nu-1)))-data.b(iu(nu-1))
  plot,data.r(0:nu-1),b_lin500, $
    /yst,xtitle='radial distance R_J',ytitle='mag. field [nT]'
  prnt,'lindetrend_data_first500'
  
  ff=fft(b_lin500)              ;fft for detrended first 500 elements
  spec1=powspec(ff,T=data.t(nu-1),frequency=f)
  plot,f,spec1,/xtype,xrange=[.001,1],xtitle='frequency',ytitle='power'
  prnt,'lindetrend_spec_first500'
  
                                ;use last 500 data points only
  nu=500
  iu=nel-1-indgen(nu)
  b_lin500=data.b(iu)- $
    (1-findgen(nu)/(nu))*(data.b(iu(0))-data.b(iu(nu-1)))-data.b(iu(nu-1))
  plot,data.r(0:nu-1),b_lin500, $
    /yst,xtitle='radial distance R_J',ytitle='mag. field [nT]'
  prnt,'lindetrend_data_last500'
  
  ff=fft(b_lin500)              ;fft for detrended first 500 elements
  spec1=powspec(ff,T=data.t(nu-1),frequency=f)
  plot,f,spec1,/xtype,xrange=[.001,1],xtitle='frequency',ytitle='power'
  prnt,'lindetrend_spec_last500'
  
  
  
  
end