df_flux_plot.pro (Kronos IDL library)

;+
; Contains the df_flux_plot procedure
;
; :Author:
;    Laurent Lamy
;
; :History:
;    2008/03/01: Created
;
;    2007/09/01: Changes
;
;    2011/04/18: Last Edit
;-
;
;+
; Procédure de tracé des cartes radio en échelle de flux (par défaut)
;
; :Uses:
;    crossp1, fill_tab_error, grille_sat, high_res_plot, horizon_radio, load_data_ephem, load_mfl, q_prod, read_antenna_set, 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
;
; :Keywords:
;    dtmin:       in, optional, type=float
;       pas de temps en minutes
;    fmin:        in, optional, type=float
;       sélection en fréquence
;    fmax:        in, optional, type=float
;       sélection en fréquence
;    vmin:        in, optional, type=float
;       sélection en taux de polarisation circulaire
;    vmax:        in, optional, type=float
;       sélection en taux de polarisation circulaire
;    snrmin:      in, optional, type=float
;       sélection en SNR 
;    betamin:     in, optional, type=float
;       sélection sur l'angle direction d'arrivée/plan des antennes
;    no_rfi:      in, optional, type=byte
;       suppresion des fréquences parasites 
;    view:        in, optional, type=string
;       champ de vue angulaire dans le plan d'observation
;       3 valeurs possibles: '45', '90' ou '120' 
;    latmin:      in, optional, type=float
;       latitude minimale pour le tracé polaire 
;    titre:       in, optional, type=string
;       titre du fichier sortant (pdf)
;    output_path: in, optional, type=string
;       chemin pour enregistrement du fichier sortant
;
;    smin:        in, optional, type=float
;       limites en flux (flux moyen, en W/m^2/Hz) correspondant à la convention A définie dans [Lamy, JGR, 2009]
;    smax:        in, optional, type=float
;       limites en flux (flux moyen, en W/m^2/Hz) correspondant à la convention A définie dans [Lamy, JGR, 2009]
;    dB:          in, optional, type=byte
;       convention d'intensité B (produit de flux, en somme de dB normalisés) définie dans [Lamy, JGR, 2009]
;    occurrence:  in, optional, type=byte
;       tracé des occurrences (somme simple pour chaque pixel des temps d'intégration)
;
;    position:    in, optional, type=byte
;       trace la position des points qui interceptent leur iso-fce avec des croix (LH en bleu, RH en rouge)
;    pts_limbe:   in, optional, type=byte
;       trace la position des points au limbe, qui n'interceptent pas leur iso-fce, avec des losanges (si /position)
;    err_bar:     in, optional, type=byte
;       trace les barres d'erreur dans la projection polaire
;    horizon:     in, optional, type=byte
;       trace les horizons radio correspondant à fmin et fmax
;    mfl:         in, optional, type=byte
;       trace deux lignes de champ magnétique typiques dans le plan d'observation (Lat=-75° & TL=TL(cassini)+/-3h
;       N.B: Décalage de 2° entre hémisphère nord et sud
;
;    spdyn_flux:  in, optional, type=byte
;       trace un spectre dynamique en flux calibré normalisé à 1UA (W/m^2/Hz) des données sélectionnées avec les critères 
;       définis plus haut
;    spdyn_polar: in, optional, type=byte
;       trace un spectre dynamique en taux de polarisation circulaire des données sélectionnées
;    high_res:    in, optional, type=byte
;       trace un spectre haute resolution avec les données [fmin=3,fmax=1500,vmin=0,vmax=2.,snrmin=0,smin=1e-24]
;    add_sls:     in, optional, type=byte
;       ajoute la référence SLS3=100° aux spectres 
;    mark_spdyn:  in, optional, type=byte
;       marque sur les spectres l'intervalle de temps-fréquence défini par dtmin et [fmin,fmax]
;-
PRO DF_FLUX_PLOT,yyyydddb,hhb,yyyyddde,hhe,dtmin=dtmin,fmin=fmin,fmax=fmax,$
                 vmin=vmin,vmax=vmax,snrmin=snrmin,betamin=betamin,no_rfi=no_rfi,$
                 view=view,latmin=latmin,titre=titre,output_path=output_path,$
                 smin=smin,smax=smax,$
                 dB=dB,occurrence=occurrence,position=position,$
                 pts_limbe=pts_limbe,err_bar=err_bar,horizon=horizon,mfl=mfl,$
                 spdyn_flux=spdyn_flux,spdyn_polar=spdyn_polar,high_res=high_res,$
                 add_sls=add_sls,mark_spdyn=mark_spdyn

if ~keyword_set(dtmin)       then dtmin=(hhe-hhb)*60.
if ~keyword_set(vmin)        then vmin=0.  else vmin=abs(float(vmin))
if ~keyword_set(vmax)        then vmax=1.1 else vmax=abs(float(vmax))
if vmax lt vmin then begin
  v1=vmax & v2=vmin
  vmin=v1 & vmax=v2
endif
if ~keyword_set(fmin)        then fmin=100. else fmin=float(fmin)
if ~keyword_set(fmax)        then fmax=1000. else fmax=float(fmax)
if ~keyword_set(dtmin)       then dtmin=60.d0 else dtmin=double(dtmin)
if ~keyword_set(snrmin)      then snrmin=20. 
if ~keyword_set(betamin)     then betamin=20.
if ~keyword_set(titre)       then titre='SKR_map'
if ~keyword_set(output_path) then output_path=''
if ~keyword_set(view)        then view='90'
if ~keyword_set(latmin)      then latmin=50.
if ~keyword_set(smin)        then smin=-230. 
if ~keyword_set(smax)        then smax=-215. 
smin1=smin & smax1=smax
smin2=smin & smax2=smax

debut = long(aj_t97(yyyydddb))
fin   = long(aj_t97(yyyyddde))
reso  = 2. ; résolution angulaire en degrés
nerr  = 8. ; nombre de points définissant chaque ellipse d'erreur

source='sq'
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

set_plot,'ps'
device,filename=output_path+titre+'.ps',/color,/landscape,bits=8
!p.font=0

; ------------------------------------------------------------------------------
; Spdyn:
; ------------------------------------------------------------------------------

if keyword_set(spdyn_flux) or keyword_set(spdyn_polar) then begin
  if keyword_set(high_res) then begin
    high_res_plot,yyyydddb,0,yyyyddde,24,ssi,vvi,tt,ff,fminmax=[3,1500],add_sls=add_sls,$
                  vmin=0.,vmax=2.,snrmin=0.,betamin=betamin,/no_rfi,/no_plot,smin=1e-23
  endif else begin
    high_res_plot,yyyydddb,hhb,yyyyddde,hhe,ssi,vvi,tt,ff,fminmax=[3,1500],add_sls=add_sls,$
                  vmin=vmin,vmax=vmax,snrmin=snrmin,betamin=betamin,no_rfi=no_rfi,/no_plot,smin=1e-23
  endelse
endif  

; ------------------------------------------------------------------------------
; Boucle jours:
; ------------------------------------------------------------------------------

for jour=debut,fin do begin

if fin*1440+hhe*60-debut*1440-hhb*60 le dtmin then begin
  jour_ajb=yyyydddb & h1=hhb
  jour_aje=yyyyddde & h2=hhe
endif else begin
  jour_ajb=t97_aj(jour) & jour_aje=jour_ajb
  h1=0.d0 & h2=24.d0
  if jour eq debut then h1=double(hhb)
  if jour eq fin   then h2=double(hhe)
endelse
hh1=long(h1) & hh2=long(h2+0.999)

; ------------------------------------------------------------------------------
; Chargement n2/n3e/loc/index:
; ------------------------------------------------------------------------------

ptrDataFileList = make_file_list(jour_ajb,hh1, jour_aje,hh2, level='n2')
if (ptr_valid(ptrDataFileList) eq 0) then begin
    n_n2DataFile=0
endif else begin
    n2DataFileList=*ptrDataFileList
    n_n2DataFile = n_elements(n2DataFileList)
    ptr_free,ptrDataFileList
endelse

ptrDataFileList = make_file_list(jour_ajb,hh1, jour_aje,hh2, level='n3e', $
source='sq', antset='d')
if (ptr_valid(ptrDataFileList) eq 0) then begin
    n_n3DataFile=0
    print,"# Warning ! No N3e data file selected..."
endif else begin
    n3eDataFileList=*ptrDataFileList
    n_n3eDataFile = n_elements(n3eDataFileList)
    ptr_free,ptrDataFileList
endelse

n2=[{data_n2}] 
n3e=[{data_n3e}] 
index=[{data_index}] 
loc=[{data_loc}]
loc_error=[{data_loc}]

for i=0l,n_n3eDataFile-1l do begin
    while strmid(n3eDataFileList(i),9,/rev) ne strmid(n2DataFileList(i),9,/rev) do begin
        if i eq 0 then n2DataFileList=n2DataFileList(1:*) else $
            n2DataFileList=[n2DataFileList(0:i-1),n2DataFileList(i+1:*)]
        n_n2DataFile = n_n2DataFile-1
    endwhile

  read_data_binary,n2DataFileList(i),n2t,level='n2',valid
  if valid eq 1 then begin
    read_data_binary,n3eDataFileList(i),n3et,level='n3e'
    n3e = [n3e,n3et]
    n2  = [n2,n2t(n3et.num)] 
  
    period_path = strmid(n3eDataFileList(i),0,strlen(n3eDataFileList(i))-21)
    str_aj=strmid(n3eDataFileList(i),9,/rev)
    dataFileIn = period_path+'index/INDEX_'+str_aj
    read_data_binary,dataFileIn,indext,level='index'
    index = [index,indext]
  
    if n_elements(n3et) eq 1 then begin
      loct=[{data_loc}] 
      loc_errort=replicate([{data_loc}],nerr)
    endif else begin
      dataFileIn = period_path+'loc/loc_SPV_'+str_aj
      read_data_binary,dataFileIn,loct,level='loc'  
      dataFileIn = period_path+'loc/loc_err_SPV_n'+string(format='(i2.2)',nerr)+'_'+str_aj
      read_data_binary,dataFileIn,loc_errort,level='loc'
    endelse
    loc = [loc,loct]
    loc_error=[loc_error,loc_errort]
  endif  
endfor

if n_elements(n2) eq 1 then return
n3e=n3e[1:*] & n2=n2[1:*] & index=index[1:*] & loc=loc[1:*]
if n_elements(loc) ne n_elements(n3e) then stop,'*** Warning *** Number of loc and n3e events incorrect'

loc_error=loc_error[1:*]
loc_error=reform(loc_error,nerr,n_elements(n3e))

wt=where(n2.t97 ge (long(n2(0).t97)+h1/24.d0) and n2.t97 lt (long(n2(n_elements(n2)-1l).t97)+h2/24.d0))
if wt(0) ne -1 then begin
  n2=n2(wt) 
  n3e=n3e(wt) 
  loc=loc(wt) 
  index=index(wt)
  loc_error=loc_error(*,wt)
endif else stop,'*** Intervalle horaire vide de donnees! ***'

; ------------------------------------------------------------------------------
; Supression fréquences parasites:
; ------------------------------------------------------------------------------

f_hf=[(findgen(15)*100+393.750),(findgen(15)*100+406.250),525,1143.75,1225,1231.25,1743.75,$
1756.25,643.750,650.000,656.250,1343.75,1425,1718.75,(findgen(11)*100+400)]
f_parasites=bytarr(n_elements(n2))
for i=0,n_elements(f_hf)-1l do f_parasites=f_parasites+(n2.f eq f_hf(i)) 
f_bf=[21.4779,25.9205,101.302,200.265,205.027]
for i=0,n_elements(f_bf)-1l do f_parasites=f_parasites+(n2.f gt (f_bf(i)-0.1) and n2.f lt (f_bf(i)+0.1)) 
wf=where(f_parasites eq 0)

if keyword_set(no_rfi) and wf(0) ne -1 then begin
  n2=n2(wf) 
  n3e=n3e(wf) 
  loc=loc(wf) 
  index=index(wf)
  loc_error=loc_error(*,wf)
endif

foot = loc.foot
temp = where(index.ind eq 0)            ; Indice 0 = flux nul
if temp(0) ne -1 then n3e(temp).s=0. 
foot(1,*) = 90.-foot(1,*)               ; Latitude du pied de la ligne de B
foot(2,*) = (foot(2,*)-180.)*24./360.   ; Temps local du pied de la ligne de B
temp=where(foot(2,*) lt 0) & if temp(0) ne -1 then foot(2,temp)=foot(2,temp)+24. 

loc_error=loc_error.foot
loc_error(1,*,*)=90.-loc_error(1,*,*) 
loc_error(2,*,*)=(loc_error(2,*,*)-180.)*24./360. 
for i=0,n_elements(loc_error(0,*,0))-1l do begin
  temp=where(loc_error(2,i,*) lt 0) 
  if temp(0) ne -1 then loc_error(2,i,temp)=loc_error(2,i,temp)+24.
endfor

; SNR:
sn=10*alog10(n3e.sn>1e-2)

; ------------------------------------------------------------------------------
; Ephémérides S/C:
; ------------------------------------------------------------------------------

rs=index.rr
tl=index.tl
lat = reform(asin(index.zz/index.rr)/!dtor)
xyz=fltarr(3,n_elements(n2.t97))
xyz(0,*)=-rs*cos(lat*!dtor)*cos(tl/24.*2*!pi)
xyz(1,*)=-rs*cos(lat*!dtor)*sin(tl/24.*2*!pi)
xyz(2,*)=rs*sin(lat*!dtor)
dr=reform(sqrt((xyz(0,*)-loc.xyz(0))^2+(xyz(1,*)-loc.xyz(1))^2+(xyz(2,*)-loc.xyz(2))^2))

; ------------------------------------------------------------------------------
; Chargement éphémérides du jour courant:
; ------------------------------------------------------------------------------

load_data_ephem,jour_ajb,h1,jour_aje,h2,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

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

if fin*1440+hhe*60-debut*1440-hhb*60 le dtmin then nn=1 else nn=round((h2-h1)*60./dtmin)

for k=0,nn-1l do begin

  if keyword_set(spdyn_flux) or keyword_set(spdyn_polar) then begin
    !p.multi=[3,1,2]  
    if (max(tt) - min(tt)) le 1. and hhb lt hhe then begin
      ti=(min(tt)-long(min(tt)))*24.  & tf=(max(tt)-long(min(tt)))*24.
      titrex = 'DOY '+strmid(t97_aj(min(tt)),7,4)+'-'+strmid(t97_aj(min(tt)),11,3)+' (h)'
    endif else begin
      ti=(min(tt)-aj_t97(long(strmid(t97_aj(min(tt)),7,4)+'000')))  & tf=(max(tt)-aj_t97(long(strmid(t97_aj(min(tt)),7,4)+'000'))) 
      titrex = 'DOY '+strmid(t97_aj(min(tt)),7,4)
    endelse

    wtt=where(tt ge (jour+(k*dtmin+h1*60.)/1440.) and tt lt (jour+((k+1)*dtmin+h1*60.)/1440.),ntt)
    wff=where(ff ge fmin and ff le fmax,nff)
    sstemp=ssi
    vvtemp=vvi

    if keyword_set(mark_spdyn) and wtt(0) ne -1 and wff(0) ne -1 then begin
      i1=wtt(0) & i2=wtt(ntt-1l)
      j1=wff(0) & j2=wff(nff-1l)
      sstemp((i1-1l)>0:i1,j1:j2)=max(ssi) & sstemp(i2:(i2+1l)<((size(ssi))(1)-1l),j1:j2)=max(ssi)
      sstemp(i1:i2,j1:j1+1l)=max(ssi)     & sstemp(i1:i2,j2-1l:j2)=max(ssi)          
      vvtemp((i1-1l)>0:i1,j1:j2)=1.1      & vvtemp(i2:(i2+1)<((size(vvi))(1)-1l),j1:j2)=1.1
      vvtemp(i1:i2,j1:j1+1l)=1.1          & vvtemp(i1:i2,j2-1l:j2)=1.1    
    endif
   
    if keyword_set(spdyn_flux) and ~keyword_set(spdyn_polar) then begin
      spdynps,sstemp,ti,tf,3,1500,titrex,'Frequency (kHz)','Flux density',0,0,0,0,0.98,/log,$
              posxmin=0.05,posxmax=0.95,posymin=0.62,posymax=1
    endif else if keyword_set(spdyn_polar) and ~keyword_set(spdyn_flux) then begin
      spdynps,vvtemp,ti,tf,3,1500,titrex,'Frequency (kHz)','Circular polarization',0,0,0,0,0.98,0,'.',/log,$
              posxmin=0.05,posxmax=0.95,posymin=0.62,posymax=1
    endif else begin
       !p.multi=[4,1,3]
       spdynps,sstemp,ti,tf,3,1500,'','','',0,0,0,0,0.98,/log,$
              posxmin=0.05,posxmax=0.95,posymin=0.805,posymax=1.
       spdynps,vvtemp,ti,tf,3,1500,'','','',0,0,0,0,0.98,0,'.',/log,$
              posxmin=0.05,posxmax=0.95,posymin=0.6,posymax=0.795
       xyouts,0.9815,0.62,'RH',charsize=.7,/normal
       xyouts,0.9815,0.761,'LH',charsize=.7,/normal,col=255
       xyouts,0.5,0.55,titrex,charsize=1.,/normal,align=0.5
       xyouts,0.,0.8,'Frequency (kHz)',charsize=1.,/normal,orient=90,align=0.5
       xyouts,0.06,0.815,'Flux density',charsize=1.,/normal
       xyouts,0.06,0.61,'Circular polarization',charsize=1.,/normal,col=255
    endelse
    
  endif else !p.multi=[0,2,1]

  ; Tracé dans le champ de vue de Cassini: 
  ww=where(n2.t97 ge (jour+(k*dtmin+h1*60.d0)/1440.) and n2.t97 lt (jour+((k+1)*dtmin+h1*60.d0)/1440.),count)
  if count gt 1 then begin

    ; Points sélectionnés: (wwj = points au limbe, wwi = croisent iso-fce , wwk = wwi+wwj)
    wwi=where(n2.f ge fmin and n2.f le fmax and abs(n3e.v) ge vmin and abs(n3e.v) le vmax and loc.dist eq 0 and foot(0,*) ne 0. and $
             sn(0,*) ge snrmin and sn(1,*) ge snrmin and abs(dthdipz) ge (betamin*!dtor) and $
             n2.t97 ge (jour+(k*dtmin+h1*60.d0)/1440.) and n2.t97 lt (jour+((k+1)*dtmin+h1*60.d0)/1440.),counti)             
    wwj=where(n2.f ge fmin and n2.f le fmax and abs(n3e.v) ge vmin and abs(n3e.v) le vmax and loc.dist gt 0 and foot(0,*) ne 0. and $
             sn(0,*) ge snrmin and sn(1,*) ge snrmin and abs(dthdipz) ge (betamin*!dtor) and $
             n2.t97 ge (jour+(k*dtmin+h1*60.d0)/1440.) and n2.t97 lt (jour+((k+1)*dtmin+h1*60.d0)/1440.),countj)    
    wwk=where(n2.f ge fmin and n2.f le fmax and abs(n3e.v) ge vmin and abs(n3e.v) le vmax and loc.dist ge 0 and foot(0,*) ne 0. and $
             sn(0,*) ge snrmin and sn(1,*) ge snrmin and abs(dthdipz) ge (betamin*!dtor) and $
             n2.t97 ge (jour+(k*dtmin+h1*60.d0)/1440.) and n2.t97 lt (jour+((k+1)*dtmin+h1*60.d0)/1440.),countk)    

    if keyword_set(spdyn_flux) or keyword_set(spdyn_polar) then begin
      !p.multi(0)=2 
      !p.multi(2)=2
    endif
    if view eq '45' then lim=mean(rs(ww))/2. else $
    if view eq '90' then lim=mean(rs(ww)) else $ 
    if view eq '120' then lim=mean(rs(ww))*tan(60.*!dtor) else $
    stop,'Choose field of view 45/90/120'
    xmin=-lim & xmax=lim
    ymin=-lim & ymax=lim
    posxmin=0.0  & posymin=posxmin 
    posxmax=0.5  & posymax=posxmax
    
    plot,[xmin,xmax],[ymin,ymax],xr=[xmin,xmax],yr=[ymin,ymax],xstyle=13,ystyle=13,$
    xtit='',ytit='',tit='',charsize=1.,$
    ticklen=0.2,/nodata,/iso,pos=[posxmin,posymin,posxmax,posymax]
    xyouts,0.185,0.51,'View from Cassini (FOV = '+view+'!eo!n)',charsize=1.3,col=0,/normal,align=0.5

    if countk gt 1 then begin  
      ; Résolution dégradée par une ellipse d'erreur:
      field_view=fltarr(1024,1024)
      field_viewdt=field_view
      fill_tab_error,field_view,index(wwk).ydf,index(wwk).zdf,$
        n3e(wwk).s*n2(wwk).dt*n2(wwk).df*(dr(wwk)/ua)^2/1064.04,reso,nerr,lim,dB=dB
      fill_tab_error,field_viewdt,index(wwk).ydf,index(wwk).zdf,n2(wwk).dt/1e3,reso,nerr,lim,/sum      
      
      ; Repassage aux Watts.m^-2:
      if ~keyword_set(dB) then field_view = field_view/total(n2(wwk).dt*n2(wwk).df)

      ; Dynamic range
      if keyword_Set(dB) then begin 
        smin1=min(field_view) 
        smax1=max(field_view)
        tab=bytscl(field_view)
      endif else begin
        ws=where(field_view gt 0)
        if ~keyword_set(set_smin) and ws(0) ne -1 then smin1=10*alog10(min(field_view(ws))) else smin1=smin
        if ~keyword_set(set_smax) and ws(0) ne -1 then smax1=10*alog10(max(field_view(ws))) else smax1=smax
        tab=bytscl(10*alog10(field_view>1e-40),min=smin1,max=smax1)
      endelse
      
      if keyword_set(occurrence) then begin
        smin1=min(field_viewdt)
        smax1=max(field_viewdt)
        tab=bytscl(field_viewdt,min=smin1,max=smax1)
      endif

      tab=255-tab

      if keyword_set(position) then begin
        wwirh=[where(n3e(wwi).v lt 0.,compl=wwilh)]
        wwjrh=[where(n3e(wwj).v lt 0.,compl=wwjlh)]
        loadct,1
        if wwilh(0) ne -1 then begin
          oplot,[index(wwi(wwilh)).ydf],[index(wwi(wwilh)).zdf],psym=1,symsize=0.5,col=140
          if keyword_set(pts_limbe) and countj gt 0 then $
            oplot,[index(wwj(wwjlh)).ydf],[index(wwj(wwjlh)).zdf],psym=4,symsize=0.5,col=140
        endif
        loadct,3
        if wwirh(0) ne -1 then begin
          oplot,[index(wwi(wwirh)).ydf],[index(wwi(wwirh)).zdf],psym=1,symsize=0.5,col=140
          if keyword_set(pts_limbe) and countj gt 0 then $
            oplot,[index(wwj(wwjrh)).ydf],[index(wwj(wwjrh)).zdf],psym=4,symsize=0.5,col=140
        endif
        loadct,0
      endif else $
        tv,tab,!x.window(0),!y.window(0),xsize=!x.window(1)-!x.window(0),ysize=!y.window(1)-!y.window(0),/normal

    endif

    ; Ligne(s) de champ moyenne(s):
    if keyword_set(mfl) and counti gt 0 then begin
      tl_mfl=[mean(tl(ww))-3,mean(tl(ww))+3]
      ;lat_mfl=[-75,-75] & if mean(lat(ww)) gt 0 then lat_mfl=abs(lat_mfl)+2
      tl_mfl=[5,7,9] & lat_mfl=[-75,-75,-75]
      for mfl=0,n_elements(tl_mfl)-1l do begin
        LOAD_MFL,mean(lat(ww)),mean(tl(ww)),lat_mfl(mfl),tl_mfl(mfl),x,y 
        oplot,x,y,line=2,thick=2,col=150
      endfor  
    endif

    axis, xaxis=0, xra=[xmin,xmax], xstyle=1, xticklen=0.015, $
          xtitle='R!iS!n', charsize = 1.1
    axis, xaxis=1, xra=[xmin,xmax], xstyle=1, xticklen=0.015, $
          xtickname=replicate(' ',9)
    axis, yaxis=0, yra=[ymin,ymax], ystyle=1, yticklen=0.015, $
          ytitle='R!iS!n', charsize = 1.1
    axis, yaxis=1, yra=[ymin,ymax], ystyle=1, yticklen=0.015, $
          ytickname=replicate(' ',9)
   
    ; Erreur typique:
    rdeg=tan(reso*!dtor)*lim
    if keyword_set(position) then $
    oplot,0.9*lim+rdeg*cos(findgen(30)/30.*2.1*!pi),0.9*lim+rdeg*sin(findgen(30)/30.*2.1*!pi),thick=2
 
    ; Grille: (fixe en TL)
    loadct,3
    grille_sat,alpha=0,beta=mean(lat(ww)),1.,0.098,!pi-mean(tl(ww))/24.*2*!pi,$
    /plot,/equ,/mer_cml,/par,/meri
    loadct,0
      
    ;-------------------------------------------------------------------------
   
    ; Tracé, polaire au pied des lignes de champ:
    if keyword_set(spdyn_flux) or keyword_set(spdyn_polar) then !p.multi(0)=1    
    rmax=sqrt(1./(cos(latmin*!dtor)^2.+(sin(latmin*!dtor)/(1-0.098))^2.))*cos(latmin*!dtor)
    xmin=-rmax & xmax=rmax
    ymin=-rmax & ymax=rmax
    posxmin=0.525  & posymin=0.025
    posxmax=0.975  & posymax=0.475
  
    plot,[xmin,xmax],[ymin,ymax],xr=[xmin,xmax],yr=[ymin,ymax],xstyle=13,ystyle=13,$
    xtit='',ytit='',tit='',charsize=1.,$
    ticklen=0.2,/nodata,/iso,pos=[posxmin,posymin,posxmax,posymax]
    xyouts,0.69,0.51,'Magnetic polar projection',charsize=1.3,col=0,/normal,align=0.5
    
    if counti gt 0 then begin
      proj=fltarr(1024,1024)  
      proj_dt=proj
      y=foot(0,wwi)*cos(foot(1,wwi)*!dtor)*cos(foot(2,wwi)/24.*2*!pi)
      x=-foot(0,wwi)*cos(foot(1,wwi)*!dtor)*sin(foot(2,wwi)/24.*2*!pi)
      err=fltarr(2,nerr,counti)
      err(0,*,*)=-loc_error(0,*,wwi)*cos(loc_error(1,*,wwi)*!dtor)*sin(loc_error(2,*,wwi)/24.*2*!pi)
      err(1,*,*)=loc_error(0,*,wwi)*cos(loc_error(1,*,wwi)*!dtor)*cos(loc_error(2,*,wwi)/24.*2*!pi)      
      
      fill_tab_error,proj,x,y,n3e(wwi).s*n2(wwi).dt*n2(wwi).df*(dr(wwi)/ua)^2/1064.04,err,nerr,1.,dB=dB
      fill_tab_error,proj_dt,x,y,n2(wwi).dt/1e3,err,nerr,1.,/sum

      ; Repassage aux Watts (J/s) .m^-2:
      if ~keyword_set(dB) then proj = proj/(total(n2(wwi).dt*n2(wwi).df))

      rmax_pix=round(rmax*512)
      for a=0l,1023l do begin
      for b=0l,1023l do begin
        if sqrt((a-512l)^2+(b-512l)^2) gt rmax_pix then proj(a,b)=0.
        if sqrt((a-512l)^2+(b-512l)^2) gt rmax_pix then proj_dt(a,b)=0.
      endfor
      endfor
      tab=proj(*,512-rmax_pix:511+rmax_pix)       & tab=tab(512-rmax_pix:511+rmax_pix,*)
      tabdt=proj_dt(*,512-rmax_pix:511+rmax_pix)  & tabdt=tabdt(512-rmax_pix:511+rmax_pix,*)

      ; Dynamic range
      if keyword_set(dB) then begin
        smin2=min(tab) 
        smax2=max(tab)
        tab=bytscl(tab,smin2,smax2)
      endif else begin
        ws=where(tab gt 0)
        if ~keyword_set(set_smin) and ws(0) ne -1 then smin2=10*alog10(min(tab(ws))) else smin2=smin-20.
        if ~keyword_set(set_smax) and ws(0) ne -1 then smax2=10*alog10(max(tab(ws))) else smax2=smax-20.
        tab=bytscl(10*alog10(tab>1e-40),min=smin2,max=smax2)
      endelse

      if keyword_set(occurrence) then begin
        smin2=min(tabdt)
        smax2=max(tabdt)
        tab=bytscl(tabdt,min=smin2,max=smax2)
      endif

      tab=255-tab
      
      if keyword_set(position) then begin
        ypol=reform(foot(0,*)*cos(foot(1,*)*!dtor)*cos(foot(2,*)/24.*2*!pi))
        xpol=reform(-foot(0,*)*cos(foot(1,*)*!dtor)*sin(foot(2,*)/24.*2*!pi))
        
        if wwilh(0) ne -1 then begin
          loadct,1 
          oplot,[xpol(wwi(wwilh))],[ypol(wwi(wwilh))],col=140,psym=1
          if keyword_set(pts_limbe) and countj gt 0 then oplot,[xpol(wwj(wwjlh))],[ypol(wwj(wwjlh))],col=140,psym=4
        endif
        if wwirh(0) ne -1 then begin
          loadct,3 
          oplot,[xpol(wwi(wwirh))],[ypol(wwi(wwirh))],col=140,psym=1
          if keyword_set(pts_limbe) and countj gt 0 then oplot,[xpol(wwj(wwjrh))],[ypol(wwj(wwjrh))],col=140,psym=4
        endif        

        if keyword_set(err_bar) then begin
          loadct,1 
          if wwilh(0) ne -1 then begin
            for tmp=0,n_elements(wwilh)-1l do begin
              y2deg=loc_error(0,indgen(nerr+1) mod nerr,wwi(wwilh(tmp)))*cos(loc_error(1,indgen(nerr+1) mod nerr,wwi(wwilh(tmp)))*!dtor)*cos(loc_error(2,indgen(nerr+1) mod nerr,wwi(wwilh(tmp)))/24.*2*!pi)
              x2deg=-loc_error(0,indgen(nerr+1) mod nerr,wwi(wwilh(tmp)))*cos(loc_error(1,indgen(nerr+1) mod nerr,wwi(wwilh(tmp)))*!dtor)*sin(loc_error(2,indgen(nerr+1) mod nerr,wwi(wwilh(tmp)))/24.*2*!pi)
              w2deg=where(x2deg ne 0. and y2deg ne 0.)
              if w2deg(0) gt 2 then begin
                x2deg=[x2deg(w2deg),x2deg(w2deg(0))] 
                y2deg=[y2deg(w2deg),y2deg(w2deg(0))]
              endif
              oplot,x2deg,y2deg,col=150,thick=2
              oplot,[xpol(wwi(wwilh(tmp)))],[ypol(wwi(wwilh(tmp)))],psym=1,col=140
            endfor
          endif
        
          loadct,3
          if wwirh(0) ne -1 then begin
            for tmp=0,n_elements(wwirh)-1l do begin
              y2deg=loc_error(0,indgen(nerr+1) mod nerr,wwi(wwirh(tmp)))*cos(loc_error(1,indgen(nerr+1) mod nerr,wwi(wwirh(tmp)))*!dtor)*cos(loc_error(2,indgen(nerr+1) mod nerr,wwi(wwirh(tmp)))/24.*2*!pi)
              x2deg=-loc_error(0,indgen(nerr+1) mod nerr,wwi(wwirh(tmp)))*cos(loc_error(1,indgen(nerr+1) mod nerr,wwi(wwirh(tmp)))*!dtor)*sin(loc_error(2,indgen(nerr+1) mod nerr,wwi(wwirh(tmp)))/24.*2*!pi)
              w2deg=where(x2deg ne 0. and y2deg ne 0.)
              if w2deg(0) gt 2 then begin
                x2deg=[x2deg(w2deg),x2deg(w2deg(0))] 
                y2deg=[y2deg(w2deg),y2deg(w2deg(0))]
              endif
              oplot,x2deg,y2deg,col=150,thick=2
              oplot,[xpol(wwi(wwirh(tmp)))],[ypol(wwi(wwirh(tmp)))],psym=1,col=140
            endfor
          endif
        endif
        loadct,0
        
      endif else tv,tab,!x.window(0),!y.window(0),xsize=!x.window(1)-!x.window(0),ysize=!y.window(1)-!y.window(0),/normal
    endif
      
    ; Grille:
    for i=0,23 do oplot,[0,rmax],([1,1]*float(i)/12)*!pi,/polar,line=1,col=100,thick=2
    for i=latmin,90,10 do oplot,fltarr(101)+sqrt(1./(cos(i*!dtor)^2.+(sin(i*!dtor)/(1-0.098))^2.))*cos(i/!radeg),findgen(101)/100*2.*!pi,/polar,line=1,col=100,thick=2
    if latmin eq 0. then begin
      for i=latmin,latmin do oplot,fltarr(101)+sqrt(1./(cos(i*!dtor)^2.+(sin(i*!dtor)/(1-0.098))^2.))*cos(i/!radeg),findgen(101)/100*2.*!pi,/polar,col=100,thick=2
    endif
    
    ; TL Cassini:
    loadct,2 & oplot,[0,-rmax*sin(tl(ww(count/2-1l))/24.*2*!pi)],[0,rmax*cos(tl(ww(count/2-1l))/24.*2*!pi)],thick=8,line=3,col=15
    arrow,-0.95*rmax*sin(tl(ww(count/2-1l))/24.*2*!pi),0.95*rmax*cos(tl(ww(count/2-1l))/24.*2*!pi),$
    -rmax*sin(tl(ww(count/2-1l))/24.*2*!pi),rmax*cos(tl(ww(count/2-1l))/24.*2*!pi),/solid,hthick=2,thick=2,col=15,/data
    loadct,0
      
    xyouts,0.69,0.482,'24:00',charsize=1.,/normal,align=0.5
    loadct,3 & xyouts,0.69,0.003,'12:00',charsize=1.,col=150,/normal,align=0.5 & loadct,0
    xyouts,0.52,0.25,'06:00',charsize=1.,/normal,align=0.5,orient=90
    xyouts,0.862,0.25,'18:00',charsize=1.,/normal,align=0.5,orient=-90
    if latmin eq 50. then begin
      xyouts,0.815,0.105,'50!eo!n',charsize=0.8,col=0,/normal
      xyouts,0.75,0.172,'70!eo!n',charsize=0.8,col=0,/normal     
    endif   
      
    ; Horizon radio:
    if keyword_set(horizon) then begin
      loadct,2
      if ~ keyword_set(ntheta) then ntheta=180
      if ~ keyword_set(nphi)   then nphi=360 
      test=horizon_radio([mean(xyz(0,ww)),mean(xyz(1,ww)),mean(xyz(2,ww))],fmin,mfl_model='SPV',horizon_foot=horizon_foot,ntheta=ntheta,nphi=nphi)
      wmfl=where(horizon_foot(1,*) gt 90.,complement=wmflcomp)
      if lat(ww(count/2-1l)) gt 0. then wmfl=wmflcomp
      x_horiz=-horizon_foot(0,wmfl)*cos((90.-horizon_foot(1,wmfl))*!dtor)*sin((horizon_foot(2,wmfl)+180.)*!dtor)
      y_horiz=horizon_foot(0,wmfl)*cos((90.-horizon_foot(1,wmfl))*!dtor)*cos((horizon_foot(2,wmfl)+180.)*!dtor)
      test_horiz=where(sqrt(x_horiz^2+y_horiz^2) le rmax)
      if test_horiz(0) ne -1 then begin
        x_horiz=x_horiz(test_horiz) & y_horiz=y_horiz(test_horiz)
      endif
      w_horiz = where(x_horiz ne 0 and y_horiz ne 0)
      if w_horiz(0) ne -1 then begin
        x_horiz = x_horiz(w_horiz)
        y_horiz = y_horiz(w_horiz)
        lon_ref = ((tl(ww(count/2-1l))/24.*2*!pi+0.5*!pi) mod (2.*!pi))
        x1_horiz=x_horiz*cos(-lon_ref)-y_horiz*sin(-lon_ref)
        y1_horiz=x_horiz*sin(-lon_ref)+y_horiz*cos(-lon_ref)
        x0_horiz=double(max(x1_horiz)+min(x1_horiz))/2.
        y0_horiz=double(max(y1_horiz)+min(y1_horiz))/2.
        test_horiz=atan((y1_horiz-y0_horiz),(x1_horiz-x0_horiz))
        wxy=where(test_horiz lt 0) & if wxy(0) ne -1 then test_horiz(wxy)+=2.*!pi
        s_horiz=sort(test_horiz)
        oplot,x_horiz(s_horiz),y_horiz(s_horiz),thick=6,col=10          
      endif
      if fmax ne fmin then begin
        test=horizon_radio([mean(xyz(0,ww)),mean(xyz(1,ww)),mean(xyz(2,ww))],fmax,mfl_model='SPV',horizon_foot=horizon_foot,ntheta=ntheta,nphi=nphi)
        wmfl=where(horizon_foot(1,*) gt 90.,complement=wmflcomp)
        if lat(ww(count/2-1l)) gt 0. then wmfl=wmflcomp
        x_horiz=-horizon_foot(0,wmfl)*cos((90.-horizon_foot(1,wmfl))*!dtor)*sin((horizon_foot(2,wmfl)+180.)*!dtor)
        y_horiz=horizon_foot(0,wmfl)*cos((90.-horizon_foot(1,wmfl))*!dtor)*cos((horizon_foot(2,wmfl)+180.)*!dtor)
        test_horiz=where(sqrt(x_horiz^2+y_horiz^2) le rmax)
        if test_horiz(0) ne -1 then begin
          x_horiz=x_horiz(test_horiz) & y_horiz=y_horiz(test_horiz)
        endif
        w_horiz = where(x_horiz ne 0 and y_horiz ne 0)
        if w_horiz(0) ne -1 then begin
          x_horiz = x_horiz(w_horiz)
          y_horiz = y_horiz(w_horiz)
          lon_ref = ((tl(ww(count/2-1l))/24.*2*!pi+0.5*!pi) mod (2.*!pi))
          x1_horiz=x_horiz*cos(-lon_ref)-y_horiz*sin(-lon_ref)
          y1_horiz=x_horiz*sin(-lon_ref)+y_horiz*cos(-lon_ref)
          x0_horiz=double(max(x1_horiz)+min(x1_horiz))/2.
          y0_horiz=double(max(y1_horiz)+min(y1_horiz))/2.
          test_horiz=atan((y1_horiz-y0_horiz),(x1_horiz-x0_horiz))
          wxy=where(test_horiz lt 0) & if wxy(0) ne -1 then test_horiz(wxy)+=2.*!pi
          s_horiz=sort(test_horiz)
          oplot,x_horiz(s_horiz),y_horiz(s_horiz),thick=6,col=20          
        endif
      endif  
      loadct,0
    endif
    
        ; Echelle d'intensité graphe 1:
    if ~keyword_set(position) then begin
      posxmin=0.38 & posymin=0
      posxmax=0.4  & posymax=0.25
      bar_titre=''
      if keyword_set(dB) then bar_tick = [' ','!9S!3 dB!Inorm!N',' '] else bar_tick = [' ','dB',' ']
      plot, [0,1], [0,1], /nodata, /noerase, xst=1, yst=12, $
      pos = [posxmin,posymin,posxmax,posymax], $
      /norm, xticklen=0, xticks=2, xtickv=[0,0.5,1], $
      xtickname=bar_tick, xtitle=bar_titre, charsize=0.8
      axis, /yaxis, yra=[smin1,smax1], ticklen=-0.02, /ysty, charsize=0.8
      bcb = 255-bytscl(replicate(1,10)#bindgen(256))
      tv,bcb,!x.window(0),!y.window(0),xsize = !x.window(1)-!x.window(0), $
      ysize = !y.window(1) - !y.window(0), /normal
      oplot, [1,0,0,1], [1,1,0,0]
   
     ; Echelle d'intensité graphe 2:
      posxmin=0.9     & posymin=0
      posxmax=0.92    & posymax=0.25
      bar_titre=''
      if keyword_set(dB) then bar_tick = [' ','!9S!3 dB!Inorm!N',' '] else bar_tick = [' ','dB',' ']
      plot, [0,1], [0,1], /nodata, /noerase, xst=1, yst=12, $
      pos = [posxmin,posymin,posxmax,posymax], $
      /norm, xticklen=0, xticks=2, xtickv=[0,0.5,1], $
      xtickname=bar_tick, xtitle=bar_titre, charsize=0.8
      axis, /yaxis, yra=[smin2,smax2], ticklen=-0.02, /ysty, charsize=0.8
      bcb = 255-bytscl(replicate(1,10)#bindgen(256))
      tv,bcb,!x.window(0),!y.window(0),xsize = !x.window(1)-!x.window(0), $
      ysize = !y.window(1) - !y.window(0), /normal
      oplot, [1,0,0,1], [1,1,0,0]
    endif 

    if keyword_set(position) then begin
      loadct,3
      posxmin=0.9     & posymin=0.1
      posxmax=0.92    & posymax=0.25
      plot, [0,1], [0,1], /nodata, /noerase, xst=1, yst=12, $
      pos = [posxmin,posymin,posxmax,posymax], $
      /norm, xticklen=0, xticks=2, xtickv=[0,0.5,1], $
      xtickname=[' ','f (kHz)',' '], xtitle='RH', charsize=0.8
      axis, /yaxis, yra=[fmin,fmax], ticklen=-0.02, /ysty, charsize=0.8
      bcb = bytscl(replicate(1,10)#(bindgen(256)/255.*150+70)) 
      tv,bcb,!x.window(0),!y.window(0),xsize = !x.window(1)-!x.window(0), $
      ysize = !y.window(1) - !y.window(0), /normal
      oplot, [1,0,0,1], [1,1,0,0]
      loadct,1
      posxmin=0.97     & posymin=0.1
      posxmax=0.99    & posymax=0.25
      plot, [0,1], [0,1], /nodata, /noerase, xst=1, yst=12, $
      pos = [posxmin,posymin,posxmax,posymax], $
      /norm, xticklen=0, xticks=2, xtickv=[0,0.5,1], $
      xtickname=[' ','f (kHz)',' '], xtitle='LH', charsize=0.8
      axis, /yaxis, yra=[fmin,fmax], ticklen=-0.02, /ysty, charsize=0.8
      bcb = bytscl(replicate(1,10)#(bindgen(256)/255.*150+70)) 
      tv,bcb,!x.window(0),!y.window(0),xsize = !x.window(1)-!x.window(0), $
      ysize = !y.window(1) - !y.window(0), /normal
      oplot, [1,0,0,1], [1,1,0,0]
    endif
    
    ; Ephems Cassini:
    xyouts,0.38,0.48,'Ephemeris:',charsize=1.1,col=0,/normal
    xyouts,0.38,0.478,'__________',charsize=1.,col=0,/normal
    aaaa=strmid(string(t97_aj(n2(ww(count/2-1l)).t97)),7,4)
    jjj=strmid(string(t97_aj(n2(ww(count/2-1l)).t97)),11,3)
    hh=strtrim(string(long(h1+k*dtmin/60.)),1)          & if (h1+k*dtmin/60.) lt 10. then hh='0'+hh
    mm=strtrim(string(long(h1*60.+k*dtmin-hh*60.)),1)   & if float(mm)        lt 10. then mm='0'+mm
    
    ;if k eq 0 then mm1=mm else if mm1 eq mm then stop
    ;mm1=mm
    xyouts,0.375,0.44,'Day',charsize=1.1,col=0,/normal
    xyouts,0.425,0.44,': '+aaaa+'-'+jjj+'   ',charsize=1.1,col=0,/normal
    xyouts,0.375,0.41,'Time',charsize=1.1,col=0,/normal
    xyouts,0.425,0.41,': '+hh+':'+mm,charsize=1.1,col=0,/normal   
    
    hhtl=strmid(strtrim(string(tl(ww(count/2-1l))),1),0,2) & if hhtl lt 10. then hhtl='0'+strmid(hhtl,0,1)
    mmtl=strtrim(string(long(60*(tl(ww(count/2-1l))-long(tl(ww(count/2-1l)))))),1) & if mmtl lt 10. then mmtl='0'+strmid(mmtl,0,1)
    xyouts,0.375,0.37,'r!iS/C!n (R!iS!n)',charsize=1.1,col=0,/normal 
    xyouts,0.435,0.37,'= '+strmid(strtrim(string(rs(ww(count/2-1l))),1),0,4),charsize=1.1,col=0,/normal 
    xyouts,0.375,0.34,'!9l!3!iS/C!n (!eo!n)',charsize=1.1,col=0,/normal
    xyouts,0.435,0.34,'= '+strmid(strtrim(string(lat(ww(count/2-1l))),1),0,5),charsize=1.1,col=0,/normal
    loadct,2 & xyouts,0.375,0.31,'TL!iS/C!n',charsize=1.,col=15,/normal 
    xyouts,0.435,0.31,'= '+hhtl+':'+mmtl,charsize=1.1,col=15,/normal & loadct,0

    ; Paramètres:
    xyouts,0.9,0.48,'Data selection:',charsize=1.1,/normal
    xyouts,0.9,0.478,'_____________',charsize=1.,/normal
    xyouts,0.9,0.44,'!9d!3t (min)',charsize=1.,/normal
    xyouts,0.97,0.44,'= '+strmid(dtmin,7,2),charsize=1.,/normal
    if keyword_set(horizon) then begin
      loadct,2
      xyouts,0.9,0.41,'f!imin!n (kHz)',charsize=1.,col=10,/normal
      xyouts,0.97,0.41,'= '+strtrim(strmid(long(fmin),8,4),1),charsize=1.,col=10,/normal
      xyouts,0.9,0.38,'f!imax!n (kHz)',charsize=1.,col=20,/normal
      xyouts,0.97,0.38,'= '+strtrim(strmid(long(fmax),8,4),1),charsize=1.,col=20,/normal
      loadct,0 
    endif else begin
      xyouts,0.9,0.41,'f!imin!n (kHz)',charsize=1.,/normal
      xyouts,0.97,0.41,'= '+strtrim(strmid(long(fmin),8,4),1),charsize=1.,col=10,/normal
      xyouts,0.9,0.38,'f!imax!n (kHz)',charsize=1.,/normal
      xyouts,0.97,0.38,'= '+strtrim(strmid(long(fmax),8,4),1),charsize=1.,col=20,/normal      
    endelse
    xyouts,0.9,0.35,'V',charsize=1.,/normal
    xyouts,0.97,0.35,'= ['+strmid(strtrim(string(vmin),1),0,4)+','+strmid(strtrim(string(vmax),1),0,4)+']',charsize=1.,/normal
    xyouts,0.9,0.32,'SN!iX,Z!n (dB)',charsize=1.,/normal
    !p.font=-1 & xyouts,0.98,0.32,'!9'+string("142b)+'!3',charsize=1.,/normal & !p.font=0
    xyouts,0.995,0.32,strmid(strtrim(string(snrmin),1),0,2),charsize=1.,/normal
    xyouts,0.9,0.29,'!9b!3 (!eo!n)',charsize=1.,/normal
    !p.font=-1 & xyouts,0.98,0.29,'!9'+string("142b)+'!3',charsize=1.,/normal & !p.font=0
    xyouts,0.995,0.29,strmid(strtrim(string(betamin),1),0,2),charsize=1.,/normal
    ;xyouts,1.02,0.29,'('+strmid(strtrim(string(abs(fix(mean(dthdipz(ww)))*!radeg)),1),0,2)+')',charsize=1.,/normal
    
  endif else begin
    if ~keyword_set(xmin) then xmin=-1.
    if ~keyword_set(xmax) then xmax=+1.
    ymin=xmin & ymax=xmax
    if keyword_set(spdyn_flux) or keyword_set(spdyn_polar) then begin
      !p.multi(0)=2 
      !p.multi(2)=2
    endif
    plot,[0],[0],xr=[xmin,xmax],yr=[ymin,ymax],xstyle=13,ystyle=13,$
    xtit='R!iS!n',ytit='R!iS!n',tit='',charsize=1.1,ticklen=0.015,/iso,pos=[0,0,0.5,0.5]   
    xyouts,0.185,0.51,'View from Cassini (FOV = '+view+'!eo!n)',charsize=1.3,col=0,/normal,align=0.5
    xyouts,0.19,0.25,'No data',charsize=1.3,col=0,/normal,align=0.5
    if keyword_set(spdyn_flux) or keyword_set(spdyn_polar) then !p.multi(0)=1    
    rmax=sqrt(1./(cos(latmin*!dtor)^2.+(sin(latmin*!dtor)/(1-0.098))^2.))*cos(latmin*!dtor)
    xmin=-rmax & xmax=rmax
    ymin=-rmax & ymax=rmax
    plot,[xmin,xmax],[ymin,ymax],xr=[xmin,xmax],yr=[ymin,ymax],xstyle=13,ystyle=13,$
    xtit='',ytit='',tit='',charsize=1.,ticklen=0.2,/nodata,/iso,pos=[0.525,0.025,0.975,0.475]
    xyouts,0.69,0.51,'Magnetic polar projection',charsize=1.3,col=0,/normal,align=0.5
    for i=0,23,1 do oplot,[0,rmax],([1,1]*float(i)/12)*!pi,/polar,line=1,col=100,thick=2
    for i=latmin,90,10 do oplot,fltarr(101)+sqrt(1./(cos(i*!dtor)^2.+(sin(i*!dtor)/(1-0.098))^2.))*cos(i/!radeg),findgen(101)/100*2.*!pi,/polar,line=1,col=100,thick=2
    if latmin eq 0. then begin
      for i=latmin,latmin do oplot,fltarr(101)+sqrt(1./(cos(i*!dtor)^2.+(sin(i*!dtor)/(1-0.098))^2.))*cos(i/!radeg),findgen(101)/100*2.*!pi,/polar,col=100,thick=2
    endif 
    xyouts,0.69,0.482,'24:00',charsize=1.,col=150,/normal,align=0.5
    loadct,3 & xyouts,0.69,0.003,'12:00',charsize=1.,col=150,/normal,align=0.5 & loadct,0
    xyouts,0.52,0.25,'06:00',charsize=1.,col=150,/normal,align=0.5,orient=90
    xyouts,0.862,0.25,'18:00',charsize=1.,col=150,/normal,align=0.5,orient=-90
    if latmin eq 50. then begin
      xyouts,0.815,0.105,'50!eo!n',charsize=0.8,col=0,/normal
      xyouts,0.75,0.172,'70!eo!n',charsize=0.8,col=0,/normal     
    endif
    xyouts,0.38,0.48,'Ephemeris:',charsize=1.1,col=0,/normal
    xyouts,0.38,0.478,'__________',charsize=1.,col=0,/normal
    aaaa=strmid(string(jour_aj),7,4)
    jjj=strmid(string(jour_aj),11,3)
    hh=strtrim(string(long(h1+k*dtmin/60.)),1)          & if (h1+k*dtmin/60.) lt 10. then hh='0'+hh
    mm=strtrim(string(long(h1*60.+k*dtmin-hh*60.)),1)   & if (h1*60.+k*dtmin-hh*60.) lt 10. then mm='0'+mm
    xyouts,0.375,0.44,'Day',charsize=1.1,col=0,/normal
    xyouts,0.425,0.44,': '+aaaa+'-'+jjj+'   ',charsize=1.1,col=0,/normal
    xyouts,0.375,0.41,'Time',charsize=1.1,col=0,/normal
    xyouts,0.425,0.41,': '+hh+':'+mm,charsize=1.1,col=0,/normal
    xyouts,0.9,0.48,'Data selection:',charsize=1.1,col=150,/normal
    xyouts,0.9,0.478,'_____________',charsize=1.,col=150,/normal
    xyouts,0.9,0.44,'No data',charsize=1.,col=150,/normal
  endelse
endfor
endfor

device,/close & set_plot,'x'
spawn,'ps2pdf '+titre+'.ps '+titre+'.pdf'
spawn,'rm -f '+titre+'.ps'


end