mithengeSirius.frink

View or download 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

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


View or download 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 17929 days, 23 hours, 27 minutes ago.