#include <iostream.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#ifndef M_PI
#define M_PI 3.1415926535
#endif

#ifndef M_2PI
#define M_2PI (2*M_PI)
#endif

#include "plgalcross.h"

//////////////////////////////////////////////////////////////////////////
//                                cmv 24/9/1999
// Pour calculer les azimuth (E of N) d'intersection d'une direction avec
// le plan galactique:
// ---- INPUT:
// TSid : temps sideral (heures decimales [0,24[)
// Lat : latitude du lieu (degres decimaux [-90,90])
// zRot : angle de rotation de l'axe du detecteur avec le zenith
//        (0 = zenith, degres decimaux [0,180]).
// alpG,decG : alpha,delta du vecteur perpendiculaire au plan considere
//             (heures et degres decimaux alpG=[0,24[, decG=[-90,90])
//     ex: Galactic plane, original IAU definition of galactic
//         North pole is : (12h 49mn    ,+27d 24'    , 1950)
//                   soit  (12h 51mn 30s,+27d 07' 42", 2000)
// ---- OUTPUT:
// azCr1,2 : azimuths (E of N) d'intersection du plan definit au dessus
//           avec le cercle balaye par l'axe du detecteur
//           (degres decimaux [0,360[).
//           Angle + sens retrograde (N>E>S>W), origine au Nord.

// ---- RETURN VALUE:
//  -1 : mauvaise valeur en entree, relisez la doc!
//   0 : OK, une ou deux intersections
//   1 : pas d'intersection
//   2 : 1 solution avec azimuth indetermine (sur l'axe +oZ ou -oZ)
//   4 : une infinite d'intersections pour tous les azimuths
//       (dans plan horizontal)
//////////////////////////////////////////////////////////////////////////

