| 1 | #include <iostream.h>
|
|---|
| 2 | #include <stdlib.h>
|
|---|
| 3 | #include <stdio.h>
|
|---|
| 4 | #include <string.h>
|
|---|
| 5 | #include <math.h>
|
|---|
| 6 |
|
|---|
| 7 | #ifndef M_PI
|
|---|
| 8 | #define M_PI 3.1415926535
|
|---|
| 9 | #endif
|
|---|
| 10 |
|
|---|
| 11 | #ifndef M_2PI
|
|---|
| 12 | #define M_2PI (2*M_PI)
|
|---|
| 13 | #endif
|
|---|
| 14 |
|
|---|
| 15 | #include "plgalcross.h"
|
|---|
| 16 |
|
|---|
| 17 | //////////////////////////////////////////////////////////////////////////
|
|---|
| 18 | // cmv 24/9/1999
|
|---|
| 19 | // Pour calculer les azimuth (E of N) d'intersection d'une direction avec
|
|---|
| 20 | // le plan galactique:
|
|---|
| 21 | // ---- INPUT:
|
|---|
| 22 | // TSid : temps sideral (heures decimales [0,24[)
|
|---|
| 23 | // Lat : latitude du lieu (degres decimaux [-90,90])
|
|---|
| 24 | // zRot : angle de rotation de l'axe du detecteur avec le zenith
|
|---|
| 25 | // (0 = zenith, degres decimaux [0,180]).
|
|---|
| 26 | // alpG,decG : alpha,delta du vecteur perpendiculaire au plan considere
|
|---|
| 27 | // (heures et degres decimaux alpG=[0,24[, decG=[-90,90])
|
|---|
| 28 | // ex: Galactic plane, original IAU definition of galactic
|
|---|
| 29 | // North pole is : (12h 49mn ,+27d 24' , 1950)
|
|---|
| 30 | // soit (12h 51mn 30s,+27d 07' 42", 2000)
|
|---|
| 31 | // ---- OUTPUT:
|
|---|
| 32 | // azCr1,2 : azimuths (E of N) d'intersection du plan definit au dessus
|
|---|
| 33 | // avec le cercle balaye par l'axe du detecteur
|
|---|
| 34 | // (degres decimaux [0,360[).
|
|---|
| 35 | // Angle + sens retrograde (N>E>S>W), origine au Nord.
|
|---|
| 36 |
|
|---|
| 37 | // ---- RETURN VALUE:
|
|---|
| 38 | // -1 : mauvaise valeur en entree, relisez la doc!
|
|---|
| 39 | // 0 : OK, une ou deux intersections
|
|---|
| 40 | // 1 : pas d'intersection
|
|---|
| 41 | // 2 : 1 solution avec azimuth indetermine (sur l'axe +oZ ou -oZ)
|
|---|
| 42 | // 4 : une infinite d'intersections pour tous les azimuths
|
|---|
| 43 | // (dans plan horizontal)
|
|---|
| 44 | //////////////////////////////////////////////////////////////////////////
|
|---|
| 45 |
|
|---|
| 46 | int PlGalCross(double TSid,double Lat,double zRot,double alpG,double decG
|
|---|
| 47 | ,double& azCr1,double& azCr2)
|
|---|
| 48 | {
|
|---|
| 49 | azCr1=0.,azCr2=0.;
|
|---|
| 50 |
|
|---|
| 51 | if(TSid<0. || TSid>=24.) return -1;
|
|---|
| 52 | if(Lat<-90. || Lat>90. ) return -1;
|
|---|
| 53 | if(zRot<0. || zRot>180.) return -1;
|
|---|
| 54 | if(alpG<0. || alpG>=24.) return -1;
|
|---|
| 55 | if(decG<-90. || decG>90. ) return -1;
|
|---|
| 56 |
|
|---|
| 57 | // (ag,hg)=(azimuth,elevation) du vecteur G perpendiculaire au plan
|
|---|
| 58 | double aG, hG;
|
|---|
| 59 | EquatToHoriz(alpG,decG,TSid,Lat,aG,hG);
|
|---|
| 60 | //printf("PlGalCross: tsid=%g lat=%g perp plan: alpha=%g dec=%g\n"
|
|---|
| 61 | // ,TSid,Lat,alpG,decG);
|
|---|
| 62 | //printf(" -> elev=%g azi=%g (+180=%g)\n"
|
|---|
| 63 | // ,hG,aG,((aG>180.)?aG-180.:aG+180.));
|
|---|
| 64 |
|
|---|
| 65 | // on convertit l'angle "zRot" en elevation des coord horizontales
|
|---|
| 66 | double hRot = 90.-zRot;
|
|---|
| 67 | //printf("PlGalCross: axe rot: elev = %g dist.zen = %g\n",hRot,zRot);
|
|---|
| 68 |
|
|---|
| 69 | // on convertit en radian
|
|---|
| 70 | hRot /= DEG_IN_RADIAN; aG /= DEG_IN_RADIAN; hG /= DEG_IN_RADIAN;
|
|---|
| 71 |
|
|---|
| 72 | // Resolution de l'intersection:
|
|---|
| 73 | // Dans le repere horizontal, on cherche le lieu des points M
|
|---|
| 74 | // du plan considere (OM.G=0) ayant l'elevation hRot.
|
|---|
| 75 | // Soient (a,h)=(azimuth,elevation) du point M (vecteur OM)
|
|---|
| 76 | // ATTENTION: l'azimuth aM,aG tourne + dans sens retrograde (N>E>S>W)!
|
|---|
| 77 | // En coordonnees cartesiennes:
|
|---|
| 78 | // | cos(h).cos(a) | | cos(hG).cos(aG) |
|
|---|
| 79 | // OM= | -cos(h).sin(a) | G= | -cos(hG).sin(aG) |
|
|---|
| 80 | // | sin(h) | | sin(hG) |
|
|---|
| 81 | // OM.G=0 avec h=hRot et ||OM||=1
|
|---|
| 82 | // - Equation du plan : a*x + b*y + c*z = 0
|
|---|
| 83 | // - Equation du cercle que decrit l'axe du detecteur sur la sphere celeste:
|
|---|
| 84 | // x^2 + y^2 = R^2 et z=sin(hRot) avec R=cos(hRot)
|
|---|
| 85 | // - Equation de la droite d'intersection du plan avec le plan du cercle:
|
|---|
| 86 | // a*x + b*y + c*z = 0 avec z=sin(hRot)=cos(zRot)
|
|---|
| 87 | double ch=cos(hRot);
|
|---|
| 88 | double z=sin(hRot), chG=cos(hG), shG=sin(hG);
|
|---|
| 89 | if(ch*ch<SMALL_ANGLE_2) {
|
|---|
| 90 | // ch=0 : le detecteur ne tourne pas (zRot==0 ou 180),
|
|---|
| 91 | // il est selon l'axe OZ ou -OZ. Une seule solution possible
|
|---|
| 92 | // si le plan passe par le pole cad si hG==0 et l'azimuth
|
|---|
| 93 | // est indetermine.
|
|---|
| 94 | if(shG*shG<SMALL_ANGLE_2) return 2; else return 1;
|
|---|
| 95 | } else if(chG*chG<SMALL_ANGLE_2) {
|
|---|
| 96 | // chG==0 : le plan est dans le plan horizontal. Seules intersections
|
|---|
| 97 | // possibles si le detecteur tourne dans le plan horizontal cad si hRot=0.
|
|---|
| 98 | // Dans ce cas une infinite d'intersections pour tous les azimuths.
|
|---|
| 99 | if(z*z<SMALL_ANGLE_2) return 3; else return 1;
|
|---|
| 100 | } else {
|
|---|
| 101 | double a=chG*cos(aG), b=-chG*sin(aG), cz=shG*z, R2=ch*ch;
|
|---|
| 102 | //printf("PlGalCross: plan a=%g b=%g c=%g R=%g\n",a,b,shG,ch);
|
|---|
| 103 | //Stabilite numerique: |b|>|a| on resoud pour "x" sinon pour "y"
|
|---|
| 104 | // **** cas 1-/ |b|>|a| : droite y = -(a/b)*x-cz/b
|
|---|
| 105 | // cas 2-/ |b|<=|a| : droite x = -(b/a)*y-cz/a
|
|---|
| 106 | // **** En remplacant dans l'equation du cercle
|
|---|
| 107 | // cas 1-/ Ax^2+2Bx+C=0
|
|---|
| 108 | // A=1+(a/b)^2; B=acz/b^2; C=(cz/b)^2-R^2
|
|---|
| 109 | // cas 2-/ Ay^2+2By+C=0
|
|---|
| 110 | // A=1+(b/a)^2; B=bcz/a^2; C=(cz/a)^2-R^2
|
|---|
| 111 | // **** Solutions: det=B^2-AC et sol=(-B +/- sqrt(det))/A
|
|---|
| 112 | // cas 1-/ x1,x2 ==> y = -(a/b)*x-cz/b ==> y1,y2
|
|---|
| 113 | // cas 2-/ y1,y2 ==> x = -(b/a)*x-cz/a ==> x1,x2
|
|---|
| 114 | // **** Azimuth: az=atan2(-y/cos(hRot),x/cos(hRot))
|
|---|
| 115 | // (car sens retrograde cf eq OM=|...|)
|
|---|
| 116 | bool calcx = (fabs(b)>fabs(a))? true: false;
|
|---|
| 117 | double A,B,C;
|
|---|
| 118 | if(calcx) {A=1+a*a/(b*b); B=a*cz/(b*b); C=cz*cz/(b*b)-R2;}
|
|---|
| 119 | else {A=1+b*b/(a*a); B=b*cz/(a*a); C=cz*cz/(a*a)-R2;}
|
|---|
| 120 | // Resolution
|
|---|
| 121 | double det = B*B-A*C;
|
|---|
| 122 | if(det<0.) return 1;
|
|---|
| 123 | double sol1 = (-B+sqrt(det))/A, sol2 = (-B-sqrt(det))/A;
|
|---|
| 124 | double x1,x2,y1,y2;
|
|---|
| 125 | if(calcx) {x1=sol1; x2=sol2; y1=-a/b*x1-cz/b; y2=-a/b*x2-cz/b;}
|
|---|
| 126 | else {y1=sol1; y2=sol2; x1=-b/a*y1-cz/a; x2=-b/a*y2-cz/a;}
|
|---|
| 127 | //printf("PlGalCross: int1=(%g,%g,%g) int2=(%g,%g,%g)\n",x1,y1,z,x2,y2,z);
|
|---|
| 128 | // Calcul des azimuths
|
|---|
| 129 | azCr1 = atan2(-y1/ch,x1/ch); azCr2 = atan2(-y2/ch,x2/ch);
|
|---|
| 130 | // on convertit en degres decimaux
|
|---|
| 131 | azCr1 *= DEG_IN_RADIAN; azCr2 *= DEG_IN_RADIAN;
|
|---|
| 132 | if(azCr1<0.) azCr1+=360.; if(azCr1>=360.) azCr1-=360.; if(azCr1<0.) azCr1=0.;
|
|---|
| 133 | if(azCr2<0.) azCr2+=360.; if(azCr2>=360.) azCr2-=360.; if(azCr2<0.) azCr2=0.;
|
|---|
| 134 | //printf("PlGalCross: azCr1 = %g, azCr2 = %g (azCr1+azCr2)/2 = %g\n"
|
|---|
| 135 | //,azCr1,azCr2,((azCr1+azCr2>=720.)?(azCr1+azCr2)/2-360.:(azCr1+azCr2)/2));
|
|---|
| 136 | return 0;
|
|---|
| 137 | }
|
|---|
| 138 |
|
|---|
| 139 | }
|
|---|
| 140 |
|
|---|
| 141 | ////////////////////////////////////////////////////////////////////////////////
|
|---|
| 142 | int EquatToHoriz(double alpha,double delta,double tsid,double latitude
|
|---|
| 143 | ,double& azimuth,double& elevation)
|
|---|
| 144 | // Pour convertir les coordonnees equatoriales en coordonnees horizontales
|
|---|
| 145 | // pour une latitude terrestre et un temps sideral donnes.
|
|---|
| 146 | // INPUT:
|
|---|
| 147 | // alpha,delta : coord equat en heures [0,24[ et degres [-90,90] decimaux
|
|---|
| 148 | // tsid : temps sideral en heures [0,24[ decimales
|
|---|
| 149 | // (ou angle horaire du point vernal ou ascension droite du zenith)
|
|---|
| 150 | // Pour une direction quelconque on a: tsid = ha+a
|
|---|
| 151 | // ou "ha" est l'angle horaire et "a" l'ascension droite
|
|---|
| 152 | // latitude : latitude du lieu en degre [-90,90] decimaux
|
|---|
| 153 | // OUTPUT: les coordonnees horizontales:
|
|---|
| 154 | // azimuth : azimuth (E of N) degre [0,360[ decimaux
|
|---|
| 155 | // tourne + sens retrograde (N>E>S>W), origine au Nord,
|
|---|
| 156 | // azimuth "des marins" (serveur BDL, USNO, xephem , skycalc...)
|
|---|
| 157 | // elevation : elevation degre [-90,90] decimaux
|
|---|
| 158 | // (distance zenithale z = 90-elevation)
|
|---|
| 159 | // RETURN: 0 si Ok, 1 si probleme
|
|---|
| 160 | // TRANSFO: (dans cette formule l'origine des azimuths est au Sud,
|
|---|
| 161 | // azimut "des astronomes" ou (E of S))
|
|---|
| 162 | // sin(azi) cos(elev) = sin(ha) cos(dec)
|
|---|
| 163 | // cos(azi) cos(elev) = cos(ha) cos(dec) sin(lat) - sin(dec) cos(lat)
|
|---|
| 164 | // sin(elev) = cos(ha) cos(dec) cos(lat) + sin(dec) sin(lat)
|
|---|
| 165 | // avec azi=azimuth (E of S), elev=elevation, lat=latitude,
|
|---|
| 166 | // dec=declinaison, ha=angle horaire=tsid-alpha
|
|---|
| 167 | // Puis azi(E of S) = azi(E of N) + Pi
|
|---|
| 168 | {
|
|---|
| 169 | azimuth = elevation = 0.;
|
|---|
| 170 | if(alpha<0. || alpha>=24 ) return 1;
|
|---|
| 171 | if(delta<-90. || delta>90. ) return 1;
|
|---|
| 172 | if(tsid<0. || tsid>=24. ) return 1;
|
|---|
| 173 | if(latitude<-90. || latitude>90.) return 1;
|
|---|
| 174 | // conversion en radians
|
|---|
| 175 | alpha /= HRS_IN_RADIAN; delta /= DEG_IN_RADIAN;
|
|---|
| 176 | tsid /= HRS_IN_RADIAN; latitude /= DEG_IN_RADIAN;
|
|---|
| 177 |
|
|---|
| 178 | double cl = cos(latitude), sl = sin(latitude);
|
|---|
| 179 | double ch = cos(tsid-alpha), sh = sin(tsid-alpha);
|
|---|
| 180 | double cd = cos(delta), sd = sin(delta);
|
|---|
| 181 |
|
|---|
| 182 | // calcul de l'elevation [-90,90]
|
|---|
| 183 | elevation = ch*cd*cl+sl*sd;
|
|---|
| 184 | if(elevation<-1.) elevation=-1.; if(elevation>1.) elevation=1.;
|
|---|
| 185 | elevation = asin(elevation); //entre [-Pi/2,Pi/2]
|
|---|
| 186 |
|
|---|
| 187 | // calcul de l'azimuth
|
|---|
| 188 | double ce = cos(elevation);
|
|---|
| 189 | if(fabs(ce)<SMALL_ANGLE) { //elevation=-Pi/2 ou Pi/2
|
|---|
| 190 | azimuth = 0.; //azimuth indifferent
|
|---|
| 191 | } else {
|
|---|
| 192 | azimuth = atan2(sh*cd/ce,(ch*cd*sl-sd*cl)/ce); //entre -Pi et +Pi
|
|---|
| 193 | azimuth += M_PI; //azimuth (E of S) -> (E of N)
|
|---|
| 194 | }
|
|---|
| 195 |
|
|---|
| 196 | azimuth *= DEG_IN_RADIAN; elevation *= DEG_IN_RADIAN;
|
|---|
| 197 | // Ca ne doit pas arriver sauf pb de precision machine
|
|---|
| 198 | if(azimuth<0. || azimuth>=360.) azimuth = 0.;
|
|---|
| 199 | if(elevation<-90.) elevation = -90.; if(elevation>90.) elevation = 90.;
|
|---|
| 200 | return 0;
|
|---|
| 201 | }
|
|---|
| 202 |
|
|---|
| 203 | ////////////////////////////////////////////////////////////////////////////////
|
|---|
| 204 | int FindPerpEquat(double a1,double d1,double a2,double d2,double& ap,double& dp)
|
|---|
| 205 | // Pour trouver les coordonnees equatoriales du vecteur perpendiculaire
|
|---|
| 206 | // a deux autres vecteurs.
|
|---|
| 207 | // INPUT:
|
|---|
| 208 | // a1,d1 : alpha,delta de la premiere direction en
|
|---|
| 209 | // (heures decimales et degres decimaux)
|
|---|
| 210 | // a2,d2 : alpha,delta de la deuxieme direction en
|
|---|
| 211 | // (heures decimales et degres decimaux)
|
|---|
| 212 | // OUTPUT:
|
|---|
| 213 | // ap,dp : alpha,delta de la direction perpendiculaire
|
|---|
| 214 | // telle que 1,2,p fasse un triedre direct
|
|---|
| 215 | // (heures decimales et degres decimaux)
|
|---|
| 216 | // RETURN:
|
|---|
| 217 | // 0 = OK, 1 = probleme
|
|---|
| 218 | {
|
|---|
| 219 | ap = dp = 0.;
|
|---|
| 220 | a1 /= HRS_IN_RADIAN; d1 /= DEG_IN_RADIAN;
|
|---|
| 221 | a2 /= HRS_IN_RADIAN; d2 /= DEG_IN_RADIAN;
|
|---|
| 222 | // coordonnees cartesiennes
|
|---|
| 223 | double v1[3],v2[3],vp[3];
|
|---|
| 224 | v1[0]=cos(d1)*cos(a1); v1[1]=cos(d1)*sin(a1); v1[2]=sin(d1);
|
|---|
| 225 | v2[0]=cos(d2)*cos(a2); v2[1]=cos(d2)*sin(a2); v2[2]=sin(d2);
|
|---|
| 226 | vp[0]=v1[1]*v2[2]-v1[2]*v2[1]; vp[1]=v1[2]*v2[0]-v1[0]*v2[2];
|
|---|
| 227 | vp[2]=v1[0]*v2[1]-v1[1]*v2[0];
|
|---|
| 228 | double nvp = vp[0]*vp[0]+vp[1]*vp[1]+vp[2]*vp[2];
|
|---|
| 229 | if(nvp<=0.) return 1;
|
|---|
| 230 | nvp = sqrt(nvp);
|
|---|
| 231 | vp[0] /= nvp; vp[1] /= nvp; vp[2] /= nvp;
|
|---|
| 232 | // | vp[0] | | cos(dp)*cos(ap) |
|
|---|
| 233 | // | vp[1] | = | cos(dp)*sin(ap) |
|
|---|
| 234 | // | vp[2] | | sin(dp) |
|
|---|
| 235 | dp = asin(vp[2]); // asin retourne entre -Pi/2 et Pi/2 ... Ok
|
|---|
| 236 | double cdp = cos(dp);
|
|---|
| 237 | if(fabs(cdp)<SMALL_ANGLE) { // cas ou dp=+Pi/2 ou -Pi/2, ap est indifferent
|
|---|
| 238 | ap = 0.;
|
|---|
| 239 | } else {
|
|---|
| 240 | ap = atan2(vp[1]/cdp,vp[0]/cdp); // renvoie entre -Pi et +Pi
|
|---|
| 241 | if(ap<0.) ap += M_2PI;
|
|---|
| 242 | }
|
|---|
| 243 |
|
|---|
| 244 | ap *= HRS_IN_RADIAN; dp *= DEG_IN_RADIAN;
|
|---|
| 245 | // Ca ne doit pas arriver sauf pb de precision machine
|
|---|
| 246 | if(ap<0. || ap>=24.) ap = 0.;
|
|---|
| 247 | if(dp<-90.) dp = -90.; if(dp>90.) dp = 90.;
|
|---|
| 248 | return 0;
|
|---|
| 249 | }
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
