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