Azimuth and Altitude of a body from earth

[alexibu]

Hi,
Ive created a method to calculate the position of a celestial body in the sky given the observers location on a planet and the time.
Im not sure how accurate it is because I don't know what to trust more (celestia and my calcs - possibly with precision probs) or other sky map programs.
This would be useful to people trying to control telescopes with celestia.

Code: Select all

void Body::getDirectionOfBodyFromSurface(const Vec3f& lonLatAlt,Body * body,double when,
                               double& Dist,double& Azimuth,double&Altitude)
{
   //get the bodies position in geographic coodinates
   Point3d bodyPointd = body->getHeliocentricPosition(when)
      *getGeographicToHeliocentric(when).inverse();

   //get the position of the observer in geographic coords
   Vec3f observerVector = planetocentricToCartesian(lonLatAlt);

   //create a plane tangential to the earths surface where the
   //observer is
   Vec3f FloorN = observerVector;
   FloorN.normalize();
   Planef FloorPlane(FloorN,observerVector.length());

   //find a point north of the observer by finding the floor planes
   //intersection with the geographical y axis
   Point3f northPoint = Planef::intersection(FloorPlane,
         Planef(Vec3f(1,0,0),0),Planef(Vec3f(0,0,1),0));
   
   //Do some precision conversions
   Point3f bodyPoint(bodyPointd.x,bodyPointd.y,bodyPointd.z);
   Vec3f bodyVector(bodyPoint.x,bodyPoint.y,bodyPoint.z);
   Vec3f northVector(northPoint.x,northPoint.y,northPoint.z);
   
   //calculate the vector from the observer to north
   northVector -= observerVector;

   //create a plane that splits the planet north south
   //from the point of view of the observer
   Vec3f NSPlaneN = northVector^observerVector;
   NSPlaneN.normalize();
   Planef NorthSouthPlane(NSPlaneN,0);
   
   //Do the same for east west
   Vec3f EWPlaneN = northVector;
   EWPlaneN.normalize();
   Planef EastWestPlane(EWPlaneN,0);


   //Calculate the Azimuth
   float bodyDistanceFromNSPlane = NorthSouthPlane.distanceTo(bodyPoint);
   float bodyDistanceFromEWPlane = EastWestPlane.distanceTo(bodyPoint);

   Azimuth = atan(bodyDistanceFromNSPlane/bodyDistanceFromEWPlane);

   if (bodyDistanceFromNSPlane > 0) {
      if (bodyDistanceFromEWPlane > 0) {
         //do nothing
      }else{
         Azimuth += PI;
      }
   }else{
      if (bodyDistanceFromEWPlane > 0) {
         Azimuth = 2*PI + Azimuth;
      }else{
         Azimuth += PI;
      }
   }

   //Calculate the Altitude
   Altitude = asin(FloorPlane.distanceTo(bodyPoint)/Dist);

   //work out the distance from the observer to the body
   bodyVector -= observerVector;
   Dist = bodyVector.length();
   Dist -= body->getRadius();
}


Alex

[alexibu]

Hi,
Ive changed the function that I posted for calculating the Alt/Azimuth of a body from the surface of a planet.
It now works for any object not just other bodies in the solar system, and is double precision for the important parts. Changed body argument to a universal coord.
Alex

Code: Select all

