AstCoordinate.cxx

Go to the documentation of this file.

//
// AstCoordinate
//
// Package: Ast (AstroUtil).
//
// S. Kasahara 05/03
//
// Purpose: Provide utility routines to convert between various astronomical
//          coordinate sytems.
//
//          Horizon coordinates (altitude,azimuth) are defined as:
//          altitude:  elevation above the horizon, measured from -90^o (nadir)
//                     through 0^o (horizon) to +90^o (zenith).
//          azimuth:  the angle in the horizon plane measured easterly from
//                    north, i.e. 0^o is North, 90^o East, 180^o South,
//                    270^o West.
//
//          Local Det. Coordinate direction cosines (dcosx,dcosy,dcosz) are 
//          defined in the local coordinate system of a specified detector
//          type (Detector::Detector_t):
//          +z: horizontal component of the beam direction.
//          +y: local vertical.
//          +x: chosen to make coordinate system right-handed.
//
//          Ideal Det. Coordinate direction cosines (dcosx_ideal,dcosy_ideal,
//          dcosz_ideal) are defined such that:
//          +z: points North.
//          +y: radius vector originating at center of earth.
//          +x: points West.   
//
//          Geocentric Equatorial coordinates (hourangle, declination),
//          are defined relative to the earth's equatorial plane.
//          hourangle:  time since the object crossed the observer's meridian 
//                      (the great circle passing through the observer's zenith
//                       and the poles).  Measured from 0 to 360^o.
//          declination: elevation above the equatorial plane, measured from
//                       -90^o to 90^o.
//
//          Geocentric Celestial coordinates (right ascension,declination),
//          are defined as in the equatorial system, but instead of rotating
//          with the earth, they are fixed on the celestial sphere.
//          right ascension: angle in the equatorial plane measured eastwards
//                           from the rising of Aries. Measured from 0 - 360^o.
//          declination: same as for Geocentric Equatorial coordinates.
//
//Panos, 10/04
//          Geocentric Ecliptic coordinates (ecliptic latitude, beta, ecliptic 
//          longitude, lamda), are defined relative to the ecliptic great circle.
//          They are also fixed on the celestial sphere 
//          latitude: elevation over the ecliptic (-90 to 90^o).
//          longitude: angle in the ecliptic plane, measured eastwards from the 
//          rising of Aries (one of two points where the equatorial and ecliptic
//          circles cross).

#include <math.h>
#include "AstroUtil/Ast.h"
#include "AstroUtil/AstCoordinate.h"
#include "AstroUtil/AstTime.h"
#include "Conventions/Mphysical.h"  // pi
#include "MessageService/MsgService.h"

CVSID("$Id: AstCoordinate.cxx,v 1.5 2005/12/13 15:57:30 rhatcher Exp $");

// Definition of methods (alphabetical order)
// ******************************************

void AstUtil::LocalToIdeal( double dcosx, double dcosy, double dcosz, 
     Detector::Detector_t dettype, double& dcosx_ideal, 
     double& dcosy_ideal, double& dcosz_ideal ) {
  // Static method to convert direction cosines from local detector
  // coordinates to ideal detector coordinates.

  const double* rotmatrix = AstUtil::GetDetRotMatrixLocalToIdeal(dettype);
  if ( !rotmatrix ) {
    MSG("Ast",Msg::kWarning) 
      << "No rotation matrix available for requested det type " << dettype
      << ".\n No rotation from local to ideal coordinates will be applied." 
      << endl;
    dcosx_ideal = dcosx;
    dcosy_ideal = dcosy;
    dcosz_ideal = dcosz;
    return;
  }

  dcosx_ideal = rotmatrix[0]*dcosx
              + rotmatrix[1]*dcosy
              + rotmatrix[2]*dcosz;
  dcosy_ideal = rotmatrix[3]*dcosx
              + rotmatrix[4]*dcosy
              + rotmatrix[5]*dcosz;
  dcosz_ideal = rotmatrix[6]*dcosx
              + rotmatrix[7]*dcosy
              + rotmatrix[8]*dcosz;
  return;

}

