行业智能>> astro>> 返回
项目作者： aligusnet

项目描述 ：
Amateur astronomical computations
高级语言： Haskell
项目地址: git://github.com/aligusnet/astro.git
创建时间： 2016-08-26T05:31:26Z
项目社区：https://github.com/aligusnet/astro
开源协议：Other
下载


Amateur astronomical computations

Build Status
Coverage Status
Documentation
Hackage

Usage

Install stack

It will take care of everything.

Get started with Stack

Build the project

  1. stack build

Run unit tests

  1. stack test

Documentation

Useful types

Decimal hours and Decimal degrees

Types to represent hours (used in celestial coordinate systems and as time zone) and degrees (used in coordinate systems).

  1. import Data.Astro.Types
  3. -- 10h 15m 19.7s
  4. dh :: DecimalHours
  5. dh = fromHMS 10 15 19.7
  6. -- DH 10.255472222222222
  8. (h, m, s) = toHMS dh
  9. -- (10,15,19.699999999999562)
  12. -- 51°28′40″
  13. dd :: DecimalDegrees
  14. dd = fromDMS 51 28 40
  15. -- DD 51.477777777777774
  17. (d, m, s) = toDMS dd
  18. -- (51,28,39.999999999987494)

Geographic Coordinates

  1. import Data.Astro.Types
  3. -- the Royal Observatory, Greenwich
  4. ro :: GeographicCoordinates
  5. ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5))
  6. -- GeoC {geoLatitude = DD 51.4778, geoLongitude = DD (-0.0014)}

Time

The main time datetime type used in the library is JulianDate defined in Data.Astro.Time.JulianDate. JulianDate is just a number of days since noon of 1 January 4713 BC:

  1. import Data.Astro.Time.JulianDate
  3. -- 2017-06-25 9:29:00 (GMT)
  4. jd :: JulianDate
  5. jd = fromYMDHMS 2017 6 25 9 29 0
  6. -- JD 2457929.895138889

LocalCiviTime and LocalCivilDate are Julian dates with time zones:

  1. import Data.Astro.Time.JulianDate
  2. import Data.Astro.Types
  4. -- 2017-06-25 10:29:00 +0100 (BST)
  5. lct :: LocalCivilTime
  6. lct = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0
  7. -- 2017-06-25 10:29:00.0000 +1.0
  9. lctJD :: JulianDate
  10. lctJD = lctUniversalTime lct
  11. -- JD 2457929.895138889
  13. lctTZ :: DecimalHours
  14. lctTZ = lctTimeZone lct
  15. -- DH 1.0
  17. lcd :: LocalCivilDate
  18. lcd = lcdFromYMD (DH 1) 2017 6 25
  20. lcdJD :: JulianDate
  21. lcdJD = lcdDate lcd
  22. -- JD 2457929.5
  24. lcdTZ :: DecimalHours
  25. lcdTZ = lcdTimeZone lcd
  26. -- DH 1.0

Celestial coordinate systems

The celestical coordinate systems are defined in Data.Astro.Coordinate.

If you would like to locate Sirius in the sky you need to know the altitude or ‘how far up’ angle in the sky and azimuth - ‘how far round’ angle from the north direction to the east. this describes the Horizontal coordinate system:

alt Horizontal coordinate system

  1. import Data.Astro.Coordinate
  2. import Data.Astro.Types
  4. hc :: HorizonCoordinates
  5. hc = HC (DD 30.5) (DD 180)
  6. -- HC {hAltitude = DD 30.0, hAzimuth = DD 180.0}

Unfortunately the Horizontal coordinate system values depend on the position of the observer. And it’s not handy when you need to share coordinates of some celestial object with your friend in Japan.

The second coordinate system is the Equatorial coordinate system. This coordinate system uses the location of the centre of the Earth as the zero point so it does not depend on the observer’s location.

