position_angle.cpp

#include <cmath>
#include <string>
#include <iostream>

#include <sash/HESSArray.hh>
#include <crash/RotatedSystem.hh>
#include <stash/Coordinate.hh>

/**
   Convert position angles between Equatorial (icrs) and Galactic 
   coordinates ( all parameters in degrees ):

   Input parameters are the original position, position angle and system:
   l = longitude coordinate
   b = latitude coordinate
   pa = position angle
   sys = system (can be "icrs" or "galactic")
   eps = step size (in deg) used to construct point 2 in the algorithm
   ( default eps should be fine if you only need a precision
   of eps and are not near a pole. )

   Output parameter is the position angle in the other system,
   restricted to be in the range 0 to 360 deg:
   new_pa

   Note: In both systems position angles are defined "east of north",
   i.e. the positive b direction has pa = 0 deg and the positive
   l direction has pa = 90 deg.

   I have tested this function for a few random values back and forth,
   it gave identical results to the NED Coordinate Transformation Calculator
   at http://nedwww.ipac.caltech.edu/forms/calculator.html
*/
double ConvertPositionAngle( double l, double b, 
			     double pa, std::string sys,
			     double eps = 1.e-3 ) {
  /**
     The following algorithm is used:
     Create a second position in the original system by
     stepping from the first position a small angle eps into the direction
     specified by the position angle.
     Then transform both positions and measure the position
     angle of the second position relative to the first position.
  */

  // Convert all angles from deg to rad
  const double RAD = M_PI / 180; // rad per deg
  //  l *= RAD; b *= RAD; pa *= RAD; eps *=RAD;

  // Construct second position
  // Note: The extra factor of cos(new_b) is necessary such that
  double l2 = l + eps * sin( RAD*pa ) / cos( RAD*b );
  double b2 = b + eps * cos( RAD*pa );
  double dist = sqrt( pow(l2-l,2) + pow(b2-b,2) );
  std::cout << "l       = " << l  << std::endl;
  std::cout << "b       = " << b  << std::endl;
  std::cout << "l2      = " << l2  << std::endl;
  std::cout << "b2      = " << b2  << std::endl;
  std::cout << "dist    = " << dist  << std::endl;
  std::cout << "pa      = " << pa  << std::endl;

  // Construct original and new coordinate systems 
  Sash::HESSArray *hess = &Sash::HESSArray::GetHESSArray();
  Stash::SystemRef* sys_ref = (Stash::SystemRef*)
    &Crash::RotatedSystem::GetGalacticSystem();
  Stash::SystemRef* new_sys_ref = (Stash::SystemRef*)
    &hess->GetRADecJ2000System();
  if ( sys == "galactic") {
  } else if ( sys == "icrs" ) {
    // switch them
     Stash::SystemRef* temp = sys_ref;
     sys_ref = new_sys_ref;
     new_sys_ref = temp;
  } else {
    std::cerr << "Coordinate system not supported: "
	      << sys << std::endl; 
    return -1;
  }
  std::cout << "sys_ref      = " << **sys_ref << endl;
  std::cout << "new_sys_ref  = " << **new_sys_ref << endl;

  // Construct coordinates
  Stash::Lambda lc ( l,  Stash::Angle::Degrees );
  Stash::Beta   bc ( b,  Stash::Angle::Degrees );
  Stash::Lambda lc2( l2, Stash::Angle::Degrees );
  Stash::Beta   bc2( b2, Stash::Angle::Degrees );
  Stash::Coordinate p( lc, bc, *sys_ref );
  Stash::Coordinate p2( lc2, bc2, *sys_ref );
  // Convert coordinates to new system
  Stash::Coordinate new_p  = p.GetCoordinate ( **new_sys_ref );
  Stash::Coordinate new_p2 = p2.GetCoordinate( **new_sys_ref );
  // Get new position angle vector as the difference
  // between the new coordinates
  double new_l  = new_p.GetLambda().GetDegrees();
  double new_l2 = new_p2.GetLambda().GetDegrees();
  double new_b  = new_p.GetBeta().GetDegrees();
  double new_b2 = new_p2.GetBeta().GetDegrees();

  // Note: The extra factor of cos(new_b) is necessary to 
  double delta_l = ( new_l2 - new_l ) * cos( RAD*new_b );
  double delta_b = ( new_b2 - new_b );
  double new_dist = sqrt( pow(delta_l,2) + pow(delta_b,2) );

  // Compute new position angle
  double new_pa = 0;
  const double almost_zero = eps * 1.e-6; // to prevent dividing by almost zero
  if ( abs(delta_b) > almost_zero ) {
    new_pa = atan( delta_l / delta_b ) / RAD; // atan returns range -pi to pi
    if ( new_pa < 0 ) new_pa += 360; // force range 0 to 360 deg
  } else {
    if ( delta_l > 0 )
      new_pa = 90;
    else
      new_pa = 270;
  }

  std::cout << "p       = " << p  << std::endl;
  std::cout << "p2      = " << p2 << std::endl;
  std::cout << "new_p   = " << new_p  << std::endl;
  std::cout << "new_p2  = " << new_p2  << std::endl;

  std::cout << "new_l   = " << new_l  << std::endl;
  std::cout << "new_b   = " << new_b  << std::endl;
  std::cout << "new_l2  = " << new_l2  << std::endl;
  std::cout << "new_b2  = " << new_b2  << std::endl;
  std::cout << "delta_l = " << delta_l << std::endl;
  std::cout << "delta_b = " << delta_b << std::endl;
  std::cout << "new_dist= " << new_dist << std::endl;

  std::cout << "new_pa  = " << new_pa << std::endl;

  return new_pa;
}