c     This program calculates the Green's function of Ting for an
c     anisotropic elliptical inclusion in an anistropic solid.
c     The Stroh roots and eigenvectors are read in as input.
c
c     Written by J. R. Berger, Colorado School of Mines
c     e-mail: jberger@mines.edu
c     Program developed:  August, 1998
c     Last revision:		11 August 1998	
     
      character*40 file1,file2
c
	complex p1(3),p2(3),sa11(3),sb11(3),sa12(3),sb12(3),
     $sa21(3),sb21(3),sa22(3),sb22(3),A1(3,3),B1(3,3),A2(3,3),
     $B2(3,3),A1inv(3,3),B1inv(3,3),A2inv(3,3),B2inv(3,3),M1(3,3)
     $,M2(3,3),ziota,q1inf(3),A1t(3,3),f(3),zi1(3,3),zi2(3,3),
     $zi3(3,3),ztemp1(3,3),ztemp2(3,3),ztemp3(3,3),zkeep1(3,3),
     $zkeep2(3,3),zvec(3),ztemp(3,3),sa13(3),sb13(3),sa23(3),sb23(3),
     $M1inv(3,3),M2inv(3,3),B1t(3,3),B2t(3,3)
c
      complex q11(3),q12(3),q13(3),q21(3),q22(3),q23(3),b2cinv(3,3),
     $c1(3),c2(3),d1(3),d2(3),t1(3),t2(3),q2(3),z1(3),xi1(3),zterm1,
     $xi1inv(3),z10(3),xi10(3),xi10inv(3),diag(3,3),diag1(3,3),
     $diag2(3,3),diag3(3,3),zvec1(3),zvec2(3),zvec3(3)
     
      complex zL1(3,3),zL2(3,3),zX(3,3),zXinv(3,3),ztau21(3),ztau22(3),
     $ztau23(3),yk(3),bigyk(3,3),czXinv(3,3),cbigyk(3,3),gk(3) 
	
	real a,b,uzip(3),u1(3)
     
c     write(6,*)'Input data filename?'
c	read(5,'(A)')file1
	open(unit=15,file='d:\programs\CuSi.txt',status='old')
	write(6,*)'Output data filename?'
	read(5,'(A)')file2
	open(unit=20,file=file2,status='unknown')
c
c
      ziota=cmplx(0.0,1.0)
	pi=acos(-1.)
c
c--------------------------------------------------------------
c     Read and echo the Stroh roots and eigenvectors for the solid
	read(15,*)p1(1),p1(2),p1(3)
c     In these variables, the first digit is the material, second
c     digit is the associated root.
	read(15,*)sa11(1),sa11(2),sa11(3)
	read(15,*)sa12(1),sa12(2),sa12(3)
	read(15,*)sa13(1),sa13(2),sa13(3)
c
	read(15,*)sb11(1),sb11(2),sb11(3)
	read(15,*)sb12(1),sb12(2),sb12(3)
	read(15,*)sb13(1),sb13(2),sb13(3)
c
c     Read and echo the Stroh roots and eigenvectors for the inclusion
      read(15,*)p2(1),p2(2),p2(3)
c
	read(15,*)sa21(1),sa21(2),sa21(3)
	read(15,*)sa22(1),sa22(2),sa22(3)
	read(15,*)sa23(1),sa23(2),sa23(3)
c
	read(15,*)sb21(1),sb21(2),sb21(3)
	read(15,*)sb22(1),sb22(2),sb22(3)
	read(15,*)sb23(1),sb23(2),sb23(3)
c
c     Read in the major and minor axis dimensions of the ellipse
      write(6,*) 'Major axis dimension for ellipse?'
	read(5,*) a
	write(6,*) 'Minor axis dimension for ellipse?'
	read(5,*) b
c
c     Read in the number of terms to take in the K series
      write(6,*) 'Number of terms to take in series?'
	read(5,*) kterms

c     Read in the location of the source point
      write(6,*) 'Source point location?'
	read(5,*)x10,x20
c

c     Read in the starting and stopping points for plotting
      write(6,*)'x1 value for plotting (x1 > a)?'
	read(5,*)x1
	write(6,*)'Starting x2 value?'
	read(5,*)x2beg
	write(6,*)'Ending x2 value?'
	read(5,*)x2end
	write(6,*)'Increment in x2 for each calculation?'
	read(5,*)xdiv
