source: Sophya/trunk/SophyaLib/NTools/datime.c@ 2808

#include "machdefs.h"
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#ifdef __MWERKS__
#include "unixmac.h"
#endif

#define DATIMEPRIVEE
#include "strutil.h"
#include "nbconst.h"
#include "datime.h"


/*       Fonctions de calcule de date et temps (Heure)       */
/*! 
  \ingroup NTools
  \file datime.c
  \brief Set of C functions and macros to manipulate date and time  + Sidereal time
  and simple astronomical functions.

  \warning  If possible, use class SOPHYA::TimeStamp in module BaseTools
  for date/time computation.

  Ce groupe de fonctions permettent de manipuler des dates et heures.
  En particulier, il est possible de calculer l'ecart (en nombre de jours
  ou de secondes separant deux dates, ou le temps sideral correspondant 
  a une date et heure legale. Deux structures simples sont definies
  afin de faciliter le passage des arguments entre differentes fonctions:
  - JMA : Jour, Mois, Annee
  - HMS : Heure, Minutes, Secondes

  \sa SOPHYA::TimeStamp
*/

/*
++
  StrgtoJMA(strg, jma)
        Decodage d'une chaine de caracteres "strg" sous forme de "12/3/97" en structure 
        JMA "jma".
  JMAtoStrg(jma, strg)
        Ecrit le contenu de la structure JMA "jma" sous forme de "jj/mm/aaaa" dans la chaine
        de caracteres "strg"
  JMAtoStrgLong(jma, strg)
        Ecriture de date en format long , par ex : "Jeudi , 19 Juin 1997"
  StrgtoHMS(strg, hms)
        Decodage d'une chaine de caracteres "strg" sous forme de "hh:mm:ss 10:43:60.5"
        en structure HMS "hms".
  HMStoStrg(hms, strg)
        Ecrit le contenu de la structure HMS "hms" sous forme de "hh:mm:ss" dans la chaine
        de caracteres "strg"
  StrtoHMS(strg, hms)
        Decodage d'une chaine de caracteres "strg" sous forme de "hh:mm:ss 10:43:60.5"
        en structure HMS "hms" avec gestion des signes (decodage correcte de -44:30:05.23)
  HMStoStr(hms, strg)
        Ecrit le contenu de la structure HMS "hms" sous forme de "hh:mm:ss" dans la chaine
        de caracteres "strg" avec gestion des signes
--
*/
/*!
  \ingroup NTools
   Decodage d'une chaine de caracteres "strg" sous forme de "hh:mm:ss 10:43:60.5"
   en structure HMS "hms" avec gestion des signes (decodage correcte de -44:30:05.23)
*/
void StrtoHMS(char *s,HMS* h)
/* On ne peut pas ecrire 1:-05:-45 pas gere et debile! (mn et sec >=0.)  cmv 12/8/97 */
{
int imoins, i2pt;
h->Heures = h->Minutes = 0; h->Secondes = 0.;
/* decodage de la chaine de characteres */
sscanf(s,"%d:%d:%lf", &(h->Heures),&(h->Minutes),&(h->Secondes));
/* cas du type 5:34:45.5 */
if( h->Heures > 0 ) return;
/* cas du type -5:34:45.5 */
if( h->Heures < 0 ) {
  h->Minutes *= -1; h->Secondes *= -1.;
  return;
}
/* Ici les Heures sont donc nulles "0:..." ou "-0:..." */
/* on n'a donne que l'heure nulle cas "0" ou "-0" sans ":" (idiot mais bon!) */
i2pt = posc(s,':');
if( i2pt < 0 ) return;
/* Ici on a au moins un ou 2 ":" */
/* cas du type 0:34:45.5 ou 0:0:45.5 sans signe - */
imoins = posc(s,'-');
if( imoins < 0 ) return;
/* cas du type -0:34:45.5 ou -0:0:45.5, le signe - n'est pas decode par sscanf */
if( imoins < i2pt ) {   /* on a donne -0:34 ou -0:34:45.5 */
  h->Minutes *= -1; h->Secondes *= -1.;
  return;
}
/* Ici il reste les cas mal codes ou le signe "-" est n'importe ou! */
/* On laisse le decodage par sscanf */
return;
}

/*!
  \ingroup NTools
  Ecrit le contenu de la structure HMS "hms" sous forme de "hh:mm:ss" dans la chaine
  de caracteres "strg" avec gestion des signes
*/
void HMStoStr(HMS h,char *s)
/* L'ecriture est forcee a h:mn:sec avec mn et sec >=0  cmv 12/8/97 */
{
double v;
v = HMStoH(h);
h =  DoubletoHMS(v);
if( v < 0. ) {
  h.Heures *= -1; h.Minutes *= -1; h.Secondes *= -1.;
  sprintf(s,"-%02d:%02d:%04.1f", h.Heures,h.Minutes,h.Secondes);
} else {
  sprintf(s,"%02d:%02d:%04.1f", h.Heures,h.Minutes,h.Secondes);
}
}

/* Nouvelle-Fonction */
/*!
  \ingroup NTools
  \brief Retourne le nombre de jours dans le mois m, annee a
*/
int NbJourMois(int a, int m)
/* Retourne le nombre de jours dans le mois m, annee a */
{
if(a<100) a+=1900;
if( m<1 || m>12 ) return(-1);
m--;
if(m!=1) return(NbJoMois[m]);
/* traitement du mois de fevrier bissextile ou non */
if(  ((a%4==0) && (a%100!=0))
  ||       (a%400==0)         ) return(NbJoMois[m]+1);
else return(NbJoMois[m]);
}

/*!
  \ingroup NTools
  Calcule le Nombre de jours ecoules depuis 0 Jan 1901. 
  Si annee < 100  On considere annee = annee+1900  (1900-1999) 
*/
/* Nouvelle-Fonction */
long JMAtoJ(JMA jma)
{
long rc,nban;
int i;

/* Protection */
if ( (jma.Mois < 1) || (jma.Mois > 12) )  jma.Mois = 1;
if ( (jma.Jour < 1) || (jma.Jour > 31) )  jma.Jour = 1;

if (jma.Annee < 100)  jma.Annee+=1900;
nban = jma.Annee - 1901;

/*  Annee Normale = 365 Jours - Bissextile = 366 Jours    */
/*  366 Jours  :  Divisible par 4 et non divisible par 100  (1972,1976)   */
/*             :  Divisible par 400 (1900 Non-Biss. et 2000 Biss.)        */

if (nban >= 0)  
  rc = nban*365 + (nban/4) - (nban/100) + ((nban+300)/400);
else 
  rc = nban*365 + (nban/4) - (nban/100) + ((nban-300)/400);

/*    Comptabilisation du nb. des jours des mois ecoule   */
for (i=1; i<jma.Mois; i++)  rc += NbJourMois(jma.Annee, i);
rc += jma.Jour;

return(rc);
}

