extract_ephem.pro (Kronos IDL library)

;+
; Contains the extract_ephem procedure
;
; :Author:
;    Laurent Lamy
;
; :History:
;    2008/11/10: Created
;
;    2011/12/16: Last Edit
;-
;
;+
; extract_ephem is a procedure that <behavior desc here>
;
; :Params:
;    jour_debut: in, required, type=sometype
;       A parameter named jour_debut
;    jour_fin:   in, required, type=sometype
;       A parameter named jour_fin
;    t:          in, required, type=sometype
;       A parameter named t
;    lat:        in, required, type=sometype
;       A parameter named lat
;    tl:         in, required, type=sometype
;       A parameter named tl
;    rs:         in, required, type=sometype
;       A parameter named rs
;
; :Keywords:
;    source:     in, optional, type=sometype
;       A keyword named source
;    dmin:       in, optional, type=sometype
;       A keyword named dmin
;    plot:       in, optional, type=sometype
;       A keyword named plot
;    ps:         in, optional, type=sometype
;       A keyword named ps
;-
PRO EXTRACT_EPHEM,jour_debut,jour_fin,t,lat,tl,rs,source=source,dmin=dmin,$
                  plot=plot,ps=ps
; ==============================================================================
; LL (10/2008)
; Version: 1.0
; Sub-routines: make_file_list.pro, read_data_binary.pro, read_src_list.pro,
;               q_vmake.pro, q_rot.pro, q_make.pro, q_prod.pro
; ==============================================================================

if ~keyword_set(source) then source = 'sq'
if ~keyword_set(dmin)   then dmin = 3.

time = dblarr(1)
lat  = fltarr(1)
tl   = fltarr(1)
rs   = fltarr(1)

read_data_src_list,src_list
src_list = (src_list(where(src_list.df_name eq source)))(0)

debut = long(aj_t97(jour_debut))
fin = long(aj_t97(jour_fin))

for courant=debut,fin do begin

  jour_courant = t97_aj(courant)  
  plist = make_file_list(jour_courant,0.,jour_courant,24., $
                         level='ephem',ephem_type=src_list.veph)
  if ptr_valid(plist) then begin
    dataFileVeph = *plist
  ptr_free,plist
  endif else begin
    print,'No Veph data... aborting'
    stop
  endelse
  
  plist = make_file_list(jour_courant,0.,jour_courant,24., $
                         level='ephem',ephem_type=src_list.qeph)
  if ptr_valid(plist) then begin
    dataFileQeph = *plist
    ptr_free,plist
  endif else begin
    print,'No Qeph data... aborting'
  endelse 
  
  read_data_binary,dataFileVeph,veph,level='veph'
  read_data_binary,dataFileQeph,qeph,level='qeph'
  
; Test et interpolation si nombres d'elements differents:
  if n_elements(veph) ne n_elements(qeph) then begin
	qeph2 = replicate({data_qeph},n_elements(veph))
	qeph2.time = veph.time
	for i = 0,3L do begin
	  qeph2.q(i) = interpol(qeph.q(i),qeph.time,veph.time)
	endfor
	qeph = qeph2
  endif  

  ; 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)
  veph.r  = qc0(1:3,*)/src_list.body_rad_value ; en Rp

  lat  = [lat,reform(atan(-veph.r(2)/sqrt(veph.r(1)^2+veph.r(0)^2))/!dtor)]
  rs   = [rs,reform(sqrt(veph.r(0)^2+veph.r(1)^2+veph.r(2)^2))]
  tl   = [tl,reform((atan(-qc0(2,*),-qc0(1,*))+!pi)/!pi*12.)]
  time = [time,veph.time]

endfor

lat = lat(1:*)
tl  = tl(1:*)
rs  = rs(1:*)
time = time(1:*)

w = where(time le 0. or abs(lat) gt 90. or tl lt 0. or tl gt 24.,comp=wc)
if (w(0) ne -1) then begin
  time = time(wc)
  lat  = lat(wc)
  rs   = rs(wc)
  tl   = tl(wc)
endif

t = dindgen((fin-debut+1l)*1440./dmin)*dmin/1440.+debut
lat = interpol(lat,time,t)
tl  = interpol(tl,time,t)
rs  = interpol(rs,time,t)


if (max(t) - min(t)) le 1. then begin
  tt=(t-long(min(t)))*24.
  xtit = 'Hours of DOY '+strmid(t97_aj(min(t)),7,7)
endif else begin
  tt=(t-aj_t97(long(strmid(t97_aj(min(t)),7,4)+'000')))  
  xtit = 'DOY '+strmid(t97_aj(min(t)),7,4)
endelse

if keyword_set(ps) then begin
  !p.multi=[0,1,3]
  titre='plot_ephem'
  set_plot,'ps'
  device,filename=titre+'.ps',/col,/landscape,bits=8
  !p.font=0
endif

if keyword_set(plot) then begin
  !p.multi=[0,1,3]
  plot,tt,lat,/xsty,/ysty,yr=[-90,90],xtit=xtit,ytit='Lat (deg)',charsize=1.8
  plot,tt,tl,/xsty,/ysty,yr=[0,24],xtit=xtit,ytit='LT (h)',charsize=1.8
  plot,tt,rs,/xsty,/ysty,yr=[0,60],xtit=xtit,ytit='r (Rs)',charsize=1.8
  !p.multi=[0,1,1]
endif

if keyword_set(ps) then begin
  device,/close & set_plot,'x'
  spawn,'ps2pdf '+titre+'.ps '+titre+'.pdf'
  spawn,'rm -f '+titre+'.ps'
  !p.multi=[0,1,1]
endif

end