high_res_plot.pro (Kronos IDL library)

;+
; Contains the high_res_plot procedure
;
; :Author:
;    Laurent Lamy
;
; :History:
;    (01/2009,05/2009,03/2012)
;
;    2010/05/18: Created
;
;    2013/07/16: Last Edit
;-
;
;+
; Procédure de création de summary plots haute résolution avec
; des flux calibré S normalisé à 1UA et degré de polarisation circulaire V
;
; :Uses:
;    crossp1, dyn, extract_ephem, load_data_ephem, load_data_n3b, load_data_n3e, q_prod, read_antenna_set, sls3, spdynps
;
; :Params:
;    yyyydddb:    in, required, type=sometype
;       A parameter named yyyydddb
;    hhb:         in, required, type=sometype
;       A parameter named hhb
;    yyyyddde:    in, required, type=sometype
;       A parameter named yyyyddde
;    hhe:         in, required, type=sometype
;       A parameter named hhe
;    image_s:     in, required, type=sometype
;       A parameter named image_s
;    image_v:     in, required, type=sometype
;       A parameter named image_v
;    t:           in, required, type=sometype
;       A parameter named t
;    U:           in, required, type=sometype
;       A parameter named U
;
; :Keywords:
;    source:      in, optional, type=sometype
;       A keyword named source
;    framp:       in, optional, type=sometype
;       A keyword named framp
;    ant_set:     in, optional, type=sometype
;       A keyword named ant_set
;    dsec:        in, optional, type=sometype
;       A keyword named dsec
;    snrmin:      in, optional, type=sometype
;       A keyword named snrmin
;    vmin:        in, optional, type=sometype
;       A keyword named vmin
;    vmax:        in, optional, type=sometype
;       A keyword named vmax
;    betamin:     in, optional, type=sometype
;       A keyword named betamin
;    no_rfi:      in, optional, type=sometype
;       A keyword named no_rfi
;    fminmax:     in, optional, type=sometype
;       A keyword named fminmax
;    smin:        in, optional, type=sometype
;       A keyword named smin
;    smax:        in, optional, type=sometype
;       A keyword named smax
;    add_sls:     in, optional, type=sometype
;       A keyword named add_sls
;    nonorm_1AU:  in, optional, type=sometype
;       A keyword named nonorm_1AU
;    type_plot:   in, optional, type=sometype
;       A keyword named type_plot
;    no_reinterp: in, optional, type=sometype
;       A keyword named no_reinterp
;    titreim:     in, optional, type=sometype
;       A keyword named titreim
;    output_path: in, optional, type=sometype
;       A keyword named output_path
;    no_plot:     in, optional, type=sometype
;       A keyword named no_plot
;-
pro high_res_plot,yyyydddb,hhb,yyyyddde,hhe,image_s,image_v,t,U,source=source,$
                  framp=framp,ant_set=ant_set,dsec=dsec,snrmin=snrmin,vmin=vmin,$
                  vmax=vmax,betamin=betamin,no_rfi=no_rfi,fminmax=fminmax,$
                  smin=smin,smax=smax,add_sls=add_sls,nonorm_1AU=nonorm_1AU,$
                  type_plot=type_plot,no_reinterp=no_reinterp,titreim=titreim,$
                  output_path=output_path,no_plot=no_plot

if ~keyword_set(source)      then source='sq'
if ~keyword_set(ant_set)     then ant_set='2-ant'
if ~keyword_set(dsec)        then dsec = 70.d0
if ~keyword_set(fminmax)     then fminmax=[3,1500]
if ~keyword_set(titreim)     then titreim='Image'
if ~keyword_set(output_path) then output_path=''
if ~keyword_set(type_plot)   then type_plot='SV' ; SV, LT

;-------------------------------------------------------------------------------
; Chargement des données:
;-------------------------------------------------------------------------------

if ant_set eq '2-ant' then load_data_n3e,yyyydddb,long(hhb),yyyyddde,long(hhe+1),n1,n2,n3,verbose=verbose,source=source else $
if ant_set eq '3-ant' then load_data_n3b,yyyydddb,long(hhb),yyyyddde,long(hhe+1),n2,n3,verbose=verbose,source=source

wt=where(n2.t97 ge (aj_t97(yyyydddb)+hhb/24.) and n2.t97 lt (aj_t97(yyyyddde)+hhe/24.))
if wt(0) eq -1 then return
n2=n2(wt) & n3=n3(wt)

if n_elements(n3) eq 1 then return

