c-----*----------------------------------------------------6---------7--
c     Subroutine for continuous cross spectrum analysis of two 
c         one-dimensional series x(i) and y(i) (i=1,...,n). 
c     Input parameters and arrays: n,m, x(n), y(n),
c       n: number of data
c       m: biggest lag time length. Generally, m is between n/3 and n/10.
c       x(n),y(n): raw series
c     Output variables:
c        ol(0:m): frequency array
c        tl(0:m): periodic array  
c        px(0:m): power spectrum of x(n)
c        py(0:m): power spectrum of y(n)
c       pxy(0:m): co-spectrum between x and y
c       qxy(0:m): quad-spectrum between x and y
c       rxy(0:m): coherence spectrum between x and y (it is between 0 and 1) 
c       cxy(0:m): phase difference spectrum between x and y
c            If cxy(i)>0, then x trails y for the wave i, i.e. y leads x. 
c       lxy(0:m): lag length spectrum between x and y
c       rxy951(0:m): 95% F-test confidence upper limit for coherence spectrum rxy
c       rxy952(0:m): 95% Goodman-test confidence upper limit for coherence spectrum rxy 
c     By Dr. Li Jianping, September 6, 2001.    
      subroutine ccrossspectrum(n,m,x,y,ol,tl,px,py,px95,py95,
     *                          rxy,cxy,lxy,rxy951,rxy952)
      dimension x(n),y(n)
      dimension px(0:m),py(0:m),ol(0:m),tl(0:m)
      real pxy(0:m),qxy(0:m),rxy(0:m),cxy(0:m),lxy(0:m)
      real rxy951(0:m),rxy952(0:m),px95(0:m),py95(0:m)
      real strwx(0:m),strwy(0:m)
      dimension r12(0:m),r21(0:m)   ! Work array
      pi=3.1415926
c---  It must be pointed out that if you use the subroutine lagcorrelation, you must
c        call the subroutine autocrrelation in cspectrum, and that if you use
c        the lagcorrelation1, the autocorrelation1 must be used in cspectrum!!!
      call lagcorrelation(n,m,x,y,r12)
      call lagcorrelation(n,m,y,x,r21)
c-----
      call cspectrum(n,m,x,ol,tl,px,px95,strwx)
      call cspectrum(n,m,y,ol,tl,py,py95,strwy)
      do 10 l=0,m
        pxy(l)=0.
        qxy(l)=0.
        do i=1,m-1 
          cc=pi*i/float(m)
          pxy(l)=pxy(l)+(r12(i)+r21(i))*(1.+cos(cc))*cos(l*cc)
          qxy(l)=qxy(l)+(r12(i)-r21(i))*(1.+cos(cc))*sin(l*cc)
        enddo
        pxy(l)=(r12(0)+pxy(l)/2.)/float(m)
        qxy(l)=(qxy(l)/2.)/float(m)
  10  continue
      pxy(0)=pxy(0)/2.
      pxy(m)=pxy(m)/2.
      qxy(0)=qxy(0)/2.
      qxy(m)=qxy(m)/2.
      do 20 l=1,m
        ol(l)=pi*l/m
        tl(l)=2.*m/l
  20  continue
      ol(0)=0
      tl(0)=100.*tl(1)         !In fact, tl(0) is infinite.
      do 30 l=0,m
        if(px(l).le.0..or.py(l).le.0.)then
          rxy(l)=0.
        else
          rxy(l)=(pxy(l)**2+qxy(l)**2)/(px(l)*py(l))
        endif
        cxy(l)=atan(qxy(l)/pxy(l))
        lxy(l)=cxy(l)*tl(l)/(2.*pi)
  30  continue
      call crossspectrumtest(n,m,rxy951,rxy952)
      return
      end
c-----*----------------------------------------------------6---------7--
c     Calculating cross lag correlation rt(0:m) between x(i) and y(i).
c     m: input parameter, maximum lag time
c     By Dr. LI Jianping, Spetember 6, 2001.
      subroutine lagcorrelation(n,m,x,y,rt)
      dimension x(n),y(n),rt(0:m)
      dimension x1(n),y1(n)               ! Work Array
c     Lagged correlations: y(i) trails x(i).
      call meanvar(n,x,ax,sx,vx)
      call meanvar(n,y,ay,sy,vy)
      call correlation(n,x,y,r)
      rt(0)=r
      do i=1,n
        x(i)=x(i)-ax
        y(i)=y(i)-ay
      enddo
      do 10 i=1,m
        sxy=0.
        do j=1,n-i
          sxy=sxy+x(j)*y(i+j)
        enddo
        sxy=sxy/float(n)
        rt(i)=sxy/(sx*sy)
  10  continue
      return
      end