c
c-----------------------------------------------------------------------
c
c     Form the A and B matrices for each material
      call zeromat(A1)
	call zeromat(B1)
	call zeromat(A2)
	call zeromat(B2)
      do 20 i=1,3
	  A1(i,1)=sa11(i)
	  A1(i,2)=sa12(i)
	  A1(i,3)=sa13(i)
	  B1(i,1)=sb11(i)
	  B1(i,2)=sb12(i)
	  B1(i,3)=sb13(i)
c
	  A2(i,1)=sa21(i)
	  A2(i,2)=sa22(i)
	  A2(i,3)=sa23(i)
	  B2(i,1)=sb21(i)
	  B2(i,2)=sb22(i)
	  B2(i,3)=sb23(i)
20    continue	  
c
c     Check to make sure the orthogonality and closure relations are
c     satisfied
      call orthog(A1,B1)
      call orthog(A2,B2)
c
c     Form the impdeance tensors
      call ctrans(B1,B1t)
	call ctrans(B2,B2t)
      call cinv(A1,A1inv)
	call cinv(B1,B1inv) 
      call cinv(A2,A2inv)
	call cinv(B2,B2inv)
	call matmult(B1,A1inv,M1)
	call matmult(B2,A2inv,M2)
	call matmult(A1,B1inv,M1inv)
	call matmult(A2,B2inv,M2inv)
	do 40 i=1,3
	do 40 j=1,3
	  M1(i,j)=-ziota*M1(i,j)
	  M2(i,j)=-ziota*M2(i,j)
	  M1inv(i,j)=ziota*M1inv(i,j)
	  M2inv(i,j)=ziota*M2inv(i,j)
40    continue
c
c--------------------------------------------------------------------
c
c     Set up the constant vector q1-infininty for a point force in the
c     x1 direction
      call ctrans(A1,A1t)
      f(1)=cmplx(1.0,0.0)
	f(2)=cmplx(0.0,0.0)
	f(3)=cmplx(0.0,0.0)
	q1inf(1)=cmplx(0.0,0.0)
	q1inf(2)=cmplx(0.0,0.0)
	q1inf(3)=cmplx(0.0,0.0)
	do 60 i=1,3
	do 60 j=1,3
	  q1inf(i)=q1inf(i)+A1t(i,j)*f(j)
60    continue
c
c----------------------------------------------------------------------
c     Calculate the 3 q1-beta and 3 q2-beta vectors.
c     First form the I-beta matrices
      call zibeta(zi1,zi2,zi3)
	call zeromat(ztemp1)
	call zeromat(ztemp2)
	call zeromat(ztemp3)
c
c     Now form the matrix terms needed from M1 and M2 for q1-beta's     
	do 80 i=1,3
	do 80 j=1,3
	  ztemp1(i,j)=M1inv(i,j)+conjg(M2inv(i,j))
	  ztemp2(i,j)=conjg(M2inv(i,j))-conjg(M1inv(i,j))
80    continue
      call cinv(ztemp1,ztemp3)
	call zeromat(ztemp1)
	call matmult(ztemp3,ztemp2,ztemp1)
	call zeromat(ztemp2)
	call zeromat(ztemp3)
	do 100 i=1,2
	do 100 j=1,2
	  ztemp2(i,j)=conjg(B1(i,j))
100   continue
      call matmult(ztemp1,ztemp2,zkeep1)
c
c     Calculate the three q1-beta vectors
c     Calculate q11
	call zerovec(q11)
      call zerovec(zvec)
	call matmult(zkeep1,zi1,ztemp1)
	do 120 i=1,3
	do 119 j=1,3
	  zvec(i)=zvec(i)+ztemp1(i,j)*q1inf(j)
119   continue
120   continue
      do 130 i=1,3
	do 129 j=1,3
	  q11(i)=q11(i)+B1inv(i,j)*zvec(j)
129   continue
130   continue
c    
c     Calculate q12
	call zerovec(q12)
      call zerovec(zvec)
	call matmult(zkeep1,zi2,ztemp1)
	do 140 i=1,3
	do 139 j=1,3
	  zvec(i)=zvec(i)+ztemp1(i,j)*q1inf(j)