/*!
  \ingroup NTools
  \brief retourne la date correspondant a \b j jours ecoules depuis le 0/1/1901
*/
/* Nouvelle-Fonction */
JMA JtoJMA(long j)
{
long i;
int inc,m;
JMA buf;

inc = (j>=0) ? 1 : -1;
buf.Jour = buf.Mois = 1; buf.Annee = 1901;

/* recherche de l'annee */
for(;;) {
  i = j-JMAtoJ(buf);
  if( i==0 ) return(buf);
  if( inc*i<0 ) {
    /* printf("...annee j=%ld scan=%ld test=%ld\n"
          ,j,JMAtoJ(buf),inc*i); */
    break;
  }
  buf.Annee += inc;
}
if(inc>0) buf.Annee--;
/* printf("Annee trouvee: %d\n",buf.Annee); */

/* recherche du mois */
for(m=1;m<=12;m++) {
  buf.Mois=m;
  i = j-JMAtoJ(buf);
  if( i==0 ) return(buf);
  if(i<0) {
    /* printf("...mois j=%ld scan=%ld test=%ld\n"
          ,j,JMAtoJ(buf),i); */
    buf.Mois--;
    break;
  }
}
/* printf("Mois trouve: %d\n",buf.Mois); */

/* recherche du jour */
for(m=1;m<=31;m++) {
  buf.Jour=m;
  i = j-JMAtoJ(buf);
  if( i==0 ) return(buf);
}
/* printf("Jour trouve: %d\n",buf.Jour); */

/* on a pas trouve !!! normallement impossible! */
printf("**** JtoJMA_Erreur: date non trouvee j=%ld ****\n",j);
buf.Jour = buf.Mois = 1; buf.Annee = 0;
return(buf);
}

/*!
  \ingroup NTools
  \brief Cette fonction calcule le numero de jour de la semaine (Lundi=1 .. Dimanche=7) 
*/
/* Nouvelle-Fonction */
int NumJour(JMA jma)
{
int l;

/*  On appelle JMAtoJ() pour avoir le nb de jours ecoule depuis 1/01/1901  */
/*  Le premier janvier 1901 etait un Mardi (JMAtoJ = 1 , NumJour = 2)      */
l = JMAtoJ(jma)%7 ;
if (l >= 0)  return(l+1);
else return(l+8);
}

/* Nouvelle-Fonction */
JMA JsmmaatoJMA(int nj, int mm, int aa)

/*   Pour "nj=kj*10+jj", calcule la date correspondant au "kj" eme 
     jour "jj" (1 .. 7 - Lundi .. Dimanche) du mois mm , annee aa
     Si "nj" est negatif ->  Dernier jour "jj". Si "kj=0", premier jour "jj".
     Si "kj" ne peut etre atteint -> Dernier jour "jj".
     Si annee "aa" <= 99   ->  annee de 1901 - 1999       */

{
JMA  jma;
int i,k;
int jj, kj;

kj = nj/10;
jj = nj%10;
if (kj < 0)  kj = 0;

if (mm > 12)  mm = 12;
if (mm < 1 ) mm = 1;
if (aa < 100)  aa += 1900;
jma.Mois = mm;
jma.Annee = aa;

if (kj == 0) {
  if (jj >= 0)  {      /*   Premier jour du mois   */    
    if (jj > 7)  jj = 7;
    if (jj < 1)  jj = 1;
    jma.Jour = 1;
    i = NumJour(jma);
    if (i > jj)  jma.Jour = (7-i)+jj+1;
    else jma.Jour = 1+(jj-i);
    }
  else   {             /*   Dernier jour du mois   */
    jj = -jj;
    if (jj > 7)  jj = 7;
    if (jj < 1)  jj = 1;
    jma.Jour = NbJourMois(aa,mm);  
    i = NumJour(jma);
    if (i >= jj)   jma.Jour -= (i-jj);
    else  jma.Jour -= (7+i-jj);
  }
}
else {  /*  On demande le kj eme jour jj du mois  */
  if (jj > 7)  jj = 7;
  if (jj < 1)  jj = 1;
  k = 0;
  for(i=1; i<=NbJourMois(aa,mm); i++) {
    jma.Jour = i;
    if (NumJour(jma) == jj)  k++;
    if (k == kj)  break;
  }
/*  On n'a pas kj jour jj ds ce mois, on renvoie le dernier jj */
  if (k<kj) return(JsmmaatoJMA(-jj, mm, aa));  
}

return(jma);
}

/*!
  \ingroup NTools
  Pour "nj=dd*10+jj", calcule la date correspondant au jour
  "jj" (1 .. 7 - Lundi .. Dimanche) du mois "mm" , annee "aa"
  se trouvant apres (ou egal a) la date "dd/mm/aa" .
  Si annee "aa" <= 99   ->  annee de 1901 - 1999
*/
JMA JApmmaatoJMA(int nj, int mm, int aa)
{
int i, dd, jj;
JMA jma;

if (nj < 0) nj = -nj;
if (nj < 1) nj = 1;
dd = nj/10;
jj = nj%10;
for(i=1; i<=6; i++) {
  jma = JsmmaatoJMA(jj+10*i, mm, aa);
  if(jma.Jour >= dd) return(jma);
}
printf("JApmmaatoJMA(%d...) Erreur / Pas de jour %d >= %d/%d/%d\n", nj, jj,dd, mm, aa);
jma.Jour = jma.Mois = jma.Annee = 0;
return(jma);
}

/*!
  \ingroup NTools
  Cette fonction calcule le nombre d'heures en decimales pour \b hms. 
  Heures/Minutes/Secondes peuvent etre +/- avec toutes les combinaisons possibles
*/
/* Nouvelle-Fonction */
double HMStoH(HMS hms)

/* Cette fonction calcule le nb d'heures en decimales pour hms */
/* Heures/Minutes/Secondes peuvent etre +/- avec toutes les combi possibles */
{
return ((double)hms.Heures + ((double)hms.Minutes)/60. + hms.Secondes/3600.);
}


/*!
  \ingroup NTools
  \brief Conversion en secondes de la structure HMS  "hms"
*/
/* Nouvelle-Fonction */
double HMStoSec(HMS hms)
/*    Calcul du nb de secondes pour hms    */
/* Heures/Minutes/Secondes peuvent etre +/- avec toutes les combi possibles */
{
return ((double)hms.Heures*3600. + (double)hms.Minutes*60. + hms.Secondes);
}

