X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=RALICE%2FAliTimestamp.cxx;h=c80833dec7fdc775c0d22bcc8b001b92ce184ab1;hb=877fbac6bfb3e4bbfb9c5d66d9ecab028c84d3f8;hp=eab54285af367bd54835e1e175df10494c91fbe8;hpb=2fa2e8164c1b8e1be4c48d80acf7071e59aed54a;p=u%2Fmrichter%2FAliRoot.git

diff --git a/RALICE/AliTimestamp.cxx b/RALICE/AliTimestamp.cxx
index eab54285af3..c80833dec7f 100644
--- a/RALICE/AliTimestamp.cxx
+++ b/RALICE/AliTimestamp.cxx
@@ -167,6 +167,7 @@
 //- Modified: NvE $Date$ Utrecht University.
 ///////////////////////////////////////////////////////////////////////////
 
+#include <cstdlib>
 #include "AliTimestamp.h"
 #include "Riostream.h"
 
@@ -216,11 +217,13 @@ void AliTimestamp::Date(Int_t mode,Double_t offset)
 // mode = 1 ==> Only the UT yy-mm-dd hh:mm:ss.sss and GMST info is printed
 //        2 ==> Only the Julian parameter info is printed
 //        3 ==> Both the UT, GMST and Julian parameter info is printed
+//       -1 ==> Only the UT yy-mm-dd hh:mm:ss.sss and GAST info is printed
+//       -3 ==> Both the UT, GAST and Julian parameter info is printed
 //
 // offset : Local time offset from UT (and also GMST) in fractional hours.
 //
 // When an offset value is specified, the corresponding local times
-// LT and LMST are printed as well.
+// LT and LMST (or LAST) are printed as well.
 //
 // The default values are mode=3 and offset=0.
 //
@@ -235,8 +238,9 @@ void AliTimestamp::Date(Int_t mode,Double_t offset)
  TString day[7]={"Mon","Tue","Wed","Thu","Fri","Sat","Sun"};
  UInt_t y,m,d,wd;
  Int_t hh,mm,ss,ns,ps;
+ Double_t gast;
  
- if (mode==1 || mode==3)
+ if (abs(mode)==1 || abs(mode)==3)
  {
   if (mjd>=40587 && (mjd<65442 || (mjd==65442 && mjsec<8047)))
   {
@@ -253,10 +257,26 @@ void AliTimestamp::Date(Int_t mode,Double_t offset)
   cout << setfill('0') << setw(2) << hh << ":"
        << setw(2) << mm << ":" << setw(2) << ss << "."
        << setw(9) << ns << setw(3) << ps << " (UT)  ";
-  GetGMST(hh,mm,ss,ns,ps);
+  if (mode>0)
+  {
+   GetGMST(hh,mm,ss,ns,ps);
+  }
+  else
+  {
+   gast=GetGAST();
+   Convert(gast,hh,mm,ss,ns,ps);
+  }
   cout << setfill('0') << setw(2) << hh << ":"
        << setw(2) << mm << ":" << setw(2) << ss << "."
-       << setw(9) << ns << setw(3) << ps << " (GMST)" << endl;
+       << setw(9) << ns << setw(3) << ps;
+  if (mode>0)
+  {
+   cout << " (GMST)" << endl;
+  }
+  else
+  {
+   cout << " (GAST)" << endl;
+  }
 
   // Local time information
   if (offset)
@@ -279,11 +299,27 @@ void AliTimestamp::Date(Int_t mode,Double_t offset)
    }
    // Determine the local time by including the offset w.r.t. the original timestamp
    Double_t hlt=GetLT(offset);
-   Double_t hlst=GetLMST(offset);
-   PrintTime(hlt,12); cout << " (LT)  "; PrintTime(hlst,12); cout << " (LMST)" << endl;
+   Double_t hlst=0;
+   if (mode>0)
+   {
+    hlst=GetLMST(offset);
+   }
+   else
+   {
+    hlst=GetLAST(offset);
+   }
+   PrintTime(hlt,12); cout << " (LT)  "; PrintTime(hlst,12);
+   if (mode>0)
+   {
+    cout << " (LMST)" << endl;
+   }
+   else
+   {
+    cout << " (LAST)" << endl;
+   }
   }
  }
