1 | //////////////////////////////// |
---|
2 | // Classe Astro // |
---|
3 | // Franck RICHARD // |
---|
4 | // franckrichard033@gmail.com // |
---|
5 | // BAORadio // |
---|
6 | // juin 2011 // |
---|
7 | //////////////////////////////// |
---|
8 | |
---|
9 | |
---|
10 | #include "astro.h" |
---|
11 | |
---|
12 | |
---|
13 | /************************************************************************************** |
---|
14 | ** Constructeur de la classe Astro |
---|
15 | ** |
---|
16 | ***************************************************************************************/ |
---|
17 | |
---|
18 | Astro::Astro() |
---|
19 | { |
---|
20 | //Initialisation des variables |
---|
21 | |
---|
22 | Longitude = 0.0; // Cordonnées du lieu d'observation |
---|
23 | Latitude = 0.0; |
---|
24 | |
---|
25 | Annee = 0.0; // Date et heure de l'observation |
---|
26 | Mois = 0.0; |
---|
27 | Jour = 0.0; |
---|
28 | Heure = 0.0; |
---|
29 | Min = 0.0; |
---|
30 | Sec = 0.0; |
---|
31 | UTCP = 0.0; // permet de gérer le fuseau horaire |
---|
32 | // mais semble inutile sous Linux |
---|
33 | |
---|
34 | hs = 0.0; // fraction du jour en cours avec hs=(heure+min/60+secs/3600)/24. |
---|
35 | ep = 0.0; // obliquité de l'écliptique |
---|
36 | tsl = 0.0; // temps sidéral local |
---|
37 | JJ = 0.0; // Jour Julien |
---|
38 | |
---|
39 | CorP = 0.0; // Nutation en longitude |
---|
40 | CorEP = 0.0; // Nutation en inclinaison |
---|
41 | |
---|
42 | Pression = 1013.0; |
---|
43 | Temp = 0.0; |
---|
44 | |
---|
45 | LongitudeSoleil = 0.0; |
---|
46 | } |
---|
47 | |
---|
48 | |
---|
49 | /************************************************************************************** |
---|
50 | ** Constructeur de la classe Astro |
---|
51 | ** |
---|
52 | ***************************************************************************************/ |
---|
53 | |
---|
54 | Astro::~Astro() |
---|
55 | { |
---|
56 | } |
---|
57 | |
---|
58 | |
---|
59 | /************************************************************************************** |
---|
60 | ** Initialisation des variables globales de la classe Astro |
---|
61 | ** |
---|
62 | ***************************************************************************************/ |
---|
63 | |
---|
64 | void Astro::DefinirDateHeure(double A, double M, double J, double H, double Mi, double S) |
---|
65 | { |
---|
66 | Annee = A; |
---|
67 | Mois = M; |
---|
68 | Jour = J; |
---|
69 | Heure = H; |
---|
70 | Min = Mi; |
---|
71 | Sec = S; |
---|
72 | } |
---|
73 | |
---|
74 | void Astro::DefinirPressionTemp(double P, double T) |
---|
75 | { |
---|
76 | Pression = P; |
---|
77 | Temp = T; |
---|
78 | } |
---|
79 | |
---|
80 | void Astro::DefinirLongitudeLatitude(double log, double lat) |
---|
81 | { |
---|
82 | Longitude = log; |
---|
83 | Latitude = lat; |
---|
84 | } |
---|
85 | |
---|
86 | |
---|
87 | |
---|
88 | /************************************************************************************** |
---|
89 | ** ParamÚtre : angle en radians |
---|
90 | ** retourne angle dans un intervalle 0 <= Angle <= 2*PI |
---|
91 | ***************************************************************************************/ |
---|
92 | |
---|
93 | double Astro::VerifAngle(double Angle) |
---|
94 | { |
---|
95 | Angle=fmod(Angle, Pi2); |
---|
96 | if (Angle<0.0) Angle+=Pi2; |
---|
97 | |
---|
98 | return(Angle); |
---|
99 | } |
---|
100 | |
---|
101 | /************************************************************************************** |
---|
102 | ** ParamÚtre : angle en radians |
---|
103 | ** retourne angle dans un intervalle 0 <= Angle <= PI |
---|
104 | ***************************************************************************************/ |
---|
105 | |
---|
106 | double Astro::VerifDistance(double Angle) |
---|
107 | { |
---|
108 | Angle=VerifAngle(Angle); |
---|
109 | if (Angle>Pi) Angle=Pi2-Angle; |
---|
110 | |
---|
111 | return(Angle); |
---|
112 | } |
---|
113 | |
---|
114 | |
---|
115 | /************************************************************************************** |
---|
116 | ** Distance entre deux points situés sur une sphÚre |
---|
117 | ** |
---|
118 | ***************************************************************************************/ |
---|
119 | |
---|
120 | double Astro::DistanceAngulaireEntre2Points(double az1, double ha1, double az2, double ha2) |
---|
121 | { |
---|
122 | return acos(sin(ha1) * sin(ha2) + cos(ha1) * cos(ha2) * cos(az2 - az1)); |
---|
123 | } |
---|
124 | |
---|
125 | |
---|
126 | /************************************************************************************** |
---|
127 | ** retourne l'arrondi du nombre donné en paramÚtre |
---|
128 | ** Exemples: Arrondi(5.7) = 6 |
---|
129 | ** Arrondi(-5.7)= -6 |
---|
130 | ***************************************************************************************/ |
---|
131 | |
---|
132 | double Astro::Arrondi(double a) |
---|
133 | { |
---|
134 | double b; |
---|
135 | |
---|
136 | if (a>=0.0) |
---|
137 | { |
---|
138 | b=a-floor(a); |
---|
139 | if (b>=0.5) return ceil(a); |
---|
140 | else return floor(a); |
---|
141 | } |
---|
142 | else |
---|
143 | { |
---|
144 | b=ceil(a)-a; |
---|
145 | if (b>=0.5) return floor(a); |
---|
146 | else return ceil(a); |
---|
147 | } |
---|
148 | } |
---|
149 | |
---|
150 | |
---|
151 | |
---|
152 | /************************************************************************************** |
---|
153 | ** converti un angle (exprimé en degrés) dans le format hh:mm:ss ou deg:mm:ss |
---|
154 | ** selon la valeur de HMS |
---|
155 | ***************************************************************************************/ |
---|
156 | |
---|
157 | string Astro::DHMS(double mema, bool HMS) |
---|
158 | { |
---|
159 | double d, m, s, a; |
---|
160 | |
---|
161 | a=mema; |
---|
162 | |
---|
163 | if (HMS) a = a / 15.0; |
---|
164 | a=fabs(a * 3600.0); |
---|
165 | s=Arrondi(fmod(a, 60.0) * 10.0) / 10.0; |
---|
166 | a=floor(a / 60.0); |
---|
167 | m=fmod(a, 60.0); |
---|
168 | d=floor(a / 60.0); |
---|
169 | |
---|
170 | if (s == 60.0) |
---|
171 | { |
---|
172 | s=0.0; |
---|
173 | m+=1.0; |
---|
174 | if (m==60.0) |
---|
175 | { |
---|
176 | m=0.0; |
---|
177 | d+=1.0; |
---|
178 | |
---|
179 | if (HMS && d>23.0) d=0.0; |
---|
180 | } |
---|
181 | } |
---|
182 | |
---|
183 | stringstream strsTemp; |
---|
184 | |
---|
185 | if (mema < 0) strsTemp << '-'; |
---|
186 | |
---|
187 | strsTemp << setfill('0') << setw(2) << d << ":" << setfill('0') << setw(2) << m << ":" << setfill('0') << setw(2) << s; |
---|
188 | |
---|
189 | return (strsTemp.str()); |
---|
190 | } |
---|
191 | |
---|
192 | /************************************************************************************** |
---|
193 | ** Décomposition et vérification des dates, heures, AD et déclinaisons |
---|
194 | ** |
---|
195 | ** type = 0 -> chaîne contient une déclinaison ou une latitude |
---|
196 | ** type = 1 -> " " une longitude |
---|
197 | ** type = 2 -> " " une AD ou une heure |
---|
198 | ** type = 3 -> " " une date |
---|
199 | ** |
---|
200 | ** Exemple de Décomposition("15:23:12", 2, a, b, c) retourne a=15, b=23, c=12 |
---|
201 | ** si la chaîne a un format incorrect, la fct retourne false |
---|
202 | ** |
---|
203 | ***************************************************************************************/ |
---|
204 | |
---|
205 | bool Astro::Decomposition(string chaine, char type, float *a1, float *a2, float *a3) |
---|
206 | { |
---|
207 | string car, s; |
---|
208 | |
---|
209 | float a, b, c; |
---|
210 | |
---|
211 | // pour les heures et les coordonnées, on attend ":" comme caractÚre séparateur, sinon |
---|
212 | // on attend d'avoir des caractÚres "/" entre chaque éléments pour une date |
---|
213 | |
---|
214 | (type == 3) ? car="/" : car=":"; |
---|
215 | |
---|
216 | // Y a-t-il 2 caractÚres ':' ou '/' dans la chaîne ? |
---|
217 | // C'est indispensable dans tous les cas |
---|
218 | |
---|
219 | int test = 0; |
---|
220 | |
---|
221 | for (size_t i=0; i<chaine.length(); i++) if (chaine[i] == car[0]) test++; |
---|
222 | |
---|
223 | if (test<2) return false; |
---|
224 | |
---|
225 | // Extraction des trois nombres |
---|
226 | |
---|
227 | s = chaine.substr(0, chaine.find(car)); |
---|
228 | |
---|
229 | a = atoi(s.c_str()); |
---|
230 | |
---|
231 | s = chaine.substr(chaine.find(car) + 1, chaine.rfind(car) - chaine.find(car) - 1); |
---|
232 | |
---|
233 | b = atoi(s.c_str()); |
---|
234 | |
---|
235 | s = chaine.substr(chaine.rfind(car)+1); |
---|
236 | |
---|
237 | // 12.125454414 sec -> 12.125 sec |
---|
238 | |
---|
239 | c = Arrondi(atof(s.c_str()) * 100.0) / 100.0; |
---|
240 | |
---|
241 | //vérification de la cohérence des infos contenues dans la chaîne |
---|
242 | |
---|
243 | if (type < 3 ) |
---|
244 | { |
---|
245 | // pour une déclinaison |
---|
246 | if (( type == 0 ) && ( a>90.0 || a<-90.0 )) return false; |
---|
247 | // pour une AD |
---|
248 | if (( type == 1 ) && ( a>360.0 || a<0.0 )) return false; |
---|
249 | // pour une heure |
---|
250 | if (( type == 2 ) && ( a>23.0 || a<0.0 )) return false; |
---|
251 | // pour les minutes |
---|
252 | if ( b < 0.0 || b> 59.0 ) return false; |
---|
253 | //pour les secondes |
---|
254 | if ( c < 0.0 || c>=60.0 ) return false; |
---|
255 | } |
---|
256 | else |
---|
257 | { |
---|
258 | //pour les jours |
---|
259 | if ( a < 0.0 || a > 31.0 ) return false; |
---|
260 | //pour les mois |
---|
261 | if ( b < 0.0 || b > 12.0 ) return false; |
---|
262 | } |
---|
263 | |
---|
264 | if ( a1 != NULL ) *a1 = a; |
---|
265 | if ( a2 != NULL ) *a2 = b; |
---|
266 | if ( a3 != NULL ) *a3 = c; |
---|
267 | |
---|
268 | return true; |
---|
269 | } |
---|
270 | |
---|
271 | /************************************************************************************** |
---|
272 | ** Calcul de la longitude du Soleil |
---|
273 | ** La quantité PS n'est pas indispensable. Elle permet d'améliorer la précision |
---|
274 | ***************************************************************************************/ |
---|
275 | |
---|
276 | double Astro::CalculLongitudeSoleil() |
---|
277 | { |
---|
278 | double T = (JJ - 2415020.0) / 36525.0; |
---|
279 | |
---|
280 | double A = 153.23 + 22518.7541 * T; |
---|
281 | double B = 216.57 + 44037.5082 * T; |
---|
282 | double C = 312.69 + 32964.3577 * T; |
---|
283 | double D = 350.74 + ( 445267.1142 - 0.00144 * T ) * T; |
---|
284 | double E = 231.19 + 20.2 * T; |
---|
285 | |
---|
286 | |
---|
287 | double L = 279.69668 + ( 36000.76892 + 0.0003025 * T ) * T; |
---|
288 | double M = 358.47583 + (((( -0.0000033 * T) -0.00015 ) * T ) + 35999.04975 ) * T; |
---|
289 | |
---|
290 | double CC = (1.919460 + ((-0.000014 * T) - 0.004789 ) * T ) * sin(M * Pidiv180) |
---|
291 | + (0.020094 - 0.0001 * T) * sin(2.0 * M * Pidiv180) |
---|
292 | + 0.000293 * sin(3.0 * M * Pidiv180); |
---|
293 | |
---|
294 | double PS = 0.00134 * cos(A * Pidiv180) |
---|
295 | +0.00154 * cos(B * Pidiv180) |
---|
296 | +0.00200 * cos(C * Pidiv180) |
---|
297 | +0.00179 * sin(D * Pidiv180) |
---|
298 | +0.00178 * sin(E * Pidiv180); |
---|
299 | |
---|
300 | |
---|
301 | double LS = VerifAngle( (L + CC + PS) * Pidiv180); |
---|
302 | |
---|
303 | return LS; |
---|
304 | } |
---|
305 | |
---|
306 | |
---|
307 | |
---|
308 | /************************************************************************************** |
---|
309 | ** Calcul de l'AD et la déclinaison du Soleil |
---|
310 | ** indispensable pour exécuter la commande "goto sun" |
---|
311 | ***************************************************************************************/ |
---|
312 | |
---|
313 | void Astro::CalculARDecSoleil(CoordonneesHorairesDouble *Soleil) |
---|
314 | { |
---|
315 | double ar; |
---|
316 | double de; |
---|
317 | |
---|
318 | ar = atan( cos(ep) * sin(LongitudeSoleil) / cos(LongitudeSoleil)); |
---|
319 | |
---|
320 | if ( cos(LongitudeSoleil) < 0.0) ar += Pi; |
---|
321 | |
---|
322 | de = asin( sin(ep) * sin(LongitudeSoleil)); |
---|
323 | |
---|
324 | Soleil->ar = VerifAngle(ar); |
---|
325 | |
---|
326 | Soleil->dec = de; |
---|
327 | } |
---|
328 | |
---|
329 | |
---|
330 | |
---|
331 | |
---|
332 | /************************************************************************************** |
---|
333 | ** Calcul de la nutation en longitude et en inclinaison |
---|
334 | ** Nécessaire pour calculer la position apparente d'un objet (étoile, galaxie etc...) |
---|
335 | ***************************************************************************************/ |
---|
336 | |
---|
337 | void Astro::Nutation() |
---|
338 | { |
---|
339 | double T = (JJ-2415020.0) / 36525.0; |
---|
340 | |
---|
341 | double L = (279.6967 + (36000.7689 + 0.000303 * T ) * T ) * Pidiv180; |
---|
342 | double LP = (270.4342 + (481267.8831 - 0.001133 * T ) * T ) * Pidiv180; |
---|
343 | double M = (358.4758 + (35999.0498 - 0.000150 * T ) * T ) * Pidiv180; |
---|
344 | double MP = (296.1046 + (477198.8491 + 0.009192 * T ) * T ) * Pidiv180; |
---|
345 | double O = (259.1833 + (-1934.1420 + 0.002078 * T ) * T ) * Pidiv180; |
---|
346 | |
---|
347 | double L2 = 2.0 * L; |
---|
348 | double LP2 = 2.0 * LP; |
---|
349 | double O2 = 2.0 * O; |
---|
350 | |
---|
351 | // Nutation en longitude |
---|
352 | |
---|
353 | CorP = -(17.2327 + 0.01737 * T) * sin(O) |
---|
354 | -(1.2729+0.00013 * T) * sin(L2) |
---|
355 | +0.2088 * sin(O2) |
---|
356 | -0.2037 * sin(LP2) |
---|
357 | +(0.1261 - 0.00031 * T) * sin(M) |
---|
358 | +0.0675 * sin(MP) |
---|
359 | -(0.0497 - 0.00012 * T) * sin(L2 + M) |
---|
360 | -0.0342 * sin(LP2 - O) |
---|
361 | -0.0261 * sin(LP2 + MP) |
---|
362 | +0.0214 * sin(L2 - M) |
---|
363 | -0.0149 * sin(L2 - LP2 + MP) |
---|
364 | +0.0124 * sin(L2 - O) |
---|
365 | +0.0114 * sin(LP2 - MP); |
---|
366 | |
---|
367 | //Nutation en inclinaison |
---|
368 | |
---|
369 | CorEP = (9.2100 + 0.00091 * T) * cos(O) |
---|
370 | +(0.5522 - 0.00029 * T) * cos(L2) |
---|
371 | -0.0904 * cos(O2) |
---|
372 | +0.0884 * cos(LP2) |
---|
373 | +0.0216 * cos(L2+M) |
---|
374 | +0.0183 * cos(LP2-O) |
---|
375 | +0.0113 * cos(LP2+MP) |
---|
376 | -0.0093 * cos(L2-M) |
---|
377 | -0.0066 * cos(L2-O); |
---|
378 | |
---|
379 | CorP = CorP / 3600.0 * Pidiv180; |
---|
380 | CorEP = CorEP / 3600.0 * Pidiv180; |
---|
381 | } |
---|
382 | |
---|
383 | |
---|
384 | |
---|
385 | /************************************************************************************** |
---|
386 | ** Nutation des étoiles - on utilise les quantités calculées précédemment |
---|
387 | ** ar et de sont exprimés en radians |
---|
388 | ***************************************************************************************/ |
---|
389 | |
---|
390 | void Astro::NutationEtoile(double *ar, double *de) |
---|
391 | { |
---|
392 | double a = (cos(ep) + sin(ep) * sin(*ar) * tan(*de)) * CorP - (cos(*ar) * tan(*de)) * CorEP; |
---|
393 | double b = (sin(ep) * cos(*ar)) * CorP + sin(*ar) * CorEP; |
---|
394 | |
---|
395 | *ar = VerifAngle(*ar+a); |
---|
396 | *de += b; |
---|
397 | } |
---|
398 | |
---|
399 | |
---|
400 | /************************************************************************************** |
---|
401 | ** Précession des équinoxes |
---|
402 | ** ar et de sont exprimés en radians |
---|
403 | ***************************************************************************************/ |
---|
404 | |
---|
405 | void Astro::Precession(double *ar, double *de) |
---|
406 | { |
---|
407 | double t = (JJ - 2451544.5) / 36524.2199; |
---|
408 | |
---|
409 | double ze = ((((0.018 * t) + 0.302) * t + 2304.948) * t) / 3600.0 * Pidiv180; |
---|
410 | double z = ((((0.019 * t) + 1.093) * t + 2304.948) * t) / 3600.0 * Pidiv180; |
---|
411 | double delta = ((((-0.042 * t) - 0.426) * t + 2004.255) * t) / 3600.0 * Pidiv180; |
---|
412 | |
---|
413 | |
---|
414 | double A = cos(*de) * sin(*ar + ze); |
---|
415 | double B = cos(delta) * cos(*de) * cos (*ar + ze) - sin(delta) * sin(*de); |
---|
416 | double C = sin(delta) * cos(*de) * cos (*ar + ze) + cos(delta) * sin(*de); |
---|
417 | |
---|
418 | double AMZ = atan(A / B); |
---|
419 | |
---|
420 | if (B<0.0) AMZ += Pi; |
---|
421 | |
---|
422 | *ar=VerifAngle(AMZ + z); |
---|
423 | *de=asin(C); |
---|
424 | } |
---|
425 | |
---|
426 | |
---|
427 | |
---|
428 | /************************************************************************************** |
---|
429 | ** Obliquité |
---|
430 | ** calcul de l'inclinaison de l'axe terrestre par rapport à au plan de l'écliptique |
---|
431 | ** Quantité indispensable pour calculer le temps sidéral local et les coordonnées |
---|
432 | ** horaires du soleil |
---|
433 | ***************************************************************************************/ |
---|
434 | |
---|
435 | void Astro::Obliquite(double JJ) |
---|
436 | { |
---|
437 | double T; |
---|
438 | |
---|
439 | T = (JJ - 2415020.0) / 36525.0; |
---|
440 | |
---|
441 | ep = (23.452294+ (((0.000000503 * T ) - 0.00000164 ) * T - 0.0130125 ) * T ) * Pidiv180; |
---|
442 | } |
---|
443 | |
---|
444 | |
---|
445 | |
---|
446 | |
---|
447 | /************************************************************************************** |
---|
448 | ** Aberration annuelle des étoiles |
---|
449 | ** ar et de sont exprimés en radians |
---|
450 | ***************************************************************************************/ |
---|
451 | |
---|
452 | void Astro::AberrationAnnuelle(double *ar, double *de) |
---|
453 | { |
---|
454 | double c = -20.49 / 3600.0 * Pidiv180; |
---|
455 | double a = c * ( cos(*ar) * cos(LongitudeSoleil) * cos(ep) |
---|
456 | + sin(*ar) * sin(LongitudeSoleil) |
---|
457 | ) / cos(*de); |
---|
458 | |
---|
459 | double b = c * ( cos(LongitudeSoleil) * cos(ep) * (tan(ep) * cos(*de) - sin(*ar) * sin(*de)) |
---|
460 | + cos(*ar) * sin(*de) * sin(LongitudeSoleil)); |
---|
461 | |
---|
462 | *ar = VerifAngle(*ar + a); |
---|
463 | *de += b; |
---|
464 | } |
---|
465 | |
---|
466 | |
---|
467 | /************************************************************************************** |
---|
468 | ** calcul Jour julien |
---|
469 | ***************************************************************************************/ |
---|
470 | |
---|
471 | void Astro::CalculJJ(double Heure) |
---|
472 | { |
---|
473 | JJ = CalculJJ(Annee, Mois, Jour, Heure); |
---|
474 | } |
---|
475 | |
---|
476 | double Astro::CalculJJ(double A, double M, double J, double Heure) |
---|
477 | { |
---|
478 | long y, a, b, c, e, m; |
---|
479 | |
---|
480 | long year = (long)A; |
---|
481 | int month = (int)M; |
---|
482 | double day = J + Heure; |
---|
483 | |
---|
484 | y = year + 4800; |
---|
485 | |
---|
486 | m = month; |
---|
487 | if ( m <= 2 ) |
---|
488 | { |
---|
489 | m += 12; |
---|
490 | y -= 1; |
---|
491 | } |
---|
492 | e = (306 * (m+1))/10; |
---|
493 | |
---|
494 | a = y/100; |
---|
495 | if ( year <= 1582L ) |
---|
496 | { |
---|
497 | if ( year == 1582L ) |
---|
498 | { |
---|
499 | if ( month < 10 ) |
---|
500 | goto julius; |
---|
501 | if ( month > 10) |
---|
502 | goto gregor; |
---|
503 | if ( day >= 15 ) |
---|
504 | goto gregor; |
---|
505 | } |
---|
506 | julius: |
---|
507 | |
---|
508 | b = -38; |
---|
509 | } |
---|
510 | else |
---|
511 | { |
---|
512 | |
---|
513 | gregor: |
---|
514 | b = (a / 4) - a; |
---|
515 | } |
---|
516 | |
---|
517 | c = (36525L * y) / 100; |
---|
518 | |
---|
519 | return b + c + e + day - 32167.5; |
---|
520 | } |
---|
521 | |
---|
522 | |
---|
523 | |
---|
524 | /************************************************************************************** |
---|
525 | ** Temps sidéral local |
---|
526 | ***************************************************************************************/ |
---|
527 | |
---|
528 | double Astro::TSL(double JJ, double HeureSiderale, double Longitude) |
---|
529 | { |
---|
530 | double T = (JJ - ET_UT / 3600.0 / 24.0 - 2415020.0 ) / 36525.0; |
---|
531 | double rd = 0.276919398 + ( 100.0021359 + 0.000001075 * T ) * T; |
---|
532 | rd += HeureSiderale; |
---|
533 | rd *= Pi2; |
---|
534 | // temps sidéral apparent |
---|
535 | rd += CorP * cos(ep); |
---|
536 | rd -= Longitude; |
---|
537 | |
---|
538 | return VerifAngle(rd); |
---|
539 | } |
---|
540 | |
---|
541 | |
---|
542 | /************************************************************************************** |
---|
543 | ** routine principale des calculs |
---|
544 | ***************************************************************************************/ |
---|
545 | |
---|
546 | void Astro::CalculTSL() |
---|
547 | { |
---|
548 | hs = (Heure + Min / 60.0 + Sec / 3600.0) / 24.0; |
---|
549 | |
---|
550 | hs -= UTCP / 24.0; |
---|
551 | |
---|
552 | CalculJJ(hs); |
---|
553 | |
---|
554 | Obliquite(JJ); |
---|
555 | |
---|
556 | Nutation(); |
---|
557 | |
---|
558 | tsl = TSL(JJ, hs, Longitude); |
---|
559 | |
---|
560 | LongitudeSoleil = CalculLongitudeSoleil(); |
---|
561 | } |
---|
562 | |
---|
563 | |
---|
564 | |
---|
565 | /************************************************************************************** |
---|
566 | ** Calcule la hauteur et l'azimut des étoiles en fct du lieu d'observation |
---|
567 | ***************************************************************************************/ |
---|
568 | |
---|
569 | void Astro::Azimut(double Ar, double De, double *azi, double *hau) |
---|
570 | { |
---|
571 | double ah = tsl - Ar; |
---|
572 | |
---|
573 | double cosLatitude = cos(Latitude); |
---|
574 | double sinLatitude = sin(Latitude); |
---|
575 | double cosDe_cosAh = cos(De) * cos(ah); |
---|
576 | double sinDe = sin(De); |
---|
577 | |
---|
578 | |
---|
579 | double zc = sinLatitude * sinDe + cosLatitude * cosDe_cosAh; |
---|
580 | double ht = atan(zc / sqrt((-zc * zc) + 1.0)); |
---|
581 | |
---|
582 | double cosHt = cos(ht); |
---|
583 | |
---|
584 | double a1 = cos(De) * sin(ah) / cosHt; |
---|
585 | double ac = (-cosLatitude * sinDe + sinLatitude * cosDe_cosAh) / cosHt; |
---|
586 | double az = atan(a1 / sqrt((-a1*a1)+1.0)); |
---|
587 | |
---|
588 | if (ac<0.0) az = Pi - az; |
---|
589 | |
---|
590 | *azi = VerifAngle( Pi + az ); |
---|
591 | *hau = ht; |
---|
592 | } |
---|
593 | |
---|
594 | |
---|
595 | /************************************************************************************** |
---|
596 | ** Calcule l'ascension droite et la déclinaison d'un objet situé à l'azimut az |
---|
597 | ** et à la déclinaison de |
---|
598 | ***************************************************************************************/ |
---|
599 | |
---|
600 | void Astro::AzHa2ADDe(double az, double ha, double *AD, double *Dec) |
---|
601 | { |
---|
602 | double sindel, del, habar, cosha, ha10; |
---|
603 | |
---|
604 | sindel = sin(ha) * sin(Latitude) + cos(ha) * cos(Latitude) * cos(az); |
---|
605 | |
---|
606 | del = asin(sindel); |
---|
607 | |
---|
608 | if ( (sindel == 1.0) || ( sindel == -1.0) ) |
---|
609 | { |
---|
610 | habar = 0.0; |
---|
611 | } |
---|
612 | else |
---|
613 | { |
---|
614 | if ( (sin(Latitude) == 1.0) || (sin(Latitude) == -1.0)) |
---|
615 | { |
---|
616 | habar = az; |
---|
617 | } |
---|
618 | else |
---|
619 | { |
---|
620 | cosha = ( sin(ha) - sin(Latitude) * sindel ) / ( cos(Latitude) * cos(del) ); |
---|
621 | if (cosha > 1.0 ) cosha = 1.0; |
---|
622 | if (cosha < -1.0 ) cosha = -1.0; |
---|
623 | |
---|
624 | habar = acos(cosha); |
---|
625 | } |
---|
626 | } |
---|
627 | |
---|
628 | if (sin(az) < 0.0) |
---|
629 | { |
---|
630 | ha10 = habar; |
---|
631 | } |
---|
632 | else |
---|
633 | { |
---|
634 | ha10 = Pi2 - habar ; |
---|
635 | } |
---|
636 | |
---|
637 | *AD = VerifAngle( tsl - ha10 ); |
---|
638 | |
---|
639 | // if (*AD < 0.0) *AD += Pi2; |
---|
640 | |
---|
641 | *Dec = del; |
---|
642 | } |
---|
643 | |
---|
644 | |
---|
645 | /************************************************************************************** |
---|
646 | ** Recherche l'azimut de l'objet suivi au moment où celui-ci passera sous les 30° |
---|
647 | ** et deviendra donc hors d'atteinte de l'instrument. |
---|
648 | ***************************************************************************************/ |
---|
649 | |
---|
650 | void Astro::RechercheAzimutFinSuivi(double AD, double Dec, int posazactuelle, int *azmin, int *azmax) |
---|
651 | { |
---|
652 | double az, ha; |
---|
653 | int min = 30000; |
---|
654 | int max = -30000; |
---|
655 | int delta; |
---|
656 | int codeuraz; |
---|
657 | |
---|
658 | CalculTSL(); |
---|
659 | |
---|
660 | |
---|
661 | posazactuelle = 0; |
---|
662 | |
---|
663 | // la recherche se fait pour les 12 prochaines heures (12 * 60 = 720) |
---|
664 | |
---|
665 | for (int i = 0; i < 720; i++) |
---|
666 | { |
---|
667 | // calcul de l'azimut de l'objet i minutes plus tard |
---|
668 | |
---|
669 | Azimut( AD, Dec, &az, &ha); |
---|
670 | |
---|
671 | if ( ha < HAUTMIN * Pidiv180 ) break; |
---|
672 | |
---|
673 | az = VerifAngle(az + Pi); |
---|
674 | |
---|
675 | codeuraz = (int) Arrondi( az * (double)NBREPASCODEURSAZ / Pi2 ); |
---|
676 | |
---|
677 | delta = codeuraz - posazactuelle; |
---|
678 | |
---|
679 | if ( delta >= NBREPASCODEURSAZ / 2 ) delta -= NBREPASCODEURSAZ; |
---|
680 | |
---|
681 | if ( delta <= -NBREPASCODEURSAZ / 2 ) delta += NBREPASCODEURSAZ; |
---|
682 | |
---|
683 | posazactuelle += delta; |
---|
684 | |
---|
685 | if ( posazactuelle < min) min = posazactuelle; |
---|
686 | |
---|
687 | if ( posazactuelle > max) max = posazactuelle; |
---|
688 | |
---|
689 | tsl += 1.0 / 60.0 / 24.0 * Pi2 ; |
---|
690 | } |
---|
691 | |
---|
692 | *azmin = min - 300; |
---|
693 | |
---|
694 | *azmax = max + 300; |
---|
695 | |
---|
696 | CalculTSL(); |
---|
697 | } |
---|
698 | |
---|
699 | |
---|
700 | |
---|
701 | /************************************************************************************** |
---|
702 | ** réfraction atmosphérique : formule simple de Jean Meeus |
---|
703 | ** la hauteur ht est exprimée en radians |
---|
704 | ***************************************************************************************/ |
---|
705 | |
---|
706 | double Astro::RefractionAtmospherique(double ht) |
---|
707 | { |
---|
708 | double gamma, R0, a, b, PressionTempCalc, tanHT; |
---|
709 | |
---|
710 | if (Pression != 0.0 && ht > 0.0) |
---|
711 | { |
---|
712 | PressionTempCalc = (Pression / 1013.25) * 273.0 / (273.0 + Temp); |
---|
713 | |
---|
714 | if (ht <= 15.0 * Pidiv180) |
---|
715 | { |
---|
716 | gamma = 2.6; |
---|
717 | a = 7.5262; |
---|
718 | b = -2.2204; |
---|
719 | |
---|
720 | if (ht >= 4.0 * Pidiv180) |
---|
721 | { |
---|
722 | gamma = -1.1; |
---|
723 | a = 4.4010; |
---|
724 | b = -0.9603; |
---|
725 | } |
---|
726 | |
---|
727 | R0 = (pow(( a + b * log( ht * N180divPi + gamma)), 2.0)) / 60.0 * Pidiv180; |
---|
728 | ht += R0 * PressionTempCalc; |
---|
729 | } |
---|
730 | else |
---|
731 | { |
---|
732 | tanHT = fabs(tan(Pidiv2 - ht)); |
---|
733 | R0 = ( 0.0161877777777777777777 - 0.000022888888888888888 * tanHT * tanHT ) |
---|
734 | * tanHT * Pidiv180; |
---|
735 | ht += R0 * PressionTempCalc; |
---|
736 | } |
---|
737 | } |
---|
738 | |
---|
739 | return ht; |
---|
740 | } |
---|
741 | |
---|
742 | |
---|
743 | |
---|
744 | |
---|
745 | double Astro::slaDrange(double angle ) |
---|
746 | { |
---|
747 | return fmod(angle, Pi); |
---|
748 | } |
---|
749 | /* |
---|
750 | float Astro::slaRange(float angle ) |
---|
751 | { |
---|
752 | return fmodf(angle, Pi); |
---|
753 | } |
---|
754 | */ |
---|
755 | void Astro::slaRefro ( double zobs, double hm, double tdk, double pmb, |
---|
756 | double rh, double wl, double phi, double tlr, |
---|
757 | double eps, double *ref ) |
---|
758 | /* |
---|
759 | ** - - - - - - - - - |
---|
760 | ** s l a R e f r o |
---|
761 | ** - - - - - - - - - |
---|
762 | ** |
---|
763 | ** Atmospheric refraction for radio and optical wavelengths. |
---|
764 | ** |
---|
765 | ** Given: |
---|
766 | ** zobs double observed zenith distance of the source (radian) |
---|
767 | ** hm double height of the observer above sea level (metre) |
---|
768 | ** tdk double ambient temperature at the observer (deg K) |
---|
769 | ** pmb double pressure at the observer (millibar) |
---|
770 | ** rh double relative humidity at the observer (range 0-1) |
---|
771 | ** wl double effective wavelength of the source (micrometre) |
---|
772 | ** phi double latitude of the observer (radian, astronomical) |
---|
773 | ** tlr double temperature lapse rate in the troposphere (degK/met |
---|
774 | ** eps double precision required to terminate iteration (radian) |
---|
775 | ** |
---|
776 | ** Returned: |
---|
777 | ** ref double refraction: in vacuo ZD minus observed ZD (radian) |
---|
778 | ** |
---|
779 | ** Notes: |
---|
780 | ** |
---|
781 | ** 1 A suggested value for the tlr argument is 0.0065D0. The |
---|
782 | ** refraction is significantly affected by tlr, and if studies |
---|
783 | ** of the local atmosphere have been carried out a better tlr |
---|
784 | ** value may be available. |
---|
785 | ** |
---|
786 | ** 2 A suggested value for the eps argument is 1e-8. The result is |
---|
787 | ** usually at least two orders of magnitude more computationally |
---|
788 | ** precise than the supplied eps value. |
---|
789 | ** |
---|
790 | ** 3 The routine computes the refraction for zenith distances up |
---|
791 | ** to and a little beyond 90 deg using the method of Hohenkerk |
---|
792 | ** and Sinclair (NAO Technical Notes 59 and 63, subsequently adopted |
---|
793 | ** in the Explanatory Supplement, 1992 edition - see section 3.281). |
---|
794 | ** |
---|
795 | ** 4 The C code is an adaptation of the Fortran optical refraction |
---|
796 | ** subroutine AREF of C.Hohenkerk (HMNAO, September 1984), with |
---|
797 | ** extensions to support the radio case. The following modifications |
---|
798 | ** to the original HMNAO optical refraction algorithm have been made: |
---|
799 | ** |
---|
800 | ** . The angle arguments have been changed to radians. |
---|
801 | ** |
---|
802 | ** . Any value of zobs is allowed (see note 6, below). |
---|
803 | ** |
---|
804 | ** . Other argument values have been limited to safe values. |
---|
805 | ** |
---|
806 | ** . Murray's values for the gas constants have been used |
---|
807 | ** (Vectorial Astrometry, Adam Hilger, 1983). |
---|
808 | ** |
---|
809 | ** . The numerical integration phase has been rearranged for |
---|
810 | ** extra clarity. |
---|
811 | ** |
---|
812 | ** . A better model for Ps(T) has been adopted (taken from |
---|
813 | ** Gill, Atmosphere-Ocean Dynamics, Academic Press, 1982). |
---|
814 | ** |
---|
815 | ** . More accurate expressions for Pwo have been adopted |
---|
816 | ** (again from Gill 1982). |
---|
817 | ** |
---|
818 | ** . Provision for radio wavelengths has been added using |
---|
819 | ** expressions devised by A.T.Sinclair, RGO (private |
---|
820 | ** communication 1989), based on the Essen & Froome |
---|
821 | ** refractivity formula adopted in Resolution 1 of the |
---|
822 | ** 13th International Geodesy Association General Assembly |
---|
823 | ** (Bulletin Geodesique 1963 p390). |
---|
824 | ** |
---|
825 | ** . Various small changes have been made to gain speed. |
---|
826 | ** |
---|
827 | ** None of the changes significantly affects the optical results |
---|
828 | ** with respect to the algorithm given in the 1992 Explanatory |
---|
829 | ** Supplement. For example, at 70 deg zenith distance the present |
---|
830 | ** routine agrees with the ES algorithm to better than 0.05 arcsec |
---|
831 | ** for any reasonable combination of parameters. However, the |
---|
832 | ** improved water-vapour expressions do make a significant difference |
---|
833 | ** in the radio band, at 70 deg zenith distance reaching almost |
---|
834 | ** 4 arcsec for a hot, humid, low-altitude site during a period of |
---|
835 | ** low pressure. |
---|
836 | ** |
---|
837 | ** 5 The radio refraction is chosen by specifying wl > 100 micrometres. |
---|
838 | ** Because the algorithm takes no account of the ionosphere, the |
---|
839 | ** accuracy deteriorates at low frequencies, below about 30 MHz. |
---|
840 | ** |
---|
841 | ** 6 Before use, the value of zobs is expressed in the range +/- pi. |
---|
842 | ** If this ranged zobs is -ve, the result ref is computed from its |
---|
843 | ** absolute value before being made -ve to match. In addition, if |
---|
844 | ** it has an absolute value greater than 93 deg, a fixed ref value |
---|
845 | ** equal to the result for zobs = 93 deg is returned, appropriately |
---|
846 | ** signed. |
---|
847 | ** |
---|
848 | ** 7 As in the original Hohenkerk and Sinclair algorithm, fixed values |
---|
849 | ** of the water vapour polytrope exponent, the height of the |
---|
850 | ** tropopause, and the height at which refraction is negligible are |
---|
851 | ** used. |
---|
852 | ** |
---|
853 | ** 8 The radio refraction has been tested against work done by |
---|
854 | ** Iain Coulson, JACH, (private communication 1995) for the |
---|
855 | ** James Clerk Maxwell Telescope, Mauna Kea. For typical conditions, |
---|
856 | ** agreement at the 0.1 arcsec level is achieved for moderate ZD, |
---|
857 | ** worsening to perhaps 0.5-1.0 arcsec at ZD 80 deg. At hot and |
---|
858 | ** humid sea-level sites the accuracy will not be as good. |
---|
859 | ** |
---|
860 | ** 9 It should be noted that the relative humidity rh is formally |
---|
861 | ** defined in terms of "mixing ratio" rather than pressures or |
---|
862 | ** densities as is often stated. It is the mass of water per unit |
---|
863 | ** mass of dry air divided by that for saturated air at the same |
---|
864 | ** temperature and pressure (see Gill 1982). |
---|
865 | ** |
---|
866 | ** Called: slaDrange, atmt, atms |
---|
867 | ** |
---|
868 | ** Defined in slamac.h: TRUE, FALSE |
---|
869 | ** |
---|
870 | ** Last revision: 30 January 1997 |
---|
871 | ** |
---|
872 | ** Copyright P.T.Wallace. All rights reserved. |
---|
873 | */ |
---|
874 | { |
---|
875 | /* Fixed parameters */ |
---|
876 | |
---|
877 | static double d93 = 1.623156204; /* 93 degrees in radians */ |
---|
878 | static double gcr = 8314.32; /* Universal gas constant */ |
---|
879 | static double dmd = 28.9644; /* Molecular weight of dry air */ |
---|
880 | static double dmw = 18.0152; /* Molecular weight of water |
---|
881 | vapour */ |
---|
882 | static double s = 6378120.0; /* Mean Earth radius (metre) */ |
---|
883 | static double delta = 18.36; /* Exponent of temperature |
---|
884 | dependence of water vapour |
---|
885 | pressure */ |
---|
886 | static double ht = 11000.0; /* Height of tropopause (metre) */ |
---|
887 | static double hs = 80000.0; /* Upper limit for refractive |
---|
888 | effects (metre) */ |
---|
889 | |
---|
890 | /* Variables used when calling the internal routine atmt */ |
---|
891 | double robs; /* height of observer from centre of Earth (metre) */ |
---|
892 | double tdkok; /* temperature at the observer (deg K) */ |
---|
893 | double alpha; /* alpha | */ |
---|
894 | double gamm2; /* gamma minus 2 | see ES */ |
---|
895 | double delm2; /* delta minus 2 | */ |
---|
896 | double c1,c2,c3,c4,c5,c6; /* various */ |
---|
897 | |
---|
898 | /* Variables used when calling the internal routine atms */ |
---|
899 | double rt; /* height of tropopause from centre of Earth (metre) */ |
---|
900 | double tt; /* temperature at the tropopause (deg k) */ |
---|
901 | double dnt; /* refractive index at the tropopause */ |
---|
902 | double gamal; /* constant of the atmospheric model = g*md/r */ |
---|
903 | |
---|
904 | int is, k, n, i, j, optic; |
---|
905 | double zobs1, zobs2, hmok, pmbok, rhok, wlok, tol, wlsq, gb, |
---|
906 | a, gamma, tdc, psat, pwo, w, tempo, dn0, rdndr0, sk0, |
---|
907 | f0, rdndrt, zt, ft, dnts, rdndrp, zts, fts, rs, |
---|
908 | dns, rdndrs, zs, fs, refold, z0, zrange, fb, ff, fo, |
---|
909 | fe, h, r, sz, rg, dr, tg, dn, rdndr, t, f, refp, reft; |
---|
910 | |
---|
911 | /* The refraction integrand */ |
---|
912 | #define refi(R,DN,RDNDR) ((RDNDR)/(DN+RDNDR)); |
---|
913 | |
---|
914 | |
---|
915 | |
---|
916 | /* Transform zobs into the normal range */ |
---|
917 | zobs1 = slaDrange ( zobs ); |
---|
918 | zobs2 = fabs ( zobs1 ); |
---|
919 | zobs2 = gmin ( zobs2, d93 ); |
---|
920 | |
---|
921 | /* Keep other arguments within safe bounds */ |
---|
922 | hmok = gmax ( hm, -1000.0 ); |
---|
923 | hmok = gmin ( hmok, 10000.0 ); |
---|
924 | tdkok = gmax ( tdk, 100.0 ); |
---|
925 | tdkok = gmin ( tdkok, 500.0 ); |
---|
926 | pmbok = gmax ( pmb, 0.0 ); |
---|
927 | pmbok = gmin ( pmbok, 10000.0 ); |
---|
928 | rhok = gmax ( rh, 0.0 ); |
---|
929 | rhok = gmin ( rhok, 1.0 ); |
---|
930 | wlok = gmax ( wl, 0.1 ); |
---|
931 | alpha = fabs ( tlr ); |
---|
932 | alpha = gmax ( alpha, 0.001 ); |
---|
933 | alpha = gmin ( alpha, 0.01 ); |
---|
934 | |
---|
935 | /* Tolerance for iteration */ |
---|
936 | w = fabs ( eps ); |
---|
937 | tol = gmin ( w, 0.1 ) / 2.0; |
---|
938 | |
---|
939 | /* Decide whether optical or radio case - switch at 100 micron */ |
---|
940 | optic = ( wlok <= 100.0 ); |
---|
941 | |
---|
942 | /* Set up model atmosphere parameters defined at the observer */ |
---|
943 | wlsq = wlok * wlok; |
---|
944 | gb = 9.784 * ( 1.0 - 0.0026 * cos ( 2.0 * phi ) - 2.8e-7 * hmok ); |
---|
945 | a = ( optic ) ? |
---|
946 | ( ( 287.604 + 1.6288 / wlsq + 0.0136 / ( wlsq * wlsq ) ) |
---|
947 | * 273.15 / 1013.25 ) * 1e-6 |
---|
948 | : |
---|
949 | 77.624e-6; |
---|
950 | gamal = gb * dmd / gcr; |
---|
951 | gamma = gamal / alpha; |
---|
952 | gamm2 = gamma - 2.0; |
---|
953 | delm2 = delta - 2.0; |
---|
954 | tdc = tdkok - 273.15; |
---|
955 | psat = pow ( 10.0, ( 0.7859 + 0.03477 * tdc ) / |
---|
956 | ( 1.0 + 0.00412 * tdc ) ) * |
---|
957 | ( 1.0 + pmbok * ( 4.5e-6 + 6e-10 * tdc * tdc ) ); |
---|
958 | pwo = ( pmbok > 0.0 ) ? |
---|
959 | rhok * psat / ( 1.0 - ( 1.0 - rhok ) * psat / pmbok ) : |
---|
960 | 0.0; |
---|
961 | w = pwo * ( 1.0 - dmw / dmd ) * gamma / ( delta - gamma ); |
---|
962 | c1 = a * ( pmbok + w ) / tdkok; |
---|
963 | c2 = ( a * w + ( optic ? 11.2684e-6 : 12.92e-6 ) * pwo ) / tdkok; |
---|
964 | c3 = ( gamma - 1.0 ) * alpha * c1 / tdkok; |
---|
965 | c4 = ( delta - 1.0 ) * alpha * c2 / tdkok; |
---|
966 | c5 = optic ? 0.0 : 371897e-6 * pwo / tdkok; |
---|
967 | c6 = c5 * delm2 * alpha / ( tdkok * tdkok ); |
---|
968 | |
---|
969 | /* Conditions at the observer */ |
---|
970 | robs = s + hmok; |
---|
971 | atmt ( robs, tdkok, alpha, gamm2, delm2, c1, c2, c3, c4, c5, c6, robs, |
---|
972 | &tempo, &dn0, &rdndr0 ); |
---|
973 | sk0 = dn0 * robs * sin ( zobs2 ); |
---|
974 | f0 = refi ( robs, dn0, rdndr0 ); |
---|
975 | |
---|
976 | /* Conditions at the tropopause in the troposphere */ |
---|
977 | rt = s + ht; |
---|
978 | atmt ( robs, tdkok, alpha, gamm2, delm2, c1, c2, c3, c4, c5, c6, rt, |
---|
979 | &tt, &dnt, &rdndrt ); |
---|
980 | zt = asin ( sk0 / ( rt * dnt ) ); |
---|
981 | ft = refi ( rt, dnt, rdndrt ); |
---|
982 | |
---|
983 | /* Conditions at the tropopause in the stratosphere */ |
---|
984 | atms ( rt, tt, dnt, gamal, rt, &dnts, &rdndrp ); |
---|
985 | zts = asin ( sk0 / ( rt * dnts ) ); |
---|
986 | fts = refi ( rt, dnts, rdndrp ); |
---|
987 | |
---|
988 | /* At the stratosphere limit */ |
---|
989 | rs = s + hs; |
---|
990 | atms ( rt, tt, dnt, gamal, rs, &dns, &rdndrs ); |
---|
991 | zs = asin ( sk0 / ( rs * dns ) ); |
---|
992 | fs = refi ( rs, dns, rdndrs ); |
---|
993 | |
---|
994 | /* |
---|
995 | ** Integrate the refraction integral in two parts; first in the |
---|
996 | ** troposphere (k=1), then in the stratosphere (k=2). |
---|
997 | */ |
---|
998 | |
---|
999 | /* Initialize previous refraction to ensure at least two iterations */ |
---|
1000 | refold = 1e6; |
---|
1001 | |
---|
1002 | /* |
---|
1003 | ** Start off with 8 strips for the troposphere integration, and then |
---|
1004 | ** use the final troposphere value for the stratosphere integration, |
---|
1005 | ** which tends to need more strips. |
---|
1006 | */ |
---|
1007 | is = 8; |
---|
1008 | |
---|
1009 | /* Troposphere then stratosphere */ |
---|
1010 | for ( k = 1; k <= 2; k++ ) { |
---|
1011 | |
---|
1012 | /* Start z, z range, and start and end values */ |
---|
1013 | if ( k == 1 ) { |
---|
1014 | z0 = zobs2; |
---|
1015 | zrange = zt - z0; |
---|
1016 | fb = f0; |
---|
1017 | ff = ft; |
---|
1018 | } else { |
---|
1019 | z0 = zts; |
---|
1020 | zrange = zs - z0; |
---|
1021 | fb = fts; |
---|
1022 | ff = fs; |
---|
1023 | } |
---|
1024 | |
---|
1025 | /* Sums of odd and even values */ |
---|
1026 | fo = 0.0; |
---|
1027 | fe = 0.0; |
---|
1028 | |
---|
1029 | /* First time through loop we have to do every point */ |
---|
1030 | n = 1; |
---|
1031 | |
---|
1032 | /* Start of iteration loop (terminates at specified precision) */ |
---|
1033 | for ( ; ; ) { |
---|
1034 | |
---|
1035 | /* Strip width */ |
---|
1036 | h = zrange / (double) is; |
---|
1037 | |
---|
1038 | /* Initialize distance from Earth centre for quadrature pass */ |
---|
1039 | r = ( k == 1 ) ? robs : rt; |
---|
1040 | |
---|
1041 | /* One pass (no need to compute evens after first time) */ |
---|
1042 | for ( i = 1; i < is; i += n ) { |
---|
1043 | |
---|
1044 | /* Sine of observed zenith distance */ |
---|
1045 | sz = sin ( z0 + h * (double) i ); |
---|
1046 | |
---|
1047 | /* Find r (to nearest metre, maximum four iterations) */ |
---|
1048 | if ( sz > 1e-20 ) { |
---|
1049 | w = sk0 / sz; |
---|
1050 | rg = r; |
---|
1051 | j = 0; |
---|
1052 | do { |
---|
1053 | if ( k == 1 ) { |
---|
1054 | atmt ( robs, tdkok, alpha, gamm2, delm2, |
---|
1055 | c1, c2, c3, c4, c5, c6, rg, |
---|
1056 | &tg, &dn, &rdndr ); |
---|
1057 | } else { |
---|
1058 | atms ( rt, tt, dnt, gamal, rg, &dn, &rdndr ); |
---|
1059 | } |
---|
1060 | dr = ( rg * dn - w ) / ( dn + rdndr ); |
---|
1061 | rg -= dr; |
---|
1062 | } while ( fabs ( dr ) > 1.0 && j++ <= 4 ); |
---|
1063 | r = rg; |
---|
1064 | } |
---|
1065 | |
---|
1066 | /* Find refractive index and integrand at r */ |
---|
1067 | if ( k == 1 ) { |
---|
1068 | atmt ( robs, tdkok, alpha, gamm2, delm2, |
---|
1069 | c1, c2, c3, c4, c5, c6, r, |
---|
1070 | &t, &dn, &rdndr ); |
---|
1071 | } else { |
---|
1072 | atms ( rt, tt, dnt, gamal, r, &dn, &rdndr ); |
---|
1073 | } |
---|
1074 | f = refi ( r, dn, rdndr ); |
---|
1075 | |
---|
1076 | /* Accumulate odd and (first time only) even values */ |
---|
1077 | if ( n == 1 && i % 2 == 0 ) { |
---|
1078 | fe += f; |
---|
1079 | } else { |
---|
1080 | fo += f; |
---|
1081 | } |
---|
1082 | } |
---|
1083 | |
---|
1084 | /* Evaluate the integrand using Simpson's Rule */ |
---|
1085 | refp = h * ( fb + 4.0 * fo + 2.0 * fe + ff ) / 3.0; |
---|
1086 | |
---|
1087 | /* Has the required precision been reached? */ |
---|
1088 | if ( fabs ( refp - refold ) > tol ) { |
---|
1089 | |
---|
1090 | /* No: prepare for next iteration */ |
---|
1091 | refold = refp; /* Save current value for convergence test */ |
---|
1092 | is += is; /* Double the number of strips */ |
---|
1093 | fe += fo; /* Sum of all = sum of evens next time */ |
---|
1094 | fo = 0.0; /* Reset odds accumulator */ |
---|
1095 | n = 2; /* Skip even values next time */ |
---|
1096 | |
---|
1097 | } else { |
---|
1098 | |
---|
1099 | /* Yes: save troposphere component and terminate loop */ |
---|
1100 | if ( k == 1 ) reft = refp; |
---|
1101 | break; |
---|
1102 | } |
---|
1103 | } |
---|
1104 | } |
---|
1105 | |
---|
1106 | /* Result */ |
---|
1107 | *ref = reft + refp; |
---|
1108 | if ( zobs1 < 0.0 ) *ref = - ( *ref ); |
---|
1109 | |
---|
1110 | *ref+=Pidiv2-zobs1; |
---|
1111 | } |
---|
1112 | |
---|
1113 | /*--------------------------------------------------------------------------*/ |
---|
1114 | |
---|
1115 | void Astro::atmt ( double robs, double tdkok, double alpha, double gamm2, |
---|
1116 | double delm2, double c1, double c2, double c3, |
---|
1117 | double c4, double c5, double c6, double r, |
---|
1118 | double *t, double *dn, double *rdndr ) |
---|
1119 | /* |
---|
1120 | ** - - - - - |
---|
1121 | ** a t m t |
---|
1122 | ** - - - - - |
---|
1123 | ** |
---|
1124 | ** Internal routine used by slaRefro: |
---|
1125 | ** |
---|
1126 | ** refractive index and derivative with respect to height for the |
---|
1127 | ** troposphere. |
---|
1128 | ** |
---|
1129 | ** Given: |
---|
1130 | ** robs double height of observer from centre of the Earth (metre) |
---|
1131 | ** tdkok double temperature at the observer (deg K) |
---|
1132 | ** alpha double alpha ) |
---|
1133 | ** gamm2 double gamma minus 2 ) see ES |
---|
1134 | ** delm2 double delta minus 2 ) |
---|
1135 | ** c1 double useful term ) |
---|
1136 | ** c2 double useful term ) |
---|
1137 | ** c3 double useful term ) see source of |
---|
1138 | ** c4 double useful term ) slaRefro main routine |
---|
1139 | ** c5 double useful term ) |
---|
1140 | ** c6 double useful term ) |
---|
1141 | ** r double current distance from the centre of the Earth (metre) |
---|
1142 | ** |
---|
1143 | ** Returned: |
---|
1144 | ** *t double temperature at r (deg K) |
---|
1145 | ** *dn double refractive index at r |
---|
1146 | ** *rdndr double r * rate the refractive index is changing at r |
---|
1147 | ** |
---|
1148 | ** This routine is derived from the ATMOSTRO routine (C.Hohenkerk, |
---|
1149 | ** HMNAO), with enhancements specified by A.T.Sinclair (RGO) to |
---|
1150 | ** handle the radio case. |
---|
1151 | ** |
---|
1152 | ** Note that in the optical case c5 and c6 are zero. |
---|
1153 | */ |
---|
1154 | { |
---|
1155 | double w, tt0, tt0gm2, tt0dm2; |
---|
1156 | |
---|
1157 | w = tdkok - alpha * ( r - robs ); |
---|
1158 | w = gmin ( w, 320.0 ); |
---|
1159 | w = gmax ( w, 200.0 ); |
---|
1160 | tt0 = w / tdkok; |
---|
1161 | tt0gm2 = pow ( tt0, gamm2 ); |
---|
1162 | tt0dm2 = pow ( tt0, delm2 ); |
---|
1163 | *t = w; |
---|
1164 | *dn = 1.0 + ( c1 * tt0gm2 - ( c2 - c5 / w ) * tt0dm2 ) * tt0; |
---|
1165 | *rdndr = r * ( - c3 * tt0gm2 + ( c4 - c6 / tt0 ) * tt0dm2 ); |
---|
1166 | } |
---|
1167 | |
---|
1168 | /*--------------------------------------------------------------------------*/ |
---|
1169 | |
---|
1170 | void Astro::atms ( double rt, double tt, double dnt, double gamal, double r, |
---|
1171 | double *dn, double *rdndr ) |
---|
1172 | /* |
---|
1173 | ** - - - - - |
---|
1174 | ** a t m s |
---|
1175 | ** - - - - - |
---|
1176 | ** |
---|
1177 | ** Internal routine used by slaRefro: |
---|
1178 | ** |
---|
1179 | ** refractive index and derivative with respect to height for the |
---|
1180 | ** stratosphere. |
---|
1181 | ** |
---|
1182 | ** Given: |
---|
1183 | ** rt double height of tropopause from centre of the Earth (metre) |
---|
1184 | ** tt double temperature at the tropopause (deg k) |
---|
1185 | ** dnt double refractive index at the tropopause |
---|
1186 | ** gamal double constant of the atmospheric model = g*md/r |
---|
1187 | ** r double current distance from the centre of the Earth (metre) |
---|
1188 | ** |
---|
1189 | ** Returned: |
---|
1190 | ** *dn double refractive index at r |
---|
1191 | ** *rdndr double r * rate the refractive index is changing at r |
---|
1192 | ** |
---|
1193 | ** This routine is derived from the ATMOSSTR routine (C.Hohenkerk, HMNAO). |
---|
1194 | */ |
---|
1195 | { |
---|
1196 | double b, w; |
---|
1197 | |
---|
1198 | b = gamal / tt; |
---|
1199 | w = ( dnt - 1.0 ) * exp ( - b * ( r - rt ) ); |
---|
1200 | *dn = 1.0 + w; |
---|
1201 | *rdndr = - r * b * w; |
---|
1202 | } |
---|
1203 | |
---|
1204 | |
---|
1205 | |
---|