diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 8c4ccd47..66aaee0c 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -3046,3 +3046,26 @@ def RotateVector(rotation, vector): rotation.rot[0][2]*vector.x + rotation.rot[1][2]*vector.y + rotation.rot[2][2]*vector.z, vector.t ) + + +def Rotation_EQJ_ECL(): + """Calculates a rotation matrix from equatorial J2000 (EQJ) to ecliptic J2000 (ECL). + + This is one of the family of functions that returns a rotation matrix + for converting from one orientation to another. + Source: EQJ = equatorial system, using equator at J2000 epoch. + Target: ECL = ecliptic system, using equator at J2000 epoch. + + Returns + ------- + RotationMatrix + A rotation matrix that converts EQJ to ECL. + """ + # ob = mean obliquity of the J2000 ecliptic = 0.40909260059599012 radians. + c = 0.9174821430670688 # cos(ob) + s = 0.3977769691083922 # sin(ob) + return RotationMatrix([ + [ 1, 0, 0], + [ 0, +c, -s], + [ 0, +s, +c] + ]) diff --git a/generate/test.py b/generate/test.py index 1f080f6e..c485d88d 100755 --- a/generate/test.py +++ b/generate/test.py @@ -45,7 +45,7 @@ def Test_GeoMoon(): cx, cy, cz = 0.002674036155459549, -0.0001531716308218381, -0.0003150201604895409 dx, dy, dz = vec.x - cx, vec.y - cy, vec.z - cz diff = math.sqrt(dx*dx + dy*dy + dz*dz) - print('Test_GeoMoon: diff = {}'.format(diff)) + print('Test_GeoMoon: diff = {}'.format(diff)) if diff > 4.34e-19: print('Test_GeoMoon: EXCESSIVE ERROR') sys.exit(1) @@ -80,8 +80,8 @@ def Test_AstroCheck(printflag): print('o {:0.6f} {:0.6f} {:0.6f}'.format(observer.latitude, observer.longitude, observer.height)) dt = 10 + math.pi/100 bodylist = [ - astronomy.Body.Sun, astronomy.Body.Moon, astronomy.Body.Mercury, astronomy.Body.Venus, - astronomy.Body.Earth, astronomy.Body.Mars, astronomy.Body.Jupiter, astronomy.Body.Saturn, + astronomy.Body.Sun, astronomy.Body.Moon, astronomy.Body.Mercury, astronomy.Body.Venus, + astronomy.Body.Earth, astronomy.Body.Mars, astronomy.Body.Jupiter, astronomy.Body.Saturn, astronomy.Body.Uranus, astronomy.Body.Neptune, astronomy.Body.Pluto ] @@ -89,7 +89,7 @@ def Test_AstroCheck(printflag): for body in bodylist: name = body.name if body != astronomy.Body.Moon: - pos = astronomy.HelioVector(body, time) + pos = astronomy.HelioVector(body, time) if printflag: print('v {} {:0.16f} {:0.16f} {:0.16f} {:0.16f}'.format(name, pos.t.tt, pos.x, pos.y, pos.z)) if body != astronomy.Body.Earth: @@ -165,7 +165,7 @@ def Test_Seasons(filename): diff_minutes = (24.0 * 60.0) * abs(calc_time.tt - correct_time.tt) if diff_minutes > max_minutes: max_minutes = diff_minutes - + if diff_minutes > 1.7: print('Test_Seasons: {} line {}: excessive error ({}): {} minutes.'.format(filename, lnum, name, diff_minutes)) return 1 @@ -538,7 +538,7 @@ def CheckSaturn(): def TestMaxMag(body, filename): # Example of input data: # - # 2001-02-21T08:00Z 2001-02-27T08:00Z 23.17 19.53 -4.84 + # 2001-02-21T08:00Z 2001-02-27T08:00Z 23.17 19.53 -4.84 # # JPL Horizons test data has limited floating point precision in the magnitude values. # There is a pair of dates for the beginning and end of the max magnitude period, @@ -716,7 +716,7 @@ def LunarApsis(filename): max_km = max(max_km, diff_km) print('LunarApsis: found {} events, max time error = {:0.3f} minutes, max distance error = {:0.3f} km.'.format(lnum, max_minutes, max_km)) return 0 - + def Test_Apsis(): if 0 != LunarApsis('apsides/moon.txt'): @@ -725,6 +725,18 @@ def Test_Apsis(): #----------------------------------------------------------------------------------------------------------- +def Test_Rotation(): + print('Python Test_Rotation: PASS') + return 0 + +#----------------------------------------------------------------------------------------------------------- + +def Test_Refraction(): + print('Python Test_Refraction: PASS') + return 0 + +#----------------------------------------------------------------------------------------------------------- + if __name__ == '__main__': if len(sys.argv) == 2: if sys.argv[1] == 'time': @@ -744,6 +756,10 @@ if __name__ == '__main__': sys.exit(Test_Apsis()) if sys.argv[1] == 'issue46': sys.exit(Test_Issue46()) + if sys.argv[1] == 'rotation': + sys.exit(Test_Rotation()) + if sys.argv[1] == 'refraction': + sys.exit(Test_Refraction()) if len(sys.argv) == 3: if sys.argv[1] == 'seasons': diff --git a/generate/unit_test_python b/generate/unit_test_python index 8805a860..4533b957 100755 --- a/generate/unit_test_python +++ b/generate/unit_test_python @@ -6,6 +6,8 @@ Fail() } python3 --version || Fail "Cannot print python version" +python3 test.py rotation || Fail "Failed Python rotation tests." +python3 test.py refraction || Fail "Failed Python refraction tests." python3 test.py time || Fail "Failure reported by test.py (time)" python3 test.py apsis || Fail "Failed Python apsis tests." python3 test.py magnitude || Fail "Failed Python magnitude tests." diff --git a/source/python/README.md b/source/python/README.md index 61fac607..808eaeca 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -911,6 +911,21 @@ A vector in the orientation specified by `rotation`. --- + +### Rotation_EQJ_ECL() + +**Calculates a rotation matrix from equatorial J2000 (EQJ) to ecliptic J2000 (ECL).** + +This is one of the family of functions that returns a rotation matrix +for converting from one orientation to another. +Source: EQJ = equatorial system, using equator at J2000 epoch. +Target: ECL = ecliptic system, using equator at J2000 epoch. + +### Returns: RotationMatrix +A rotation matrix that converts EQJ to ECL. + +--- + ### Search(func, context, t1, t2, dt_tolerance_seconds) diff --git a/source/python/astronomy.py b/source/python/astronomy.py index c0c26b22..b9ec400a 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -5107,3 +5107,26 @@ def RotateVector(rotation, vector): rotation.rot[0][2]*vector.x + rotation.rot[1][2]*vector.y + rotation.rot[2][2]*vector.z, vector.t ) + + +def Rotation_EQJ_ECL(): + """Calculates a rotation matrix from equatorial J2000 (EQJ) to ecliptic J2000 (ECL). + + This is one of the family of functions that returns a rotation matrix + for converting from one orientation to another. + Source: EQJ = equatorial system, using equator at J2000 epoch. + Target: ECL = ecliptic system, using equator at J2000 epoch. + + Returns + ------- + RotationMatrix + A rotation matrix that converts EQJ to ECL. + """ + # ob = mean obliquity of the J2000 ecliptic = 0.40909260059599012 radians. + c = 0.9174821430670688; # cos(ob) + s = 0.3977769691083922; # sin(ob) + return RotationMatrix([ + [ 1, 0, 0], + [ 0, +c, -s], + [ 0, +s, +c] + ])