mirror/astronomy
1
0
Fork 0
mirror of https://github.com/cosinekitty/astronomy.git synced 2026-05-25 01:09:09 -04:00
Code Issues Packages Projects Releases Wiki Activity

C: Eliminated duplicate code - nutation() is based on nutation_rot().

Browse Source 
Instead of having the same calculations duplicated in both
nutation() and nutation_rot(), I reworked nutation() in terms of
nutation_rot(). Use nutation_rot() to calculate the rotation matrix,
then multiply that matrix by the input vector to produce the output vector.
This commit is contained in:
Don Cross
2019-12-14 18:15:13 -05:00
parent c4a537a236
commit e7d8804bf1
2 changed files with 16 additions and 76 deletions
Download Patch File Download Diff File Expand all files Collapse all files

46
generate/template/astronomy.c
View File

@@ -942,44 +942,6 @@ static astro_equatorial_t vector2radec(const double pos[3])
return equ;
}
static void nutation(astro_time_t *time, int direction, const double inpos[3], double outpos[3])
{
earth_tilt_t tilt = e_tilt(time);
double oblm = tilt.mobl * DEG2RAD;
double oblt = tilt.tobl * DEG2RAD;
double psi = tilt.dpsi * ASEC2RAD;
double cobm = cos(oblm);
double sobm = sin(oblm);
double cobt = cos(oblt);
double sobt = sin(oblt);
double cpsi = cos(psi);
double spsi = sin(psi);
double xx = cpsi;
double yx = -spsi * cobm;
double zx = -spsi * sobm;
double xy = spsi * cobt;
double yy = cpsi * cobm * cobt + sobm * sobt;
double zy = cpsi * sobm * cobt - cobm * sobt;
double xz = spsi * sobt;
double yz = cpsi * cobm * sobt - sobm * cobt;
double zz = cpsi * sobm * sobt + cobm * cobt;
if (direction == 0)
{
/* forward rotation */
outpos[0] = xx * inpos[0] + yx * inpos[1] + zx * inpos[2];
outpos[1] = xy * inpos[0] + yy * inpos[1] + zy * inpos[2];
outpos[2] = xz * inpos[0] + yz * inpos[1] + zz * inpos[2];
}
else
{
/* inverse rotation */
outpos[0] = xx * inpos[0] + xy * inpos[1] + xz * inpos[2];
outpos[1] = yx * inpos[0] + yy * inpos[1] + yz * inpos[2];
outpos[2] = zx * inpos[0] + zy * inpos[1] + zz * inpos[2];
}
}
static astro_rotation_t nutation_rot(astro_time_t *time, int direction)
{
@@ -1036,6 +998,14 @@ static astro_rotation_t nutation_rot(astro_time_t *time, int direction)
return rotation;
}
static void nutation(astro_time_t *time, int direction, const double inpos[3], double outpos[3])
{
astro_rotation_t r = nutation_rot(time, direction);
outpos[0] = r.rot[0][0]*inpos[0] + r.rot[1][0]*inpos[1] + r.rot[2][0]*inpos[2];
outpos[1] = r.rot[0][1]*inpos[0] + r.rot[1][1]*inpos[1] + r.rot[2][1]*inpos[2];
outpos[2] = r.rot[0][2]*inpos[0] + r.rot[1][2]*inpos[1] + r.rot[2][2]*inpos[2];
}
static double era(double ut) /* Earth Rotation Angle */
{
double thet1 = 0.7790572732640 + 0.00273781191135448 * ut;

46
source/c/astronomy.c
View File

@@ -1111,44 +1111,6 @@ static astro_equatorial_t vector2radec(const double pos[3])
return equ;
}
static void nutation(astro_time_t *time, int direction, const double inpos[3], double outpos[3])
{
earth_tilt_t tilt = e_tilt(time);
double oblm = tilt.mobl * DEG2RAD;
double oblt = tilt.tobl * DEG2RAD;
double psi = tilt.dpsi * ASEC2RAD;
double cobm = cos(oblm);
double sobm = sin(oblm);
double cobt = cos(oblt);
double sobt = sin(oblt);
double cpsi = cos(psi);
double spsi = sin(psi);
double xx = cpsi;
double yx = -spsi * cobm;
double zx = -spsi * sobm;
double xy = spsi * cobt;
double yy = cpsi * cobm * cobt + sobm * sobt;
double zy = cpsi * sobm * cobt - cobm * sobt;
double xz = spsi * sobt;
double yz = cpsi * cobm * sobt - sobm * cobt;
double zz = cpsi * sobm * sobt + cobm * cobt;
if (direction == 0)
{
/* forward rotation */
outpos[0] = xx * inpos[0] + yx * inpos[1] + zx * inpos[2];
outpos[1] = xy * inpos[0] + yy * inpos[1] + zy * inpos[2];
outpos[2] = xz * inpos[0] + yz * inpos[1] + zz * inpos[2];
}
else
{
/* inverse rotation */
outpos[0] = xx * inpos[0] + xy * inpos[1] + xz * inpos[2];
outpos[1] = yx * inpos[0] + yy * inpos[1] + yz * inpos[2];
outpos[2] = zx * inpos[0] + zy * inpos[1] + zz * inpos[2];
}
}
static astro_rotation_t nutation_rot(astro_time_t *time, int direction)
{
@@ -1205,6 +1167,14 @@ static astro_rotation_t nutation_rot(astro_time_t *time, int direction)
return rotation;
}
static void nutation(astro_time_t *time, int direction, const double inpos[3], double outpos[3])
{
astro_rotation_t r = nutation_rot(time, direction);
outpos[0] = r.rot[0][0]*inpos[0] + r.rot[1][0]*inpos[1] + r.rot[2][0]*inpos[2];
outpos[1] = r.rot[0][1]*inpos[0] + r.rot[1][1]*inpos[1] + r.rot[2][1]*inpos[2];
outpos[2] = r.rot[0][2]*inpos[0] + r.rot[1][2]*inpos[1] + r.rot[2][2]*inpos[2];
}
static double era(double ut) /* Earth Rotation Angle */
{
double thet1 = 0.7790572732640 + 0.00273781191135448 * ut;