module Data.Astro.CoordinateTest ( tests , testEC1 , testEcC ) where import Test.Framework (testGroup) import Test.Framework.Providers.HUnit import Test.Framework.Providers.QuickCheck2 (testProperty) import Test.HUnit import Test.HUnit.Approx import Test.QuickCheck import Control.Monad (unless) import Data.Astro.TypesTest (testDecimalDegrees, testDecimalHours) import Data.Astro.Time.JulianDate (JulianDate(..), LocalCivilTime(..) , fromYMDHMS, lctFromYMDHMS) import Data.Astro.Coordinate import Data.Astro.Types import Data.Astro.Utils (reduceToZeroRange) ro :: GeographicCoordinates ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5)) dt :: LocalCivilTime dt = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0 jd :: JulianDate jd = lctUniversalTime dt sunHC :: HorizonCoordinates sunHC = HC (fromDMS 49 18 21.77) (fromDMS 118 55 19.53) sunEC1 :: EquatorialCoordinates1 sunEC1 = EC1 {e1Declination = DD 23.378295912623855, e1RightAscension = DH 6.29383725890224} tests = [testGroup "RA <-> HA" [ testDecimalHours "RA 18.53 in 1980-04-22 18:36:51.67 to HA" 0.0000001 (DH 9.873237) (raToHA (DH 18.539167) (DD (-64.0)) $ fromYMDHMS 1980 4 22 18 36 51.67) , testDecimalHours "HA 9.87 in 1980-04-22 18:36:51.67 to RA" 0.0000001 (DH 9.873237) (raToHA (DH 18.539167) (DD (-64.0)) $ fromYMDHMS 1980 4 22 18 36 51.67) , testProperty "RA <-> HA property for Novosibirsk" $ prop_HARAConv (DD 83) (JD 2457607.97281) , testProperty "RA <-> HA property for Rio de Janeiro" $ prop_HARAConv (DD (-43)) (JD 2457607.97281) ] , testGroup "EC2 <-> HC" [ testHC "EC2 23.219 5.862 -> HC 19.334 283.271" 0.00000001 (HC (DD 19.334345224) (DD 283.27102726)) (equatorialToHorizon (DD 52) (EC2 (DD 23.219444444) (DH 5.862222222222222))) , testEC2 "HC 19.334 283.271 -> EC2 23.219 5.862" 0.00000001 (EC2 (DD 23.219444444) (DH 5.862222222222222)) (horizonToEquatorial (DD 52) (HC (DD 19.334345224) (DD 283.27102726))) , testProperty "property" prop_EC2HCConv ] , testGroup "EC1 <-> HC" [ testHC "sunEC1 -> sunHC" 0.00000001 sunHC (ec1ToHC ro jd sunEC1) , testEC1 "sunHC -> sunEC1" 0.00000001 sunEC1 (hcToEC1 ro jd sunHC) ] , testGroup "obliquity" [ testDecimalDegrees "2009-07-06" 0.000001 (DD 23.43925285) (obliquity $ JD 2455018.5) , testDecimalDegrees "2016-08-12" 0.000001 (DD 23.434706359) (obliquity $ JD 2457612.5) ] , testGroup "Ecliptic <-> Equatorial" [ testEC1 "EcC 4.875 139.6861 on 2009-07-06 to EC1" 0.000001 (EC1 (DD 19.535712) (DH 9.581501)) (eclipticToEquatorial (EcC (DD 4.875277777777778) (DD 139.6861111111111)) (JD 2455018.5)) , testEcC "EC1 19.535 9.581 to EcC on 2009-07-06" 0.000001 (EcC (DD 4.875278) (DD 139.686111)) (equatorialToEcliptic (EC1 (DD 19.535712) (DH 9.581501)) (JD 2455018.5)) , testProperty "property" prop_EC1EcCConv ] , testGroup "Galactic <-> Equatorial" [ testGC "EC 10.053 10.35 -> GC" 0.000001 (GC (DD 51.12226779935766) (DD 232.24788348766026)) (equatorialToGalactic (EC1 (DD 10.053055555555556) (DH 10.35))) , testEC1 "GC 51.122 232.247 -> EC" 0.000001 (EC1 (DD 10.053055555555556) (DH 10.35)) (galacticToEquatorial (GC (DD 51.12226779935766) (DD 232.24788348766026))) , testProperty "property" prop_EC1GCConv ] ] prop_HARAConv longitude ut dh = let dh' = reduceToZeroRange 24 dh DH ra = haToRA (raToHA (DH dh') longitude ut) longitude ut DH ha = raToHA (haToRA (DH dh') longitude ut) longitude ut eps = 0.00000001 in abs (ra - dh') < eps && (ha - dh') < eps where types = (dh::Double) testHC msg eps expected actual = testCase msg $ assertHC eps expected actual assertHC eps expected@(HC (DD eAlt) (DD eAz)) actual@(HC (DD aAlt) (DD aAz)) = unless (abs(eAlt-aAlt) <= eps && abs(eAz-aAz) <= eps) (assertFailure msg) where msg = "expected: " ++ show expected ++ "\n but got: " ++ show actual ++ "\n (maximum margin of error: " ++ show eps ++ ")" testEC2 msg eps expected actual = testCase msg $ assertEC2 eps expected actual assertEC2 eps expected@(EC2 (DD eDec) (DH eHa)) actual@(EC2 (DD aDec) (DH aHa)) = unless (abs(eDec-aDec) <= eps && abs(eHa-aHa) <= eps) (assertFailure msg) where msg = "expected: " ++ show expected ++ "\n but got: " ++ show actual ++ "\n (maximum margin of error: " ++ show eps ++ ")" prop_EC2HCConv (latitude, up, round) = let latitude' = DD $ reduceToZeroRange 180 latitude up' = (reduceToZeroRange 89.9 up) + 0.1 round' = (reduceToZeroRange 359.9 round) + 0.1 HC (DD up'') (DD round'') = equatorialToHorizon latitude' $ horizonToEquatorial latitude' (HC (DD up') (DD round')) eps = 0.00001 in abs(up'-up'') < eps && abs(round'-round'') < eps where types = ((latitude, up, round)::(Double, Double, Double)) testEC1 msg eps expected actual = testCase msg $ assertEC1 eps expected actual assertEC1 eps expected@(EC1 (DD e1) (DH e2)) actual@(EC1 (DD a1) (DH a2)) = unless (abs(e1-a1) <= eps && abs(e2-a2) <= eps) (assertFailure msg) where msg = "expected: " ++ show expected ++ "\n but got: " ++ show actual ++ "\n (maximum margin of error: " ++ show eps ++ ")" testEcC msg eps expected actual = testCase msg $ assertEcC eps expected actual assertEcC eps expected@(EcC (DD eLat) (DD eLon)) actual@(EcC (DD aLat) (DD aLon)) = unless (abs(eLat-aLat) <= eps && abs(eLon-aLon) <= eps) (assertFailure msg) where msg = "expected: " ++ show expected ++ "\n but got: " ++ show actual ++ "\n (maximum margin of error: " ++ show eps ++ ")" prop_EC1EcCConv (jd, up, round) = let jd' = JD $ (reduceToZeroRange 10000 jd) + 2455018.5 up' = (reduceToZeroRange 90 up) round' = (reduceToZeroRange 359 round) EcC (DD up'') (DD round'') = equatorialToEcliptic (eclipticToEquatorial (EcC (DD up') (DD round')) jd') jd' eps = 0.000001 in abs(up'-up'') < eps && abs(round'-round'') < eps where types = ((jd, up, round)::(Double, Double, Double)) testGC msg eps expected actual = testCase msg $ assertGC eps expected actual assertGC eps expected@(GC (DD eLat) (DD eLon)) actual@(GC (DD aLat) (DD aLon)) = unless (abs(eLat-aLat) <= eps && abs(eLon-aLon) <= eps) (assertFailure msg) where msg = "expected: " ++ show expected ++ "\n but got: " ++ show actual ++ "\n (maximum margin of error: " ++ show eps ++ ")" prop_EC1GCConv (up, round) = let up' = (reduceToZeroRange 90 up) round' = (reduceToZeroRange 359 round) GC (DD up'') (DD round'') = equatorialToGalactic $ galacticToEquatorial (GC (DD up') (DD round')) eps = 0.000001 in abs(up'-up'') < eps && abs(round'-round'') < eps where types = ((up, round)::(Double, Double))