/*!
  \ingroup NTools
  Calcule le nombre de secondes correspondant a "date,heure" 
  a partir de l'origine 1er Janv 1990  0H00 ( Date entre 1930 - 2050 )
*/
/* Nouvelle-Fonction */
long DatetoSec(char const* date, char const* heure)
/*   Calcul du nombre de secondes correspondant a date/heure    */
/*  a partir de l'origine 1er Janv 1990  0H00                   */
/*  Date entre 1930 - 2050                                      */
{
JMA jma,jma0;
HMS hms;
long dj;

StrgtoJMA(date,jma);   StrgtoHMS(heure,hms);
StrgtoJMA("01/01/1990",jma0);
dj = JMAtoJ(jma) - JMAtoJ(jma0);
if (dj < -24000)  return(-0x7FFFFFFF);
if (dj > 24000)  return(0x7FFFFFFF);
return( dj*3600*24 + (int)(HMStoSec(hms)) );    
}

/*!
  \ingroup NTools
  Calcule le nombre de secondes correspondant a "date,heure" 
  a partir de l'origine 1er Janv 1990  0H00 ( Date entre 1930 - 2050 )
  en tenant compte des decalages des heures legales.
*/
/* Nouvelle-Fonction */
long DatetoSecOff(char const* date, char const* heure)
{
JMA jma,jma0;
HMS hms, hms0;
long dj,tdec;

StrgtoJMA(date,jma);   StrgtoHMS(heure,hms);
StrgtoJMA("01/01/1990",jma0);   StrgtoHMS("00:00:00",hms0);
dj = JMAtoJ(jma) - JMAtoJ(jma0);
if (dj < -24000)  return(-0x7FFFFFFF);
if (dj > 24000)  return(0x7FFFFFFF);

/* H_vrai+TOff=H_leg : H_vrai-H_vrai0 = (H_leg-TOff)-(H_leg0-TOff0) */
tdec  = 3600*(TLegOffset(jma, hms)-TLegOffset(jma0, hms0));

dj = dj*3600*24 + (int)(HMStoSec(hms)) - tdec;
return(dj);    
}

/*!
  \ingroup NTools
  \brief Conversion heures en decimales en structure HMS
*/
/* Nouvelle-Fonction */
HMS HtoHMS(double h)
/* Heures decimales en H , M , Sec  Il y a un modulo 24  */
/* la sortie est telle que 0 <= heures < 24 */
{
HMS hms; 

while (h < 0.) h += 24.;
  
hms.Heures = (int) h;
h = (h-hms.Heures)*60.;
hms.Minutes = (int) h;
hms.Secondes = (h-hms.Minutes)*60.;
hms.Heures = hms.Heures % 24;

return(hms);
}

/*!
  \ingroup NTools
  Conversion degres decimaux en structure HMS (deg,min,sec)
  Degres decimaux en D , M , Sec  Il y a un modulo 360.
  la sortie est telle que -180 < deg <= 180.
*/
/* Nouvelle-Fonction */
HMS DtoDMS(double h)
{
HMS hms;
int sgn = 1;

while (h > 180.) h -= 360.;
while (h <= -180.) h += 360.;
if( h < 0. ) {h *= -1.; sgn = -1;}
  
hms.Heures = (int) h;
h = (h-hms.Heures)*60.;
hms.Minutes = (int) h;
hms.Secondes = (h-hms.Minutes)*60.;
hms.Heures = hms.Heures % 360;

if( sgn<0 ) {hms.Heures *= -1; hms.Minutes *= -1; hms.Secondes *= -1.;}
return(hms);
}

/*!
  \ingroup NTools
  Conversion degres/heures decimaux en structure HMS (deg/heure,min,sec)
  SANS modulo 360/24
*/
/* Nouvelle-Fonction */
HMS DoubletoHMS(double h)
/*    Degres/Heures decimaux en D/H , M , Sec  SANS modulo 360/24  */
{
HMS hms;

int sgn = 1;
if( h < 0. ) {h *= -1.; sgn = -1;}

hms.Heures = (int) h;
h = (h-hms.Heures)*60.;
hms.Minutes = (int) h;
hms.Secondes = (h-hms.Minutes)*60.;
hms.Heures = hms.Heures;

if( sgn<0 ) {hms.Heures *= -1; hms.Minutes *= -1; hms.Secondes *= -1.;}
return(hms);
}

/*
++
 Titre  Calcule sur dates et heures
--
*/

/*
++
int NbJourMois(int a, int m)
        Retourne le nombre de jours dans le mois "m", annee "a".
long JMAtoJ(JMA jma)
        Calcule le Nb. de jours ecoules depuis 0 Jan 1901, 
        Si annee < 100  On considere annee = annee+1900  (1900-1999)
JMA JtoJMA(long j)
        Retourne la date correspondant a un nombre de jours "j"
        ecoules depuis le 0 Jan 1901.
int NumJour(JMA jma)
        Calcule le numero du jour de la semaine correspondant a la date indiquee
        par la structure JMA "jma". (Lundi=1 .. Dimanche=7)
JMA JsmmaatoJMA(int nj, int mm, int aa)
        Pour "nj=kj*10+jj", calcule la date correspondant au "kj" eme 
        jour "jj" (1 .. 7 - Lundi .. Dimanche) du mois mm , annee aa
        Si "nj" est negatif ->  Dernier jour "jj". Si "kj=0", premier jour "jj".
        Si "kj" ne peut etre atteint -> Dernier jour "jj".
        Si annee "aa" <= 99   ->  annee de 1901 - 1999       
JMA JApmmaatoJMA(int nj, int mm, int aa)
        Pour "nj=dd*10+jj", calcule la date correspondant au jour
        "jj" (1 .. 7 - Lundi .. Dimanche) du mois "mm" , annee "aa"
        se trouvant apres (ou egal a) la date "dd/mm/aa" .
        Si annee "aa" <= 99   ->  annee de 1901 - 1999       
double HMStoH(HMS hms)
        Conversion en heures decimales de la structure  "hms"
double HMStoSec(HMS hms)
        Conversion en secondes de la structure HMS  "hms"
HMS HtoHMS(double h)
        Conversion heures en decimales en structure HMS
HMS DtoDMS(double h)
        Conversion degres decimaux en structure HMS (deg,min,sec)
HMS DoubletoHMS(double h)
        Conversion degres/heures decimaux en structure HMS (deg/heure,min,sec)
                   SANS modulo 360/24
long DatetoSec(char const* date, char const* heure)
        Calcule le nombre de secondes correspondant a "date,heure" 
        a partir de l'origine 1er Janv 1990  0H00 ( Date entre 1930 - 2050 )            
long DatetoSecOff(char const* date, char const* heure)
        Calcule le nombre de secondes correspondant a "date,heure" 
        a partir de l'origine 1er Janv 1990  0H00 ( Date entre 1930 - 2050 )
        en tenant compte des decalages des heures legales.
--
*/



