]> git.uio.no Git - u/mrichter/AliRoot.git/blame - RALICE/AliAstrolab.cxx
Updated calibration object
[u/mrichter/AliRoot.git] / RALICE / AliAstrolab.cxx
CommitLineData
8bde545d 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16// $Id$
17
18///////////////////////////////////////////////////////////////////////////
19// Class AliAstrolab
20// Virtual lab to correlate measurements with astrophysical phenomena.
21//
22// The lab can be given a terrestrial location via the usual longitude
23// and latitude specifications.
24// Since this class is derived from AliTimestamp, a lab can also be
25// given a specific timestamp. Together with the terrestrial location
26// this provides access to local (sidereal) times etc...
27// In addition to the usual astronomical reference frames, a local
28// lab reference frame can also be specified. Together with the lab's
29// timestamp this uniquely defines all the coordinate transformations
30// between the various reference frames.
31// These lab characteristics allow a space-time correlation of lab
32// observations with external (astrophysical) phenomena.
33//
34// Observations are entered as generic signals containing a position,
35// reference frame specification and a timestamp.
36// These observations can then be analysed in various reference frames
37// via the available GET functions.
38//
39// Various external (astrophysical) phenomena may be entered as
40// so-called reference signals.
41// This class provides facilities (e.g. MatchRefSignal) to check
42// correlations of the stored measurement with these reference signals.
43//
44// Coding example :
45// ----------------
46// gSystem->Load("ralice");
47//
48// AliAstrolab lab("IceCube","The South Pole Neutrino Observatory");
49// lab.SetLabPosition(0,-90,"deg"); // South Pole
50//
51// lab.SetLocalFrame(90,180,90,270,0,0); // Local frame has X-axis to the North
52//
53// lab.Data(1,"dms"); // Print laboratory parameters
54//
55// // Enter some observed event to be investigated
56// AliTimestamp ts;
57// ts.SetUT(1989,7,30,8,14,23,738504,0);
58// Float_t vec[3]={1,23.8,118.65};
59// Ali3Vector r;
60// r.SetVector(vec,"sph","deg");
61// lab.SetSignal(&r,"loc","M",&ts,0,"Event10372");
62//
63// // Enter some reference signals
64// Float_t alpha=194818.0;
65// Float_t delta=84400.;
66// lab.SetSignal(alpha,delta,"B",1950,"M",-1,"Altair");
67// alpha=124900.0;
68// delta=272400.;
69// lab.SetSignal(alpha,delta,"B",1950,"M",-1,"NGP");
70// alpha=64508.917;
71// delta=-164258.02;
72// lab.SetSignal(alpha,delta,"J",2000,"M",-1,"Sirius");
73// alpha=23149.08;
74// delta=891550.8;
75// lab.SetSignal(alpha,delta,"J",2000,"M",-1,"Polaris");
76// alpha=43600.;
77// delta=163100.;
78// lab.SetSignal(alpha,delta,"J",2000,"M",-1,"Aldebaran");
79// Float_t l=327.531;
80// Float_t b=-35.8903;
81// Float_t pos[3]={1,90.-b,l};
82// r.SetVector(pos,"sph","deg");
83// lab.SetUT(1989,7,30,8,14,16,0,0);
84// lab->SetSignal(&r,"gal","M",0,-1,"GRB890730");
85//
86// // List all stored objects
87// lab.ListSignals("equ","M",5);
88//
89// // Determine minimal space and time differences with reference signals
90// Double_t da,dt;
91// Int_t ia,it;
92// da=lab.GetDifference(0,"deg",dt,"s",1,&ia,&it);
93// cout << " Minimal differences damin (deg) : " << da << " dtmin (s) : " << dt
94// << " damin-index : " << ia << " dtmin-index : " << it << endl;
95// cout << " damin for "; lab->PrintSignal("equ","T",&ts,5,ia); cout << endl;
96// cout << " dtmin for "; lab->PrintSignal("equ","T",&ts,5,it); cout << endl;
97//
98// // Search for space and time match with the reference signals
99// da=5;
100// dt=10;
101// TArrayI* arr=lab.MatchRefSignal(da,"deg",dt,"s");
102// Int_t index=0;
103// if (arr)
104// {
105// for (Int_t i=0; i<arr->GetSize(); i++)
106// {
107// index=arr->At(i);
108// cout << " Match found for index : " << index << endl;
109// cout << " Corresponding ref. object "; lab->PrintSignal("equ","T",&ts,5,index); cout << endl;
110// }
111// }
112//
113//
114//--- Author: Nick van Eijndhoven 15-mar-2007 Utrecht University
115//- Modified: NvE $Date$ Utrecht University
116///////////////////////////////////////////////////////////////////////////
117
118#include "AliAstrolab.h"
119#include "Riostream.h"
120
121ClassImp(AliAstrolab) // Class implementation to enable ROOT I/O
122
123AliAstrolab::AliAstrolab(const char* name,const char* title) : TNamed(name,title),AliTimestamp()
124{
125// Default constructor
126
127 fToffset=0;
128 fXsig=0;
129 fRefs=0;
130 fBias=0;
131 fGal=0;
132 fIndices=0;
133}
134///////////////////////////////////////////////////////////////////////////
135AliAstrolab::~AliAstrolab()
136{
137// Destructor to delete all allocated memory.
138
139 if (fXsig)
140 {
141 delete fXsig;
142 fXsig=0;
143 }
144 if (fRefs)
145 {
146 delete fRefs;
147 fRefs=0;
148 }
149 if (fIndices)
150 {
151 delete fIndices;
152 fIndices=0;
153 }
154}
155///////////////////////////////////////////////////////////////////////////
156AliAstrolab::AliAstrolab(const AliAstrolab& t) : TNamed(t),AliTimestamp(t)
157{
158// Copy constructor
159
160 fToffset=t.fToffset;
161 fLabPos=t.fLabPos;
162 fXsig=0;
163 if (t.fXsig) fXsig=new AliSignal(*(t.fXsig));
164 fRefs=0;
165 if (t.fRefs)
166 {
167 Int_t size=t.fRefs->GetSize();
168 fRefs=new TObjArray(size);
169 for (Int_t i=0; i<size; i++)
170 {
171 AliSignal* sx=(AliSignal*)t.fRefs->At(i);
172 if (sx) fRefs->AddAt(sx->Clone(),i);
173 }
174 }
175 fBias=0;
176 fGal=0;
177 fIndices=0;
178}
179///////////////////////////////////////////////////////////////////////////
180void AliAstrolab::Data(Int_t mode,TString u)
181{
182// Provide lab information.
183//
184// "mode" indicates the mode of the timestamp info (see AliTimestamp::Date).
185//
186// The string argument "u" allows to choose between different angular units
187// in case e.g. a spherical frame is selected.
188// u = "rad" : angles provided in radians
189// "deg" : angles provided in degrees
190// "dms" : angles provided in ddd:mm:ss.sss
191// "hms" : angles provided in hh:mm:ss.sss
192//
193// The defaults are mode=1 and u="deg".
194
195 const char* name=GetName();
196 const char* title=GetTitle();
197 cout << " *" << ClassName() << "::Data*";
198 if (strlen(name)) cout << " Name : " << GetName();
199 if (strlen(title)) cout << " Title : " << GetTitle();
200 cout << endl;
201
202 Double_t l,b;
203 GetLabPosition(l,b,"deg");
204 cout << " Lab position longitude : "; PrintAngle(l,"deg",u,2);
205 cout << " latitude : "; PrintAngle(b,"deg",u,2);
206 cout << endl;
207 cout << " Lab time offset w.r.t. UT : "; PrintTime(fToffset,12); cout << endl;
208
209 // UT and Local time info
210 Date(mode,fToffset);
211}
212///////////////////////////////////////////////////////////////////////////
213void AliAstrolab::SetLabPosition(Ali3Vector& p)
214{
215// Set the lab position in the terrestrial coordinates.
216// The right handed reference frame is defined such that the North Pole
217// corresponds to a polar angle theta=0 and the Greenwich meridian corresponds
218// to an azimuth angle phi=0, with phi increasing eastwards.
219
220 fLabPos.SetPosition(p);
221
222 // Determine local time offset in fractional hours w.r.t. UT.
223 Double_t vec[3];
224 p.GetVector(vec,"sph","deg");
225 Double_t l=vec[2];
226 fToffset=l/15.;
227}
228///////////////////////////////////////////////////////////////////////////
229void AliAstrolab::SetLabPosition(Double_t l,Double_t b,TString u)
230{
231// Set the lab position in the terrestrial longitude (l) and latitude (b).
232// Positions north of the equator have b>0, whereas b<0 indicates
233// locations south of the equator.
234// Positions east of Greenwich have l>0, whereas l<0 indicates
235// locations west of Greenwich.
236//
237// The string argument "u" allows to choose between different angular units
238// u = "rad" : angles provided in radians
239// "deg" : angles provided in degrees
240// "dms" : angles provided in dddmmss.sss
241// "hms" : angles provided in hhmmss.sss
242//
243// The default is u="deg".
244
245 Double_t r=1,theta=0,phi=0;
246
247 l=ConvertAngle(l,u,"deg");
248 b=ConvertAngle(b,u,"deg");
249
250 Double_t offset=90.;
251
252 theta=offset-b;
253 phi=l;
254
255 Double_t p[3]={r,theta,phi};
256 fLabPos.SetPosition(p,"sph","deg");
257
258 // Local time offset in fractional hours w.r.t. UT.
259 fToffset=l/15.;
260}
261///////////////////////////////////////////////////////////////////////////
262AliPosition AliAstrolab::GetLabPosition() const
263{
264// Provide the lab position in the terrestrial coordinates.
265// The right handed reference frame is defined such that the North Pole
266// corresponds to a polar angle theta=0 and the Greenwich meridian corresponds
267// to an azimuth angle phi=0, with phi increasing eastwards.
268
269 return fLabPos;
270}
271///////////////////////////////////////////////////////////////////////////
272void AliAstrolab::GetLabPosition(Double_t& l,Double_t& b,TString u) const
273{
274// Provide the lab position in the terrestrial longitude (l) and latitude (b).
275// Positions north of the equator have b>0, whereas b<0 indicates
276// locations south of the equator.
277// Positions east of Greenwich have l>0, whereas l<0 indicates
278// locations west of Greenwich.
279//
280// The string argument "u" allows to choose between different angular units
281// u = "rad" : angles provided in radians
282// "deg" : angles provided in degrees
283//
284// The default is u="deg".
285
286 Double_t pi=acos(-1.);
287
288 Double_t offset=90.;
289 if (u=="rad") offset=pi/2.;
290
291 Double_t p[3];
292 fLabPos.GetPosition(p,"sph",u);
293 b=offset-p[1];
294 l=p[2];
295}
296///////////////////////////////////////////////////////////////////////////
297Double_t AliAstrolab::GetLT()
298{
299// Provide the Lab's local time in fractional hours.
300// A mean solar day lasts 24h (i.e. 86400s).
301//
302// In case a hh:mm:ss format is needed, please use the Convert() facility.
303
304 Double_t h=GetLT(fToffset);
305 return h;
306}
307///////////////////////////////////////////////////////////////////////////
308Double_t AliAstrolab::GetLMST()
309{
310// Provide the Lab's Local Mean Sidereal Time (LMST) in fractional hours.
311// A sidereal day corresponds to 23h 56m 04.09s (i.e. 86164.09s) mean solar time.
312// The definition of GMST is such that a sidereal clock corresponds with
313// 24 sidereal hours per revolution of the Earth.
314// As such, local time offsets w.r.t. UT and GMST can be treated similarly.
315//
316// In case a hh:mm:ss format is needed, please use the Convert() facility.
317
318 Double_t h=GetLMST(fToffset);
319 return h;
320}
321///////////////////////////////////////////////////////////////////////////
322Double_t AliAstrolab::GetLAST()
323{
324// Provide the Lab's Local Apparent Sidereal Time (LAST) in fractional hours.
325// A sidereal day corresponds to 23h 56m 04.09s (i.e. 86164.09s) mean solar time.
326// The definition of GMST and GAST is such that a sidereal clock corresponds with
327// 24 sidereal hours per revolution of the Earth.
328// As such, local time offsets w.r.t. UT, GMST and GAST can be treated similarly.
329//
330// In case a hh:mm:ss format is needed, please use the Convert() facility.
331
332 Double_t h=GetLAST(fToffset);
333 return h;
334}
335///////////////////////////////////////////////////////////////////////////
336void AliAstrolab::PrintAngle(Double_t a,TString in,TString out,Int_t ndig) const
337{
338// Printing of angles in various formats.
339//
340// The input argument "a" denotes the angle to be printed.
341// The string arguments "in" and "out" specify the angular I/O formats.
342//
343// in = "rad" : input angle provided in radians
344// "deg" : input angle provided in degrees
345// "dms" : input angle provided in dddmmss.sss
346// "hms" : input angle provided in hhmmss.sss
347//
348// out = "rad" : output angle provided in radians
349// "deg" : output angle provided in degrees
350// "dms" : output angle provided in dddmmss.sss
351// "hms" : output angle provided in hhmmss.sss
352//
353// The argument "ndig" specifies the number of digits for the fractional
354// part (e.g. ndig=6 for "dms" corresponds to micro-arcsecond precision).
355// No rounding will be performed, so an arcsecond count of 3.473 with ndig=1
356// will appear as 03.4 on the output.
357// Due to computer accuracy, precision on the pico-arcsecond level may get lost.
358//
359// The default is ndig=1.
360//
361// Note : The angle info is printed without additional spaces or "endline".
362// This allows the print to be included in various composite output formats.
363
364 Double_t b=ConvertAngle(a,in,out);
365
366 if (out=="deg" || out=="rad")
367 {
368 cout.setf(ios::fixed,ios::floatfield);
369 cout << setprecision(ndig) << b << " " << out.Data();
370 cout.unsetf(ios::fixed);
371 return;
372 }
373
374 Double_t epsilon=1.e-12; // Accuracy in (arc)seconds
375 Int_t word=0,ddd=0,hh=0,mm=0,ss=0;
376 Double_t s;
377 ULong64_t sfrac=0;
378
379 if (out=="dms")
380 {
381 word=Int_t(b);
382 word=abs(word);
383 ddd=word/10000;
384 word=word%10000;
385 mm=word/100;
386 ss=word%100;
387 s=fabs(b)-Double_t(ddd*10000+mm*100+ss);
388 if (s>(1.-epsilon))
389 {
390 s=0.;
391 ss++;
392 }
393 while (ss>=60)
394 {
395 ss-=60;
396 mm++;
397 }
398 while (mm>=60)
399 {
400 mm-=60;
401 ddd++;
402 }
403 while (ddd>=360)
404 {
405 ddd-=360;
406 }
407 s*=pow(10.,ndig);
408 sfrac=ULong64_t(s);
409 if (b<0) cout << "-";
410 cout << ddd << "d " << mm << "' " << ss << "."
411 << setfill('0') << setw(ndig) << sfrac << "\"";
412 return;
413 }
414
415 if (out=="hms")
416 {
417 word=Int_t(b);
418 word=abs(word);
419 hh=word/10000;
420 word=word%10000;
421 mm=word/100;
422 ss=word%100;
423 s=fabs(b)-Double_t(hh*10000+mm*100+ss);
424 if (s>(1.-epsilon))
425 {
426 s=0.;
427 ss++;
428 }
429 while (ss>=60)
430 {
431 ss-=60;
432 mm++;
433 }
434 while (mm>=60)
435 {
436 mm-=60;
437 hh++;
438 }
439 while (hh>=24)
440 {
441 hh-=24;
442 }
443 s*=pow(10.,ndig);
444 sfrac=ULong64_t(s);
445 if (b<0) cout << "-";
446 cout << hh << "h " << mm << "m " << ss << "."
447 << setfill('0') << setw(ndig) << sfrac << "s";
448 return;
449 }
450}
451///////////////////////////////////////////////////////////////////////////
452void AliAstrolab::SetSignal(Ali3Vector* r,TString frame,TString mode,AliTimestamp* ts,Int_t jref,TString name)
453{
454// Store a signal as specified by the position r and the timestamp ts.
455// The position is stored in International Celestial Reference System (ICRS) coordinates.
456// The ICRS is a fixed, time independent frame and as such provides a unique reference
457// frame without the need of specifying any epoch etc...
458// The ICRS coordinate definitions match within 20 mas with the mean ones of the J2000.0
459// equatorial system. Nevertheless, to obtain the highest accuracy, the slight
460// coordinate correction between J2000 and ICRS is performed here via the
461// so-called frame bias matrix.
462// For further details see the U.S. Naval Observatory (USNO) circular 179 (2005),
463// which is available on http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
464//
465// The input parameter "frame" allows the user to specify the frame to which
466// the components of r refer. Available options are :
467//
468// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d),
469// where the "sph" components of r correspond to theta=(pi/2)-d and phi=a.
470// "gal" ==> Galactic coordinates with longitude (l) and lattitude (b).
471// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
472// "ecl" ==> Ecliptic coordinates with longitude (l) and lattitude (b),
473// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
474// "hor" ==> Horizontal coordinates at the AliAstrolab location, where the "sph"
475// components of r correspond to theta=zenith angle and phi=pi-azimuth.
476// "icr" ==> ICRS coordinates with longitude (l) and lattitude (b),
477// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
478// "loc" ==> Local coordinates at the AliAstrolab location, where the "sph"
479// components of r correspond to the usual theta and phi angles.
480//
481// In case the coordinates are the equatorial right ascension and declination (a,d),
482// they can represent so-called "mean" and "true" values.
483// The distinction between these two representations is the following :
484//
485// mean values : (a,d) are only corrected for precession and not for nutation
486// true values : (a,d) are corrected for both precession and nutation
487//
488// The input parameter "mode" allows the user to specifiy either "mean" or "true"
489// values for the input in case of equatorial (a,d) coordinates.
490//
491// mode = "M" --> Input coordinates are the mean values
492// "T" --> Input coordinates are the true values
493//
494// The input parameter "jref" allows the user to store so-called "reference" signals.
495// These reference signals may be used to check space-time event coincidences with the
496// stored measurement (e.g. coincidence of the measurement with transient phenomena).
497//
498// jref = 0 --> Storage of the measurement
499// j --> Storage of a reference signal at the j-th position (j=1 is first)
500// < 0 --> Add a reference signal at the next available position
501//
502// Via the input argument "name" the user can give the stored signal also a name.
503//
504// The default values are jref=0 and name="".
505//
506// Note : In case ts=0 the current timestamp of the lab will be taken.
507
508 if (!r)
509 {
510 if (!jref && fXsig)
511 {
512 delete fXsig;
513 fXsig=0;
514 }
515 return;
516 }
517
518 if (frame!="equ" && frame!="gal" && frame!="ecl" && frame!="hor" && frame!="icr" && frame!="loc")
519 {
520 if (!jref && fXsig)
521 {
522 delete fXsig;
523 fXsig=0;
524 }
525 return;
526 }
527
528 if (frame=="equ" && mode!="M" && mode!="m" && mode!="T" && mode!="t")
529 {
530 if (!jref && fXsig)
531 {
532 delete fXsig;
533 fXsig=0;
534 }
535 return;
536 }
537
538 if (!ts) ts=(AliTimestamp*)this;
539
540 Double_t vec[3]={1,0,0};
541 vec[1]=r->GetX(2,"sph","rad");
542 vec[2]=r->GetX(3,"sph","rad");
543 Ali3Vector q;
544 q.SetVector(vec,"sph","rad");
545
546 AliSignal* sxref=0;
547
548 if (!jref) // Storage of the measurement
549 {
550 if (!fXsig)
551 {
552 fXsig=new AliSignal();
553 }
554 else
555 {
556 fXsig->Reset(1);
557 }
558 if (name != "") fXsig->SetName(name);
559 fXsig->SetTitle("Event in ICRS coordinates");
560 fXsig->SetTimestamp(*ts);
561 }
562 else // Storage of a reference signal
563 {
564 if (!fRefs)
565 {
566 fRefs=new TObjArray();
567 fRefs->SetOwner();
568 }
569 // Expand array size if needed
570 if (jref>0 && jref>=fRefs->GetSize()) fRefs->Expand(jref+1);
571 sxref=new AliSignal();
572 if (name != "") sxref->SetName(name);
573 sxref->SetTitle("Reference event in ICRS coordinates");
574 sxref->SetTimestamp(*ts);
575 }
576
577 if (frame=="loc")
578 {
579 // Convert to horizontal coordinates
580 q=q.GetUnprimed(&fL);
581
582 // Store the signal
583 SetSignal(&q,"hor",mode,ts,jref);
584 return;
585 }
586
587 if (frame=="equ")
588 {
589 // Convert to "mean" values at specified epoch
590 if (mode=="T" || mode=="t")
591 {
592 SetNmatrix(ts);
593 q=q.GetUnprimed(&fN);
594 }
595
596 // Convert to "mean" values at J2000
597 SetPmatrix(ts);
598 q=q.GetUnprimed(&fP);
599
600 // Convert to ICRS values
601 if (!fBias) SetBmatrix();
602 q=q.GetUnprimed(&fB);
603 }
604
605 if (frame=="gal")
606 {
607 // Convert to J2000 equatorial mean coordinates
608 if (fGal != 2) SetGmatrix("J");
609 q=q.GetUnprimed(&fG);
610
611 // Convert to ICRS values
612 if (!fBias) SetBmatrix();
613 q=q.GetUnprimed(&fB);
614 }
615
616 if (frame=="ecl")
617 {
618 // Convert to mean equatorial values at specified epoch
619 SetEmatrix(ts);
620 q=q.GetUnprimed(&fE);
621
622 // Convert to "mean" values at J2000
623 SetPmatrix(ts);
624 q=q.GetUnprimed(&fP);
625
626 // Convert to ICRS values
627 if (!fBias) SetBmatrix();
628 q=q.GetUnprimed(&fB);
629 }
630
631 if (frame=="hor")
632 {
633 // Convert to "true" equatorial values at the specified timestamp
634 SetHmatrix(ts);
635 q=q.GetUnprimed(&fH);
636
637 // Convert to "mean" values at specified timestamp
638 SetNmatrix(ts);
639 q=q.GetUnprimed(&fN);
640
641 // Convert to "mean" values at J2000
642 SetPmatrix(ts);
643 q=q.GetUnprimed(&fP);
644
645 // Convert to ICRS values
646 if (!fBias) SetBmatrix();
647 q=q.GetUnprimed(&fB);
648 }
649
650 // Store the signal in ICRS coordinates
651 if (!jref) // Storage of a regular signal
652 {
653 fXsig->SetPosition(q);
654 }
655 else // Storage of a reference signal
656 {
657 sxref->SetPosition(q);
658 if (jref<0)
659 {
660 fRefs->Add(sxref);
661 }
662 else
663 {
664 fRefs->AddAt(sxref,jref-1);
665 }
666 }
667}
668///////////////////////////////////////////////////////////////////////////
669void AliAstrolab::SetSignal(Double_t a,Double_t d,TString s,Double_t e,TString mode,Int_t jref,TString name)
670{
671// Store a signal with right ascension (a) and declination (d) given for epoch e.
672// The position is stored in International Celestial Reference System (ICRS) coordinates.
673// The ICRS is a fixed, time independent frame and as such provides a unique reference
674// frame without the need of specifying any epoch etc...
675// The ICRS coordinate definitions match within 20 mas the mean ones of the J2000.0
676// equatorial system. Nevertheless, to obtain the highest accuracy, the slight
677// coordinate correction between J2000 and ICRS is performed here via the
678// so-called frame bias matrix.
679// For further details see the U.S. Naval Observatory (USNO) circular 179 (2005),
680// which is available on http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
681//
682// The coordinates (a,d) can represent so-called "mean" and "true" values.
683// The distinction between these two representations is the following :
684//
685// mean values : (a,d) are only corrected for precession and not for nutation
686// true values : (a,d) are corrected for both precession and nutation
687//
688// The input parameter "mode" allows the user to specifiy either "mean" or "true"
689// values for the input (a,d) coordinates.
690//
691// a : Right ascension in hhmmss.sss
692// d : Declination in dddmmss.sss
693// s = "B" --> Besselian reference epoch.
694// "J" --> Julian reference epoch.
695// e : Reference epoch for the input coordinates (e.g. 1900, 1950, 2000,...)
696// mode = "M" --> Input coordinates are the mean values
697// "T" --> Input coordinates are the true values
698//
699// The input parameter "jref" allows the user to store so-called "reference" signals.
700// These reference signals may be used to check space-time event coincidences with the
701// stored measurement (e.g. coincidence of the measurement with transient phenomena).
702//
703// jref = 0 --> Storage of the measurement
704// j --> Storage of a reference signal at the j-th position (j=1 is first)
705// < 0 --> Add a reference signal at the next available position
706//
707// Via the input argument "name" the user can give the stored signal also a name.
708//
709// The default values are jref=0 and name="".
710
711 if (s!="B" && s!="b" && s!="J" && s!="j")
712 {
713 if (!jref && fXsig)
714 {
715 delete fXsig;
716 fXsig=0;
717 }
718 return;
719 }
720
721 if (mode!="M" && mode!="m" && mode!="T" && mode!="t")
722 {
723 if (!jref && fXsig)
724 {
725 delete fXsig;
726 fXsig=0;
727 }
728 return;
729 }
730
731 // Convert coordinates to fractional degrees.
732 a=ConvertAngle(a,"hms","deg");
733 d=ConvertAngle(d,"dms","deg");
734
735
736 AliTimestamp tx;
737 tx.SetEpoch(e,s);
738
739 Ali3Vector r;
740 Double_t vec[3]={1.,90.-d,a};
741 r.SetVector(vec,"sph","deg");
742
743 SetSignal(&r,"equ",mode,&tx,jref,name);
744}
745///////////////////////////////////////////////////////////////////////////
746void AliAstrolab::SetSignal(Double_t a,Double_t d,TString mode,AliTimestamp* ts,Int_t jref,TString name)
747{
748// Store a signal with right ascension (a) and declination (d) given for timestamp ts.
749// The position is stored in International Celestial Reference System (ICRS) coordinates.
750// The ICRS is a fixed, time independent frame and as such provides a unique reference
751// frame without the need of specifying any epoch etc...
752// The ICRS coordinate definitions match within 20 mas the mean ones of the J2000.0
753// equatorial system. Nevertheless, to obtain the highest accuracy, the slight
754// coordinate correction between J2000 and ICRS is performed here via the
755// so-called frame bias matrix.
756// For further details see the U.S. Naval Observatory (USNO) circular 179 (2005),
757// which is available on http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
758//
759// The coordinates (a,d) can represent so-called "mean" and "true" values.
760// The distinction between these two representations is the following :
761//
762// mean values : (a,d) are only corrected for precession and not for nutation
763// true values : (a,d) are corrected for both precession and nutation
764//
765// The input parameter "mode" allows the user to specifiy either "mean" or "true"
766// values for the input (a,d) coordinates.
767//
768// a : Right ascension in hhmmss.sss
769// d : Declination in dddmmss.sss
770// mode = "M" --> Input coordinates are the mean values
771// "T" --> Input coordinates are the true values
772//
773// The input parameter "jref" allows the user to store so-called "reference" signals.
774// These reference signals may be used to check space-time event coincidences with the
775// stored measurement (e.g. coincidence of the measurement with transient phenomena).
776//
777// jref = 0 --> Storage of the measurement
778// j --> Storage of a reference signal at the j-th position (j=1 is first)
779// < 0 --> Add a reference signal at the next available position
780//
781// Via the input argument "name" the user can give the stored signal also a name.
782//
783// The default values are jref=0 and name="".
784//
785// Note : In case ts=0 the current timestamp of the lab will be taken.
786
787 // Convert coordinates to fractional degrees.
788 a=ConvertAngle(a,"hms","deg");
789 d=ConvertAngle(d,"dms","deg");
790
791 Ali3Vector r;
792 Double_t vec[3]={1.,90.-d,a};
793 r.SetVector(vec,"sph","deg");
794
795 SetSignal(&r,"equ",mode,ts,jref,name);
796}
797///////////////////////////////////////////////////////////////////////////
798AliSignal* AliAstrolab::GetSignal(Ali3Vector& r,TString frame,TString mode,AliTimestamp* ts,Int_t jref)
799{
800// Provide the user specified coordinates of a signal at the specific timestamp ts.
801// The coordinates are returned via the vector argument "r".
802// In addition also a pointer to the stored signal object is provided.
803// In case no stored signal was available or one of the input arguments was
804// invalid, the returned pointer will be 0.
805//
806// The input parameter "frame" allows the user to specify the frame to which
807// the components of r refer. Available options are :
808//
809// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d),
810// where the "sph" components of r correspond to theta=(pi/2)-d and phi=a.
811// "gal" ==> Galactic coordinates with longitude (l) and lattitude (b).
812// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
813// "ecl" ==> Ecliptic coordinates with longitude (l) and lattitude (b),
814// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
815// "hor" ==> Horizontal coordinates at the AliAstrolab location, where the "sph"
816// components of r correspond to theta=zenith angle and phi=pi-azimuth.
817// "icr" ==> ICRS coordinates with longitude (l) and lattitude (b),
818// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
819// "loc" ==> Local coordinates at the AliAstrolab location, where the "sph"
820// components of r correspond to the usual theta and phi angles.
821//
822// In case the coordinates are the equatorial right ascension and declination (a,d),
823// they can represent so-called "mean" and "true" values.
824// The distinction between these two representations is the following :
825//
826// mean values : (a,d) are only corrected for precession and not for nutation
827// true values : (a,d) are corrected for both precession and nutation
828//
829// The input parameter "mode" allows the user to specifiy either "mean" or "true"
830// values for the input in case of equatorial (a,d) coordinates.
831//
832// mode = "M" --> Input coordinates are the mean values
833// "T" --> Input coordinates are the true values
834//
835// The input parameter "jref" allows the user to access so-called "reference" signals.
836// These reference signals may be used to check space-time event coincidences with the
837// stored measurement (e.g. coincidence of the measurement with transient phenomena).
838//
839// jref = 0 --> Access to the measurement
840// j --> Access to the reference signal at the j-th position (j=1 is first)
841//
842// Default value is jref=0.
843//
844// Note : In case ts=0 the current timestamp of the lab will be taken.
845
846 r.SetZero();
847
848 if (frame!="equ" && frame!="gal" && frame!="ecl" && frame!="hor" && frame!="icr" && frame!="loc") return 0;
849
850 if (frame=="equ" && mode!="M" && mode!="m" && mode!="T" && mode!="t") return 0;
851
852 AliSignal* sx=GetSignal(jref);
853
854 if (!sx) return 0;
855
856 if (!ts) ts=(AliTimestamp*)this;
857
858 Double_t vec[3];
859 sx->GetPosition(vec,"sph","rad");
860 Ali3Vector q;
861 q.SetVector(vec,"sph","rad");
862
863 if (frame=="icr")
864 {
865 r.Load(q);
866 return sx;
867 }
868
869 // Convert from ICRS to equatorial J2000 coordinates
870 if (!fBias) SetBmatrix();
871 q=q.GetPrimed(&fB);
872
873 if (frame=="equ")
874 {
875 // Precess to specified timestamp
876 AliTimestamp ts1;
877 ts1.SetEpoch(2000,"J");
878 Precess(q,&ts1,ts);
879
880 // Nutation correction if requested
881 if (mode=="T" || mode=="t") Nutate(q,ts);
882 }
883
884 if (frame=="gal")
885 {
886 // Convert from equatorial J2000 to galactic
887 if (fGal != 2) SetGmatrix("J");
888 q=q.GetPrimed(&fG);
889 }
890
891 if (frame=="ecl")
892 {
893 // Precess to specified timestamp
894 AliTimestamp ts1;
895 ts1.SetEpoch(2000,"J");
896 Precess(q,&ts1,ts);
897
898 // Convert from equatorial to ecliptic coordinates
899 SetEmatrix(ts);
900 q=q.GetPrimed(&fE);
901 }
902
903 if (frame=="hor")
904 {
905 // Precess to specified timestamp
906 AliTimestamp ts1;
907 ts1.SetEpoch(2000,"J");
908 Precess(q,&ts1,ts);
909
910 // Nutation correction
911 Nutate(q,ts);
912
913 // Convert from equatorial to horizontal coordinates
914 SetHmatrix(ts);
915 q=q.GetPrimed(&fH);
916 }
917
918 if (frame=="loc")
919 {
920 // Get the signal in horizontal coordinates
921 GetSignal(q,"hor",mode,ts);
922
923 // Convert from horizontal local-frame coordinates
924 q=q.GetPrimed(&fL);
925 }
926
927 r.Load(q);
928 return sx;
929}
930///////////////////////////////////////////////////////////////////////////
931AliSignal* AliAstrolab::GetSignal(Ali3Vector& r,TString frame,TString mode,AliTimestamp* ts,TString name)
932{
933// Provide the user specified coordinates of the signal with the specified
934// name at the specific timestamp ts.
935// The coordinates are returned via the vector argument "r".
936// In addition also a pointer to the stored signal object is provided.
937// In case no such stored signal was available or one of the input arguments was
938// invalid, the returned pointer will be 0.
939//
940// The input parameter "frame" allows the user to specify the frame to which
941// the components of r refer. Available options are :
942//
943// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d),
944// where the "sph" components of r correspond to theta=(pi/2)-d and phi=a.
945// "gal" ==> Galactic coordinates with longitude (l) and lattitude (b).
946// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
947// "ecl" ==> Ecliptic coordinates with longitude (l) and lattitude (b),
948// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
949// "hor" ==> Horizontal coordinates at the AliAstrolab location, where the "sph"
950// components of r correspond to theta=zenith angle and phi=pi-azimuth.
951// "icr" ==> ICRS coordinates with longitude (l) and lattitude (b),
952// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
953// "loc" ==> Local coordinates at the AliAstrolab location, where the "sph"
954// components of r correspond to the usual theta and phi angles.
955//
956// In case the coordinates are the equatorial right ascension and declination (a,d),
957// they can represent so-called "mean" and "true" values.
958// The distinction between these two representations is the following :
959//
960// mean values : (a,d) are only corrected for precession and not for nutation
961// true values : (a,d) are corrected for both precession and nutation
962//
963// The input parameter "mode" allows the user to specifiy either "mean" or "true"
964// values for the input in case of equatorial (a,d) coordinates.
965//
966// mode = "M" --> Input coordinates are the mean values
967// "T" --> Input coordinates are the true values
968//
969// Note : In case ts=0 the current timestamp of the lab will be taken.
970
971 AliSignal* sx=0;
972 Int_t j=GetSignalIndex(name);
973 if (j>=0) sx=GetSignal(r,frame,mode,ts,j);
974 return sx;
975}
976///////////////////////////////////////////////////////////////////////////
977AliSignal* AliAstrolab::GetSignal(Double_t& a,Double_t& d,TString mode,AliTimestamp* ts,Int_t jref)
978{
979// Provide precession (and nutation) corrected right ascension (a) and
980// declination (d) of the stored signal object at the specified timestamp.
981// In addition also a pointer to the stored signal object is provided.
982// In case no stored signal was available or one of the input arguments was
983// invalid, the returned pointer will be 0.
984//
985// The coordinates (a,d) can represent so-called "mean" and "true" values.
986// The distinction between these two representations is the following :
987//
988// mean values : (a,d) are only corrected for precession and not for nutation
989// true values : (a,d) are corrected for both precession and nutation
990//
991// The input parameter "mode" allows the user to select either
992// "mean" or "true" values for (a,d).
993//
994// The correction methods used are the new IAU 2000 ones as described in the
995// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
996// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
997//
998// a : Right ascension in hhmmss.sss
999// d : Declination in dddmmss.sss
1000// mode = "M" --> Output coordinates are the mean values
1001// "T" --> Output coordinates are the true values
1002// ts : Timestamp at which the corrected coordinate values are requested.
1003//
1004// The input parameter "jref" allows the user to access so-called "reference" signals.
1005// These reference signals may be used to check space-time event coincidences with the
1006// stored regular signal (e.g. coincidence of the stored signal with transient phenomena).
1007//
1008// jref = 0 --> Access to the measurement
1009// j --> Access to the reference signal at the j-th position (j=1 is first)
1010//
1011// Default value is jref=0.
1012//
1013// Note : In case ts=0 the current timestamp of the lab will be taken.
1014
1015 a=0;
1016 d=0;
1017
1018 Ali3Vector r;
1019 AliSignal* sx=GetSignal(r,"equ",mode,ts,jref);
1020
1021 if (!sx) return 0;
1022
1023 // Retrieve the requested (a,d) values
1024 Double_t vec[3];
1025 r.GetVector(vec,"sph","deg");
1026 d=90.-vec[1];
1027 a=vec[2];
1028
1029 while (a<-360.)
1030 {
1031 a+=360.;
1032 }
1033 while (a>360.)
1034 {
1035 a-=360.;
1036 }
1037 while (d<-90.)
1038 {
1039 d+=90.;
1040 }
1041 while (d>90.)
1042 {
1043 d-=90.;
1044 }
1045
1046 // Convert coordinates to appropriate format
1047 a=ConvertAngle(a,"deg","hms");
1048 d=ConvertAngle(d,"deg","dms");
1049
1050 return sx;
1051}
1052///////////////////////////////////////////////////////////////////////////
1053AliSignal* AliAstrolab::GetSignal(Double_t& a,Double_t& d,TString mode,AliTimestamp* ts,TString name)
1054{
1055// Provide precession (and nutation) corrected right ascension (a) and
1056// declination (d) of the stored signal object with the specified name
1057// at the specific timestamp ts.
1058// In addition also a pointer to the stored signal object is provided.
1059// In case no stored signal was available or one of the input arguments was
1060// invalid, the returned pointer will be 0.
1061//
1062// The coordinates (a,d) can represent so-called "mean" and "true" values.
1063// The distinction between these two representations is the following :
1064//
1065// mean values : (a,d) are only corrected for precession and not for nutation
1066// true values : (a,d) are corrected for both precession and nutation
1067//
1068// The input parameter "mode" allows the user to select either
1069// "mean" or "true" values for (a,d).
1070//
1071// The correction methods used are the new IAU 2000 ones as described in the
1072// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1073// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1074//
1075// a : Right ascension in hhmmss.sss
1076// d : Declination in dddmmss.sss
1077// mode = "M" --> Output coordinates are the mean values
1078// "T" --> Output coordinates are the true values
1079// ts : Timestamp at which the corrected coordinate values are requested.
1080// name : Name of the requested signal object
1081//
1082// Note : In case ts=0 the current timestamp of the lab will be taken.
1083
1084 AliSignal* sx=0;
1085 Int_t j=GetSignalIndex(name);
1086 if (j>=0) sx=GetSignal(a,d,mode,ts,j);
1087 return sx;
1088}
1089///////////////////////////////////////////////////////////////////////////
1090AliSignal* AliAstrolab::GetSignal(Double_t& a,Double_t& d,TString s,Double_t e,TString mode,Int_t jref)
1091{
1092// Provide precession (and nutation) corrected right ascension (a) and
1093// declination (d) of the stored signal object at the specified epoch e.
1094// In addition also a pointer to the stored signal object is provided.
1095// In case no stored signal was available or one of the input arguments was
1096// invalid, the returned pointer will be 0.
1097//
1098// The coordinates (a,d) can represent so-called "mean" and "true" values.
1099// The distinction between these two representations is the following :
1100//
1101// mean values : (a,d) are only corrected for precession and not for nutation
1102// true values : (a,d) are corrected for both precession and nutation
1103//
1104// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1105// values for the input (a,d) coordinates.
1106//
1107// a : Right ascension in hhmmss.sss
1108// d : Declination in dddmmss.sss
1109// s = "B" --> Besselian reference epoch.
1110// "J" --> Julian reference epoch.
1111// e : Reference epoch for the input coordinates (e.g. 1900, 1950, 2000,...)
1112// mode = "M" --> Input coordinates are the mean values
1113// "T" --> Input coordinates are the true values
1114//
1115// The input parameter "jref" allows the user to access so-called "reference" signals.
1116// These reference signals may be used to check space-time event coincidences with the
1117// stored measurement (e.g. coincidence of the measurement with transient phenomena).
1118//
1119// jref = 0 --> Access to the measurement
1120// j --> Access to the reference signal at the j-th position (j=1 is first)
1121//
1122// Default value is jref=0.
1123
1124 a=0;
1125 d=0;
1126
1127 if (s!="B" && s!="b" && s!="J" && s!="j") return 0;
1128
1129 if (mode!="M" && mode!="m" && mode!="T" && mode!="t") return 0;
1130
1131 // Convert coordinates to fractional degrees.
1132 a=ConvertAngle(a,"hms","deg");
1133 d=ConvertAngle(d,"dms","deg");
1134
1135
1136 AliTimestamp tx;
1137 tx.SetEpoch(e,s);
1138
1139 AliSignal* sx=GetSignal(a,d,mode,&tx,jref);
1140 return sx;
1141}
1142///////////////////////////////////////////////////////////////////////////
1143AliSignal* AliAstrolab::GetSignal(Double_t& a,Double_t& d,TString s,Double_t e,TString mode,TString name)
1144{
1145// Provide precession (and nutation) corrected right ascension (a) and
1146// declination (d) of the stored signal object with the specified name
1147// at the specific epoch e.
1148// In addition also a pointer to the stored signal object is provided.
1149// In case no stored signal was available or one of the input arguments was
1150// invalid, the returned pointer will be 0.
1151//
1152// The coordinates (a,d) can represent so-called "mean" and "true" values.
1153// The distinction between these two representations is the following :
1154//
1155// mean values : (a,d) are only corrected for precession and not for nutation
1156// true values : (a,d) are corrected for both precession and nutation
1157//
1158// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1159// values for the input (a,d) coordinates.
1160//
1161// a : Right ascension in hhmmss.sss
1162// d : Declination in dddmmss.sss
1163// s = "B" --> Besselian reference epoch.
1164// "J" --> Julian reference epoch.
1165// e : Reference epoch for the input coordinates (e.g. 1900, 1950, 2000,...)
1166// mode = "M" --> Input coordinates are the mean values
1167// "T" --> Input coordinates are the true values
1168// name : Name of the requested signal object
1169
1170 AliSignal* sx=0;
1171 Int_t j=GetSignalIndex(name);
1172 if (j>=0) sx=GetSignal(a,d,s,e,mode,j);
1173 return sx;
1174}
1175///////////////////////////////////////////////////////////////////////////
1176AliSignal* AliAstrolab::GetSignal(Int_t jref)
1177{
1178// Provide the pointer to a stored signal object.
1179//
1180// The input parameter "jref" allows the user to access so-called "reference" signals.
1181// These reference signals may be used to check space-time event coincidences with the
1182// stored measurement (e.g. coincidence of the measurement with transient phenomena).
1183//
1184// jref = 0 --> Access to the measurement
1185// j --> Access to the reference signal at the j-th position (j=1 is first)
1186//
1187// Default value is jref=0.
1188
1189 AliSignal* sx=0;
1190 if (!jref)
1191 {
1192 sx=fXsig;
1193 }
1194 else
1195 {
1196 if (jref>0 && jref<fRefs->GetSize()) sx=(AliSignal*)fRefs->At(jref-1);
1197 }
1198 return sx;
1199}
1200///////////////////////////////////////////////////////////////////////////
1201AliSignal* AliAstrolab::GetSignal(TString name)
1202{
1203// Provide the pointer to the stored signal object with the specified name.
1204
1205 AliSignal* sx=0;
1206 Int_t j=GetSignalIndex(name);
1207 if (j>=0) sx=GetSignal(j);
1208 return sx;
1209}
1210///////////////////////////////////////////////////////////////////////////
1211void AliAstrolab::RemoveRefSignal(Int_t j,Int_t compress)
1212{
1213// Remove the reference signal which was stored at the j-th position (j=1 is first).
1214// Note : j=0 means that all stored reference signals will be removed.
1215// j<0 allows array compression (see below) without removing any signals.
1216//
1217// The "compress" parameter allows compression of the ref. signal storage array.
1218//
1219// compress = 1 --> Array will be compressed
1220// 0 --> Array will not be compressed
1221//
1222// Note : Compression of the storage array means that the indices of the
1223// reference signals in the storage array will change.
1224
1225 if (!fRefs) return;
1226
1227 // Clearing of the complete storage
1228 if (!j)
1229 {
1230 delete fRefs;
1231 fRefs=0;
1232 return;
1233 }
1234
1235 // Removing a specific reference signal
1236 if (j>0 && j<fRefs->GetSize())
1237 {
1238 TObject* obj=fRefs->RemoveAt(j-1);
1239 if (obj) delete obj;
1240 }
1241
1242 // Compression of the storage array
1243 if (compress) fRefs->Compress();
1244}
1245///////////////////////////////////////////////////////////////////////////
1246void AliAstrolab::RemoveRefSignal(TString name,Int_t compress)
1247{
1248// Remove the reference signal with the specified name.
1249//
1250// The "compress" parameter allows compression of the ref. signal storage array.
1251//
1252// compress = 1 --> Array will be compressed
1253// 0 --> Array will not be compressed
1254//
1255// Note : Compression of the storage array means that the indices of the
1256// reference signals in the storage array will change.
1257
1258 Int_t j=GetSignalIndex(name);
1259 if (j>0) RemoveRefSignal(j,compress);
1260}
1261///////////////////////////////////////////////////////////////////////////
1262Int_t AliAstrolab::GetSignalIndex(TString name)
1263{
1264// Provide storage index of the signal with the specified name.
1265// In case the name matches with the stored measurement,
1266// the value 0 is returned.
1267// In case no signal with the specified name was found, the value -1 is returned.
1268
1269 Int_t index=-1;
1270
1271 if (fXsig)
1272 {
1273 if (name==fXsig->GetName()) return 0;
1274 }
1275
1276 if (!fRefs) return -1;
1277
1278 for (Int_t i=0; i<fRefs->GetSize(); i++)
1279 {
1280 AliSignal* sx=(AliSignal*)fRefs->At(i);
1281 if (!sx) continue;
1282
1283 if (name==sx->GetName())
1284 {
1285 index=i+1;
1286 break;
1287 }
1288 }
1289
1290 return index;
1291}
1292///////////////////////////////////////////////////////////////////////////
1293void AliAstrolab::PrintSignal(TString frame,TString mode,AliTimestamp* ts,Int_t ndig,Int_t jref)
1294{
1295// Print data of a stored signal in user specified coordinates at the specific timestamp ts.
1296// In case no stored signal was available or one of the input arguments was
1297// invalid, no printout is produced.
1298//
1299// The argument "ndig" specifies the number of digits for the fractional
1300// part (e.g. ndig=6 for "dms" corresponds to micro-arcsecond precision).
1301// No rounding will be performed, so an arcsecond count of 3.473 with ndig=1
1302// will appear as 03.4 on the output.
1303// Due to computer accuracy, precision on the pico-arcsecond level may get lost.
1304//
1305// Note : The angle info is printed without additional spaces or "endline".
1306// This allows the print to be included in various composite output formats.
1307//
1308// The input parameter "frame" allows the user to specify the frame to which
1309// the coordinates refer. Available options are :
1310//
1311// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d).
1312//
1313// "gal" ==> Galactic coordinates with longitude (l) and lattitude (b).
1314//
1315// "ecl" ==> Ecliptic coordinates with longitude (l) and lattitude (b).
1316//
1317// "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
1318//
1319// "icr" ==> ICRS coordinates with longitude (l) and lattitude (b).
1320//
1321// "loc" ==> Local spherical angles theta and phi at the AliAstrolab location.
1322//
1323// In case the coordinates are the equatorial right ascension and declination (a,d),
1324// they can represent so-called "mean" and "true" values.
1325// The distinction between these two representations is the following :
1326//
1327// mean values : (a,d) are only corrected for precession and not for nutation
1328// true values : (a,d) are corrected for both precession and nutation
1329//
1330// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1331// values for the input in case of equatorial (a,d) coordinates.
1332//
1333// mode = "M" --> Input coordinates are the mean values
1334// "T" --> Input coordinates are the true values
1335//
1336// The input parameter "jref" allows printing of a so-called "reference" signal.
1337// These reference signals may serve to check space-time event coincidences with the
1338// stored measurement (e.g. coincidence of the measurement with transient phenomena).
1339//
1340// jref = 0 --> Printing of the measurement
1341// j --> Printing of the j-th reference signal
1342//
1343// Default value is jref=0.
1344//
1345// Note : In case ts=0 the current timestamp of the lab will be taken.
1346
1347 Ali3Vector r;
1348 AliSignal* sx=GetSignal(r,frame,mode,ts,jref);
1349
1350 if (!sx) return;
1351
1352 TString name=sx->GetName();
1353 if (name != "") cout << name.Data() << " ";
1354
1355 if (frame=="equ")
1356 {
1357 Double_t a,d;
1358 d=90.-r.GetX(2,"sph","deg");
1359 a=r.GetX(3,"sph","rad");
1360 cout << "Equatorial (" << mode.Data() <<") a : "; PrintAngle(a,"rad","hms",ndig);
1361 cout << " d : "; PrintAngle(d,"deg","dms",ndig);
1362 return;
1363 }
1364
1365 if (frame=="gal")
1366 {
1367 Double_t l,b;
1368 b=90.-r.GetX(2,"sph","deg");
1369 l=r.GetX(3,"sph","deg");
1370 cout << "Galactic l : "; PrintAngle(l,"deg","deg",ndig);
1371 cout << " b : "; PrintAngle(b,"deg","deg",ndig);
1372 return;
1373 }
1374
1375 if (frame=="icr")
1376 {
1377 Double_t a,d;
1378 d=90.-r.GetX(2,"sph","deg");
1379 a=r.GetX(3,"sph","rad");
1380 cout << "ICRS l : "; PrintAngle(a,"rad","hms",ndig);
1381 cout << " b : "; PrintAngle(d,"deg","dms",ndig);
1382 return;
1383 }
1384
1385 if (frame=="ecl")
1386 {
1387 Double_t a,d;
1388 d=90.-r.GetX(2,"sph","deg");
1389 a=r.GetX(3,"sph","deg");
1390 cout << "Ecliptic l : "; PrintAngle(a,"deg","deg",ndig);
1391 cout << " b : "; PrintAngle(d,"deg","deg",ndig);
1392 return;
1393 }
1394
1395 if (frame=="hor")
1396 {
1397 Double_t alt=90.-r.GetX(2,"sph","deg");
1398 Double_t azi=180.-r.GetX(3,"sph","deg");
1399 while (azi>360)
1400 {
1401 azi-=360.;
1402 }
1403 while (azi<0)
1404 {
1405 azi+=360.;
1406 }
1407 cout << "Horizontal azi : "; PrintAngle(azi,"deg","deg",ndig);
1408 cout << " alt : "; PrintAngle(alt,"deg","deg",ndig);
1409 return;
1410 }
1411
1412 if (frame=="loc")
1413 {
1414 Double_t theta=r.GetX(2,"sph","deg");
1415 Double_t phi=r.GetX(3,"sph","deg");
1416 cout << "Local-frame phi : "; PrintAngle(phi,"deg","deg",ndig);
1417 cout << " theta : "; PrintAngle(theta,"deg","deg",ndig);
1418 return;
1419 }
1420}
1421///////////////////////////////////////////////////////////////////////////
1422void AliAstrolab::PrintSignal(TString frame,TString mode,AliTimestamp* ts,Int_t ndig,TString name)
1423{
1424// Print data of the stored signal with the specified name in user specified coordinates
1425// at the specific timestamp ts.
1426// In case such stored signal was available or one of the input arguments was
1427// invalid, no printout is produced.
1428//
1429// The argument "ndig" specifies the number of digits for the fractional
1430// part (e.g. ndig=6 for "dms" corresponds to micro-arcsecond precision).
1431// No rounding will be performed, so an arcsecond count of 3.473 with ndig=1
1432// will appear as 03.4 on the output.
1433// Due to computer accuracy, precision on the pico-arcsecond level may get lost.
1434//
1435// Note : The angle info is printed without additional spaces or "endline".
1436// This allows the print to be included in various composite output formats.
1437//
1438// The input parameter "frame" allows the user to specify the frame to which
1439// the coordinates refer. Available options are :
1440//
1441// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d).
1442//
1443// "gal" ==> Galactic coordinates with longitude (l) and lattitude (b).
1444//
1445// "ecl" ==> Ecliptic coordinates with longitude (l) and lattitude (b).
1446//
1447// "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
1448//
1449// "icr" ==> ICRS coordinates with longitude (l) and lattitude (b).
1450//
1451// "loc" ==> Local spherical angles theta and phi at the AliAstrolab location.
1452//
1453// In case the coordinates are the equatorial right ascension and declination (a,d),
1454// they can represent so-called "mean" and "true" values.
1455// The distinction between these two representations is the following :
1456//
1457// mean values : (a,d) are only corrected for precession and not for nutation
1458// true values : (a,d) are corrected for both precession and nutation
1459//
1460// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1461// values for the input in case of equatorial (a,d) coordinates.
1462//
1463// mode = "M" --> Input coordinates are the mean values
1464// "T" --> Input coordinates are the true values
1465//
1466// Note : In case ts=0 the current timestamp of the lab will be taken.
1467
1468 Int_t j=GetSignalIndex(name);
1469 if (j>=0) PrintSignal(frame,mode,ts,ndig,j);
1470}
1471///////////////////////////////////////////////////////////////////////////
1472void AliAstrolab::ListSignals(TString frame,TString mode,Int_t ndig)
1473{
1474// List all stored signals in user specified coordinates at the timestamp
1475// of the actual recording of the stored measurement under investigation.
1476// In case no (timestamp of the) actual measurement is available,
1477// the current timestamp of the lab will be taken.
1478// In case no stored signal is available or one of the input arguments is
1479// invalid, no printout is produced.
1480//
1481// The argument "ndig" specifies the number of digits for the fractional
1482// part (e.g. ndig=6 for "dms" corresponds to micro-arcsecond precision).
1483// No rounding will be performed, so an arcsecond count of 3.473 with ndig=1
1484// will appear as 03.4 on the output.
1485// Due to computer accuracy, precision on the pico-arcsecond level may get lost.
1486//
1487// The default value is ndig=1;
1488//
1489// Note : The angle info is printed without additional spaces or "endline".
1490// This allows the print to be included in various composite output formats.
1491//
1492// The input parameter "frame" allows the user to specify the frame to which
1493// the coordinates refer. Available options are :
1494//
1495// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d).
1496//
1497// "gal" ==> Galactic coordinates with longitude (l) and lattitude (b).
1498//
1499// "ecl" ==> Ecliptic coordinates with longitude (l) and lattitude (b).
1500//
1501// "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
1502//
1503// "icr" ==> ICRS coordinates with longitude (l) and lattitude (b).
1504//
1505// "loc" ==> Local spherical angles theta and phi at the AliAstrolab location.
1506//
1507// In case the coordinates are the equatorial right ascension and declination (a,d),
1508// they can represent so-called "mean" and "true" values.
1509// The distinction between these two representations is the following :
1510//
1511// mean values : (a,d) are only corrected for precession and not for nutation
1512// true values : (a,d) are corrected for both precession and nutation
1513//
1514// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1515// values for the input in case of equatorial (a,d) coordinates.
1516//
1517// mode = "M" --> Input coordinates are the mean values
1518// "T" --> Input coordinates are the true values
1519//
1520// Note : In case ts=0 the current timestamp of the lab will be taken.
1521
1522 Int_t iprint=0;
1523
1524 AliTimestamp* tx=0;
1525
1526 if (fXsig)
1527 {
1528 AliTimestamp* tx=fXsig->GetTimestamp();
1529 if (!tx) tx=(AliTimestamp*)this;
1530 cout << " *AliAstrolab::ListSignals* List of all stored signals." << endl;
1531 cout << " The measurement under investigation." << endl;
1532 cout << " --- Timestamp of the actual observation ---" << endl;
1533 cout << " "; tx->Date(1);
1534 cout << " --- Location of the actual observation ---" << endl;
1535 cout << " "; PrintSignal(frame,mode,tx,ndig); cout << endl;
1536 iprint=1;
1537 }
1538
1539 if (!fRefs) return;
1540
1541 for (Int_t i=1; i<=fRefs->GetSize(); i++)
1542 {
1543 AliSignal* sx=GetSignal(i);
1544 if (!sx) continue;
1545
1546 if (!iprint)
1547 {
1548 cout << " *AliAstrolab::ListRefSignals* List of all stored signals." << endl;
1549 tx=(AliTimestamp*)this;
1550 cout << " --- Current timestamp of the laboratory ---" << endl;
1551 cout << " "; tx->Date(1);
1552 iprint=1;
1553 }
1554 if (iprint==1)
1555 {
1556 cout << " --- All stored reference signals according to the above timestamp ---" << endl;
1557 iprint=2;
1558 }
1559 cout << " Index : " << i << " "; PrintSignal(frame,mode,tx,ndig,i); cout << endl;
1560 }
1561}
1562///////////////////////////////////////////////////////////////////////////
1563void AliAstrolab::Precess(Ali3Vector& r,AliTimestamp* ts1,AliTimestamp* ts2)
1564{
1565// Correct mean right ascension and declination given for Julian date "jd"
1566// for the earth's precession corresponding to the specified timestamp.
1567// The results are the so-called "mean" (i.e. precession corrected) values,
1568// corresponding to the specified timestamp.
1569// The method used is the new IAU 2000 one as described in the
1570// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1571// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1572// Since the standard reference epoch is J2000, this implies that all
1573// input (a,d) coordinates will be first internally converted to the
1574// corresponding J2000 values before the precession correction w.r.t. the
1575// specified lab timestamp will be applied.
1576//
1577// r : Input vector containing the right ascension and declination information
1578// in the form of standard Ali3Vector coordinates.
1579// In spherical coordinates the phi corresponds to the right ascension,
1580// whereas the declination corresponds to (pi/2)-theta.
1581// jd : Julian date corresponding to the input coordinate values.
1582// ts : Timestamp corresponding to the requested corrected coordinate values.
1583//
1584// Note : In case ts=0 the current timestamp of the lab will be taken.
1585
1586 // Convert back to J2000 values
1587 Ali3Vector r0;
1588 SetPmatrix(ts1);
1589 r0=r.GetUnprimed(&fP);
1590
1591 // Precess to the specified timestamp
1592 if (!ts2) ts2=(AliTimestamp*)this;
1593 SetPmatrix(ts2);
1594 r=r0.GetPrimed(&fP);
1595}
1596///////////////////////////////////////////////////////////////////////////
1597void AliAstrolab::Nutate(Ali3Vector& r,AliTimestamp* ts)
1598{
1599// Correct mean right ascension and declination for the earth's nutation
1600// corresponding to the specified timestamp.
1601// The results are the so-called "true" (i.e. nutation corrected) values,
1602// corresponding to the specified timestamp.
1603// The method used is the new IAU 2000 one as described in the
1604// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1605// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1606//
1607// r : Input vector containing the right ascension and declination information
1608// in the form of standard Ali3Vector coordinates.
1609// In spherical coordinates the phi corresponds to the right ascension,
1610// whereas the declination corresponds to (pi/2)-theta.
1611// ts : Timestamp for which the corrected coordinate values are requested.
1612//
1613// Note : In case ts=0 the current timestamp of the lab will be taken.
1614
1615 // Nutation correction for the specified timestamp
1616 if (!ts) ts=(AliTimestamp*)this;
1617 SetNmatrix(ts);
1618 r=r.GetPrimed(&fN);
1619}
1620///////////////////////////////////////////////////////////////////////////
1621void AliAstrolab::SetBmatrix()
1622{
1623// Set the frame bias matrix elements.
1624// The formulas and numerical constants used are the ones from the
1625// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1626// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1627
1628 Double_t pi=acos(-1.);
1629
1630 // Parameters in mas
1631 Double_t a=-14.6;
1632 Double_t x=-16.6170;
1633 Double_t e=-6.8192;
1634
1635 // Convert to radians
1636 a*=pi/(180.*3600.*1000.);
1637 x*=pi/(180.*3600.*1000.);
1638 e*=pi/(180.*3600.*1000.);
1639
1640 Double_t mat[9];
1641 mat[0]=1.-0.5*(a*a+x*x);
1642 mat[1]=a;
1643 mat[2]=-x;
1644 mat[3]=-a-e*x;
1645 mat[4]=1.-0.5*(a*a+e*e);
1646 mat[5]=-e;
1647 mat[6]=x-e*a;
1648 mat[7]=e+x*a;
1649 mat[8]=1.-0.5*(e*e+x*x);
1650
1651 fB.SetMatrix(mat);
1652 fBias=1;
1653}
1654///////////////////////////////////////////////////////////////////////////
1655void AliAstrolab::SetPmatrix(AliTimestamp* ts)
1656{
1657// Set precession matrix elements for Julian date jd w.r.t. J2000.
1658// The formulas and numerical constants used are the ones from the
1659// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1660// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1661// All numerical constants refer to the standard reference epoch J2000.
1662
1663 Double_t mat[9]={0,0,0,0,0,0,0,0,0};
1664 if (!ts)
1665 {
1666 fP.SetMatrix(mat);
1667 return;
1668 }
1669
1670 Double_t pi=acos(-1.);
1671
1672 Double_t t=(ts->GetJD()-2451545.0)/36525.; // Julian centuries since J2000.0
1673
1674 // Parameters for the precession matrix in arcseconds
1675 Double_t eps0=84381.406; // Mean ecliptic obliquity at J2000.0
1676 Double_t psi=5038.481507*t-1.0790069*pow(t,2)-0.00114045*pow(t,3)+0.000132851*pow(t,4)
1677 -0.0000000951*pow(t,4);
1678 Double_t om=eps0-0.025754*t+0.0512623*pow(t,2)-0.00772503*pow(t,3)-0.000000467*pow(t,4)
1679 +0.0000003337*pow(t,5);
1680 Double_t chi=10.556403*t-2.3814292*pow(t,2)-0.00121197*pow(t,3)+0.000170663*pow(t,4)
1681 -0.0000000560*pow(t,5);
1682
1683 // Convert to radians
1684 eps0*=pi/(180.*3600.);
1685 psi*=pi/(180.*3600.);
1686 om*=pi/(180.*3600.);
1687 chi*=pi/(180.*3600.);
1688
1689 Double_t s1=sin(eps0);
1690 Double_t s2=sin(-psi);
1691 Double_t s3=sin(-om);
1692 Double_t s4=sin(chi);
1693 Double_t c1=cos(eps0);
1694 Double_t c2=cos(-psi);
1695 Double_t c3=cos(-om);
1696 Double_t c4=cos(chi);
1697
1698 mat[0]=c4*c2-s2*s4*c3;
1699 mat[1]=c4*s2*c1+s4*c3*c2*c1-s1*s4*s3;
1700 mat[2]=c4*s2*s1+s4*c3*c2*s1+c1*s4*s3;
1701 mat[3]=-s4*c2-s2*c4*c3;
1702 mat[4]=-s4*s2*c1+c4*c3*c2*c1-s1*c4*s3;
1703 mat[5]=-s4*s2*s1+c4*c3*c2*s1+c1*c4*s3;
1704 mat[6]=s2*s3;
1705 mat[7]=-s3*c2*c1-s1*c3;
1706 mat[8]=-s3*c2*s1+c3*c1;
1707
1708 fP.SetMatrix(mat);
1709}
1710///////////////////////////////////////////////////////////////////////////
1711void AliAstrolab::SetNmatrix(AliTimestamp* ts)
1712{
1713// Set nutation matrix elements for the specified Julian date jd.
1714// The formulas and numerical constants used are the ones from the
1715// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1716// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1717
1718 Double_t mat[9]={0,0,0,0,0,0,0,0,0};
1719 if (!ts)
1720 {
1721 fN.SetMatrix(mat);
1722 return;
1723 }
1724
1725 Double_t pi=acos(-1.);
1726
1727 Double_t dpsi,deps,eps;
1728 ts->Almanac(&dpsi,&deps,&eps);
1729
1730 // Convert to radians
1731 dpsi*=pi/(180.*3600.);
1732 deps*=pi/(180.*3600.);
1733 eps*=pi/(180.*3600.);
1734
1735 Double_t s1=sin(eps);
1736 Double_t s2=sin(-dpsi);
1737 Double_t s3=sin(-(eps+deps));
1738 Double_t c1=cos(eps);
1739 Double_t c2=cos(-dpsi);
1740 Double_t c3=cos(-(eps+deps));
1741
1742 mat[0]=c2;
1743 mat[1]=s2*c1;
1744 mat[2]=s2*s1;
1745 mat[3]=-s2*c3;
1746 mat[4]=c3*c2*c1-s1*s3;
1747 mat[5]=c3*c2*s1+c1*s3;
1748 mat[6]=s2*s3;
1749 mat[7]=-s3*c2*c1-s1*c3;
1750 mat[8]=-s3*c2*s1+c3*c1;
1751
1752 fN.SetMatrix(mat);
1753}
1754///////////////////////////////////////////////////////////////////////////
1755void AliAstrolab::SetGmatrix(TString mode)
1756{
1757// Set the mean equatorial to galactic coordinate conversion matrix.
1758// The B1950 parameters were taken from section 22.3 of the book
1759// "An Introduction to Modern Astrophysics" by Carrol and Ostlie (1996).
1760// The J2000 parameters are obtained by precession of the B1950 values.
1761//
1762// Via the input argument "mode" the required epoch can be selected
1763// mode = "B" ==> B1950
1764// "J" ==> J2000
1765
1766 Ali3Vector x; // The Galactic x-axis in the equatorial frame
1767 Ali3Vector y; // The Galactic y-axis in the equatorial frame
1768 Ali3Vector z; // The Galactic z-axis in the equatorial frame
1769
1770 Double_t a,d;
1771 Double_t vec[3]={1,0,0};
1772
1773 fGal=1; // Set flag to indicate B1950 matrix values
1774
1775 // B1950 equatorial coordinates of the North Galactic Pole (NGP)
1776 a=124900.;
1777 d=272400.;
1778 a=ConvertAngle(a,"hms","deg");
1779 d=ConvertAngle(d,"dms","deg");
1780 vec[1]=90.-d;
1781 vec[2]=a;
1782 z.SetVector(vec,"sph","deg");
1783
1784 // B1950 equatorial coordinates of the Galactic l=b=0 point
1785 a=174224.;
1786 d=-285500.;
1787 a=ConvertAngle(a,"hms","deg");
1788 d=ConvertAngle(d,"dms","deg");
1789 vec[1]=90.-d;
1790 vec[2]=a;
1791 x.SetVector(vec,"sph","deg");
1792
1793 // Precess to the corresponding J2000 values if requested
1794 if (mode=="J")
1795 {
1796 fGal=2; // Set flag to indicate J2000 matrix values
1797 AliTimestamp t1;
1798 t1.SetEpoch(1950,"B");
1799 AliTimestamp t2;
1800 t2.SetEpoch(2000,"J");
1801 Precess(z,&t1,&t2);
1802 Precess(x,&t1,&t2);
1803 }
1804
1805 // The Galactic y-axis is determined for the right handed frame
1806 y=z.Cross(x);
1807
1808 fG.SetAngles(x.GetX(2,"sph","deg"),x.GetX(3,"sph","deg"),
1809 y.GetX(2,"sph","deg"),y.GetX(3,"sph","deg"),
1810 z.GetX(2,"sph","deg"),z.GetX(3,"sph","deg"));
1811}
1812///////////////////////////////////////////////////////////////////////////
1813void AliAstrolab::SetEmatrix(AliTimestamp* ts)
1814{
1815// Set the mean equatorial to ecliptic coordinate conversion matrix
1816// for the specified timestamp.
1817// A nice sketch and explanation of the two frames can be found
1818// in chapter 3 of the book "Astronomy Methods" by Hale Bradt (2004).
1819
1820 Double_t dpsi,deps,eps;
1821 ts->Almanac(&dpsi,&deps,&eps);
1822
1823 // Convert to degrees
1824 eps/=3600.;
1825
1826 // Positions of the ecliptic axes w.r.t. the equatorial ones
1827 // at the moment of the specified timestamp
1828 Double_t theta1=90; // Ecliptic x-axis
1829 Double_t phi1=0;
1830 Double_t theta2=90.-eps; //Ecliptic y-axis
1831 Double_t phi2=90;
1832 Double_t theta3=eps; // Ecliptic z-axis
1833 Double_t phi3=270;
1834
1835 fE.SetAngles(theta1,phi1,theta2,phi2,theta3,phi3);
1836}
1837///////////////////////////////////////////////////////////////////////////
1838void AliAstrolab::SetHmatrix(AliTimestamp* ts)
1839{
1840// Set the mean equatorial to horizontal coordinate conversion matrix
1841// for the specified timestamp.
1842// A nice sketch and explanation of the two frames can be found
1843// in chapter 3 of the book "Astronomy Methods" by Hale Bradt (2004).
1844//
1845// Note : In order to simplify the calculations, we use here a
1846// right-handed horizontal frame.
1847
1848 Ali3Vector x; // The (South pointing) horizontal x-axis in the equatorial frame
1849 Ali3Vector y; // The (East pointing) horizontal y-axis in the equatorial frame
1850 Ali3Vector z; // The (Zenith pointing) horizontal z-axis in the equatorial frame
1851
1852 Double_t l,b;
1853 GetLabPosition(l,b,"deg");
1854
1855 Double_t a;
1856 Double_t vec[3]={1,0,0};
1857
1858 // Equatorial coordinates of the horizontal z-axis
1859 // at the moment of the specified timestamp
1860 a=ts->GetLAST(fToffset);
1861 a*=15.; // Convert fractional hours to degrees
1862 vec[1]=90.-b;
1863 vec[2]=a;
1864 z.SetVector(vec,"sph","deg");
1865
1866 // Equatorial coordinates of the horizontal x-axis
1867 // at the moment of the specified timestamp
1868 vec[1]=180.-b;
1869 vec[2]=a;
1870 x.SetVector(vec,"sph","deg");
1871
1872 // The horizontal y-axis is determined for the right handed frame
1873 y=z.Cross(x);
1874
1875 fH.SetAngles(x.GetX(2,"sph","deg"),x.GetX(3,"sph","deg"),
1876 y.GetX(2,"sph","deg"),y.GetX(3,"sph","deg"),
1877 z.GetX(2,"sph","deg"),z.GetX(3,"sph","deg"));
1878}
1879///////////////////////////////////////////////////////////////////////////
1880void AliAstrolab::SetLocalFrame(Double_t t1,Double_t p1,Double_t t2,Double_t p2,Double_t t3,Double_t p3)
1881{
1882// Specification of the orientations of the local-frame axes.
1883// The input arguments represent the angles (in degrees) of the local-frame axes
1884// w.r.t. a so called Master Reference Frame (MRF), with the same convention
1885// as the input arguments of TRrotMatix::SetAngles.
1886//
1887// The right handed Master Reference Frame is defined as follows :
1888// Z-axis : Points to the local Zenith
1889// X-axis : Makes an angle of 90 degrees with the Z-axis and points South
1890// Y-axis : Makes an angle of 90 degrees with the Z-axis and points East
1891//
1892// Once the user has specified the local reference frame, any observed event
1893// can be related to astronomical space-time locations via the SetSignal
1894// and GetSignal memberfunctions.
1895
1896 // Set the matrix for the conversion of our reference frame coordinates
1897 // into the local-frame ones.
1898
1899 fL.SetAngles(t1,p1,t2,p2,t3,p3);
1900}
1901///////////////////////////////////////////////////////////////////////////
1902Double_t AliAstrolab::ConvertAngle(Double_t a,TString in,TString out) const
1903{
1904// Conversion of various angular formats.
1905//
1906// The input argument "a" denotes the angle to be converted.
1907// The string arguments "in" and "out" specify the angular I/O formats.
1908//
1909// in = "rad" : input angle provided in radians
1910// "deg" : input angle provided in degrees
1911// "dms" : input angle provided in dddmmss.sss
1912// "hms" : input angle provided in hhmmss.sss
1913//
1914// out = "rad" : output angle provided in radians
1915// "deg" : output angle provided in degrees
1916// "dms" : output angle provided in dddmmss.sss
1917// "hms" : output angle provided in hhmmss.sss
1918
1919 if (in==out) return a;
1920
1921 // Convert input to its absolute value in (fractional) degrees.
1922 Double_t pi=acos(-1.);
1923 Double_t epsilon=1.e-12; // Accuracy in (arc)seconds
1924 Int_t word=0,ddd=0,hh=0,mm=0,ss=0;
1925 Double_t s=0;
1926
1927 Double_t b=fabs(a);
1928
1929 if (in=="rad") b*=180./pi;
1930
1931 if (in=="dms")
1932 {
1933 word=Int_t(b);
1934 ddd=word/10000;
1935 word=word%10000;
1936 mm=word/100;
1937 ss=word%100;
1938 s=b-Double_t(ddd*10000+mm*100+ss);
1939 b=Double_t(ddd)+Double_t(mm)/60.+(Double_t(ss)+s)/3600.;
1940 }
1941
1942 if (in=="hms")
1943 {
1944 word=Int_t(b);
1945 hh=word/10000;
1946 word=word%10000;
1947 mm=word/100;
1948 ss=word%100;
1949 s=b-Double_t(hh*10000+mm*100+ss);
1950 b=15.*(Double_t(hh)+Double_t(mm)/60.+(Double_t(ss)+s)/3600.);
1951 }
1952
1953 while (b>360)
1954 {
1955 b-=360.;
1956 }
1957
1958 if (out=="rad") b*=pi/180.;
1959
1960 if (out=="dms")
1961 {
1962 ddd=Int_t(b);
1963 b=b-Double_t(ddd);
1964 b*=60.;
1965 mm=Int_t(b);
1966 b=b-Double_t(mm);
1967 b*=60.;
1968 ss=Int_t(b);
1969 s=b-Double_t(ss);
1970 if (s>(1.-epsilon))
1971 {
1972 s=0.;
1973 ss++;
1974 }
1975 while (ss>=60)
1976 {
1977 ss-=60;
1978 mm++;
1979 }
1980 while (mm>=60)
1981 {
1982 mm-=60;
1983 ddd++;
1984 }
1985 while (ddd>=360)
1986 {
1987 ddd-=360;
1988 }
1989 b=Double_t(10000*ddd+100*mm+ss)+s;
1990 }
1991
1992 if (out=="hms")
1993 {
1994 b/=15.;
1995 hh=Int_t(b);
1996 b=b-Double_t(hh);
1997 b*=60.;
1998 mm=Int_t(b);
1999 b=b-Double_t(mm);
2000 b*=60.;
2001 ss=Int_t(b);
2002 s=b-Double_t(ss);
2003 if (s>(1.-epsilon))
2004 {
2005 s=0.;
2006 ss++;
2007 }
2008 while (ss>=60)
2009 {
2010 ss-=60;
2011 mm++;
2012 }
2013 while (mm>=60)
2014 {
2015 mm-=60;
2016 hh++;
2017 }
2018 while (hh>=24)
2019 {
2020 hh-=24;
2021 }
2022 b=Double_t(10000*hh+100*mm+ss)+s;
2023 }
2024
2025 if (a<0) b=-b;
2026
2027 return b;
2028}
2029///////////////////////////////////////////////////////////////////////////
2030Double_t AliAstrolab::GetHourAngle(TString mode,AliTimestamp* ts,Int_t jref)
2031{
2032// Provide the Local Hour Angle (in fractional degrees) of a stored signal
2033// object at the specified timestamp.
2034//
2035// The input parameter "mode" allows the user to select either the
2036// "mean" or "apparent" value for the returned Hour Angle.
2037//
2038// mode = "M" --> Output is the Mean Hour Angle
2039// "A" --> Output is the Apparent Hour Angle
2040// ts : Timestamp at which the hour angle is requested.
2041//
2042// The input parameter "jref" allows the user to specify so-called "reference" signals.
2043// These reference signals may be used to check space-time event coincidences with the
2044// stored regular signal (e.g. coincidence of the stored signal with transient phenomena).
2045//
2046// jref = 0 --> Use the regular signal
2047// j --> Use the reference signal at the j-th position (j=1 is first)
2048//
2049// Default value is jref=0.
2050//
2051// Note : In case ts=0 the current timestamp of the lab will be taken.
2052
2053 if (!ts) ts=(AliTimestamp*)this;
2054
2055 // Get corrected right ascension and declination for the specified timestamp.
2056 Double_t a,d;
2057 GetSignal(a,d,"T",ts,jref);
2058
2059 // Convert coordinates to fractional degrees.
2060 a=ConvertAngle(a,"hms","deg");
2061 d=ConvertAngle(d,"dms","deg");
2062
2063 a/=15.; // Convert a to fractional hours
2064 Double_t ha=0;
2065 if (mode=="M" || mode=="m") ha=ts->GetLMST(fToffset)-a;
2066 if (mode=="A" || mode=="a") ha=ts->GetLAST(fToffset)-a;
2067 ha*=15.; // Convert to (fractional) degrees
2068
2069 return ha;
2070}
2071///////////////////////////////////////////////////////////////////////////
2072void AliAstrolab::SetLT(Int_t y,Int_t m,Int_t d,Int_t hh,Int_t mm,Int_t ss,Int_t ns,Int_t ps)
2073{
2074// Set the AliTimestamp parameters corresponding to the LT date and time
2075// in the Gregorian calendar as specified by the input arguments.
2076//
2077// The input arguments represent the following :
2078// y : year in LT (e.g. 1952, 2003 etc...)
2079// m : month in LT (1=jan 2=feb etc...)
2080// d : day in LT (1-31)
2081// hh : elapsed hours in LT (0-23)
2082// mm : elapsed minutes in LT (0-59)
2083// ss : elapsed seconds in LT (0-59)
2084// ns : remaining fractional elapsed second of LT in nanosecond
2085// ps : remaining fractional elapsed nanosecond of LT in picosecond
2086//
2087// Note : ns=0 and ps=0 are the default values.
2088//
2089
2090 SetLT(fToffset,y,m,d,hh,mm,ss,ns,ps);
2091}
2092///////////////////////////////////////////////////////////////////////////
2093void AliAstrolab::SetLT(Int_t y,Int_t d,Int_t s,Int_t ns,Int_t ps)
2094{
2095// Set the AliTimestamp parameters corresponding to the specified elapsed
2096// timespan since the beginning of the new LT year.
2097//
2098// The LT year and elapsed time span is entered via the following input arguments :
2099//
2100// y : year in LT (e.g. 1952, 2003 etc...)
2101// d : elapsed number of days
2102// s : (remaining) elapsed number of seconds
2103// ns : (remaining) elapsed number of nanoseconds
2104// ps : (remaining) elapsed number of picoseconds
2105//
2106// The specified d, s, ns and ps values will be used in an additive
2107// way to determine the elapsed timespan.
2108// So, specification of d=1, s=100, ns=0, ps=0 will result in the
2109// same elapsed time span as d=0, s=24*3600+100, ns=0, ps=0.
2110// However, by making use of the latter the user should take care
2111// of possible integer overflow problems in the input arguments,
2112// which obviously will provide incorrect results.
2113//
2114// Note : ns=0 and ps=0 are the default values.
2115
2116 SetLT(fToffset,y,d,s,ns,ps);
2117}
2118///////////////////////////////////////////////////////////////////////////
2119Double_t AliAstrolab::GetDifference(Int_t j,TString au,Double_t& dt,TString tu,Int_t mode,Int_t* ia,Int_t* it)
2120{
2121// Provide space and time difference between the j-th reference signal
2122// (j=1 indicates first) and the stored measurement.
2123//
2124// The return value of this memberfunction provides the positional angular
2125// difference, whereas the output argument "dt" provides the time difference.
2126//
2127// The units of the angular difference can be specified via the the "au"
2128// input argument, where
2129//
2130// au = "rad" --> Angular difference in (fractional) radians
2131// "deg" --> Angular difference in (fractional) degrees
2132//
2133// The units of the time difference can be specified via the "tu" and "mode"
2134// input arguments. For details please refer to AliTimestamp::GetDifference().
2135// Also here mode=1 is the default value.
2136//
2137// For the time difference the reference signal is used as the standard.
2138// This means that in case of a positive time difference, the stored
2139// measurement occurred later than the reference signal.
2140//
2141// In case j=0, the stored measurement will be compared with each
2142// reference signal and the returned angular and time differences are
2143// the minimal differences which were encountered.
2144// In this case the user may obtain the indices of the two stored reference signals
2145// which had the minimal angular and minimal time difference via the output
2146// arguments "ia" and "it" as follows :
2147//
2148// ia = Index of the stored reference signal with minimial angular difference
2149// it = Index of the stored reference signal with minimial time difference
2150//
2151// In case these indices are the same, there obviously was 1 single reference signal
2152// which showed both the minimal angular and time difference.
2153//
2154// The default values are mode=1, ia=0 and it=0;
2155
2156 Double_t dang=999;
2157 dt=1.e20;
2158 if (ia) *ia=0;
2159 if (it) *it=0;
2160
2161 if (!fRefs) return dang;
2162
2163 Ali3Vector rx; // Position of the measurement
2164 Ali3Vector r0; // Position of the reference signal
2165
2166 AliSignal* sx=GetSignal(rx,"icr","M",0);
2167 if (!sx) return dang;
2168
2169 AliTimestamp* tx=sx->GetTimestamp();
2170 if (!tx) return dang;
2171
2172 // Space and time difference w.r.t. a specific reference signal
2173 if (j>0)
2174 {
2175 AliSignal* s0=GetSignal(r0,"icr","M",0,j);
2176 if (!s0) return dang;
2177 AliTimestamp* t0=s0->GetTimestamp();
2178
2179 if (!t0) return dang;
2180
2181 dang=r0.GetOpeningAngle(rx,au);
2182 dt=t0->GetDifference(tx,tu,mode);
2183 return dang;
2184 }
2185
2186 // Minimal space and time difference encountered over all reference signals
2187 Double_t dangmin=dang;
2188 Double_t dtmin=dt;
2189 for (Int_t i=1; i<=fRefs->GetSize(); i++)
2190 {
2191 AliSignal* s0=GetSignal(r0,"icr","M",0,i);
2192 if (!s0) continue;
2193 AliTimestamp* t0=s0->GetTimestamp();
2194 if (!t0) continue;
2195 dang=r0.GetOpeningAngle(rx,au);
2196 dt=t0->GetDifference(tx,tu,mode);
2197 if (fabs(dang)<dangmin)
2198 {
2199 dangmin=fabs(dang);
2200 *ia=i;
2201 }
2202 if (fabs(dt)<dtmin)
2203 {
2204 dtmin=fabs(dt);
2205 *it=i;
2206 }
2207 }
2208
2209 dt=dtmin;
2210 return dangmin;
2211}
2212///////////////////////////////////////////////////////////////////////////
2213Double_t AliAstrolab::GetDifference(TString name,TString au,Double_t& dt,TString tu,Int_t mode)
2214{
2215// Provide space and time difference between the reference signal
2216// with the specified name and the stored measurement.
2217//
2218// The return value of this memberfunction provides the positional angular
2219// difference, whereas the output argument "dt" provides the time difference.
2220//
2221// The units of the angular difference can be specified via the the "au"
2222// input argument, where
2223//
2224// au = "rad" --> Angular difference in (fractional) radians
2225// "deg" --> Angular difference in (fractional) degrees
2226//
2227// The units of the time difference can be specified via the "tu" and "mode"
2228// input arguments. For details please refer to AliTimestamp::GetDifference().
2229// Also here mode=1 is the default value.
2230//
2231// For the time difference the reference signal is used as the standard.
2232// This means that in case of a positive time difference, the stored
2233// measurement occurred later than the reference signal.
2234
2235 Double_t dang=999;
2236 dt=1.e20;
2237
2238 Int_t j=GetSignalIndex(name);
2239 if (j>0) dang=GetDifference(j,au,dt,tu,mode);
2240 return dang;
2241}
2242///////////////////////////////////////////////////////////////////////////
2243TArrayI* AliAstrolab::MatchRefSignal(Double_t da,TString au,Double_t dt,TString tu,Int_t mode)
2244{
2245// Provide the storage indices of the reference signals which match in space
2246// and time with the stored measurement.
2247// The indices are returned via a pointer to a TArrayI object.
2248// In case no matches were found, the null pointer is returned.
2249// A reference signal is regarded as matching with the stored measurement
2250// if the positional angular difference doesn't exceed "da" and the absolute
2251// value of the time difference doesn't exceed "dt".
2252//
2253// The units of the angular difference "da" can be specified via the the "au"
2254// input argument, where
2255//
2256// au = "rad" --> Angular difference in (fractional) radians
2257// "deg" --> Angular difference in (fractional) degrees
2258//
2259// The units of the time difference "dt" can be specified via the "tu" and "mode"
2260// input arguments. For details please refer to AliTimestamp::GetDifference().
2261// Also here mode=1 is the default value.
2262
2263 if (!fXsig || !fRefs) return 0;
2264
2265 Int_t nrefs=fRefs->GetEntries();
2266
2267 if (fIndices) delete fIndices;
2268 fIndices=new TArrayI(nrefs);
2269
2270 Double_t dang,dtime;
2271 Int_t jfill=0;
2272 for (Int_t i=1; i<=fRefs->GetSize(); i++)
2273 {
2274 dang=GetDifference(i,au,dtime,tu,mode);
2275 if (fabs(dang)<=da && fabs(dtime)<=dt)
2276 {
2277 fIndices->AddAt(i,jfill);
2278 jfill++;
2279 }
2280 }
2281
2282 fIndices->Set(jfill);
2283 return fIndices;
2284}
2285///////////////////////////////////////////////////////////////////////////
2286TObject* AliAstrolab::Clone(const char* name) const
2287{
2288// Make a deep copy of the current object and provide the pointer to the copy.
2289// This memberfunction enables automatic creation of new objects of the
2290// correct type depending on the object type, a feature which may be very useful
2291// for containers when adding objects in case the container owns the objects.
2292
2293 AliAstrolab* lab=new AliAstrolab(*this);
2294 if (name)
2295 {
2296 if (strlen(name)) lab->SetName(name);
2297 }
2298 return lab;
2299}
2300///////////////////////////////////////////////////////////////////////////