cahsl.f 4.26 KB
      subroutine cahsl (alfag,alfas,deltas,rre,thets,phis,hsl,ihsl,ier)
c*
c***********************************************************************
c*
c*          "Copyright [c] CNES 98 - tous droits reserves"
c*          **********************************************
c*
c*PRO MAGLIB
c*
c*VER 99.03.31 - V 1.0
c*VER 01.06.01 - V 2.0
c*VER 03.01.06 - V 2.1
c*
c*AUT spec. CNES - JC KOSIK - janvier 1991
c*AUT port. CISI
c*
c*ROL Theme : Frontieres et regions
c*ROL         Calcul de la position relative du satellite par
c*ROL         rapport au cone d'ombre de la terre.
c*ROL         La position du satellite est donnee par ses
c*ROL         coordonnees geocentriques.
c*
c*PAR alfag  (I) : ascension droite de Greenwich (radians)
c*PAR alfas  (I) : ascension droite du soleil (radians)
c*PAR deltas (I) : declinaison du soleil (radians)
c*
c*PAR rre    (I) : distance radiale geocentrique du satellite
C*               : (rayons terrestres)
c*PAR thets  (I) : colatitude geocentrique du satellite (radians)
c*PAR phis   (I) : longitude geocentrique du satellite (radians)
c*
c*PAR hsl    (O) : distance satellite - cylindre d'ombre (km)
c*
c*PAR ihsl   (O) : position du satellite par rapport a l'ombre
c*
c*PAR ier    (O) : code de retour
c*
c*NOT ihsl       : ihsl = 1 : le satellite est dans l'ombre
c*NOT ihsl       : ihsl = 0 : le satellite est hors de l'ombre
c*
c*NOT ier        : sans objet
c*
c*NOT common     : util
c*
c*INF utilise    : sans objet
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 2.0 - 01.06.01 - correction de commentaires de code
c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
c*
c***********************************************************************
c*
      implicit none
c
c     ---------------------------------
c*FON Declaration identificateur rcs_id
c     ---------------------------------
c
      character rcs_id*100
c
c     --------------------------
c*FON Declaration des parametres
c     --------------------------
c
      double precision alfag, alfas, deltas
      double precision rre, thets, phis
      double precision hsl
      integer ihsl
      integer ier
c
c     ----------------------------------
c*FON Declaration des variables communes
c     ----------------------------------
c
      double precision pi,dpi,rad,deg,pid,xmu,rayt
c
c*COM pi   : constante pi (obtenue a partir de acos(-1.))
c*COM dpi  : constante 2 * pi
c*COM pid  : constante pi / 2
c*COM rad  : facteur de conversion degres  ----> radians
c*COM deg  : facteur de conversion radians ----> degres
c*COM xmu  : constante de gravitation terrestre (km**3/sec**2)
c*COM rayt : rayon equatorial terrestre (km)
c
      common/util/pi,dpi,rad,deg,pid,xmu,rayt
c
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
      double precision rs,delsa,alsa,delh
      double precision alh,cu,ur,hsat,hshade,diff
c*LOC Variables de travail intermediaires
c
      SAVE
c
c     ---------------------------------
c*FON Affectation identificateur rcs_id
c     ---------------------------------
c
      data rcs_id /"
     >$Id$"/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier  = 0
c
c     -----------------------------------------------
c*FON Calcul de la position geocentrique du satellite
c     -----------------------------------------------
c
      rs = rre * rayt
c
c     ----------------------------------------------------
c*FON Calcul de l'ascension et de la declinaison associees
c     ----------------------------------------------------
c
      delsa =  pid - thets
      alsa  =  alfag + phis
      delh  = -deltas
      alh   =  alfas + pi
      alh   =  mod(alh,dpi)
      cu    =  sin(delsa) * sin(delh) +
     >         cos(delsa) * cos(delh) * cos(alh - alsa)
      ur    =  acos(cu)
c
c     -------------------------------------
c*FON Distance du satellite au cone d'ombre
c     -------------------------------------
c
      ihsl = 0
      hsl  = 999.0d0
c
      if (ur .lt. pid) then
         hsat   = rs * sin(ur)
         hshade = 6370.d0
         diff   = hshade - hsat
c
         if (diff .ge. 0.d0) then
            hsl  = diff
            ihsl = 1
         endif
      endif
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end