]> git.uio.no Git - u/mrichter/AliRoot.git/blame - RALICE/AliAstrolab.cxx
Modifying logging output of components.
[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//
661385a7 1336// The input parameter "mode" also determines which type of time and
1337// local hour angle will appear in the printout.
1338//
1339// mode = "M" --> Mean Siderial Time (MST) and Local Mean Hour Angle (LMHA)
1340// "T" --> Apparent Siderial Time (AST) and Local Apparent Hour Angle (LAHA)
1341//
8bde545d 1342// The input parameter "jref" allows printing of a so-called "reference" signal.
1343// These reference signals may serve to check space-time event coincidences with the
1344// stored measurement (e.g. coincidence of the measurement with transient phenomena).
1345//
1346// jref = 0 --> Printing of the measurement
1347// j --> Printing of the j-th reference signal
1348//
1349// Default value is jref=0.
1350//
1351// Note : In case ts=0 the current timestamp of the lab will be taken.
1352
1353 Ali3Vector r;
1354 AliSignal* sx=GetSignal(r,frame,mode,ts,jref);
1355
1356 if (!sx) return;
1357
661385a7 1358 // Local Hour Angle of the signal
1359 Double_t lha=GetHourAngle("A",ts,jref);
1360 TString slha="LAHA";
1361 if (mode=="M" || mode=="m")
1362 {
1363 lha=GetHourAngle("M",ts,jref);
1364 slha="LMHA";
1365 }
1366
8bde545d 1367 TString name=sx->GetName();
1368 if (name != "") cout << name.Data() << " ";
1369
1370 if (frame=="equ")
1371 {
1372 Double_t a,d;
1373 d=90.-r.GetX(2,"sph","deg");
1374 a=r.GetX(3,"sph","rad");
1375 cout << "Equatorial (" << mode.Data() <<") a : "; PrintAngle(a,"rad","hms",ndig);
1376 cout << " d : "; PrintAngle(d,"deg","dms",ndig);
661385a7 1377 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1378 return;
1379 }
1380
1381 if (frame=="gal")
1382 {
1383 Double_t l,b;
1384 b=90.-r.GetX(2,"sph","deg");
1385 l=r.GetX(3,"sph","deg");
1386 cout << "Galactic l : "; PrintAngle(l,"deg","deg",ndig);
1387 cout << " b : "; PrintAngle(b,"deg","deg",ndig);
661385a7 1388 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1389 return;
1390 }
1391
1392 if (frame=="icr")
1393 {
1394 Double_t a,d;
1395 d=90.-r.GetX(2,"sph","deg");
1396 a=r.GetX(3,"sph","rad");
1397 cout << "ICRS l : "; PrintAngle(a,"rad","hms",ndig);
1398 cout << " b : "; PrintAngle(d,"deg","dms",ndig);
661385a7 1399 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1400 return;
1401 }
1402
1403 if (frame=="ecl")
1404 {
1405 Double_t a,d;
1406 d=90.-r.GetX(2,"sph","deg");
1407 a=r.GetX(3,"sph","deg");
1408 cout << "Ecliptic l : "; PrintAngle(a,"deg","deg",ndig);
1409 cout << " b : "; PrintAngle(d,"deg","deg",ndig);
661385a7 1410 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1411 return;
1412 }
1413
1414 if (frame=="hor")
1415 {
1416 Double_t alt=90.-r.GetX(2,"sph","deg");
1417 Double_t azi=180.-r.GetX(3,"sph","deg");
1418 while (azi>360)
1419 {
1420 azi-=360.;
1421 }
1422 while (azi<0)
1423 {
1424 azi+=360.;
1425 }
1426 cout << "Horizontal azi : "; PrintAngle(azi,"deg","deg",ndig);
1427 cout << " alt : "; PrintAngle(alt,"deg","deg",ndig);
661385a7 1428 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1429 return;
1430 }
1431
1432 if (frame=="loc")
1433 {
1434 Double_t theta=r.GetX(2,"sph","deg");
1435 Double_t phi=r.GetX(3,"sph","deg");
1436 cout << "Local-frame phi : "; PrintAngle(phi,"deg","deg",ndig);
1437 cout << " theta : "; PrintAngle(theta,"deg","deg",ndig);
661385a7 1438 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1439 return;
1440 }
1441}
1442///////////////////////////////////////////////////////////////////////////
1443void AliAstrolab::PrintSignal(TString frame,TString mode,AliTimestamp* ts,Int_t ndig,TString name)
1444{
1445// Print data of the stored signal with the specified name in user specified coordinates
1446// at the specific timestamp ts.
1447// In case such stored signal was available or one of the input arguments was
1448// invalid, no printout is produced.
1449//
1450// The argument "ndig" specifies the number of digits for the fractional
1451// part (e.g. ndig=6 for "dms" corresponds to micro-arcsecond precision).
1452// No rounding will be performed, so an arcsecond count of 3.473 with ndig=1
1453// will appear as 03.4 on the output.
1454// Due to computer accuracy, precision on the pico-arcsecond level may get lost.
1455//
1456// Note : The angle info is printed without additional spaces or "endline".
1457// This allows the print to be included in various composite output formats.
1458//
1459// The input parameter "frame" allows the user to specify the frame to which
1460// the coordinates refer. Available options are :
1461//
1462// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d).
1463//
1464// "gal" ==> Galactic coordinates with longitude (l) and lattitude (b).
1465//
1466// "ecl" ==> Ecliptic coordinates with longitude (l) and lattitude (b).
1467//
1468// "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
1469//
1470// "icr" ==> ICRS coordinates with longitude (l) and lattitude (b).
1471//
1472// "loc" ==> Local spherical angles theta and phi at the AliAstrolab location.
1473//
1474// In case the coordinates are the equatorial right ascension and declination (a,d),
1475// they can represent so-called "mean" and "true" values.
1476// The distinction between these two representations is the following :
1477//
1478// mean values : (a,d) are only corrected for precession and not for nutation
1479// true values : (a,d) are corrected for both precession and nutation
1480//
1481// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1482// values for the input in case of equatorial (a,d) coordinates.
1483//
1484// mode = "M" --> Input coordinates are the mean values
1485// "T" --> Input coordinates are the true values
1486//
661385a7 1487// The input parameter "mode" also determines which type of time and
1488// local hour angle will appear in the printout.
1489//
1490// mode = "M" --> Mean Siderial Time (MST) and Local Mean Hour Angle (LMHA)
1491// "T" --> Apparent Siderial Time (AST) and Local Apparent Hour Angle (LAHA)
1492//
8bde545d 1493// Note : In case ts=0 the current timestamp of the lab will be taken.
1494
1495 Int_t j=GetSignalIndex(name);
1496 if (j>=0) PrintSignal(frame,mode,ts,ndig,j);
1497}
1498///////////////////////////////////////////////////////////////////////////
1499void AliAstrolab::ListSignals(TString frame,TString mode,Int_t ndig)
1500{
1501// List all stored signals in user specified coordinates at the timestamp
1502// of the actual recording of the stored measurement under investigation.
1503// In case no (timestamp of the) actual measurement is available,
1504// the current timestamp of the lab will be taken.
1505// In case no stored signal is available or one of the input arguments is
1506// invalid, no printout is produced.
1507//
1508// The argument "ndig" specifies the number of digits for the fractional
1509// part (e.g. ndig=6 for "dms" corresponds to micro-arcsecond precision).
1510// No rounding will be performed, so an arcsecond count of 3.473 with ndig=1
1511// will appear as 03.4 on the output.
1512// Due to computer accuracy, precision on the pico-arcsecond level may get lost.
1513//
661385a7 1514// The default value is ndig=1.
8bde545d 1515//
1516// Note : The angle info is printed without additional spaces or "endline".
1517// This allows the print to be included in various composite output formats.
1518//
1519// The input parameter "frame" allows the user to specify the frame to which
1520// the coordinates refer. Available options are :
1521//
1522// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d).
1523//
1524// "gal" ==> Galactic coordinates with longitude (l) and lattitude (b).
1525//
1526// "ecl" ==> Ecliptic coordinates with longitude (l) and lattitude (b).
1527//
1528// "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
1529//
1530// "icr" ==> ICRS coordinates with longitude (l) and lattitude (b).
1531//
1532// "loc" ==> Local spherical angles theta and phi at the AliAstrolab location.
1533//
1534// In case the coordinates are the equatorial right ascension and declination (a,d),
1535// they can represent so-called "mean" and "true" values.
1536// The distinction between these two representations is the following :
1537//
1538// mean values : (a,d) are only corrected for precession and not for nutation
1539// true values : (a,d) are corrected for both precession and nutation
1540//
1541// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1542// values for the input in case of equatorial (a,d) coordinates.
1543//
1544// mode = "M" --> Input coordinates are the mean values
1545// "T" --> Input coordinates are the true values
1546//
661385a7 1547// The input parameter "mode" also determines which type of time and
1548// local hour angle will appear in the listing.
1549//
1550// mode = "M" --> Mean Siderial Time (MST) and Local Mean Hour Angle (LMHA)
1551// "T" --> Apparent Siderial Time (AST) and Local Apparent Hour Angle (LAHA)
1552//
8bde545d 1553// Note : In case ts=0 the current timestamp of the lab will be taken.
1554
1555 Int_t iprint=0;
1556
1557 AliTimestamp* tx=0;
1558
661385a7 1559 Int_t dform=1;
1560 if (mode=="T" || mode=="t") dform=-1;
1561
8bde545d 1562 if (fXsig)
1563 {
661385a7 1564 tx=fXsig->GetTimestamp();
8bde545d 1565 if (!tx) tx=(AliTimestamp*)this;
1566 cout << " *AliAstrolab::ListSignals* List of all stored signals." << endl;
661385a7 1567 cout << " === The measurement under investigation ===" << endl;
1568 cout << " Timestamp of the actual observation" << endl;
1569 tx->Date(dform,fToffset);
1570 cout << " Location of the actual observation" << endl;
1571 cout << " "; PrintSignal(frame,mode,tx,ndig); cout << endl;
8bde545d 1572 iprint=1;
1573 }
1574
1575 if (!fRefs) return;
1576
1577 for (Int_t i=1; i<=fRefs->GetSize(); i++)
1578 {
1579 AliSignal* sx=GetSignal(i);
1580 if (!sx) continue;
1581
1582 if (!iprint)
1583 {
1584 cout << " *AliAstrolab::ListRefSignals* List of all stored signals." << endl;
1585 tx=(AliTimestamp*)this;
661385a7 1586 cout << " Current timestamp of the laboratory" << endl;
1587 tx->Date(dform,fToffset);
8bde545d 1588 iprint=1;
1589 }
1590 if (iprint==1)
1591 {
661385a7 1592 cout << " === All stored reference signals according to the above timestamp ===" << endl;
8bde545d 1593 iprint=2;
1594 }
661385a7 1595 cout << " Index : " << i << " "; PrintSignal(frame,mode,tx,ndig,i); cout << endl;
8bde545d 1596 }
1597}
1598///////////////////////////////////////////////////////////////////////////
1599void AliAstrolab::Precess(Ali3Vector& r,AliTimestamp* ts1,AliTimestamp* ts2)
1600{
1601// Correct mean right ascension and declination given for Julian date "jd"
1602// for the earth's precession corresponding to the specified timestamp.
1603// The results are the so-called "mean" (i.e. precession corrected) values,
1604// corresponding to the specified timestamp.
1605// The method used is the new IAU 2000 one as described in the
1606// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1607// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1608// Since the standard reference epoch is J2000, this implies that all
1609// input (a,d) coordinates will be first internally converted to the
1610// corresponding J2000 values before the precession correction w.r.t. the
1611// specified lab timestamp will be applied.
1612//
1613// r : Input vector containing the right ascension and declination information
1614// in the form of standard Ali3Vector coordinates.
1615// In spherical coordinates the phi corresponds to the right ascension,
1616// whereas the declination corresponds to (pi/2)-theta.
1617// jd : Julian date corresponding to the input coordinate values.
1618// ts : Timestamp corresponding to the requested corrected coordinate values.
1619//
1620// Note : In case ts=0 the current timestamp of the lab will be taken.
1621
1622 // Convert back to J2000 values
1623 Ali3Vector r0;
1624 SetPmatrix(ts1);
1625 r0=r.GetUnprimed(&fP);
1626
1627 // Precess to the specified timestamp
1628 if (!ts2) ts2=(AliTimestamp*)this;
1629 SetPmatrix(ts2);
1630 r=r0.GetPrimed(&fP);
1631}
1632///////////////////////////////////////////////////////////////////////////
1633void AliAstrolab::Nutate(Ali3Vector& r,AliTimestamp* ts)
1634{
1635// Correct mean right ascension and declination for the earth's nutation
1636// corresponding to the specified timestamp.
1637// The results are the so-called "true" (i.e. nutation corrected) values,
1638// corresponding to the specified timestamp.
1639// The method used is the new IAU 2000 one as described in the
1640// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1641// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1642//
1643// r : Input vector containing the right ascension and declination information
1644// in the form of standard Ali3Vector coordinates.
1645// In spherical coordinates the phi corresponds to the right ascension,
1646// whereas the declination corresponds to (pi/2)-theta.
1647// ts : Timestamp for which the corrected coordinate values are requested.
1648//
1649// Note : In case ts=0 the current timestamp of the lab will be taken.
1650
1651 // Nutation correction for the specified timestamp
1652 if (!ts) ts=(AliTimestamp*)this;
1653 SetNmatrix(ts);
1654 r=r.GetPrimed(&fN);
1655}
1656///////////////////////////////////////////////////////////////////////////
1657void AliAstrolab::SetBmatrix()
1658{
1659// Set the frame bias matrix elements.
1660// The formulas and numerical constants used are the ones from the
1661// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1662// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1663
1664 Double_t pi=acos(-1.);
1665
1666 // Parameters in mas
1667 Double_t a=-14.6;
1668 Double_t x=-16.6170;
1669 Double_t e=-6.8192;
1670
1671 // Convert to radians
1672 a*=pi/(180.*3600.*1000.);
1673 x*=pi/(180.*3600.*1000.);
1674 e*=pi/(180.*3600.*1000.);
1675
1676 Double_t mat[9];
1677 mat[0]=1.-0.5*(a*a+x*x);
1678 mat[1]=a;
1679 mat[2]=-x;
1680 mat[3]=-a-e*x;
1681 mat[4]=1.-0.5*(a*a+e*e);
1682 mat[5]=-e;
1683 mat[6]=x-e*a;
1684 mat[7]=e+x*a;
1685 mat[8]=1.-0.5*(e*e+x*x);
1686
1687 fB.SetMatrix(mat);
1688 fBias=1;
1689}
1690///////////////////////////////////////////////////////////////////////////
1691void AliAstrolab::SetPmatrix(AliTimestamp* ts)
1692{
1693// Set precession matrix elements for Julian date jd w.r.t. J2000.
1694// The formulas and numerical constants used are the ones from the
1695// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1696// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1697// All numerical constants refer to the standard reference epoch J2000.
1698
1699 Double_t mat[9]={0,0,0,0,0,0,0,0,0};
1700 if (!ts)
1701 {
1702 fP.SetMatrix(mat);
1703 return;
1704 }
1705
1706 Double_t pi=acos(-1.);
1707
1708 Double_t t=(ts->GetJD()-2451545.0)/36525.; // Julian centuries since J2000.0
1709
1710 // Parameters for the precession matrix in arcseconds
1711 Double_t eps0=84381.406; // Mean ecliptic obliquity at J2000.0
1712 Double_t psi=5038.481507*t-1.0790069*pow(t,2)-0.00114045*pow(t,3)+0.000132851*pow(t,4)
1713 -0.0000000951*pow(t,4);
1714 Double_t om=eps0-0.025754*t+0.0512623*pow(t,2)-0.00772503*pow(t,3)-0.000000467*pow(t,4)
1715 +0.0000003337*pow(t,5);
1716 Double_t chi=10.556403*t-2.3814292*pow(t,2)-0.00121197*pow(t,3)+0.000170663*pow(t,4)
1717 -0.0000000560*pow(t,5);
1718
1719 // Convert to radians
1720 eps0*=pi/(180.*3600.);
1721 psi*=pi/(180.*3600.);
1722 om*=pi/(180.*3600.);
1723 chi*=pi/(180.*3600.);
1724
1725 Double_t s1=sin(eps0);
1726 Double_t s2=sin(-psi);
1727 Double_t s3=sin(-om);
1728 Double_t s4=sin(chi);
1729 Double_t c1=cos(eps0);
1730 Double_t c2=cos(-psi);
1731 Double_t c3=cos(-om);
1732 Double_t c4=cos(chi);
1733
1734 mat[0]=c4*c2-s2*s4*c3;
1735 mat[1]=c4*s2*c1+s4*c3*c2*c1-s1*s4*s3;
1736 mat[2]=c4*s2*s1+s4*c3*c2*s1+c1*s4*s3;
1737 mat[3]=-s4*c2-s2*c4*c3;
1738 mat[4]=-s4*s2*c1+c4*c3*c2*c1-s1*c4*s3;
1739 mat[5]=-s4*s2*s1+c4*c3*c2*s1+c1*c4*s3;
1740 mat[6]=s2*s3;
1741 mat[7]=-s3*c2*c1-s1*c3;
1742 mat[8]=-s3*c2*s1+c3*c1;
1743
1744 fP.SetMatrix(mat);
1745}
1746///////////////////////////////////////////////////////////////////////////
1747void AliAstrolab::SetNmatrix(AliTimestamp* ts)
1748{
1749// Set nutation matrix elements for the specified Julian date jd.
1750// The formulas and numerical constants used are the ones from the
1751// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1752// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1753
1754 Double_t mat[9]={0,0,0,0,0,0,0,0,0};
1755 if (!ts)
1756 {
1757 fN.SetMatrix(mat);
1758 return;
1759 }
1760
1761 Double_t pi=acos(-1.);
1762
1763 Double_t dpsi,deps,eps;
1764 ts->Almanac(&dpsi,&deps,&eps);
1765
1766 // Convert to radians
1767 dpsi*=pi/(180.*3600.);
1768 deps*=pi/(180.*3600.);
1769 eps*=pi/(180.*3600.);
1770
1771 Double_t s1=sin(eps);
1772 Double_t s2=sin(-dpsi);
1773 Double_t s3=sin(-(eps+deps));
1774 Double_t c1=cos(eps);
1775 Double_t c2=cos(-dpsi);
1776 Double_t c3=cos(-(eps+deps));
1777
1778 mat[0]=c2;
1779 mat[1]=s2*c1;
1780 mat[2]=s2*s1;
1781 mat[3]=-s2*c3;
1782 mat[4]=c3*c2*c1-s1*s3;
1783 mat[5]=c3*c2*s1+c1*s3;
1784 mat[6]=s2*s3;
1785 mat[7]=-s3*c2*c1-s1*c3;
1786 mat[8]=-s3*c2*s1+c3*c1;
1787
1788 fN.SetMatrix(mat);
1789}
1790///////////////////////////////////////////////////////////////////////////
1791void AliAstrolab::SetGmatrix(TString mode)
1792{
1793// Set the mean equatorial to galactic coordinate conversion matrix.
1794// The B1950 parameters were taken from section 22.3 of the book
1795// "An Introduction to Modern Astrophysics" by Carrol and Ostlie (1996).
1796// The J2000 parameters are obtained by precession of the B1950 values.
1797//
1798// Via the input argument "mode" the required epoch can be selected
1799// mode = "B" ==> B1950
1800// "J" ==> J2000
1801
1802 Ali3Vector x; // The Galactic x-axis in the equatorial frame
1803 Ali3Vector y; // The Galactic y-axis in the equatorial frame
1804 Ali3Vector z; // The Galactic z-axis in the equatorial frame
1805
1806 Double_t a,d;
1807 Double_t vec[3]={1,0,0};
1808
1809 fGal=1; // Set flag to indicate B1950 matrix values
1810
1811 // B1950 equatorial coordinates of the North Galactic Pole (NGP)
1812 a=124900.;
1813 d=272400.;
1814 a=ConvertAngle(a,"hms","deg");
1815 d=ConvertAngle(d,"dms","deg");
1816 vec[1]=90.-d;
1817 vec[2]=a;
1818 z.SetVector(vec,"sph","deg");
1819
1820 // B1950 equatorial coordinates of the Galactic l=b=0 point
1821 a=174224.;
1822 d=-285500.;
1823 a=ConvertAngle(a,"hms","deg");
1824 d=ConvertAngle(d,"dms","deg");
1825 vec[1]=90.-d;
1826 vec[2]=a;
1827 x.SetVector(vec,"sph","deg");
1828
1829 // Precess to the corresponding J2000 values if requested
1830 if (mode=="J")
1831 {
1832 fGal=2; // Set flag to indicate J2000 matrix values
1833 AliTimestamp t1;
1834 t1.SetEpoch(1950,"B");
1835 AliTimestamp t2;
1836 t2.SetEpoch(2000,"J");
1837 Precess(z,&t1,&t2);
1838 Precess(x,&t1,&t2);
1839 }
1840
1841 // The Galactic y-axis is determined for the right handed frame
1842 y=z.Cross(x);
1843
1844 fG.SetAngles(x.GetX(2,"sph","deg"),x.GetX(3,"sph","deg"),
1845 y.GetX(2,"sph","deg"),y.GetX(3,"sph","deg"),
1846 z.GetX(2,"sph","deg"),z.GetX(3,"sph","deg"));
1847}
1848///////////////////////////////////////////////////////////////////////////
1849void AliAstrolab::SetEmatrix(AliTimestamp* ts)
1850{
1851// Set the mean equatorial to ecliptic coordinate conversion matrix
1852// for the specified timestamp.
1853// A nice sketch and explanation of the two frames can be found
1854// in chapter 3 of the book "Astronomy Methods" by Hale Bradt (2004).
1855
1856 Double_t dpsi,deps,eps;
1857 ts->Almanac(&dpsi,&deps,&eps);
1858
1859 // Convert to degrees
1860 eps/=3600.;
1861
1862 // Positions of the ecliptic axes w.r.t. the equatorial ones
1863 // at the moment of the specified timestamp
1864 Double_t theta1=90; // Ecliptic x-axis
1865 Double_t phi1=0;
1866 Double_t theta2=90.-eps; //Ecliptic y-axis
1867 Double_t phi2=90;
1868 Double_t theta3=eps; // Ecliptic z-axis
1869 Double_t phi3=270;
1870
1871 fE.SetAngles(theta1,phi1,theta2,phi2,theta3,phi3);
1872}
1873///////////////////////////////////////////////////////////////////////////
1874void AliAstrolab::SetHmatrix(AliTimestamp* ts)
1875{
1876// Set the mean equatorial to horizontal coordinate conversion matrix
1877// for the specified timestamp.
1878// A nice sketch and explanation of the two frames can be found
1879// in chapter 3 of the book "Astronomy Methods" by Hale Bradt (2004).
1880//
1881// Note : In order to simplify the calculations, we use here a
1882// right-handed horizontal frame.
1883
1884 Ali3Vector x; // The (South pointing) horizontal x-axis in the equatorial frame
1885 Ali3Vector y; // The (East pointing) horizontal y-axis in the equatorial frame
1886 Ali3Vector z; // The (Zenith pointing) horizontal z-axis in the equatorial frame
1887
1888 Double_t l,b;
1889 GetLabPosition(l,b,"deg");
1890
1891 Double_t a;
1892 Double_t vec[3]={1,0,0};
1893
1894 // Equatorial coordinates of the horizontal z-axis
1895 // at the moment of the specified timestamp
1896 a=ts->GetLAST(fToffset);
1897 a*=15.; // Convert fractional hours to degrees
1898 vec[1]=90.-b;
1899 vec[2]=a;
1900 z.SetVector(vec,"sph","deg");
1901
1902 // Equatorial coordinates of the horizontal x-axis
1903 // at the moment of the specified timestamp
1904 vec[1]=180.-b;
1905 vec[2]=a;
1906 x.SetVector(vec,"sph","deg");
1907
1908 // The horizontal y-axis is determined for the right handed frame
1909 y=z.Cross(x);
1910
1911 fH.SetAngles(x.GetX(2,"sph","deg"),x.GetX(3,"sph","deg"),
1912 y.GetX(2,"sph","deg"),y.GetX(3,"sph","deg"),
1913 z.GetX(2,"sph","deg"),z.GetX(3,"sph","deg"));
1914}
1915///////////////////////////////////////////////////////////////////////////
1916void AliAstrolab::SetLocalFrame(Double_t t1,Double_t p1,Double_t t2,Double_t p2,Double_t t3,Double_t p3)
1917{
1918// Specification of the orientations of the local-frame axes.
1919// The input arguments represent the angles (in degrees) of the local-frame axes
1920// w.r.t. a so called Master Reference Frame (MRF), with the same convention
1921// as the input arguments of TRrotMatix::SetAngles.
1922//
1923// The right handed Master Reference Frame is defined as follows :
1924// Z-axis : Points to the local Zenith
1925// X-axis : Makes an angle of 90 degrees with the Z-axis and points South
1926// Y-axis : Makes an angle of 90 degrees with the Z-axis and points East
1927//
1928// Once the user has specified the local reference frame, any observed event
1929// can be related to astronomical space-time locations via the SetSignal
1930// and GetSignal memberfunctions.
1931
1932 // Set the matrix for the conversion of our reference frame coordinates
1933 // into the local-frame ones.
1934
1935 fL.SetAngles(t1,p1,t2,p2,t3,p3);
1936}
1937///////////////////////////////////////////////////////////////////////////
1938Double_t AliAstrolab::ConvertAngle(Double_t a,TString in,TString out) const
1939{
1940// Conversion of various angular formats.
1941//
1942// The input argument "a" denotes the angle to be converted.
1943// The string arguments "in" and "out" specify the angular I/O formats.
1944//
1945// in = "rad" : input angle provided in radians
1946// "deg" : input angle provided in degrees
1947// "dms" : input angle provided in dddmmss.sss
1948// "hms" : input angle provided in hhmmss.sss
1949//
1950// out = "rad" : output angle provided in radians
1951// "deg" : output angle provided in degrees
1952// "dms" : output angle provided in dddmmss.sss
1953// "hms" : output angle provided in hhmmss.sss
1954
1955 if (in==out) return a;
1956
1957 // Convert input to its absolute value in (fractional) degrees.
1958 Double_t pi=acos(-1.);
1959 Double_t epsilon=1.e-12; // Accuracy in (arc)seconds
1960 Int_t word=0,ddd=0,hh=0,mm=0,ss=0;
1961 Double_t s=0;
1962
1963 Double_t b=fabs(a);
1964
1965 if (in=="rad") b*=180./pi;
1966
1967 if (in=="dms")
1968 {
1969 word=Int_t(b);
1970 ddd=word/10000;
1971 word=word%10000;
1972 mm=word/100;
1973 ss=word%100;
1974 s=b-Double_t(ddd*10000+mm*100+ss);
1975 b=Double_t(ddd)+Double_t(mm)/60.+(Double_t(ss)+s)/3600.;
1976 }
1977
1978 if (in=="hms")
1979 {
1980 word=Int_t(b);
1981 hh=word/10000;
1982 word=word%10000;
1983 mm=word/100;
1984 ss=word%100;
1985 s=b-Double_t(hh*10000+mm*100+ss);
1986 b=15.*(Double_t(hh)+Double_t(mm)/60.+(Double_t(ss)+s)/3600.);
1987 }
1988
1989 while (b>360)
1990 {
1991 b-=360.;
1992 }
1993
1994 if (out=="rad") b*=pi/180.;
1995
1996 if (out=="dms")
1997 {
1998 ddd=Int_t(b);
1999 b=b-Double_t(ddd);
2000 b*=60.;
2001 mm=Int_t(b);
2002 b=b-Double_t(mm);
2003 b*=60.;
2004 ss=Int_t(b);
2005 s=b-Double_t(ss);
2006 if (s>(1.-epsilon))
2007 {
2008 s=0.;
2009 ss++;
2010 }
2011 while (ss>=60)
2012 {
2013 ss-=60;
2014 mm++;
2015 }
2016 while (mm>=60)
2017 {
2018 mm-=60;
2019 ddd++;
2020 }
2021 while (ddd>=360)
2022 {
2023 ddd-=360;
2024 }
2025 b=Double_t(10000*ddd+100*mm+ss)+s;
2026 }
2027
2028 if (out=="hms")
2029 {
2030 b/=15.;
2031 hh=Int_t(b);
2032 b=b-Double_t(hh);
2033 b*=60.;
2034 mm=Int_t(b);
2035 b=b-Double_t(mm);
2036 b*=60.;
2037 ss=Int_t(b);
2038 s=b-Double_t(ss);
2039 if (s>(1.-epsilon))
2040 {
2041 s=0.;
2042 ss++;
2043 }
2044 while (ss>=60)
2045 {
2046 ss-=60;
2047 mm++;
2048 }
2049 while (mm>=60)
2050 {
2051 mm-=60;
2052 hh++;
2053 }
2054 while (hh>=24)
2055 {
2056 hh-=24;
2057 }
2058 b=Double_t(10000*hh+100*mm+ss)+s;
2059 }
2060
2061 if (a<0) b=-b;
2062
2063 return b;
2064}
2065///////////////////////////////////////////////////////////////////////////
2066Double_t AliAstrolab::GetHourAngle(TString mode,AliTimestamp* ts,Int_t jref)
2067{
2068// Provide the Local Hour Angle (in fractional degrees) of a stored signal
2069// object at the specified timestamp.
2070//
2071// The input parameter "mode" allows the user to select either the
2072// "mean" or "apparent" value for the returned Hour Angle.
2073//
2074// mode = "M" --> Output is the Mean Hour Angle
2075// "A" --> Output is the Apparent Hour Angle
2076// ts : Timestamp at which the hour angle is requested.
2077//
2078// The input parameter "jref" allows the user to specify so-called "reference" signals.
2079// These reference signals may be used to check space-time event coincidences with the
661385a7 2080// stored measurement (e.g. coincidence of the measurement with transient phenomena).
8bde545d 2081//
661385a7 2082// jref = 0 --> Use the stored measurement
8bde545d 2083// j --> Use the reference signal at the j-th position (j=1 is first)
2084//
2085// Default value is jref=0.
2086//
2087// Note : In case ts=0 the current timestamp of the lab will be taken.
2088
2089 if (!ts) ts=(AliTimestamp*)this;
2090
2091 // Get corrected right ascension and declination for the specified timestamp.
2092 Double_t a,d;
661385a7 2093 if (mode=="M" || mode=="m") GetSignal(a,d,"M",ts,jref);
2094 if (mode=="A" || mode=="a") GetSignal(a,d,"T",ts,jref);
8bde545d 2095
2096 // Convert coordinates to fractional degrees.
2097 a=ConvertAngle(a,"hms","deg");
2098 d=ConvertAngle(d,"dms","deg");
2099
2100 a/=15.; // Convert a to fractional hours
2101 Double_t ha=0;
2102 if (mode=="M" || mode=="m") ha=ts->GetLMST(fToffset)-a;
2103 if (mode=="A" || mode=="a") ha=ts->GetLAST(fToffset)-a;
2104 ha*=15.; // Convert to (fractional) degrees
2105
2106 return ha;
2107}
2108///////////////////////////////////////////////////////////////////////////
2109void 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)
2110{
2111// Set the AliTimestamp parameters corresponding to the LT date and time
2112// in the Gregorian calendar as specified by the input arguments.
2113//
2114// The input arguments represent the following :
2115// y : year in LT (e.g. 1952, 2003 etc...)
2116// m : month in LT (1=jan 2=feb etc...)
2117// d : day in LT (1-31)
2118// hh : elapsed hours in LT (0-23)
2119// mm : elapsed minutes in LT (0-59)
2120// ss : elapsed seconds in LT (0-59)
2121// ns : remaining fractional elapsed second of LT in nanosecond
2122// ps : remaining fractional elapsed nanosecond of LT in picosecond
2123//
2124// Note : ns=0 and ps=0 are the default values.
2125//
2126
2127 SetLT(fToffset,y,m,d,hh,mm,ss,ns,ps);
2128}
2129///////////////////////////////////////////////////////////////////////////
2130void AliAstrolab::SetLT(Int_t y,Int_t d,Int_t s,Int_t ns,Int_t ps)
2131{
2132// Set the AliTimestamp parameters corresponding to the specified elapsed
2133// timespan since the beginning of the new LT year.
2134//
2135// The LT year and elapsed time span is entered via the following input arguments :
2136//
2137// y : year in LT (e.g. 1952, 2003 etc...)
2138// d : elapsed number of days
2139// s : (remaining) elapsed number of seconds
2140// ns : (remaining) elapsed number of nanoseconds
2141// ps : (remaining) elapsed number of picoseconds
2142//
2143// The specified d, s, ns and ps values will be used in an additive
2144// way to determine the elapsed timespan.
2145// So, specification of d=1, s=100, ns=0, ps=0 will result in the
2146// same elapsed time span as d=0, s=24*3600+100, ns=0, ps=0.
2147// However, by making use of the latter the user should take care
2148// of possible integer overflow problems in the input arguments,
2149// which obviously will provide incorrect results.
2150//
2151// Note : ns=0 and ps=0 are the default values.
2152
2153 SetLT(fToffset,y,d,s,ns,ps);
2154}
2155///////////////////////////////////////////////////////////////////////////
2156Double_t AliAstrolab::GetDifference(Int_t j,TString au,Double_t& dt,TString tu,Int_t mode,Int_t* ia,Int_t* it)
2157{
2158// Provide space and time difference between the j-th reference signal
2159// (j=1 indicates first) and the stored measurement.
2160//
2161// The return value of this memberfunction provides the positional angular
2162// difference, whereas the output argument "dt" provides the time difference.
2163//
2164// The units of the angular difference can be specified via the the "au"
2165// input argument, where
2166//
2167// au = "rad" --> Angular difference in (fractional) radians
2168// "deg" --> Angular difference in (fractional) degrees
2169//
2170// The units of the time difference can be specified via the "tu" and "mode"
2171// input arguments. For details please refer to AliTimestamp::GetDifference().
2172// Also here mode=1 is the default value.
2173//
2174// For the time difference the reference signal is used as the standard.
2175// This means that in case of a positive time difference, the stored
2176// measurement occurred later than the reference signal.
2177//
2178// In case j=0, the stored measurement will be compared with each
2179// reference signal and the returned angular and time differences are
2180// the minimal differences which were encountered.
2181// In this case the user may obtain the indices of the two stored reference signals
2182// which had the minimal angular and minimal time difference via the output
2183// arguments "ia" and "it" as follows :
2184//
2185// ia = Index of the stored reference signal with minimial angular difference
2186// it = Index of the stored reference signal with minimial time difference
2187//
2188// In case these indices are the same, there obviously was 1 single reference signal
2189// which showed both the minimal angular and time difference.
2190//
2191// The default values are mode=1, ia=0 and it=0;
2192
2193 Double_t dang=999;
2194 dt=1.e20;
2195 if (ia) *ia=0;
2196 if (it) *it=0;
2197
2198 if (!fRefs) return dang;
2199
2200 Ali3Vector rx; // Position of the measurement
2201 Ali3Vector r0; // Position of the reference signal
2202
2203 AliSignal* sx=GetSignal(rx,"icr","M",0);
2204 if (!sx) return dang;
2205
2206 AliTimestamp* tx=sx->GetTimestamp();
2207 if (!tx) return dang;
2208
2209 // Space and time difference w.r.t. a specific reference signal
2210 if (j>0)
2211 {
2212 AliSignal* s0=GetSignal(r0,"icr","M",0,j);
2213 if (!s0) return dang;
2214 AliTimestamp* t0=s0->GetTimestamp();
2215
2216 if (!t0) return dang;
2217
2218 dang=r0.GetOpeningAngle(rx,au);
2219 dt=t0->GetDifference(tx,tu,mode);
2220 return dang;
2221 }
2222
2223 // Minimal space and time difference encountered over all reference signals
2224 Double_t dangmin=dang;
2225 Double_t dtmin=dt;
2226 for (Int_t i=1; i<=fRefs->GetSize(); i++)
2227 {
2228 AliSignal* s0=GetSignal(r0,"icr","M",0,i);
2229 if (!s0) continue;
2230 AliTimestamp* t0=s0->GetTimestamp();
2231 if (!t0) continue;
2232 dang=r0.GetOpeningAngle(rx,au);
2233 dt=t0->GetDifference(tx,tu,mode);
2234 if (fabs(dang)<dangmin)
2235 {
2236 dangmin=fabs(dang);
2237 *ia=i;
2238 }
2239 if (fabs(dt)<dtmin)
2240 {
2241 dtmin=fabs(dt);
2242 *it=i;
2243 }
2244 }
2245
2246 dt=dtmin;
2247 return dangmin;
2248}
2249///////////////////////////////////////////////////////////////////////////
2250Double_t AliAstrolab::GetDifference(TString name,TString au,Double_t& dt,TString tu,Int_t mode)
2251{
2252// Provide space and time difference between the reference signal
2253// with the specified name and the stored measurement.
2254//
2255// The return value of this memberfunction provides the positional angular
2256// difference, whereas the output argument "dt" provides the time difference.
2257//
2258// The units of the angular difference can be specified via the the "au"
2259// input argument, where
2260//
2261// au = "rad" --> Angular difference in (fractional) radians
2262// "deg" --> Angular difference in (fractional) degrees
2263//
2264// The units of the time difference can be specified via the "tu" and "mode"
2265// input arguments. For details please refer to AliTimestamp::GetDifference().
2266// Also here mode=1 is the default value.
2267//
2268// For the time difference the reference signal is used as the standard.
2269// This means that in case of a positive time difference, the stored
2270// measurement occurred later than the reference signal.
2271
2272 Double_t dang=999;
2273 dt=1.e20;
2274
2275 Int_t j=GetSignalIndex(name);
2276 if (j>0) dang=GetDifference(j,au,dt,tu,mode);
2277 return dang;
2278}
2279///////////////////////////////////////////////////////////////////////////
2280TArrayI* AliAstrolab::MatchRefSignal(Double_t da,TString au,Double_t dt,TString tu,Int_t mode)
2281{
2282// Provide the storage indices of the reference signals which match in space
2283// and time with the stored measurement.
2284// The indices are returned via a pointer to a TArrayI object.
2285// In case no matches were found, the null pointer is returned.
2286// A reference signal is regarded as matching with the stored measurement
2287// if the positional angular difference doesn't exceed "da" and the absolute
2288// value of the time difference doesn't exceed "dt".
2289//
2290// The units of the angular difference "da" can be specified via the the "au"
2291// input argument, where
2292//
2293// au = "rad" --> Angular difference in (fractional) radians
2294// "deg" --> Angular difference in (fractional) degrees
2295//
2296// The units of the time difference "dt" can be specified via the "tu" and "mode"
2297// input arguments. For details please refer to AliTimestamp::GetDifference().
2298// Also here mode=1 is the default value.
2299
2300 if (!fXsig || !fRefs) return 0;
2301
2302 Int_t nrefs=fRefs->GetEntries();
2303
2304 if (fIndices) delete fIndices;
2305 fIndices=new TArrayI(nrefs);
2306
2307 Double_t dang,dtime;
2308 Int_t jfill=0;
2309 for (Int_t i=1; i<=fRefs->GetSize(); i++)
2310 {
2311 dang=GetDifference(i,au,dtime,tu,mode);
2312 if (fabs(dang)<=da && fabs(dtime)<=dt)
2313 {
2314 fIndices->AddAt(i,jfill);
2315 jfill++;
2316 }
2317 }
2318
2319 fIndices->Set(jfill);
2320 return fIndices;
2321}
2322///////////////////////////////////////////////////////////////////////////
2323TObject* AliAstrolab::Clone(const char* name) const
2324{
2325// Make a deep copy of the current object and provide the pointer to the copy.
2326// This memberfunction enables automatic creation of new objects of the
2327// correct type depending on the object type, a feature which may be very useful
2328// for containers when adding objects in case the container owns the objects.
2329
2330 AliAstrolab* lab=new AliAstrolab(*this);
2331 if (name)
2332 {
2333 if (strlen(name)) lab->SetName(name);
2334 }
2335 return lab;
2336}
2337///////////////////////////////////////////////////////////////////////////