139   continue
140   continue
	do 150 i=1,3
	do 149 j=1,3
	  q12(i)=q12(i)+B1inv(i,j)*zvec(j)
149   continue
150   continue
c    
c     Calculate q13
      call zerovec(q13)
	call zerovec(zvec)
	call matmult(zkeep1,zi3,ztemp1)
	do 160 i=1,3
	do 159 j=1,3
	  zvec(i)=zvec(i)+ztemp1(i,j)*q1inf(j)
159   continue
160   continue
	do 170 i=1,3
	do 169 j=1,3
	  q13(i)=q13(i)+B1inv(i,j)*zvec(j)
169   continue
170   continue
c
c
c
c     Now form the matrix terms needed from M1 and M2 for q2-beta's
      call zeromat(ztemp1)
	call zeromat(ztemp2)
	call zeromat(ztemp3)
c    
      do 171 i=1,3
	do 171 j=1,3
	   ztemp1(i,j)=conjg(B2(i,j))
171   continue
      call cinv(ztemp1,B2cinv)
c
	do 180 i=1,3
	do 180 j=1,3
	  ztemp1(i,j)=M1inv(i,j)+conjg(M2inv(i,j))
	  ztemp2(i,j)=conjg(M2inv(i,j))+conjg(M1inv(i,j))
180   continue
      call cinv(ztemp1,ztemp3)
	call zeromat(ztemp1)
	call matmult(ztemp3,ztemp2,ztemp1)
	call zeromat(ztemp2)
	call zeromat(ztemp3)
	do 200 i=1,2
	do 200 j=1,2
	  ztemp2(i,j)=conjg(B1(i,j))
200   continue
      call matmult(ztemp1,ztemp2,zkeep2)
c
c     Calculate q21
      call zerovec(q21)
	call zerovec(zvec)
	call matmult(zkeep2,zi1,ztemp1)
	do 220 i=1,3
	do 219 j=1,3
	  zvec(i)=zvec(i)+ztemp1(i,j)*q1inf(j)
219   continue
220   continue
      do 230 i=1,3
	do 229 j=1,3
	  q21(i)=q21(i)+B2cinv(i,j)*zvec(j)
229   continue
230   continue
c 
c     Calculate q22
      call zerovec(q22)
	call zerovec(zvec)
	call matmult(zkeep2,zi2,ztemp1)
	do 240 i=1,3
	do 239 j=1,3
	  zvec(i)=zvec(i)+ztemp1(i,j)*q1inf(j)
239   continue
240   continue
      do 250 i=1,3
	do 249 j=1,3
	  q22(i)=q22(i)+B2cinv(i,j)*zvec(j)
249   continue
250   continue
c
c     Calculate q23
      call zerovec(q23)
	call zerovec(zvec)
	call matmult(zkeep2,zi3,ztemp1)
	do 260 i=1,3
	do 259 j=1,3
	  zvec(i)=zvec(i)+ztemp1(i,j)*q1inf(j)
259   continue
260   continue
      do 270 i=1,3
	do 269 j=1,3
	  q23(i)=q23(i)+B2cinv(i,j)*zvec(j)
269   continue
270   continue
c--------------------------------------------------------------------
c     Calculate the Hirth-Lothe L tensor, and the auxillary matrices
c     needed to get the gk vector.
      call matmult(B1,B1t,ztemp1)
      do 600 i=1,3
	do 600 j=1,3
	  zL1(i,j)=-2.*ziota*ztemp(i,j)
600   continue
      call matmult(B2,B2t,ztemp1)
	do 620 i=1,3
	do 620 j=1,3
	 zL2(i,j)=-2.*ziota*ztemp1(i,j)
620   continue
      do 640 i=1,3
	do 640 j=1,3
	  ztemp1(i,j)=M1inv(i,j)+conjg(M2inv(i,j))
640   continue
      call matmult(zL2,ztemp1,zX)
	do 660 i=1,3
	do 660 j=1,3
	  zX(i,j)=0.5*zX(i,j)
660   continue
      call cinv(zX,zXinv)
	do 661 i=1,3
	do 661 j=1,3
	  czXinv(i,j)=conjg(zXinv(i,j))