We have two flavours of equatorial coordinates:

  • the first one uses the vernal equinox as a starting direction for the ‘how far round’ coordinate (right ascension, α),

  • the second one uses the meridian instead of the vernal equinox (hour angle).

We can consider the second one as a transition coordinate system between the horizontal one and the ‘true’ equatorial one.

  1. import Data.Astro.Coordinate
  2. import Data.Astro.Types
  4. ec1 :: EquatorialCoordinates1
  5. ec1 = EC1 (DD 71.7) (DH 8)
  6. -- EC1 {e1Declination = DD 71.7, e1RightAscension = DH 8.0}
  8. ec2 :: EquatorialCoordinates2
  9. ec2 = EC1 (DD 77.7) (DH 11)
  10. -- EC2 {e2Declination = DD 77.7, e2HoursAngle = DH 11.0}

Transformations

Say, now is 2017-06-25 10:29 BST and we are somewhere near the Royal Observatory, Greenwich.

Let convert the current location of the Sun in horizon coordinates (altitude: 49°18′21.77″, azimuth: 118°55′19.53″) to equatorial coordinates and back to horizon ones:

  1. import Data.Astro.Time.JulianDate
  2. import Data.Astro.Coordinate
  3. import Data.Astro.Types
  5. ro :: GeographicCoordinates
  6. ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5))
  8. dt :: LocalCivilTime
  9. dt = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0
  11. sunHC :: HorizonCoordinates
  12. sunHC = HC (fromDMS 49 18 21.77) (fromDMS 118 55 19.53)
  13. -- HC {hAltitude = DD 49.30604722222222, hAzimuth = DD 118.92209166666666}
  15. sunEC2 :: EquatorialCoordinates2
  16. sunEC2 = horizonToEquatorial (geoLatitude ro) sunHC
  17. -- EC2 {e2Declination = DD 23.378295912623855, e2HoursAngle = DH 21.437117068873537}
  19. sunEC1 :: EquatorialCoordinates1
  20. sunEC1 = EC1 (e2Declination sunEC2) (haToRA (e2HoursAngle sunEC2) (geoLongitude ro) (lctUniversalTime dt))
  21. -- EC1 {e1Declination = DD 23.378295912623855, e1RightAscension = DH 6.29383725890224}
  24. sunEC2' :: EquatorialCoordinates2
  25. sunEC2' = EC2 (e1Declination sunEC1) (raToHA (e1RightAscension sunEC1) (geoLongitude ro) (lctUniversalTime dt))
  26. -- EC2 {e2Declination = DD 23.378295912623855, e2HoursAngle = DH 21.437117068873537}
  28. sunHC' :: HorizonCoordinates
  29. sunHC' = equatorialToHorizon (geoLatitude ro) sunEC2'
  30. -- HC {hAltitude = DD 49.30604722222222, hAzimuth = DD 118.92209166666666}

You can use function-shortcuts to simplify transformation EquatorialCoordinates1 <-> HorizonCoordinates: ec1ToHC and hcToEC1:

  1. import Data.Astro.Time.JulianDate
  2. import Data.Astro.Coordinate
  3. import Data.Astro.Types
  5. ro :: GeographicCoordinates
  6. ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5))
  8. dt :: LocalCivilTime
  9. dt = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0
  11. sunHC :: HorizonCoordinates
  12. sunHC = HC (fromDMS 49 18 21.77) (fromDMS 118 55 19.53)
  13. -- HC {hAltitude = DD 49.30604722222222, hAzimuth = DD 118.92209166666666}
  15. sunEC1 :: EquatorialCoordinates1
  16. sunEC1 = hcToEC1 ro (lctUniversalTime dt) sunHC
  17. -- EC1 {e1Declination = DD 23.378295912623855, e1RightAscension = DH 6.29383725890224}
  19. sunHC' :: HorizonCoordinates
  20. sunHC' = ec1ToHC ro (lctUniversalTime dt) sunEC1
  21. -- HC {hAltitude = DD 49.30604722222222, hAzimuth = DD 118.92209166666666}

