Horizontal sundial calculations: Difference between revisions

Content deleted Content added
No edit summary
m →‎{{header|REXX}}: added/changed whitespace and comments, simplified some functions, corrected some typos.
Line 1,590: Line 1,590:


No attempt was made to explain the inner workings of the trigonometric functions.
No attempt was made to explain the inner workings of the trigonometric functions.
<lang rexx>/*REXX program shows hour/sun hour angle/dial hour line angle, 6am─►6pm.*/
<lang rexx>/*REXX pgm shows: hour, sun hour angle, dial hour line angle, 6am ───► 6pm*/
numeric digits 60 /*better overkill then underkill. */
numeric digits 60 /*better overkill then underkill for dig*/


parse arg lat lng mer . /*get the arguments (if any). */
parse arg lat lng mer . /*get the optional arguments from the CL*/
/*If none specified, then use the */
/*None specified? Then use the default of Jules */
/*default of Jules Verne's Lincoln*/
/*Verne's Lincoln Island, aka Ernest Legouve Reef.*/
/*Island, aka Ernest Legouve Reef.*/


if lat=='' | lat==',' then lat=-4.95 /*No argument? Then use default.*/
if lat=='' | lat==',' then lat=-4.95 /*Not specified? Then use the default.*/
if lng=='' | lng==',' then lng=-150.5 /*No argument? Then use default.*/
if lng=='' | lng==',' then lng=-150.5 /* " " " " " " */
if mer=='' | mer==',' then mer=-150 /*No argument? Then use default.*/
if mer=='' | mer==',' then mer=-150 /* " " " " " " */
L=max(length(lat),length(lng),length(mer))
L=max(length(lat), length(lng), length(mer))
say ' latitude:' right(lat,L)
say ' latitude:' right(lat,L)
say ' longitude:' right(lng,L)
say ' longitude:' right(lng,L)
say ' legal meridian:' right(mer,L)
say ' legal meridian:' right(mer,L)
sineLat=sin(d2r(lat))
sineLat=sin(d2r(lat))
w1=max(length('hour'),length('midnight'))+2
w1=max(length('hour') ,length('midnight'))+2
w2=max(length('sun hour') ,length('angle'))+2
w2=max(length('sun hour') ,length('angle'))+2
w3=max(length('dial hour'),length('line angle'))+2
w3=max(length('dial hour'),length('line angle'))+2
indent=left('',30) /*make presentation prettier. */
indent=left('',30) /*make the presentation prettier. */
say indent center(' ',w1) center('sun hour',w2) center('dial hour' ,w3)
say indent center(' ',w1) center('sun hour',w2) center('dial hour' ,w3)
say indent center('hour',w1) center('angle' ,w2) center('line angle',w3)
say indent center('hour',w1) center('angle' ,w2) center('line angle',w3)
call sep /*add separator line for eyeballs*/
call sep /*add a separator line for the eyeballs*/


do h=-6 to 6 /*Okey dokey then, let's get busy*/
do h=-6 to 6 /*Okey dokey then, let's get busy. */
select
select
when abs(h)==12 then hc='midnight' /*above artic circle ? */
when abs(h)==12 then hc='midnight' /*above the arctic circle? */
when h<0 then hc=-h 'am' /*convert hour for human beans. */
when h<0 then hc=-h 'am' /*convert the hour for human beans. */
when h==0 then hc='noon' /* ... easy to understand now. */
when h==0 then hc='noon' /* ... easier to understand now. */
when h>0 then hc=h 'pm' /* ... even meaningfull. */
when h>0 then hc=h 'pm' /* ... even more meaningful. */
end /*select*/
end /*select*/
hra=15*h-lng+mer
hra=15*h-lng+mer
Line 1,627: Line 1,626:


call sep
call sep
exit /*stick a fork in it, we're done.*/
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────subroutines───────────────────────────────*/

