mirror of
https://github.com/cosinekitty/astronomy.git
synced 2025-12-27 01:30:50 -05:00
33de407723
See the following discussion for context: https://github.com/cosinekitty/astronomy/discussions/308 Added a demo program that shows how to search for when a body enters a window defined in terms of an observer's horizontal frame of reference, given a range of altitudes and a range of azimuths.
263 lines
8.1 KiB
C++
263 lines
8.1 KiB
C++
/*
|
|
altazsearch.cpp - Don Cross - 2023-06-17
|
|
|
|
https://github.com/cosinekitty/astronomy/discussions/308
|
|
|
|
Problem: given a range of altitudes and azimuths that
|
|
form a "window" on the sky, search for when the Moon
|
|
enters/exits that window starting from a given search time.
|
|
*/
|
|
|
|
#include <cmath>
|
|
#include <cstdio>
|
|
#include <string>
|
|
#include <stdexcept>
|
|
#include "astronomy.h"
|
|
|
|
|
|
struct Event
|
|
{
|
|
astro_time_t time;
|
|
double azimuth = NAN;
|
|
double altitude = NAN;
|
|
|
|
Event()
|
|
{
|
|
time.tt = time.ut = time.eps = time.psi = time.st = NAN;
|
|
}
|
|
|
|
void Print() const
|
|
{
|
|
char text[TIME_TEXT_BYTES];
|
|
Astronomy_FormatTime(time, TIME_FORMAT_SECOND, text, sizeof(text));
|
|
printf("%s az=%0.2lf alt=%0.2lf", text, azimuth, altitude);
|
|
}
|
|
};
|
|
|
|
|
|
struct Solution
|
|
{
|
|
bool valid = false;
|
|
Event start;
|
|
Event finish;
|
|
|
|
void Print() const
|
|
{
|
|
if (valid)
|
|
{
|
|
printf("Start: ");
|
|
start.Print();
|
|
printf("; Finish: ");
|
|
finish.Print();
|
|
printf(".\n");
|
|
}
|
|
else
|
|
{
|
|
printf("No solution.\n");
|
|
}
|
|
}
|
|
};
|
|
|
|
|
|
void Verify(astro_status_t status, const char *message)
|
|
{
|
|
if (status != ASTRO_SUCCESS)
|
|
throw std::logic_error(std::string(message) + ": error " + std::to_string(static_cast<int>(status)));
|
|
}
|
|
|
|
|
|
class SearchProblem
|
|
{
|
|
private:
|
|
astro_body_t body;
|
|
astro_observer_t observer;
|
|
double az1;
|
|
double az2;
|
|
double alt1;
|
|
double alt2;
|
|
astro_vector_t center;
|
|
|
|
astro_horizon_t Position(astro_time_t time) const
|
|
{
|
|
// Get topocentric equatorial coordinates of body, using the Earth's equator of date.
|
|
astro_equatorial_t equ = Astronomy_Equator(body, &time, observer, EQUATOR_OF_DATE, ABERRATION);
|
|
Verify(equ.status, "Equator");
|
|
|
|
// Convert to observer's horizontal coordinates, correcting for atmospheric refraction.
|
|
return Astronomy_Horizon(&time, observer, equ.ra, equ.dec, REFRACTION_NORMAL);
|
|
}
|
|
|
|
double AngularDistance(astro_time_t time) const
|
|
{
|
|
astro_horizon_t hor = Position(time);
|
|
|
|
// Translate angular horizontal coordinates to a vector.
|
|
// Do NOT remove the refraction from the vector (REFRACTION_NONE).
|
|
astro_spherical_t sphere;
|
|
sphere.status = ASTRO_SUCCESS;
|
|
sphere.dist = 1.0;
|
|
sphere.lat = hor.altitude;
|
|
sphere.lon = hor.azimuth;
|
|
astro_vector_t vec = Astronomy_VectorFromHorizon(sphere, time, REFRACTION_NONE);
|
|
Verify(vec.status, "VectorFromHorizon");
|
|
|
|
// Calculate the angle in degrees between the body and the center of the window.
|
|
astro_angle_result_t result = Astronomy_AngleBetween(center, vec);
|
|
Verify(result.status, "AngleBetween");
|
|
|
|
return result.angle;
|
|
}
|
|
|
|
static astro_func_result_t DistanceSlopeCallback(void *context, astro_time_t time)
|
|
{
|
|
const SearchProblem& p = *static_cast<const SearchProblem *>(context);
|
|
astro_func_result_t result;
|
|
result.value = p.DistanceSlope(time);
|
|
result.status = ASTRO_SUCCESS;
|
|
return result;
|
|
}
|
|
|
|
double DistanceSlope(astro_time_t time) const
|
|
{
|
|
const double dt = 0.1 / 86400.0;
|
|
astro_time_t t1 = Astronomy_AddDays(time, -dt);
|
|
astro_time_t t2 = Astronomy_AddDays(time, +dt);
|
|
double a1 = AngularDistance(t1);
|
|
double a2 = AngularDistance(t2);
|
|
return (a2 - a1) / (2 * dt);
|
|
}
|
|
|
|
bool IsInsideWindow(astro_time_t time) const
|
|
{
|
|
astro_horizon_t hor = Position(time);
|
|
bool insideAzimuthLimits;
|
|
if (az1 <= az2)
|
|
insideAzimuthLimits = (az1 <= hor.azimuth && hor.azimuth <= az2);
|
|
else
|
|
insideAzimuthLimits = (az1 <= hor.azimuth || hor.azimuth <= az2);
|
|
|
|
return insideAzimuthLimits && alt1 <= hor.altitude && hor.altitude <= alt2;
|
|
}
|
|
|
|
Solution FindBracket(astro_time_t closestTime) const
|
|
{
|
|
Solution solution;
|
|
|
|
// If the closestTime is inside the window, we can find a bracket.
|
|
// Otherwise, there is no bracket here.
|
|
if (IsInsideWindow(closestTime))
|
|
{
|
|
const double dt = 10.0 / (24.0 * 60.0); // 10 minutes, converted to days
|
|
|
|
// Look backwards until we find a time before entering the window.
|
|
// Do a binary search to find the moment when we enter the window.
|
|
solution.start = FindTransition(closestTime, -dt);
|
|
|
|
// Look forward until we find a time after leaving the window.
|
|
// Do a binary search to find the moment we leave the window.
|
|
solution.finish = FindTransition(closestTime, +dt);
|
|
|
|
solution.valid = true;
|
|
}
|
|
|
|
return solution;
|
|
}
|
|
|
|
Event FindTransition(astro_time_t closestTime, double dt) const
|
|
{
|
|
// Find a bracket [t1, t2] that straddles being inside/outside the window.
|
|
astro_time_t t1 = closestTime;
|
|
astro_time_t t2 = Astronomy_AddDays(closestTime, dt);
|
|
while (IsInsideWindow(t2))
|
|
{
|
|
t1 = t2;
|
|
t2 = Astronomy_AddDays(t2, dt);
|
|
}
|
|
|
|
// Do a binary search to find the moment of transition, within tolerance.
|
|
const double tolerance = 0.1 / (3600.0 * 24.0); // one tenth of a second, expressed in days
|
|
astro_time_t tm = Astronomy_TimeFromDays((t1.ut + t2.ut) / 2);
|
|
while (fabs(t2.ut - t1.ut) > tolerance)
|
|
{
|
|
if (IsInsideWindow(tm))
|
|
t1 = tm;
|
|
else
|
|
t2 = tm;
|
|
tm = Astronomy_TimeFromDays((t1.ut + t2.ut) / 2);
|
|
}
|
|
|
|
astro_horizon_t hor = Position(tm);
|
|
|
|
Event event;
|
|
event.time = tm;
|
|
event.altitude = hor.altitude;
|
|
event.azimuth = hor.azimuth;
|
|
return event;
|
|
}
|
|
|
|
public:
|
|
SearchProblem(astro_body_t _body, astro_observer_t _observer, double _az1, double _az2, double _alt1, double _alt2)
|
|
: body(_body)
|
|
, observer(_observer)
|
|
, az1(_az1)
|
|
, az2(_az2)
|
|
, alt1(_alt1)
|
|
, alt2(_alt2)
|
|
{
|
|
astro_time_t dummyTime = Astronomy_TimeFromDays(0.0);
|
|
astro_spherical_t sphere;
|
|
sphere.status = ASTRO_SUCCESS;
|
|
sphere.dist = 1.0;
|
|
sphere.lat = (_alt1 + _alt2) / 2;
|
|
sphere.lon = (_az1 + _az2) / 2;
|
|
center = Astronomy_VectorFromHorizon(sphere, dummyTime, REFRACTION_NONE);
|
|
}
|
|
|
|
Solution FindNext(astro_time_t startTime, double limitDays)
|
|
{
|
|
astro_time_t stopTime = Astronomy_AddDays(startTime, limitDays);
|
|
const double stepDays = 1.0 / 24.0; // one hour
|
|
astro_time_t t1 = startTime;
|
|
double m1 = DistanceSlope(t1);
|
|
|
|
while (t1.ut < stopTime.ut)
|
|
{
|
|
astro_time_t t2 = Astronomy_AddDays(t1, stepDays);
|
|
double m2 = DistanceSlope(t2);
|
|
if (m1 <= 0.0 && m2 >= 0.0)
|
|
{
|
|
astro_search_result_t result = Astronomy_Search(DistanceSlopeCallback, this, t1, t2, 0.1);
|
|
if (result.status == ASTRO_SUCCESS)
|
|
{
|
|
// We found a time bracket [t1, t2] where the body passes closest to
|
|
// the center of the target window. Now search nearby for when the
|
|
// body enters and exits the window.
|
|
Solution solution = FindBracket(result.time);
|
|
if (solution.valid)
|
|
return solution;
|
|
}
|
|
}
|
|
t1 = t2;
|
|
m1 = m2;
|
|
}
|
|
|
|
// We could not find a solution.
|
|
// Return a default-constructed Solution, which indicates failure.
|
|
return Solution();
|
|
}
|
|
};
|
|
|
|
|
|
int main(int argc, const char *argv[])
|
|
{
|
|
// 2023-06-17T12:00:00Z Moon AZ = 74.03 ALT = 27.44
|
|
// 2023-06-17T13:00:00Z Moon AZ = 79.22 ALT = 39.74
|
|
astro_observer_t observer = Astronomy_MakeObserver(30.0, -80.0, 0.0);
|
|
SearchProblem problem(BODY_MOON, observer, 74.0, 78.0, 25.0, 40.0);
|
|
astro_time_t startTime = Astronomy_MakeTime(2023, 6, 17, 0, 0, 0.0);
|
|
Solution solution = problem.FindNext(startTime, 2.0);
|
|
solution.Print();
|
|
return 0;
|
|
}
|
|
|