If you want a stable method of computing geodesic distances, I recommend Richie Carmichael's wrapper for ESRI's Projection Engine.
Update: I just tried Richie's code with ArcGIS 10.0 on Vista64 and get an exception after calling LoadLibrary. I'll look into that more later.
For now though, here is some code in response to questions in the comments of another answer.
The code compares IProximityOperator for points with and without spatial references. Then it shows how to use an azimuthal equidistant projection (with first point being the point of tangency) to find the great circle distance.
private void Test()
{
IPoint p1 = new PointClass();
p1.PutCoords(-98.0, 28.0);
IPoint p2 = new PointClass();
p2.PutCoords(-78.0, 28.0);
Debug.Print("Euclidian Distance {0}", EuclidianDistance(p1, p2));
Debug.Print("Distance with no spatialref {0}", GetDistance(p1, p2));
ISpatialReferenceFactory srf = new SpatialReferenceEnvironmentClass();
IGeographicCoordinateSystem gcs =
srf.CreateGeographicCoordinateSystem((int)esriSRGeoCSType.esriSRGeoCS_WGS1984);
p1.SpatialReference = gcs;
p2.SpatialReference = gcs;
Debug.Print("Distance with spatialref {0}", GetDistance(p1, p2));
Debug.Print("Great Circle Distance {0}", GreatCircleDist(p1, p2));
}
private double GetDistance(IPoint p1, IPoint p2)
{
return ((IProximityOperator)p1).ReturnDistance(p2);
}
private double EuclidianDistance(IPoint p1, IPoint p2)
{
return Math.Sqrt(Math.Pow((p2.X - p1.X),2.0) + Math.Pow((p2.Y - p1.Y), 2.0));
}
private double GreatCircleDist(IPoint p1, IPoint p2)
{
ISpatialReferenceFactory srf = new SpatialReferenceEnvironmentClass();
IProjectedCoordinateSystem pcs =
srf.CreateProjectedCoordinateSystem((int)esriSRProjCSType.esriSRProjCS_WGS1984N_PoleAziEqui);
pcs.set_CentralMeridian(true, p1.X);
((IProjectedCoordinateSystem2)pcs).LatitudeOfOrigin = p1.Y;
p1.SpatialReference = pcs.GeographicCoordinateSystem;
p1.Project(pcs);
p2.SpatialReference = pcs.GeographicCoordinateSystem;
p2.Project(pcs);
return EuclidianDistance(p1, p2);
}
Here's the output:
Euclidian Distance 20
Distance with no spatialref 20
Distance with spatialref 20
Great Circle Distance 1965015.61318737
I think it would be interesting to test this against the projection engine dll (pe.dll). Will post results if I ever get Richie's code to work.
Update:
Once I changed Richies code to compile for x86, I got it to run. Interesting ... the great circle distance it give me is 1960273.80162999 - a significant difference from that returned from the azimuthal equidistant method above.
Given that your inputs are in degrees, trigonometric operations are unavoidable. At a minimum you need to account for the distortion in latitude relative to longitude at all points not on the equator.
The simplest approximation (based on a spherical model of the earth and suitable only for small distances) I can think of projects the two points (lat1, lon1), (lat2, lon2) to the points
(x1, y1) = (lon1 * cos(lat1), lat1) and
(x2, y2) = (lon2 * cos(lat2), lat2)
in the Euclidean plane. (This is a Sinusoidal projection.) The direction vector from the viewpoint (#1) to the target point (#2) therefore is
(u, v) = (x2-x1, y2-y1)
with length n = Sqrt(u^2 + v^2), whence the unit direction vector equals
(u/n, v/n).
If you can specify the direction you're facing as a vector, rather than as an angle, you can avoid some more trig. Otherwise, a direction t degrees east of north has to be converted to the unit direction vector
(a, b) = (sin(t), cos(t)).
Finally, visibility is tested by comparing the inner product
(u/n, v/n) . (a, b) = (a*u + v*b)/n
to the cosine of half the horizontal field of view (another trig value, but it can be precomputed once and for all). The inner product is less than this cosine only for the invisible points.
This method requires computing one cosine for each base point and target point, as well as a sine and a cosine for the direction of view.
For example, with a HFOV of 30 degrees you can look 15 degrees to the right and left. The cosine of 15 degrees equals 0.965926. If you are standing at lon = -1.080 and lat = 51.40 and the target point is located at lon = -1.085 and lat = 51.20, then
(x1, y1) = (-1.080 / cos(51.4), 51.4) = (-0.67397, 51.4)
(x2, y2) = (-1.085 / cos(51.2), 51.2) = (-0.67987, 51.2)
(u, v) = (-0.00608, -0.200)
n = 0.200092
(u/n, v/n) = (-0.03036, -0.99954).
The direction vector for due north is (sin(0), cos(0)) = (0, 1). Its inner product with (u/n, v/n) equals -0.99954. Because this is (a lot) less than 0.965926, the point is invisible. (This, I hope, was obvious at the outset, because the target point is south of--that is, behind--the point of view.) If you happened to be looking due south, your direction vector would be (sin(180), cos(180)) = (0, -1), the inner product would equal +0.99954, and because this exceeds 0.965926, you would conclude that the target point is visible.
In some cases you can avoid a lot of trig. If you have a single base point and lots of target points nearby, all the cosines of the latitudes will be approximately the same and equal to the cosine of the base point's latitude. In this case you need compute it only once and use that approximation for all the other calculations. In this fashion, no matter how many target points there are, you need only four trig evaluations for each base point, HFOV, and direction of view combination. Of course, because this is an approximation (albeit a good one), it would be silly to worry about "edge cases" (which I presume are points right at the edge of the field of view).
Best Answer
You appear to be describing rotations of the Riemann Sphere: they are conformal and preserve the spherical metric. There are several simple ways to write them down (not involving any trigonometry!). A good one is to consider coordinates (x,y) in the plane as complex numbers z = x + yI, I^2 = -1. Given any four complex numbers a, b, c, and d, a fractional linear transformation of z is the value (az + b)/(cz + d). Thinking of z as the image of a point on the sphere via stereographic projection, it's easy to work out that for a fractional linear transformation to be a true rotation, a and d must be complex conjugates of each other and b and -c must also be mutual complex conjugates. (This defines the matrix group SU(2,C).)
Let Q be the point (in projected coordinates) you would like to base a new projection on. All we need to do is rotate Q to infinity, because infinity corresponds to the projection's base point (the north pole). This implies the denominator cz + d must equal zero, giving a solution
In other words, the change of coordinates you seek has the formula
As an example, suppose you want to make Q = (5, 5) = 5 + 5I the new origin of projection. Using the rules of complex arithmetic we can work out the transformation in terms of the coordinates (x,y):
That is, writing the new coordinates of (x,y) as (x',y'), we have
and
You can follow this up by any rotation around the origin to reorient the new map. A good choice is to multiply the result by Q/Conjugate(Q) = Q^2 / |Q|^2. This is because the transformation has a nice interpretation:
The first factor merely translates all points by the (small) amount 1/Conjugate(Q). In particular, this keeps the map oriented correctly. The second factor can be ignored when z/Q is very small. Typically, for small rotations, Q is already "near infinity": that is, it is large compared with any point on the map, because it is the projection a point close to the original north pole. This justifies approximating the change in projection by means of a translation and, in addition, it tells us how to measure the error: the size of z/Q (the first, and largest, neglected term in the series on the right) is the ratio of the sizes of z and Q (that is, the ratio of their distances from the map origin). In other words, when the original projection of Q is way beyond the extent of your map, you will likely be ok with the approximate formula.