Horizontal sundial calculations: Difference between revisions

Line 1,637:

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

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

Acos: procedure; arg x; if x<-1|x>1 then call AcosErr; return .5*pi()-Asin(x)

Atan: procedure; arg x; if abs(x)=1 then return pi()/4*sign(x)

return Asin(x/sqrt(1+x**2))

sin: procedure; arg x; x=r2r(x); numeric fuzz min(5,digits()-3)

if abs(x)=pi() then return 0; 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)

cos: procedure; arg x; x=r2r(x); a=abs(x); ~~numeric~~ ~~fuzz~~ ~~min(9,digits()-9)~~

if a=~~pi()~~ then return -1; if a=pi()/~~2 | a=pi()*2~~ 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 k=2 by 2; _=-_*x/(k*(k+i));z=z+_; if z=p then leave; p=z; end; return z

Line 1,659:

Line 1,660:

return z

pi: ~~return~~, /*a bit of an overkill, but hey !! */

pi: pi=, /*a bit of an overkill, but hey !! */

3.1415926535897932384626433832795028841971693993751058209749445923078164062862

/*Note: the real PI subroutine returns PI's accuracy that */

return pi /*Note: the real PI subroutine returns PI's accuracy that */

/*matches the current NUMERIC DIGITS, up to 1 million digits.*/

/*John Machin's formula is used for calculating more digits. */