661   continue
c----------------------------------------------------------------------
c     Calculate the mapping constants
	do 280 i=1,3
	  c1(i)=0.5*(a-ziota*p1(i)*b)
	  c2(i)=0.5*(a-ziota*p2(i)*b)
	  d1(i)=0.5*(a+ziota*p1(i)*b)
	  d2(i)=0.5*(a+ziota*p2(i)*b)
280   continue
      do 285 j=1,3
	  t1(j)=d1(j)/c1(j)
	  t2(j)=d2(j)/c2(j)
285   continue

c----------------------------------------------------------------------
c     Calculate the constant uzip
	do 290 i=1,3
	  q2(i)=q21(i)+q22(i)+q23(i)
290   continue
      call zeromat(ztemp1)
	do 300 j=1,3
	  ztemp(j,j)=clog(c2(j))
300   continue
	call matmult(A2,ztemp1,ztemp2)
	call zerovec(zvec)
	do 319 i=1,3
	do 319 j=1,3
	  zvec(i)=zvec(i)+ztemp2(i,j)*q2(j)
319   continue
      do 340 i=1,3
	  uzip(i)=-(1./pi)*aimag(zvec(i))
340   continue
c--------------------------------------------------------------------
c
c     Calculate the Green's function for the outer solid, material 1
      x2=x2beg
c     Form the complex coordinates and their mappings
500	do 520 k=1,3
        z1(k)=x1+p1(k)*x2
	  zterm1=a**2+(p1(k))**2*b**2
	  zterm2=a-ziota*p1(k)*b
	  zterm3=a+ziota*p1(k)*b
	  xi1(k)=(z1(k)+sqrt(z1(k)**2-zterm1))/zterm2
     	  xi1inv(k)=(z1(k)-sqrt(z1(k)**2-zterm1))/zterm3
c     
     	  z10(k)=x10+p1(k)*x20
	  xi10(k)=(z10(k)+sqrt(z10(k)**2-zterm1))/zterm2
     	  xi10inv(k)=(z10(k)-sqrt(z10(k)**2-zterm1))/zterm3
520   continue
c
c     Form diagonal matrices except for k-series terms
      call zeromat(diag)
	call zeromat(diag1)
	call zeromat(diag2)
	call zeromat(diag3)     
     	do 540 j=1,3
	  zterm1=clog(xi1(j)-xi10(j))
	  diag(j,j)=zterm1
     	  zterm1=clog(xi1inv(j)-conjg(xi10inv(1)))
	  diag1(j,j)=zterm1
	  zterm1=clog(xi1inv(j)-conjg(xi10inv(2)))
     	  diag2(j,j)=zterm1
     	  zterm1=clog(xi1inv(j)-conjg(xi10inv(3)))
     	  diag3(j,j)=zterm1
540   continue     
c
c     Form the u1-displacement GF without the k-series terms
      call matmult(A1,diag,ztemp)
	call matmult(A1,diag1,ztemp1)
	call matmult(A1,diag2,ztemp2)     
     	call matmult(A1,diag3,ztemp3)
	call zerovec(zvec)
	call zerovec(zvec1)
	call zerovec(zvec2)
	call zerovec(zvec3)
     	do 560 i=1,3
	do 560 j=1,3
	  zvec(i)=zvec(i)+ztemp(i,j)*q1inf(j)
     	  zvec1(i)=zvec1(i)+ztemp1(i,j)*q11(j)
     	  zvec2(i)=zvec2(i)+ztemp2(i,j)*q12(j)
     	  zvec3(i)=zvec3(i)+ztemp3(i,j)*q13(j)
560   continue
      do 580 k=1,3
	  u1(k)=(1./pi)*aimag(zvec(k)+zvec1(k)+zvec2(k)+zvec3(k))
580   continue
c
c-------------------------------------------------------------------
c     Now add the k-series terms to the GF.           
     	do 888 idx=1,kterms
c       First thing to do is to solve for the vector gk.
 	  call zeromat(diag1)
	  call zeromat(diag2)
	  call zeromat(diag3)
	  do 700 j=1,3
	    ztau21(j)=t2(j)*xi10(1)
		ztau22(j)=t2(j)*xi10(2)
		ztau23(j)=t2(j)*xi10(3)
