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 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 = "<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 17645 days, 8 hours, 31 minutes ago.