$debug c program for cross-spectrum ananysis program main parameter(nt=852,mm=160) dimension x1(nt),x2(nt) dimension rx(0:mm) dimension ry(0:mm) dimension rxy(0:mm) dimension ryx(0:mm) dimension sx1(0:mm) dimension sx2(0:mm) dimension sy1(0:mm) dimension sy2(0:mm) dimension p1(0:mm) dimension q1(0:mm) dimension p2(0:mm) dimension q2(0:mm) dimension r12(0:mm) dimension f12(0:mm) dimension zta(0:mm) c dimension ar(nt,num),ai(nt,num) c========= read POP coefficients to Num.6 device =================== c open(16,file='d:\for\pop2.time',status='old') c read(16,*) c do 79 k=1,num c read(16,*) n c79 read(16,*) (ar(i,k),i=1,nt) c write(*,*) n c read(16,*) c do 89 k=1,num c read(16,*) n c89 read(16,*) (ai(i,k),i=1,nt) c write(*,*) n c close(16) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c do 10 it=1,nt c x1(it)=ar(it,5) c10 x2(it)=ai(it,5) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c open(1,file='31r') c read(1,*) x1 c open(2,file='31i') c read(2,*) x2 c close(1) c close(2) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do 10 i=1,nt 10 x1(i)=4.0*cos((2*3.14159/100.)*i)+3.0*cos((6.28/50)*i) do 110 i=1,nt 110 x2(i)=2.0*sin((2*3.14159/50.)*i) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc call cross(x1,x2,nt,mm,rx,ry,rxy,ryx,sx1,sx2,sy1,sy2,p1,q1, $ p2,q2,r12,f12,zta) stop end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine cross(x,y,n,m,rx,ry,rxy,ryx,sx1,sx2,sy1,sy2, $ p1,q1,p2,q2,r12,f12,zta) dimension x(n),y(n) dimension rx(0:m) dimension ry(0:m) dimension rxy(0:m) dimension ryx(0:m) dimension sx1(0:m) dimension sx2(0:m) dimension sy1(0:m) dimension sy2(0:m) dimension p1(0:m) dimension q1(0:m) dimension p2(0:m) dimension q2(0:m) dimension r12(0:m) dimension f12(0:m) dimension zta(0:m) open(6,file='c:\lunwen\sp.dat') c write (6,605) c write (6,606) x c write (6,607) c write (6,606) y xm=0. ym=0. do 10 i=1,n xm=xm+x(i)/real(n) 10 ym=ym+y(i)/real(n) sx=0. sy=0. do 20 i=1,n sx=sx+(x(i)-xm)**2 20 sy=sy+(y(i)-ym)**2 sx=sqrt(sx/real(n)) sy=sqrt(sy/real(n)) do 30 i=1,n x(i)=(x(i)-xm)/sx 30 y(i)=(y(i)-ym)/sy do 40 itau=0,m rx(itau)=0. ry(itau)=0. rxy(itau)=0. ryx(itau)=0. do 45 it=1,n-itau rx(itau)=rx(itau)+x(it)*x(it+itau) ry(itau)=ry(itau)+y(it)*y(it+itau) rxy(itau)=rxy(itau)+x(it)*y(it+itau) ryx(itau)=ryx(itau)+y(it)*x(it+itau) 45 continue rx(itau)=rx(itau)/real(n-itau) ry(itau)=ry(itau)/real(n-itau) rxy(itau)=rxy(itau)/real(n-itau) ryx(itau)=ryx(itau)/real(n-itau) 40 continue c write(6,111) (rx(i),i=0,m) 111 format(3x,'rx(tau) tau=1--mm',/2x,10f8.2) c write(6,112) (ry(i),i=0,m) 112 format(3x,'ry(tau) tau=1--mm',/2x,10f6.2) do 50 j=0,m sx1(j)=rx(0)+(-1.0)**j*rx(m) sy1(j)=ry(0)+(-1.0)**j*ry(m) do 55 k=1,m-1 sx1(j)=sx1(j)+2.0*rx(k)*cos(real(j*k)*3.14159/real(m)) sy1(j)=sy1(j)+2.0*ry(k)*cos(real(j*k)*3.14159/real(m)) 55 continue sx1(j)=2.0*sx1(j) sy1(j)=2.0*sy1(j) 50 continue do 60 j=0,m p1(j)=(rxy(0)+ryx(0))*0.5+(-1.0)**j*(rxy(m)+ryx(m))*0.5 do 65 k=1,m-1 p1(j)=p1(j)+(rxy(k)+ryx(k))*cos(real(j*k)*3.14159/real(m)) 65 continue p1(j)=p1(j)*2.0 60 continue do 70 j=0,m q1(j)=0. do 75 k=1,m-1 75 q1(j)=q1(j)+2.0*(rxy(k)-ryx(k))*sin(real(k*j)*3.14159/real(m)) 70 continue ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc sx2(0)=(sx1(0)+sx1(1))/2.0 sy2(0)=(sy1(0)+sy1(1))/2.0 p2(0)=(p1(0)+p1(1))/2.0 q2(0)=(q1(0)+q1(1))/2.0 sx2(m)=(sx1(m)+sx1(m-1))/2.0 sy2(m)=(sy1(m)+sy1(m-1))/2.0 p2(m)=(p1(m)+p1(m-1))/2.0 q2(m)=(q1(m)+q1(m-1))/2.0 do 73 j=1,m-1 sx2(j)=0.25*sx1(j-1)+0.5*sx1(j)+0.25*sx1(j+1) sy2(j)=0.25*sy1(j-1)+0.5*sy1(j)+0.25*sy1(j+1) p2(j)=0.25*p1(j-1)+0.5*p1(j)+0.25*p1(j+1) q2(j)=0.25*q1(j-1)+0.5*q1(j)+0.25*q1(j+1) 73 continue cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc spx=0. spy=0. do 415 j=1,m spx=spx+sx2(j)/real(m) 415 spy=spy+sy2(j)/real(m) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc xalpha=6.01 ccccccccccccccccccccccccccccccccccccccccccccccccccccc h=(2.0*n-m/2.0)/m write(*,*) 'H= ',h do 417 j=1,m sx1(j)=xalpha*spx*(1.0-rx(1)**2)/ $ (h*(1.0+rx(1)**2-2.0*rx(1)*cos(3.14159*real(j)/real(m)))) sy1(j)=xalpha*spy*(1.0-ry(1)**2)/ $ (h*(1.0+ry(1)**2-2.0*ry(1)*cos(3.14159*real(j)/real(m)))) 417 continue rnu=(2.0*n-0.5*m)/real(m) do 80 j=1,m r12(j)=(p2(j)**2+q2(j)**2)/(sx2(j)*sy2(j)) c if(r12(j).gt.1.0) print *,'error ','j=',j,',r12(j)= ',r12(j) 80 continue do 90 j=1,m zta(j)=atan2(q2(j),p2(j)) 90 continue do 100 j=1,m 100 f12(j)=(rnu-1.0)*r12(j)/(1.0-r12(j)) 605 format(/2x,'the original series x') 607 format(/2x,'the original series y') 606 format((1x,10f7.2)) write(6,610) 610 format(/10x,'cross spectral ananysis of x and y ') write(6,620) (i,i=1,m) 620 format(/5x,' order ',10i6) 630 format(1x,' freqency ',10f6.2) write (6,640) (2.0*m/real(j),j=1,m) 640 format(4x,' period ',20f6.2) write(6,650) (sx2(j),j=1,m) 650 format(1x,' spe. of x',20f6.2) write(6,652) (sx1(j),j=1,m) 652 format(1x,'x-sp-test',20f6.2) write (6,660) (sy2(j),j=1,m) 660 format(1x,'spe. of y' ,20f6.2) write(6,662) (sy1(j),j=1,m) 662 format(1x,'y-sp-test',20f6.2) write (6,670) (p2(j),j=1,m) 670 format(1x,' real part ',20f6.2) write(6,675) (q2(j),j=1,m) 675 format(1x,'imaginary',20f6.2) write (6,680) (r12(j),j=1,m) 680 format(1x,'coherence',20f6.2) write(6,685) (zta(j)*57.3,j=1,m) 685 format(1x,'phase lag',20f6.2) write (6,690) (real(m)*zta(j)/(3.14159*j),j=1,m) 690 format(2x,' time lag ',20f6.2) write (6,695) (f12(j),j=1,m) 695 format(3x,'f-value',20f6.2) inu=int(rnu) write(*,699) inu 699 format(/5x,'degrees of freedom of f-distribution (2,',i3')') ccccccccccccccccccccccccccccccccccccccccccccccccccccccc g1=4.10/(4.10+rnu-1) g2=7.56/(7.56+rnu-1) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc write(*,*) 'g1 ',g1 write(*,*) 'g2= ',g2 open(3,file='phase.dat') do 42 j=0,m 42 write (3,*) real(j),zta(j)*57.3 close(3) open(23,file='time.dat') do 44 j=1,m 44 write (23,*) real(j),real(m)*zta(j)/(3.14159*j) close(23) open(4,file='cohere.dat') do 41 j=0,m 41 write (4,*) real(j),r12(j) close(4) open(9,file='ak.dat') do 991 i=-180,180 write(9,*) real(m)*2.0/77.,real(i) c write(9,*) 9.9,real(i) 991 continue close(9) open(10,file='g1.dat') do 992 i=0,m write(10,'(2f20.15)') real(i),g1 992 continue close(10) open(11,file='g2.dat') do 993 i=0,m write(11,'(2f20.15)') real(i),g2 993 continue close(11) return end