source: tags/Mars-V0.9.4.1/mastro/MAstro.cc

Last change on this file was 7227, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 18.5 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz, 11/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MAstro
28// ------
29//
30////////////////////////////////////////////////////////////////////////////
31#include "MAstro.h"
32
33#include <iostream>
34
35#include <TVector3.h> // TVector3
36
37#include "MTime.h" // MTime::GetGmst
38#include "MString.h"
39
40#include "MAstroCatalog.h" // FIXME: replace by MVector3!
41
42using namespace std;
43
44ClassImp(MAstro);
45
46Double_t MAstro::Trunc(Double_t val)
47{
48 // dint(A) - truncate to nearest whole number towards zero (double)
49 return val<0 ? TMath::Ceil(val) : TMath::Floor(val);
50}
51
52Double_t MAstro::Round(Double_t val)
53{
54 // dnint(A) - round to nearest whole number (double)
55 return val<0 ? TMath::Ceil(val-0.5) : TMath::Floor(val+0.5);
56}
57
58Double_t MAstro::Hms2Sec(Int_t deg, UInt_t min, Double_t sec, Char_t sgn)
59{
60 const Double_t rc = TMath::Sign((60.0 * (60.0 * (Double_t)TMath::Abs(deg) + (Double_t)min) + sec), (Double_t)deg);
61 return sgn=='-' ? -rc : rc;
62}
63
64Double_t MAstro::Dms2Rad(Int_t deg, UInt_t min, Double_t sec, Char_t sgn)
65{
66 // pi/(180*3600): arcseconds to radians
67 //#define DAS2R 4.8481368110953599358991410235794797595635330237270e-6
68 return Hms2Sec(deg, min, sec, sgn)*TMath::Pi()/(180*3600)/**DAS2R*/;
69}
70
71Double_t MAstro::Hms2Rad(Int_t hor, UInt_t min, Double_t sec, Char_t sgn)
72{
73 // pi/(12*3600): seconds of time to radians
74//#define DS2R 7.2722052166430399038487115353692196393452995355905e-5
75 return Hms2Sec(hor, min, sec, sgn)*TMath::Pi()/(12*3600)/**DS2R*/;
76}
77
78Double_t MAstro::Dms2Deg(Int_t deg, UInt_t min, Double_t sec, Char_t sgn)
79{
80 return Hms2Sec(deg, min, sec, sgn)/3600.;
81}
82
83Double_t MAstro::Hms2Deg(Int_t hor, UInt_t min, Double_t sec, Char_t sgn)
84{
85 return Hms2Sec(hor, min, sec, sgn)/240.;
86}
87
88Double_t MAstro::Dms2Hor(Int_t deg, UInt_t min, Double_t sec, Char_t sgn)
89{
90 return Hms2Sec(deg, min, sec, sgn)/54000.;
91}
92
93Double_t MAstro::Hms2Hor(Int_t hor, UInt_t min, Double_t sec, Char_t sgn)
94{
95 return Hms2Sec(hor, min, sec, sgn)/3600.;
96}
97
98void MAstro::Day2Hms(Double_t day, Char_t &sgn, UShort_t &hor, UShort_t &min, UShort_t &sec)
99{
100 /* Handle sign */
101 sgn = day<0?'-':'+';
102
103 /* Round interval and express in smallest units required */
104 Double_t a = Round(86400. * TMath::Abs(day)); // Days to seconds
105
106 /* Separate into fields */
107 const Double_t ah = Trunc(a/3600.);
108 a -= ah * 3600.;
109 const Double_t am = Trunc(a/60.);
110 a -= am * 60.;
111 const Double_t as = Trunc(a);
112
113 /* Return results */
114 hor = (UShort_t)ah;
115 min = (UShort_t)am;
116 sec = (UShort_t)as;
117}
118
119void MAstro::Rad2Hms(Double_t rad, Char_t &sgn, UShort_t &deg, UShort_t &min, UShort_t &sec)
120{
121 Day2Hms(rad/(TMath::Pi()*2), sgn, deg, min, sec);
122}
123
124void MAstro::Rad2Dms(Double_t rad, Char_t &sgn, UShort_t &deg, UShort_t &min, UShort_t &sec)
125{
126 Rad2Hms(rad*15, sgn, deg, min, sec);
127}
128
129void MAstro::Deg2Dms(Double_t d, Char_t &sgn, UShort_t &deg, UShort_t &min, UShort_t &sec)
130{
131 Day2Hms(d/24, sgn, deg, min, sec);
132}
133
134void MAstro::Deg2Hms(Double_t d, Char_t &sgn, UShort_t &deg, UShort_t &min, UShort_t &sec)
135{
136 Day2Hms(d/360, sgn, deg, min, sec);
137}
138
139void MAstro::Hor2Dms(Double_t h, Char_t &sgn, UShort_t &deg, UShort_t &min, UShort_t &sec)
140{
141 Day2Hms(h*15/24, sgn, deg, min, sec);
142}
143
144void MAstro::Hor2Hms(Double_t h, Char_t &sgn, UShort_t &deg, UShort_t &min, UShort_t &sec)
145{
146 Day2Hms(h/24, sgn, deg, min, sec);
147}
148
149void MAstro::Day2Hm(Double_t day, Char_t &sgn, UShort_t &hor, Double_t &min)
150{
151 /* Handle sign */
152 sgn = day<0?'-':'+';
153
154 /* Round interval and express in smallest units required */
155 Double_t a = Round(86400. * TMath::Abs(day)); // Days to seconds
156
157 /* Separate into fields */
158 const Double_t ah = Trunc(a/3600.);
159 a -= ah * 3600.;
160
161 /* Return results */
162 hor = (UShort_t)ah;
163 min = a/60.;
164}
165
166void MAstro::Rad2Hm(Double_t rad, Char_t &sgn, UShort_t &deg, Double_t &min)
167{
168 Day2Hm(rad/(TMath::Pi()*2), sgn, deg, min);
169}
170
171void MAstro::Rad2Dm(Double_t rad, Char_t &sgn, UShort_t &deg, Double_t &min)
172{
173 Rad2Hm(rad*15, sgn, deg, min);
174}
175
176void MAstro::Deg2Dm(Double_t d, Char_t &sgn, UShort_t &deg, Double_t &min)
177{
178 Day2Hm(d/24, sgn, deg, min);
179}
180
181void MAstro::Deg2Hm(Double_t d, Char_t &sgn, UShort_t &deg, Double_t &min)
182{
183 Rad2Hm(d/360, sgn, deg, min);
184}
185
186void MAstro::Hor2Dm(Double_t h, Char_t &sgn, UShort_t &deg, Double_t &min)
187{
188 Day2Hm(h*15/24, sgn, deg, min);
189}
190
191void MAstro::Hor2Hm(Double_t h, Char_t &sgn, UShort_t &deg, Double_t &min)
192{
193 Day2Hm(h/24, sgn, deg, min);
194}
195
196TString MAstro::GetStringDeg(Double_t deg, const char *fmt)
197{
198 Char_t sgn;
199 UShort_t d, m, s;
200 Deg2Dms(deg, sgn, d, m, s);
201
202 MString str;
203 str.Print(fmt, sgn, d, m ,s);
204 return str;
205}
206
207TString MAstro::GetStringHor(Double_t deg, const char *fmt)
208{
209 Char_t sgn;
210 UShort_t h, m, s;
211 Hor2Hms(deg, sgn, h, m, s);
212
213 MString str;
214 str.Print(fmt, sgn, h, m ,s);
215 return str;
216}
217
218// --------------------------------------------------------------------------
219//
220// Interpretes a string ' - 12 30 00.0' or '+ 12 30 00.0'
221// as floating point value -12.5 or 12.5. If interpretation is
222// successfull kTRUE is returned, otherwise kFALSE. ret is not
223// touched if interpretation was not successfull. The successfull
224// interpreted part is removed from the TString.
225//
226Bool_t MAstro::String2Angle(TString &str, Double_t &ret)
227{
228 Char_t sgn;
229 Int_t d, len;
230 UInt_t m;
231 Float_t s;
232
233 // Skip whitespaces before %c and after %f
234 int n=sscanf(str.Data(), " %c %d %d %f %n", &sgn, &d, &m, &s, &len);
235
236 if (n!=4 || (sgn!='+' && sgn!='-'))
237 return kFALSE;
238
239 str.Remove(0, len);
240
241 ret = Dms2Deg(d, m, s, sgn);
242 return kTRUE;
243}
244
245// --------------------------------------------------------------------------
246//
247// Interpretes a string '-12:30:00.0', '12:30:00.0' or '+12:30:00.0'
248// as floating point value -12.5, 12.5 or 12.5. If interpretation is
249// successfull kTRUE is returned, otherwise kFALSE. ret is not
250// touched if interpretation was not successfull.
251//
252Bool_t MAstro::Coordinate2Angle(const TString &str, Double_t &ret)
253{
254 Char_t sgn = str[0]=='-' ? '-' : '+';
255 Int_t d;
256 UInt_t m;
257 Float_t s;
258
259 const int n=sscanf(str[0]=='+'||str[0]=='-' ? str.Data()+1 : str.Data(), "%d:%d:%f", &d, &m, &s);
260
261 if (n!=3)
262 return kFALSE;
263
264 ret = Dms2Deg(d, m, s, sgn);
265 return kTRUE;
266}
267
268// --------------------------------------------------------------------------
269//
270// Returns val=-12.5 as string '-12:30:00'
271//
272TString MAstro::Angle2Coordinate(Double_t val)
273{
274 Char_t sgn;
275 UShort_t d,m,s;
276
277 Deg2Dms(val, sgn, d, m, s);
278
279 return Form("%c%02d:%02d:%02d", sgn, d, m, s);
280}
281
282// --------------------------------------------------------------------------
283//
284// Return year y, month m and day d corresponding to Mjd.
285//
286void MAstro::Mjd2Ymd(UInt_t mjd, UShort_t &y, Byte_t &m, Byte_t &d)
287{
288 // Express day in Gregorian calendar
289 const ULong_t jd = mjd + 2400001;
290 const ULong_t n4 = 4*(jd+((6*((4*jd-17918)/146097))/4+1)/2-37);
291 const ULong_t nd10 = 10*(((n4-237)%1461)/4)+5;
292
293 y = n4/1461L-4712;
294 m = ((nd10/306+2)%12)+1;
295 d = (nd10%306)/10+1;
296}
297
298// --------------------------------------------------------------------------
299//
300// Return Mjd corresponding to year y, month m and day d.
301//
302Int_t MAstro::Ymd2Mjd(UShort_t y, Byte_t m, Byte_t d)
303{
304 // Month lengths in days
305 static int months[12] = { 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 };
306
307 // Validate month
308 if (m<1 || m>12)
309 return -1;
310
311 // Allow for leap year
312 months[1] = (y%4==0 && (y%100!=0 || y%400==0)) ? 29 : 28;
313
314 // Validate day
315 if (d<1 || d>months[m-1])
316 return -1;
317
318 // Precalculate some values
319 const Byte_t lm = 12-m;
320 const ULong_t lm10 = 4712 + y - lm/10;
321
322 // Perform the conversion
323 return 1461L*lm10/4 + (306*((m+9)%12)+5)/10 - (3*((lm10+188)/100))/4 + d - 2399904;
324}
325
326// --------------------------------------------------------------------------
327//
328// theta0, phi0 [rad]: polar angle/zenith distance, azimuth of 1st object
329// theta1, phi1 [rad]: polar angle/zenith distance, azimuth of 2nd object
330// AngularDistance [rad]: Angular distance between two objects
331//
332Double_t MAstro::AngularDistance(Double_t theta0, Double_t phi0, Double_t theta1, Double_t phi1)
333{
334 TVector3 v0(1);
335 v0.Rotate(phi0, theta0);
336
337 TVector3 v1(1);
338 v1.Rotate(phi1, theta1);
339
340 return v0.Angle(v1);
341}
342
343// --------------------------------------------------------------------------
344//
345// Calls MTime::GetGmst() Better use MTime::GetGmst() directly
346//
347Double_t MAstro::UT2GMST(Double_t ut1)
348{
349 return MTime(ut1).GetGmst();
350}
351
352// --------------------------------------------------------------------------
353//
354// RotationAngle
355//
356// calculates the angle for the rotation of the sky coordinate system
357// with respect to the local coordinate system. This is identical
358// to the rotation angle of the sky image in the camera.
359//
360// sinl [rad]: sine of observers latitude
361// cosl [rad]: cosine of observers latitude
362// theta [rad]: polar angle/zenith distance
363// phi [rad]: rotation angle/azimuth
364//
365// Return sin/cos component of angle
366//
367// The convention is such, that the rotation angle is -pi/pi if
368// right ascension and local rotation angle are counted in the
369// same direction, 0 if counted in the opposite direction.
370//
371// (In other words: The rotation angle is 0 when the source culminates)
372//
373// Using vectors it can be done like:
374// TVector3 v, p;
375// v.SetMagThetaPhi(1, theta, phi);
376// p.SetMagThetaPhi(1, TMath::Pi()/2-latitude, 0);
377// v = v.Cross(l));
378// v.RotateZ(-phi);
379// v.Rotate(-theta)
380// rho = TMath::ATan2(v(2), v(1));
381//
382// For more information see TDAS 00-11, eqs. (18) and (20)
383//
384void MAstro::RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi, Double_t &sin, Double_t &cos)
385{
386 const Double_t sint = TMath::Sin(theta);
387 const Double_t cost = TMath::Cos(theta);
388
389 const Double_t snlt = sinl*sint;
390 const Double_t cslt = cosl*cost;
391
392 const Double_t sinp = TMath::Sin(phi);
393 const Double_t cosp = TMath::Cos(phi);
394
395 const Double_t v1 = sint*sinp;
396 const Double_t v2 = cslt - snlt*cosp;
397
398 const Double_t denom = TMath::Sqrt(v1*v1 + v2*v2);
399
400 sin = cosl*sinp / denom; // y-component
401 cos = (snlt-cslt*cosp) / denom; // x-component
402}
403
404// --------------------------------------------------------------------------
405//
406// RotationAngle
407//
408// calculates the angle for the rotation of the sky coordinate system
409// with respect to the local coordinate system. This is identical
410// to the rotation angle of the sky image in the camera.
411//
412// sinl [rad]: sine of observers latitude
413// cosl [rad]: cosine of observers latitude
414// theta [rad]: polar angle/zenith distance
415// phi [rad]: rotation angle/azimuth
416//
417// Return angle [rad] in the range -pi, pi
418//
419// The convention is such, that the rotation angle is -pi/pi if
420// right ascension and local rotation angle are counted in the
421// same direction, 0 if counted in the opposite direction.
422//
423// (In other words: The rotation angle is 0 when the source culminates)
424//
425// Using vectors it can be done like:
426// TVector3 v, p;
427// v.SetMagThetaPhi(1, theta, phi);
428// p.SetMagThetaPhi(1, TMath::Pi()/2-latitude, 0);
429// v = v.Cross(l));
430// v.RotateZ(-phi);
431// v.Rotate(-theta)
432// rho = TMath::ATan2(v(2), v(1));
433//
434// For more information see TDAS 00-11, eqs. (18) and (20)
435//
436Double_t MAstro::RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi)
437{
438 const Double_t snlt = sinl*TMath::Sin(theta);
439 const Double_t cslt = cosl*TMath::Cos(theta);
440
441 const Double_t sinp = TMath::Sin(phi);
442 const Double_t cosp = TMath::Cos(phi);
443
444 return TMath::ATan2(cosl*sinp, snlt-cslt*cosp);
445}
446
447
448// --------------------------------------------------------------------------
449//
450// Kepler - solve the equation of Kepler
451//
452Double_t MAstro::Kepler(Double_t m, Double_t ecc)
453{
454 m *= TMath::DegToRad();
455
456 Double_t delta = 0;
457 Double_t e = m;
458 do {
459 delta = e - ecc * sin(e) - m;
460 e -= delta / (1 - ecc * cos(e));
461 } while (fabs(delta) > 1e-6);
462
463 return e;
464}
465
466// --------------------------------------------------------------------------
467//
468// GetMoonPhase - calculate phase of moon as a fraction:
469// Returns -1 if calculation failed
470//
471Double_t MAstro::GetMoonPhase(Double_t mjd)
472{
473 /****** Calculation of the Sun's position. ******/
474
475 // date within epoch
476 const Double_t epoch = 44238; // 1980 January 0.0
477 const Double_t day = mjd - epoch;
478 if (day<0)
479 {
480 cout << "MAstro::GetMoonPhase - Day before Jan 1980" << endl;
481 return -1;
482 }
483
484 // mean anomaly of the Sun
485 const Double_t n = fmod(day*360/365.2422, 360);
486
487 const Double_t elonge = 278.833540; // ecliptic longitude of the Sun at epoch 1980.0
488 const Double_t elongp = 282.596403; // ecliptic longitude of the Sun at perigee
489
490 // convert from perigee co-ordinates to epoch 1980.0
491 const Double_t m = fmod(n + elonge - elongp + 360, 360);
492
493 // solve equation of Kepler
494 const Double_t eccent = 0.016718; // eccentricity of Earth's orbit
495 const Double_t k = Kepler(m, eccent);
496 const Double_t ec0 = sqrt((1 + eccent) / (1 - eccent)) * tan(k / 2);
497 // true anomaly
498 const Double_t ec = 2 * atan(ec0) * TMath::RadToDeg();
499
500 // Sun's geocentric ecliptic longitude
501 const Double_t lambdasun = fmod(ec + elongp + 720, 360);
502
503
504 /****** Calculation of the Moon's position. ******/
505
506 // Moon's mean longitude.
507 const Double_t mmlong = 64.975464; // moon's mean lonigitude at the epoch
508 const Double_t ml = fmod(13.1763966*day + mmlong + 360, 360);
509 // Moon's mean anomaly.
510 const Double_t mmlongp = 349.383063; // mean longitude of the perigee at the epoch
511 const Double_t mm = fmod(ml - 0.1114041*day - mmlongp + 720, 360);
512 // Evection.
513 const Double_t ev = 1.2739 * sin((2 * (ml - lambdasun) - mm)*TMath::DegToRad());
514 // Annual equation.
515 const Double_t sinm = TMath::Sin(m*TMath::DegToRad());
516 const Double_t ae = 0.1858 * sinm;
517 // Correction term.
518 const Double_t a3 = 0.37 * sinm;
519 // Corrected anomaly.
520 const Double_t mmp = (mm + ev - ae - a3)*TMath::DegToRad();
521 // Correction for the equation of the centre.
522 const Double_t mec = 6.2886 * sin(mmp);
523 // Another correction term.
524 const Double_t a4 = 0.214 * sin(2 * mmp);
525 // Corrected longitude.
526 const Double_t lp = ml + ev + mec - ae + a4;
527 // Variation.
528 const Double_t v = 0.6583 * sin(2 * (lp - lambdasun)*TMath::DegToRad());
529 // True longitude.
530 const Double_t lpp = lp + v;
531 // Age of the Moon in degrees.
532 const Double_t age = (lpp - lambdasun)*TMath::DegToRad();
533
534 // Calculation of the phase of the Moon.
535 return (1 - TMath::Cos(age)) / 2;
536}
537
538// --------------------------------------------------------------------------
539//
540// Calculate the Period to which the time belongs to. The Period is defined
541// as the number of synodic months ellapsed since the first full moon
542// after Jan 1st 1980 (which was @ MJD=44240.37917)
543//
544Double_t MAstro::GetMoonPeriod(Double_t mjd)
545{
546 const Double_t synmonth = 29.53058868; // synodic month (new Moon to new Moon)
547 const Double_t epoch0 = 44240.37917; // First full moon after 1980/1/1
548
549 const Double_t et = mjd-epoch0; // Ellapsed time
550 return et/synmonth;
551}
552
553// --------------------------------------------------------------------------
554//
555// To get the moon period as defined for MAGIC observation we take the
556// nearest integer mjd, eg:
557// 53257.8 --> 53258
558// 53258.3 --> 53258
559// Which is the time between 13h and 12:59h of the following day. To
560// this day-period we assign the moon-period at midnight. To get
561// the MAGIC definition we now substract 284.
562//
563// For MAGIC observation period do eg:
564// GetMagicPeriod(53257.91042)
565// or
566// MTime t;
567// t.SetMjd(53257.91042);
568// GetMagicPeriod(t.GetMjd());
569// or
570// MTime t;
571// t.Set(2004, 1, 1, 12, 32, 11);
572// GetMagicPeriod(t.GetMjd());
573//
574Int_t MAstro::GetMagicPeriod(Double_t mjd)
575{
576 const Double_t mmjd = (Double_t)TMath::Nint(mjd);
577 const Double_t period = GetMoonPeriod(mmjd);
578
579 return (Int_t)TMath::Floor(period)-284;
580}
581
582// --------------------------------------------------------------------------
583//
584// Returns the distance in x,y between two polar-vectors (eg. Alt/Az, Ra/Dec)
585// projected on aplain in a distance dist. For Magic this this the distance
586// of the camera plain (1700mm) dist also determins the unit in which
587// the TVector2 is returned.
588//
589// v0 is the reference vector (eg. the vector to the center of the camera)
590// v1 is the vector to which we determin the distance on the plain
591//
592// (see also MStarCamTrans::Loc0LocToCam())
593//
594TVector2 MAstro::GetDistOnPlain(const TVector3 &v0, TVector3 v1, Double_t dist)
595{
596 v1.RotateZ(-v0.Phi());
597 v1.RotateY(-v0.Theta());
598 v1.RotateZ(-TMath::Pi()/2); // exchange x and y
599 v1 *= dist/v1.Z();
600
601 return v1.XYvector(); //TVector2(v1.Y(), -v1.X());//v1.XYvector();
602}
603
604// --------------------------------------------------------------------------
605//
606// Calculate the absolute misspointing from the nominal zenith angle nomzd
607// and the deviations in zd (devzd) and az (devaz).
608// All values given in deg, the return value, too.
609//
610Double_t MAstro::GetDevAbs(Double_t nomzd, Double_t devzd, Double_t devaz)
611{
612 const Double_t pzd = nomzd * TMath::DegToRad();
613 const Double_t azd = devzd * TMath::DegToRad();
614 const Double_t aaz = devaz * TMath::DegToRad();
615
616 const double el = TMath::Pi()/2-pzd;
617
618 const double dphi2 = aaz/2.;
619 const double cos2 = TMath::Cos(dphi2)*TMath::Cos(dphi2);
620 const double sin2 = TMath::Sin(dphi2)*TMath::Sin(dphi2);
621 const double d = TMath::Cos(azd)*cos2 - TMath::Cos(2*el)*sin2;
622
623 return TMath::ACos(d)*TMath::RadToDeg();
624}
Note: See TracBrowser for help on using the repository browser.