MODULE  SINCOS;  { Copyright (c)  by  T. W. Lougheed  24 April 1981 }

{ This procedure provides a replacement for the SINE function currently
available to the pascal system.  It is faster and more accurate. }

{ First version   5 February 1981

  By    T. W. Lougheed
        Dept. T. & A. Mechanics
        Thurston Hall, Cornell U.
        Ithaca,  NY  14853

  Last version  23 February 1981

  This software is in the public domain, and may not be sold by any
  person or corperation without permission of the author.  }



FUNCTION SIN( Z :REAL) :REAL;

CONST  TWO_PI  = 6.283185307179586;
       PI      = 3.141592653589793;
       HALF_PI = 1.570796326794897;

{ The following is from the HANDBOOK OF MATHEMATICAL FUNCTIONS, by
Abramowicz and Stegun, tenth printing.  Formula 4.3.97;  its error
is less than  |X| 10^8  for  |X| < PI/2.  }

FUNCTION  TCHEBYSHEV( X :REAL) :REAL;

   CONST A2  = -0.1666666664;
         A4  =  0.0083333315;
         A6  = -0.0001984090;
         A8  =  0.0000027526;
         A10 = -0.0000000239;

   VAR  S  :REAL;
   BEGIN
   S := SQR( X );
   TCHEBYSHEV := (((((A10*S + A8)*S + A6)*S + A4)*S + A2)*S + 1)*X;
   END;


BEGIN

{ Map the argument  Z  onto the interval  -PI..PI, rounding
towards the smallest absolute value. }

IF Z > 0 THEN  Z := Z - TWO_PI*ROUND(  Z/TWO_PI )
         ELSE  Z := Z + TWO_PI*ROUND( -Z/TWO_PI );

{ Always ask for  TCHEBYSHEV  of a number in the interval -PI/2..PI/2. }

IF  Z >  HALF_PI  THEN  SIN := TCHEBYSHEV( PI - Z )
ELSE IF  Z > -HALF_PI  THEN  SIN := TCHEBYSHEV( Z )
ELSE  SIN :=  -TCHEBYSHEV( PI + Z );

END;



{ Lazy cosine function based on the identity

    cos x  =  ( 1 - (sin x/2)^2 )/2

note that this has been re-arranged somewhat to
improve accuracy. }

FUNCTION COS( Z :REAL) :REAL;
   CONST SQRT_2 = 1.4142135623730950;
   VAR   S  :REAL;
   BEGIN
   S := SQRT_2*SIN( Z/2 );
   COS := (1 - S)*(1 + S);
   END;


MODEND
.