- if (mode==2 || mode==3)
+ if (abs(mode)==2 || abs(mode)==3)
  {
   Int_t jd,jsec,jns;
   GetJD(jd,jsec,jns);
@@ -914,7 +950,7 @@ void AliTimestamp::SetMJD(Int_t mjd,Int_t sec,Int_t ns,Int_t ps)
  Int_t limit=65442; // MJD of the latest possible TTimeStamp date/time
  
  Int_t date,time;
- if (mjd<epoch || (mjd>=limit && sec>=8047))
+ if (mjd<epoch || mjd>limit || (mjd==limit && sec>=8047))
  {
   Set(0,kFALSE,0,kFALSE);
   date=GetDate();
@@ -1712,66 +1748,18 @@ Double_t AliTimestamp::GetGAST()
 //
 //      GAST = GMST + (equation of the equinoxes)
 //
-// With the right ascension and declination values to be (0,0) for the
-// vernal equinox, we can workout the shift in right ascension of the vernal
-// equinox due to nutation.
-// The nutation model used is the new one as documented in :
-//   "The IAU Resolutions on Astronomical Reference Systems,
-//    Time Scales and Earth Rotation Models".
-// This document is freely available as Circular 179 (2005) of the
-// United States Naval Observatory (USNO).
-// (See : http://aa.usno.navy.mil/publications/docs).
-// The new expression for the equation of the equinoxes is based on a series
-// expansion and is the most accurate one known to date.
-// The components are documented on p.17 of the USNO Circular 179.
+// The equation of the equinoxes is determined via the Almanac() memberfunction.
 // 
-// Since GMST is based on the MJD, the TTimeStamp limitations
-// do not apply here.
-
- Double_t pi=acos(-1.);
-
- Double_t gmst=GetGMST();
-
- Double_t days;    // Time difference in fractional Julian days w.r.t. the start of J2000.
- Double_t t;       // Time difference in fractional Julian centuries w.r.t. the start of J2000.
- Double_t epsilon; // Mean obliquity of the ecliptic
- Double_t omega;   // Mean longitude of the Moon's mean ascending mode
- Double_t f;       // Mean argument of latitude of the moon
- Double_t d;       // Mean elongation of the Moon from the Sun
-
- days=GetJD()-2451545.0;
- t=days/36525.;
-
- // Values in arcseconds
- epsilon=84381.406-46.836769*t-0.0001831*pow(t,2)+0.00200340*pow(t,3)
-         -0.000000576*pow(t,4)-0.0000000434*pow(t,5);
- omega=450160.398036-6962890.5431*t+7.4722*pow(t,2)+0.007702*pow(t,3)-0.00005939*pow(t,4);
- f=335779.526232+1739527262.8478*t-12.7512*pow(t,2)-0.001037*pow(t,3)+0.00000417*pow(t,4);
- d=1072260.70369+1602961601.2090*t-6.3706*pow(t,2)+0.006593*pow(t,3)-0.00003169*pow(t,4);
-
- // Convert to radians
- epsilon=epsilon*pi/(180.*3600.);
- omega=omega*pi/(180.*3600.);
- f=f*pi/(180.*3600.);
- d=d*pi/(180.*3600.);
+// Since GMST is based on the MJD, the TTimeStamp limitations do not apply here.
 
- Double_t dpsi; // Nutation in longitude in arcseconds
- Double_t beta=(125.04-0.052954*days)*pi/180.;
- Double_t gamma=(200.94+1.97130*days)*pi/180.;
- dpsi=-17.226*sin(beta)-1.296*sin(gamma);
-
- Double_t da; // Right ascension shift of the vernal equinox in arcseconds
- da=dpsi*cos(epsilon)+0.00264096*sin(omega)+0.00006352*sin(2.*omega)
-    +0.00001175*sin(2.*f-2.*d+3.*omega)+0.00001121*sin(2.*f-2.*d+omega)
-    -0.00000455*sin(2.*f-2.*d+2.*omega)+0.00000202*sin(2.*f+3.*omega)+0.00000198*sin(2.*f+omega)
-    -0.00000172*sin(3.*omega)-0.00000087*t*sin(omega);
+ Double_t da=Almanac();
 
  // Convert to fractional hours
- da=da/(3600.*15.);
+ da/=3600.;
 
- Double_t gast=gmst+da;
+ Double_t gast=GetGMST()+da;
 
- while (gast<-24.)
+ while (gast<0)
  {
   gast+=24.;
  }
@@ -1796,7 +1784,7 @@ Double_t AliTimestamp::GetLT(Double_t offset)
  
  h+=offset;
 
- while (h<-24)
+ while (h<0)
  {
   h+=24.;
  }
@@ -1824,7 +1812,7 @@ Double_t AliTimestamp::GetLMST(Double_t offset)
 
  h+=offset;
 
- while (h<-24)
+ while (h<0)
  {
   h+=24.;
  }
@@ -1852,7 +1840,7 @@ Double_t AliTimestamp::GetLAST(Double_t offset)
 
  h+=offset;
 
- while (h<-24)
+ while (h<0)
  {
   h+=24.;
  }
@@ -1985,3 +1973,157 @@ Double_t AliTimestamp::GetTJD(Double_t e,TString mode) const
  return tjd;
 }
 ///////////////////////////////////////////////////////////////////////////
+Double_t AliTimestamp::Almanac(Double_t* dpsi,Double_t* deps,Double_t* eps)
+{
+// Determination of some astronomical observables which may be needed
+// for further calculations like e.g. precession of coordinates.
+//
+// The standard returned value is the "equation of the equinoxes" 
+// (i.e. the nutational shift of the RA of the vernal equinox) in seconds.
+// The memberfunction arguments provide the possibility of retrieving
+// optional returned values. The corresponding observables are :
+//
+// dpsi : Nutational shift in ecliptic longitude in arcseconds
+// deps : Nutational shift in ecliptic obliquity in arcseconds
+// eps  : Mean obliquity of the ecliptic in arcseconds
+//
+// All shifts are determined for the current timestamp with
+// J2000.0 (i.e. 01-jan-2000 12:00:00 UT) as the reference epoch.
+//
+// Invokation example :
+// --------------------
+// AliTimestamp t;
+// Double_t da,dpsi,deps,eps;
+// da=t.Almanac(&dpsi,&deps,&eps); 
+//
+// The nutation model used is the new one as documented in :
+//   "The IAU Resolutions on Astronomical Reference Systems,
+//    Time Scales and Earth Rotation Models".
+// This document is freely available as Circular 179 (2005) of the
+// United States Naval Observatory (USNO).
+// (See : http://aa.usno.navy.mil/publications/docs).
+//
+// The change in ecliptic longitude (dpsi) and ecliptic obliquity (deps)
+// are evaluated using the IAU 2000A nutation series expansion
+// as provided in the USNO Circular 179.
+// The new expression for the equation of the equinoxes is based on a series
+// expansion and is the most accurate one known to date.
+// The components are documented on p.17 of the USNO Circular 179.
+//
+// In the current implementation only the first 28 terms of the nutation series
+// are used. This provides an accuracy of about 0.01 arcsec corresponding to 0.001 sec.
+// In case a better accuracy is required, the series can be extended.
+// The total series expansion consists of 1365 terms.
+// 
+// Since all calculations are based on the JD, the TTimeStamp limitations
+// do not apply here.
+
+ Double_t pi=acos(-1.);
+
+ Double_t t;       // Time difference in fractional Julian centuries w.r.t. the start of J2000.
+ Double_t epsilon; // Mean obliquity of the ecliptic
+ Double_t l;       // Mean anomaly of the Moon
+ Double_t lp;      // Mean anomaly of the Sun
+ Double_t f;       // Mean argument of latitude of the moon
+ Double_t d;       // Mean elongation of the Moon from the Sun
+ Double_t om;      // Mean longitude of the Moon's mean ascending mode
+
+ t=(GetJD()-2451545.0)/36525.;
+
+ // Values of epsilon and the fundamental luni-solar arguments in arcseconds
+ epsilon=84381.406-46.836769*t-0.0001831*pow(t,2)+0.00200340*pow(t,3)
+         -0.000000576*pow(t,4)-0.0000000434*pow(t,5);
+ l=485868.249036+1717915923.2178*t+31.8792*pow(t,2)+0.051635*pow(t,3)-0.00024470*pow(t,4);
+ lp=1287104.79305+129596581.0481*t-0.5532*pow(t,2)+0.000136*pow(t,3)-0.00001149*pow(t,4);
+ f=335779.526232+1739527262.8478*t-12.7512*pow(t,2)-0.001037*pow(t,3)+0.00000417*pow(t,4);
+ d=1072260.70369+1602961601.2090*t-6.3706*pow(t,2)+0.006593*pow(t,3)-0.00003169*pow(t,4);
+ om=450160.398036-6962890.5431*t+7.4722*pow(t,2)+0.007702*pow(t,3)-0.00005939*pow(t,4);
+
+ if (eps) *eps=epsilon;
+
+ // Convert to radians
+ epsilon=epsilon*pi/(180.*3600.);
+ f=f*pi/(180.*3600.);
+ d=d*pi/(180.*3600.);
+ l=l*pi/(180.*3600.);
+ lp=lp*pi/(180.*3600.);
+ om=om*pi/(180.*3600.);
+
+ //The IAU 2000A nutation series expansion.
+ Double_t phi[28]={om,2.*(f-d+om),2.*(f+om),2.*om,lp,lp+2.*(f-d+om),l,
+                   2.*f+om,l+2.*(f+om),2.*(f-d+om)-lp,2.*(f-d)+om,2.*(f+om)-l,2.*d-l,l+om,
+                   om-l,2.*(f+d+om)-l,l+2.*f+om,2.*(f-l)+om,2.*d,2.*(f+d+om),2.*(f-d+om-lp),
+                   2.*(d-l),2.*(l+d+om),l+2.*(f-d+om),2.*f+om-l,2.*l,2.*f,lp+om};
+ Double_t s[28]={-17.2064161,-1.3170907,-0.2276413, 0.2074554, 0.1475877,-0.0516821, 0.0711159,
+                  -0.0387298,-0.0301461, 0.0215829, 0.0128227, 0.0123457, 0.0156994, 0.0063110,
+                  -0.0057976,-0.0059641,-0.0051613, 0.0045893, 0.0063384,-0.0038571, 0.0032481,
+                  -0.0047722,-0.0031046, 0.0028593, 0.0020441, 0.0029243, 0.0025887,-0.0014053};
+ Double_t sd[28]={-0.0174666,-0.0001675,-0.0000234, 0.0000207,-0.0003633, 0.0001226, 0.0000073,
+                  -0.0000367,-0.0000036,-0.0000494, 0.0000137, 0.0000011, 0.0000010, 0.0000063,
+                  -0.0000063,-0.0000011,-0.0000042, 0.0000050, 0.0000011,-0.0000001, 0.0000000,
+                   0.0000000,-0.0000001, 0.0000000, 0.0000021, 0.0000000, 0.0000000,-0.0000025};
+ Double_t cp[28]={ 0.0033386,-0.0013696, 0.0002796,-0.0000698, 0.0011817,-0.0000524,-0.0000872,
+                   0.0000380, 0.0000816, 0.0000111, 0.0000181, 0.0000019,-0.0000168, 0.0000027,
+                  -0.0000189, 0.0000149, 0.0000129, 0.0000031,-0.0000150, 0.0000158, 0.0000000,
+                  -0.0000018, 0.0000131,-0.0000001, 0.0000010,-0.0000074,-0.0000066, 0.0000079};
+ Double_t c[28]= { 9.2052331, 0.5730336, 0.0978459,-0.0897492, 0.0073871, 0.0224386,-0.0006750,
+                   0.0200728, 0.0129025,-0.0095929,-0.0068982,-0.0053311,-0.0001235,-0.0033228,
+                   0.0031429, 0.0025543, 0.0026366,-0.0024236,-0.0001220, 0.0016452,-0.0013870,
+                   0.0000477, 0.0013238,-0.0012338,-0.0010758,-0.0000609,-0.0000550, 0.0008551};
+ Double_t cd[28]={ 0.0009086,-0.0003015,-0.0000485, 0.0000470,-0.0000184,-0.0000677, 0.0000000,
+                   0.0000018,-0.0000063, 0.0000299,-0.0000009, 0.0000032, 0.0000000, 0.0000000,
+                   0.0000000,-0.0000011, 0.0000000,-0.0000010, 0.0000000,-0.0000011, 0.0000000,
+                   0.0000000,-0.0000011, 0.0000010, 0.0000000, 0.0000000, 0.0000000,-0.0000002};
+ Double_t sp[28]={ 0.0015377,-0.0004587, 0.0001374,-0.0000291,-0.0001924,-0.0000174, 0.0000358,
+                   0.0000318, 0.0000367, 0.0000132, 0.0000039,-0.0000004, 0.0000082,-0.0000009,
+                  -0.0000075, 0.0000066, 0.0000078, 0.0000020, 0.0000029, 0.0000068, 0.0000000,
+                  -0.0000025, 0.0000059,-0.0000003,-0.0000003, 0.0000013, 0.0000011,-0.0000045};
+ 
+ Double_t dp=0,de=0,da=0;
+ for (Int_t i=0; i<28; i++)
+ {
+  dp+=(s[i]+sd[i]*t)*sin(phi[i])+cp[i]*cos(phi[i]);
+  de+=(c[i]+cd[i]*t)*cos(phi[i])+sp[i]*sin(phi[i]);
+ }
+
+ da=dp*cos(epsilon)+0.00264096*sin(om)+0.00006352*sin(2.*om)
+     +0.00001175*sin(2.*f-2.*d+3.*om)+0.00001121*sin(2.*f-2.*d+om)
+     -0.00000455*sin(2.*f-2.*d+2.*om)+0.00000202*sin(2.*f+3.*om)+0.00000198*sin(2.*f+om)
+     -0.00000172*sin(3.*om)-0.00000087*t*sin(om);
+
+ if (dpsi) *dpsi=dp;
+ if (deps) *deps=de;
+
+ // Convert to seconds
+ da/=15.;
+
+ return da;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliTimestamp::SetEpoch(Double_t e,TString mode)
+{
+// Set the timestamp parameters according to the epoch as specified by
+// the input argument "e".
+// Via the input argument "mode" the user can specify the type of epoch
+//
+// mode = "B" ==> Besselian epoch
+//        "J" ==> Julian epoch
+
+ Double_t jd=GetJD(e,mode);
+ SetJD(jd);
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliTimestamp::GetEpoch(TString mode)
+{
+// Provide the corresponding epoch value.
+// Via the input argument "mode" the user can specify the type of epoch
+//
+// mode = "B" ==> Besselian epoch
+//        "J" ==> Julian epoch
+
+ Double_t e=0;
+ if (mode=="B" || mode=="b") e=GetBE();
+ if (mode=="J" || mode=="j") e=GetJE();
+ return e;
+}
+///////////////////////////////////////////////////////////////////////////