Stars

The ancient astronomers noted that there were 2 types of stars: some of them were fixed, travelling the same way across the sky every sidereal day and another were wanderers (planetai in ancient Greek).

Of course, stars are not fixed, they are travelling with high speeds but distances to them are so high that their movement is very difficult to note. So we can assume that they are fixed for our purposes.

Given the “fixed” equatorial coordinates of the star we only need to transform them to the horizon coordinates to find out where the star in the sky.

In the example below we will use Data.Astro.Star module which defines equatorial coordinates of some stars:

  1. import Data.Astro.Time.JulianDate
  2. import Data.Astro.Coordinate
  3. import Data.Astro.Types
  4. import Data.Astro.Star
  7. ro :: GeographicCoordinates
  8. ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5))
  10. dt :: LocalCivilTime
  11. dt = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0
  13. -- Calculate location of Betelgeuse
  15. betelgeuseEC1 :: EquatorialCoordinates1
  16. betelgeuseEC1 = starCoordinates Betelgeuse
  17. -- EC1 {e1Declination = DD 7.407064, e1RightAscension = DH 5.919529}
  19. betelgeuseHC :: HorizonCoordinates
  20. betelgeuseHC = ec1ToHC ro (lctUniversalTime dt) betelgeuseEC1
  21. -- HC {hAltitude = DD 38.30483892505852, hAzimuth = DD 136.75755644642248}

Rise and Set

Data.Astro.CelestialObject.RiseSet module defines RiseSet type to represent time and azimuth of rise and set.

Let calculate rise and set time of Rigel:

  1. import Data.Astro.Time.JulianDate
  2. import Data.Astro.Coordinate
  3. import Data.Astro.Types
  4. import Data.Astro.Effects
  5. import Data.Astro.CelestialObject.RiseSet
  6. import Data.Astro.Star
  9. ro :: GeographicCoordinates
  10. ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5))
  12. today :: LocalCivilDate
  13. today = lcdFromYMD (DH 1) 2017 6 25
  15. -- Calculate location of Betelgeuse
  17. rigelEC1 :: EquatorialCoordinates1
  18. rigelEC1 = starCoordinates Rigel
  20. verticalShift :: DecimalDegrees
  21. verticalShift = refract (DD 0) 12 1012
  22. -- DD 0.5660098245614035
  24. rigelRiseSet :: RiseSetLCT
  25. rigelRiseSet = riseAndSetLCT ro today verticalShift rigelEC1
  26. -- RiseSet (2017-06-25 06:38:18.4713 +1.0,DD 102.51249855335433) (2017-06-25 17:20:33.4902 +1.0,DD 257.48750144664564)

As we can see Rigel rose today at 06:38:18 and will set at 17:20:33, azimuths of rise and set 102.51° and 257.49° correspondingly.

We used refract function of Data.Astro.Effects module with reasonable default parameters to calculate vertical shift.

Planets

The planets is completely different story. We cannot assume that the planets have “fixed” location in equatorial coordinates like stars.

What we can do is to describe details of the planets’ orbit and calculate their positions at any given moment.

Planets and planet details are defined in Data.Astro.Planet module. j2010PlanetDetails returns details for the given planet.
This module also defines planetPosition, planetDistance1 and planetAngularDiameter to calculate position of the given planet, distance to the planet and angular size of the planet correspondingly.

1 at the end of the planetDistance1 means that this function uses not very precise method to do calculations. Sometimes there are 2-methods available in the library, but not always.

Let us do some planets-related calculations.

