Я не знаю, является ли это простейшим решением, и предполагается, что мир - это сфера.
Определение:
R - радиус сферы (то есть земли).
r - радиус круга (в тех же единицах).
t - это угол, образованный дугой большого круга длины r в центре сферы, поэтому t = r / R радиан.
Теперь предположим, что сфера имеет радиус 1 и центрирована в начале координат.
C - единичный вектор, представляющий центр круга.
Представьте себе круг вокруг северного полюса и рассмотрите точку, где плоскость круга пересекает линию от центра земли к северному полюсу. Очевидно, что эта точка будет где-то ниже северного полюса.
K - это соответствующая точка "ниже" C (то есть, где плоскость вашего круга пересекает C), поэтому K = cos (t) C
s - радиус окружности, измеренный в трехмерном пространстве (т.е. не на сфере), поэтому s = sin (t)
Теперь нам нужны точки на окружности в трехмерном пространстве с центром K, радиусом s и лежащим в плоскости, проходящей и перпендикулярной к K.
В этом ответе (игнорируя вращение) объясняется, как найти базисный вектор для плоскости (то есть вектор, ортогональный к нормали K или C). Используйте перекрестное произведение, чтобы найти секунду.
Назовите эти базисные векторы U и V.
// Pseudo-code to calculate 20 points on the circle
for (a = 0; a != 360; a += 18)
{
// A point on the circle and the unit sphere
P = K + s * (U * sin(a) + V * cos(a))
}
Преобразуйте каждую точку в сферические координаты, и все готово.
Мне было скучно, я написал это на C #. Результаты правдоподобны: они находятся по кругу и лежат на сфере. Большая часть кода реализует struct, представляющий вектор. Фактический расчет очень прост.
using System;
namespace gpsCircle
{
struct Gps
{
// In degrees
public readonly double Latitude;
public readonly double Longtitude;
public Gps(double latitude, double longtitude)
{
Latitude = latitude;
Longtitude = longtitude;
}
public override string ToString()
{
return string.Format("({0},{1})", Latitude, Longtitude);
}
public Vector ToUnitVector()
{
double lat = Latitude / 180 * Math.PI;
double lng = Longtitude / 180 * Math.PI;
// Z is North
// X points at the Greenwich meridian
return new Vector(Math.Cos(lng) * Math.Cos(lat), Math.Sin(lng) * Math.Cos(lat), Math.Sin(lat));
}
}
struct Vector
{
public readonly double X;
public readonly double Y;
public readonly double Z;
public Vector(double x, double y, double z)
{
X = x;
Y = y;
Z = z;
}
public double MagnitudeSquared()
{
return X * X + Y * Y + Z * Z;
}
public double Magnitude()
{
return Math.Sqrt(MagnitudeSquared());
}
public Vector ToUnit()
{
double m = Magnitude();
return new Vector(X / m, Y / m, Z / m);
}
public Gps ToGps()
{
Vector unit = ToUnit();
// Rounding errors
double z = unit.Z;
if (z > 1)
z = 1;
double lat = Math.Asin(z);
double lng = Math.Atan2(unit.Y, unit.X);
return new Gps(lat * 180 / Math.PI, lng * 180 / Math.PI);
}
public static Vector operator*(double m, Vector v)
{
return new Vector(m * v.X, m * v.Y, m * v.Z);
}
public static Vector operator-(Vector a, Vector b)
{
return new Vector(a.X - b.X, a.Y - b.Y, a.Z - b.Z);
}
public static Vector operator+(Vector a, Vector b)
{
return new Vector(a.X + b.X, a.Y + b.Y, a.Z + b.Z);
}
public override string ToString()
{
return string.Format("({0},{1},{2})", X, Y, Z);
}
public double Dot(Vector that)
{
return X * that.X + Y * that.Y + Z * that.Z;
}
public Vector Cross(Vector that)
{
return new Vector(Y * that.Z - Z * that.Y, Z * that.X - X * that.Z, X * that.Y - Y * that.X);
}
// Pick a random orthogonal vector
public Vector Orthogonal()
{
double minNormal = Math.Abs(X);
int minIndex = 0;
if (Math.Abs(Y) < minNormal)
{
minNormal = Math.Abs(Y);
minIndex = 1;
}
if (Math.Abs(Z) < minNormal)
{
minNormal = Math.Abs(Z);
minIndex = 2;
}
Vector B;
switch (minIndex)
{
case 0:
B = new Vector(1, 0, 0);
break;
case 1:
B = new Vector(0, 1, 0);
break;
default:
B = new Vector(0, 0, 1);
break;
}
return (B - minNormal * this).ToUnit();
}
}
class Program
{
static void Main(string[] args)
{
// Phnom Penh
Gps centre = new Gps(11.55, 104.916667);
// In metres
double worldRadius = 6371000;
// In metres
double circleRadius = 1000;
// Points representing circle of radius circleRadius round centre.
Gps[] points = new Gps[20];
CirclePoints(points, centre, worldRadius, circleRadius);
}
static void CirclePoints(Gps[] points, Gps centre, double R, double r)
{
int count = points.Length;
Vector C = centre.ToUnitVector();
double t = r / R;
Vector K = Math.Cos(t) * C;
double s = Math.Sin(t);
Vector U = K.Orthogonal();
Vector V = K.Cross(U);
// Improve orthogonality
U = K.Cross(V);
for (int point = 0; point != count; ++point)
{
double a = 2 * Math.PI * point / count;
Vector P = K + s * (Math.Sin(a) * U + Math.Cos(a) * V);
points[point] = P.ToGps();
}
}
}
}