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

/* ======================================================================== *\
!
! *
! * 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
// "-"
//
void MVPObject::SetObjectByName(char* object)
{
  fObjectName=object;
  fGotName=kTRUE;
  
//   cout<<"OBJ:"<<object<<endl;

  unsigned int delim=0;
  for (unsigned int i=0; i<strlen(object); i++) 
    if ((object[i]=='+')||(object[i]=='-'))
      delim=i;
  
  char ra[6];
  char de[6];
  
  unsigned int i;
  for (i=0;  i<=delim; i++)
    ra[i]=object[i];
  ra[i-1]=0;

  for (i=delim+1;  i<strlen(object); i++)
    de[i-delim-1]=object[i];
  de[i-delim-1]=0;

  Float_t RA, Dec;

  sscanf(ra,"%f",&RA);
  sscanf(de,"%f",&Dec);

//   cout<<"OBJd:"<<Dec<<endl; //220
//   cout<<"OBJr:"<<RA<<endl; //1959

  if (object[delim]=='-') Dec*=-1;

  fRA=(Double_t)( fgHrsToRad*  ((Int_t)(RA/100) + ( RA-(Int_t)(RA/100)*100)/60        ));
  fDec=(Double_t)( fgDegToRad* ((Int_t)(Dec/10) + (Dec-(Int_t)(Dec/10)*10 )/10        ));

 //  fRA=(Double_t)( fgHrsToRad*  ((Int_t)(RA/100)   + ((RA / 100)-(Int_t)(RA/100))/60   ));
//   fDec=(Double_t)( fgDegToRad* ((Int_t)(Dec/10)  + ((Dec / 10)-(Int_t)(Dec/100))/10 ));

//     cout<<"OBJd:"<<fDec/fgDegToRad<<endl;
//     cout<<"OBJr:"<<fRA/fgHrsToRad<<endl;

  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; 
}