From e7d8804bf1fa26e2fd6334eeb45202eec324b927 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Sat, 14 Dec 2019 18:15:13 -0500 Subject: [PATCH] C: Eliminated duplicate code - nutation() is based on nutation_rot(). 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. --- generate/template/astronomy.c | 46 ++++++----------------------------- source/c/astronomy.c | 46 ++++++----------------------------- 2 files changed, 16 insertions(+), 76 deletions(-) diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index 1c086bf0..d0e29b96 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -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; diff --git a/source/c/astronomy.c b/source/c/astronomy.c index 57c20637..d968106c 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -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;