
C	circrest

C	Sven Ivansson 1990




	SUBROUTINE CIRCREST (Z, CCIRCR,STCIRCR,SDCIRCR)
C	to compute  ccircr  = cos(z)
C	            stcircr = sin(z) * z
C	            sdcircr = sin(z) / z
C	accurately, even for small z

	IMPLICIT NONE

	COMPLEX*16 Z
	COMPLEX*16 CCIRCR,STCIRCR,SDCIRCR

	INTEGER N, I
	REAL*8 EPS
	COMPLEX*16 Z2, CSLASK, CREL(20), SREL(20)

	DATA EPS /1.0D-17/  ! desired relative accuracy (approximately), 
C	strict analysis possible with MacLaurin remainder term

	IF ( DMAX1(DABS(DREAL(Z)),DABS(DIMAG(Z))) .GT. 0.5D0 ) THEN
C	  more economical by using Euler and "cdexp(-)=1/cdexp(+)"
C	  if z were very small, num. canc. would occur for "cdsin(z)" since 
C	  this seems to be computed (on Vax) using Euler's formulas when z 
C	  is nonreal
	  CCIRCR = CDCOS(Z)
	  CSLASK = CDSIN(Z)
	  STCIRCR = CSLASK * Z
	  SDCIRCR = CSLASK / Z
	ELSE
C	  ccircr  =            1 - z**2/2 + z**4/24  + ...
C	  stcircr = (z**2) * ( 1 - z**2/6 + z**4/120 + ... )
C	  sdcircr =            1 - z**2/6 + z**4/120 + ...
	  Z2 = Z*Z
	  N = 1
	  I = 2
	  CREL(1) = -Z2/2
	  SREL(1) = -Z2/6
	  DO WHILE ( DMAX1(DABS(DREAL(CREL(N))),DABS(DIMAG(CREL(N)))) 
     $	  .GT. EPS )
	    N = N + 1
	    I = I + 2
	    CREL(N) = -CREL(N-1) * Z2 / ((I-1)*I)
C	                                     !  (-1)**n * z**(2*n) / ((2*n)!)

	    SREL(N) = CREL(N) / (I+1)        !  (-1)**n * z**(2*n) / ((2*n+1)!)
	  END DO
	  CCIRCR = (0.0D0,0.0D0)
	  DO I = N, 1, -1
	    CCIRCR = CCIRCR + CREL(I)
	  END DO
	  CCIRCR = (1.0D0,0.0D0) + CCIRCR
	  SDCIRCR = (0.0D0,0.0D0)
	  DO I = N, 1, -1
	    SDCIRCR = SDCIRCR + SREL(I)
	  END DO
	  SDCIRCR = (1.0D0,0.0D0) + SDCIRCR

	  STCIRCR = Z2 * SDCIRCR
	END IF

	RETURN
	END
