Index: /trunk/Poubelle/archTOI.old/plgalcross.cc
===================================================================
--- /trunk/Poubelle/archTOI.old/plgalcross.cc	(revision 423)
+++ /trunk/Poubelle/archTOI.old/plgalcross.cc	(revision 423)
@@ -0,0 +1,232 @@
+#include <iostream.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+
+#ifndef M_PI 
+#define M_PI 3.14159265358979
+#endif
+
+#ifndef M_2PI
+#define M_2PI (2*M_PI)
+#endif
+
+#include "plgalcross.h"
+
+//////////////////////////////////////////////////////////////////////////
+// Pour calculer les azimuth 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 (galactique) 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 d'intersection du plan definit au dessus
+//           avec le cercle balaye par l'axe du detecteur
+//           (degres decimaux [0,360[)
+// ---- 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);
+
+// on convertit l'angle "zRot" en elevation des coord horizontales
+double hRot = 90.-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 le sens retrograde!
+// 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)
+double R2=cos(hRot); R2*=R2;
+double z=sin(hRot), chG=cos(hG), shG=sin(hG);
+if(R2<SMALL_ANGLE_2) {
+  // R2=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;
+  //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,x) (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;}
+  // Calcul des azimuths
+  azCr1 = atan2(-y1,x1); azCr2 = atan2(-y2,x2);
+  // 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.;
+  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 equatoriales 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 degre [0,360[ decimaux
+//   elevation : elevation degre [-90,90] decimaux
+//               (distance zenithale z = 90-elevation)
+// RETURN: 0 si Ok, 1 si probleme
+// TRANSFO:
+// 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, elev=elevation, lat=latitude,
+//      dec=declinaison, ha=angle horaire=tsid-alpha
+{
+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,Pi]
+
+// calcul de l'azimuth
+if(fabs(cos(elevation))<SMALL_ANGLE) { //elevation=-Pi/2 ou Pi/2
+  azimuth = 0.; //azimuth indifferent
+} else {
+  azimuth = atan2(sh*cd,ch*cd*sl-sd*cl); //entre -Pi et +Pi
+  if(azimuth<0.) azimuth += M_2PI; //remise azimuth entre 0 et 2Pi
+}
+
+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
+nvp = cos(dp);
+if(fabs(nvp)<1.e-30) { // cas ou dp=+Pi/2 ou -Pi/2, ap est indifferent
+  ap = 0.;
+} else {
+  ap = atan2(vp[1],vp[0]); // 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;
+}
Index: /trunk/Poubelle/archTOI.old/plgalcross.h
===================================================================
--- /trunk/Poubelle/archTOI.old/plgalcross.h	(revision 423)
+++ /trunk/Poubelle/archTOI.old/plgalcross.h	(revision 423)
@@ -0,0 +1,11 @@
+#define DEG_IN_RADIAN 57.2957795130823
+#define HRS_IN_RADIAN 3.819718634205
+// SMALL_ANGLE 1.e-8 radian = 1e-3 ArcSec
+#define SMALL_ANGLE 1.e-8
+#define SMALL_ANGLE_2 (SMALL_ANGLE*SMALL_ANGLE)
+
+int PlGalCross(double TSid,double Lat,double zRot,double alpG,double decG
+              ,double& azCr1,double& azCr2);
+int EquatToHoriz(double alpha,double delta,double TSid,double latitude
+                ,double& azimuth,double& hauteur);
+int FindPerpEquat(double a1,double d1,double a2,double d2,double& ap,double& dp);
Index: /trunk/Poubelle/archTOI.old/templocator.cc
===================================================================
--- /trunk/Poubelle/archTOI.old/templocator.cc	(revision 422)
+++ /trunk/Poubelle/archTOI.old/templocator.cc	(revision 423)
@@ -4,4 +4,7 @@
 #include "aa_hadec.h"
 }
+
+#include "fitsio.h"
+#include "plgalcross.h"
 
 #ifndef M_PI
@@ -14,29 +17,78 @@
 {
   lon = lat = ts = 0;
+  ra  = dec = -99999;
+  xSampleNum = -99999;
+  fitsfile* fptr;
+  int status=0;
+  fits_open_file(&fptr, "samplenum_gal_cross.fits", READONLY, &status);
+  int simple, bitpix, naxis;
+  long naxes;
+  long pcount, gcount;
+  int extend;
+  fits_read_imghdr(fptr, 1, &simple, &bitpix, &naxis, &naxes, &pcount, &gcount, &extend, &status);
+  nGalCross = naxes;
+  crossings = new long[nGalCross];
+  int anynul;
+  fits_read_img_lng(fptr, 0, 1, nGalCross, 0, crossings, &anynul, &status);
+  fits_close_file(fptr, &status);
+  fits_report_error(stderr, status);  /* print out any error messages */
+}
+           
+           
+void TempLocator::setEarthPos(double lon, double lat) {
+  if (this->lon == lon && this->lat == lat) return;
+  this->lon = lon;
+  this->lat = lat;
+  ra  = dec = -99999; xSampleNum = -99999;
 }
 
-void TempLocator::setEarthPos(double lon, double lat) {
-  this->lon = lon;
-  this->lat = lat;
+void TempLocator::setTSid(double ts) {
+  if (this->ts == ts) return;
+  this->ts  = ts;
+  ra  = dec = -99999; xSampleNum = -99999;
 }
-void TempLocator::setTSid(double ts) {
-  this->ts  = ts;
+  
+void   TempLocator::ComputeZenith() {
+  double ha;
+  aa_hadec (lat * M_PI/180, .5 * M_PI, 0, &ha, &dec);
+  ra = - (ha * 180. / M_PI / 15) + (ts/3600.);
+  dec = dec * 180. / M_PI;
 }
   
 double TempLocator::getAlphaZenith() {
-  double ha, dec;
-  aa_hadec (lat * M_PI/180, .5 * M_PI, 0, &ha, &dec);
-  double ra = - (ha * 180. / M_PI / 15) + (ts/3600.);
-  //dec = dec * 180. / M_PI;
+  if (ra < -100) ComputeZenith();
   return ra;
 }
 
 double TempLocator::getDeltaZenith() {
-  double ha, dec;
-  aa_hadec (lat * M_PI/180, .5 * M_PI, 0, &ha, &dec);
-  //double ra = - (ha * 180. / M_PI / 15) + (ts/3600.);
-  dec = dec * 180. / M_PI;
+  if (dec < -100) ComputeZenith();
   return dec;
 }
 
+void TempLocator::findGeomFromGC(int sampleNum) // pour le bolo qui voit les xing
+{
+  if (sampleNum == xSampleNum) return;
+  if (dec < -100) ComputeZenith();
+  
+  azimBolGC = -9999;
+  
+  // On trouve les croisements juste avant et juste apres notre sampleNum
+  int icross;
+  for (icross=0; icross<nGalCross; icross++) {
+    if (crossings[icross] > sampleNum) break;
+  }
+  if (icross == 0 || icross >= nGalCross) return;
+  
+  
+}
 
+
+ 			//	 ;
+  			//	 ;                            /\         4     2
+  			//	 ;                  elevation ||                       positive scanning, clockwise
+  			//	 ;                            ||       6    1    5      ---------->
+  			//	 ;                            ||                      == positive azimut
+  			//	 ;                            ||         x     3
+            //                                         |----|      = 0.78 deg / cos(elev)
+            // bol 1 = 41° elevation
+            
Index: /trunk/Poubelle/archTOI.old/templocator.h
===================================================================
--- /trunk/Poubelle/archTOI.old/templocator.h	(revision 422)
+++ /trunk/Poubelle/archTOI.old/templocator.h	(revision 423)
@@ -14,5 +14,9 @@
   double getDeltaZenith();
   
-
+      
+  double getAzimutBolo(int sampleNum, int ibolo); // ibolo, comme dans les TOI
+  double getElvBolo(int sampleNum, int ibolo);
+  double getAlphaBolo(int sampleNum, int ibolo);
+  double getDeltaBolo(int sampleNum, int ibolo);
 
 protected:
@@ -20,4 +24,16 @@
   double lat; // degrees  +=NORD
   double ts;  // secondes 
+  
+  int    nGalCross;
+  long*   crossings; // samplenum des croisements galactiques
+  
+  double ra;
+  double dec;
+  int xSampleNum; // pour le dernier calcul
+  double azimBolGC;
+  
+  void findGeomFromGC(int sampleNum); // pour le bolo qui voit les xing
+    // il faut avoir positionne lon, lat et ts.
+  void ComputeZenith();
 };
 