if keyword_set(betamin) then begin
  load_data_ephem,yyyydddb,long(hhb),yyyyddde,long(hhe+1),veph,qeph,veph='vsat',qeph='qssq',t97=n2.t97
  ; qcenter: quaternion containing the body position (with respect to S/C) in S/C frame
  qcenter = q_vmake(veph.r)
  ; qc0    : quaternion containing the body position (with respect to S/C) in body frame
  qc0     = q_rot(qeph.q,qcenter)
  ; rr0    : distance S/C to body (km)
  rr0      = sqrt(total(qc0(1:3,*)^2.,1))
  ; qr1    : rotational quaternion to put the body direction in (XoZ) plane
  ;          i.e. rotation around Z of angle - azimuth(body)
  qr1     = q_make(reform(-atan(qc0(2,*),qc0(1,*))),[0,0,1])
  ; qr2    : rotational quaternion to put the body direction in the X direction
  ;          i.e. rotation around Y of angle pi/2 - colatitude
  qr2     = q_make(reform(!pi/2-acos(qc0(3,*)/rr0)),[0,1,0])
  ; qrot0  : rotational quaternion that merges the qr2, qr1 and qeph.q
  ;          in the resulting frame:
  ;             - Z is pointing into North
  ;             - X is pointing to the body center  
  ;
  ; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  ;  NB:  in the new frame X pointing to the body center, Z is up along the North direction, 
  ;      and Y toward left.
  ; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  qrot0   = q_prod(q_prod(qr2,qr1),qeph.q)
  ; qcenter: quaternion containing the body position (with respect to S/C) in the new frame
  qcenter2= q_rot(qrot0,qcenter)
  ; veph0  : just to keep a link to the original ephemeris data
  veph0   = veph
  ; veph   : updating the body position in the new frame (veph.r should be rr*[1,0,0])
  veph.r  = qcenter2(1:3,*)

  read_antenna_set,ant,'calDec04_H12',/rad
  dthdip = angular_distance(veph0.r,make_vect_sph(1,ant.dip.al,ant.dip.be))
  dthz = angular_distance(veph0.r,make_vect_sph(1,ant.z.al,ant.z.be))
  dthxp = angular_distance(veph0.r,make_vect_sph(1,ant.xp.al,ant.xp.be))
  dthxm = angular_distance(veph0.r,make_vect_sph(1,ant.xm.al,ant.xm.be))
  ; vv = perp au plan des antennes

  dthdipz = fltarr(n_elements(veph0.r(0))) 
  wwdz = where(n2.ant eq 3)
  if wwdz(0) ne -1 then begin
    vv = crossP1(make_vect_sph(1,ant.dip.al,ant.dip.be),make_vect_sph(1,ant.z.al,ant.z.be),/norm)
    dthdipz(wwdz) = !pi/2.-angular_distance(veph0(wwdz).r,vv) 
  endif
  wwxpz = where(n2.ant mod 10 eq 1)
  if wwxpz(0) ne -1 then begin
    vv = crossP1(make_vect_sph(1,ant.xp.al,ant.xp.be),make_vect_sph(1,ant.z.al,ant.z.be),/norm)
    dthdipz(wwxpz) = !pi/2.-angular_distance(veph0(wwxpz).r,vv) 
  endif
  wwxmz = where(n2.ant mod 10 eq 2)
  if wwxmz(0) ne -1 then begin
    vv = crossP1(make_vect_sph(1,ant.xm.al,ant.xm.be),make_vect_sph(1,ant.z.al,ant.z.be),/norm)
    dthdipz(wwxmz) = !pi/2.-angular_distance(veph0(wwxmz).r,vv) 
  endif
endif

;-------------------------------------------------------------------------------
; Sélection:
;-------------------------------------------------------------------------------

; Suppression freqs parasitées:
if keyword_set(no_rfi) then begin
  f_hf=[(findgen(15)*100+393.750),(findgen(15)*100+400),(findgen(15)*100+406.250),525,$
  findgen(3)*100+643.75,findgen(3)*100+650,findgen(5)*100+656.25,781.25,818.75,881.25,918.75,931.25,$
  981.25,1018.75,1081.25,1118.75,1131.25,1050,1143.75,1150,1181.25,1218.75,$
  1225,1231.25,1343.75,1350,1425,1431.25,1550,1650,1718.75,1725,1743.75,1750,1756.25,1825]
  f_parasites=bytarr(n_elements(n2))
  for i=0,n_elements(f_hf)-1l do f_parasites=f_parasites+(abs(n2.f-f_hf(i)) le 0.1) 
  f_bf=[4.0489,4.2437,4.448,5.62630,5.8971,6.7901,7.11690,8.1946,8.589,21.4779,25.9205,48.8898,51.2426,85.9355,90.0711,$
  98.9490,101.302,103.711,191.070,195.613,200.265,205.027,209.902,214.894,284.906,291.681,305.718,312.988,318.75]  
  for i=0,n_elements(f_bf)-1l do f_parasites=f_parasites+(abs(n2.f-f_bf(i)) le 0.1) 
  wf=where(f_parasites eq 0)
  if wf(0) ne -1 then begin
    n2=n2(wf) 
    n3=n3(wf)
    if keyword_set(betamin) then dthdipz=dthdipz(wf)
  endif