/*   Variables contenant les informations pour conversion  */
/*  Temps Legal <> TU <> Temps solaire moyen               */

static double  TSolMOff = 0.;     /*  TSol.Moyen - TU  (en heures) */
static char SiteName[128];           /*  Nom de site  */

static struct
  {
  char fg;         /*  Flag (=A ou B) Type de specification  */
  int jj[2];       /*  Numero ou date du jour                */
  int mm[2];       /*  Mois    */
  int Off[3];      /*  Decalage = TLegal - TU    */
  }  TLegInfo = {'\0',{0,0},{0,0},{0,0,0}} ;

static HMS HChgH = {0,0,0.} ; /*  Heure de changement d'heure */

/* 
++ 
  Module        Temps Sideral, Universel (C)
  Lib   LibsUtil
  include       datime.h

        Ce groupe de fonctions permet convertir differentes heures:
        les temps local, le temps universel (TU), le temps legal (TLeg). Il est 
        aussi possible de calculer le temps sideral (TS). 
--
*/
/*
++
  Links Voir aussi:
  Dates (C)
--
*/
/*
++
 Titre  T-Legal, TU
        Conversion entre l'heure legale, le temps universel (TU) et le temps solaire
        moyen
--
*/

/*
++
SetTSolMOff(double tsmoff, char *name)
        Definit pour un site donne, le decalage T Sol. Moyen - TU (en heures) "tmsoff"
        et le nom du site 
PrtTSolMOff()
        Imprime les infos definis par "SetTSolMOff()"
--
*/

/* Nouvelle-Fonction */
void SetTSolMOff(double tsmoff, char *name)

/*     Offset  T Sol. Moyen - TU (en hh:mm:ss) ,  Nom de site name   */
{
TSolMOff = tsmoff;
strcpy(SiteName,name);
return;
}


/* Nouvelle-Fonction */
void PrtTSolMOff(void)
{
printf("  **Site : %s , T Sol.Moy.-TU  %g  (heures) \n",SiteName,TSolMOff);
return;
}

/*
++
int SetTLegOff(char* dtz, char* dhz)
        Initialise la structure TLegInfo. L'argument *dtz peut etre 
        sous 2 forme differentes et permet de  definir 3 plages 
        de temps avec 3 decalages horaires (nb entier). 
        Les dates pivots sont les jours de changement d'heure a 
        l'heure indiquee par dhz ("hh:mm:ss"), si "dhz==NULL ou 
        chaine vide, alors changement a 0 heure.  
|   1) dtz = "A Off1 j1/m1 Off2 j2/m2 Off3"
|      jx/mx Jour/Mois - Offx : Decalage TLeg - TU                  
|      pour 01/01 <= jj/mm <  j1/m1   ->  Off1 
|           j1/m1 <= jj/mm <  j2/m2   ->  Off2 
|           j2/m2 <= jj/mm <= 31/12   ->  off3 
|   2) dtz = "B off1 Js1/m1 Off2 Js2/m2 Off3"
|      Jsx/mx  Numero du jour de la semaine (Js) et numero du mois
|      Js =  1 ..  7  (Lundi .. Dimanche)  :  1er Jour Js du mois    
|         = -1 .. -7   Dernier jour Js du mois                  
|         = k1 .. k7   k eme jour Js du mois                  
|   3) dtz = "C off1 Js1/m1 Off2 Js2/m2 Off3"
|      Jsx/mx  Numero du jour de la semaine (Js) et numero du mois avec une 
|      precision d'une date butee inferieurs    
|      Js =  ddj ..  (j=1..7  Lundi .. Dimanche)  :  1er Jour j du mois
|      apres (ou egal a) la date dd/mm    
        Retour : 0 si Ok -1 Sinon                                          
void PrtTLegInfo(int aa)        
        Imprime les informations de decalage TLeg ,  TU - Annee "aa"
int  TLegOffset(JMA date, HMS heure)
        Indique le decalage (en nombre d'heures) TLeg-TU        
--
*/

/* Nouvelle-Fonction */
int SetTLegOff(char* dtz, char* dhz)

                        
{
int is;

HChgH.Heures = HChgH.Minutes = 0;
HChgH.Secondes = 0.;

if (dhz) StrgtoHMS(dhz, HChgH);
if ( (HChgH.Heures < 0) || (HChgH.Heures > 23) )
  { HChgH.Heures = HChgH.Minutes = 0;  HChgH.Secondes = 0.; }

is = sscanf(dtz,"%c %d %d/%d %d %d/%d %d", &(TLegInfo.fg), &(TLegInfo.Off[0]), &(TLegInfo.jj[0]),  &(TLegInfo.mm[0]), &(TLegInfo.Off[1]), &(TLegInfo.jj[1]),  &(TLegInfo.mm[1]), &(TLegInfo.Off[2]) );

if ( (TLegInfo.fg != 'A') && (TLegInfo.fg != 'B')
     && (TLegInfo.fg != 'C') )  TLegInfo.fg = '\0';
if (is != 8) TLegInfo.fg = '\0';
if (TLegInfo.fg != '\0') return(0);
else { printf("SetTLegOff()/Erreur %s \n", dtz); return(-1); }
}

/* Nouvelle-Fonction */
void PrtTLegInfo(int aa)

