source: trunk/MagicSoft/Mars/mastro/MAstro.cc@ 3830

/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz, 11/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MAstro
// ------
//
////////////////////////////////////////////////////////////////////////////
#include "MAstro.h"

#include <iostream>

#include <TVector3.h> // TVector3

#include "MTime.h"    // MTime::GetGmst

using namespace std;

ClassImp(MAstro);

Double_t MAstro::Trunc(Double_t val)
{
    // dint(A) - truncate to nearest whole number towards zero (double)
    return val<0 ? TMath::Ceil(val) : TMath::Floor(val);
}

Double_t MAstro::Round(Double_t val)
{
    // dnint(A) - round to nearest whole number (double)
    return val<0 ? TMath::Ceil(val-0.5) : TMath::Floor(val+0.5);
}

Double_t MAstro::Hms2Sec(Int_t deg, UInt_t min, Double_t sec, Char_t sgn)
{
    const Double_t rc = TMath::Sign((60.0 * (60.0 * (Double_t)TMath::Abs(deg) + (Double_t)min) + sec), (Double_t)deg);
    return sgn=='-' ? -rc : rc;
}

Double_t MAstro::Dms2Rad(Int_t deg, UInt_t min, Double_t sec, Char_t sgn)
{
    // pi/(180*3600):  arcseconds to radians
    //#define DAS2R 4.8481368110953599358991410235794797595635330237270e-6
    return Hms2Sec(deg, min, sec, sgn)*TMath::Pi()/(180*3600)/**DAS2R*/;
}

Double_t MAstro::Hms2Rad(Int_t hor, UInt_t min, Double_t sec, Char_t sgn)
{
    // pi/(12*3600):  seconds of time to radians
//#define DS2R 7.2722052166430399038487115353692196393452995355905e-5
    return Hms2Sec(hor, min, sec, sgn)*TMath::Pi()/(12*3600)/**DS2R*/;
}

Double_t MAstro::Dms2Deg(Int_t deg, UInt_t min, Double_t sec, Char_t sgn)
{
    return Hms2Sec(deg, min, sec, sgn)/3600.;
}

Double_t MAstro::Hms2Deg(Int_t hor, UInt_t min, Double_t sec, Char_t sgn)
{
    return Hms2Sec(hor, min, sec, sgn)/240.;
}

Double_t MAstro::Dms2Hor(Int_t deg, UInt_t min, Double_t sec, Char_t sgn)
{
    return Hms2Sec(deg, min, sec, sgn)/54000.;
}

Double_t MAstro::Hms2Hor(Int_t hor, UInt_t min, Double_t sec, Char_t sgn)
{
    return Hms2Sec(hor, min, sec, sgn)/3600.;
}

void MAstro::Day2Hms(Double_t day, Char_t &sgn, UShort_t &hor, UShort_t &min, UShort_t &sec)
{
    /* Handle sign */
    sgn = day<0?'-':'+';

    /* Round interval and express in smallest units required */
    Double_t a = Round(86400. * TMath::Abs(day)); // Days to seconds

    /* Separate into fields */
    const Double_t ah = Trunc(a/3600.);
    a -= ah * 3600.;
    const Double_t am = Trunc(a/60.);
    a -= am * 60.;
    const Double_t as = Trunc(a);

    /* Return results */
    hor = (UShort_t)ah;
    min = (UShort_t)am;
    sec = (UShort_t)as;
}

void MAstro::Rad2Hms(Double_t rad, Char_t &sgn, UShort_t &deg, UShort_t &min, UShort_t &sec)
{
    Day2Hms(rad/(TMath::Pi()*2), sgn, deg, min, sec);
}

void MAstro::Rad2Dms(Double_t rad, Char_t &sgn, UShort_t &deg, UShort_t &min, UShort_t &sec)
{
    Rad2Hms(rad*15, sgn, deg, min, sec);
}

void MAstro::Deg2Dms(Double_t d, Char_t &sgn, UShort_t &deg, UShort_t &min, UShort_t &sec)
{
    Day2Hms(d/24, sgn, deg, min, sec);
}

void MAstro::Deg2Hms(Double_t d, Char_t &sgn, UShort_t &deg, UShort_t &min, UShort_t &sec)
{
    Rad2Hms(d/360, sgn, deg, min, sec);
}

void MAstro::Hor2Dms(Double_t h, Char_t &sgn, UShort_t &deg, UShort_t &min, UShort_t &sec)
{
    Day2Hms(h*15/24, sgn, deg, min, sec);
}

