// 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 http://futureboy.us/mithenge/ // // Alan Eliasen, eliasen@mindspring.com use mithengecorridor.frink use cambridgetempFourier.frink use sun.frink // 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 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] while (date <= #2018#) { // 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] 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] [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 } }