Do some initialisation:

  1. import Data.Astro.Time.JulianDate
  2. import Data.Astro.Coordinate
  3. import Data.Astro.Types
  4. import Data.Astro.Effects
  5. import Data.Astro.CelestialObject.RiseSet
  6. import Data.Astro.Planet
  8. ro :: GeographicCoordinates
  9. ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5))
  11. dt :: LocalCivilTime
  12. dt = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0
  14. today :: LocalCivilDate
  15. today = lcdFromYMD (DH 1) 2017 6 25
  17. jupiterDetails :: PlanetDetails
  18. jupiterDetails = j2010PlanetDetails Jupiter
  20. earthDetails :: PlanetDetails
  21. earthDetails = j2010PlanetDetails Earth
  23. jupiterPosition :: JulianDate -> EquatorialCoordinates1
  24. jupiterPosition = planetPosition planetTrueAnomaly1 jupiterDetails earthDetails

Calculate Jupiter’s coordinates:

  1. jupiterEC1 :: EquatorialCoordinates1
  2. jupiterEC1 = jupiterPosition (lctUniversalTime dt)
  3. -- EC1 {e1Declination = DD (-4.104626810672402), e1RightAscension = DH 12.863365504382228}
  5. jupiterHC :: HorizonCoordinates
  6. jupiterHC = ec1ToHC ro (lctUniversalTime dt) jupiterEC1
  7. -- HC {hAltitude = DD (-30.67914598469227), hAzimuth = DD 52.29376845044007}

As be can see Jupiter is below the horizon now (the altitude is negative), that’s unfortunate.

Now let us calculate distance to Jupiter:

  1. jupiterDistance :: AstronomicalUnits
  2. jupiterDistance = planetDistance1 jupiterDetails earthDetails (lctUniversalTime dt)
  3. -- AU 5.193435872521039

1 Astronomical Unit is an average distance from the Earth to the Sun.

and calculate an angular size now:

  1. jupiterAngularSize :: DecimalDegrees
  2. jupiterAngularSize = planetAngularDiameter jupiterDetails jupiterDistance
  3. -- DD 1.052289877865987e-2
  5. toDMS jupiterAngularSize
  6. -- (0,0,37.88243560317554)

Rise and Set

Calculate rise and set times of planets are not easy task, because planets change their equatorial coordinates during the day.

riseAndSet2 function of Data.Astro.CelestialObject.RiseSet module applies iterative approach: calculates rise and set date for midday coordinates and then recalculates rise time for rise coordinates and set for set coordinates obtained from the previous step:

  1. verticalShift :: DecimalDegrees
  2. verticalShift = refract (DD 0) 12 1012
  3. -- DD 0.5660098245614035
  5. jupiterRiseSet :: RiseSetMB
  6. jupiterRiseSet = riseAndSet2 0.000001 jupiterPosition ro verticalShift today
  7. -- RiseSet
  8. --    (Just (2017-06-25 13:53:27.3109 +1.0,DD 95.88943953535569))
  9. --    (Just (2017-06-25 01:21:23.5835 +1.0,DD 264.1289033612776))

We can see now why at 10 am Jupiter is below horizon because it will rise only at 1:53 pm.

Sun

Some examples of doing the Sun’s related calculations:

  1. import Data.Astro.Time.JulianDate
  2. import Data.Astro.Coordinate
  3. import Data.Astro.Types
  4. import Data.Astro.Sun
  6. ro :: GeographicCoordinates
  7. ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5))
  9. dt :: LocalCivilTime
  10. dt = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0
  12. today :: LocalCivilDate
  13. today = lcdFromYMD (DH 1) 2017 6 25
  15. jd :: JulianDate
  16. jd = lctUniversalTime dt
  18. verticalShift :: DecimalDegrees
  19. verticalShift = refract (DD 0) 12 1012
  21. -- distance from the Earth to the Sun in kilometres
  22. distance :: Double
  23. distance = sunDistance jd
  24. -- 1.5206375976421073e8
  26. -- Angular Size
  27. angularSize :: DecimalDegrees
  28. angularSize = sunAngularSize jd
  29. -- DD 0.5244849215333616
  31. -- The Sun's coordinates
  32. ec1 :: EquatorialCoordinates1
  33. ec1 = sunPosition2 jd
  34. -- EC1 {e1Declination = DD 23.37339098989099, e1RightAscension = DH 6.29262026252748}
  36. hc :: HorizonCoordinates
  37. hc = ec1ToHC ro jd ec1
  38. -- HC {hAltitude = DD 49.312050979507404, hAzimuth = DD 118.94723825710143}
  41. -- Rise and Set
  42. riseSet :: RiseSetMB
  43. riseSet = sunRiseAndSet ro 0.833333 today
  44. -- RiseSet
  45. --    (Just (2017-06-25 04:44:04.3304 +1.0,DD 49.043237261724215))
  46. --    (Just (2017-06-25 21:21:14.4565 +1.0,DD 310.91655607595595))