700     continue 
        do 720 i=1,3
	    diag1(i,i)=(ztau21(i))**idx
	    diag2(i,i)=(ztau22(i))**idx
	    diag3(i,i)=(ztau23(i))**idx
720     continue
	  call matmult(B2,diag1,ztemp1)
	  call matmult(B2,diag1,ztemp2)
	  call matmult(B2,diag1,ztemp3)
	  call zerovec(yk)
	  do 740 j=1,3
	  do 740 k=1,3
	    yk(j)=yk(j)+ztemp1(j,k)*q21(k)
740	  continue
        do 760 j=1,3
	  do 760 k=1,3
	    yk(j)=yk(j)+ztemp2(j,k)*q22(k)
760	  continue
        do 780 j=1,3
	  do 780 k=1,3
	    yk(j)=yk(j)+ztemp3(j,k)*q23(k)
780	  continue
        do 790 i=1,3
	  do 790 j=1,3
	    ztemp1(i,j)=conjg(M1inv(i,j))-conjg(M2inv(i,j))
790     continue
        call matmult(B2t,ztemp1,ztemp2)
	  call zeromat(diag)
	  do 800 i=1,3
	    diag(i,i)=(t2(i))**idx
800     continue
        call matmult(diag,ztemp2,ztemp3)
	  call matmult(B2,ztemp3,bigyk)
	  call matmult(bigyk,czXinv,ztemp1)
	  call zerovec(zvec)
	  do 805 j=1,3
	  do 805 k=1,3
	    zvec(j)=zvec(j)+ztemp1(j,k)*conjg(yk(k))
805     continue
        do 810 j=1,3
	  do 810 k=1,3
	    zvec(j)=ziota*zvec(j)-yk(j)
810     continue
	  do 815 i=1,3
	  do 815 j=1,3
	    cbigyk(i,j)=conjg(bigyk(i,j))
815     continue
        call matmult(czXinv,cbigyk,ztemp1)
	  call matmult(bigyk,ztemp1,ztemp2)
	  do 820 i=1,3
	  do 820 j=1,3
	    ztemp3(i,j)=zX(i,j)-ztemp2(i,j)
820     continue
        call matmult(ztemp3,B1,ztemp1)
	  call cinv(ztemp1,ztemp2)
	  call zerovec(gk)
	  do 825 i=1,3
	  do 825 j=1,3
	    gk(i)=gk(i)+ztemp2(i,j)*zvec(j)
825     continue
c
c       Now the additional terms to the displacement GF can be determined:
	  call zeromat(diag)
	  xk=float(idx)
	  do 900 i=1,3
	    zterm1=xi1(i)
	    diag(i,i)=1.0/(zterm1**xk) 
900     continue
        call matmult(A1,diag,ztemp)
	  call zerovec(zvec)
	  do 920 j=1,3
	  do 920 k=1,3
	    zvec(j)=zvec(j)+ztemp(j,k)*gk(k)
920     continue
        do 940 m=1,3
	    u1(m)=u1(m)+(1.0/pi)*aimag(zvec(m)/xk)
940     continue
c
c     Go get next term in k-series.
888   continue
c
c     Write out displacment components, then increment to next x2 value.
      if(u1(3).ne.0.0) then
	  write(6,*)'Problem!  u3 not equal to zero!'
	  write(6,*)'***Ending exectution***'
	  goto 999
	endif
     	write(20,*)x1,x2,u1(1),u1(2)
      if(x2.ge.x2end) goto 999
	x2=x2+xdiv
	goto 500     
c
c
999   close(20)
      stop
	end
c

c *********************************************************************
c *===================================================================*
c *========================= SUBROUTINES =============================*
c *===================================================================*
c *********************************************************************
c
      subroutine zeromat(cc)
c     This subroutines zeroes all elements of the 3x3 matrix cc
      complex cc(3,3)
	  do 10 i=1,3
	  do 10 j=1,3
	    cc(i,j)=cmplx(0.0,0.0)
10    continue
      return
	  end
c********************************************************************
c
      subroutine matmult(aa,bb,cc)