void Body::getDirectionFromSurface(const Vec3f& lonLatAlt,const UniversalCoord & TargetPos,
                           const double & TargetRadius,const double when,
                               double& Dist,double& Azimuth,double&Altitude)
{
   //get the position of the observers sun in universal coords
   Point3f sunPos(0.0f, 0.0f, 0.0f);
   const Star* sun = getSystem()->getStar();
   if (sun != NULL) sunPos = sun->getPosition();

   //get the targets position in geographic coodinates
   Point3d bodyPoint = astro::heliocentricPosition(TargetPos,sunPos)
      *getGeographicToHeliocentric(when).inverse();

   //get the position of the observer in geographic coords
   Vec3f observerVectorf = planetocentricToCartesian(lonLatAlt);
   Vec3d observerVector(observerVectorf.x,observerVectorf.y,observerVectorf.z);

   //create a plane tangential to the earths surface where the
   //observer is
   Vec3d FloorN = observerVector;
   FloorN.normalize();
   Planed FloorPlane(FloorN,-observerVector.length());

   //find a point north of the observer by finding the floor planes
   //intersection with the geographical y axis
   Point3d northPoint = Planed::intersection(FloorPlane,
         Planed(Vec3d(1,0,0),0),Planed(Vec3d(0,0,1),0));
   
   //Do some conversions
   Vec3d northVector(northPoint.x,northPoint.y,northPoint.z);
   
   //calculate the vector from the observer to north
   northVector -= observerVector;

   //create a plane that splits the planet north south
   //from the point of view of the observer
   Vec3d NSPlaneN = northVector^observerVector;
   NSPlaneN.normalize();
   Planed NorthSouthPlane(NSPlaneN,0);
   
   //Do the same for east west
   Vec3d EWPlaneN = northVector;
   EWPlaneN.normalize();
   Planed EastWestPlane(EWPlaneN,0);

   //Calculate the Azimuth
   Azimuth = atan2(NorthSouthPlane.distanceTo(bodyPoint),
               EastWestPlane.distanceTo(bodyPoint));

   //work out the distance from the observer to the body
   Point3d Temp = bodyPoint - observerVector;
   Dist = Temp.distanceFromOrigin();

   //Calculate the Altitude
   Altitude = asin(FloorPlane.distanceTo(bodyPoint)/Dist);

   //Subtract the radius of the target
   Dist -= TargetRadius;
   
   //The distance to the planet the observer is on is sometimes slightly negative
   //due to precision errors - this looks bad so make sure that distance is positive
   Dist = fabs(Dist);
}



[Christophe]

Hi Alex,

I haven't tried your code yet, but I think it's a nice addition. Once 1.3.1 is released, I'll try and add it for display in the popup menu.

Something that would be really useful too would be Alt/Az cursor coordinates. Do you think you could do that as well?

[alexibu]

Err,
The code I posted doesn't work,
Ill post a new one soon
sorry
Alex

[alexibu]

Hi,
Ive now fixed this code.

Code: Select all

void Body::getDirectionFromSurface(const Vec3f& lonLatAlt,const UniversalCoord & TargetPos,
                           const double & TargetRadius,const double when,
                               double& Dist,double& Azimuth,double&Altitude)
{
   //get the position of the observers sun in universal coords
   Point3f sunPos(0.0f, 0.0f, 0.0f);
   const Star* sun = getSystem()->getStar();
   if (sun != NULL) sunPos = sun->getPosition();

   //get the targets position in geographic coodinates
   Point3d bodyPoint = astro::heliocentricPosition(TargetPos,sunPos)
      *getGeographicToHeliocentric(when).inverse();

   //get the position of the observer in geographic coords
   Vec3f observerVectorf = planetocentricToCartesian(lonLatAlt);
   Vec3d observerVector(observerVectorf.x,observerVectorf.y,observerVectorf.z);

   //create a plane tangential to the earths surface where the
   //observer is
   Vec3d FloorN = observerVector;
   FloorN.normalize();
   Planed FloorPlane(FloorN,-observerVector.length());

   //find a point north of the observer by finding the floor planes
   //intersection with the geographical y axis
   Point3d northPoint = Planed::intersection(FloorPlane,
         Planed(Vec3d(1,0,0),0),Planed(Vec3d(0,0,1),0));
   
   //Do some conversions
   Vec3d northVector(northPoint.x,northPoint.y,northPoint.z);
   
   //calculate the vector from the observer to north
   northVector -= observerVector;

   //create a plane that splits the planet north south
   //from the point of view of the observer
   Vec3d NSPlaneN = northVector^observerVector;
   NSPlaneN.normalize();
   Planed NorthSouthPlane(NSPlaneN,0);
   
   //Do the same for east west
   Vec3d EWPlaneN = northVector;
   EWPlaneN.normalize();
   Planed EastWestPlane(EWPlaneN,0);

   //Calculate the Azimuth
   Azimuth = atan2(NorthSouthPlane.distanceTo(bodyPoint),
               EastWestPlane.distanceTo(bodyPoint));

   //work out the distance from the observer to the body
   Point3d Temp = bodyPoint - observerVector;
   Dist = Temp.distanceFromOrigin();

   //Calculate the Altitude
   Altitude = asin(FloorPlane.distanceTo(bodyPoint)/Dist);

   //Subtract the radius of the target
   Dist -= TargetRadius;
   
   //The distance to the planet the observer is on is sometimes slightly negative
   //due to precision errors - this looks bad so make sure that distance is positive
   Dist = fabs(Dist);
}