void AstUtil::IdealToLocal( double dcosx_ideal, double dcosy_ideal,
      double dcosz_ideal, Detector::Detector_t dettype,
      double& dcosx, double& dcosy, double& dcosz) {
  // Static method to convert direction cosines from ideal detector 
  // coordinates to local detector coordinates.

  const double* rotmatrix = AstUtil::GetDetRotMatrixIdealToLocal(dettype);
  if ( !rotmatrix ) {
    MSG("Ast",Msg::kWarning) 
      << "No rotation matrix available for requested det type " << dettype
      << ".\n No rotation from ideal to local coordinates  will be applied." 
      << endl;
    dcosx = dcosx_ideal;
    dcosy = dcosy_ideal;
    dcosz = dcosz_ideal;
    return;
  }

  dcosx = rotmatrix[0]*dcosx_ideal
        + rotmatrix[1]*dcosy_ideal
        + rotmatrix[2]*dcosz_ideal;
  dcosy = rotmatrix[3]*dcosx_ideal
        + rotmatrix[4]*dcosy_ideal
        + rotmatrix[5]*dcosz_ideal;
  dcosz = rotmatrix[6]*dcosx_ideal
        + rotmatrix[7]*dcosy_ideal
        + rotmatrix[8]*dcosz_ideal;
  return;

}

void AstUtil::LocalToHorizon(double dcosx,double dcosy,double dcosz, 
     Detector::Detector_t dettype, double& altitude, double& azimuth){
  // Static method to convert direction cosines to horizon coordinates
  // The direction cosines are local to the detector type specified by dettype.


  double dcosx_ideal,dcosy_ideal,dcosz_ideal;

  AstUtil::LocalToIdeal(dcosx,dcosy,dcosz,dettype,dcosx_ideal,
                         dcosy_ideal,dcosz_ideal);
  AstUtil::IdealToHorizon(dcosx_ideal,dcosy_ideal,dcosz_ideal,altitude,
                             azimuth);

}

void AstUtil::IdealToHorizon(double dcosx_ideal,double dcosy_ideal,
                 double dcosz_ideal,double& altitude, double& azimuth){
  // Static method to convert ideal coordinate system direction cosines to 
  // horizon coordinates.

  double zenith = acos(dcosy_ideal)*180./Mphysical::pi;
  altitude = 90. - zenith;
  azimuth = 180.*atan2(-dcosx_ideal,dcosz_ideal)/Mphysical::pi;
  if ( azimuth < 0. ) azimuth += 360.;

}

void AstUtil::HorizonToIdeal(double altitude, double azimuth,
    double& dcosx_ideal, double& dcosy_ideal, double& dcosz_ideal) {
  // Static method to convert horizon coordinates to ideal detector
  // coordinates.
  double zenrad = (90. - altitude)*Mphysical::pi/180.; 
  double azirad = azimuth*Mphysical::pi/180.;
  dcosy_ideal =  cos(zenrad);
  dcosx_ideal = -sin(zenrad)*sin(azirad);
  dcosz_ideal =  sin(zenrad)*cos(azirad);
  return;

}

void AstUtil::HorizonToLocal(double altitude, double azimuth,
          Detector::Detector_t dettype, double& dcosx, double& dcosy, 
          double& dcosz) {
  // Static method to convert horizon coordinates to direction cosines

  double dcosx_ideal,dcosy_ideal,dcosz_ideal;
  AstUtil::HorizonToIdeal(altitude,azimuth,dcosx_ideal,dcosy_ideal,
                                                       dcosz_ideal);
  AstUtil::IdealToLocal(dcosx_ideal,dcosy_ideal,dcosz_ideal,dettype,
                        dcosx,dcosy,dcosz);
  return;
}

