Download or view mithengeSirius.frink in plain text format
// Program to calculate a star crossing of the Infinite Corridor at MIT
// known as "MIThenge"
//
// More info at also http://web.mit.edu/planning/www/mithenge.html
// Thanks to Keith Winstein and Ken Olum for various data.
//
// For worked predictions, see https://futureboy.us/mithenge/
//
// Alan Eliasen, eliasen@mindspring.com
use mithengecorridor.frink
use cambridgetempFourier.frink
use sun.frink
/* I attempted to find the most-accurate positions for Sirius from the GAIA
Data Release 2, but found a few things:
Sirius A evidently doesn't have an entry in the GAIA DR2 release, maybe
because it's so bright. Sirius B does have an entry, source_id
2947050466531873024
GAIA DR2 search engine:
https://gea.esac.esa.int/archive/
Sirius A is HIP32349
http://simbad.u-strasbg.fr/simbad/sim-basic?Ident=Sirius+A&submit=SIMBAD+search
Sirius B is Gaia DR2 source_id 2947050466531873024
http://simbad.u-strasbg.fr/simbad/sim-basic?Ident=Sirius+B&submit=SIMBAD+search
The GAIA DR2 data release currently uses epoch 2015.5 (column ref_epoch)
for all stars, but this may change in future releases. Thus,
right ascension and declination coordinates may be different than those
in J2000 epoch if you're spot-checking them.
Sirius B has a parallax of 376.68011318544626 +/- 0.4525695839121846 mas.
Using:
distance = au/parallax
This gives a distance of
8.6587097725589463725 lightyears. With the above error, this ranges from
8.648319082 ly to
8.669125460 ly
GAIA DR2 apparently does not have a radial_velocity value for Sirius B!
I even hand-searched through the gaia_source_with_rv directories.
*/
// Sirius, in J2000 coordinates. The calculations below will correct these
// positions for the epoch of the actual date.
//
starRA = HMS[6, 45, 08.91728]
starDecl = DMS[-16, 42, 58.0171]
raV = -546.01 mas/yr
declV = -1223.07 mas/yr
distance = 8.60 lightyear
radialVelocity = -5.50 km/s
sep = "\t"
preamble = ""
html = false
if length[ARGS] > 0 && ARGS@0 == "--html"
{
sep = "<TD>"
preamble = "<TR><TD>"
html = true
}
out = ["degrees", "arcmin", "arcsec"]
dateOut = ### yyyy-MM-dd hh:mm:ss a zzz ###
date = #2016#
temperature = F[59]
date = beginningOfYear[now[]]
end = beginningOfYearPlus[now[], 2]
while (date <= end)
{
// Correct the star's right ascension and declination for that date
// (this applies corrections for proper motion, precession, nutation,
// and aberration.)
[ra, decl] = apparentPosition[date, starRA, starDecl, raV, declV, distance, radialVelocity]
date = starSecantAzimuth[date, ra, decl, lat, long, corridorAzimuthMeeus, temperature, pressure]
// Now refine with temperature for that time of day.
temperature = cambridgeTemp[date]
date = starSecantAzimuth[date, ra, decl, lat, long, corridorAzimuthMeeus, temperature, pressure]
// This correction again is probably unnecessary but improves accuracy
[ra, decl] = apparentPosition[date, starRA, starDecl, raV, declV, distance, radialVelocity]
[azimuth, altitude] = refractedAzimuthAltitude[date, ra, decl, lat, long, distance, temperature, pressure]
print[preamble]
print[(date -> [dateOut, "Eastern"]) + sep]
print[format[JD[date],day,5] + sep]
print[format[altitude,degrees,2] + sep]
print[format[F[temperature],1,0] + sep]
print[skyDarkness[date, lat, long] + sep]
println[]
date = date + 1 siderealday
}
// Secant method solver to find when the a fixed star is at a given azimuth to
// an observer at the specified latitude and longitude on the Earth.
// It takes a first guess at date1 and adds an hour to that to find
// its two starting points.
//
// This function corrects for refraction.
//
// This currently doesn't handle the case where the star doesn't pass through
// the desired azimuth on that date and may give unpredictable results in
// this case or the case where the initial guess is far off.
//
// desiredAzimuth should be specified in the coordinate system specified by
// Meeus, which is angle west of south. To convert from normal true compass
// bearings, use (trueAzimuth + 180 degrees) mod circle
starSecantAzimuth[date1, starRA, starDecl, lat, long, desiredAzimuth, temperature = 283 K,
pressure=1010 millibars ] :=
{
date2 = date1 + 5 minutes
[azimuth1, altitude1] = refractedAzimuthAltitude[date1, starRA, starDecl, lat, long, undef, temperature, pressure]
[azimuth2, altitude2] = refractedAzimuthAltitude[date2, starRA, starDecl, lat, long, undef, temperature, pressure]
azimuth1 = azimuth1 - desiredAzimuth
azimuth2 = azimuth2 - desiredAzimuth
while (true)
{
// Don't overcorrect... if we're correcting more than half a phase
// from the initial guess, then go back a phase.
if azimuth1 > 180 degrees
azimuth1 = azimuth1 - circle
if azimuth1 < -180 degrees
azimuth1 = azimuth1 + circle
if azimuth2 > 180 degrees
azimuth2 = azimuth2 - circle
if azimuth2 < -180 degrees
azimuth2 = azimuth2 + circle
// did we wrap around the circle?
if abs[azimuth1 - azimuth2] > 180 degrees
if (azimuth1 < azimuth2)
azimuth2 = azimuth2 - circle
else
azimuth1 = azimuth1 - circle
// println["Azimuth 1: $azimuth1"]
// println["Azimuth 2: $azimuth2"]
// println["Altitude 1: " + (altitude1 -> "degrees")]
correction = (azimuth1 * (date1 - date2)) / (azimuth1-azimuth2)
// println["Correction: $correction"]
date = date1 - correction
// println[date]
if abs[correction] < 1 second
return date
date2 = date1
date1 = date
azimuth2 = azimuth1
[azimuth1, altitude1] = refractedAzimuthAltitude[date, starRA, starDecl, lat, long, undef, temperature, pressure]
azimuth1 = azimuth1 - desiredAzimuth
}
}
Download or view mithengeSirius.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 20217 days, 15 hours, 58 minutes ago.