// 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 = "" preamble = "" 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 } }