c-----*----------------------------------------------------6---------7--
c     For the correlation coefficient r between two series 
c       x(i) and y(i), where i=1,...,n.
c     input: n,x(n),y(n)
c       n: number of time series
c       x(n): raw series 
c       y(n): raw series 
c     output: r
c       r: correlation coefficient between x and y
c     By Dr. LI Jianping, January 5, 2000.
      subroutine correlation(n,x,y,r)
      dimension x(n),y(n)
      call meanvar(n,x,ax,sx,vx)
      call meanvar(n,y,ay,sy,vy)
      sxy=0.
      do 10 i=1,n
        sxy=sxy+(x(i)-ax)*(y(i)-ay)
  10  continue
      sxy=sxy/float(n)
      r=sxy/(sx*sy)
      return
      end
c-----*----------------------------------------------------6---------7--
c     Subroutine for continuous spectrum analysis of an one-dimensional series
c         x(i) (i=1,...,n). 
c     Input parameters and arrays: n,m,lw, x(n),
c       n: number of data
c       m: biggest lag time length. Generally, m is between n/3 and n/10.
c       lw: maximum wave number (=[n/2]).
c       x(n): raw series
c     Output variables:
c        ol(0:m): frequency array
c        tl(0:m): periodic array  
c        sl(0:m): power spectrum. Its sum from 0 to m is 1.
c      st95(0:m): 95% confidence upper limit of red or white noise spectrum
c      strw(0:m): spectrum density of red or white noise
c     By Dr. Li Jianping, May 12,1998.    
      subroutine cspectrum(n,m,x,ol,tl,sl,st95,strw)
      dimension x(n)
      dimension sl(0:m),ol(0:m),tl(0:m),st95(0:m),strw(0:m)
      dimension r(0:m)            ! Work array
      pi=3.1415926
      call autocorrelation(n,m,x,r)
      do 10 l=0,m
        sl(l)=r(0)
        do i=1,m-1 
          cc=pi*i/float(m)
          sl(l)=sl(l)+r(i)*(1.+cos(cc))*cos(l*cc)
        enddo
  10  continue
      sl(0)=sl(0)/2.
      sl(m)=sl(m)/2.
      do 20 l=0,m
        sl(l)=sl(l)/m
  20  continue
      do 30 l=1,m
        ol(l)=pi*l/m
        tl(l)=2.*m/l
  30  continue
      ol(0)=0
      tl(0)=100.*tl(1)
      call cnoisetest(n,m,r(1),sl,st95,strw)
      return
      end
c-----*----------------------------------------------------6---------7--
c     This is a subroutine for calculating the lag autocorrelations
c        r(m) of a series x(i).
c     input: n,m,and x(n)
c        n: number of the data x(n)
c        x: the raw series
c        m: maximum lagged length;  m: maximum lag time
c     output: r(0:m)
c        r: lag autocorrelations from 0 to m.
c     By Dr. LI Jianping, April 3, 1998.
      subroutine autocorrelation(n,m,x,r)
      dimension x(n),r(0:m)
      call meanvar(n,x,ax,sx,vx)
      r(0)=1.
      do i=1,n
        x(i)=x(i)-ax
      enddo
      do 10 i=1,m
        sxy=0.
        do j=1,n-i
          sxy=sxy+x(j)*x(i+j)
        enddo
        sxy=sxy/float(n)
        r(i)=sxy/vx
  10  continue
      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
c-----*----------------------------------------------------6---------7--
c     Noise test for continuous spectrum (95% confidence level).
c     strw: spectrum density of red or white noise
c     chi**2 distribution with degrees of freedom v=(2*n-m/2)/m is used.
c     By Dr. Li Jianping, May 12,1998.    
      subroutine cnoisetest(n,m,r1,sl,st95,strw)
      dimension sl(0:m),st95(0:m),strw(0:m)
      dimension chi95(30),chi99(30)
      call chi_table(chi95,chi98,chi99)
      pi=3.1415926
      as=0.
      do 10 l=0,m
        as=as+sl(l)
  10  continue
      as=as/float(m+1)
      r2=r1*r1
      v=(2.*n-m/2.)/float(m)
      if(v.gt.30.)then  !! invaliable from mathmatic table for v>30
        write(*,*)'Unable continue... Please reinput m !'
        stop
      endif
      iv=int(v)
      c95=chi95(iv)+(v-iv)*(chi95(iv+1)-chi95(iv))
      if(r1.gt.0.)then
c       tname='red noise test'
        do l=0,m
          r3=1+r2-2*r1*cos(l*pi/m)
          strw(l)=as*(1.-r2)/r3
          st95(l)=strw(l)*c95/v
        enddo
      else
c       tname='white noise test'
        do l=0,m
          strw(l)=as
          st95(l)=as*c95/v
        enddo
      endif
      return
      end 