/*    Cette fonction  print les infos de TLegInfo     */
{
int i,k,j;
char ForL[2][128];
JMA jma;

if (TLegInfo.fg == '\0')    printf ("  !!! TlegInfo not set or Error \n");


HMStoStrg(HChgH, ForL[0]);
printf("  ==== TLegInfo : Changement d'heure a %s \n", ForL[0]);
ForL[0][0] = '\0';
if (TLegInfo.fg == 'A')  
  {
  printf("  ---- TLegInfo : Selection par date (Flag = A)  ---- \n");
  printf("   -> 01 Janvier - %d %s   :  %d heures \n",TLegInfo.jj[0], NomMo[TLegInfo.mm[0]-1], TLegInfo.Off[0]);
  printf("   -> %d %s - %d %s  :  %d heures \n",TLegInfo.jj[0], NomMo[TLegInfo.mm[0]-1],TLegInfo.jj[1], NomMo[TLegInfo.mm[1]-1],TLegInfo.Off[1]);
  printf("   -> %d %s - 31 Decembre  :  %d heures \n",TLegInfo.jj[1], NomMo[TLegInfo.mm[1]-1], TLegInfo.Off[2]);
  }

if (TLegInfo.fg == 'B')  
  {
  printf("  --- TLegInfo : Selection par Jour de semaine et mois (Flag = B)  --- \n");
  for (i=0;i<2;i++) {
    k = TLegInfo.jj[i]/10;
    if (k < 1)  k = 1;
    j = TLegInfo.jj[i]%10;
    if (j < 0)  j = -j;
    if (j < 1) j = 1;
    if (j > 7) j = 7;
    if (TLegInfo.jj[i] < 0)  sprintf(ForL[i]," Dernier  %s  de  %s  (KJ=%d)", NomJo[j-1], 
                                     NomMo[TLegInfo.mm[i]-1], TLegInfo.jj[i]);
    else  sprintf(ForL[i]," %d eme  %s  de  %s (KJ=%d)", k, NomJo[j-1], 
                  NomMo[TLegInfo.mm[i]-1], TLegInfo.jj[i]);
    }
  printf(" Butees: %s  -  %s \n", ForL[0], ForL[1]);
  for (i=0;i<2;i++) {
    jma = JsmmaatoJMA(TLegInfo.jj[i], TLegInfo.mm[i], aa);
    JMAtoStrg(jma, ForL[i]);
  }
  printf("   -> 01 Janvier %d - %s  :  %d heures \n",aa, ForL[0], TLegInfo.Off[0]);
  printf("   -> %s - %s  :  %d heures \n",ForL[0], ForL[1], TLegInfo.Off[1]);
  printf("   -> %s - 31 Decembre %d :  %d heures \n",ForL[1], aa, TLegInfo.Off[2]);
 }

if (TLegInfo.fg == 'C')  
  {
  printf("  --- TLegInfo : Selection par Jour de semaine >= dd et mois (Flag = C)  --- \n");
  for (i=0;i<2;i++) {
    k = TLegInfo.jj[i]/10;
    j = TLegInfo.jj[i]%10;
    if (j < 0)  j = -j;
    if (j < 1) j = 1;
    if (j > 7) j = 7;
    sprintf(ForL[i],"1er %s  apres le %d %s ", NomJo[j-1], k, NomMo[TLegInfo.mm[i]-1]);
    }
  printf(" Butees: %s - %s \n", ForL[0], ForL[1]);
  for (i=0;i<2;i++) {
    jma = JApmmaatoJMA(TLegInfo.jj[i], TLegInfo.mm[i], aa);
    JMAtoStrg(jma, ForL[i]);
  }
  printf("   -> 01 Janvier %d - %s  :  %d heures \n",aa, ForL[0], TLegInfo.Off[0]);
  printf("   -> %s - %s  :  %d heures \n",ForL[0], ForL[1], TLegInfo.Off[1]);
  printf("   -> %s - 31 Decembre %d :  %d heures \n",ForL[1], aa, TLegInfo.Off[2]);
 }

return;
}


/* Nouvelle-Fonction */
int  TLegOffset(JMA date, HMS heure)

/*     Cette fonction calcule l'offset   TLeg - TU  pour la date donnee date  */
{
JMA jma1,jma2;
long nbj1,nbj2,nbj;
double decj1, decj2, decj;
int Off;

/*   On recherche l'offset correspondant a la date   */
jma1.Annee = jma2.Annee = date.Annee;
if (TLegInfo.fg == 'A')
/*    Periode defini par des dates ds les mois   */ 
  {
  jma1.Mois = TLegInfo.mm[0];
  jma1.Jour = TLegInfo.jj[0];
  jma2.Mois = TLegInfo.mm[1];
  jma2.Jour = TLegInfo.jj[1];
  }
if (TLegInfo.fg == 'B')
/*    Periodes definies par jour de semaine et mois   */
  {
  jma1 = JsmmaatoJMA(TLegInfo.jj[0],TLegInfo.mm[0],date.Annee);
  jma2 = JsmmaatoJMA(TLegInfo.jj[1],TLegInfo.mm[1],date.Annee);
  }

if (TLegInfo.fg == 'C')
/*    Periodes definies par jour de semaine et mois  avec date inferieure */
  {
  jma1 = JApmmaatoJMA(TLegInfo.jj[0],TLegInfo.mm[0],date.Annee);
  jma2 = JApmmaatoJMA(TLegInfo.jj[1],TLegInfo.mm[1],date.Annee);
  }

nbj = JMAtoJ(date);
nbj1 = JMAtoJ(jma1);
nbj2 = JMAtoJ(jma2);

/*  Les dates pivots sont les jours de changement d'heure a Heure= HChgH */
decj = (double)nbj + HMStoH(heure)/24.;
decj1 = (double)nbj1 + HMStoH(HChgH)/24.;
decj2 = (double)nbj2 + HMStoH(HChgH)/24.;

if ( decj < decj1 )  Off = TLegInfo.Off[0];
else if (decj < decj2)  Off = TLegInfo.Off[1];
     else Off = TLegInfo.Off[2];

/*     Impression de debugging   */
/*Debug--  printf ("   @@@ Debug_TLegOffset (date= %d/%d/%d)  :  %d H. \n",date.Jour,date.Mois,date.Annee,Off);
printf("  @-1 %s , %d %s %d \n", NomJo[NumJour(jma1)-1], jma1.Jour, NomMo[jma1.Mois-1], jma1.Annee) ;
printf("  @-2 %s , %d %s %d \n", NomJo[NumJour(jma2)-1], jma2.Jour, NomMo[jma2.Mois-1], jma2.Annee) ;  --*/


return (Off);
}

/*
++
HMS TLegtoTU(JMA date, HMS TLeg)
        Conversion TLegal en TU pour une date donnée
HMS TUtoTLeg(JMA date, HMS TLeg)
        Conversion TU en TLegal pour une date donnée. L'heure légale calculée
        peut être erronnée autour des dates de changement d'heure.
HMS TUtoTSolM(HMS TU)
        Conversion TU en temps solaire moyen
--
*/

/* Nouvelle-Fonction */
HMS TLegtoTU(JMA date, HMS hmsz)

/*   Cette fonction transforme le temps legal en temps universel   */
{
HMS  TU;
int Off;

Off = TLegOffset(date, hmsz);
TU = hmsz;
TU.Heures -= Off;
if (TU.Heures >= 24)  TU.Heures -= 24;
if (TU.Heures < 0)  TU.Heures += 24;
return(TU);
}



/* Nouvelle-Fonction */
HMS TUtoTLeg(JMA date, HMS hms)