void MAstro::Hor2Hms(Double_t h, Char_t &sgn, UShort_t &deg, UShort_t &min, UShort_t &sec)
{
    Day2Hms(h/24, sgn, deg, min, sec);
}

void MAstro::Day2Hm(Double_t day, Char_t &sgn, UShort_t &hor, Double_t &min)
{
    /* Handle sign */
    sgn = day<0?'-':'+';

    /* Round interval and express in smallest units required */
    Double_t a = Round(86400. * TMath::Abs(day)); // Days to seconds

    /* Separate into fields */
    const Double_t ah = Trunc(a/3600.);
    a -= ah * 3600.;

    /* Return results */
    hor = (UShort_t)ah;
    min = a/60.;
}

void MAstro::Rad2Hm(Double_t rad, Char_t &sgn, UShort_t &deg, Double_t &min)
{
    Day2Hm(rad/(TMath::Pi()*2), sgn, deg, min);
}

void MAstro::Rad2Dm(Double_t rad, Char_t &sgn, UShort_t &deg, Double_t &min)
{
    Rad2Hm(rad*15, sgn, deg, min);
}

void MAstro::Deg2Dm(Double_t d, Char_t &sgn, UShort_t &deg, Double_t &min)
{
    Day2Hm(d/24, sgn, deg, min);
}

void MAstro::Deg2Hm(Double_t d, Char_t &sgn, UShort_t &deg, Double_t &min)
{
    Rad2Hm(d/360, sgn, deg, min);
}

void MAstro::Hor2Dm(Double_t h, Char_t &sgn, UShort_t &deg, Double_t &min)
{
    Day2Hm(h*15/24, sgn, deg, min);
}

void MAstro::Hor2Hm(Double_t h, Char_t &sgn, UShort_t &deg, Double_t &min)
{
    Day2Hm(h/24, sgn, deg, min);
}

