//---------------------------------------------------------------------------
#include <vcl\vcl.h>
#pragma hdrstop

#include "calcul_maree.h"
#include "table2NC.h"
#include "trigo.h"

extern double uneminute;		    // une minute exprimée en heures
extern constituentsNC table2NC[nb_const];   // table des constantes harmoniques
extern double ampli[nb_const];		// amplitude de la constante
extern double phase[nb_const];		// phase     de la constante

extern double s;       	            // mean longitude of moon
extern double h;		    	    // mean longitude of sun
extern double p;			        // longitude of lunar perigee
extern double p1;			        // longitude of solar perigee
extern double N;			        // longitude of moon's node

double equilbrm[nb_const];		    // equilibrium argument (theorical phase)
double nodefctr[nb_const];		    // node factor

void   equi_tide();
double heure_maree(double t0);
bool   signe_derivee(double t);
double amplitude(double t);
double amplimax();

//---------------------------------------------------------------------------
void equi_tide()
{ double T=180.;
  String symbole;
  for (int i=0;i<nb_const;i++)
    { equilbrm[i]=table2NC[i].T*T     + table2NC[i].s*s   + table2NC[i].h*h
                 +table2NC[i].p*p     + table2NC[i].p1*p1
                 +(double)table2NC[i].deg ;
      equilbrm[i]+=(*table2NC[i].fnu)();
      equilbrm[i]=reduc360(equilbrm[i]);
      nodefctr[i]=(*table2NC[i].fnf)();// appel fonction pointée par table2NC
    }
}
//---------------------------------------------------------------------------
double heure_maree(double t0)	    // heure de la marée en heures
{ /* la hauteur de la maree est la fonction somme des composants,
     chaque composant étant de la forme acos(vt-p).
     La fonction dérivée =-vsin(vt-p) s'annule lorsque la fonction
     passe par un maxi (pleine mer) ou un mini (basse mer). Elle est
     positive lorsque la marée monte, et négative lorsque la marée
     descend. Plutôt que de rechercher les racines de la fonction dérivée,
     on a trouvé plus simple de procéder par approches successives.
  */

  double t, dh;
  // sens de la maree 1=maree montante 0=maree descendante
  bool sens, sens0;

  dh=0.5;                       // une demi-heure
  sens0=signe_derivee(t0);	    // sens de la marée a l'instant initial
  // si 30 secondes avant, la marée n'était pas dans le même sens, c'est
  // qu' elle a changé de sens entre-temps ! heure marée = instant initial
  if (signe_derivee(t0-0.5*uneminute)!= sens0) return t0;

  t=t0;			                // on part de l'instant initial
  do                            // on va regarder
    { t=t+dh;		            // les demi-heures suivantes
      sens=signe_derivee(t);    // quel est le sens de la marée
    }
  while (sens==sens0);	        // jusqu'à ce qu'elle change de sens
  // elle a changé de sens !
  sens0=sens;	    	        // on note le nouveau sens
  do			                // on revient en arriere
    { t=t-uneminute;	        // minute par minute
      sens=signe_derivee(t);	// tant que le sens
    }
  while (sens==sens0);	        // est toujours le même
  // on a atteint le sens précédent ! t est l'heure de la marée

  // si 30 secondes après, la marée était toujours dans le même sens
  // on prend la minute suivante comme heure de la marée
  if (signe_derivee(t+.5*uneminute)==sens) t=t+uneminute;

  return t;
}
//---------------------------------------------------------------------------
bool signe_derivee(double t)
// la derivée est de la forme -vsin(vt-p)
{ double sens=0.;
  double var;

  for (int i=0; i<nb_const;i++)
    if (ampli[i]>0.)
      { var = reduc360(table2NC[i].speed*t+equilbrm[i]-phase[i]);
        sens -= nodefctr[i]*ampli[i]*table2NC[i].speed*sin_deg(var);
      }
  // return 1 si étale ou marée montante, 0 si marée descendante
  return (sens>=0.);
}
//---------------------------------------------------------------------------
double amplitude(double t)
{ double amplitude=0.;
  double var;
  for (int i=0; i<nb_const;i++)
    if (ampli[i]>0.)
       { var = reduc360(table2NC[i].speed*t+equilbrm[i]-phase[i]);
         amplitude += nodefctr[i]*ampli[i]*cos_deg(var);
       }
  return amplitude;
}
//---------------------------------------------------------------------------
double amplimax()
{ double ampmax=0.;
  for (int i=0;i<nb_const;i++)
    ampmax = ampmax + ampli[i];
  return ampmax;
}
//---------------------------------------------------------------------------



</pre>
</body>
</html>