/*   Cette fonction transforme le temps universel en temps legal */
{
HMS  TLeg;
int Off;

Off = TLegOffset(date, hms);    /*  $CHECK$  */
TLeg = hms;
TLeg.Heures += Off;
TLeg.Heures = TLeg.Heures%24;
if (TLeg.Heures < 0)  TLeg.Heures += 24;
return(TLeg);
}


/* Nouvelle-Fonction  */
HMS TUtoTSolM(HMS TU)

/*   Passage  TU  ->  Temps local    */
{
double t;
t = HMStoH(TU) + TSolMOff;
return(HtoHMS(t));
}

/*
++
 Titre  Temps Sideral
|      Jour Sideral = 86400 Sec. Sid.  ( = 24 H Sid. )    
|      Annee Solaire = 365 J 5H 48' 46"  =~  365.2422 J                       
|      Annee Solaire = 365.2422 Jours Solaire Moyen = 366.2422 Jours Sideraux 
|      Jour Sol. Moyem = 24 H Sol. = 24 H Sid. * (366.2422 / 365.2422)        
|      Jour Solaire Moyen = 86636.555359 Sec. Sid. = 24 H 3' 56.555359"       
|      Heure Solaire Moyen = 3609.85647 Sec. Sid.                             
|      Sec. Sol. = 1.0027379 Sec. Sid.                                        
|      Le 1er Janvier 1989 0H TU Il etait 6H 42' 30" GMT_Sideral 
--
*/

/*
++
HMS TLegtoTSid(JMA date, HMS TZ)
        Cette fonction calcule le temps sideral a partir du temps
        legal et de la date
GMST_at_0h_UT(JMA date)
        Cette fonction calcule le temps sideral moyen a greenwich a 0h UT
--
*/

/* Nouvelle-Fonction  */
HMS TLegtoTSid(JMA date, HMS TZ)

/*  Cette fonction  donne le temps sideral a partir du temps legal  */
/*  et la date                                                      */
{
HMS  Orgh,TU,TSid;
JMA  Orgd;
int off;
long lj;
double ts,ratio;

/* Le 1er Janvier 1989 0H TU Il etait 6H 42' 30" GMT_Sideral */
StrgtoJMA("1/1/1989",Orgd);
StrgtoHMS("6:42:30",Orgh);



/*    Nb de jours / Date Origine   */
lj = JMAtoJ(date) - JMAtoJ(Orgd);

/*  On calcule l'heure TU et on tient compte du decalage eventuel de + ou - 1 jour  */
off = TLegOffset(date, TZ);
TU = TZ;
TU.Heures -= off;
if (TU.Heures < 0)  
  {
  lj -= 1;
  TU.Heures += 24;
  }
if (TU.Heures >= 24)  
  {
  lj += 1;
  TU.Heures -= 24;
  }



ratio = ((365.*24.) + 5. + (48. / 60.) + (46. / 3600.)) / 24.;
ratio = (ratio+1) / ratio ;

/*   Calcul du decalage en heures du temps Sideral par rapport a          */
/*   l'origine choisie :  ts                                              */

ts = lj * (24. * (ratio - 1) );
ts += ( TU.Heures * ratio );
ts += ( ( TU.Minutes / 60. ) * ratio );
ts += ( (TU.Secondes / 3600. ) * ratio );

ts += HMStoH(Orgh);

/*   On tient compte de l'offset du temps solaire moyen  */
ts += TSolMOff;

TSid = HtoHMS(ts);
TSid.Heures = TSid.Heures % 24;
if (TSid.Heures < 0)  TSid.Heures += 24;

return(TSid);
}

/*!
  \ingroup NTools
  \brief Cette fonction  donne le temps sideral moyen a greenwich a 0h UT 
*/
/* Nouvelle-Fonction  */
HMS GMST_at_0h_UT (JMA date)

/*  Cette fonction  donne le temps sideral moyen a greenwich a 0h UT */
/*                                                       CMV 5/11/93 */
{
HMS Orgh, TSid;
JMA  Orgd;
long lj, ejs;
double js, hs, hs0;

/*  J H ' " sont les vraies valeurs indiscutables (cf metre etalon etc..)  */
/*  Annee Solaire  = 365 J 5H 48' 45.2" =~  365.2421898148 J               */
/*  Annee Siderale = 365 J 6H  9'  9.8" =~  365.2563634259 J               */
/*  Le 1er Janvier 1994 0H TU Il sera 6h 41' 40.4305" GMT_Sideral          */
StrgtoJMA("1/1/1994",Orgd);
StrgtoHMS("6:41:40.4305",Orgh);
hs0 = HMStoH(Orgh);

/*   Nb de jours ecoules depuis la Date Origine */
lj = JMAtoJ(date) - JMAtoJ(Orgd);
/* printf("nombre de jours ecoules depuis 1/1/94 %d\n",lj); */

/*   Nb de jours sideraux depuis la Date Origine */
js = (double) lj * 366.2421898148 / 365.2421898148;
/*   Nb d'heures siderales depuis Date Origine (modulo un jour) */
ejs = (int) js;
hs = js - (double) ejs;
hs *= 24.;
hs += hs0;

TSid = HtoHMS(hs);
/* printf("nombre d'heures siderales depuis la reference (modulo un jour) %d %d %f\n"
         ,TSid.Heures,TSid.Minutes,TSid.Secondes); */

return(TSid);
}

/*
++
HMS TUtoTSid(JMA date, HMS TU)
        Cette fonction calcule le temps sideral GMT a partir du temps
        universel et de la date GMT.
int TSidtoTU(JMA date, HMS TS, HMS *TU1, HMS *TU2)
        Cette fonction calcule les temps universels a partir du temps
        sideral GMT et de la date GMT. Selon la valeur du TS, il peut y avoir
        une ou deux possibilites pour le TU (return code).
--
*/
/*!
  \ingroup NTools
  Cette fonction calcule le temps sideral GMT a partir du temps
  universel et de la date GMT.
  \verbatim
  Temps sideral
  Jour Sideral = 86400 Sec. Sid.  ( = 24 H Sid. )    
  Annee Solaire = 365 J 5H 48' 46"  =~  365.2422 J                       
  Annee Solaire = 365.2422 Jours Solaire Moyen = 366.2422 Jours Sideraux 
  Jour Sol. Moyem = 24 H Sol. = 24 H Sid. * (366.2422 / 365.2422)        
  Jour Solaire Moyen = 86636.555359 Sec. Sid. = 24 H 3' 56.555359"       
  Heure Solaire Moyen = 3609.85647 Sec. Sid.                             
  Sec. Sol. = 1.0027379 Sec. Sid.                                        
  Le 1er Janvier 1989 0H TU Il etait 6H 42' 30" GMT_Sideral 
  \endverbatim
*/
/* Nouvelle-Fonction  */
HMS TUtoTSid(JMA date, HMS TU)

