source: trunk/MagicSoft/Mars/mtemp/MVPObject.cc@ 4308

/* ======================================================================== *\
!
! *
! * 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): Robert Wagner <mailto:magicsoft@rwagner.de> 10/2002
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//  MVPObject                                                              //
//                                                                         //
//  Class used by the visibility plotter to convert RA/Dec to Alt/Az       //
//                                                                         //
//  This class represents an object and is used with the Visibility        //
//  macro. It must be provided with its RA/Dec coordinates and an          //
//  object name (cf. MVPObject::SetRA, MVPObject::SetDec, and              //
//  MVPObject::SetName). Alternatively, you can require the MVPObject      //
//  to be a solar system object like the Sun, Mars or the Moon             //
//  (cf. MVPObject::SetObject).                                            //
//                                                                         //
//  MVPObject is ready to be used in a Mars Eventloop. You must provide    //
//  an Observatory Location as well as a time at which the position        //
//  of the MVPObject is to be calculated. MVPObject::PreProcess            //
//  checks the existence of the required containers and also makes sure    //
//  all necessary setters have been called. MVPObject::Process             //
//  then calculates the Alt/Az position of the object, as well as the      //
//  Zenith angle and the object diameter (Solar system objects).           //
//                                                                         //
//  The astronomical algorithms used are taken from SLALIB 2.4-8.          //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MVPObject.h"

#include <TMath.h>

#include "MLog.h"
#include "MLogManip.h"
#include "MParList.h"

#include "../../slalib/slalib.h"

ClassImp(MVPObject);

// --------------------------------------------------------------------------
//
// Default constructor. 
//
MVPObject::MVPObject(const char *name, const char *title) : fDiameter(0), fCalcEc(kFALSE), fUT1(52000), fBody(10), fGotRA(kFALSE), fGotDec(kFALSE), fGotName(kFALSE)
{
  fName  = name  ? name  : "MVPObject";
  fTitle = title ? title : "Task to calculate Alt, Az of a given object";
  
  fgDegToRad=2*TMath::Pi()/360;
  fgHrsToRad=2*TMath::Pi()/24;  
}

MVPObject::~MVPObject()
{
  //Destructor: nothing special yet.
}

// --------------------------------------------------------------------------
//
//  Check if necessary containers exist in the parameter list already.
//  We need an ObservatoryLocation and a MVPTime object.
//
Bool_t MVPObject::PreProcess(MParList *pList)
{
  fObservatory = (MObservatoryLocation*)pList->FindObject("MObservatoryLocation");
  if (!fObservatory)
    {
      *fLog << dbginf << "MObservatoryLocation not found... aborting." << endl;
      return kFALSE;
    }
  
  fTime = (MVPTime*)pList->FindObject("MVPTime");
  if (!fTime)
    {
      *fLog << dbginf << "MVPTime not found... aborting." << endl;
      return kFALSE;
    }

  if (!fGotRA || !fGotDec || !fGotName)
    {
      *fLog << dbginf << "Object information is not complete." << endl;
      return kFALSE;
    }

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Sets coordinates from object name. Instead of providing RA, Dec and Name
// of an object, you may also just provide the object name in the from
// HHMMsDDT, where RA is given in hours and minutes and Declination is
// given by degrees DD and tenths of degrees T. "s" may be "+" or
// "-", eg. "1959+650"
//
void MVPObject::SetObjectByName(const char *object)
{
  fObjectName=object;
  fGotName=kTRUE;
  
  Int_t ra, dec;
  sscanf(object, "%d%d", &ra, &dec);

  fRA  = fgHrsToRad * (0.01*ra + (r%100)/60.);
  fDec = fgDegToRad * 0.1 * dec;

  fGotRA=kTRUE;
  fGotDec=kTRUE;
}


// --------------------------------------------------------------------------
//
// Sets RA position of object. Position is to be provided in hours, minutes,
// seconds, and microseconds (if needed)
//
void MVPObject::SetRA(Int_t rh, Int_t rm, Int_t rs, Int_t ru)
{
  // Rect is a timelike value...
  fRA = fgHrsToRad*((Double_t)rh + (Double_t)rm/60 + (Double_t)rs/(60*60) + (Double_t)ru/(36000));
  fBody = 10;
  fGotRA = kTRUE;
}

// --------------------------------------------------------------------------
//
// Sets Dec position of object. Position is to be provided in degrees, 
// minutes, seconds, and microseconds (if needed)
//
void MVPObject::SetDec(Int_t dh, Int_t dm, Int_t ds, Int_t du)
{
  // Dec is an anglelike value
  fDec = fgDegToRad*((Double_t)dh + (Double_t)dm/60 + (Double_t)ds/(60*60) + (Double_t)du/(36000));
  fBody = 10;
  fGotDec = kTRUE;
}

// --------------------------------------------------------------------------
//
// Alternatively to providing RA, Dec and Name of an object, you may provide 
// a solar system object (which has no fixed RA, Dec, by the way!) with
// MVPObject::SetObject.
// -
// UInt_t body | Object      Sun and Moon will be objects needed at most,
//           0 | Sun         presumably.
//           1 | Mercury
//           2 | Venus
//           3 | Moon
//           4 | Mars
//           5 | Jupiter
//           6 | Saturn
//           7 | Uranus
//           8 | Neptune
//           9 | Pluto
//
Bool_t MVPObject::SetObject(UInt_t body)
{
  if (body > 9) 
    {
      *fLog << dbginf << "No solar system object associated with value " << body <<"! Ignoring request." << endl;
      return kFALSE;
    }
  else  // We are working on a solar system body.
    {                   
      switch (body) 
        {
        case 1: fObjectName="Mercury"; break;
        case 2: fObjectName="Venus"; break;
        case 3: fObjectName="Moon"; break;
        case 4: fObjectName="Mars"; break;
        case 5: fObjectName="Jupiter"; break;
        case 6: fObjectName="Saturn"; break;
        case 7: fObjectName="Uranus"; break;
        case 8: fObjectName="Neptune"; break;
        case 9: fObjectName="Pluto"; break;
        default: fObjectName="Sun"; 
        }            
    }
  
  fBody = body; 
  fGotRA = fGotDec = fGotName = kTRUE;
  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Given RA, Dec or a solar system object as well as an observatory
// location and a MVPTime, MVPObject::Process() calculates
// Alt, Az, ZA and (in the case of solar system objects) the apparent
// object diameter
//
Bool_t MVPObject::Process()
{
  Double_t diameter = 0.0;

  if (fBody < 10) // We are working on a solar system body.
    {      
      slaRdplan(fTime->GetMJD(), fBody, fObservatory->GetLongitudeRad(), fObservatory->GetLatitudeRad(), &fRA, &fDec, &diameter);         
    }
  
  if (fCalcEc) slaEqecl(fRA, fDec, fTime->GetMJD(), &fEcLong, &fEcLat);
  
  Float_t azimuth;
  Float_t elevation;
  
  Float_t hourAngle = (Float_t)UT1ToGMST(fTime->GetMJD()) - fRA;

  //  cout << "ha: " << hourAngle  << " ra: " << fRA << " dec " << fDec <<endl;

  slaE2h (hourAngle, (Float_t)fDec, (Float_t)fObservatory->GetLatitudeRad(), &azimuth, &elevation);
   
  fZA  = slaZd(hourAngle, fDec, fObservatory->GetLatitudeRad());
  fAlt = (Double_t)elevation;
  fAz  = (Double_t)azimuth; 
  fDiameter = diameter;

  return kTRUE;
}


// --------------------------------------------------------------------------
//
// Returns distance of given object to this object in degrees
//
Double_t MVPObject::GetDistance(MVPObject* object)
{
  return slaSep(fRA, fDec, object->GetRARad(), object->GetDecRad())/fgDegToRad;
}

// --------------------------------------------------------------------------
//
// Returns distance of given object to this object in radians
//
Double_t MVPObject::GetDistanceRad(MVPObject* object)
{
  return slaSep(fRA, fDec, object->GetRARad(), object->GetDecRad());
}


// --------------------------------------------------------------------------
//
// Converts UT1 (given as MJD) to Greenwich mean star time in radians
//
Double_t MVPObject::UT1ToGMST(Double_t ut1)
{
  return slaGmst(ut1);
}

void MVPObject::Print(Option_t *) const
{
  *fLog << all;
  *fLog << "Position of "<< fObjectName << 
    ": Dec " << fDec/fgDegToRad << " deg, " << 
    "RA  " << fRA/fgHrsToRad << " hrs" << endl;
  if (fCalcEc) *fLog << "Ecliptic Long: " << fEcLong/fgDegToRad << " deg, " << 
                  "Ecliptic Lat: " << fEcLat/fgDegToRad << " deg, " << endl; 
}