void AstUtil::HorizonToEquatorial(double altitude, double azimuth, 
              double latitude, double& hourangle, double& declination) {
  // Static method to convert from horizon to equatorial coordinates.
  // From Practical Ephemeris Calculations by O. Montenbruck:
  // 1)sin(dec)         = sin(lat)*sin(alt) + cos(lat)*cos(alt)*cos(azi)
  // 2)cos(dec)*sin(ha) = -cos(alt)*sin(azi)
  // 3)cos(dec)*cos(ha) = cos(lat)*sin(alt)-sin(lat)*cos(alt)*cos(azi)
  // This implementation is adapted from an M. Thomson Soudan2 fortran
  // routine.

  double latrad = latitude*Mphysical::pi/180.;
  double altrad = altitude*Mphysical::pi/180.;
  double azirad = azimuth*Mphysical::pi/180.;

  double sinlat = sin(latrad);
  double sinalt = sin(altrad);
  double sinazi = sin(azirad);
  double coslat = cos(latrad);
  double cosalt = cos(altrad);
  double cosazi = cos(azirad);

  // solve equation 1 to determine declination
  double sindec = sinalt*sinlat + cosalt*cosazi*coslat;
  if ( sindec >= 1 ) declination = Mphysical::pi/2.;
  else if ( sindec <= -1 ) declination = -Mphysical::pi/2.;
  else declination = asin(sindec);

  // solve equation 2 to determine hourangle in range -pi/2 to pi/2
  double du = 90.0001*Mphysical::pi/180.;
  double dl = 89.9999*Mphysical::pi/180.;
  if ( fabs(declination) > dl && fabs(declination) < du ) {
    hourangle = 0.;
  }
  else {
   double sinhourangle = -cosalt*sinazi/cos(declination);
   if ( sinhourangle >= 1 ) hourangle = Mphysical::pi/2.;
   else if ( sinhourangle <= -1 ) hourangle = -Mphysical::pi/2.;
   else hourangle = asin(sinhourangle);
  }

  // use 3rd equation to fully determine hourangle
  double lhs = cos(declination)*cos(hourangle);
  double rhs = sinalt*coslat - cosalt*cosazi*sinlat;

  double prod = lhs*rhs;
  if ( prod != fabs(prod) ) {
    if ( hourangle > 0. ) hourangle = Mphysical::pi - hourangle;
    else hourangle = -Mphysical::pi - hourangle;
  }

  if ( hourangle < 0. ) hourangle += 2.*Mphysical::pi;

  hourangle   = hourangle*180./Mphysical::pi; // degrees
  declination = declination*180./Mphysical::pi; // degrees

  return;

}

void AstUtil::EquatorialToHorizon(double hourangle,double declination, 
    double latitude, double& altitude, double& azimuth) {
  // Static method to convert from equatorial to horizon coordinates
  // From Practical Ephemeris Calculations by O. Montenbruck:
  // 1)sin(alt)          =  sin(lat)*sin(dec) + cos(lat)*cos(dec)*cos(ha)
  // 2)cos(alt)*sin(azi) = -cos(dec)*sin(ha)
  // 3)cos(alt)*cos(azi) =  cos(lat)*sin(dec) - sin(lat)*cos(dec)*cos(ha)
  // This implementation is adapted from an M. Thomson Soudan2 fortran
  // routine.

  double latrad = latitude*Mphysical::pi/180.;
  double harad  = hourangle*Mphysical::pi/180.;
  double decrad = declination*Mphysical::pi/180.;

  double sinlat = sin(latrad);
  double sinha  = sin(harad);
  double sindec = sin(decrad);
  double coslat = cos(latrad);
  double cosha  = cos(harad);
  double cosdec = cos(decrad);

  // Solve equation 1 for altitude
  double sinalt = sindec*sinlat + cosdec*cosha*coslat;
  if ( sinalt >= 1 ) altitude = Mphysical::pi/2.;
  else if ( sinalt <= -1 ) altitude = -Mphysical::pi/2.;
  else altitude = asin(sinalt);

  // Solve equation 2 to determine azimuth in range -pi/2 to pi/2.
  double du = 90.0001*Mphysical::pi/180.;
  double dl = 89.9999*Mphysical::pi/180.;
  if (fabs(altitude) > dl && fabs(altitude) < du ) {
    azimuth = 0;
  }
  else {
    double sinazi = -cosdec*sinha/cos(altitude);
    if ( sinazi >= 1 ) azimuth = Mphysical::pi/2.;
    else if (sinazi <= -1 ) azimuth = -Mphysical::pi/2.; 
    else azimuth = asin(sinazi);
  }

  // Use 3rd equation to fully determine azimuth
  double lhs = cos(altitude)*cos(azimuth);
  double rhs = sindec*coslat - cosdec*cosha*sinlat;

  double prod = lhs*rhs;
  if ( prod != fabs(prod) ) {
    if ( azimuth > 0. ) azimuth = Mphysical::pi - azimuth;
    else azimuth = -Mphysical::pi - azimuth;
  }

  if ( azimuth < 0. ) azimuth += 2.*Mphysical::pi;

  altitude = altitude*180./Mphysical::pi;  // degrees
  azimuth  = azimuth*180./Mphysical::pi;   // degrees

  return;

}

