solterv.f 6.27 KB
      subroutine solterv (iyear,imonth,iday,ihour,imin,isec,year,
     >                    alfag,alfas,deltas,obliq,tetdip,phidip,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.05 - V 2.0
c*VER 03.01.06 - V 2.1
c*VER 01.01.07 - V 3.0
c*
c*AUT spec. CNES - JC KOSIK - janvier 1998
c*AUT port. CISI
c*
c*ROL Theme : Astronomie et calculs d'orbite
c*ROL         Calcul de l'ascension droite de Greenwich.
c*ROL         Calcul de l'ascension droite et de la declinaison
c*ROL         du soleil.
c*ROL         Calcul de la position geographique du point ou le dipole
c*ROL         coupe l'hemisphere nord du 1er janvier 70 au 31 decembre 00. 
c*
c*PAR iyear  (I) : annee de la date (annee >= 1970 et < 2000)
c*PAR imonth (I) : mois de la date
c*PAR iday   (I) : jour de la date
c*PAR ihour  (I) : heures de la date
c*PAR imin   (I) : minutes de la date
c*PAR isec   (I) : secondes de la date
c*
c*PAR year   (O) : annee fractionnaire (>= 1970. et < 2000.)
c*
c*PAR alfag  (O) : ascension droite de Greenwich (radians)
c*
c*PAR alfas  (O) : ascension droite du soleil (radians)
c*PAR deltas (O) : declinaison du soleil (radians)
c*
c*PAR obliq  (O) : obliquite de l'ecliptique (radians)
c*
c*PAR tetdip (O) : colatitude geocentrique du dipole (radians)
c*PAR phidip (O) : longitude geocentrique du dipole (radians)
c*
c*PAR ier    (O) : code de retour
c*
c*NOT ier        : = 1 si iyear < 1970 ou iyear >= 2000
c*
c*INF utilise    : julg, sun
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 2.0 - 01.06.05 - correction de commentaires de code
c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
c*HST version 3.0 - 01.01.07 - extension jusqu'a l'annee 2000
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
      integer iyear, imonth, iday, ihour, imin, isec
      double precision year
      double precision alfag, alfas, deltas, obliq
      double precision tetdip, phidip
      integer ier
c
c     ----------------------------------
c*FON Declaration des fonctions externes
c     ----------------------------------
c
      external julg
      integer julg
c
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
      integer glg(3,7)
      integer julc,jul1,nujour,imin1,imax
      integer i,k
c*LOC i,k : indices de boucles et de tableaux
c
      double precision delt
c*LOC delt : intervalle entre l'annee courante et l'annee de reference
c
      double precision tzero(7)
c*LOC tzero : tableau des annees de reference de 1970 a 2000
c
      double precision fjour,g10,g11,c11,h11,gdt
      double precision gg(3,7)
      double precision tm, tmold
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*FON Affectation des constantes
c     --------------------------
c
      data (glg(i,1),i=1,3) /-30220.d0, -2068.d0, 5737.d0/
      data (glg(i,2),i=1,3) /-30100.d0, -2013.d0, 5675.d0/
      data (glg(i,3),i=1,3) /-29992.d0, -1956.d0, 5604.d0/
      data (glg(i,4),i=1,3) /-29873.d0, -1905.d0, 5500.d0/
      data (glg(i,5),i=1,3) /-29775.d0, -1848.d0, 5406.d0/
      data (glg(i,6),i=1,3) /-29682.d0, -1789.d0, 5318.d0/
      data (glg(i,7),i=1,3) /-29619.4d0, -1728.2d0, 5186.1d0/
c
      data tzero /1970.d0, 1975.d0, 1980.d0, 1985.d0,
     >            1990.d0, 1995.d0, 2000.d0/
      data tmold /0.d0/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier = 0
c
      if (iyear .lt. 2005 .and. iyear .ge. 1970) then
c
c     -----------------------------
c*FON Calcul de la date en secondes
c     -----------------------------
c
         fjour = dble ((ihour * 60 + imin) * 60 + isec) / 86400.d0
c
c     -------------------------------
c*FON Calcul du jour julien CNES 1950
c     -------------------------------
c
         julc  = julg (iday,imonth,iyear,ier)
c
         jul1 = julg(1,1,iyear,ier)
c
c   ------------------------------------------------------------
c*FON Calcul du numero du jour compte a partir du 1er janvier de  
c*FON l'annee courante
c     ----------------------------------------------------------
c
         nujour = julc - jul1 + 1
c
c     ---------------------------------------------------------
c*FON Calcul de la position du soleil dans le repere equatorial
c*FON a gamma 50
c     ---------------------------------------------------------
c
         call sun (iyear,nujour,fjour,alfag,alfas,deltas,obliq,ier)
         year = dble(iyear) + dble(nujour-1)/365.d0 
c
c     -------------------------------------------------------------
c*FON Calculs intermediaires pour obtenir les coordonnees du dipole
c     --------------------------------------------------------------
c 
         tm = dble(iyear) + dble(imonth) / 12.d0
         if (tm .ne. tmold) then
c
            tmold = tm
            imin1 = idint((tm - 1970.d0) / 5.d0) + 1
            imax  = imin1 + 1
            delt  = (tm - tzero(imin1)) / 5.d0
c
            do 30 k = 1, 3
                gg(k,imin1) = glg(k,imin1)
     >                        + (glg(k,imax)
     >                        - glg(k,imin1)) * delt
30          continue
c
            g10 = - gg(1,imin1)
            g11 = - gg(2,imin1)
            h11 = - gg(3,imin1)
c
            c11 = sqrt(g11**2 + h11**2)
            gdt = sqrt(g11**2 + h11**2 + g10**2)
c
c     ---------------------------------------------------------------
c*FON Calcul de la longitude et de la colatitude du dipole en radians
c     ---------------------------------------------------------------
c
            phidip = atan(h11 / g11)
            tetdip = asin(c11 / gdt)
c
         endif
      else
         ier = 1
      endif
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end