! Trig package
! JWM August 1982

%option "-nocheck-nodiag-noline"

%const %real PIO4   = 0.7853981633974483096156608
%const %real PIO2   = 1.570796326794896619231322
%const %real PI     = 3.141592653589793238462643
%const %real PI2    = 6.283185307179586476925287
%const %real E      = 2.718281828

%real %function NORMALISE(%real x)
{brings result into range  -pi <= result <= +pi}
  x = x-int pt(x/pi2)*pi2
  %result = pi2+x %if x <= -pi
  %result = x-pi2 %if x >= pi
  %result = x
%end {of normalise}

%real %function RAW SINE(%real x)
{parameter must be between 0 and pi/4}
%real y;  y = x*x
  %result = (((0.0032811761*y-0.1335639326)*y+1.0)*x)/ %c
            ((0.0004649838*y+0.0331027317)*y+1.0)
%end {of raw sine}

%real %function RAW COSINE(%real x)
{parameter must be between 0 and pi/4}
%real y;  y = x*x
  %result = ((0.0205121130*y-0.4558922221)*y+0.9999999992)/ %c
            ((0.0008996261*y+0.0441077396)*y+1.0)
%end {of raw cosine}

%external %real %function SIN(%real x)
  %result = -sin(-x) %if x < 0
  x = normalise(x)
  x = pi-x %if x > pio2
  %result = raw cosine(pio2-x) %if x > pio4
  %result = raw sine(x)
%end {of sin}

%external %real %function COS(%real x)
  x = normalise(x)
  x = -x %if x < 0
  %if x > pio2 %start
    x = pi-x
    %result = -raw sine(pio2-x) %if x > pio4
    %result = -raw cosine(x)
  %finish %else %start
    %result = raw sine(pio2-x) %if x > pio4
    %result = raw cosine(x)
  %finish
%end {of cos}

%external %real %function TAN(%real x)
{%real s
{  s = sin(x)
{  %result = s/sqrt(1.0-s*s)
  %result = sin(x)/cos(x)
%end {of tan}

%real %function A TAN(%real x)
! parameter must be in range 0 <= x <= 1
%real y;  y = x*x
  %result = %c
  x*(((0.0089472229*y+0.2870044785)*y+1.1303754276)*y+0.9999999992)/ %c
    (((0.0506770959*y+0.5749098994)*y+1.4637086946)*y+1.0)
%end

%external %real %function P ARC TAN(%real x)
  %result = -p arc tan(-x) %if x < 0
  %result = a tan(x) %if 0 <= x <= 1.0
  %result = pio2 - a tan(1.0/x)
%end {of arc tan1}

%external %real %function ARC TAN(%real x,y)
! Returns the value, in radians, of the angle whose tangent
! is specified by "y/x". This value is between +&- pi. If
! x<0 then it is in the second or third quadrants otherwise
! it is in the first or fourth quadrants.
%real z
   z = p arc tan(y/x)       {don't have to worry about -ve case}
   %result = z %if x >= 0
   %result = pi-z %if z > 0
   %result = -pi-z
%end {of arc tan}

%external %real %function ARC SIN(%real x)
  %if x >= 1.0 %start
    %signal 10 %if x > 1.0
    %result = pio2
  %finish
  %if x <= -1.0 %start
    %signal 10 %if x < 1.0
    %result = -pio2
  %finish
  %result = p arc tan(x/(sqrt(1.0-x*x)))
%end {of arc sin}

%external %real %function ARC COS(%real x)
%real w
  %if x > 0.999 %start
    %signal 10 %if x > 1.0
    %result = 0
  %finish
  %if x < -0.999 %start
    %signal 10 %if x < -1.0
    %result = pi
  %finish
  %result = pio2 %if -0.001 < x < 0.001
  w = p arc tan(sqrt(1.0/(x*x)-1.0))
  {this next because arc tan is in range between +&-pi}
  %result = w %if x >= 0
  %result = -w+pi
%end {of arc cos}

%real %function LOG E(%real z)
! z is such that 0.5 <= z <= 1
%real x,y
  x = z+z-1.0                   {nasty VAX compiler bug}
  y = x*x
  %result = ((0.0956558162*x+0.5297501385)*y+0.0677412133*x-0.6931471773) %c
              /((0.0286818192*x+0.45477291277)*y+1.3449644663*x+1.0)
%end {of log e}

%external %real %function LOG(%real x)
! logarithm to base e
%const %real log of two = 0.6931471806
%integer factor;  factor = 0
  %signal 10 %if x <= 0
  %if x > 1.0 %start
    %cycle
      factor = factor+1;  x = x/2
    %repeat %until 0.5 <= x <= 1.0
    %result = log e(x) + factor*log of two
  %finish %else %start
    %result = log e(x) %if x >= 0.5
    %cycle
      factor = factor+1;  x = x+x
    %repeat %until 0.5 <= x <= 1.0
    %result = log e(x) - factor*log of two
  %finish
%end {of log}

%external %real %function LOG TEN(%real x)
! logarithm to base 10
  %result = log(x)/log(10.0)
%end {of log10}

%real %function RAW EXP(%real x)
! x must be in the range 0 <= x <= 1 (to give e to the power x)
%real y;  y = x*x
  %result = ((0.0106337905*x+0.1125548636)*y+0.5240642207*x+1.0)/ %c
            ((0.0884921370-0.0065658101*x)*y+1.0000000007-0.4759358618*x)
%end {of raw exp}

%external %real %function EXP(%real x)
! e to the power x
%real temp
  %if -1.0 <= x <= 1.0 %start
    %result = raw exp(x) %if x >= 0
    %result = 1.0/raw exp(-x)
  %finish %else %start
    temp = exp(0.5*x);  %result = temp*temp
  %finish
%end {of exp}

%external %real %function SIN H(%real x)
  %result = 0.5*(exp(x)-exp(-x))
%end {of sin h}

%external %real %function COS H(%real x)
  %result = 0.5*(exp(x)+exp(-x))
%end {of cos h}

%external %real %function TAN H(%real x)
  %result = sin h(x)/cos h(x)
%end {of tan h}

%end %of %file