
C	cossinc

C	Sven Ivansson 2020




	COMPLEX*16 FUNCTION COSSINC (Z)
C	to compute  cossinc = ( cos(z) - sin(z)/z ) / z**2
C	accurately, even for small z

	IMPLICIT NONE

	COMPLEX*16 Z

	INTEGER N, I
	REAL*8 EPS
	COMPLEX*16 Z2, CREL(20), CSLASK1,CSLASK2, EIMULT

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

	Z2 = Z*Z

	IF ( MAX(ABS(REAL(Z)),ABS(AIMAG(Z))) .GT. 0.5D0 ) THEN
	  CSLASK1 = EIMULT(Z)
	  CSLASK1 = EXP(CSLASK1)
	  CSLASK2 = 1/CSLASK1
	  COSSINC = ( (CSLASK1+CSLASK2) - (CSLASK1-CSLASK2)/EIMULT(Z) ) / 
     &	   (2*Z2)
	ELSE
C	  cossinc = -1/3 + z**2/30 - ...
	  N = 1
	  CREL(1) = -1.0D0/3
	  DO WHILE ( MAX(ABS(REAL(CREL(N))),ABS(AIMAG(CREL(N)))) .
     $	  GT. EPS )
	    N = N + 1
	    IF (N.GT.20) STOP '%%% cossinc, crel dim'
	    I = 2*N  ! crel(n) = (-1)**n * z**(2*n-2) / ((2*n-1)!*(2*n+1))
	    CREL(N) = -CREL(N-1) * Z2 / ((I-2)*(I+1))
	  END DO
	  COSSINC = (0.0D0,0.0D0)
	  DO I = N, 1, -1
	    COSSINC = COSSINC + CREL(I)
	  END DO
	END IF

	RETURN
	END