void AstUtil::CelestialToEquatorial(double ra, double gst,
                          double longitude,double& hourangle) {
  // Static method to convert from celestial to equatorial coordinates

  double lst;
  AstUtil::GSTToLST(gst,longitude,lst);
  lst *= 15.; // convert from hours to degrees
  hourangle = lst - ra;
  if ( hourangle < 0. ) hourangle += 360.;
 
  return;

}

void AstUtil::EquatorialToCelestial(double hourangle, double gst,
                                    double longitude, double& ra) {
  // Static method to convert from equatorial to celestial coordinates

  double lst; 
  AstUtil::GSTToLST(gst,longitude,lst);
  lst *= 15.; // convert from hours to degrees
  ra = lst - hourangle; // in degrees
  if ( ra < 0. ) ra += 360.;

  return;
}

void AstUtil::CelestialToEcliptic(double declination, double ra,
                double& beta, double& lamda)
{ // Static method to convert from celestial to ecliptic coordinates

    const double epsilon  = 23.43;
    
    double sindelta = sin(declination);
    double sinalpha = sin(ra);
    double cosdelta = cos(declination);
    double cosalpha = cos(ra);
    double cosbeta  = cos(beta);
    
    //We will need both trigonometric numbers of epsilon, the constant angle between
    //the ecliptic and the celestial equator great circles. 
    const double Alpha = cos(Mphysical::pi*epsilon/180.);
    const double Beta  = sin(Mphysical::pi*epsilon/180.);
    
    beta  =  asin(sindelta*Alpha-cosdelta*Beta*sinalpha);
    if (cosbeta!=0){
        lamda =  acos(cosalpha*cosdelta/cosbeta); 
    } else { //the point lies on pole of ecliptic, and lamda has no meaning
        lamda = 0;
    }
    if (lamda < 0.) lamda += 360.;
    return;
}


void AstUtil::EclipticToCelestial(double beta, double lamda,
                double& declination, double& ra)
{ // Static method to convert from ecliptic to celestial coordinates

    const double epsilon  = 23.43;
    
    double sinlamda = sin(lamda);
    double coslamda = cos(lamda);
    double cosbeta  = cos(beta);
    double sinbeta  = sin(beta);
    double cosdelta = cos(declination);
    
    //We will need both trigonometric numbers of epsilon, the constant angle between
    //the ecliptic and the celestial equator great circles. 
    const double Alpha = cos(Mphysical::pi*epsilon/180.);
    const double Beta  = sin(Mphysical::pi*epsilon/180.);
    
    declination  =  asin(sinbeta*Alpha+cosbeta*Beta*sinlamda);
    if (cosdelta!=0){
        ra =  acos(coslamda*cosbeta/cosdelta); 
    } else { //the point lies on pole of equator, and ra has no meaning
        ra = 0;
    }
    if (ra <0.) ra += 360.;
    return;
}

---