int PlGalCross(double TSid,double Lat,double zRot,double alpG,double decG
              ,double& azCr1,double& azCr2)
{
azCr1=0.,azCr2=0.;

if(TSid<0.   || TSid>=24.) return -1;
if(Lat<-90.  || Lat>90.  ) return -1;
if(zRot<0.   || zRot>180.) return -1;
if(alpG<0.   || alpG>=24.) return -1;
if(decG<-90. || decG>90. ) return -1;

// (ag,hg)=(azimuth,elevation) du vecteur G perpendiculaire au plan
double aG, hG;
EquatToHoriz(alpG,decG,TSid,Lat,aG,hG);
//printf("PlGalCross: tsid=%g lat=%g perp plan: alpha=%g dec=%g\n"
//      ,TSid,Lat,alpG,decG);
//printf("            -> elev=%g azi=%g (+180=%g)\n"
//      ,hG,aG,((aG>180.)?aG-180.:aG+180.));

// on convertit l'angle "zRot" en elevation des coord horizontales
double hRot = 90.-zRot;
//printf("PlGalCross: axe rot: elev = %g  dist.zen = %g\n",hRot,zRot);

// on convertit en radian
hRot /= DEG_IN_RADIAN; aG /= DEG_IN_RADIAN; hG /= DEG_IN_RADIAN;

// Resolution de l'intersection:
// Dans le repere horizontal, on cherche le lieu des points M
// du plan considere (OM.G=0) ayant l'elevation hRot.
// Soient (a,h)=(azimuth,elevation) du point M (vecteur OM)
// ATTENTION: l'azimuth aM,aG tourne + dans sens retrograde (N>E>S>W)!
// En coordonnees cartesiennes:
//     |  cos(h).cos(a) |      |  cos(hG).cos(aG) |
// OM= | -cos(h).sin(a) |   G= | -cos(hG).sin(aG) |
//     |  sin(h)        |      |  sin(hG)         |
// OM.G=0 avec h=hRot et ||OM||=1
// - Equation du plan : a*x + b*y + c*z = 0
// - Equation du cercle que decrit l'axe du detecteur sur la sphere celeste:
//            x^2 + y^2 = R^2  et z=sin(hRot) avec R=cos(hRot)
// - Equation de la droite d'intersection du plan avec le plan du cercle:
//            a*x + b*y + c*z = 0 avec z=sin(hRot)=cos(zRot)
double ch=cos(hRot);
double z=sin(hRot), chG=cos(hG), shG=sin(hG);
if(ch*ch<SMALL_ANGLE_2) {
  // ch=0 : le detecteur ne tourne pas (zRot==0 ou 180),
  // il est selon l'axe OZ ou -OZ. Une seule solution possible
  // si le plan passe par le pole cad si hG==0 et l'azimuth
  // est indetermine.
  if(shG*shG<SMALL_ANGLE_2) return 2; else return 1;
} else if(chG*chG<SMALL_ANGLE_2) {
  // chG==0 : le plan est dans le plan horizontal. Seules intersections
  // possibles si le detecteur tourne dans le plan horizontal cad si hRot=0.
  // Dans ce cas une infinite d'intersections pour tous les azimuths.
  if(z*z<SMALL_ANGLE_2) return 3; else return 1;
} else {
  double a=chG*cos(aG), b=-chG*sin(aG), cz=shG*z, R2=ch*ch;
  //printf("PlGalCross: plan a=%g b=%g c=%g R=%g\n",a,b,shG,ch);
  //Stabilite numerique: |b|>|a| on resoud pour "x" sinon pour "y"
  // **** cas 1-/ |b|>|a|  :  droite y = -(a/b)*x-cz/b
  //      cas 2-/ |b|<=|a| :  droite x = -(b/a)*y-cz/a
  // **** En remplacant dans l'equation du cercle
  //      cas 1-/ Ax^2+2Bx+C=0
  //              A=1+(a/b)^2; B=acz/b^2; C=(cz/b)^2-R^2
  //      cas 2-/ Ay^2+2By+C=0
  //              A=1+(b/a)^2; B=bcz/a^2; C=(cz/a)^2-R^2
  // **** Solutions: det=B^2-AC et sol=(-B +/- sqrt(det))/A
  //      cas 1-/ x1,x2 ==> y = -(a/b)*x-cz/b ==> y1,y2
  //      cas 2-/ y1,y2 ==> x = -(b/a)*x-cz/a ==> x1,x2
  // **** Azimuth: az=atan2(-y/cos(hRot),x/cos(hRot))
  //               (car sens retrograde cf eq OM=|...|)
  bool calcx = (fabs(b)>fabs(a))? true: false;
  double A,B,C;
  if(calcx) {A=1+a*a/(b*b); B=a*cz/(b*b); C=cz*cz/(b*b)-R2;}
    else    {A=1+b*b/(a*a); B=b*cz/(a*a); C=cz*cz/(a*a)-R2;}
  // Resolution
  double det = B*B-A*C;
  if(det<0.) return 1;
  double sol1 = (-B+sqrt(det))/A, sol2 = (-B-sqrt(det))/A;
  double x1,x2,y1,y2;
  if(calcx) {x1=sol1; x2=sol2; y1=-a/b*x1-cz/b; y2=-a/b*x2-cz/b;}
     else   {y1=sol1; y2=sol2; x1=-b/a*y1-cz/a; x2=-b/a*y2-cz/a;}
  //printf("PlGalCross: int1=(%g,%g,%g) int2=(%g,%g,%g)\n",x1,y1,z,x2,y2,z);
  // Calcul des azimuths
  azCr1 = atan2(-y1/ch,x1/ch); azCr2 = atan2(-y2/ch,x2/ch);
  // on convertit en degres decimaux
  azCr1 *= DEG_IN_RADIAN; azCr2 *= DEG_IN_RADIAN;
  if(azCr1<0.) azCr1+=360.; if(azCr1>=360.) azCr1-=360.; if(azCr1<0.) azCr1=0.;
  if(azCr2<0.) azCr2+=360.; if(azCr2>=360.) azCr2-=360.; if(azCr2<0.) azCr2=0.;
  //printf("PlGalCross: azCr1 = %g,    azCr2 = %g  (azCr1+azCr2)/2 = %g\n"
  //,azCr1,azCr2,((azCr1+azCr2>=720.)?(azCr1+azCr2)/2-360.:(azCr1+azCr2)/2));
  return 0;
}

}

////////////////////////////////////////////////////////////////////////////////
int EquatToHoriz(double alpha,double delta,double tsid,double latitude
                 ,double& azimuth,double& elevation)
