function b_dipol,lat,l_shell=l_shell
  common adiabatic,kev,par_dipol,r_planet
  
  return,par_dipol/(r_planet*l_shell*(cos(lat)^2))^3*sqrt(4-3*cos(lat)^2)
end

function t_gyro,energy=e,mass=m,l_shell=l,pa=pa,charge=charge,r_gyro=r_gyro, $
                bf=bf
  common adiabatic,kev,par_dipol,r_planet
  
  if keyword_set(charge) eq 0 then charge=1
  v=calc_v(e,m)
  mred=!amu*m/sqrt(1-v^2/!c_light^2)
  bf=b_dipol(0.,l_shell=l)
  omega_gyro=!e_charge*charge*bf/mred
  t_gyro=2*!pi/omega_gyro
  if keyword_set(pa) or keyword_set(r_gyro) then $
    r_gyro=abs(mred*v*sin(pa)/bf/(!e_charge*charge))
  return,t_gyro
end

function psd,ft
  
  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
  return,p
end
  
  

pro make_data
  common adiabatic,kev,par_dipol,r_planet
  
  physical_constants
  r_planet=!RJ
  par_dipol=!PAR_DIPOL_JUP
  kev=1d3*!e_charge
  userlct,/nologo & erase
  
  lat=0.
  l_shell=[9.,9.00005]
  l_shell=[8.5,9.5]
  speed=100000.
  
  dkm=(l_shell(1)-l_shell(0))*r_planet
  sec=dkm/speed
;  sec=1000.
  res=.33
  nel=fix(abs(sec)/res/2)*2
  tsec=findgen(nel)/(nel-1)*sec
  T=max(tsec)
  
  lvec=findgen(nel)/(nel+1)*(l_shell(1)-l_shell(0))+l_shell(0)
  bdip=b_dipol(lat,l_shell=lvec)*1e9
  
                                ;add compression of magnetosphere
  cpr=-20*exp(-(lvec-total(l_shell)/2.)^2/.1^2)
  bdip=bdip+cpr
  
  mass=[64.,32.,16.,1]
  ampl=[2.,6.,3.,3.]
  nm=n_elements(mass)
  fg=fltarr(nm,nel)
  bsig=bdip
  for i=0,nm-1 do begin
    tg=t_gyro(energy=20.,mass=mass(i),pa=90,charge=1,l_shell=lvec,r_gyro=rg)
    freq=1./tg
    bsig=bsig+ampl(i)*sin(2*!pi*freq*tsec+randomu(s)*2*!pi)
print,mass(i),tg(0)    
  endfor
  
  fftb=fft(bsig-bdip)
  x=tsec
  dx=x(1)-x(0)
  n21=nel/2+1
  f=findgen(nel)
  f[n21]=n21-nel+findgen(n21-2)
  f=f/nel/dx
  sf=shift(f,-n21)
  sffy=shift(fftb,-n21)
  
                                ;convolute signal with response
                                ;function
  fft_rf=fft(exp(-f^2/.000001))
  fft_rf=fft_rf/max(fft_rf)
  fftb_rf=fftb*fft_rf
  b_rf=fft(fftb_rf,1)+bdip
  plot,f,fft_rf & wait,1 &  plot,bsig-bsig(0) & oplot,b_rf-bsig(0),color=1
  
                                ;add noise spectrum
  fft_noise=fltarr(nel)+complex(1,1)
  noise=fft(fft_noise,1)
  
  noise=(randomu(asd,nel)*2-1)*5.
  b_rfn=b_rf+noise
  plot,bsig-bsig(0) & oplot,b_rfn-bsig(0),color=1
  
  
  fftb_rfn=fft(b_rfn-bdip)
  sfftb_rfn=shift(fftb_rfn,-n21)
  
  wait,1
  psd=psd(fftb_rfn)
  np=n_elements(psd)
  fpsd=findgen(np)/nel/dx
  plot,fpsd,psd,xrange=[.02,10],/xtype
  
  data={r:lvec,t:tsec,b:b_rfn}
  save,file='gll_data.sav',data
end