select_data_n3b.pro (Kronos IDL library)

;+
; Contains the select_data_n3b procedure
;
; :Author:
;    Laurent Lamy
;
; :History:
;    2011/01/20: Created
;
;    2011/01/20: Last Edit
;-
;
;+
; select_data_n3b is a procedure that <behavior desc here>
;
; :Uses:
;    q_prod
;
; :Params:
;    n2:     in, required, type=sometype
;       A parameter named n2
;    n3b:    in, required, type=sometype
;       A parameter named n3b
;    source: in, required, type=sometype
;       A parameter named source
;    snmin:  in, required, type=sometype
;       A parameter named snmin
;    altmin: in, required, type=sometype
;       A parameter named altmin
;    altmax: in, required, type=sometype
;       A parameter named altmax
;    index:  in, required, type=sometype
;       A parameter named index
;    rr:     in, required, type=sometype
;       A parameter named rr
;    tl:     in, required, type=sometype
;       A parameter named tl
;    zz:     in, required, type=sometype
;       A parameter named zz
;    ydf:    in, required, type=sometype
;       A parameter named ydf
;    zdf:    in, required, type=sometype
;       A parameter named zdf
;    azi:    in, required, type=sometype
;       A parameter named azi
;-
pro SELECT_DATA_N3B, n2, n3b, source, $
                 snmin,altmin,altmax, $
                 index,rr,tl,zz,$
                 ydf,zdf,azi
;======================================================================================

; 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(n3b)

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

jour_debut = t97_aj(long(n2(0).t97))
jour_fin   = t97_aj(long(n2(nn-1L).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.,n3b(ww).th,n3b(ww).ph)
         qdf = q_vmake(vdf)


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

         qrot = fltarr(4,nww)
         if nww gt 2 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 2 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 = n3b(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( (sn(0,*) le snmin) or (sn(1,*) le snmin), $
                    count0)

         ; Nord
         ww1 = where( (sn(0,*) gt snmin and sn(1,*) gt snmin) and $
                      ( (sqrt(ydf^2+zdf^2) le dmax) and (zdf ge dmin) ), $
                    count1 )
                    
         ; Sud               
         ww2 = where( (sn(0,*) gt snmin and sn(1,*) gt snmin) and $
                      ( (sqrt(ydf^2+zdf^2) le dmax) and (zdf le -dmin) ), $
                    count2 )
                       
         ; Centre 
         ww3 = where( (sn(0,*) gt snmin and sn(1,*) gt snmin) and $
                      ( (abs(zdf) lt dmin) and (sqrt(ydf^2+zdf^2) le dmax) ), $
                    count3 )

         ; "grand" nord
         ww4 = where( (sn(0,*) gt snmin and sn(1,*) gt snmin) and $
                      ( (sqrt(ydf^2+zdf^2) gt dmax) and (zdf ge 0.) ), $
                    count4 )                  
   
         ; "grand" sud
         ww5 = where( (sn(0,*) gt snmin and sn(1,*) gt snmin) and $
                      ( (sqrt(ydf^2+zdf^2) gt dmax) and (zdf lt 0.) ), $
                    count5 ) 

         n_ind=(count0+count1+count2+count3+count4+count5)
         if n_ind ne uniq_list(n_ind) then stop,'Warning: error in index procedure'
         ; et pas n_elements(ww) car parfois zdf = nan
         
         ; Remplissage du tableau d'indices:
         ;----------------------------------

         if count0  ne 0 then index(ww(ww0))  = 0b
         if count1  ne 0 then index(ww(ww1))  = 11b
         if count2  ne 0 then index(ww(ww2))  = 12b
         if count3  ne 0 then index(ww(ww3))  = 13b
         if count4  ne 0 then index(ww(ww4))  = 14b
         if count5  ne 0 then index(ww(ww5))  = 15b

       endif
endfor

end