// --------------------------------------------------------------------------
//
// Interpretes a string ' - 12 30 00.0' or '+ 12 30 00.0'
// as floating point value -12.5 or 12.5. If interpretation is
// successfull kTRUE is returned, otherwise kFALSE. ret is not
// touched if interpretation was not successfull. The successfull
// interpreted part is removed from the TString.
//
Bool_t MAstro::String2Angle(TString &str, Double_t &ret)
{
    Char_t  sgn;
    Int_t   d, len;
    UInt_t  m;
    Float_t s;

    // Skip whitespaces before %c and after %f
    int n=sscanf(str.Data(), " %c %d %d %f %n", &sgn, &d, &m, &s, &len);

    if (n!=4 || (sgn!='+' && sgn!='-'))
        return kFALSE;

    str.Remove(0, len);

    ret = Dms2Deg(d, m, s, sgn);
    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Interpretes a string '-12:30:00.0', '12:30:00.0' or '+12:30:00.0'
// as floating point value -12.5, 12.5 or 12.5. If interpretation is
// successfull kTRUE is returned, otherwise kFALSE. ret is not
// touched if interpretation was not successfull.
//
Bool_t MAstro::Coordinate2Angle(const TString &str, Double_t &ret)
{
    Char_t  sgn = str[0]=='-' ? '-' : '+';
    Int_t   d;
    UInt_t  m;
    Float_t s;

    const int n=sscanf(str[0]=='+'||str[0]=='-' ? str.Data()+1 : str.Data(), "%d:%d:%f", &d, &m, &s);

    if (n!=3)
        return kFALSE;

    ret = Dms2Deg(d, m, s, sgn);
    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Return year y, month m and day d corresponding to Mjd.
//
void MAstro::Mjd2Ymd(UInt_t mjd, UShort_t &y, Byte_t &m, Byte_t &d)
{
    // Express day in Gregorian calendar
    const ULong_t jd   = mjd + 2400001;
    const ULong_t n4   = 4*(jd+((6*((4*jd-17918)/146097))/4+1)/2-37);
    const ULong_t nd10 = 10*(((n4-237)%1461)/4)+5;

    y = n4/1461L-4712;
    m = ((nd10/306+2)%12)+1;
    d = (nd10%306)/10+1;
}

// --------------------------------------------------------------------------
//
//  Return Mjd corresponding to year y, month m and day d.
//
Int_t MAstro::Ymd2Mjd(UShort_t y, Byte_t m, Byte_t d)
{
    // Month lengths in days
    static int months[12] = { 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 };

    // Validate month
    if (m<1 || m>12)
        return -1;

    // Allow for leap year
    months[1] = (y%4==0 && (y%100!=0 || y%400==0)) ? 29 : 28;

    // Validate day
    if (d<1 || d>months[m-1])
        return -1;

    // Precalculate some values
    const Byte_t  lm = 12-m;
    const ULong_t lm10 = 4712 + y - lm/10;

    // Perform the conversion
    return 1461L*lm10/4 + (306*((m+9)%12)+5)/10 - (3*((lm10+188)/100))/4 + d - 2399904;
}

// --------------------------------------------------------------------------
//
//  theta0, phi0    [rad]: polar angle/zenith distance, azimuth of 1st object
//  theta1, phi1    [rad]: polar angle/zenith distance, azimuth of 2nd object
//  AngularDistance [rad]: Angular distance between two objects
//
Double_t MAstro::AngularDistance(Double_t theta0, Double_t phi0, Double_t theta1, Double_t phi1)
{
    TVector3 v0(1);
    v0.Rotate(phi0, theta0);

    TVector3 v1(1);
    v1.Rotate(phi1, theta1);

    return v0.Angle(v1);
}

// --------------------------------------------------------------------------
//
// Calls MTime::GetGmst() Better use MTime::GetGmst() directly
//
Double_t MAstro::UT2GMST(Double_t ut1)
{
    return MTime(ut1).GetGmst();
}

// --------------------------------------------------------------------------
//
// RotationAngle 
// 
// calculates the angle for the rotation of the sky coordinate system
// with respect to the local coordinate system. This is identical
// to the rotation angle of the sky image in the camera.
//
//  sinl  [rad]: sine of observers latitude
//  cosl  [rad]: cosine of observers latitude
//  theta [rad]: polar angle/zenith distance
//  phi   [rad]: rotation angle/azimuth
//
// Return sin/cos component of angle
//
// The convention is such, that the rotation angle is -pi/pi if
// right ascension and local rotation angle are counted in the
// same direction, 0 if counted in the opposite direction.
//
// (In other words: The rotation angle is 0 when the source culminates)
//
// Using vectors it can be done like:
//    TVector3 v, p;
//    v.SetMagThetaPhi(1, theta, phi);
//    p.SetMagThetaPhi(1, TMath::Pi()/2-latitude, 0);
//    v = v.Cross(l)); 
//    v.RotateZ(-phi);
//    v.Rotate(-theta)
//    rho = TMath::ATan2(v(2), v(1));
//
// For more information see TDAS 00-11, eqs. (18) and (20)
//
void MAstro::RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi, Double_t &sin, Double_t &cos)
{
    const Double_t sint = TMath::Sin(theta);
    const Double_t cost = TMath::Cos(theta);

    const Double_t snlt = sinl*sint;
    const Double_t cslt = cosl*cost;

    const Double_t sinp = TMath::Sin(phi);
    const Double_t cosp = TMath::Cos(phi);

    const Double_t v1 = sint*sinp;
    const Double_t v2 = cslt - snlt*cosp;

    const Double_t denom = TMath::Sqrt(v1*v1 + v2*v2);

    sin =   cosl*sinp      / denom; // y-component
    cos = (snlt-cslt*cosp) / denom; // x-component
}

// --------------------------------------------------------------------------
//
// RotationAngle 
// 
// calculates the angle for the rotation of the sky coordinate system
// with respect to the local coordinate system. This is identical
// to the rotation angle of the sky image in the camera.
//
//  sinl  [rad]: sine of observers latitude
//  cosl  [rad]: cosine of observers latitude
//  theta [rad]: polar angle/zenith distance
//  phi   [rad]: rotation angle/azimuth
//
// Return angle [rad] in the range -pi, pi
//
// The convention is such, that the rotation angle is -pi/pi if
// right ascension and local rotation angle are counted in the
// same direction, 0 if counted in the opposite direction.
//
// (In other words: The rotation angle is 0 when the source culminates)
//
// Using vectors it can be done like:
//    TVector3 v, p;
//    v.SetMagThetaPhi(1, theta, phi);
//    p.SetMagThetaPhi(1, TMath::Pi()/2-latitude, 0);
//    v = v.Cross(l)); 
//    v.RotateZ(-phi);
//    v.Rotate(-theta)
//    rho = TMath::ATan2(v(2), v(1));
//
// For more information see TDAS 00-11, eqs. (18) and (20)
//
Double_t MAstro::RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi)
{
    const Double_t snlt = sinl*TMath::Sin(theta);
    const Double_t cslt = cosl*TMath::Cos(theta);

    const Double_t sinp = TMath::Sin(phi);
    const Double_t cosp = TMath::Cos(phi);

    return TMath::ATan2(cosl*sinp, snlt-cslt*cosp);
}