endif

; Sélection gamme de fréquences
wf = where(n2.f ge fminmax(0) and n2.f le fminmax(1))
if wf(0) ne -1 then begin
  n3=n3(wf)
  n2=n2(wf)
  if keyword_set(betamin) then dthdipz=dthdipz(wf)
endif

; Selection d'une paire d'antennes pour le mode 3 antennes:
if ant_set eq '3-ant' then begin
  ww=where(n2.ant eq 11b,compl=wc)
  if wc(0) ne -1 then n3(wc).s=0.
endif

; Sélection SNR:
if keyword_set(snrmin) then begin
  ww=where(n3.sn(0) ge snrmin and n3.sn(1) ge snrmin,compl=wc)
  if wc(0) ne -1 then n3(wc).s=0.
endif

; Sélection V:
if keyword_set(vmin) or keyword_set(vmax) then begin
  ww=where(abs(n3.v) ge vmin and abs(n3.v) le vmax,compl=wc)
  if wc(0) ne -1 then n3(wc).s=0.
endif

;-------------------------------------------------------------------------------
; Création des tableaux:
;-------------------------------------------------------------------------------

tb = aj_t97(yyyydddb)+hhb/24.d0
te = aj_t97(yyyyddde)+hhe/24.d0
nt = long(((te-tb)*86400.d0/dsec)+0.99999)
t = dindgen(nt)*dsec/86400.d0 + tb
indt = long((n2.t97 - tb)/(dsec/86400.d0))>0

n2t=n2
if ~keyword_set(framp) then begin
  fs = n2t(sort(n2t.f)).f
  f = fs(uniq(fs))
endif else f=framp(sort(framp))
nf=n_elements(f)
fmin = [min(f)-1.,(f(0:nf-2)+f(1:nf-1))/2.]
fmax = [(f(0:nf-2)+f(1:nf-1))/2.,max(f)+25.]
if fmin(0)     gt fminmax(0) then fmin(0)     = fminmax(0)-1.
if fmax(nf-1l) lt fminmax(1) then fmax(nf-1l) = fminmax(1)+25.

; Correction du temps de propagation Sat-S/C:
read_data_src_list,src_list
src_list = (src_list(where(src_list.df_name eq source)))(0)
ua = 150.e6/src_list.body_rad_value
rp = src_list.body_rad_value*1e3

extract_ephem,yyyydddb,yyyyddde,teph,lat,tl,rs,source=source
wwt=where(teph ge min(t) and teph le max(t))
ti = n2.t97-interpol(rs(wwt),teph(wwt),n2.t97)*rp/3.e8/86400.

; Remplissage des tableaux:
s  = fltarr(nt,nf)-1. 
v  = fltarr(nt,nf)
ns = fltarr(2,nt,nf)    
if ant_set eq '3-ant' and type_plot eq 'LT' then begin
  lin = fltarr(nt,nf)    
  tot = fltarr(nt,nf)    
endif

