c c........|.........|.........|.........|.........|.........|.........|.. c........|.........|.........|.........|.........|.........|.........|.. c c Rmatrix and associated programs may be used freely by University c researchers and teachers. Commercial parties interested in using c the program should contact C. J. Thomson. c c........|.........|.........|.........|.........|.........|.........|.. c........|.........|.........|.........|.........|.........|.........|.. c c subroutine aniso6c1(a,rho,px,py,eval,evec) implicit real*8 (a-h,o-z) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c Finds e-vals and e-vecs for given layer. c Based on Sean Guest's routines for finding vertical slowness given c horizontal slownesses by solving the 6X6 matrix problem (using c public domain library routines - Eispack? They look like it.). c dimension a(3,3,3,3),p(3) complex*16 px,py complex*16 eval(6),evec(6,6),xnorm logical iso common/info/lay,id,iso iso=.false. p(1)=dble(px) p(2)=dble(py) call eig (rho,a,p,eval,evec) c Normalise these evecs to unit displacement do 100 j=1,6 xnorm=cdsqrt(evec(1,j)**2+evec(2,j)**2+evec(3,j)**2) do 101 i=1,6 evec(i,j)=evec(i,j)/xnorm 101 continue 100 continue return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine eig (den,ae,p,wc,zc) c Copyright (c) 1996 W. S. Guest and C. J. Thomson. c All rights reserved by the authors. c c W.S. Guest October 1988 c c This subroutine finds the eigenvalues (vertical slowness) and c eigenvectors (displacement-stress) when given the 81 component c elastic constants AE, the density DEN, and the horizontal slownesses c p(1) and p(2). If the medium is isotropic, then a call is made c to ISOEIG which provide a rapid analytic solution (Fryer & Frazer). c THe anisotropic formulation is after Woodhouse (1974) and then the c system matrix is reduced using the `numerical recipes' routine DEIGV. c The vectors must then be sorted into c Pu,SVu,SHu,Pd,SVd,SHd c This is done in SORT which bases the S orientation on polarization c or in SORT_AMP which sorts the eigenvalues from smallest to largest. c implicit real*8(a-h,o-z) c integer i,j,k,l,m,n,lay,id logical iso dimension c(3,3,3,3),ae(3,3,3,3),cm(3,3),ci(3,3),t(3,3),tt(3,3) dimension ct2(3,3,3),ct(3,3),qc(3,3,3,3),s(3,3),a(6,6),p(3) dimension zi(6,6),wr(6),wi(6),zr(6,6) complex*16 wc(6),zc(6,6) c Common Block common/info/lay,id,iso c c c write (*,*) 'ineig',iso do 10 i=1,3 do 20 j=1,3 do 30 k=1,3 do 40 l=1,3 c(i,j,k,l) = ae(i,j,k,l) * den 40 continue 30 continue 20 continue 10 continue c if (iso) then c write (*,*)'IsoEig used to solve "A"' c33=c(3,3,3,3) c44=c(2,3,2,3) c write (*,*) c33,c44 call isoeig(den,p,c33,c44,wc,zc) goto 999 endif c c Calculate the c33 inverse matrix c cm(1,1)=c(2,3,2,3)*c(3,3,3,3) - c(2,3,3,3)**2.d0 cm(2,2)=c(1,3,1,3)*c(3,3,3,3) - c(1,3,3,3)**2.d0 cm(3,3)=c(1,3,1,3)*c(2,3,2,3) - c(1,3,2,3)**2.d0 cm(1,2)=-1.d0*(c(2,3,1,3)*c(3,3,3,3) - c(2,3,3,3)*c(3,3,1,3)) cm(1,3)=c(2,3,1,3)*c(3,3,2,3) - c(2,3,2,3)*c(3,3,1,3) cm(2,3)=-1.d0*(c(1,3,1,3)*c(3,3,2,3) - c(1,3,2,3)*c(3,3,1,3)) cm(2,1)=cm(1,2) cm(3,1)=cm(1,3) cm(3,2)=cm(2,3) det=c(1,3,1,3)*cm(1,1) + c(2,3,1,3)*cm(2,1) + c(3,3,1,3)*cm(3,1) c do 50 i=1,3 do 55 j=1,3 ci(i,j)=cm(i,j)/det 55 continue 50 continue c c Calculate the 'T' and 'T' transpose matrices c do 60 m=1,3 do 65 n=1,3 c1=0.d0 c2=0.d0 do 70 k=1,3 c1 = c1 + c(m,1,k,3)*ci(k,n) c2 = c2 + c(m,2,k,3)*ci(k,n) 70 continue tt(n,m) = -1.d0*(p(1)*c1 + p(2)*c2) t(m,n) = tt(n,m) 65 continue 60 continue c c Calculate the 'Q' matrix c do 75 i=1,2 do 80 k=1,3 do 85 l=1,3 ct(k,l)=0.d0 do 90 m=1,3 ct(k,l) = ct(k,l) + c(k,i,m,3)*ci(m,l) 90 continue 85 continue 80 continue c do 95 j=1,2 do 100 k=1,3 do 105 l=1,3 ct2(j,k,l)=0.d0 do 110 m=1,3 ct2(j,k,l)= ct2(j,k,l) + ct(k,m)*c(m,3,l,j) 110 continue qc(i,j,k,l)=c(k,i,l,j) - ct2(j,k,l) 105 continue 100 continue 95 continue 75 continue c c Calculate the 'S' matrix c do 115 k=1,3 do 120 l=1,3 s(k,l)=0.d0 do 125 i=1,2 do 130 j=1,2 s(k,l) = s(k,l) + p(i)*p(j)*qc(i,j,k,l) 130 continue 125 continue s(k,l) = -1.d0 *s(k,l) if (k.eq.l) s(k,l)=s(k,l)+den 120 continue 115 continue c c Form the complete system matrix 'A' c do 135 i=1,3 j=i+3 do 140 k=1,3 l=k+3 a(i,k) = tt(i,k) a(i,l) = ci(i,k) a(j,k) = s(i,k) a(j,l) = t(i,k) 140 continue 135 continue c c Find the eigenvalues and eigenvectors c call deigv(6,a,wr,wi,zr,zi) c call sort_pol(wr,wi,zr,zi,wc,zc,ae,p) c call sort_amp(wr,wi,zr,zi,wc,zc,ae,p) c 999 continue c c write (19,*) c write (19,'(6e12.3)') wc(1),wc(2),wc(3),wc(4),wc(5),wc(6) c write (19,*) c write (19,'(12d10.2)') zc(1,1),zc(1,2),zc(1,3), c * zc(1,4),zc(1,5),zc(1,6) c write (19,'(12d10.2)') zc(2,1),zc(2,2),zc(2,3), c * zc(2,4),zc(2,5),zc(2,6) c write (19,'(12d10.2)') zc(3,1),zc(3,2),zc(3,3), c * zc(3,4),zc(3,5),zc(3,6) c return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine sort_amp (wr,wi,zr,zi,wc,zc,ae,p) c Copyright (c) 1996 W. S. Guest and C. J. Thomson. c All rights reserved by the authors. c c same as SORT except that the shear eigenvalue separation is based c soley on magnitudes. If all are real then the eigenvalues are arranged c as 1, the smallest magnitude, to 3, the largest. This means that 3 will c be the last to go evanescant. When all three waves are imaginary then the c magnitude of the imaginary components are monitored and they are arranged c from 1, the largest, to 3, the smallest. If there are no intersections then c this will guarentee smooth flowing curves. c implicit real*8 (a-h,o-z) c integer np(6),neg(3),pos(3),neg_ev(3),pos_ev(3),in_ev,ip_ev, + nvmax,pvmax,nvmin,pvmin,in,ip,i,j,k dimension wr(6),wi(6),zr(6,6),zi(6,6),p(3),v(3),ae(3,3,3,3) complex*16 wc(6),zc(6,6) c Common Blocks common/info/lay,id,iso c c numreal=0 in=0 ip=0 in_ev=0 ip_ev=0 velmaxn=0.d0 velmaxp=0.d0 velminn=10000.d0 velminp=10000.d0 holdn=-10.d0 holdp=10.d0 c do i=1,6 np(i)=0 enddo c c write (19,'(6e10.3)') wr c write (19,'(6e10.3)') wi c c c do i=1,6 if (dabs(wi(i)) .lt. 1.d-9)then numreal=numreal+1 p(3)=wr(i) call vtot(ae,p,v) velocity=dsqrt(v(1)**2. + v(2)**2. + v(3)**2.) c write (19,'(i4,4f10.1)')i,v,velocity c if (v(3) .lt. 0.d0)then in=in+1 neg(in)=i if (velocity .lt. velminn)then velminn=velocity nvmin=i endif if (velocity .gt. velmaxn)then velmaxn=velocity nvmax=i endif else ip=ip+1 pos(ip)=i if (velocity .lt. velminp)then velminp=velocity pvmin=i endif if (velocity .gt. velmaxp)then velmaxp=velocity pvmax=i endif endif elseif (wi(i) .lt. 0.d0)then in=in+1 in_ev=in_ev+1 neg_ev(in_ev)=i neg(in)=i else ip=ip+1 ip_ev=ip_ev+1 pos_ev(ip_ev)=i pos(ip)=i endif enddo c c write (19,*) ip,in,pos,neg c write (19,*) 'number of real eigenvalues',numreal c write (19,*) pos_ev,neg_ev c write (19,*) 'pvmax,pvmin',pvmax,pvmin c write (19,*) 'nvmax,nvmin',nvmax,nvmin c j=0 k=0 c if (numreal .ge. 4)then c np(3)=pvmin np(6)=nvmin c if (numreal .eq. 6)then np(1)=pvmax np(4)=nvmax else np(1)=pos_ev(1) np(4)=neg_ev(1) endif c do i=1,3 if (pos(i).ne.np(1) .and. pos(i).ne.np(3)) np(2)=pos(i) if (neg(i).ne.np(4) .and. neg(i).ne.np(6)) np(5)=neg(i) enddo c elseif (numreal .ge. 2)then c np(3)=pvmin np(6)=nvmin c if (wi(pos_ev(1)) .gt. wi(pos_ev(2)))then np(1)=pos_ev(1) np(2)=pos_ev(2) else np(2)=pos_ev(1) np(1)=pos_ev(2) endif if (wi(neg_ev(1)) .lt. wi(neg_ev(2)))then np(4)=neg_ev(1) np(5)=neg_ev(2) else np(5)=neg_ev(1) np(4)=neg_ev(2) endif c else c do 200 i=1,3 if (wi(pos_ev(1)) .gt. wi(pos_ev(2))) then j=j+1 if (j .eq. 1)then holdp=wi(pos(i)) np(2)=pos(i) elseif (wi(pos(i)) .lt. holdp) then np(3)=np(2) np(2)=pos(i) else np(3)= pos(i) endif endif c if (neg(i) .ne. np(4)) then k=k+1 if (j .eq. 1)then holdn=wi(neg(i)) np(5)=neg(i) elseif (wi(neg(i)) .gt. holdn) then np(6)=np(5) np(5)=neg(i) else np(6)= neg(i) endif endif 200 continue c endif c c write (19,*) 'np',np c c Switch the eigenvalues and eigenvectors to the order P,S1,S2 c do i=1,6 wc(i)=dcmplx(wr(np(i)),wi(np(i))) do j=1,3 zc(j,i)=dcmplx(zr(j,np(i)),zi(j,np(i))) c zc(j+3,i)=dcmplx(zr(j+3,np(i)),zi(j+3,np(i)))*(0.d0,1.d0) zc(j+3,i)=dcmplx(zr(j+3,np(i)),zi(j+3,np(i))) enddo enddo c return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine sort_pol (wr,wi,zr,zi,wc,zc,ae,p) c Copyright (c) 1996 W. S. Guest and C. J. Thomson. c All rights reserved by the authors. c c W.S. Guest July 1989 c c This subroutine takes the input of the eigenvalues and eigenvectors c calculated in eig and orders them into pu,svu,shu,pd,svd,shd. a c check is first made to see if any of the eigenvalues are complex c representing evanescant waves. separation of the up and downgoing c waves is based on th sign of the eigenvalue but a check is made to c make sure there are three of each sign, if not then the total vel. c are calculated to get energy flux direction. c if all of the eigenvalues are real, c then the p-wave is taken as the smallest eigenvalues, both in the c up and downgoing direction. the s-waves are not sorted, but c always come out as sv and then sh. s-waves can be separated based c on displacement, but problem occur for vertically travelling waves c implicit real*8 (a-h,o-z) c integer i,j,k,np(6),neg(3),npos(3) dimension wr(6),wi(6),zr(6,6),zi(6,6), + p(3),v(3),ae(3,3,3,3) complex*16 wc(6),zc(6,6) c Common Blocks common/info/lay,id,iso c Parameter defintions parameter (pi = 3.1415926535897932d0) c j=0 k=0 holdn=-10.d0 holdp=10.d0 numev=0 hldevn=0.d0 hldevp=0.d0 do 5 i=1,6 np(i)=0 5 continue c c Check to see is any of the waves are evanesant c do 10 i=1,6 if (wi(i) .ne. 0.d0)then if (wi(i) .lt. 0.d0)then j=j+1 neg(j)=i numev=numev+1 if (wi(i) .lt. hldevn) then np(4)=i hldevn=wi(i) endif else numev=numev+1 k=k+1 npos(k)=i if (wi(i) .gt. hldevp) then np(1)=i hldevp=wi(i) endif endif else if (wr(i) .lt. 0.d0)then j=j+1 neg(j)=i if (numev .eq. 0) then if (wr(i) .gt. holdn)then holdn=wr(i) np(4)=i endif endif else k=k+1 npos(k)=i if (numev .eq. 0) then if (wr(i) .lt. holdp)then holdp=wr(i) np(1)=i endif endif endif endif 10 continue c c write (19,*) 'numev',numev,' j',j,' k',k c if(j .gt. 3 .or. k .gt. 3)then j=0 k=0 do 20 i=1,6 if (dabs(wi(i)) .lt. 1.d-9)then p(3)=-wr(i) call vtot(ae,p,v) if (v(3) .gt. 0.d0)then j=j+1 neg(j)=i else k=k+1 npos(k)=i endif elseif (wi(i) .lt. 0.d0)then j=j+1 neg(j)=i else k=k+1 npos(k)=i endif 20 continue endif c j=0 k=0 if (numev .le. 2)then do 755 i=1,3 if (npos(i) .ne. np(1)) then j=j+1 hfrac=dsqrt(dble(zr(1,npos(i))**2.d0)+dble * (zr(2,npos(i)))**2.d0)/dsqrt(dble(zr(1,npos(i)))**2.d0 * +dble(zr(2,npos(i)))**2.d0+dble(zr(3,npos(i)))**2.d0) if (j .eq. 1) then hfrac1=hfrac ihld=i elseif (dabs(hfrac1) .gt. dabs(hfrac))then np(3)= npos(ihld) np(2)= npos(i) else np(3)= npos(i) np(2)= npos(ihld) endif endif if (neg(i) .ne. np(4)) then k=k+1 hfracn=dsqrt(dble(zr(1,neg(i)))**2.d0+ * dble(zr(2,neg(i)))**2.d0)/dsqrt(dble(zr(1,neg(i)))**2.d0 * +dble(zr(2,neg(i)))**2.d0+dble(zr(3,neg(i)))**2.d0) if (k .eq. 1) then hfracn1=hfracn ihld=i elseif (dabs(hfracn1) .gt. dabs(hfracn))then np(6)= neg(ihld) np(5)= neg(i) else np(6)= neg(i) np(5)= neg(ihld) endif endif 755 continue c elseif (numev .le. 4)then c do 756 i=1,3 c if (npos(i) .ne. np(1)) then c j=j+1 c if (j .eq. 1) then c if (wi(npos(i)) .gt. 0.d0)then c ihld=i c else c if (dabs(dble(zr(2,npos(i)))) .gt. c * dabs(dble(zr(3,npos(i)))))then c np(3)=npos(i) c else c np(2)=npos(i) c endif c endif c else c if (wi(npos(i)) .gt. 0.d0)then c if (np(3) .eq. 0)then c np(3)=npos(i) c else c np(2)=npos(i) c endif c else c if (dabs(dble(zr(2,npos(i)))) .gt. c * dabs(dble(zr(3,npos(i)))))then c np(2)=npos(ihld) c np(3)=npos(i) c else c c np(2)=npos(i) c np(3)=npos(ihld) c endif c endif c endif c endif if (npos(i) .ne. np(1)) then if (wi(npos(i)) .gt. 0.d0)then np(5-id)=npos(i) c np(2)=neg(i) else np(id)=npos(i) c np(3)=neg(i) endif endif c if (neg(i) .ne. np(4)) then c k=k+1 c if (k .eq. 1) then c if (wi(neg(i)) .gt. 0.d0)then c jhld=i c else c if (dabs(dble(zr(2,neg(i)))) .gt. c * dabs(dble(zr(3,neg(i)))))then c np(6)=neg(i) c else c c np(5)=neg(i) c endif c endif c else c if (wi(neg(i)) .gt. 0.d0)then c if (np(6) .eq. 0)then c np(6)=neg(i) c else c np(5)=neg(i) c endif c else c if (dabs(dble(zr(2,neg(i)))) .gt. c * dabs(dble(zr(3,neg(i)))))then c np(5)=neg(jhld) c np(6)=neg(i) c else c np(5)=neg(i) c np(6)=neg(jhld) c endif c endif c endif c endif c if (neg(i) .ne. np(4)) then if (wi(neg(i)) .lt. 0.d0)then np(8-id)=neg(i) c np(5)=neg(i) else c np(6)=neg(i) np(3+id)=neg(i) endif endif 756 continue c elseif (numev .le. 6)then c do 757 i=1,3 if (npos(i) .ne. np(1)) then if (np(2) .eq. 0.d0)then np(2)=npos(i) else np(3)=npos(i) endif endif if (neg(i) .ne. np(4)) then if (np(5) .eq. 0.d0)then np(5)=neg(i) else np(6)=neg(i) endif endif 757 continue endif c c write (19,*) 'np',np c c Switch the eigenvalues and eigenvectors to the order P,S1,S2 c do i=1,6 wc(i)=dcmplx(wr(np(i)),wi(np(i))) do j=1,3 zc(j,i)=dcmplx(zr(j,np(i)),zi(j,np(i))) c zc(j+3,i)=dcmplx(zr(j+3,np(i)),zi(j+3,np(i)))*(0.d0,1.d0) zc(j+3,i)=dcmplx(zr(j+3,np(i)),zi(j+3,np(i))) enddo enddo c return end c c........|.........|.........|.........|.........|.........|.........|.. c SUBROUTINE VTOT(CC,P,V) C IMPLICIT REAL*8 (A-H,O-Z) integer lay,id logical iso DIMENSION CC(3,3,3,3),A(3,3),CP(3,3),DCP(3,3,3),V(3),P(3),DET(3) COMMON/info/lay,id,iso C if (iso)then c pmag=p(1)**2.d0+p(2)**2.d0+p(3)**2.d0 do 10 i=1,3 v(i)=p(i)/pmag 10 continue c else c DO 28 I=1,3 DO 27 L=1,3 PSUM = 0.0D0 DO 26 J=1,3 DO 25 K=1,3 PSUM = PSUM+CC(I,J,K,L)*P(J)*P(K) 25 CONTINUE 26 CONTINUE CP(I,L) = PSUM CP(L,I) = PSUM 27 CONTINUE 28 CONTINUE C C DO 34 M=1,3 do 33 i=1,3 DO 32 L=1,3 X = 0.0D0 Y = 0.0D0 DO 30 J=1,3 X = X+CC(I,J,M,L)*P(J) 30 CONTINUE DO 31 K=1,3 Y = Y+CC(I,M,K,L)*P(K) 31 CONTINUE DCP(M,I,L) = X+Y 32 CONTINUE 33 CONTINUE 34 CONTINUE C C DO 45 I=1,3 DET(I)=0.0D0 45 CONTINUE C C DO 55 N=1,3 DO 54 M=1,3 DO 51 I=1,3 DO 50 J=1,3 A(I,J)=CP(I,J) 50 CONTINUE 51 CONTINUE DO 52 I=1,3 A(I,I) = CP(I,I)-1.D0 52 CONTINUE DO 53 I=1,3 A(I,M) = DCP(N,I,M) 53 CONTINUE Z1 = A(1,1)*A(2,2)*A(3,3)-A(1,1)*A(2,3)*A(3,2) Z2 = -A(1,2)*A(2,1)*A(3,3)+A(1,2)*A(2,3)*A(3,1) Z3 = A(1,3)*A(2,1)*A(3,2)-A(1,3)*A(2,2)*A(3,1) DET(N) = DET(N)+Z1+Z2+Z3 54 CONTINUE 55 CONTINUE C DEN = 0.0D0 DO 65 I=1,3 DEN = DEN + P(I)*DET(I) 65 CONTINUE C DO 71 I=1,3 V(I) = DET(I)/DEN 71 CONTINUE c endif C RETURN END c c........|.........|.........|.........|.........|.........|.........|.. c subroutine isoeig(den,p,c33,c44,wc,zc) c Copyright (c) 1996 W. S. Guest and C. J. Thomson. c All rights reserved by the authors. c c Original version: W.S. Guest November 1988 c c Modified/Corrected: C. Thomson October 1993 c Sean's original lines commented out for c comparison. c NOTE: c33 = lambda + 2mu and c44 = mu in Lame notation c c c Sean's comments follow... c c ``This subroutine calculates the displacement-stress eigenvectors c and the eigenvalues in the case of an isotropic medium. The coding c follows Fryer and Frazer (1987) with the obvious errors corrected. c The vectors are choosen to be in a right handed system. P is along c the slowness vector, SV is pointed upwards, and SH is horizontal c but its direction depends on the value of p(x) (It goes to the c left of the p(x) direction). All of the vectors are normalized so c that the displacements are of unit magnitude. c c Orientations rechecked 07/92 Sean Guest'' c implicit real*8(a-h,o-z) c integer i dimension p(3) complex*16 wc(6),zc(6,6),norm,acomp,ccomp,bcomp c p2=p(1)**2.d0 +p(2)**2.d0 fac1=den/c33-p2 fac2=den/c44-p2 c write (19,*) p2,fac1,fac2 c c Define the eigenvalues c wc(1)=cdsqrt(dcmplx(fac1,0.d0)) wc(2)=cdsqrt(dcmplx(fac2,0.d0)) wc(3)=wc(2) wc(4)=-wc(1) wc(5)=-wc(2) wc(6)=-wc(3) c c write (19,'(12e10.3)') wc c c Evaluation of the Eigenvectors c Solution of the Eigenvectors c do 10 i=0,3,3 c zc(1,1+i)=p(1) zc(2,1+i)=p(2) zc(3,1+i)=wc(1+i) zc(4,1+i)=2.d0*p(1)*c44*wc(1+i) zc(5,1+i)=2.d0*p(2)*c44*wc(1+i) zc(6,1+i)=den-2.d0*p2*c44 c c zc(4,1+i)=-2.d0*p(1)*c44*wc(1+i) c zc(5,1+i)=-2.d0*p(2)*c44*wc(1+i) c zc(6,1+i)=2.d0*p2*c44-den c zc(1,2+i)=p(1) zc(2,2+i)=p(2) zc(3,2+i)=-p2/wc(2+i) zc(4,2+i)=p(1)*zc(6,1+i)/wc(2+i) zc(5,2+i)=p(2)*zc(6,1+i)/wc(2+i) zc(6,2+i)=-2.d0*p2*c44 c c zc(4,2+i)=p(1)*(2.d0*p2*c44-den)/wc(2+i) c zc(5,2+i)=p(2)*(2.d0*p2*c44-den)/wc(2+i) c zc(6,2+i)=2.d0*p2*c44 c if (dble(wc(2+i)) .gt. 0.d0)then c do j=1,6 c zc(j,2+i) = -zc(j,2+i) c enddo c endif c zc(1,3+i)=-p(2) zc(2,3+i)=p(1) zc(3,3+i)=0.d0 zc(4,3+i)=-p(2)*c44*wc(3+i) zc(5,3+i)=p(1)*c44*wc(3+i) zc(6,3+i)=0.d0 c 10 continue c cs do 40 j=1,6 cs do 30 i=4,6 cs zc(i,j)=zc(i,j)*(0.d0,-1.d0) cs 30 continue cs 40 continue c c Normalize the vectors based on displacement c do 50 j=1,6 acomp=zc(1,j)*zc(1,j) bcomp=zc(2,j)*zc(2,j) ccomp=zc(3,j)*zc(3,j) norm=cdsqrt(acomp+bcomp+ccomp) do 60 i=1,6 zc(i,j)=zc(i,j)/norm 60 continue 50 continue c return end c c........|.........|.........|.........|.........|.........|.........|.. c