Drag Function Subroutines (csnx, frdrag)

SUBROUTINE csnx(x,cosx,sinx,nmax)
C
C:  Iterative computation of cos(nx), sin(nx)
C
C   cos(nx) = cos[(n-1)x+x] = cos(n-1)x*cosx - sin(n-1)x*sinx
C   sin(nx) = sin[(n-1)x+x] = sin(n-1)x*cosx + cos(n-1)x*sinx
C
C   Parameter description:
C
C   x ...... variable
C   cosx ... vector of cosnx
C   sinx ... vector of sinnx
C   nmax ... max. evaluation degree
C
C   Author: Zongping Chen, GFZ/D-PAF   95/03/30
C	    with some modifications by Andrew Sinclair, RGO
C
parameter (n0=30)
implicit DOUBLE PRECISION (A-H,O-Z)
dimension cosx(n0),sinx(n0)

if (nmax.gt.n0) then
print *,' !! Fatal !! SUBR CSNX: parameter N0 too small',
+		   ' dimensioned. stop!'
stop ' ++ Error ++ '
endif
cosx(1)=cos(x)
sinx(1)=sin(x)
if (nmax.eq.1) goto 11
do 10 ii=2,nmax
cosx(ii)=cosx(ii-1)*cosx(1)-sinx(ii-1)*sinx(1)
sinx(ii)=sinx(ii-1)*cosx(1)+cosx(ii-1)*sinx(1)
10	 continue
11	 continue
return
end

Subroutine frdrag(t,x1,x2,x3,nmax,v)
C
C: Evaluation of the drag function
C
C   Drag Function:
C      v = x1 + x2*sum[(-1)**k coskX/k**2] + x3*sum[(-1)**k sinkX/k]
C	 where X = 2*pi*(t-0.5)
C
C   Parameter description:
C
C   t ........... time in days from epoch
C   x1,x2,x3 .... DRAG FRCO coefficients
C   nmax ........ degree of the Fourier expansion
C   v ........... Value after the Fourier expansion
C
C   Author: Zongping Chen, GFZ/D-PAF   95/03/30
C	    with some modifications by Andrew Sinclair, RGO
C
parameter (n0=30)
implicit DOUBLE PRECISION (A-H,O-Z)
dimension cosx(n0),sinx(n0)

if (nmax.gt.n0) then
print *,' !! Fatal !! SUBR FRDRAG: parameter N0 too small',
+		 ' dimensioned. stop!'
stop ' ++ Error ++ '
endif

pi=4.d0*atan(1.d0)
gx=2.d0*pi*(t-0.5d0)
call csnx (gx,cosx,sinx,nmax)
k=-1
sumr2i=0.d0
sumr3i=0.d0
do 20 ii=1,nmax
cosxii=cosx(ii)*dble(k)/dble(ii*ii)
sinxii=sinx(ii)*dble(k)/dble(ii)
sumr2i=sumr2i+cosxii
sumr3i=sumr3i+sinxii
k=-k
20	 continue
v=x1+x2*sumr2i+x3*sumr3i

return
end