/*  Cette fonction  donne le temps sideral GMT a partir du temps universel  */
/*  et la date  GMT                                                         */
{
HMS  Orgh,TSid;
JMA  Orgd;
long lj;
double ts,ratio;

/*   Le 1er Janvier 1989 0H TU Il etait 6H 42' 30" GMT_Sideral */
StrgtoJMA("1/1/1989",Orgd);
StrgtoHMS("6:42:30",Orgh);

/*    Nb de jours / Date Origine   */
lj = JMAtoJ(date) - JMAtoJ(Orgd);

ratio = ((365.*24.) + 5. + (48. / 60.) + (46. / 3600.)) / 24.;
ratio = (ratio+1) / ratio ;

/*   Calcul du decalage en heures du temps Sideral par rapport a          */
/*   l'origine choisie :  ts                                              */

ts = lj * (24. * (ratio - 1) );
ts += ( TU.Heures * ratio );
ts += ( ( TU.Minutes / 60. ) * ratio );
ts += ( (TU.Secondes / 3600. ) * ratio );

ts += HMStoH(Orgh);

TSid = HtoHMS(ts);
TSid.Heures = TSid.Heures % 24;
if (TSid.Heures < 0)  TSid.Heures += 24;

return(TSid);
}

/*!
  \ingroup NTools
  Cette fonction calcule les temps universels a partir du temps
  sideral GMT et de la date GMT. Selon la valeur du TS, il peut y avoir
  une ou deux possibilites pour le TU (return code).
*/
/* Nouvelle-Fonction  */
int TSidtoTU(JMA date, HMS TS, HMS *TU1, HMS *TU2)

/*  Cette fonction  donne les temps universels a partir du temps sideral GMT */
/*  et la date GMT                                                           */
{
int rc = 0;
HMS  Orgh,T0;
double ratio,delt,deltt,jsecsid;
char strd[16],strh[16];
/* char strh1[16]; */

TU1->Heures = TU2->Heures = -1;
TU1->Minutes = TU2->Minutes = 0;
TU1->Secondes = TU2->Secondes = 0.;

ratio = ((365.*24.) + 5. + (48. / 60.) + (46. / 3600.)) / 24.;
ratio = (ratio / (ratio+1));
jsecsid = 86400./ratio;

/* temps sideral pour la date a 0h TU */
StrgtoHMS("00:00:00",T0);
Orgh =  TUtoTSid(date,T0);
/* JMAtoStrg(date,strd); HMStoStrg(T0,strh); HMStoStrg(Orgh,strh1);
printf("... date=%s tu=%s -> org_ts=%s\n",strd,strh,strh1); */

/* difference de temps sideral / a 0h TU */
delt = HMStoSec(TS) - HMStoSec(Orgh);
/* printf("... ts=%f org_ts=%f en sec -> delt=%f (jsecsid=%g)\n"
      ,HMStoSec(TS),HMStoSec(Orgh),delt,jsecsid); */

if(fabs(delt)<1.e-10) delt=0.;
if(delt<0.) delt += 86400.;
if( fabs(delt-86400.)<1.e-10 ) delt = 86400. - 1.e-10;
if(delt<0. || delt>=86400.+1.e-10) {
  HMStoStrg(TS,strh);
  JMAtoStrg(date,strd);
  printf("*** BUG *** TSidtoTU_erreur: delt=%g<0. (d=%s,ts=%s) *******\n"
        ,delt,strd,strh);
  return(0);
}

/* difference de temps / a 0h TU 1ere solution */
/* printf("... delt=%f en sec sid",delt); */
deltt = delt*ratio;
*TU1 = SectoHMS(deltt);
*TU2 = *TU1;
rc++;
/* HMStoStrg((*TU1),strh1);
printf("... deltt=%f en sec -> TU1=%s\n",deltt,strh1); */

/* Y a t-il une seconde solution ? */
if(delt>=jsecsid-86400.-1.e-10) return(rc);
delt += 86400.;
deltt = delt*ratio;
*TU2 = SectoHMS(deltt);
/* HMStoStrg((*TU2),strh1);
printf("... deltt=%f en sec -> TU2=%s\n",deltt,strh1); */
rc++;

return (rc);
}

/*!
  \ingroup NTools
  Fonction d'initialisation de conversion  de temps sideral et legal 
  pour la ESO-Silla "SetTSolMOff() , SetTLegOff()".
  \verbatim
  Longitude = 70 deg 43.8 min ouest
  Latitude  = 29 deg 15.4 min Sud    
  T Sol. Moyen - TU = -4.715333 Heures
  TLegal-TU : -4 Heures de Mars a Octobre, -3 Heures sinon
  Changement le 1er dimanche >= 8/03  8/10 a midi (12:00:00)
  \endverbatim
*/
/* Fonction d'initialisation de calcul de temps sideral pour la Silla */
void TSidSetupLaSilla()
{
SetTSolMOff((double)(-4.715333),"La Silla (Chili)");
if (SetTLegOff("C -3 87/3 -4 87/10 -3","12:00:00") != 0)
  printf("  TSidSetupLaSilla / Erreur SetTLegOff \n"); 
/*  Position de l'observatoire */
/* SetObsPos("La Silla (Chili)","O 70:43:48","S 29:15:24",2347);  */
return;
}



/*!
  \ingroup NTools
  Calcul du jour Julien pour une date TU.
  Uniquement valable a partir du 15/10/1582 00:00:00.
  
*/
double ToJulianDay(JMA dateTU, HMS hmsTU)
/* Cf Fundamental astronomie, Springer Verlag 2sd ed.  cmv 4/12/98 */
{
double j_dat_0_ut;
long int f,g,a;
f = (dateTU.Mois>=3) ? dateTU.Annee: dateTU.Annee-1;
g = (dateTU.Mois>=3) ? dateTU.Mois: dateTU.Mois+12;
a = 2 - (int)((double)f/100.) + (int)((double)f/400.);
/* Julian date at 0h UT */
j_dat_0_ut = (int)(365.25*(double)f)+(int)(30.6001*(double)(g+1))+dateTU.Jour+a+1720994.5;
return j_dat_0_ut + HMStoH(hmsTU)/24.;
}

