/** This program calculates the coordinates of the sun as seen from earth
    in a variety of coordinate systems. */

use sun.frink

lat = 40 deg North
long = 105 deg West

d = now[]

g = new graphics

//for date = beginningOfDayPlus[d, 0] + 6 hours to beginningOfDayPlus[d, 0] + 18 hours step 1 hour
for hour = 6 hours to 18 hours step 1 hour
for date = beginningOfYearPlus[d, 0] + hour to beginningOfYearPlus[d, 1] step 1 day
{
   // Azimuth-altitude
   [az,alt] = refractedSunAzimuthAltitude[date, lat, long]

   // Change Meeus convention to have 180 deg as south.
   az = (az+180 deg) mod circle
   
   if (alt > 0 deg)
   {
      print["$date\t"]
      //   print[format[ra, deg, 5] + "\t" + format[decl, deg, 5] + "\t"]

      print[format[az mod circle, deg, 5] + "\t" + format[alt, deg, 5] + "\t"]

      // x y z angle toward the sun
      // The y axis is positive toward south.
      x = -cos[alt] sin[az]
      y = -cos[alt] cos[az]
      z = sin[alt]

      print[format[x, 1, 5] + "\t" + format[y, 1, 5] + "\t" + format[z, 1, 5]]
      
      //   print[format[lambda, deg, 5] + "\t" + format[beta, deg, 5]]
      println[]
      g.fillEllipseCenter[(az-180 deg) cos[-alt], -alt, .5 deg, .5 deg]
   }
}

g.show[]