mirror/astronomy
1
0
Fork 0
mirror of https://github.com/cosinekitty/astronomy.git synced 2026-05-19 14:27:52 -04:00
Code Issues Packages Projects Releases Wiki Activity

C RotationAxis: calculate north pole vector.

Browse Source 
Calculate the vector that points in the direction
of the body's north pole.
The unit test now checks for excessive angle
between the expected north pole vector and the
calculated north pole vector.
This commit is contained in:
Don Cross
2021-11-28 21:51:47 -05:00
parent ec35f21827
commit 20ff46bb27
5 changed files with 60 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files

39
generate/ctest.c
View File

@@ -4984,12 +4984,15 @@ fail:
/*-----------------------------------------------------------------------------------------------------------*/
static int AxisTestBody(astro_body_t body, const char *filename)
static int AxisTestBody(astro_body_t body, const char *filename, double arcmin_tolerance)
{
int error, lnum, nscanned, found_data;
int error, lnum, nscanned, found_data, count;
astro_axis_t axis;
astro_time_t time;
double jd, ra, dec;
astro_spherical_t sphere;
astro_vector_t north;
astro_angle_result_t diff;
double jd, ra, dec, arcmin, max_arcmin;
FILE *infile;
char line[100];
@@ -4997,8 +5000,10 @@ static int AxisTestBody(astro_body_t body, const char *filename)
if (infile == NULL)
FAIL("AxisTestBody: cannot open input file: %s\n", filename);
count = 0;
lnum = 0;
found_data = 0;
max_arcmin = 0.0;
while (ReadLine(line, sizeof(line), infile, filename, lnum))
{
++lnum;
@@ -5020,18 +5025,36 @@ static int AxisTestBody(astro_body_t body, const char *filename)
if (nscanned != 3)
FAIL("C AxisTestBody(%s line %d): could not scan data.\n", filename, lnum);
ra /= 15.0; /* convert degrees to sidereal hours */
time = Astronomy_TimeFromDays(jd - 2451545.0);
axis = Astronomy_RotationAxis(body, time);
if (axis.status != ASTRO_SUCCESS)
FAIL("C AxisTestBody(%s line %d): Astronomy_Axis returned error %d\n", filename, lnum, axis.status);
/* Convert the reference angles to a reference north pole vector. */
sphere.status = ASTRO_SUCCESS;
sphere.dist = 1.0;
sphere.lat = dec;
sphere.lon = ra; /* tricky: RA is in degrees, not sidereal hours */
north = Astronomy_VectorFromSphere(sphere, time);
if (north.status != ASTRO_SUCCESS)
FAIL("C AxisTestBody(%s line %d): VectorFromSphere error %d\n", filename, lnum, north.status);
/* Find angle between two versions of the north pole. Use that as the measure of error. */
DEBUG("C AxisTestBody(%s): correct (ra=%lf, dec=%lf) calc(ra=%lf, dec=%lf)\n", filename, ra, dec, axis.ra, axis.dec);
diff = Astronomy_AngleBetween(north, axis.zdir);
if (diff.status != ASTRO_SUCCESS)
FAIL("C AxisTestBody(%s line %d): AngleBetween error %d\n", filename, lnum, diff.status);
arcmin = diff.angle * 60.0; /* convert error degrees to arcminutes */
if (arcmin > max_arcmin)
max_arcmin = arcmin;
++count;
}
}
DEBUG("C AxisTestBody(%s): %d test cases, max arcmin error = %0.6lf.\n", filename, count, max_arcmin);
if (max_arcmin > arcmin_tolerance)
FAIL("C AxisTestBody(%s): EXCESSIVE ERROR = %lf arcmin\n", filename, max_arcmin);
error = 0;
fail:
if (infile != NULL) fclose(infile);
@@ -5041,8 +5064,8 @@ fail:
static int AxisTest(void)
{
int error;
CHECK(AxisTestBody(BODY_SUN, "axis/Sun.txt"));
CHECK(AxisTestBody(BODY_MERCURY, "axis/Mercury.txt"));
CHECK(AxisTestBody(BODY_SUN, "axis/Sun.txt", 0.000000));
CHECK(AxisTestBody(BODY_MERCURY, "axis/Mercury.txt", 0.074340));
printf("C AxisBody: PASS\n");
fail:
return error;

15
generate/template/astronomy.c
View File

@@ -8738,12 +8738,13 @@ astro_axis_t Astronomy_RotationAxis(astro_body_t body, astro_time_t time)
double d = time.tt;
double T = d / 36525.0;
double M1, M2, M3, M4, M5;
double radlat, radlon, rcoslat;
switch (body)
{
case BODY_SUN:
ra = 286.13;
dec = 68.87;
dec = 63.87;
w = 84.176 + (14.1844 * d);
break;
@@ -8773,8 +8774,18 @@ astro_axis_t Astronomy_RotationAxis(astro_body_t body, astro_time_t time)
axis.ra = ra / 15.0; /* convert degrees to sidereal hours */
axis.dec = dec;
axis.spin = w;
axis.status = ASTRO_SUCCESS;
/* calculate the north pole vector using the given angles. */
radlat = dec * DEG2RAD;
radlon = ra * DEG2RAD;
rcoslat = cos(radlat);
axis.zdir.x = rcoslat * cos(radlon);
axis.zdir.y = rcoslat * sin(radlon);
axis.zdir.z = sin(radlat);
axis.zdir.t = time;
axis.zdir.status = ASTRO_SUCCESS;
axis.status = ASTRO_SUCCESS;
return axis;
}

1
source/c/README.md
View File

@@ -3220,6 +3220,7 @@ The fields `ra`, `dec`, and `spin` correspond to the variables α0, δ0, and W,
| `double` | `ra` | The J2000 right ascension of the body's north pole direction, in sidereal hours. |
| `double` | `dec` | The J2000 declination of the body's north pole direction, in degrees. |
| `double` | `spin` | Rotation angle of the body's prime meridian, in degrees. |
| [`astro_vector_t`](#astro_vector_t) | `zdir` | A J2000 dimensionless unit vector pointing in the direction of the body's north pole. |
---

15
source/c/astronomy.c
View File

@@ -10432,12 +10432,13 @@ astro_axis_t Astronomy_RotationAxis(astro_body_t body, astro_time_t time)
double d = time.tt;
double T = d / 36525.0;
double M1, M2, M3, M4, M5;
double radlat, radlon, rcoslat;
switch (body)
{
case BODY_SUN:
ra = 286.13;
dec = 68.87;
dec = 63.87;
w = 84.176 + (14.1844 * d);
break;
@@ -10467,8 +10468,18 @@ astro_axis_t Astronomy_RotationAxis(astro_body_t body, astro_time_t time)
axis.ra = ra / 15.0; /* convert degrees to sidereal hours */
axis.dec = dec;
axis.spin = w;
axis.status = ASTRO_SUCCESS;
/* calculate the north pole vector using the given angles. */
radlat = dec * DEG2RAD;
radlon = ra * DEG2RAD;
rcoslat = cos(radlat);
axis.zdir.x = rcoslat * cos(radlon);
axis.zdir.y = rcoslat * sin(radlon);
axis.zdir.z = sin(radlat);
axis.zdir.t = time;
axis.zdir.status = ASTRO_SUCCESS;
axis.status = ASTRO_SUCCESS;
return axis;
}

2
source/c/astronomy.h
View File

@@ -935,6 +935,8 @@ typedef struct
double ra; /**< The J2000 right ascension of the body's north pole direction, in sidereal hours. */
double dec; /**< The J2000 declination of the body's north pole direction, in degrees. */
double spin; /**< Rotation angle of the body's prime meridian, in degrees. */
astro_vector_t zdir; /**< A J2000 dimensionless unit vector pointing in the direction of the body's north pole. */
}
astro_axis_t;