for ifreq = 0l,nf-1l do begin
  wf = where(n2.f eq f(ifreq),countf) 
  iff = where(fmin le f(ifreq) and fmax gt f(ifreq)) ; case fréquence tableau final
  if n_elements(iff) ne 1 then stop,'Erreur rampe de freq relle differente de rampe de freq tableau'
  iff = iff(0)
  
  if countf gt 0 then begin
    indt = long((ti(wf) - tb)/(dsec/86400.d0))>0;<(nt-1l)
    for itt = min(indt),max(indt) do begin    
      wt = where(indt eq itt, countt) 
      
      if countt gt 0 then begin       
        ns(0,itt,iff) = total(n2(wf(wt)).dt) ; deltat en ms
        ns(1,itt,iff) = mean(n2(wf(wt)).df)  ; deltaf en kHz 
        if (ns(0,itt,iff) eq 0) or (ns(1,itt,iff) eq 0) then stop,'Temps integration ou bande passante nulle'
        if s(itt,iff) eq -1 then s(itt,iff) = 0. else stop,'Probleme remplissage tableaux'
          s(itt,iff)   = total(n3(wf(wt)).s*n2(wf(wt)).df*n2(wf(wt)).dt) $
                         /(ns(0,itt,iff)*ns(1,itt,iff))      
          v(itt,iff)   = total(n3(wf(wt)).v*n3(wf(wt)).s*n2(wf(wt)).df*n2(wf(wt)).dt) $
                         /((s(itt,iff)>1e-30)*ns(0,itt,iff)*ns(1,itt,iff))
            if ant_set eq '3-ant' and type_plot eq 'LT' then begin
          lin(itt,iff) = total(sqrt(n3(wf(wt)).q^2+n3(wf(wt)).u^2)*n3(wf(wt)).s*n2(wf(wt)).df*n2(wf(wt)).dt) $
                         /((s(itt,iff)>1e-30)*ns(0,itt,iff)*ns(1,itt,iff))
          tot(itt,iff) = total(sqrt(n3(wf(wt)).q^2+n3(wf(wt)).u^2+n3(wf(wt)).v^2)*n3(wf(wt)).s*n2(wf(wt)).df*n2(wf(wt)).dt) $
                         /((s(itt,iff)>1e-30)*ns(0,itt,iff)*ns(1,itt,iff))   
            endif
      endif
    endfor   
  endif  
endfor

; Conversion des flux en W.m^-2.Hz-1 :
Z0 = 377. 
LCaCb=1.68
w = where(s ge 0,count)
if count gt 0 then s(w) = s(w)/(Z0*LCaCb^2)

; Réinterpolation (différents modes):
if ~keyword_set(no_reinterp) then begin
  for i=0l,nt-1l do begin
    wgap = where(s(i,*) ge 0,ng)
    if ng gt 1 then begin
      s(i,*)   = interpol(s(i,wgap),f(wgap),f)
      v(i,*)   = interpol(v(i,wgap),f(wgap),f)
        if ant_set eq '3-ant' and type_plot eq 'LT' then begin
      lin(i,*) = interpol(lin(i,wgap),f(wgap),f)
      tot(i,*) = interpol(tot(i,wgap),f(wgap),f)
        endif
    endif
  endfor
endif

; Flux reçu normé à 1 ua:
if ~keyword_set(nonorm_1AU) then begin
  rs=interpol(rs(wwt),teph(wwt),t)
  r=rebin(reform(rs,nt,1),nt,nf)
  s = (s)*(r^2)/(ua^2) 
endif

; Ajout fc:
;restore,'../Saves/Champs/fcyclo/fcyclo_2008291.sav'
;fc=interpol(fc,tfc,t)
;for i=0l,nt-1l do begin
;  w=where(abs((f-fc(i))) eq min(abs((f-fc(i)))))
;  if w(0) ne -1 then begin
;    v(i,w(0))=1.
;    s(i,w(0))=max(s)
;  endif
;endfor


;-------------------------------------------------------------------------------
; Tracé:
;-------------------------------------------------------------------------------

if keyword_set(autocorr)     then s=autoz
if ~keyword_set(smin)        then smin=1e-24
if ~keyword_set(smax)        then smax=max(s)<1e-19

ws = where(s lt smin) & if ws(0) ne -1 then v(ws) = 0.
v  = v>(-1.)<(1.)
s  = alog10(s>smin<smax)

if ant_set eq '3-ant' and type_plot eq 'LT' then begin
  lin = lin>0.<1.
  tot = tot>0.<1.
endif

