c-----*----------------------------------------------------6---------7--
c     The discrete Fourier spectrum of one-dimensional series x(n).
c     input: n,x(n)
c       n: number of data
c       x(n): raw series
c     output: a(0:m),b(0:m),c(0:m),cta(0:m) and s(0:m)
c       a(0:m): the Fourier coefficient of sine component, m=[n/2.] 
c       b(0:m): the Fourier coefficient of cosine component, m=[n/2.] 
c       c(0:m): the amplitude spectrum = sqrt(a**2+b**2)/2
c       cta(0:m): the phase angle spectrum =atan(b/a)
c       s(0:m): the discrete power spectrum = (a**2+b**2)/2
c     By Dr. Li Jianping, May 10, 1998.    
      subroutine fourier(n,x,a,b,c,s,cta)
      dimension x(n),a(0:n),b(0:n),c(0:n),cta(0:n),s(0:n)
      pi=3.1415926
      om=2.*pi/n
      m=int(n/2.)
      an=2./n
      do 5 k=0,n
        a(k)=0.
        b(k)=0.
        c(k)=0.
        s(k)=0.
        cta(k)=0.
  5   continue
      call meanvar(n,x,ax,sx,vx)
      a(0)=ax
      b(0)=0.
      do 10 k=1,m
        a(k)=0.
        b(k)=0.
        omg=om*k
        do i=1,n
          i1=i-1
          a(k)=a(k)+x(i)*cos(omg*i1)
          b(k)=b(k)+x(i)*sin(omg*i1)
        enddo
        a(k)=a(k)*an
        b(k)=b(k)*an
  10  continue
      do 20 k=0,m
        sk=a(k)*a(k)+b(k)*b(k)
        c(k)=sqrt(sk)/2.
        s(k)=sk/2.
  20  continue
      cta(0)=0.
      do k=1,m
        cta(k)=atan(b(k)/a(k))
      enddo
      return
      end
c-----*----------------------------------------------------6---------7--
c     Computing the mean ax, standard deviation sx
c       and variance vx of a series x(i) (i=1,...,n).
c     input: n and x(n)
c       n: number of raw series
c       x(n): raw series
c     output: ax, sx and vx
c       ax: the mean value of x(n)
c       sx: the standard deviation of x(n)
c       vx: the variance of x(n)
c     By Dr. LI Jianping, May 6, 1998.
      subroutine meanvar(n,x,ax,sx,vx)
      dimension x(n)
      ax=0.
      vx=0.
      sx=0.
      do 10 i=1,n
        ax=ax+x(i)
  10  continue
      ax=ax/float(n)
      do 20 i=1,n
        vx=vx+(x(i)-ax)**2
  20  continue
      vx=vx/float(n)
      sx=sqrt(vx)
      return
      end