C
C     One-purpose code to compare the finite-difference and analytical
C     solutions in the twisted-crystal model.
C
C     The program reads the output files TCFDX='tc-fdx.dat' and
C     TCFDY='tc-fdx.dat' of Matlab scripts 'tc-fdx.m' and 'tc-fdy.m',
C     calculated for relatively simple and well-defined initial
C     conditions by a finite difference method.
C     The program reads the output file TCSTRESS='stress.out' of program
C     'tcgreen.for', containing the frequency-dependent initial stress
C     corresponding to the analytical solution (initial displacement
C     is constant).
C     The program then converts the given finite-difference solution to
C     the finite-difference solution with the initial conditions
C     corresponding to the analytical solution, and writes it into
C     output file TCFD='tc-fd.out'.
C     The program also reads analytical solution GREENTC='tc-anal.out'
C     calculated by program 'tcgreen.for', and writes the relative
C     difference between finite-difference solution TCFD='tc-fd.out' and
C     analytical solution GREENTC='tc-anal.out' into output file
C     TCCOMP='tc-fderr.out'.
C
C     Date: 2003, May 26
C     Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      CHARACTER*80 FILSEP,FILFDX,FILFDY,FILSTR,FILFD,FILTC,FILERR
      PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4,LU5=5,LU6=6)
C
      READ(*,*) FILSEP
      CALL RSEP1(1,FILSEP)
      CALL RSEP3T('TCFDX'   ,FILFDX,'tc-fdx.dat')
      CALL RSEP3T('TCFDY'   ,FILFDY,'tc-fdy.dat')
      CALL RSEP3T('TCSTRESS',FILSTR,'stress.out')
      CALL RSEP3T('TCFD'    ,FILFD ,'tc-fd.out')
      CALL RSEP3T('GREENTC' ,FILTC ,'tc-anal.out')
      CALL RSEP3T('TCCOMP'  ,FILERR,'tc-fderr.out')
C
      OPEN(LU1,FILE=FILFDX)
      OPEN(LU2,FILE=FILFDY)
      OPEN(LU3,FILE=FILSTR)
      OPEN(LU4,FILE=FILFD )
      OPEN(LU5,FILE=FILTC )
      OPEN(LU6,FILE=FILERR)
C
   10 CONTINUE
        READ(LU1,*,END=90) F1,G11R,G11Y,G21R,G21Y
        READ(LU2,*,END=90) F2,G12R,G12Y,G22R,G22Y
        READ(LU3,*,END=90) F3,S11R,S11Y,S12R,S12Y,S22R,S22Y
        READ(LU5,*,END=90) F3,U11R,U11Y,U21R,U21Y,U31R,U31Y,
     *                        U12R,U12Y,U22R,U22Y
        IF(F1.NE.F3) THEN
          WRITE(*,*) F1,F3
          STOP 1
        END IF
        IF(F2.NE.F3) THEN
          WRITE(*,*) F2,F3
          STOP 2
        END IF
        V11R=G11R+G11Y*S11R+G12Y*S12R
        V11Y=     G11Y*S11Y+G12Y*S12Y
        V21R=G21R+G21Y*S11R+G22Y*S12R
        V21Y=     G21Y*S11Y+G22Y*S12Y
        V12R=G12R+G11Y*S12R+G12Y*S22R
        V12Y=     G11Y*S12Y+G12Y*S22Y
        V22R=G22R+G21Y*S12R+G22Y*S22R
        V22Y=     G21Y*S12Y+G22Y*S22Y
        WRITE(LU4,'(9(F9.6,A))')
     *    F3,' ',V11R,' ',V11Y,' ',V21R,' ',V21Y,' 0 0 '
     *          ,V12R,' ',V12Y,' ',V22R,' ',V22Y,' 0 0  0 0  0 0  0 0'
        DIF=(V11R-U11R)**2+(V21R-U21R)**2+(V12R-U12R)**2+(V22R-U22R)**2
     *     +(V11R-U11R)**2+(V21R-U21R)**2+(V12R-U12R)**2+(V22R-U22R)**2
        DIF=SQRT(0.5*DIF)
        WRITE(LU6,'(9(F9.6,1X))') F3,DIF
      GO TO 10
C
   90 CONTINUE
      CLOSE(LU1)
      CLOSE(LU2)
      CLOSE(LU3)
      CLOSE(LU4)
      CLOSE(LU5)
      CLOSE(LU6)
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'sep.for'
      INCLUDE 'error.for'
      INCLUDE 'length.for'
C
C=======================================================================
C