Changeset 427 in Sophya

Timestamp:

Sep 24, 1999, 9:16:54 AM (26 years ago)

Author:

ansari

Message:

correction cmv azimut galaxie

Location:

trunk/Poubelle/archTOI.old

Files:

• plgalcross.cc (modified) (9 diffs)
• plgalcross.h (modified) (1 diff)

trunk/Poubelle/archTOI.old/plgalcross.cc

@@ -5 +5 @@
 #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
+// Pour calculer les azimuth (E of N) d'intersection d'une direction avec
 // le plan galactique:
 // ---- INPUT:
@@ -30 +22 @@
 //                   soit  (12h 51mn 30s, +27d 07' 42", 2000)
 // ---- OUTPUT:
-// azCr1,2 : azimuths d'intersection du plan definit au dessus
+// 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[)
@@ -78 +70 @@
 // - 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 ch=cos(hRot);
 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),
+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
@@ -92 +84 @@
   if(z*z<SMALL_ANGLE_2) return 3; else return 1;
 } else {
-  double a=chG*cos(aG), b=-chG*sin(aG), cz=shG*z;
+  double a=chG*cos(aG), b=-chG*sin(aG), cz=shG*z, R2=ch*ch;
   //Stabilite numerique: |b|>|a| on resoud pour "x" sinon pour "y"
   // **** cas 1-/ |b|>|a|  :  droite y = -(a/b)*x-cz/b
@@ -104 +96 @@
   //      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=|...|)
+  // **** 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;
@@ -117 +110 @@
      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);
+  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;
@@ -140 +133 @@
 //   latitude : latitude du lieu en degre [-90,90] decimaux
 // OUTPUT: les coordonnees horizontales:
-//   azimuth : azimuth degre [0,360[ decimaux
+//   azimuth : azimuth (E of N) degre [0,360[ decimaux
+//             ici compte + dans sens retrograde avec 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:
+// 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, elev=elevation, lat=latitude,
+// 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.;
@@ -167 +164 @@
 elevation = ch*cd*cl+sl*sd;
 if(elevation<-1.) elevation=-1.; if(elevation>1.) elevation=1.;
-elevation = asin(elevation); //entre [-Pi,Pi]
+elevation = asin(elevation); //entre [-Pi/2,Pi/2]
 
 // calcul de l'azimuth
-if(fabs(cos(elevation))<SMALL_ANGLE) { //elevation=-Pi/2 ou Pi/2
+double ce = cos(elevation);
+if(fabs(ce)<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 = atan2(sh*cd/ce,(ch*cd*sl-sd*cl)/ce); //entre -Pi et +Pi
+  azimuth += M_PI; //azimuth (E of S) -> (E of N)
 }
 
@@ -217 +215 @@
 // | 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
+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],vp[0]); // renvoie entre -Pi et +Pi
+  ap = atan2(vp[1]/cdp,vp[0]/cdp); // renvoie entre -Pi et +Pi
   if(ap<0.) ap += M_2PI;
 }

trunk/Poubelle/archTOI.old/plgalcross.h

@@ -8 +8 @@
               ,double& azCr1,double& azCr2);
 int EquatToHoriz(double alpha,double delta,double TSid,double latitude
-                ,double& azimuth,double& hauteur);
+                ,double& azimuth,double& elevation);
 int FindPerpEquat(double a1,double d1,double a2,double d2,double& ap,double& dp);