source: Sophya/trunk/Poubelle/archTOI.old/templocator.cc@ 635

Last change on this file since 635 was 555, checked in by ansari, 26 years ago

portage cxx en cours
File size: 6.9 KB
Line 
1#include "templocator.h"
2#include "gondolageom.h"
3#include <math.h>
4extern "C" {
5#include "aa_hadec.h"
6}
7
8#include "fitsio.h"
9#include "plgalcross.h"
10#include "archparam.h"
11
12#ifndef M_PI
13#define M_PI 3.14159265358979323846
14#endif
15
16TempLocator tempLocator;
17
18TempLocator::TempLocator()
19{
20 lon = lat = ts = 0;
21 raZ = decZ = -99999;
22 xSampleNum = -99999;
23 fitsfile* fptr;
24 int status=0;
25 fits_open_file(&fptr, "samplenum_gal_cross.fits", READONLY, &status);
26 int simple, bitpix, naxis;
27 long naxes;
28 long pcount, gcount;
29 int extend;
30 fits_read_imghdr(fptr, 1, &simple, &bitpix, &naxis, &naxes, &pcount, &gcount, &extend, &status);
31 nGalCross = naxes;
32 crossings = new long[nGalCross];
33 int anynul;
34 fits_read_img_lng(fptr, 0, 1, nGalCross, 0, crossings, &anynul, &status);
35 fits_close_file(fptr, &status);
36 fits_report_error(stderr, status); /* print out any error messages */
37}
38
39
40void TempLocator::setEarthPos(double lon, double lat) {
41 if (this->lon == lon && this->lat == lat) return;
42 this->lon = lon;
43 this->lat = lat;
44 raZ = decZ = -99999; xSampleNum = -99999;
45}
46
47void TempLocator::setTSid(double ts) {
48 if (this->ts == ts) return;
49 this->ts = ts;
50 raZ = decZ = -99999; xSampleNum = -99999;
51}
52
53void TempLocator::ComputeZenith() {
54 double ha;
55 aa_hadec (lat * M_PI/180, .5 * M_PI, 0, &ha, &decZ);
56 raZ = - (ha * 180. / M_PI / 15) + (ts/3600.);
57 if (raZ>24) raZ -= 24;
58 if (raZ<0) raZ += 24;
59 decZ = decZ * 180. / M_PI;
60}
61
62double TempLocator::getAlphaZenith() {
63 if (raZ < -100) ComputeZenith();
64 return raZ;
65}
66
67double TempLocator::getDeltaZenith() {
68 if (decZ < -100) ComputeZenith();
69 return decZ;
70}
71
72#define altbolo1 41.5
73
74void TempLocator::findGeomFromGC(int sampleNum) // pour le bolo qui voit les xing
75{
76 if (sampleNum == xSampleNum) return;
77 if (decZ < -100) ComputeZenith();
78
79 azimBolGC = -9999;
80
81 // On trouve les croisements juste avant et juste apres notre sampleNum
82 for (icross=0; icross<nGalCross; icross++) {
83 if (crossings[icross] > sampleNum) break;
84 }
85 if (icross == 0 || icross >= nGalCross) return;
86
87 // On trouve l'azimut du croisement principal pour notre position actuelle
88 double alpG = 12. + 51./60. + 30./3600.;
89 double delG = 27. + 07./60. + 42./3600.;
90 double azCr1, azCr2;
91 int rc = PlGalCross(ts/3600., lat, (90. - altbolo1), alpG, delG, azCr1, azCr2);
92 if (rc != 0) return; // pas deux points d'intersection
93
94 // Il faut determiner le croisement principal, ie le plus proche
95 // du centre galactique. Pendant le vol de Trapani, c'etait celui
96 // le plus proche de 220° d'azimut.
97
98 double azCross = azCr1;
99 if (fabs(azCr2-220) < fabs(azCr1-220)) azCross = azCr2;
100
101 rotSpeed = 360./(crossings[icross] - crossings[icross-1]); // °/sample
102
103 azimBolGC = azCross - (sampleNum - (crossings[icross-1]+12))*rotSpeed;
104 // azimut bolo 1 from crossing, for sampleNum
105 if (azimBolGC > 360) azimBolGC -= 360;
106}
107
108double TempLocator::getRotSpeed(int sampleNum) {
109 findGeomFromGC(sampleNum);
110 return rotSpeed / archParam.acq.perEch;
111}
112
113int TempLocator::getCrossSamples(int sampleNum, int& SN1, int& SN2) {
114 findGeomFromGC(sampleNum);
115 if (icross == 0 || icross >= nGalCross) return -1;
116 SN1 = crossings[icross-1];
117 SN2 = crossings[icross];
118 return 0;
119}
120
121
122void TempLocator::getAltAzBolo(int sampleNum, int ibolo, double& elv, double& az) {
123 findGeomFromGC(sampleNum);
124 double delElv = 0;
125 double delAz = 0; // relative to ch1-bolo1
126 elv = -99999;
127 az = -99999;
128 switch (ibolo) {
129 case 11:
130 delElv = 0;
131 delAz = 0;
132 break;
133 case 8:
134 delElv = 0.78 * sqrt(3.)/2.;
135 delAz = 0.78 * 1./2.;
136 break;
137 case 13:
138 delElv = - 0.78 * sqrt(3.)/2.;
139 delAz = 0.78 * 1./2.;
140 break;
141 case 9:
142 delElv = 0.78 * sqrt(3.)/2.;
143 delAz = - 0.78 * 1./2.;
144 break;
145 case 4:
146 delElv = 0.;
147 delAz = 0.78;
148 break;
149 case 15:
150 delElv = 0.;
151 delAz = - 0.78;
152 break;
153 default:
154 return;
155 }
156 delAz /= cos(41 * M_PI/180);
157 elv = altbolo1 + delElv;
158 az = azimBolGC + delAz;
159 if (az>360) az -= 360;
160 if (az<0) az += 360;
161 return;
162}
163
164double TempLocator::getAzimutSST(int sampleNum) {
165 findGeomFromGC(sampleNum);
166 double az = azimBolGC;
167 if (az>360) az -= 360;
168 if (az<0) az += 360;
169 return az;
170}
171
172double TempLocator::getElvSST(int /*sampleNum*/) {
173 return GondolaGeom::elevSST0;
174}
175
176double TempLocator::getAlphaSST(int sampleNum) {
177 double elv = getElvSST(sampleNum);
178 double az = getAzimutSST(sampleNum);
179 double ha;
180 double ra,dec;
181 aa_hadec (lat * M_PI/180, elv * M_PI/180, az * M_PI/180, &ha, &dec);
182 ra = - (ha * 180. / M_PI / 15) + (ts/3600.);
183 while (ra>24) ra -= 24;
184 while (ra<0) ra += 24;
185 //dec = dec * 180. / M_PI;
186 return ra;
187}
188
189double TempLocator::getDeltaSST(int sampleNum) {
190 double elv = getElvSST(sampleNum);
191 double az = getAzimutSST(sampleNum);
192 double ha;
193 //double ra;
194 double dec;
195 aa_hadec (lat * M_PI/180, elv * M_PI/180, az * M_PI/180, &ha, &dec);
196 //ra = - (ha * 180. / M_PI / 15) + (ts/3600.);
197 //while (ra>24) ra -= 24;
198 //while (ra<0) ra += 24;
199 dec = dec * 180. / M_PI;
200 return dec;
201}
202
203
204
205double TempLocator::getAzimutBolo(int sampleNum, int ibolo) {
206 double elv, az;
207 getAltAzBolo(sampleNum, ibolo, elv, az);
208 return az;
209}
210
211double TempLocator::getElvBolo(int sampleNum, int ibolo) {
212 double elv, az;
213 getAltAzBolo(sampleNum, ibolo, elv, az);
214 return elv;
215}
216
217double TempLocator::getAlphaBolo(int sampleNum, int ibolo) {
218 double elv, az;
219 getAltAzBolo(sampleNum, ibolo, elv, az);
220 double ha;
221 double ra,dec;
222 aa_hadec (lat * M_PI/180, elv * M_PI/180, az * M_PI/180, &ha, &dec);
223 ra = - (ha * 180. / M_PI / 15) + (ts/3600.);
224 while (ra>24) ra -= 24;
225 while (ra<0) ra += 24;
226 //dec = dec * 180. / M_PI;
227 return ra;
228}
229
230double TempLocator::getDeltaBolo(int sampleNum, int ibolo) {
231 double elv, az;
232 getAltAzBolo(sampleNum, ibolo, elv, az);
233 double ha;
234 //double ra;
235 double dec;
236 aa_hadec (lat * M_PI/180, elv * M_PI/180, az * M_PI/180, &ha, &dec);
237 //ra = - (ha * 180. / M_PI / 15) + (ts/3600.);
238 //while (ra>24) ra -= 24;
239 //while (ra<0) ra += 24;
240 dec = dec * 180. / M_PI;
241 return dec;
242}
243
244double TempLocator::getAlphaCenter(int sampleNum) {
245 return getAlphaBolo(sampleNum, 11);
246}
247
248double TempLocator::getDeltaCenter(int sampleNum) {
249 return getDeltaBolo(sampleNum, 11);
250}
251
252double TempLocator::getAzimutCenter(int sampleNum) {
253 return getAzimutBolo(sampleNum, 11);
254}
255
256double TempLocator::getElvCenter(int sampleNum) {
257 return getElvBolo(sampleNum, 11);
258}
259
260
261
262
263 // ;
264 // ; /\ 4 2
265 // ; elevation || positive scanning, clockwise
266 // ; || 6 1 5 ---------->
267 // ; || == positive azimut
268 // ; || x 3
269 // |----| = 0.78 deg / cos(elev)
270 // bol 1 = 41° elevation
271
Note: See TracBrowser for help on using the repository browser.