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

/*  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]

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, 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
   }
}


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 18014 days, 13 hours, 47 minutes ago.