Task

Create a program that calculates the hour, sun hour angle, dial hour line angle from 6am to 6pm for an operator entered location.



For example, the user is prompted for a location and inputs the latitude and longitude 4°57′S 150°30′W (4.95°S 150.5°W of Jules Verne's Lincoln Island, aka Ernest Legouve Reef), with a legal meridian of 150°W.

(Note: the "meridian" is approximately the same concept as the "longitude" - the distinction is that the meridian is used to determine when it is "noon" for official purposes. This will typically be slightly different from when the sun appears at its highest location, because of the structure of time zones. For most, but not all, time zones (hour wide zones with hour zero centred on Greenwich), the legal meridian will be an even multiple of 15 degrees.)

Wikipedia: A sundial is a device that measures time by the position of the Sun. In common designs such as the horizontal sundial, the sun casts a shadow from its style (also called its Gnomon, a thin rod or a sharp, straight edge) onto a flat surface marked with lines indicating the hours of the day (also called the dial face or dial plate). As the sun moves across the sky, the shadow-edge progressively aligns with different hour-lines on the plate. Such designs rely on the style being aligned with the axis of the Earth's rotation. Hence, if such a sundial is to tell the correct time, the style must point towards true north (not the north or south magnetic pole) and the style's angle with horizontal must equal the sundial's geographical latitude.




/*REXX program displays:   hour,  sun hour angle,  dial hour line angle,  6am ───► 6pm. */
numeric digits length( pi() )    -    length(.)  /*in case sundial is in polar regions. */
parse arg lat lng .                              /*obtain optional arguments from the CL*/
                              /*     ┌───────────◄ None specified?  Then use the default*/
                              /*     │             of Jules Verne's Lincoln Island,     */
$= left('', 30)               /*     ↓             aka      Ernest Legouve Reef.        */
if lat=='' | lat==","  then lat=   -4.95         /*Not specified?  Then use the default.*/
if lng=='' | lng==","  then lng= -150.5          /* "      "         "   "   "     "    */
mer= format(lng/15, , 0)  *  15                  /*calculate legal meridian longitude.  */
sineLat= sin( d2r(lat) )                         /*calculate sine of (radian) latitude. */
w1= max( length('hour'     ), length("midnight"  ))  + 2   /*compute the max hour  width*/
w2= max( length('sun hour' ), length("angle"     ))  + 2   /*   "     "   " angle    "  */
w3= max( length('dial hour'), length("line angle"))  + 2   /*   "     "   " lineº    "  */
L=  max( length(lat), length(lng), length(mer) ) /*find the maximum length of 3 numbers.*/
     say '       latitude:'    right(lat, L)     /*display the  latitude to the terminal*/
     say '      longitude:'    right(lng, L)     /*   "     "  longitude  "  "     "    */
     say ' legal meridian:'    right(mer, L)     /*   "    legal meridian "  "     "    */
     say $  center('    ', w1)    center("sun hour", w2)      center('dial hour' , w3)
     say $  center('hour', w1)    center(" angle"  , w2)      center('line angle', w3)
call sep                                         /*to help a one─eyed pirate's eyeball. */
        do h=5  to 19  by .25                    /*Okey dokey then, now let's show stuff*/
        hra= 15 * h   -   lng  +  mer - 180      /*calculate sun hour angle (in degrees)*/
        hla= r2d( Atan(sineLat * tan(d2r(hra)))) /*this is the heavy lifting calculation*/
        if abs(hra)>90  then hla= hla + 180 * sign(hra*lat)  /*adjust for negative angle*/
        call civil                               /*convert the time─of─day to civil time*/
        say $ center(hc, w1)   right(format(hra,,1), w2-2)@   right(format(hla,,1), w3-5)
        end        /*h*/
call sep                                         /*to help a one─eyed pirate's eyeball. */
exit 0                                           /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
pi:   pi= 3.1415926535897932384626433832795028841971693993751058209749445923078; return pi
d2d:  return arg(1)              // 360          /*normalize degrees ──► a unit circle. */
d2r:  return r2r( arg(1) * pi()   / 180)         /*convert degrees   ──► radians.       */
r2d:  return d2d( (arg(1) * 180   / pi() )  )    /*convert radians   ──► degrees.       */
r2r:  return arg(1)              //(pi() * 2)    /*normalize radians ──► a unit circle. */
sep:  say $  copies('═', w1)  copies("═", w2)   copies('═', w3);   @= left('',3);   return
tan:  procedure; parse arg x;  _= cos(x);   if _=0  then call tanErr;      return sin(x)/_
tellErr: say;    say '*** error ***';       say;        say arg(1);      say;      exit 13
AsinErr: call tellErr 'Asin(x),  X  must be in the range of  -1 ──► +1,  X='   x
AcosErr: call tellErr 'Acos(x),  X  must be in the range of  -1 ──► +1,  X='   x
tanErr:  call tellErr 'tan(' || x") causes division by zero,  X="              x
Acos: procedure; arg x;  if x<-1 | x>1  then call AcosErr;      return .5 * pi() - Asin(x)
Atan: procedure; parse arg x; if abs(x)=1 then return pi()*x/4; return Asin(x/sqrt(1+x*x))
/*──────────────────────────────────────────────────────────────────────────────────────*/
Asin: procedure; parse arg x;   if x<-1 | x>1  then call AsinErr;    s= x*x
      if abs(x)>=sqrt(2)*.5  then return sign(x) * Acos(sqrt(1-s));  z= x;  o= x;     p= z
        do j=2 by 2; o= o*s*(j-1)/j; z=z+o/(j+1); if z=p  then leave; p= z;  end; return z
/*──────────────────────────────────────────────────────────────────────────────────────*/
sin:  procedure; parse arg x;        x= r2r(x);          numeric fuzz min(5, digits() - 3)
                 if abs(x)=pi  then return;              return .sinCos(x, x, 1)
/*──────────────────────────────────────────────────────────────────────────────────────*/
cos:  procedure; parse arg x;        x= r2r(x);          a= abs(x);          hpi= pi * .5
                 numeric fuzz min(6, digits() - 3);      if a=pi      then return -1
                 if a=hpi | a=hpi*3  then return 0;      if a=pi/3    then return .5
                 if a=pi * 2 / 3     then return '-.5';  return .sinCos(1, 1, '-1')
/*──────────────────────────────────────────────────────────────────────────────────────*/
.sinCos: parse arg z,!,i;                  x= x*x;                     p= z
           do #=2 by 2; != -!*x/(#*(#+i)); z= z+!; if z=p  then leave; p= z; end; return z
/*──────────────────────────────────────────────────────────────────────────────────────*/
civil:      select                               /* [↓]  Maybe below  Antarctic circle, */
            when h== 0  then hc= 'midnight'      /*        or  above     Arctic    "    */
            when h <12  then hc=    h  'am'      /*convert da hour for human beans (sic)*/
            when h==12  then hc= 'noon'          /*   ···  easier to understand now.    */
            when h >0   then hc=    h  'pm'      /*   ···  even more meaningful.        */
            end   /*select*/
      parse var hc  hh ampm .;    if \datatype(hh, 'N')  then return     /*not numeric? */
      hh= hh / 1;   if hh>12  then hh= hh - 12                           /*civil time ? */
      if pos(., hh)==0  then do; hc= right(hh, 2)'   ' ampm; return; end /*exact hour ? */
      parse var  hh   hr '.' -0 mn;  if hr==0  then hr= 12               /*get MN; noon?*/
      mn= mn * 60 / 1;  hc= right(hr, 2)":"right(mn, 2, 0) ampm; return  /*reformat time*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
sqrt: procedure; parse arg x;  if x=0  then return 0;  d=digits();  numeric digits;  h=d+6
      m.=9; numeric form; parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g=g*.5'e'_ % 2
        do j=0  while h>9;      m.j= h;             h= h % 2 + 1;    end  /*j*/
        do k=j+5  to 0  by -1;  numeric digits m.k; g= (g+x/g) *.5;  end  /*k*/;  return g