I'm following @whuber's post for calculating azimuth and elevation angles based on two vector3's. The elevation equation works great and gives me numbers that makes sense. I'm running into issues with the azimuth equation though and its way above my head on how it fully works. My end goal is to calculate az/el angles from a point on earth to a satellite. Also, while this post includes JavaScript, I'm not looking for coding help but how to solve the problem mathematically.
I've created this little 3D simulation for calculating these to try and wrap my head around it. I'm not looking for a perfect simulation so the fact that the earth is perfectly round is okay. Also the "earth" has a radius of 1 and the satellite a radius of 2 just to make the math simpler.
https://jsfiddle.net/51yerkvu/16/
By changing the longitude slider for either the base station or satellite you can see the El angle properly calculates. As it goes below the horizon, the numbers become negative and directly above its 90°.
The az calculation seems to be correct in some scenarios, but not others. For example, the default position when loading the jsfiddle it calculates an azimuth of -63.43 degrees. Moving its longitude -1° causes it to calculate -90 degrees, and moving +1° causes it to calculate +90 degrees. The azimuth angle should be the same when changing the longitude as its on the same plane as the satellite and is just moving closer or further away. Moving its latitude causes it to calculate really wrong numbers.
Cos(azimuth) = (-z*x*dx - z*y*dy + (x^2+y^2)*dz) / Sqrt((x^2+y^2)(x^2+y^2+z^2)(dx^2+dy^2+dz^2))
Sin(azimuth) = (-y*dx + x*dy) / Sqrt((x^2+y^2)(dx^2+dy^2+dz^2))
By plugging in whubers example points (1285410, -4797210, 3994830) and (1202990, -4824940, 3999870) into my program, I do get roughly the same value as him. He calculated .331 for his azimuth where I'm getting .336. I found this to be a discrepancy in the Cos(azimuth) function providing a slightly different number. I don't think this is the issue though and I've rechecked that I put the expression in correctly a dozen times.
Maybe I'm misunderstanding and this approach and this isn't a good way to calculate azimuth for what I'm trying to do. Elevation works perfect, and I've tried a few other approaches leveraging the 3d library I'm using but haven't been able to get something reliable.
Converting whubers equation to JavaScript I get this:
function calculateAzimuth(p1, p2) {
const d = p2.clone().sub(p1.clone());
// Cos(azimuth) = (-z*x*dx - z*y*dy + (x^2+y^2)*dz) / Sqrt((x^2+y^2)(x^2+y^2+z^2)(dx^2+dy^2+dz^2))
let azCos = (-p1.z*p1.x*d.x - p1.z*p1.y*d.y + (p1.x**2+p1.y**2)*d.z);
azCos /= Math.sqrt((p1.x**2+p1.y**2) * (p1.x**2+p1.y**2+p1.z**2) * (d.x**2+d.y**2+d.z**2));
// Sin(azimuth) = (-y*dx + x*dy) / Sqrt((x^2+y^2)(dx^2+dy^2+dz^2))
let azSin = (-p1.y*d.x + p1.x*d.y);
azSin /= Math.sqrt((p1.x**2+p1.y**2) * (d.x**2+d.y**2+d.z**2));
return Math.atan(azCos / azSin);
}
I've checked it over a dozen times to make sure its the same, but I'm hoping I just have an equation error. I'm also open to trying other ways of calculating azimuth, but I would like my solution to be based on positions relative to the earths center.
If I look down the blue axis of the point on the earth and project the satellite onto the same plane, its then just a simple 2D problem. I'm not sure how to go about that mathematically however.
