PROGRAM fco
C     compare CP and SOS tensors
      IMPLICIT INTEGER*4 (I-N) 
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (NT0=500,MX3=3*NT0)
      DIMENSION alp(MX3,3,3),alc(MX3,3,3)
      DIMENSION avp(MX3,3,3),avc(MX3,3,3)
      DIMENSION gp(MX3,3,3),gc(MX3,3,3)
      DIMENSION ap(MX3,3,3,3),ac(MX3,3,3,3)
c
      open(34,file='X.TTT',status='old')
      read(34,*)
      read(34,*)nat
c
      read(34,*)
      read(34,*)
      do 10 i=1,3
      read(34,*)
      do 10 l=1,nat
      do 10 ix=1,3
      IIND=3*(l-1)+ix
10    read(34,*)idum,idum,(alp(IIND,i,j),j=1,3)
c
      read(34,*)
      read(34,*)
      do 101 i=1,3
      read(34,*)
      do 101 l=1,nat
      do 101 ix=1,3
      IIND=3*(l-1)+ix
101   read(34,*)idum,idum,(gp(IIND,i,j),j=1,3)
c
      read(34,*)
      read(34,*)
      do 102 k=1,3
      do 102 i=1,3
      read(34,*)
      do 102 l=1,nat
      do 102 ix=1,3
      IIND=3*(l-1)+ix
102   read(34,*)idum,idum,(ap(IIND,i,j,k),j=1,3)
c
      read(34,*)
      read(34,*)
      do 103 i=1,3
      read(34,*)
      do 103 l=1,nat
      do 103 ix=1,3
      IIND=3*(l-1)+ix
103   read(34,*)idum,idum,(avp(IIND,i,j),j=1,3)
      close(34)
      write(6,*)'X.TTT read',nat,' atoms'
c
      open(34,file='Y.TTT',status='old')
      read(34,*)
      read(34,*)nat
c
      read(34,*)
      read(34,*)
      do 11 i=1,3
      read(34,*)
      do 11 l=1,nat
      do 11 ix=1,3
      IIND=3*(l-1)+ix
11    read(34,*)idum,idum,(alc(IIND,i,j),j=1,3)
c
      read(34,*)
      read(34,*)
      do 111 i=1,3
      read(34,*)
      do 111 l=1,nat
      do 111 ix=1,3
      IIND=3*(l-1)+ix
111   read(34,*)idum,idum,(gc(IIND,i,j),j=1,3)
c
      read(34,*)
      read(34,*)
      do 112 k=1,3
      do 112 i=1,3
      read(34,*)
      do 112 l=1,nat
      do 112 ix=1,3
      IIND=3*(l-1)+ix
112   read(34,*)idum,idum,(ac(IIND,i,j,k),j=1,3)
c
      read(34,*)
      read(34,*)
      do 113 i=1,3
      read(34,*)
      do 113 l=1,nat
      do 113 ix=1,3
      IIND=3*(l-1)+ix
113   read(34,*)idum,idum,(avc(IIND,i,j),j=1,3)
      close(34)
      write(6,*)'Y.TTT read',nat,' atoms'
      n=3*nat*3*3
      a1=0.0d0
      a2=0.0d0
      a3=0.0d0
      a4=0.0d0
      a1=0.0d0
      d2=0.0d0
      d3=0.0d0
      d4=0.0d0
      xiyi=0.0d0
      yiyi=0.0d0
      xixi=0.0d0
      do 12 i=1,3
      do 12 l=1,nat
      do 12 ix=1,3
      IIND=3*(l-1)+ix
      do 12 j=1,3
      xixi=xixi+alp(IIND,i,j)*alp(IIND,i,j)
      yiyi=yiyi+alc(IIND,i,j)*alc(IIND,i,j)
12    xiyi=xiyi+alc(IIND,i,j)*alp(IIND,i,j)
      if(n.gt.0.and.abs(xixi).gt.1.0d-10)then
       a1=xiyi/xixi
       d1=sqrt((yiyi-2.0d0*a1*xiyi+a1*a1*xixi)/real(n))
      endif
      xiyi=0.0d0
      yiyi=0.0d0
      xixi=0.0d0
      do 13 i=1,3
      do 13 l=1,nat
      do 13 ix=1,3
      IIND=3*(l-1)+ix
      do 13 j=1,3
      xixi=xixi+gp(IIND,i,j)*gp(IIND,i,j)
      yiyi=yiyi+gc(IIND,i,j)*gc(IIND,i,j)
13    xiyi=xiyi+gc(IIND,i,j)*gp(IIND,i,j)
      if(n.gt.0.and.abs(xixi).gt.1.0d-10)then
       a2=xiyi/xixi
       d2=sqrt((yiyi-2.0d0*a2*xiyi+a2*a2*xixi)/real(n))
      endif
      xiyi=0.0d0
      yiyi=0.0d0
      xixi=0.0d0
      do 14 k=1,3
      do 14 i=1,3
      do 14 l=1,nat
      do 14 ix=1,3
      IIND=3*(l-1)+ix
      do 14 j=1,3
      xixi=xixi+ap(IIND,i,j,k)*ap(IIND,i,j,k)
      yiyi=yiyi+ac(IIND,i,j,k)*ac(IIND,i,j,k)
14    xiyi=xiyi+ac(IIND,i,j,k)*ap(IIND,i,j,k)
      na=3*nat*3*3*3
      if(na.gt.0.and.abs(xixi).gt.1.0d-10)then
       a3=xiyi/xixi
       d3=sqrt((yiyi-2.0d0*a3*xiyi+a3*a3*xixi)/real(na))
      endif
      xiyi=0.0d0
      yiyi=0.0d0
      xixi=0.0d0
      do 15 i=1,3
      do 15 l=1,nat
      do 15 ix=1,3
      IIND=3*(l-1)+ix
      do 15 j=1,3
      xixi=xixi+avp(IIND,i,j)*avp(IIND,i,j)
      yiyi=yiyi+avc(IIND,i,j)*avc(IIND,i,j)
15    xiyi=xiyi+avc(IIND,i,j)*avp(IIND,i,j)
      if(n.gt.0.and.abs(xixi).gt.1.0d-10)then
       a4=xiyi/xixi
       d4=sqrt((yiyi-2.0d0*a4*xiyi+a4*a4*xixi)/real(n))
      endif
      write(6,6000)a1,d1,a2,d2,a3,d3,a4,d4
6000  format('al y = a x; a = ',f10.5,' (del = ',f10.5,')',/,
     1       ' g y = a x; a = ',f10.5,' (del = ',f10.5,')',/,
     1       ' a y = a x; a = ',f10.5,' (del = ',f10.5,')',/,
     1       'av y = a x; a = ',f10.5,' (del = ',f10.5,')',/)
      stop
      end