// Pour convertir les coordonnees equatoriales en coordonnees horizontales
// pour une latitude terrestre et un temps sideral donnes.
// INPUT:
//   alpha,delta : coord equat en heures [0,24[ et degres [-90,90] decimaux
//   tsid : temps sideral en heures [0,24[ decimales
//          (ou angle horaire du point vernal ou ascension droite du zenith)
//          Pour une direction quelconque on a: tsid = ha+a 
//          ou "ha" est l'angle horaire et "a" l'ascension droite
//   latitude : latitude du lieu en degre [-90,90] decimaux
// OUTPUT: les coordonnees horizontales:
//   azimuth : azimuth (E of N) degre [0,360[ decimaux
//             tourne +  sens retrograde (N>E>S>W), origine au Nord,
//             azimuth "des marins" (serveur BDL, USNO, xephem , skycalc...)
//   elevation : elevation degre [-90,90] decimaux
//               (distance zenithale z = 90-elevation)
// RETURN: 0 si Ok, 1 si probleme
// TRANSFO: (dans cette formule l'origine des azimuths est au Sud,
//           azimut "des astronomes" ou  (E of S))
// sin(azi) cos(elev) = sin(ha) cos(dec)
// cos(azi) cos(elev) = cos(ha) cos(dec) sin(lat) - sin(dec) cos(lat)
//          sin(elev) = cos(ha) cos(dec) cos(lat) + sin(dec) sin(lat)
// avec azi=azimuth (E of S), elev=elevation, lat=latitude,
//      dec=declinaison, ha=angle horaire=tsid-alpha
// Puis azi(E of S) = azi(E of N) + Pi
{
azimuth = elevation = 0.;
if(alpha<0.      || alpha>=24   )   return 1;
if(delta<-90.    || delta>90.   ) return 1;
if(tsid<0.       || tsid>=24.   ) return 1;
if(latitude<-90. || latitude>90.) return 1;
// conversion en radians
alpha /= HRS_IN_RADIAN; delta    /= DEG_IN_RADIAN;
tsid  /= HRS_IN_RADIAN; latitude /= DEG_IN_RADIAN;

double cl = cos(latitude),   sl = sin(latitude);
double ch = cos(tsid-alpha), sh = sin(tsid-alpha);
double cd = cos(delta),      sd = sin(delta);

// calcul de l'elevation [-90,90]
elevation = ch*cd*cl+sl*sd;
if(elevation<-1.) elevation=-1.; if(elevation>1.) elevation=1.;
elevation = asin(elevation); //entre [-Pi/2,Pi/2]

// calcul de l'azimuth
double ce = cos(elevation);
if(fabs(ce)<SMALL_ANGLE) { //elevation=-Pi/2 ou Pi/2
  azimuth = 0.; //azimuth indifferent
} else {
  azimuth = atan2(sh*cd/ce,(ch*cd*sl-sd*cl)/ce); //entre -Pi et +Pi
  azimuth += M_PI; //azimuth (E of S) -> (E of N)
}

azimuth *= DEG_IN_RADIAN; elevation *= DEG_IN_RADIAN;
// Ca ne doit pas arriver sauf pb de precision machine
if(azimuth<0. || azimuth>=360.) azimuth = 0.;
if(elevation<-90.) elevation = -90.; if(elevation>90.) elevation = 90.;
return 0;
}

////////////////////////////////////////////////////////////////////////////////
int FindPerpEquat(double a1,double d1,double a2,double d2,double& ap,double& dp)
// Pour trouver les coordonnees equatoriales du vecteur perpendiculaire
// a deux autres vecteurs.
// INPUT:
//  a1,d1 : alpha,delta de la premiere direction en
//          (heures decimales et degres decimaux)
//  a2,d2 : alpha,delta de la deuxieme direction en
//          (heures decimales et degres decimaux)
// OUTPUT:
//  ap,dp : alpha,delta de la direction perpendiculaire
//          telle que 1,2,p fasse un triedre direct
//          (heures decimales et degres decimaux)
// RETURN:
//  0 = OK, 1 = probleme
{
ap = dp = 0.;
a1 /= HRS_IN_RADIAN; d1 /= DEG_IN_RADIAN;
a2 /= HRS_IN_RADIAN; d2 /= DEG_IN_RADIAN;
// coordonnees cartesiennes
double v1[3],v2[3],vp[3];
v1[0]=cos(d1)*cos(a1); v1[1]=cos(d1)*sin(a1); v1[2]=sin(d1);
v2[0]=cos(d2)*cos(a2); v2[1]=cos(d2)*sin(a2); v2[2]=sin(d2);
vp[0]=v1[1]*v2[2]-v1[2]*v2[1]; vp[1]=v1[2]*v2[0]-v1[0]*v2[2];
                               vp[2]=v1[0]*v2[1]-v1[1]*v2[0];
double nvp = vp[0]*vp[0]+vp[1]*vp[1]+vp[2]*vp[2];
if(nvp<=0.) return 1;
nvp = sqrt(nvp);
vp[0] /= nvp; vp[1] /= nvp; vp[2] /= nvp;
// | vp[0] |   | cos(dp)*cos(ap) |
// | vp[1] | = | cos(dp)*sin(ap) |
// | vp[2] |   | sin(dp)         |
dp = asin(vp[2]); // asin retourne entre -Pi/2 et Pi/2 ... Ok
double cdp = cos(dp);
if(fabs(cdp)<SMALL_ANGLE) { // cas ou dp=+Pi/2 ou -Pi/2, ap est indifferent
  ap = 0.;
} else {
  ap = atan2(vp[1]/cdp,vp[0]/cdp); // renvoie entre -Pi et +Pi
  if(ap<0.) ap += M_2PI;
}

ap *= HRS_IN_RADIAN; dp *= DEG_IN_RADIAN;
// Ca ne doit pas arriver sauf pb de precision machine
if(ap<0. || ap>=24.) ap = 0.;
if(dp<-90.) dp = -90.; if(dp>90.) dp = 90.;
return 0;
}