if keyword_set(add_sls) then begin
  ; SLS3:
  if (tb ge aj_t97(2004001)) and (te lt aj_t97(2007223)) then begin
    sls3,yyyydddb,yyyyddde,tsls,sls,casslong,psls,driftsls
    sls=interpol(sls,tsls,t)
    ww=where(abs(sls-100.) le 1.)
    wwf=where(f gt 800.,nnf)
    if ww(0) ne -1 then begin
      s(ww,n_elements(f)-nnf:*)=max(s)
      v(ww,n_elements(f)-nnf:*)=1.2
      if ant_set eq '3-ant' and type_plot eq 'LT' then begin
        lin(ww,n_elements(f)-nnf:*)=1.2
        tot(ww,n_elements(f)-nnf:*)=1.2
      endif
    endif
  endif
  
  ; SLS meudon:
  ;restore,'/Volumes/HomeDirs/laurent/Documents/RPWS/Pskr/Saves/Periodes/variations_pskr_2004001_2010065_40-500kHz_LS_96000pts_pts>0.sav'
  restore,'/Volumes/HomeDirs/laurent/Documents/RPWS/Saves/Phases/SKR_periods_phases_200d_2004-001_2010-193.sav'
  ;sls=interpol(sls_rh,t_rh+2556.,t)
  sls=interpol(Phi_N,time+2556.,t)
  ww=where(abs(sls) le 4.)
  wwf=where(f lt 6.,nnf)
    if ww(0) ne -1 then begin
      s(ww,0:nnf)=max(s)
      v(ww,0:nnf)=-1.2
      if ant_set eq '3-ant' and type_plot eq 'LT' then begin
        lin(ww,n_elements(f)-nnf:*)=1.2
        tot(ww,n_elements(f)-nnf:*)=1.2
      endif
  endif
  ;sls=interpol(sls_lh,t_lh+2556.,t)
  sls=interpol(Phi_S,time+2556.,t)
  ww=where(abs(sls) le 4.)
  wwf=where(f lt 6.,nnf)
    if ww(0) ne -1 then begin
      s(ww,0:nnf)=max(s)
      v(ww,0:nnf)=1.2
      if ant_set eq '3-ant' and type_plot eq 'LT' then begin
        lin(ww,n_elements(f)-nnf:*)=1.2
        tot(ww,n_elements(f)-nnf:*)=1.2
      endif
  endif
endif

if keyword_set(betamin) then begin
  beta=interpol(abs(dthdipz)*!radeg,n2.t97,t)
  wbeta=where(beta lt betamin)
  DYN,s,0,0.98,tabmin,tabmax
  if wbeta(0) ne -1 then s(wbeta,*)>=((tabmax-tabmin)*0.1+tabmin)
endif


; Bornes:
if (max(t) - min(t)) le 1. and hhb lt hhe then begin
  tt=(t-long(min(t)))*24.
  titrex = 'DOY '+strmid(t97_aj(min(t)),7,4)+'-'+strmid(t97_aj(min(t)),11,3)+' (h)'
endif else begin
  tt=(t-aj_t97(long(strmid(t97_aj(min(t)),7,4)+'000')))
  titrex = 'DOY '+strmid(t97_aj(min(t)),7,4)
endelse
titref = 'Frequency (kHz)'

; Interpolation:
nf=n_elements(f)*4
U=10^(alog10(min(f))+(findgen(nf)/(nf-1L)*(alog10(max(f))-alog10(min(f)))))
ntm=nt
image_s=fltarr(nt,nf) & image_v=image_s
image_lin=image_s     & image_tot=image_s
for ittm=0l,ntm-1l do begin
  image_s(ittm,*)=interpol(s(ittm,*),f,U)
  image_v(ittm,*)=interpol(v(ittm,*),f,U)
  if ant_set eq '3-ant' and type_plot eq 'LT' then begin
    image_lin(ittm,*)=interpol(lin(ittm,*),f,U)
    image_tot(ittm,*)=interpol(tot(ittm,*),f,U)
  endif
endfor

; Tracé:
if ~keyword_set(no_plot) then begin
  !p.multi=[0,1,2] 
  !p.font=0
  set_plot,'ps'
  device,filename=output_path+titreim+'.ps',/landscape,bits=8,/col
  if type_plot eq 'SV' then begin
    if ~keyword_set(nonorm_1AU) then titre1='Flux density @ 1AU (Log W.m!e-2!n.Hz!e-1!n)' else $
    titre1='Flux density (Log W.m!e-2!n.Hz!e-1!n)'
    titre2='Circular polarization'
    spdynps,image_s,min(tt),max(tt),min(f),max(f),'',titref,titre1,0,0,0,alog10(smin),alog10(smax),1,'.',/log,/col_scale 
    spdynps,image_v,min(tt),max(tt),min(f),max(f),titrex,titref,titre2,0,0,0,-1,1,0,'.',/log
    xyouts,0.913,0.17,'RH',/normal
    xyouts,0.913,0.445,'LH',/normal,col=255
  endif
  if type_plot eq 'LT' then begin
    titre1='Linear polarization'
    titre2='Total polarization'
    spdynps,image_lin,min(tt),max(tt),min(f),max(f),'',titref,titre1,0,0,0,0,1,1,'.',/log,/col_scale 
    spdynps,image_tot,min(tt),max(tt),min(f),max(f),titrex,titref,titre2,0,0,0,0,1,1,'.',/log,/col_scale 
  endif
  device,/close & set_plot,'x'
  print,"Writing "+output_path+titreim+'.pdf'
  spawn,'ps2pdf '+output_path+titreim+'.ps '+output_path+titreim+'.pdf'
  spawn,'rm -f '+output_path+titreim+'.ps'
endif

return
end