/** This program generates a digital sundial. */

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/2 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 = r * cos[alt] sin[az]
      // y = r * cos[alt] cos[az]
      // z = r *  sin[alt]
      //
      // Inverse:
      // r   = sqrt[x^2 + y^2 + z^2]
      // alt = arcsin[z / r]    (if r != 0)
      // az  = arctan[x,y]
      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, 1 deg, 1 deg]
   }
}

g.show[]