/*!
  \ingroup NTools
  Calcul de la date date et l'heure TU a partir d'un jour Julien.
  Retourne 0 si succes.
*/
int FromJulianDay(double JD,JMA* dateTU, HMS* hmsTU)
/* Cf Fundamental astronomie, Springer Verlag 2sd ed.  cmv 4/12/98 */
{
double jd,i,f,a,x,b,c,d,e,h;
if(JD<0.) return 1;

i = JD - (int)(JD);
if(i<0.5) jd = (int)(JD)-0.5; else jd = (int)(JD)+0.5;

i = (int)(jd+0.5);
f = jd+0.5 - i;
if(i<2299161.0) a = i;
else {x = (int)((i-1867216.25)/36524.25); a = i+1.+x-(int)(x/4.);}
b = a+1524.;
c = (int)((b-122.1)/365.25);
d = (int)(365.25*c);
e = (int)((b-d)/30.6001);

x = b-d-(int)(30.6001*e)+f;
dateTU->Jour = (int)(x);

if(e<13.5) dateTU->Mois = (int)(e-1.); else dateTU->Mois = (int)(e-13.);

if(dateTU->Mois>2.5) dateTU->Annee = (int)(c-4716.); else dateTU->Annee = (int)(c-4715.);

hmsTU->Heures = hmsTU->Minutes = 0; hmsTU->Secondes = 0.;
h = (JD-ToJulianDay(*dateTU,*hmsTU)) * 24.;
*hmsTU = HtoHMS(h);

return 0;
}


/*!
  \ingroup NTools
  Cette fonction renvoie la valeur decimale en heures
  d'un angle specifie sous forme de [-][+]dd:mm:ss.
*/
double StrgtoDegDec(char *strg)
/*         cmv 4/12/98 */
{
int dd,mm;
float ss;
char sgn;
double deg;

dd = mm = 0;
ss = 0.;
sgn = *strg;
if (isdigit(sgn))
  sscanf(strg,"%d:%d:%f",&dd,&mm,&ss);
else
  sscanf(strg+1,"%d:%d:%f",&dd,&mm,&ss);
deg = (double)(dd) + (double)(mm)/60. + (double)(ss)/3600.;
if (sgn == '-')  deg = (-deg);
return(deg);
}

/*!
  \ingroup NTools
  Cette fonction ecrit une valeur en degre decimal sous 
  forme de [-]dd:mm:ss.
*/

char * DegDectoStrg(double deg, char *strg)
/*         cmv 4/12/98 */
{
int dd,mm;
float ss;
char sgn;

if (deg < (double)(0.))
  { deg = (-deg) ; sgn = '-' ; }
else sgn = '+';
dd = (int)deg;
deg = (deg-(double)(dd)) * (double)60.;
mm = (int)deg;
ss = (deg-(double)(mm)) * (double)60.;
if (sgn == '-')
  sprintf(strg,"-%d:%d:%f",dd,mm,ss);
else
  sprintf(strg,"%d:%d:%f",dd,mm,ss);
return(strg);
}

/*!
  \ingroup NTools
  Valeur de l'eccentricite de l'orbite terrestre
  a la date TU ``dateTU''.
  Cf Fundamental astronomie, Springer Verlag 2sd ed p477 table E10
*/
double EccEarth(JMA dateTU)
/* Cf Fundamental astronomie, Springer Verlag 2sd ed p477 table E10.  cmv 9/12/98 */
{
double JD,T;
HMS hmsTU;
hmsTU.Heures = 12; hmsTU.Minutes = 0; hmsTU.Secondes = 0.;
JD = ToJulianDay(dateTU,hmsTU);
/* T en siecles juliens depuis 31/12/1899 a 12h TU --> JD=2415020. */
T = (JD-2415020.)/36525;
return(0.0167498 - 0.00004258*T - 0.000000137*T*T);
}

/*!
  \ingroup NTools
  Valeur de l'obliquite de l'orbite terrestre
  a la date TU ``dateTU'' (en degres decimaux).
*/
double ObliqEarth(JMA dateTU)
/* Cf Fundamental astronomie, Springer Verlag 2sd ed.  cmv 9/12/98 */
/* le 1/1/1900 a 0h00 TU --> jd = 2415020.5 */
{
double jd;
HMS hmsTU;
hmsTU.Heures = hmsTU.Minutes = 0; hmsTU.Secondes = 0.;
jd = ToJulianDay(dateTU,hmsTU) - 2415020.5;
return(23.452294 - 0.00013 *jd/365.25);
}

/*! 
  \ingroup NTools
  Retourne la Longitude Ecliptique du perihelie de
  l'orbite terrestre a la date TU ``dateTU'' (en degres decimaux).
*/
double LongEcPerihelie(JMA dateTU)
/* Cf Fundamental astronomie, Springer Verlag 2sd ed p477 table E10.  cmv 9/12/98 */
{
double JD,T;
HMS hmsTU;
hmsTU.Heures = 12; hmsTU.Minutes = 0; hmsTU.Secondes = 0.;
JD = ToJulianDay(dateTU,hmsTU);
/* T en siecles juliens depuis 31/12/1899 a 12h TU --> JD=2415020. */
T = (JD-2415020.)/36525;
/* l = 101d 13' 07.15" + 6171.77" *T + 1.823" *T^2 */
return(101.218652778 + 1.71438055556*T + 0.000506388889*T*T);
}

/*!
  \ingroup NTools
  Renvoie les coordonnees ecliptiques \b (l,b) a partir de
  coordonnees equatoriales \b (a,d) a la date TU \b dateTU.
  La date permet de corriger la variation annuelle
  de l'obliquite terrestre.
*/
void EquatToEclip(double a,double d,double* l,double *b,JMA dateTU)
/* Cf Fundamental astronomie, Springer Verlag 2sd ed p477 table E10.  cmv 9/12/98 */
/* l=[0,360[, b=[-90,90] */
{
double e,sinb,cosb,sinl,cosl;

e = ObliqEarth(dateTU)*DegenRad;
a *= DegenRad; d *= DegenRad;

sinb = sin(d)*cos(e)-cos(d)*sin(e)*sin(a);
if(sinb>1.) sinb=1.; if(sinb<-1.) sinb=-1.;
*b = asin(sinb);
cosb = cos(*b);

*l = 0.;
if(cosb!=0.) {
  sinl = (sin(d)*sin(e)+cos(d)*cos(e)*sin(a))/cosb;
  cosl = cos(d)*cos(a)/cosb;
  if(cosl>1.) cosl=1.; if(cosl<-1.) cosl=-1.;
  *l = acos(cosl);
  if(sinl<0.) *l += M_PI;
  if(*l>2.*M_PI) *l -= 2.*M_PI;
}

*b /= DegenRad;
*l /= DegenRad;
}