soltervo.f 6.16 KB
      subroutine soltervo (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*
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 45 au 31 decembre 69.
c*
c*PAR iyear  (I) : annee de la date (annee >= 1945 et < 1970)
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 (>= 1945. et < 1970.)
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 < 1945 ou iyear >= 1970
c*
c*NOT L'annee 1990 represente l'origine des temps pour les
c*NOT calculs du champ magnetique terrestre, il faut donc
c*NOT que le parametre iyear soit superieur ou egal a 1990.
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*
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
      double precision obliq, 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 lg(3,6)
      integer julc, jul1, nujour, imin1, imax
      integer i,k
c*LOC i,k : indices de boucles et de tableaux
c
      double precision gg(3,6)
      double precision delt
c*LOC delt : intervalle entre l'annee courante et l'annee de reference
c
      double precision tzero(6)
c*LOC tzero : tableau des annees de reference de 1945 a 1970
c
      double precision fjour,g10,g11,c11,h11,gdt
      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 (lg(i,1),i=1,3) /-30594, -2285, 5810/
      data (lg(i,2),i=1,3) /-30554, -2250, 5815/
      data (lg(i,3),i=1,3) /-30500, -2215, 5820/
      data (lg(i,4),i=1,3) /-30421, -2169, 5791/
      data (lg(i,5),i=1,3) /-30334, -2119, 5776/
      data (lg(i,6),i=1,3) /-30220, -2068, 5737/
c
      data tzero /1945.d0, 1950.d0, 1955.d0, 1960.d0,
     >            1965.d0, 1970.d0/
      data tmold /0.d0/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier = 0
c
      if (iyear .lt. 1970 .and. iyear .ge. 1945) 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
c*FON janvier de l'annee courante
c     -----------------------------------------------
c
         nujour = julc - jul1 + 1

c     ----------------------------------------------
c*FON Calcul de la position du soleil dans le repere
c*FON equatorial 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
c*FON du dipole
c     ---------------------------------------------------
c
         tm = dble(iyear) + dble(imonth) / 12.d0
         if (tm .ne. tmold) then
c
            tmold = tm
            imin1 = idint((tm - 1945.d0) / 5.d0) + 1
            imax  = imin1 + 1
            delt = (tm - tzero(imin1)) / 5.d0
c
            do 30 k = 1, 3
               gg(k,imin1) = dble(lg(k,imin1))
     >                       + dble(lg(k,imax)
     >                        - lg(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
c*FON 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