Moon

The Moon’s related calculations. Data.Astro.Moon module defines 2 new types of functions we haven’t seen before: moonPhase and moonBrightLimbPositionAngle which calculate the phase (the area of the visible segment expressed as a fraction of the whole disk) and the position-angle which is the angle of the midpoint of the illuminated limb measured eastwards from the north point of the disk.

  1. import Data.Astro.Time.JulianDate
  2. import Data.Astro.Coordinate
  3. import Data.Astro.Types
  4. import Data.Astro.Effects
  5. import Data.Astro.CelestialObject.RiseSet
  6. import Data.Astro.Moon
  8. ro :: GeographicCoordinates
  9. ro = GeoC (fromDMS 51 28 40) (-(fromDMS 0 0 5))
  11. dt :: LocalCivilTime
  12. dt = lctFromYMDHMS (DH 1) 2017 6 25 10 29 0
  14. today :: LocalCivilDate
  15. today = lcdFromYMD (DH 1) 2017 6 25
  17. jd :: JulianDate
  18. jd = lctUniversalTime dt
  20. -- distance from the Earth to the Moon in kilometres
  21. mdu :: MoonDistanceUnits
  22. mdu = moonDistance1 j2010MoonDetails jd
  23. -- MDU 0.9550170577020396
  25. distance :: Double
  26. distance = mduToKm mdu
  27. -- 367109.51199772174
  29. -- Angular Size
  30. angularSize :: DecimalDegrees
  31. angularSize = moonAngularSize mdu
  32. -- DD 0.5425033990980761
  34. -- The Moon's coordinates
  35. position :: JulianDate -> EquatorialCoordinates1
  36. position = moonPosition1 j2010MoonDetails
  38. ec1 :: EquatorialCoordinates1
  39. ec1 = position jd
  40. -- EC1 {e1Declination = DD 18.706180658927323, e1RightAscension = DH 7.56710547682055}
  42. hc :: HorizonCoordinates
  43. hc = ec1ToHC ro jd ec1
  44. -- HC {hAltitude = DD 34.57694951316064, hAzimuth = DD 103.91119101451832}
  46. -- Rise and Set
  47. riseSet :: RiseSetMB
  48. riseSet = riseAndSet2 0.000001 position ro verticalShift today
  49. -- RiseSet
  50. --    (Just (2017-06-25 06:22:51.4858 +1.0,DD 57.81458864497365))
  51. --    (Just (2017-06-25 22:28:20.3023 +1.0,DD 300.4168238905249))
  53. -- Phase
  54. phase :: Double
  55. phase = moonPhase j2010MoonDetails jd
  56. -- 2.4716141948212922e-2
  59. sunEC1 :: EquatorialCoordinates1
  60. sunEC1 = sunPosition2 jd
  61. -- EC1 {e1Declination = DD 23.37339098989099, e1RightAscension = DH 6.29262026252748}
  63. limbAngle :: DecimalDegrees
  64. limbAngle = moonBrightLimbPositionAngle ec1 sunEC1
  65. -- DD 287.9869373767473