The observations Newton used to calculate the orbit of Kirch’s comet can be found in the tables in http://www.newtonproject.ox.ac.uk/view/texts/diplomatic/NATP00301 page 491. Observations up to 5 February are by Flamsteed; later observations are by Newton.
Different versions of these tables exist. This version is sufficient to demonstrate the calculation.
I have not found the Newton Project’s licensing terms; their terms will also extend to this post.
I used the code below and manual editing to put the data in a consistent programmable format.
I converted Newton’s fractions to seconds in the latitude column, and omitted them elsewhere
I converted Newton’s times to solar longitude using https://clearskytonight.com/projects/astronomycalculator/sun/sunlongitude.html.
Update: I have posted this code at https://bitbucket.org/AlistairWall/blog/src/master/comets/. There are two files: preparation.hs is the code to regularise the data and Observations.hs is the library that will be used in the rest of the project. (the Main.* files are still a test.)
import Data.List.Split (splitOn)
import Data.Char (ord)
{- "♈♓" = "\9800\9811" -}
{- observations from http://www.newtonproject.ox.ac.uk/view/texts/diplomatic/NATP00301 page 491: the Newton Project's licencing terms apply. -}
-- dummy value ♈0.0.0 added for Newton's solar longitudes
observations'::[String]
observations' = [
"1680 December 12,4.46,4.46.00,♑1.53.2,♑6.33.0,8.26.0",
"1680 December 21,6.32,6.36.59,♑11.8.10,♒5.7.38,21.45.30",
"1680 December 24,6.12,6.17.52,♑14.10.49,♒18.49.10,25.23.24",
"1680 December 26,5.14,5.20.44,♑16.10.38,♒28.24.6,27.0.57",
"1680 December 29,7.55,8.03.2,♑19.20.56,♓13.11.45,28.10.05",
"1680 December 30,8.2,8.10.26,♑20.22.20,♓17.37.5,28.11.12",
"1681 January 5,5.51,6.1.38,♑26.23.19,♈8.49.10,26.15.26",
"1681 January 9,6.49,7.0.53,♒0.29.54,♈18.43.18,24.12.42",
"1681 January 10,5.54,6.6.10,♒1.28.34,♈20.40.57,23.44.0",
"1681 January 13,6.56,7.8.55,♒4.34.6,♈25.59.34,22.17.36",
"1681 January 25,7.44,7.58.42,♒16.45.58,♉9.55.48,17.56.54",
"1681 January 30,8.07,8.21.53,♒21.50.9,♉13.19.36,16.40.57",
"1681 February 2,6.20,6.34.51,♒24.47.41,♉15.13.48,16.02.02",
"1681 February 5,6.50,7.4.41,♒27.49.51,♉16.59.52,15.27.23",
"1681 February 25,8h.30,,♈0.0.0,♉26.19.22,12.46.43",
"1681 February 27,8.15,,♈0.0.0,♉27.4.28,12.36.0",
"1681 March 1,11.0,,♈0.0.0,♉27.53.8,12.24.26",
"1681 March 2,8.0,,♈0.0.0,♉28.12.29,12.19.30",
"1681 March 5,11.30,,♈0.0.0,♉29.20.51,12.2.40",
"1681 March 9,8.30,,♈0.0.0,♊0.43.2,11.44.36"]
observations'' = map (splitOn",") observations'
data Angle = Angle {deg::Int, min::Int, sec::Int}
instance Show Angle where
show = show.degrees
-- I am not writing a parser now because I want to get on to the geometry.
parts::String->[Int]
parts s = map read $ splitOn "." s
-- read angle written as deg.min.sec
readAngle::String->Angle
readAngle s = readAngle' $ parts s
readAngle'::[Int]->Angle
readAngle' [] = Angle 0 0 0 -- dummy value which I will replace with external values for Newton's solar longitudes.
readAngle' (d:m:s:[]) = Angle d m s
readAngle' _ = undefined -- all the observations are valid, so no need to handle errors
-- read angle written with constellation symbol
readAngleAstro::String->Angle
readAngleAstro (x:xs) =
let stars = ((ord x) - (ord '♈')) * 30
(Angle d m s) = readAngle xs
in Angle (d + stars) m s
readAngleAstro _ = undefined -- all the observations are valid, so no need to handle errors
degrees::Angle->Float
degrees (Angle d m s) = fromIntegral d + fromIntegral m /60.0 + fromIntegral s/3600.0
data Observation = Observation {
date::String,
time::String,
trueLongitude::String,
solarLongitude::Float,
cometLongitude::Float,
cometAscension::Float } deriving (Read, Show)
fromList::[String]->Observation
fromList (dt:tm:trueL:solarL:cometLong:cometLat:[]) =
Observation dt tm trueL (degrees $ readAngleAstro solarL) (degrees $ readAngleAstro cometLong) (degrees $ readAngle cometLat)
observations''' = map fromList observations''
main = mapM_ putStrLn $ map show observations'''
-- manually adding Newton's solar longitudes (https://clearskytonight.com/projects/astronomycalculator/sun/sunlongitude.html) to the output of main gives
observations::[Observation]
observations = [
Observation {date = "1680 December 12", time = "4.46", trueLongitude = "4.46.00", solarLongitude = 271.88388, cometLongitude = 276.55, cometAscension = 8.433333},
Observation {date = "1680 December 21", time = "6.32", trueLongitude = "6.36.59", solarLongitude = 281.1361, cometLongitude = 305.12723, cometAscension = 21.758333},
Observation {date = "1680 December 24", time = "6.12", trueLongitude = "6.17.52", solarLongitude = 284.18027, cometLongitude = 318.81946, cometAscension = 25.39},
Observation {date = "1680 December 26", time = "5.14", trueLongitude = "5.20.44", solarLongitude = 286.17722, cometLongitude = 328.40167, cometAscension = 27.015833},
Observation {date = "1680 December 29", time = "7.55", trueLongitude = "8.03.2", solarLongitude = 289.3489, cometLongitude = 343.19583, cometAscension = 28.168055},
Observation {date = "1680 December 30", time = "8.2", trueLongitude = "8.10.26", solarLongitude = 290.37222, cometLongitude = 347.61807, cometAscension = 28.186666},
Observation {date = "1681 January 5", time = "5.51", trueLongitude = "6.1.38", solarLongitude = 296.3886, cometLongitude = 8.819445, cometAscension = 26.257223},
Observation {date = "1681 January 9", time = "6.49", trueLongitude = "7.0.53", solarLongitude = 300.49835, cometLongitude = 18.721666, cometAscension = 24.211668},
Observation {date = "1681 January 10", time = "5.54", trueLongitude = "6.6.10", solarLongitude = 301.4761, cometLongitude = 20.682499, cometAscension = 23.733334},
Observation {date = "1681 January 13", time = "6.56", trueLongitude = "7.8.55", solarLongitude = 304.56836, cometLongitude = 25.992779, cometAscension = 22.293333},
Observation {date = "1681 January 25", time = "7.44", trueLongitude = "7.58.42", solarLongitude = 316.7661, cometLongitude = 39.93, cometAscension = 17.948332},
Observation {date = "1681 January 30", time = "8.07", trueLongitude = "8.21.53", solarLongitude = 321.83585, cometLongitude = 43.326664, cometAscension = 16.682499},
Observation {date = "1681 February 2", time = "6.20", trueLongitude = "6.34.51", solarLongitude = 324.7947, cometLongitude = 45.23, cometAscension = 16.033888},
Observation {date = "1681 February 5", time = "6.50", trueLongitude = "7.4.41", solarLongitude = 327.83084, cometLongitude = 46.99778, cometAscension = 15.456388},
Observation {date = "1681 February 25", time = "8h.30", trueLongitude = "", solarLongitude = 337.43662733763, cometLongitude = 56.322777, cometAscension = 12.778611},
Observation {date = "1681 February 27", time = "8.15", trueLongitude = "", solarLongitude = 339.43057972743, cometLongitude = 57.074444, cometAscension = 12.6},
Observation {date = "1681 March 1", time = "11.0", trueLongitude = "", solarLongitude = 341.54725185297, cometLongitude = 57.88556, cometAscension = 12.407222},
Observation {date = "1681 March 2", time = "8.0", trueLongitude = "", solarLongitude = 342.422357371082, cometLongitude = 58.208057, cometAscension = 12.325},
Observation {date = "1681 March 5", time = "11.30", trueLongitude = "", solarLongitude = 345.565036319878, cometLongitude = 59.3475, cometAscension = 12.044445},
Observation {date = "1681 March 9", time = "8.30", trueLongitude = "", solarLongitude = 349.429015768727, cometLongitude = 60.717224, cometAscension = 11.743334}]
Great Wall 