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 19966 days, 2 hours, 55 minutes ago.