Spheric trigonometry - large circle on earth (sphere)

Question

I'm trying to draw a very large circle on the earth surface and my formula fails either when the circle goes above the north pole or when some component reaches 90 degree length (I'm not sure which one). The circle should represent the footprint of visible satellites in a given height. Works fine with a height up to 1000km, but for larger values the upper part of the circle flips. Instead of going to the other side of the world it stays on the same side.

I assume the polar triangle exceeds its definition, but I'm unable to find a solution.

Here is a test code. It calculates the values, and creates a KML which can be used to show it on a map (there are multiple web services and programs available to do so). With a satheight of 20000km it is extremely obvious.

I'd appreciate any hints how to fix this issue. Some simple case handling or do I need a totally different formula?

#!/usr/bin/perl -w
use strict; # enforce variable declaration
use Math::Complex qw(pi);
use POSIX;  # ceil function

my $rhodeg = 180.0 / pi;
my $earthradius = 6365000.0; #[m]
my $reflat = 52;
my $reflon = 13;
#my $satheight = 1300 * 1000.0;
my $satheight = 5000 * 1000.0;

open FILE,">","out.kml" or die;

print FILE "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n"
. "<kml xmlns=\"http://www.opengis.net/kml/2.2\">\n"
. "<Document>\n";

foreach my $elev (0)
{
  my $arc = -asin($earthradius * cos($elev/$rhodeg)/($earthradius+$satheight)) + (90-$elev)/$rhodeg;

  print FILE "<Placemark>\n<LineString>\n<coordinates>\n";
  for(my $az = 0; $az <= 360; $az += 1)
  {
    my $a = sin($reflat/$rhodeg);
    my $b = cos($arc);
    my $c = cos($reflat/$rhodeg);
    my $d = sin($arc);
    my $e = cos($az/$rhodeg);
    my $f = sin($az/$rhodeg);
    my $g = $reflon/$rhodeg;

    my $lat = asin($a*$b+$c*$d*$e);
    my $lon = $g+asin($d*$f/cos($lat));
    $lat *= $rhodeg;
    $lon *= $rhodeg;
    printf("Elevation $elev Azimuth $az Latitude $lat Longitude $lon : a %f b %f c %f d %f e %f f %f g %f ab %f cde %f\n", $a, $b, $c, $d, $e, $f, $g, $a*$b, $c*$d*$e);

    if($lat > 90)
    {
      $lat = 180 - $lat;
      $lon = 180 + $lon;

      $lon -= 360 if($lon >= 180);
      $lon += 360 if($lon <= -180);
    }
    print FILE "$lon,$lat,0\n";
  }
  print FILE "</coordinates>\n</LineString>\n</Placemark>\n";
}
print FILE "</Document>\n"
. "</kml>\n";

close FILE;

P.S. It is clear, that the circle is not really a visible circle with most projections, except for very small height values. It will always be a somewhat distorted thing :-)

Answer 1

Using another method to calculate the longitude works as expected:

my $lon = ($g + ($d*$f/cos($lat) <= 0.0)?-1:1) * acos(($b*$c-$d*$a*$e)/cos($lat)));