sep: say indent copies('═',w1) copies('═',w2) copies('═',w3); return
/*──────────────────────────────────subroutines─────────────────────────*/
d2d: return arg(1) // 360 /*normalize degrees ──► a unit circle. */
/*looking at subroutines is like looking at saugages being made. Don't.*/
d2r: return r2r(arg(1)*pi() / 180) /*convert degrees ──► radians. */
sep: say indent copies('═',w1) copies('═',w2) copies('═',w3); return
pi: return, /*a bit of overkill, but hey !! */
r2d: return d2d((arg(1)*180 / pi())) /*convert radians ──► degrees. */
r2r: return arg(1) // (2*pi()) /*normalize radians ──► a unit circle. */
3.1415926535897932384626433832795028841971693993751058209749445923078164062862
tan: procedure; parse arg x; _=cos(x); if _=0 then call tanErr; return sin(x)/_
/*Note: the real PI subroutine returns PI's accuracy that */
tanErr: call tellErr 'tan(' || x") causes division by zero, X=" || x
/*matches the current NUMERIC DIGITS, up to 1 million digits.*/
/*John Machin's formula is used for calculating more digits. */

d2d: return arg(1)//360 /*normalize degrees►1 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)//(2*pi()) /*normalize radians►1 unit circle*/
tan: procedure; arg x; _=cos(x); if _=0 then call tanErr; return sin(x)/_
tanErr: call tellErr 'tan('||x") causes division by zero, X="||x
tellErr: say; say '*** error! ***'; say; say arg(1); say; exit 13
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
AsinErr: call tellErr 'Asin(x), X must be in the range of -1 ──► +1, X=' ||x
sqrtErr: call tellErr "sqrt(x), X can't be negative, X="||x
sqrtErr: call tellErr "sqrt(x), X can't be negative, X="||x
AcosErr: call tellErr 'Acos(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
Acos: procedure; arg x; if x<-1|x>1 then call AcosErr; return .5*pi()-Asin(x)
Acos: procedure; arg x; if x<-1|x>1 then call AcosErr; return .5*pi()-Asin(x)


Line 1,661: Line 1,652:
return .sincos(1,1,-1)
return .sincos(1,1,-1)
.sinCos: parse arg z,_,i; x=x*x; p=z
.sinCos: parse arg z,_,i; x=x*x; p=z
do k=2 by 2; _=-_*x/(k*(k+i));z=z+_; if z=p then leave;p=z;end; return z
do k=2 by 2; _=-_*x/(k*(k+i));z=z+_; if z=p then leave; p=z; end; return z


Asin: procedure; parse arg x; if x<-1 | x>1 then call AsinErr; s=x*x
Asin: procedure; parse arg x; if x<-1 | x>1 then call AsinErr; s=x*x
Line 1,668: Line 1,659:
return z
return z


pi: return, /*a bit of an overkill, but hey !! */
sqrt: procedure; parse arg x; if x=0 then return 0;d=digits();numeric digits 11
3.1415926535897932384626433832795028841971693993751058209749445923078164062862
g=.sqrtGuess(); do j=0 while p>9; m.j=p; p=p%2+1; end
/*Note: the real PI subroutine returns PI's accuracy that */
do k=j+5 to 0 by -1; if m.k>11 then numeric digits m.k; g=.5*(g+x/g); end
/*matches the current NUMERIC DIGITS, up to 1 million digits.*/
numeric digits d; return g/1
/*John Machin's formula is used for calculating more digits. */
.sqrtGuess: if x<0 then call sqrtErr; numeric form scientific; m.=11; p=d+d%4+2

parse value format(x,2,1,,0) 'E0' with g 'E' _ .; return g*.5'E'_%2</lang>
sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); i=; m.=9
{{out}}
numeric digits 9; if x<0 then do; x=-x; i='i'; end; numeric form; h=d+2
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=.5*(g+x/g); end /*k*/
numeric digits d; return (g/1)i</lang>
'''pit[it'''
<pre>
<pre>
latitude: -4.95
latitude: -4.95