! copyright (c) 2008 Swedish Defence Research Agency (FOI)
! may not be copied or distributed without written permission from FOI



C	rhypsmall

C	Sven Ivansson 1992




	SUBROUTINE RHYPSMALL (Y, RCOSHYP,RSINHYP)
C	to compute  rcoshyp = coshyp(y)
C	            rsinhyp = sinhyp(y)
C	accurately, even for small y

	IMPLICIT NONE

	REAL*16 Y
	REAL*16 RCOSHYP,RSINHYP

	INTEGER N, I
	REAL*16 EPS
	REAL*16 Y2, SLASK1,SLASK2, RCREL(50), RSREL(50)

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

	IF (ABS(Y).GT.0.1Q0) THEN
C	  if y were very small, num. canc. would occur for rsinhyp here
	  SLASK1 = EXP(Y)
	  SLASK2 = 1/SLASK1
	  RCOSHYP = 0.5Q0 * ( SLASK1 + SLASK2 )
	  RSINHYP = 0.5Q0 * ( SLASK1 - SLASK2 )
	ELSE
C	  rcoshyp  =       1 + y**2/2 + y**4/24  + ...
C	  rsinhyp  = y * ( 1 + y**2/6 + y**4/120 + ... )
	  Y2 = Y*Y
	  N = 1
	  I = 2
	  RCREL(1) = Y2/2
	  RSREL(1) = Y2/6
	  DO WHILE (RCREL(N).GT.EPS)
	    N = N + 1
	    I = I + 2
	    IF (N.GT.50) STOP '%%% rhypsmall, dim'
	    RCREL(N) = RCREL(N-1) * Y2 / ((I-1)*I)   !  y**(2*n) / ((2*n)!)
	    RSREL(N) = RCREL(N) / (I+1)              !  y**(2*n) / ((2*n+1)!)
	  END DO
	  RCOSHYP = 0.0Q0
	  DO I = N, 1, -1
	    RCOSHYP = RCOSHYP + RCREL(I)
	  END DO
	  RCOSHYP = 1.0Q0 + RCOSHYP
	  RSINHYP = 0.0Q0
	  DO I = N, 1, -1
	    RSINHYP = RSINHYP + RSREL(I)
	  END DO
	  RSINHYP = Y * ( 1.0Q0 + RSINHYP )

	END IF

	RETURN
	END