c     This subroutines multiplies two 3x3 matrices, cc=aa*bb
      complex aa(3,3),bb(3,3),cc(3,3)
	call zeromat(cc)
	do 10 i=1,3
	do 10 j=1,3
	do 10 k=1,3
	  cc(i,j)=cc(i,j)+aa(i,k)*bb(k,j)
10    continue
      return
	end
c*******************************************************************
c
      subroutine ctrans(aa,aat)
c     This subroutines computes the transpose of a 3x3 matrix aa
      complex aa(3,3),aat(3,3)
	do 10 i=1,3
	do 10 j=1,3
	  aat(i,j)=aa(j,i)
10    continue
      return
	end
c*********************************************************************
c
      subroutine cinv(a,ainv)
c     This subroutine calculates the inverse of a 3x3 complex matrix aa
      complex a(3,3),ainv(3,3),cdet,aa,bb,cc,cof(3,3)
	aa=a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2))
      bb=a(2,1)*(a(1,2)*a(3,3)-a(1,3)*a(3,2))
      cc=a(3,1)*(a(1,2)*a(2,3)-a(1,3)*a(2,2))
	cdet=aa-bb+cc
c	
	cof(1,1)=a(2,2)*a(3,3)-a(2,3)*a(3,2)
	cof(1,2)=a(2,1)*a(3,3)-a(2,3)*a(3,1)
	cof(1,3)=a(2,1)*a(3,2)-a(2,2)*a(3,1)
c
	cof(2,1)=a(1,2)*a(3,3)-a(1,3)*a(3,2)
	cof(2,2)=a(1,1)*a(3,3)-a(1,3)*a(3,1)
	cof(2,3)=a(1,1)*a(3,2)-a(1,2)*a(3,1)
c
      cof(3,1)=a(1,2)*a(2,3)-a(1,3)*a(2,2)
	cof(3,2)=a(1,1)*a(2,3)-a(1,3)*a(2,1)
	cof(3,3)=a(1,1)*a(2,2)-a(1,2)*a(2,1)
c
      do 100 i=1,3
	do 100 j=1,3
	  ainv(i,j)=(-1.0)**(i+j)*cof(j,i)
100   continue
      do 150 i=1,3
	do 150 j=1,3
	  ainv(i,j)=ainv(i,j)/cdet
150   continue
	return
	end
c*********************************************************************
c
      subroutine zibeta(zi1,zi2,zi3)
c     This subroutine forms the I-beta matrices
      complex zi1(3,3),zi2(3,3),zi3(3,3)
	call zeromat(zi1)
	call zeromat(zi2)
	call zeromat(zi3)
	zi1(1,1)=cmplx(1.0,0.0)
	zi2(2,2)=cmplx(1.0,0.0)
	zi3(3,3)=cmplx(1.0,0.0)
	return
	end
c*********************************************************************
c
      subroutine zerovec(bb)
c     This subroutine zeroes out a complex 3-vector bb
      complex bb(3)
	do 10 i=1,3
	  bb(i)=cmplx(0.0,0.0)
10    continue
      return
	end
c**********************************************************************
c
      subroutine orthog(a,b)
c     This subroutine checks the orthogonality and closure relations
c     on the matrices A and B.  The results of the check are printed
c     to the file orthog.chk
      complex a(3,3),b(3,3),at(3,3),bt(3,3)
	complex ca(3,3),cb(3,3),cat(3,3),cbt(3,3)
	complex ztemp1(3,3),ztemp2(3,3),check1(3,3),check2(3,3)
c
	open(unit=33,file='orthog.chk',status='unknown')
	write(33,*)'**************************************************'
c
c     Write the A and B matrices to the output file
      write(33,*) 'A Matrix'
	write(33,*) '========'
	write(33,*) a(1,1), a(1,2), a(1,3)
	write(33,*) a(2,1), a(2,2), a(2,3)
	write(33,*) a(3,1), a(3,2), a(3,3)
	write(33,*)
	write(33,*) 'B Matrix'
	write(33,*) '========'
	write(33,*) b(1,1), b(1,2), b(1,3)
	write(33,*) b(2,1), b(2,2), b(2,3)
	write(33,*) b(3,1), b(3,2), b(3,3)
	write(33,*)
	write(33,*)
      call ctrans(a,at)
	call ctrans(b,bt)
	call zconj(a,ca)
	call zconj(b,cb)
	call zconj(at,cat)
	call zconj(bt,cbt)