c-----*----------------------------------------------------6---------7--
c     chi-square distribution
c     Chi-table with right tail probabilities
c       P(X**2>=Xa**2)=a
c       where a (alpha) significance level and (1-a)*100% is confindence level.
c     chi95: chi-square distribution at 95% confidence level.
c     chi98: chi-square distribution at 98% confidence level.
c     chi99: chi-square distribution at 99% confidence level.
c     By Dr. Li Jianping, May 12,1998.    
      subroutine chi_table(chi95,chi98,chi99)
      dimension chi95(30),chi98(30),chi99(30)
      dimension c95(30),c98(30),c99(30) !work array
      data c95/ 3.841, 5.991, 7.815, 9.488,11.070,
     *         12.592,14.067,15.507,16.919,18.307,
     *         19.675,21.026,22.362,23.685,24.996,
     *         26.296,27.587,28.869,30.144,31.410,
     *         32.671,33.924,35.172,36.415,37.652,
     *         38.885,40.113,41.337,42.557,43.773/
      data c98/ 5.412, 7.824, 9.837,11.668,13.388,
     *         15.033,16.622,18.168,19.679,21.161,
     *         22.618,24.054,25.472,26.873,28.259,
     *         29.633,30.995,32.346,33.687,35.020,
     *         36.343,37.659,38.968,40.270,41.566,
     *         42.856,44.140,45.419,46.693,47.662/
      data c99/ 6.635, 9.210,11.345,13.277,15.086,
     *         16.812,18.475,20.090,21.666,23.209,
     *         24.725,26.217,27.688,29.141,30.578,
     *         32.000,33.409,34.805,36.191,37.566,
     *         38.932,40.289,41.638,42.980,44.314,
     *         45.642,46.963,48.278,49.588,50.892/
      do i=1,30
        chi95(i)=c95(i)
        chi98(i)=c98(i)
        chi99(i)=c99(i)
      enddo
      return
      end
c-----*----------------------------------------------------6---------7--
c     Significant test for continuous cross spectrum (95% confidence level).
c     There are two test methods.
c     rxy951: F-test (=sqrt(F/(F+v-1)), where degree of freedom v=(2*n-m/2)/m).
c             F distribution here is with degrees of freedom numerator 2 and 
c               degrees of freedom denominator 2*(v-1).
c     rxy952: Goodman-test (=sqrt(1-a**(1/(v-1)), where a=0.05 is significant level)
c             The degrees of freedom v=(2*n-m/2)/m).
c     By Dr. Li Jianping, September 6, 2001.    
      subroutine crossspectrumtest(n,m,rxy951,rxy952)
      dimension rxy951(0:m),rxy952(0:m)
      dimension fsd95(1000)
      call F_table(fsd95)
      pi=3.1415926
      v=(2.*n-m/2.)/float(m)
      v1=2.*(v-1)
      iv=int(v1)
      if(iv.gt.1000)then
        f95=fsd95(1000)
      else
        f95=fsd95(iv)+(v1-iv)*(fsd95(iv+1)-fsd95(iv))
      endif
      r951=f95/(f95+v-1)
      r952=1.-0.05**(1./(v-1.))
      do 10 i=0,m
        rxy951(i)=r951
        rxy952(i)=r952
 10   continue
      return
      end
c-----*----------------------------------------------------6---------7--
c     The F distribution is a right-skewed distribution used most commonly in 
c         Analysis of Variance. The F distribution is a ration of two Chi-square
c         distributions, and a specific F distribution is denoted by the degrees
c         of freedom for thw numerator Chi-square and the degrees of freedom
c         for the denominator Chi-square.
c         F=S1**2/S2**2,
c         P(F>=Fa)=0.05
c     F(n1,n2): F distribution at 95% confidence level.
c       n1 the numerator degrees of freedom
c       n2 the denominator degrees of freedom
c     Here fsd95(n)=F(2,n) for discrete power spectrum test.
      subroutine F_table(fsd95)
      dimension fsd95(1000)
      dimension fsd(30)
      data fsd/200.00,19.00,9.55,6.94,5.79,5.14,4.74,4.46,4.26,4.10,
     *           3.98, 3.89,3.81,3.74,3.68,3.63,3.59,3.55,3.52,3.49,
     *           3.47, 3.44,3.42,3.40,3.39,3.37,3.35,3.34,3.33,3.32/
      do i=1,30
        fsd95(i)=fsd(i)
      enddo
      do i=31,40
        fsd95(i)=3.32-(3.32-3.23)*(i-30.)/10.
      enddo
      do i=41,50
        fsd95(i)=3.23-(3.23-3.18)*(i-40.)/10.
      enddo
      do i=51,60
        fsd95(i)=3.18-(3.18-3.15)*(i-50.)/10.
      enddo
      do i=61,70
        fsd95(i)=3.15-(3.15-3.13)*(i-60.)/10.
      enddo
      do i=71,80
        fsd95(i)=3.13-(3.13-3.11)*(i-70.)/10.
      enddo
      do i=81,90
        fsd95(i)=3.11-(3.11-3.10)*(i-80.)/10.
      enddo
      do i=91,100
        fsd95(i)=3.10-(3.10-3.09)*(i-90.)/10.
      enddo
      do i=101,500
        fsd95(i)=3.09-(3.09-3.01)*(i-100.)/400.
      enddo
      do i=501,1000
        fsd95(i)=3.01-(3.01-3.00)*(i-500.)/500.
      enddo
      return
      end