select_data_n3e_n3d.pro (Kronos IDL library)

;+
; Contains the select_data_n3e_n3d procedure
;
; :Author:
;    Laurent Lamy
;
; :History:
;    2006: Created
;
;    2009/05/18: File Created
;
;    2009/11/24: Last Edit
;-
;
;+
; Attribution d'un indice de -1 à 10 à chaque fichier n3e 
;
; :Uses:
;    q_prod, make_file_list, read_data_binary, read_src_list,
;    q_vmake, q_rot
;
; :Params:
;    n2:     in, required, type=sometype
;       donnees niveau 2  (chargees avec load_data_n3e)
;    n3e:    in, required, type=sometype
;       n3e: donnees niveau 3e (chargees avec load_data_n3e_n3d)
;    n3d:    in, required, type=sometype
;       donnees niveau 3d (chargees avec load_data_n3e_n3d)
;    source: in, required, type=sometype
;       source initial guess DF (code en 2 caracteres)
;    vmin:   in, required, type=sometype
;       Taux de polarisation circulaire min en valeur absolue
;    vmax:   in, required, type=sometype
;       Taux de polarisation circulaire max en valeur absolue
;    snmin:  in, required, type=sometype
;       S/N minimum (typiquement 3 ou 5 <=> auto-fond = 3 ou 5*sigma)
;    altmin: in, required, type=sometype
;       Altitude min par rapport au disque (valeur en Rp)
;    altmax: in, required, type=sometype
;       Altitude max par rapport au disque (valeur en Rp)
;    index:  out, required, type=int
;         -1  -> donnees non selectionnees
;          0  -> donnees nulles (S~0)
;          1  -> RH,Nord
;          2  -> LH,Nord
;          3  -> RH,Sud
;          4  -> LH,Sud
;          5  -> RH,centre
;          6  -> LH,centre
;          7  -> RH,grand nord
;          8  -> LH,grand nord
;          9  -> RH,grand sud
;          10 -> LH,grand sud
;    rr:     out, required, type=sometype
;       Coordonnees (en Rs) de Cassini dans le repere Saturn Solar
;       Equatorial (SSQ): distance (en Rs)
;    tl:     out, required, type=sometype
;       Coordonnees (en Rs) de Cassini dans le repere Saturn Solar
;       Equatorial (SSQ): temps local (0->24h)
;    zz:     out, required, type=sometype
;       Coordonnees (en Rs) de Cassini dans le repere Saturn Solar
;       Equatorial (SSQ): altitude au dessus du plan des anneaux (en Rs)
;    ydf:    in, required, type=sometype
;       resultats DF dans le repere final (obtenu par qrot0): -tan(azi)*rr
;       NB: dans ce repere, Saturne est a lat=0 et azi=0, 
;                       lat=0 <-> plan des anneaux, lat>0 <-> hemisphere N
;    zdf:    in, required, type=sometype
;       resultats DF dans le repere final (obtenu par qrot0): tan(lat)*rr
;       NB: dans ce repere, Saturne est a lat=0 et azi=0, 
;                       lat=0 <-> plan des anneaux, lat>0 <-> hemisphere N
;-
pro SELECT_DATA_N3E_N3D, n2, n3e, n3d, source, $
                 vmin,vmax,snmin,altmin,altmax, $
                 index,rr,tl,zz,$
                 ydf,zdf

; Chargement du type d'ephemerides:
; ---------------------------------

;    (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.

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)

; ua = 1 Unite astronomique en Rp

ua = 150.e6/src_list.body_rad_value

; Creation tableau de sortie:
; ---------------------------

nn = n_elements(n3e)

index = bytarr(nn)-1b
rr = fltarr(nn)
tl = fltarr(nn) 
zz = fltarr(nn)
ydf = fltarr(nn)
zdf = fltarr(nn)

jour_debut = t97_aj(long(min(n2.t97)))
jour_fin   = t97_aj(long(max(n2.t97)))

for jour_courant=jour_debut,jour_fin do begin
  
