Coverage for tests / test_moid.py: 25%

1# This file is part of pipe_tasks. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21 

22import unittest 

23import lsst.utils.tests 

24from lsst.pipe.tasks.ssp.moid import MOIDSolver, earth_orbit 

25 

26# J2000.0 epoch as MJD 

27MJD_J2000 = 51544.5 

28 

29 

30class TestEarthOrbit(lsst.utils.tests.TestCase): 

31 

32 def testEarthOrbitJ2000(self): 

33 """Verify J2000 values match JPL coefficients.""" 

34 el = earth_orbit(MJD_J2000) 

35 self.assertAlmostEqual(el.a_AU, 1.00000261, places=8) 

36 self.assertAlmostEqual(el.e, 0.01671123, places=8) 

37 self.assertAlmostEqual(el.omega_deg, 102.93768193, places=5) 

38 self.assertEqual(el.Omega_deg, 0.0) 

39 self.assertAlmostEqual(el.inc_deg, 0.00005, places=6) 

40 

41 def testEarthOrbitSecularEvolution(self): 

42 """Verify elements evolve at the expected rate.""" 

43 el_j2000 = earth_orbit(MJD_J2000) 

44 # One century later 

45 el_later = earth_orbit(MJD_J2000 + 36525.0) 

46 

47 # omega should increase by ~0.323 deg/century 

48 domega = el_later.omega_deg - el_j2000.omega_deg 

49 self.assertAlmostEqual(domega, 0.32327364, places=5) 

50 

51 # eccentricity should decrease by ~0.0000439/century 

52 de = el_later.e - el_j2000.e 

53 self.assertAlmostEqual(de, -0.00004392, places=8) 

54 

55 # semi-major axis should increase slightly 

56 da = el_later.a_AU - el_j2000.a_AU 

57 self.assertAlmostEqual(da, 0.00000562, places=8) 

58 

59 def testEarthOrbitEclipticFrame(self): 

60 """Verify Omega=0 and inc~0 at all epochs.""" 

61 for mjd in [40000.0, MJD_J2000, 60000.0, 80000.0]: 

62 el = earth_orbit(mjd) 

63 self.assertEqual(el.Omega_deg, 0.0) 

64 self.assertAlmostEqual(el.inc_deg, 0.00005, places=6) 

65 

66 

67class TestMOIDSolver(lsst.utils.tests.TestCase): 

68 

69 def setUp(self): 

70 self.solver = MOIDSolver() 

71 

72 def testIdenticalOrbits(self): 

73 """Two identical circular orbits should have MOID = 0.""" 

74 el = (1.0, 0.0, 0.0, 0.0, 0.0) 

75 result = self.solver.compute(el, el) 

76 self.assertAlmostEqual(result.MOID_AU, 0.0, places=6) 

77 

78 def testCoplanarCircular(self): 

79 """Coplanar circular orbits a=1, a=2 should have MOID=1.0.""" 

80 el1 = (1.0, 0.0, 0.0, 0.0, 0.0) 

81 el2 = (2.0, 0.0, 0.0, 0.0, 0.0) 

82 result = self.solver.compute(el1, el2) 

83 self.assertAlmostEqual(result.MOID_AU, 1.0, places=4) 

84 

85 def testCoplanarCrossing(self): 

86 """Eccentric orbit crossing a circular one should have 

87 MOID near zero. 

88 """ 

89 # a=1 e=0 (circular) vs a=1.5 e=0.5 

90 # (perihelion=0.75, aphelion=2.25 — crosses the r=1 circle) 

91 el1 = (1.0, 0.0, 0.0, 0.0, 0.0) 

92 el2 = (1.5, 0.5, 0.0, 0.0, 0.0) 

93 result = self.solver.compute(el1, el2) 

94 self.assertLess(result.MOID_AU, 0.001) 

95 

96 def testPerpendicularOrbit(self): 

97 """90-degree inclination orbit should have a finite, 

98 positive MOID. 

99 """ 

100 earth = earth_orbit(60000.0) 

101 perp = (1.5, 0.3, 90.0, 0.0, 0.0) 

102 result = self.solver.compute(earth, perp) 

103 self.assertGreater(result.MOID_AU, 0.0) 

104 self.assertLess(result.MOID_AU, 2.0) 

105 

106 def testKnownNEA(self): 

107 """Apophis-like elements should give MOID ~ 0.0002 AU.""" 

108 # Approximate elements for (99942) Apophis 

109 apophis = (0.9224, 0.1914, 3.339, 204.43, 126.39) 

110 earth = earth_orbit(60000.0) 

111 result = self.solver.compute(earth, apophis) 

112 # MPC gives ~0.00025 AU; we allow some tolerance for 

113 # element epoch differences 

114 self.assertLess(result.MOID_AU, 0.005) 

115 self.assertGreater(result.MOID_AU, 0.0) 

116 

117 def testDeltaVPositive(self): 

118 """DeltaV should always be positive.""" 

119 earth = earth_orbit(60000.0) 

120 for el in [ 

121 (1.5, 0.3, 10.0, 30.0, 45.0), 

122 (2.5, 0.1, 5.0, 120.0, 200.0), 

123 (0.9224, 0.1914, 3.339, 204.43, 126.39), 

124 ]: 

125 result = self.solver.compute(earth, el) 

126 self.assertGreater(result.DeltaV_kms, 0.0) 

127 

128 def testResultFields(self): 

129 """MOIDResult should have all expected fields.""" 

130 el1 = (1.0, 0.1, 5.0, 30.0, 45.0) 

131 el2 = (1.5, 0.2, 15.0, 60.0, 10.0) 

132 result = self.solver.compute(el1, el2) 

133 

134 self.assertTrue(hasattr(result, 'MOID_AU')) 

135 self.assertTrue(hasattr(result, 'DeltaV_kms')) 

136 self.assertTrue(hasattr(result, 'EclipticLongitude_deg')) 

137 self.assertTrue(hasattr(result, 'TrueAnomaly1_deg')) 

138 self.assertTrue(hasattr(result, 'TrueAnomaly2_deg')) 

139 

140 # Angles should be in [0, 360) 

141 self.assertGreaterEqual(result.EclipticLongitude_deg, 0.0) 

142 self.assertLess(result.EclipticLongitude_deg, 360.0) 

143 self.assertGreaterEqual(result.TrueAnomaly1_deg, 0.0) 

144 self.assertLess(result.TrueAnomaly1_deg, 360.0) 

145 self.assertGreaterEqual(result.TrueAnomaly2_deg, 0.0) 

146 self.assertLess(result.TrueAnomaly2_deg, 360.0) 

147 

148 

149class MemoryTester(lsst.utils.tests.MemoryTestCase): 

150 pass 

151 

152 

153def setup_module(module): 

154 lsst.utils.tests.init() 

155 

156 

157if __name__ == "__main__": 157 ↛ 158line 157 didn't jump to line 158 because the condition on line 157 was never true

158 lsst.utils.tests.init() 

159 unittest.main()