c
c     Check the orthogonality relations
      write(33,*)'Orthogonality relations check'
	write(33,*)'============= ========= ====='
	write(33,*)'The next two matrices should equal the identity'
	write(33,*)
	call matmult(bt,a,ztemp1)
	call matmult(at,b,ztemp2)
	do 100 i=1,3
	do 100 k=1,3
	  check1(i,k)=ztemp1(i,k)+ztemp2(i,k)
100   continue
      do 105 i=1,3
	   write(33,*)check1(i,1),check1(i,2),check1(i,3)
105   continue
      write(33,*)
	write(33,*)
      call matmult(cbt,ca,ztemp1)
	call matmult(cat,cb,ztemp2)
	do 120 i=1,3
	do 120 k=1,3
	  check2(i,k)=ztemp1(i,k)+ztemp2(i,k)
120   continue  	
      do 125 i=1,3
	   write(33,*)check2(i,1),check2(i,2),check2(i,3)
125   continue
      write(33,*)
	write(33,*)
	write(33,*)
	write(33,*)'The next two matrices should be null'
 	call matmult(bt,ca,ztemp1)
	call matmult(at,cb,ztemp2)
	do 130 i=1,3
	do 130 k=1,3
	  check1(i,k)=ztemp1(i,k)+ztemp2(i,k)
130   continue
      do 135 i=1,3
	   write(33,*)check1(i,1),check1(i,2),check1(i,3)
135   continue
      write(33,*)
	write(33,*)
      call matmult(cbt,a,ztemp1)
	call matmult(cat,b,ztemp2)
	do 140 i=1,3
	do 140 k=1,3
	  check2(i,k)=ztemp1(i,k)+ztemp2(i,k)
140   continue 	
      do 145 i=1,3
	   write(33,*)check2(i,1),check2(i,2),check2(i,3)
145   continue
      write(33,*)
	write(33,*)
	write(33,*)
	write(33,*)
	write(33,*)'Closure relations check'
	write(33,*)'======= ========= ====='
	write(33,*)'The next two matrices should equal the identity'
	write(33,*)
	call matmult(a,bt,ztemp1)
	call matmult(ca,cbt,ztemp2)
	do 200 i=1,3
	do 200 k=1,3
	  check1(i,k)=ztemp1(i,k)+ztemp2(i,k)
200   continue
      do 205 i=1,3
	   write(33,*)check1(i,1),check1(i,2),check1(i,3)
205   continue
      write(33,*)
      call matmult(b,at,ztemp1)
	call matmult(cb,cat,ztemp2)
	do 220 i=1,3
	do 220 k=1,3
	  check2(i,k)=ztemp1(i,k)+ztemp2(i,k)
220   continue 	
      do 225 i=1,3
	   write(33,*)check2(i,1),check2(i,2),check2(i,3)
225   continue
      write(33,*)
	write(33,*)
	write(33,*)
	write(33,*)'The next two matrices should be null'
 	call matmult(a,at,ztemp1)
	call matmult(ca,cat,ztemp2)
	do 230 i=1,3
	do 230 k=1,3
	  check1(i,k)=ztemp1(i,k)+ztemp2(i,k)
230   continue
      do 235 i=1,3
	   write(33,*)check1(i,1),check1(i,2),check1(i,3)
235   continue
      write(33,*)
      call matmult(b,bt,ztemp1)
	call matmult(cb,cbt,ztemp2)
	do 240 i=1,3
	do 240 k=1,3
	  check2(i,k)=ztemp1(i,k)+ztemp2(i,k)
240   continue 	
      do 245 i=1,3
	   write(33,*)check2(i,1),check2(i,2),check2(i,3)
245   continue
      write(33,*)
	write(33,*)
	write(33,*)'****************************************************'
	return
	end
c********************************************************************
c
      subroutine zconj(a,ca)
c     This subroutine forms the conjugate of the matrix a
      complex a(3,3),ca(3,3)
	do 100 i=1,3
	do 100 j=1,3
	  ca(i,j)=conjg(a(i,j))
100   continue
      return
	end
c********************************************************************
