#! /bin/sh
#
#    Last change: Norel Rimbu 2004
#
if [ $# -lt 5 ]
then
  echo " "
  echo "           *****     PARTIAL-CORWINDOW     *****      "
  echo " "
  echo "      The program performs window correlations"
  echo " "
  echo " Syntax: corwindow <file1> <file2> <file3> <file4> <win_length>"
  echo " "
  echo "    <file1>      - input series 1 (one record)"
  echo "    <file2>      - input series 2 (one record)"
  echo "    <file3>      - input series 3 (one record)"
  echo "    <file4>      - output series"
  echo "    <win_length> - window length"
  echo " "
  exit 0
else
 echo " "
 echo "     *****     CORWINDOW     *****      "
 echo " "
 echo "file1="$1
 echo "file2="$2
 echo "file3="$3
 echo "file4="$4
 echo "Window length="$5
 cp $1 fort.1
 cp $2 fort.2
 cp $3 fort.3
fi
cat > corwindow.f <<EOF
C
C      Correlates two time series in moving WINDOW 
C
      program corwindow
      dimension x(20000),y(20000),z(20000)
       dimension rxy(20000),rxz(20000),ryz(20000),rpar(20000)
      integer ia(4),ib(4)

      read(*,'(i10)')nw

C  I read files of an input

      read(1)(ia(l),l=1,4)
      n=ia(4)
      read(1)(x(l),l=1,n)
      read(2)(ib(l),l=1,4)
      n1=ib(4)
      read(2)(y(l),l=1,n1)
         read(3)(ib(l),l=1,4)
      n1=ib(4)
        read(3)(z(l),l=1,n1)
c        write(*,*)(z(l),l=1,n1)
C   Has stopped to read files of an input 

C  The algorithm of calculation of correlation 1 with 2 follows 

             
c   I calculate average 

      nf=n-nw
C      write(*,*)nf
C      read(*,*)

C  The algorithm of calculation of correlation 1 with 2 and 1 with 3 follows 

      do 99 k=1,nf
        s1=0.0
        s2=0.0
        s3=0.0
        nnw=k+nw-1
        do 1 l=k,nnw
          s1=s1+x(l)
          s2=s2+y(l)
          s3=s3+z(l)
    1   continue
        a=nw
        xm=s1/a
        ym=s2/a
        zm=s3/a      
c I calculate STANDARD DEVIATION

        s1=0.0
        s2=0.0
        s3=0.0 
        do 2 l=k,nnw
          s1=s1+(x(l)-xm)*(x(l)-xm)
          s2=s2+(y(l)-ym)*(y(l)-ym)
          s3=s3+(z(l)-zm)*(z(l)-zm)
    2   continue
        xd=sqrt(s1/a)
        yd=sqrt(s2/a)
        zd=sqrt(s3/a)

c I calculate the cross sum 

        s1=0.0
        s2=0.0
        s3=0.0
        do 3 l=k,nnw
          s1=s1+(x(l)-xm)*(y(l)-ym)
          s2=s2+(x(l)-xm)*(z(l)-zm)
          s3=s3+(y(l)-ym)*(z(l)-zm)
    3   continue
        sxy=s1/a
        sxz=s2/a
        syz=s3/a  

c I calculate factor of correlation 

        rxy(k)=sxy/(xd*yd)
        rxz(k)=sxz/(xd*zd)
        ryz(k)=syz/(yd*zd)
        rpar(k)=(rxy(k)-rxz(k)*ryz(k))/
     *sqrt( (1-rxz(k)*rxz(k))*(1-ryz(k)*ryz(k)))
c        write(*,*)ryz(k)
   99 continue

C I write results 

      open(4,form='unformatted')

      ia(4)=nf
      write(4)(ia(l),l=1,4)
      write(4)(rpar(l),l=1,nf)
                         
      close(1)
      close(2)
      close(3)
      close(4)
      print *,'REGULAR END!'
      end
EOF
f77 corwindow.f -o corwindow.x
corwindow.x $1 $2 $3 $4<< M
$5
M
cp fort.4 $4
rm fort*
rm corwindow.x corwindow.f
exit