; Chargement des fichiers d'ephemerides du jour 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
	   
	   
       ; Preparation des Ephemerides:
       ;-----------------------------

       ; 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,*)

       ; Selection des donnees correspondant au jour courant:
       ;-----------------------------------------------------

       ww = where(long(t97_aj(n2.t97)) eq jour_courant, count)
       if count ne 0 then begin
         nww = n_elements(ww)
                  
         ; Creation du quaternion vecteur des resultats du DF:
         ;----------------------------------------------------

         vdf = make_vect_sph(1.,n3e(ww).th,n3e(ww).ph)
         qdf = q_vmake(vdf)

         ; Interpolation des Ephemerides:
         ;-------------------------------

         qrot = fltarr(4,nww)
         if nww gt 1 then begin
           qrot(0,*) = interpol(qrot0(0,*),qeph.time,n2(ww).t97)
           qrot(1,*) = interpol(qrot0(1,*),qeph.time,n2(ww).t97)
           qrot(2,*) = interpol(qrot0(2,*),qeph.time,n2(ww).t97)
           qrot(3,*) = interpol(qrot0(3,*),qeph.time,n2(ww).t97)
         endif
         
         qdf = q_rot(qrot,qdf)

         ; Direction finding results in new frame:
         ;----------------------------------------

         lat  = asin(qdf(3,*))
         azi  = atan(qdf(2,*),qdf(1,*))

         ; lat = latitude (angle XoY->DF)
         ; azi = azimuth  (angle X->DFsin(th2))
         
         ; Remplissage des coordonnees rr,lt,zz:
         ;--------------------------------------
         
         if nww gt 1 then begin
           rr(ww) = interpol(rr0,veph.time,n2(ww).t97)/src_list.body_rad_value
             qc0_1 = interpol(qc0(1,*),qeph.time,n2(ww).t97)
             qc0_2 = interpol(qc0(2,*),qeph.time,n2(ww).t97)
           tl(ww) = (atan(-qc0_2,-qc0_1)+!pi)/!pi*12.     
           zz(ww) = interpol(-qc0(3,*)/src_list.body_rad_value,qeph.time,n2(ww).t97)
         endif

         zdf(ww) =  tan(lat)*rr(ww)
         ydf(ww) = -tan(azi)*rr(ww)        
         
         ; Preparation des tableaux d'indices:
         ;------------------------------------
         sn = n3e(ww).sn
         
         ; Erreur max pour 99% des points sur lat/azi avec SN>23dB:
         ;---------------------------------------------------------
         teta = 5.*!dtor 
         dmin = rr(ww)*tan(atan(altmin/rr(ww))+teta)
         dmax = rr(ww)*tan(atan(altmax/rr(ww))+teta)
         
         ;"nul" = données non sélectionnées
         ww0 = where(n3e(ww).s eq 0,count0)

         ;RH nord
         ww1 = where( (sn(0,ww) gt snmin and sn(1,ww) gt snmin) and $
                      ( (sqrt(ydf^2+zdf^2) le dmax) and (zdf ge dmin) ) and $
                      (n3e(ww).v le -vmin) and (n3e(ww).v ge -vmax) and $
                      (n3d(ww).v le -vmin) and (n3d(ww).v ge -vmax), $
                    count1 )

         ;LH nord
         ww2 = where( (sn(0,ww) gt snmin and sn(1,ww) gt snmin) and $
                      ( (sqrt(ydf^2+zdf^2) le dmax) and (zdf ge dmin) ) and $
                      (n3e(ww).v ge vmin) and (n3e(ww).v le vmax) and $
                      (n3d(ww).v ge vmin) and (n3d(ww).v le vmax), $
                    count2 )   
                    
         ;RH sud               
         ww3 = where( (sn(0,ww) gt snmin and sn(1,ww) gt snmin) and $
                      ( (sqrt(ydf^2+zdf^2) le dmax) and (zdf le -dmin) ) and $
                      (n3e(ww).v le -vmin) and (n3e(ww).v ge -vmax) and $
                      (n3d(ww).v le -vmin) and (n3d(ww).v ge -vmax), $
                    count3 )
         ;LH sud
         ww4 = where( (sn(0,ww) gt snmin and sn(1,ww) gt snmin) and $
                      ( (sqrt(ydf^2+zdf^2) le dmax) and (zdf le -dmin) ) and $
                      (n3e(ww).v ge vmin) and (n3e(ww).v le vmax) and $
                      (n3d(ww).v ge vmin) and (n3d(ww).v le vmax), $
                    count4 )       
                    
         ;RH centre 
         ww5 = where( (sn(0,ww) gt snmin and sn(1,ww) gt snmin) and $
                      ( (abs(zdf) lt dmin) and (sqrt(ydf^2+zdf^2) le dmax) ) and $
                      (n3e(ww).v le -vmin) and (n3e(ww).v ge -vmax) and $
                      (n3d(ww).v le -vmin) and (n3d(ww).v ge -vmax), $
                    count5 )
                    
         ;LH centre 
         ww6 = where( (sn(0,ww) gt snmin and sn(1,ww) gt snmin) and $
                      ( (abs(zdf) lt dmin) and (sqrt(ydf^2+zdf^2) le dmax) ) and $
                      (n3e(ww).v ge vmin) and (n3e(ww).v le vmax) and $
                      (n3d(ww).v ge vmin) and (n3d(ww).v le vmax), $
                    count6 )

         ;RH "grand" nord
         ww7 = where( (sn(0,ww) gt snmin and sn(1,ww) gt snmin) and $
                      ( (sqrt(ydf^2+zdf^2) gt dmax) and (zdf ge 0.) ) and $
                      (n3e(ww).v le -vmin) and (n3e(ww).v ge -vmax) and $
                      (n3d(ww).v le -vmin) and (n3d(ww).v ge -vmax), $
                    count7 )                  
         
         ;LH "grand" nord
         ww8 = where( (sn(0,ww) gt snmin and sn(1,ww) gt snmin) and $
                      ( (sqrt(ydf^2+zdf^2) gt dmax) and (zdf ge 0.) ) and $
                      (n3e(ww).v ge vmin) and (n3e(ww).v le vmax) and $
                      (n3d(ww).v ge vmin) and (n3d(ww).v le vmax), $
                    count8 )  

         ;RH "grand" sud
         ww9 = where( (sn(0,ww) gt snmin and sn(1,ww) gt snmin) and $
                      ( (sqrt(ydf^2+zdf^2) gt dmax) and (zdf lt 0.) ) and $
                      (n3e(ww).v le -vmin) and (n3e(ww).v ge -vmax) and $
                      (n3d(ww).v le -vmin) and (n3d(ww).v ge -vmax), $
                    count9 ) 

         ;LH "grand" sud
         ww10 = where( (sn(0,ww) gt snmin and sn(1,ww) gt snmin) and $
                      ( (sqrt(ydf^2+zdf^2) gt dmax) and (zdf lt 0.) ) and $
                      (n3e(ww).v ge vmin) and (n3e(ww).v le vmax) and $
                      (n3d(ww).v ge vmin) and (n3d(ww).v le vmax), $
                    count10 ) 
         

         ; Remplissage du tableau d'indices:
         ;----------------------------------

         if count0  ne 0 then index(ww(ww0))  = 0b
         if count1  ne 0 then index(ww(ww1))  = 1b
         if count2  ne 0 then index(ww(ww2))  = 2b
         if count3  ne 0 then index(ww(ww3))  = 3b
         if count4  ne 0 then index(ww(ww4))  = 4b
         if count5  ne 0 then index(ww(ww5))  = 5b
         if count6  ne 0 then index(ww(ww6))  = 6b
         if count7  ne 0 then index(ww(ww7))  = 7b
         if count8  ne 0 then index(ww(ww8))  = 8b
         if count9  ne 0 then index(ww(ww9))  = 9b
         if count10 ne 0 then index(ww(ww10)) = 10b
       endif
endfor

end