grille_sat.pro (Kronos IDL library)

;+
; Contains the grille_sat procedure
;
; :Author:
;    Laurent Lamy
;
; :History:
;    1992: Created by Guy Munhoven (Liege)
;
;    1993/08/**: Modified by Vincent Dols
;
;    1994/07/**: Modified by Renee Prange juillet  avril 1995, correction des latitudes (août 1993)
;
;    2005/02/**: Adapte de elisa de thierry Fouchet par Laurent Pallier
;
;    2006/01/**: Corrige par Renee Prange
;
;    2006/02/**: Corrige et modifie par Laurent Lamy de fevrier 2006 à Janvier 2008
;
;    2009/06/01: Last Edit
;-
;
;+
; Procédure de tracé d'une grille de longitudes/parallèles/anneaux
; pour Saturne quelle que soit la position de l'observateur
;
; :Params:
;    req:      in, required, type=sometype
;       Rayon equatorial
;    e:        in, required, type=sometype
;       Applatissement de la planete
;    CML:      in, required, type=sometype
;       Méridien Subterrestre
;    pointa:   out, required, type=sometype
;       Pointage HST exprime en latitude, longitude
;    grille1:  out, required, type=sometype
;       Matrice des parallèles
;
;       grille1(a,b,c) est le tableau des coordonnées y,z des différents points définissant
;       la grille de parallèles projetée sur le plan de l'image. a=nbre de parallèles, b=nbre 
;       de points par parallèles, c=2 pour les vecteurs y et z
;    grille2:  out, required, type=sometype
;       Matrice des méridiens
;
;       grille2(a,b,c) est le tableau des coordonnées y,z des différents points définissant 
;       la grille de méridiens projetée sur le plan de l'image. a=nbre de méridiens, b=nbre 
;       de points par méridien, c=2 pour les vecteurs y et z (correspond a xx2,xx3)
;    feature:  out, required, type=sometype
;       Matrice de l'équateur, du méridien central, du contour et du méridien 180°
;
;       feature(a,b,c) est le tableau des coordonnées y,z des différents points définissant 
;       certaines courbes particulières projetée sur le plan de l'image. a=nbre de courbes 
;       spéciales, de 0 a 3 correspond à contour, equateur, meridien central et meridien a 180°,
;       b=nbre de points par courbes, c=2 pour les vecteurs y et z
;    lieu:     out, required, type=sometype
;       Matrice des aurores
;
;       lieu(a,b,c) est le tableau des coordonnées y,z des différents points définissant 
;       les aurores éventuelles
;    anneau:   out, required, type=sometype
;       Matrice des anneaux
;
;       anneau(a,b,c) est le tableau des coordonnées y,z des différents points définissant 
;       les anneaux projetée sur le plan de l'image. a=nbre de frontières d'anneaux, b=nbre 
;       de points par frontières d'anneaux, c=2 pour les vecteurs y et z
;    fov:      out, required, type=sometype
;       Pointage HST (field of view)
;
; :Keywords:
;    alpha:    in, optional, type=sometype
;       Angle de Pitch ou Position Angle
;    beta:     in, optional, type=sometype
;       Angle de tilt
;    plot:     in, optional, type=sometype
;       A keyword named plot
;    meri:     in, optional, type=sometype
;       A keyword named meri
;    par:      in, optional, type=sometype
;       A keyword named par
;    equ:      in, optional, type=sometype
;       A keyword named equ
;    mer_cml:  in, optional, type=sometype
;       A keyword named mer_cml
;    anti_cml: in, optional, type=sometype
;       A keyword named anti_cml
;-
pro grille_sat,alpha=alp,beta=bet,req,e,CML,pointa,grille1,grille2,feature,lieu,anneau,$
           fov,plot=plot,meri=meri,par=par,equ=equ,mer_cml=mer_cml,anti_cml=anti_cml

;-------------------------------------------------------------------------------
; Matrice des lieux:
;-------------------------------------------------------------------------------

pointa=fltarr(2)
grille1=fltarr(19,720,2)
grille2=fltarr(24,337,2)
feature=fltarr(4,1440,2)
anneau=fltarr(5,1440,2)
lieu=fltarr(5,180,2)
fov=fltarr(2)

;-------------------------------------------------------------------------------
; Matrice de changement de repere:
;-------------------------------------------------------------------------------

sa=double(sin(alp*!dtor))
ca=double(cos(alp*!dtor))
sb=double(sin(bet*!dtor))
cb=double(cos(bet*!dtor))
tb=double(tan(bet*!dtor))

;matrice de changement de repère de GM (que RP a changee en Aout 93)
;a11=cb      &      a12=sa*sb      &      a13=ca*sb
;a21=0.      &      a22=ca         &      a23=-sa
;a31=-sb     &      a32=sa*cb      &      a33=ca*cb

;nouvelle matrice de RP
a11=cb        &      a12=0.       &      a13=sb
a21=sa*sb     &      a22=ca       &      a23=-sa*cb
a31=-sb*ca    &      a32=sa       &      a33=ca*cb

;-------------------------------------------------------------------------------
; PARAMETRES de la planete:
;-------------------------------------------------------------------------------

e1=1./(1-e)
req2=req^2
rp=(1-e)*req
; Ancien rayon polaire apparent: rpa=req*(1-e*(a11)^2) (visiblement faux, cf calcul limbe)

;-------------------------------------------------------------------------------
; CONTOUR du limbe:
;-------------------------------------------------------------------------------

; Calcul préalable du rayon polaire apparent:
; Equation de l'ellipse après rotation
;    x1=req*cb*cos(zeta)+rp*sb*sin(zeta)
;    x3=-req*sb*cos(zeta)+rp*cb*sin(zeta)
;    rpa=x3_max
zeta=atan(-rp/(tb*req)) 
rpa=abs(-req*sb*cos(zeta)+rp*cb*sin(zeta)) 

l=findgen(1440)
lambda=l*!dtor/4.
d2=req*sin(lambda)  ; Coordonnees dans le plan (req,rpa)
d3=rpa*cos(lambda)  ; Idem
x2=a22*d2+a32*d3    ; Coordonnees dans le plan (req,nord) Terre
x3=a32*d2-a22*d3    ; Idem

feature(0,l,0)=x2  
feature(0,l,1)=x3  
rlimbe=sqrt(x2^2+x3^2) & tanlimbe=atan(x3,x2)
x2test=x2 & x3test=x3

;-------------------------------------------------------------------------------
; ANNEAUX: 
;-------------------------------------------------------------------------------

; Trace les limites interieure et exterieure des anneaux A,B et C:
for an=0,4 do begin
  case an of
    0: ran=1.23
    1: ran=1.52
    2: ran=1.95
    3: ran=2.02
    4: ran=2.27
  endcase 

; Sur 1 anneau point tous les 1/4 deg:
  p=findgen(1440)
  psi=p*!dtor/4.
  xx1=req*ran*cos(psi)
  xx2=req*ran*sin(psi)
  xx3=0.
  x1=a11*xx1+a12*xx2+a13*xx3
  x2=a21*xx1+a22*xx2+a23*xx3
  x3=a31*xx1+a32*xx2+a33*xx3
  w=where(abs(x2) le 1e-30) & if w(0) ne -1 then x2(w)=1e-30
  rl=sqrt(x2^2+x3^2) & tananneau=atan(x3,x2)

  
; Sélection du rlimbe correspondant à la condition de visibilité ci-après:
  if min(rl) le max(rlimbe) then begin
    previs=fltarr(1440)
    for i=0,1439l do begin
      test=abs(tanlimbe-tananneau(i))
      retest=where(test eq min(test),count_test)
      if count_test ne 0 then previs(i)=retest(0)
    endfor
  endif else previs=findgen(1440)
  
; Condition de visibilite des anneaux: 
  vis=where(rl ge rlimbe(previs) or x1 ge 0,count_vis)
  if count_vis ne 0 then begin
    anneau(an,vis,0)=x2(vis)
    anneau(an,vis,1)=x3(vis)
  endif
endfor

; Condition de visibilité du limbe au travers des anneaux
tanan_int=atan(anneau(0,*,1),anneau(0,*,0)) & rl_int=sqrt(anneau(0,*,0)^2+anneau(0,*,1)^2) 
tanan_ext=atan(anneau(4,*,1),anneau(4,*,0)) & rl_ext=sqrt(anneau(4,*,0)^2+anneau(4,*,1)^2)
for i=0l,1439l do begin
  test_int=abs(tanan_int-tanlimbe(i)) & test_ext=abs(tanan_ext-tanlimbe(i))
  if abs(min(test_int))*!radeg le 1. and abs(min(test_ext))*!radeg le 1. then begin
    retest_int=where(test_int eq min(test_int)) & retest_ext=where(test_ext eq min(test_ext))
    if retest_int(0) ne -1 and retest_ext(0) ne -1 then begin
      if rlimbe(i) ge rl_int(retest_int(0)) and rlimbe(i) le rl_ext(retest_ext(0)) then begin
        feature(0,i,0) = 0.
        feature(0,i,1) = 0.
      endif
    endif
  endif  
endfor

if keyword_set(plot) then begin
  w=where(feature(0,*,0) ne 0 and feature(0,*,1) ne 0)
  if w(0) ne -1 then oplot,feature(0,w,0),feature(0,w,1),psym=3
  for i=0,4 do begin
    wi=where(anneau(i,*,0) ne 0 or anneau(i,*,1) ne 0)
    if wi(0) ne -1 then oplot,anneau(i,wi,0),anneau(i,wi,1),psym=3
  endfor
endif

;-------------------------------------------------------------------------------
; PARALLELES:
;-------------------------------------------------------------------------------

; Parallele tous les 10 deg:
for l=-9,9 do begin

; Calcul du rayon local rl a la latitude planetocentrique lambda:
  lambda=10.*l*!dtor
  rl=sqrt(req2/(cos(lambda)^2+(sin(lambda)*e1)^2)) ; cf eq de l'ellipse = 1/sqtr((cos(lambda)/req)^2+(sin(lambda)/rp)^2)
  d12=rl*cos(lambda)

; Sur 1 parallele 1 point tous les 1/2 deg:
  p=findgen(720) 
  psi=p*!dtor/2.
  xx1=d12*cos(psi)
  xx2=d12*sin(psi)
  xx3=rl*sin(lambda)
  x1=a11*xx1+a12*xx2+a13*xx3
  x2=a21*xx1+a22*xx2+a23*xx3
  x3=a31*xx1+a32*xx2+a33*xx3

; Recherche condition de visibilité vraie: (Calcul à 2D)
;   Equation de l'ellipse après rotation
;     x1=req*cb*cos(zeta)+rp*sb*sin(zeta)
;     x3=-req*sb*cos(zeta)+rp*cb*sin(zeta)
;   Recherche de limite:
;     zeta=atan(-rp/(tb*req)) 
;     x3_max=abs(-req*sb*cos(zeta)+rp*cb*sin(zeta))  
;     x1_lim=abs(req*cb*cos(zeta)+rp*sb*sin(zeta))
  
; Recherche condition de visibilité vraie: (2 Rotations successives à 2D)
; 1ère rotation dans le plan (x1,x3):
  zeta=atan(-rp/(tb*req))
  x1_lim=abs(req*cb*cos(zeta)*cos(psi)+rp*sb*sin(zeta)*cos(psi))
;   N.B: Les ellipses dans (x1,x3) ont pour côtés req*cos(psi) et rp*cos(psi)
;        On a ainsi x1_lim pour chacun des points des parallèles

; Condition de visibilité finale:
;   On compare x1 à x1_lim (pas de changement par la 2ème rotation)
;   On compare xx3*cb à 0 car les points concernés ne changent pas avec la 2ème rotation 
  if ((bet ge 0. and bet lt 90.) or (bet ge 180. and bet lt 270.)) then begin
    vis=where((xx3*cb ge 0. and x1 ge -x1_lim) or (xx3*cb lt 0. and x1 ge x1_lim),count_vis)
  endif else begin
    vis=where((xx3*cb ge 0 and x1 ge x1_lim) or (xx3*cb lt 0 and x1 ge -x1_lim),count_vis)
  endelse  

  if count_vis ne 0 then begin
    grille1(l+9,vis,0)=x2(vis)
    grille1(l+9,vis,1)=x3(vis)
  endif
  

  ; Condition de visibilité au travers des anneaux
  tan_par=atan(grille1(l+9,*,1),grille1(l+9,*,0)) & rl_par=sqrt(grille1(l+9,*,0)^2+grille1(l+9,*,1)^2) 
  for i=0l,n_elements(grille1(l+9,*,0))-1l do begin
    test_int=abs(tanan_int-tan_par(i)) & test_ext=abs(tanan_ext-tan_par(i))
    if abs(min(test_int))*!radeg le 1. and abs(min(test_ext))*!radeg le 1. then begin
      retest_int=where(test_int eq min(test_int)) & retest_ext=where(test_ext eq min(test_ext))
      if retest_int(0) ne -1 and retest_ext(0) ne -1 then begin
        if rl_par(i) ge rl_int(retest_int(0)) and rl_par(i) le rl_ext(retest_ext(0)) then begin
          grille1(l+9,i,0) = 0.
          grille1(l+9,i,1) = 0.
        endif
      endif
    endif  
  endfor
    
  if keyword_set(plot) and keyword_set(par) then begin
    w=where(grille1(l+9,*,0) ne 0 or grille1(l+9,*,1) ne 0)
    if w(0) ne -1 then oplot,grille1(l+9,w,0),grille1(l+9,w,1),psym=3
  endif

endfor

;-------------------------------------------------------------------------------
; EQUATEUR:
;-------------------------------------------------------------------------------

p=findgen(1440) 
psi=p*!dtor/4.    
xx1=req*cos(psi)
xx2=req*sin(psi)
xx3=0.
x1=a11*xx1+a12*xx2+a13*xx3
x2=a21*xx1+a22*xx2+a23*xx3
x3=a31*xx1+a32*xx2+a33*xx3

; Recherche condition de visibilité vraie:
vis=where(x1 ge 0,count_vis) 
if count_vis ne 0 then begin
  feature(1,vis,0)=x2(vis)
  feature(1,vis,1)=x3(vis)
endif

; Condition de visibilité au travers des anneaux
tan_eq=atan(feature(1,*,1),feature(1,*,0)) & rl_eq=sqrt(feature(1,*,0)^2+feature(1,*,1)^2) 
for i=0l,n_elements(feature(1,*,0))-1l do begin
  test_int=abs(tanan_int-tan_eq(i)) & test_ext=abs(tanan_ext-tan_eq(i))
  if abs(min(test_int))*!radeg le 1. and abs(min(test_ext))*!radeg le 1. then begin
    retest_int=where(test_int eq min(test_int)) & retest_ext=where(test_ext eq min(test_ext))
    if retest_int(0) ne -1 and retest_ext(0) ne -1 then begin
      if rl_eq(i) ge rl_int(retest_int(0)) and rl_eq(i) le rl_ext(retest_ext(0)) then begin
        feature(1,i,0) = 0.
        feature(1,i,1) = 0.
      endif
    endif
  endif  
endfor

if keyword_set(plot) and keyword_set(equ) then begin
  w=where(feature(1,*,0) ne 0 or feature(1,*,1) ne 0)
  if w(0) ne -1 then oplot,feature(1,w,0),feature(1,w,1),psym=3
endif

;-------------------------------------------------------------------------------
; MERIDIENS:
;-------------------------------------------------------------------------------

; Meridien tous les 15 deg:
for p=0l,23l do begin 
  psi=15.*p*!dtor

; Point tous les 1/2 deg de -84 a +84°:
  l=findgen(337)
  lambda=(l-168.)*!dtor/2.
  rl=sqrt(req2/(cos(lambda)^2+(sin(lambda)*e1)^2))
  xx1=rl*cos(psi+CML)*cos(lambda)
  xx2=rl*sin(psi+CML)*cos(lambda)
  xx3=rl*sin(lambda)
  x1=a11*xx1+a12*xx2+a13*xx3
  x2=a21*xx1+a22*xx2+a23*xx3
  x3=a31*xx1+a32*xx2+a33*xx3
                          
; Recherche condition de visibilité vraie: (cf parallèles pour détail)
  zeta=atan(-rp/(tb*req))
  x1_lim=abs(req*cb*cos(zeta)*cos(psi)+rp*sb*sin(zeta)*cos(psi))
; N.B: On a ainsi x1_lim pour chacun des méridiens

; Condition de visibilité finale:
;   On compare x1 à x1_lim (pas de changement par la 2ème rotation)
;   On compare xx3*cb à 0 car les points concernés ne changent pas avec la 2ème rotation 
  if ((bet ge 0. and bet lt 90.) or (bet ge 180. and bet lt 270.)) then begin
    vis=where((xx3*cb ge 0 and x1 ge -x1_lim) or (xx3*cb lt 0 and x1 ge x1_lim),count_vis)
  endif else begin
    vis=where((xx3*cb ge 0 and x1 ge x1_lim) or (xx3*cb lt 0 and x1 ge -x1_lim),count_vis)
  endelse  
    
  if count_vis ne 0 then begin
    grille2(p,vis,0)=x2(vis)
    grille2(p,vis,1)=x3(vis)
  endif
  
  ; Condition de visibilité au travers des anneaux
  tan_mer=atan(grille2(p,*,1),grille2(p,*,0)) & rl_mer=sqrt(grille2(p,*,0)^2+grille2(p,*,1)^2) 
  for i=0l,n_elements(grille2(p,*,0))-1l do begin
    test_int=abs(tanan_int-tan_mer(i)) & test_ext=abs(tanan_ext-tan_mer(i))
    if abs(min(test_int))*!radeg le 1. and abs(min(test_ext))*!radeg le 1. then begin
      retest_int=where(test_int eq min(test_int)) & retest_ext=where(test_ext eq min(test_ext))
      if retest_int(0) ne -1 and retest_ext(0) ne -1 then begin
        if rl_mer(i) ge rl_int(retest_int(0)) and rl_mer(i) le rl_ext(retest_ext(0)) then begin
          grille2(p,i,0) = 0.
          grille2(p,i,1) = 0.
        endif
      endif
    endif  
  endfor
    
  if keyword_set(plot) and keyword_set(meri) then begin
    w=where(grille2(p,*,0) ne 0 or grille2(p,*,1) ne 0)
    if w(0) ne -1 then oplot,grille2(p,w,0),grille2(p,w,1),psym=3
  endif

endfor

;-------------------------------------------------------------------------------
; MERIDIEN CENTRAL:
;-------------------------------------------------------------------------------

; En gras le meridien 0 deg:
l=findgen(720) 
lambda=(l-360)*!dtor/4.
rl=sqrt(req2/(cos(lambda)^2+(sin(lambda)*e1)^2))
xx1=rl*cos(CML)*cos(lambda)
xx2=rl*sin(CML)*cos(lambda)
xx3=rl*sin(lambda)
x1=a11*xx1+a12*xx2+a13*xx3
x2=a21*xx1+a22*xx2+a23*xx3
x3=a31*xx1+a32*xx2+a33*xx3

; Recherche condition de visibilité vraie:
zeta=atan(rp/(tan(bet*!dtor)*req)) 
;x3_max=abs(req*sin(bet*!dtor)*cos(zeta)+rp*cos(bet*!dtor)*sin(zeta))
x1_lim=abs(req*cos(bet*!dtor)*cos(zeta)-rp*sin(bet*!dtor)*sin(zeta))
if ((bet ge 0. and bet lt 45.) or (bet ge 180. and bet lt 270.)) then begin
  vis=where((x3 ge 0 and x1 ge -x1_lim) or (x3 lt 0 and x1 ge x1_lim),count_vis)
endif else begin
  vis=where((x3 ge 0 and x1 ge x1_lim) or (x3 lt 0 and x1 ge -x1_lim),count_vis)
endelse 

if count_vis ne 0 then begin
  feature(2,l(vis),0)=x2(vis)
  feature(2,l(vis),1)=x3(vis)
endif

; Condition de visibilité au travers des anneaux
tan_cml=atan(feature(2,*,1),feature(2,*,0)) & rl_cml=sqrt(feature(2,*,0)^2+feature(2,*,1)^2) 
for i=0l,n_elements(feature(1,*,0))-1l do begin
  test_int=abs(tanan_int-tan_cml(i)) & test_ext=abs(tanan_ext-tan_cml(i))
  if abs(min(test_int))*!radeg le 1. and abs(min(test_ext))*!radeg le 1. then begin
    retest_int=where(test_int eq min(test_int)) & retest_ext=where(test_ext eq min(test_ext))
    if retest_int(0) ne -1 and retest_ext(0) ne -1 then begin
      if rl_cml(i) ge rl_int(retest_int(0)) and rl_cml(i) le rl_ext(retest_ext(0)) then begin
        feature(2,i,0) = 0.
        feature(2,i,1) = 0.
      endif
    endif
  endif  
endfor

if keyword_set(plot) and keyword_set(mer_cml) then begin
  w=where(feature(2,*,0) ne 0 or feature(2,*,1) ne 0)
  if w(0) ne -1 then oplot,feature(2,w,0),feature(2,w,1),psym=3,col=150
endif

;-------------------------------------------------------------------------------
; MERIDIEN A 180 DEGRES (en pointilles): 
;-------------------------------------------------------------------------------

; Pointille = 1 point sur 4 par rapport au meridien central:
l=findgen(720) 
lambda=(l-360)*!dtor/4+180.*!dtor
rl=sqrt(req2/(cos(lambda)^2+(sin(lambda)*e1)^2))
xx1=rl*cos(CML+!pi)*cos(lambda)
xx2=rl*sin(CML+!pi)*cos(lambda)
xx3=rl*sin(lambda)
x1=a11*xx1+a12*xx2+a13*xx3
x2=a21*xx1+a22*xx2+a23*xx3
x3=a31*xx1+a32*xx2+a33*xx3

vis=where(xx1 lt 0. or x1 ge 0,count_vis) 
if count_vis ne 0 then begin
  feature(3,l(vis),0)=x2(vis)
  feature(3,l(vis),1)=x3(vis)
endif

; Condition de visibilité au travers des anneaux
tan_180=atan(feature(3,*,1),feature(3,*,0)) & rl_180=sqrt(feature(3,*,0)^2+feature(3,*,1)^2) 
for i=0l,n_elements(feature(1,*,0))-1l do begin
  test_int=abs(tanan_int-tan_180(i)) & test_ext=abs(tanan_ext-tan_180(i))
  if abs(min(test_int))*!radeg le 1. and abs(min(test_ext))*!radeg le 1. then begin
    retest_int=where(test_int eq min(test_int)) & retest_ext=where(test_ext eq min(test_ext))
    if retest_int(0) ne -1 and retest_ext(0) ne -1 then begin
      if rl_180(i) ge rl_int(retest_int(0)) and rl_180(i) le rl_ext(retest_ext(0)) then begin
        feature(3,i,0) = 0.
        feature(3,i,1) = 0.
      endif
    endif
  endif  
endfor

if keyword_set(plot) and keyword_set(anti_cml) then begin
  w=where(feature(3,*,0) ne 0 or feature(3,*,1) ne 0)
  if w(0) ne -1 then oplot,feature(3,w,0),feature(3,w,1),psym=3,col=150
endif

return
end