ephemeris.pro (Kronos IDL library)

;+
; Contains the ephemeris procedure
;
; :Author:
;    ??
;
; :History:
;    2007/01/04: Created
;
;    2007/01/04: Last Edit
;-
;
;+
; Ephemeris : Procédure de visualisation des éphémérides de Cassini depuis 
; Saturne dans différents repères
;
; :Params:
;    jour_debut: in, required, type=int
;       aaaajjj
;    jour_fin:   in, required, type=int
;       aaaajjj
;    source:     in, required, type=string
;       'sc' pour écliptique / 'sq' pour équatorial
;    xmax:       in, required, type=sometype
;       calibration des cadrans de visualisation des coordonnees en Rsat
;
; :Keywords:
;    lat:        in, optional, type=byte
;       pour tracer la latitude / t2004 dans la dernier cadran, sinon tracé par 
;       défaut de la distance à saturne
;    ps:         in, optional, type=byte
;       création image.ps
;-
PRO EPHEMERIS,jour_debut,jour_fin,source,xmax,lat=lat,ps=ps

; ==============================================================================
;    (Planetary Solar Equatorial, Planet centered Z=Omega, V2=Planet-Sun vector)
;    X      - In the Planet->Sun plane (at noon), positive towards the Sun.
;    Y      - Z x X
;    Z      - Northward spin axis of the Planet.
; ==============================================================================
;    (Planetary Solar Ecliptic, Planet centered X=Sun, V2=Orbital normal)
;    X      - Along the Planet->Sun line (S), positive towards the Sun.
;    Y      - V2 x X
;    Z      - Parallel to the Planetary orbital plane upward normal.
; ==============================================================================
; Dernière modif par LL le 04/01/07
; ==============================================================================


!p.charsize=1
!p.multi=[0,2,2] & !p.font=0
device,decompose=0

if keyword_set(ps) then begin
  set_plot,'ps'
  device,filename='image.ps',/landscape,bits=8
endif else window,1,xs=600,ys=600

; ------------------------------------------------------------------------------
; Paramètres:
; ------------------------------------------------------------------------------

hd=0.
hf=24.

xmin = -xmax 
ymax = xmax & ymin = xmin 
zmax = xmax & zmin = xmin

; ------------------------------------------------------------------------------
; Titre:
; ------------------------------------------------------------------------------

debut = strtrim(string(jour_debut),1) 
fin = strtrim(string(jour_fin),1)
if strmid(debut,0,4) ne strmid(fin,0,4) then begin
  titre =strmid(debut,0,4)+strmid(debut,4,3)+'_'+strmid(fin,0,4)+strmid(fin,4,3)+' '+source
endif else begin
  titre =strmid(debut,0,4)+'_'+strmid(debut,4,3)+'_'+strmid(fin,4,3)+' '+source
endelse

; ------------------------------------------------------------------------------
; Définition d'une structure type veph pour la concaténation finale:
; ------------------------------------------------------------------------------

ephem={data_veph}
time = fltarr(1)
read_data_src_list,src_list ;,file='/c/rpws/pro/src_list.txt'
src_list = (src_list(where(src_list.df_name eq source)))(0)

; ------------------------------------------------------------------------------
; Boucle de chargement journalier:
; ------------------------------------------------------------------------------

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

for courant=debut,fin do begin
  
; Chargement des fichiers d'ephemerides du jour courant:
; ------------------------------------------------------
  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
  
; Restriction à la période sélectionnée:
; --------------------------------------

  date_debut = aj_t97(double(jour_debut)+0.1*long(hd*10./24.))
  date_fin   = aj_t97(double(jour_fin)+0.1*long(hf*10./24.))
  where_t = where(veph.time ge date_debut and veph.time le date_fin,count)
  if count gt 0 then begin
    veph = veph(where_t)
    qeph = qeph(where_t)

; Calcul des éphémérides:
; -----------------------
  
    ; qpos: quaternion containing the S/C position (with respect to body) in S/C frame
    qpos = q_vmake(-veph.r)
    ; temp: quaternion containing the S/C position (with respect to body) in body frame
    temp = q_rot(qeph.q,qpos)
    veph.r = temp(1:3,*)/src_list.body_rad_value ;(Donne les coordonnées en Rp)
    ;veph.r=temp(1:3,*)/1e5                      ;(Donne les coordonnées en millions de km)
    ephem = [ephem,veph]
    time = [time,veph.time]
  endif  
endfor

; Ephémérides finales:
ephem = ephem(1:*)
pos = ephem.r
time = time(1:*)

; Distance Cassini/Saturne:
rr= sqrt(total(pos(0:2,*)^2.,1))

; Latitude S/C:
if source eq 'sq' then lat = reform(atan(pos(2,*)/sqrt(pos(1,*)^2+pos(0,*)^2))/!dtor)
if source eq 'sc' and keyword_set(lat) then stop,$
'Attention! La latitude doit etre calculee en equatorial, passez en source = sq'

; ------------------------------------------------------------------------------
; Trace final 1:
; ------------------------------------------------------------------------------

plot,pos(0,*),pos(1,*),psym=3,col=0,back=255,xran=[xmin,xmax],yran=[ymin,ymax],$
tit='Cassini ephemeris: '+titre,xtit='X (Rsat)',ytit='Y';,/iso

; Trace additionnel des heures locales:
;--------------------------------------

p=findgen(10)*20
for k=0,23 do begin
  oplot,p*cos(15*k*!dtor),p*sin(15*k*!dtor),col=100,linestyle=2
endfor

; Trace additionnel de Saturne:
;------------------------------

beta=findgen(360)*!pi/180
oplot,cos(beta),sin(beta),col=0

; ------------------------------------------------------------------------------
; Trace final 2:
; ------------------------------------------------------------------------------

plot,pos(2,*),pos(1,*),psym=3,col=0,back=255,xran=[zmin,zmax],yran=[ymin,ymax],$
xtit='Z (Rsat)',ytit='Y',tit='Y(Z)';,/iso

; Trace additionnel de Saturne:
;------------------------------

beta=findgen(360)*!pi/180
oplot,cos(beta),sin(beta),col=0
a=(findgen(10)-5)*10
b=fltarr(10)
oplot,a,b,linestyle=2,col=0

; ------------------------------------------------------------------------------
; Trace final 3:
; ------------------------------------------------------------------------------

plot,pos(0,*),pos(2,*),psym=3,col=0,back=255,xran=[xmin,xmax],yran=[zmin,zmax],$
xtit='X (Rsat)',ytit='Z',tit='Z(X)';,/iso

; Trace additionnel de Saturne:
;------------------------------

beta=findgen(360)*!pi/180
oplot,cos(beta),sin(beta),col=0
a=(findgen(10)-5)*10
b=fltarr(10)
oplot,a,b,linestyle=2,col=0

; ------------------------------------------------------------------------------
; Trace final 4:
; ------------------------------------------------------------------------------

if keyword_set(lat) then begin
  plot,time-2556.,lat,psym=3,col=0,back=255,xtit='DOY 2004',ytit='Lat ',$
  tit='S/C Latitude (deg)'
endif else begin
  plot,time-2556.,rr,psym=3,col=0,back=255,xtit='DOY 2004 ',ytit='r',$
  tit='Distance S/C to Saturn (Rsat)'
endelse

if keyword_set(ps) then begin
  device,/close
  set_plot,'x'
endif

!p.multi=[0,1,1]
end