C                                                                            C
C       P R O G R A M    M I S O R I E N T A T I O N   M A T R I X           C
C                                                                            C
C       Written by P. Heilmann                          03-Feb-82            C
C                  Ohio State University                                     C
C                  Met. Engineering                                          C
C                                                                            C
C       Loosely based on                                                     C
C                  C. T. Young                          20-DEC-73            C
C                  Virginia Polytechnic Institute and State University       C
C                  Met. Engineering                                          C
C                                                                            C
C       Published by                                                         C
C                  C.T. Young, J.H. Steele,Jr., and J.L. Lytton              C
C                  Met. Trans., 1973, Vol.4, 2081-2089                       C
C                                                                            C
C       INPUT:                                                               C
C               Input from keyboard                                          C
C                 Beam direction (real)                                      C
C                 Kikuchi line (vector pointing away from the origin         C
C                 Rot. Angle gamma from pattern to reference frame           C
C                      (in deg, counterclockwise = positive)                 C
C               or                                                           C
C               Input from File, written by program ORIENT                   C
C                 Rotation matrix from crystal to reference frame            C
C                      (two files can be open at the same time)              C
C                                                                            C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C234567890223456789032345678904234567890523456789062345678907234567890
 
        PROGRAM MISMAT
        DIMENSION NPAT(8,2),DM(3,3,24),CR(3,3,2),RMAT(3,3),
     1    RMATS(3,3)
	  Double Precision  range(10,2)
	  Real              rr(10,3,3)
        CHARACTER*64 infile,outfile
	  character    new,old,answer
	  character*7  progsource
	  character*1  cracked
	  integer      Eof, status, grainA, grainB
	  logical      bothfound
C
	Parameter   (Eof=-1)
C
	  Common /ExtRots/ rr(10,3,3),range(10,2),numex
c
  552   FORMAT (1H1)
  553   FORMAT (/'  PROGRAM TO CALCULATE THE MISORIENTATION',
     1          'MATRIX BETWEEN TWO CRYSTALS'/)
  555   FORMAT (' Enter name of input file containing beam
     1 directions and rotation matrices')
  556   FORMAT (' Enter name of transfer output file')
  557   FORMAT (' Continue with min. angle (sym. op. # ',
     1         I2,')   (Y) ? ',\)
  558   FORMAT (' Which other symm. op.     (none = 1) ? ',\)
 560  format(8(2A))
  561   FORMAT (' Input from keyboard              (Y) ? ',\)
C
	open (8,file = 'pairfile',status = old)
	CALL SYMMAT(DM)
	NU=2
	NUM=NU
C
	WRITE (9,553)
C
	WRITE (9,555)
	read(9,560) infile
   3  write(9,*) ' enter name of transfer output file '
	read(9,5) outfile
   5  format(8(2A))
	write(9,*) ' new transfer output file ? (n)   or'
	write(9,*) '             or'
	write(9,*) '  append an existing file ? (a)'
	read(9,*) answer
	if(answer.eq.'n') then
	  Open(unit=4,file=outfile,status=new,position='asis')
	else if(answer.eq.'a') then
	  Open(unit=4,file=outfile,status=old,position='append')
	else
	  goto 3
	End if
c
	answer = 'none'
 115  write(9,*) ' Detailed output will go into a file named',
     x           ' MisMat.out.'
	write(9,*) ' Is this to be a new output file ? (n)   or'
	write(9,*) '    should we'
	write(9,*) 'append an existing MisMat.out file ? (a)'
	read(9,*) answer
	if(answer.eq.'n') then
	  Open(unit=2,file='MisMat.out',status='new ')
	else if(answer.eq.'a') then
	  Open(unit=2,file='MisMat.out',status='old',
     x              position='append')
	else
	  goto 115
	End if
c
c	   prompt user for any relative rotations that should 
c	   be performed on grains within any ranges
c        of grain numbers
c
	numex = 0
	answer = 'none'
	write(9,*)
	write(9,*) ' Do any of the rotation matrices need',
     1           ' further rotation ? (y/n)'
	read(9,*) answer
	if (answer.eq.'y') call extrot(infile)
c
	Do While(status.ne.Eof)
	  cracked = 'u'                ! for unknown
	  read(8,*,iostat=status) grainA,grainB,cracked
	  write(9,*) ' grainA = ',grainA,' grainB = ',grainB
	  call infile(CR,grainA,grainB,NU,infile,
     1              progsource,bothfound)

c	  write(9,*) ' both grain numbers found : ',bothfound
	  if (bothfound) then
C
C       Invert rot matrix of second crystal
C
         CALL INVERT (CR,2)
C
C       Multiply both matrices
C
         CALL MMPRD (CR,2,CR,1,RMAT,1)
C
C       Apply all 24 symm. operations to find
C       minimum rot angle
C
	   NUM = 9
 31      CALL DEGMAT (RMAT,DM,LDM,NUM)
c        WRITE (*,557) LDM
c        if(IANS('Y').eq.1) goto 35
C
C       Apply min. ang. symm. op. (LDM)
C
 35      CALL MMPRD (DM,LDM,RMAT,1,RMATS,1)
C
C       Print results
C
         CALL OUTPUT (RMATS,LDM,DM,NU,grainA,grainB,
     1               progsource,cracked)

	  End if            ! bothfound block
  100 repeat
C
C
      close(2)
	close(4)
	close(8)
	STOP
      END
C
C
C-------DATA INPUT FROM KEYBOARD------------------------------------
C
        SUBROUTINE INKBRD(CR,NPAT,NCRY,NU)
        DIMENSION NPAT(8,2),CR(3,3,2),IHKL(3),BEAMD(3)
C
  502   FORMAT (/' Input data of',I2,'. crystal :')
  511   FORMAT (' Plate number         (16 char. max.) ? ',\)
  513   FORMAT (' Beam direction  X,Y,Z         (real) ? ',\)
  514   FORMAT (' K. - line  U,V,W           (integer) ? ',\)
  517   FORMAT (' Angle betw. K.-line a. x-axis  (deg) ? ',\)
  518   FORMAT (/' Do you want to make any changes  (N) ?',\)
  520   FORMAT (8A2)
  609   FORMAT (/' Input data of ',I1,'. crystal:'/)
  610   FORMAT (' Plate #',T21,'Beam direction',T58,'K.-line',
     1         T71,'Gamma')
  611   FORMAT (' ',8A2,3F12.7,1X,3I4,2X,F7.2/)
  612   FORMAT (T21,' Rotation matrix crystal - ref. frame')
  613   FORMAT (3(T18,3F12.7,/))
C
  1     WRITE (*,502) NCRY
        WRITE (*,511)
        READ (*,520) (NPAT(I,NCRY),I=1,8)
        WRITE (*,513)
        READ (*,*) (BEAMD(I),I=1,3)
        WRITE (*,514)
        READ (*,*) (IHKL(I),I=1,3)
        WRITE (*,517)
        READ (*,*) GAMMA
        WRITE (*,518)
        IF (IANSW('N') .EQ. 0) GOTO 1
C
C
        CALL ROTMAT (BEAMD,IHKL,GAMMA*0.1745329E-1,CR,NCRY)
C
C      Print out input and rot. matrices (ref., spec. frames)
C
        WRITE (NU,609) NCRY
        WRITE (NU,610)
        WRITE (NU,611) (NPAT(I,NCRY),I=1,8),BEAMD,IHKL,GAMMA
        WRITE (NU,612)
        WRITE (NU,613) ((CR(K,I,NCRY),I=1,3),K=1,3)
C
        RETURN
        END
C
C-------DATA INPUT FROM FILE-----------------------------------
C
	  subroutine infile(CR,grainA,grainB,nu,infile,
     1                    progsource,bothfound)
        Real              CR(3,3,2)
	  logical           Afound, Bfound, bothfound
	  integer           Eof3, nu, status3, ne, numex
	  integer           gn, grainA, grainB, gn1, gn2
	  character*64      infile
	  character         old
	  character*7       psA,psB,ps,progsource
	  Double Precision  grainNum,range(10,2),grainNumA,grainNumB
	  Real              rr(10,3,3),temp(3,3)
C
	  Parameter    (Eof3 = -1)
C
	  Common /ExtRots/ rr(10,3,3),range(10,2),numex
c
  13	format(f8.3,7a)
  14	format(' grainNum = ',f7.3)
  507	FORMAT (' '/' No combination for grain ',i4,
     1     ' and grain ',i4/' ')
C
	ne = numex
	status3 = 0
	Bfound = .false.
	Afound = .false.
	bothfound = .false.
	open (3,file = infile,status = old)
	Do While(status3.ne.Eof3)
  1	  Read(3,13,iostat=status3) grainNum,ps
	  write(9,*) ' ps = ',ps
	  if ((ps.ne.'Orient').and.(ps.ne.'Polser')) ps = 'unknown'
c	  write(9,14) grainNum
	  gn1 = grainNum
	  gn2 = nint(grainNum*1000)
	  gn = gn2 -1000*gn1
	  write(9,*) ' gn,a,b ',gn,grainA,grainB
	  if(gn.eq.grainA) then
	    grainNumA = grainNum
	    read(3,*,iostat=status3) camcon,caml
	    Do 10 jj = 1,3
  10        read(3,*,iostat=status3) (CR(jj,kk,1),kk=1,3)
	      if ((ne).ne.0) then
	        Do (jj = 1,ne)
		    if ((grainNum.ge.range(jj,1)).and.(grainNum
     1             .le.range(jj,2))) then
	           call MatMult(rr,jj,CR,1,temp)
		     CR(1,1,1) = temp(1,1)
		     CR(1,2,1) = temp(1,2)
		     CR(1,3,1) = temp(1,3)
		     CR(2,1,1) = temp(2,1)
		     CR(2,2,1) = temp(2,2)
		     CR(2,3,1) = temp(2,3)
		     CR(3,1,1) = temp(3,1)
		     CR(3,2,1) = temp(3,2)
		     CR(3,3,1) = temp(3,3)
		     call NORMAT(CR,1)
		    end if
	        End do
		End if
	    psA = ps
	    Afound = .true.
	  End if
	  if(gn.eq.grainB) then
	    grainNumB = grainNum
	    read(3,*,iostat=status3) camcon,caml
	    Do 12 jj = 1,3
  12        read(3,*,iostat=status3) (CR(jj,kk,2),kk=1,3)
	      if ((ne).ne.0) then
	        Do (jj = 1,ne)
		    if ((grainNum.ge.range(jj,1)).and.(grainNum
     1             .le.range(jj,2))) then
	           call MatMult(rr,jj,CR,2,temp)
		     CR(1,1,2) = temp(1,1)
		     CR(1,2,2) = temp(1,2)
		     CR(1,3,2) = temp(1,3)
		     CR(2,1,2) = temp(2,1)
		     CR(2,2,2) = temp(2,2)
		     CR(2,3,2) = temp(2,3)
		     CR(3,1,2) = temp(3,1)
		     CR(3,2,2) = temp(3,2)
		     CR(3,3,2) = temp(3,3)
		     call NORMAT(CR,2)
		    end if
	        End do
		End if
	    psB = ps
	    Bfound = .true.
	  End if
	  if ((Afound).and.(Bfound)) then
  21	    format(f7.3)
	    write(4,21) grainNumA
	    write(4,21) grainNumB
	    if ((psA.eq.'Orient').and.(psB.eq.'Orient')) then
	      progsource = 'Orient'
	    else if ((psA.eq.'Polser').or.(psB.eq.'Polser')) then
	      progsource = 'Polser'
	    else
	      progsource = 'unknown'
	    End if
	    bothfound = .true.
	    close(3)
	    return
	  else
	    if ((gn.ne.grainA).and.(gn.ne.grainB)) then
	      read(3,*,iostat=status3) a,b
	      read(3,*,iostat=status3) a,b,c
	      read(3,*,iostat=status3) a,b,c
	      read(3,*,iostat=status3) a,b,c
	    End if
	  End if
	repeat
	bothfound = .false.     ! redundant, but just to be certain
	write(2,507) grainA,grainB
	write(9,507) grainA,grainB
	close(3)
	return
	end
C 
C-------CALCULATE ALL DEGENERATE MISORIENTATIONS-----------------
C
        SUBROUTINE DEGMAT(RMAT,DM,LDM,NU)
        DIMENSION DM(3,3,24),RMAT(3,3),AMM(3,3),AXMS(3)
C
  500   FORMAT (/' Degenerate misorientation of the two ',
     1           'crystals:'//4(' Sym.Op. Mis.Angle '))
  501   FORMAT (I4,F12.3,4X,\)
  502   FORMAT (' ')
  503   FORMAT (/' Minimum angle is ',F7.3,' deg.')
  622   FORMAT (//' DEGENERATE MISORIENTATIONS  (24 symm.',
     1            ' operations used)'//,2X,'#',T12,
     2            'Mis. Angle (deg)',T33,'Misorientation',
     3            ' axis in Crystal frame 1',/)
  623   FORMAT (1X,I2,T12,F9.2,T33,3F10.6)
  624   FORMAT ('1')
C
        CF = 3.141593/180.
        IF (NU .EQ. 2) WRITE (NU,622)
        ANGMIN = 180.0
        WRITE (*,500)
        DO 1 N=1,24
        CALL MMPRD  (DM,N,RMAT,1,AMM,1)
	   arg = 0.5*(AMM(1,1)+AMM(2,2)+AMM(3,3)-1.)
	   if ((abs(arg).gt.1.00).and.(abs(arg).lt.1.005)) then
	     arg = arg / abs(arg)
	   end if
         ANGM=ACOS(arg)/CF
c       ANGM=ACOS(0.5*(AMM(1,1)+AMM(2,2)+AMM(3,3)-1.))/CF
        WRITE (*,501) N,ANGM
        IF (MOD(N,4) .EQ. 0) WRITE (*,502)
        IF (ANGM .GE. ANGMIN) GOTO 4
        ANGMIN = ANGM
        LDM = N
 4      CALL AXIS (AMM,1,AXMS)
        WRITE (NU,623) N,ANGM,AXMS
        IF (MOD(N,4) .EQ. 0) WRITE (NU,502)
 1      CONTINUE
        WRITE (NU,624)
 3      WRITE (*,503) ANGMIN
        RETURN
        END
CCC
C-------CALCULATE MISORIENTATION PARAMETERS-------------------
C
        SUBROUTINE OUTPUT (RMAT,LDM,DM,NU,grainA,grainB,
     1                     progsource,cracked)
        DIMENSION RMAT(3,3),AXMIS(3),DM(3,3,24)
C
	  character*1    cracked
	  character*7    progsource

 611    FORMAT (/' MisO Pars for boundary ',i7,' -',i7,
     1         T45,'Symm. operation',T65,'|',3I3,' |')
 612    FORMAT (T45,'# ',I2,' used',T65,'|',3I3,' |')
 613    FORMAT (T65,'|',3I3,' |')
 617    FORMAT (5X,'Mis. Angle  (deg.)       ',F10.3)
 618    FORMAT (/5X,'Mis. Axis   in Cryst.Fr.1',3F13.6)
 619    FORMAT (/5X,'Mis. Matrix in Cryst.Fr.1',3F13.6,
     1         2(/30X,3F13.6))
 620    FORMAT (5X,'Equiv. to succ. rot.'/7X,'around x,y,z',
     1         ' axis by',2X,3(F8.3,6X),'deg.')
 621    FORMAT (7X,'around z,y,x axis by',2X,3(F8.3,6X),
     1          'deg.')
 622    FORMAT (7X,'around x,z,y axis by',2X,3(F8.3,6X),
     1          'deg.')
C
        CF = 180./3.141593
        CALL AXIS  (RMAT,1,AXMIS)
        ANGMIS = ACOS(0.5*(RMAT(1,1)+RMAT(2,2)+
     1                RMAT(3,3)-1.))*CF
C
        WRITE (NU,611) grainA,grainB,(IRND(DM(I,1,LDM)),
     1                                            I=1,3)
        WRITE (NU,612) LDM,(IRND(DM(I,2,LDM)),I=1,3)
        WRITE (NU,613) (IRND(DM(I,3,LDM)),I=1,3)
        WRITE (NU,617) ANGMIS
        WRITE (NU,618) AXMIS
        WRITE (NU,619) ((RMAT(I,J),J=1,3),I=1,3)
c
c  send output to transfer output file for later use
c
 710  format(3(f12.7,x),f7.2)
 711	format(3(f12.7,x))
 712	format(3(f12.7,x),7a)
 713	format(3(f12.7,x),a)
	  write(4,710) AXMIS(1),AXMIS(2),AXMIS(3),ANGMIS
	  write(4,711) (RMAT(1,J),J=1,3)
	  write(4,712) (RMAT(2,J),J=1,3),progsource
	  write(4,713) (RMAT(3,J),J=1,3),cracked
c
c	  Call CSLbdry(ANGMIS,AXMIS,RMAT)
c
        CALL PROTAX(RMAT(1,1),RMAT(2,1),RMAT(3,1),
     1             RMAT(3,2),RMAT(3,3),AL,BE,GA)
        WRITE (NU,620) AL,BE,GA
        CALL PROTAX(RMAT(1,1),-RMAT(1,2),-RMAT(1,3),
     1             -RMAT(2,3),RMAT(3,3),AL,BE,GA)
        WRITE (NU,621) GA,BE,AL
        CALL PROTAX(RMAT(1,1),-RMAT(3,1),-RMAT(2,1),
     1             -RMAT(2,3),RMAT(2,2),AL,GA,BE)
        WRITE (NU,622) AL,GA,BE
        RETURN
        END
C
C-------Get answer from keyboard, 1 if answer=default(NDEFLT) or <CR>
C
        FUNCTION IANSW(DEFLT)
C
        CHARACTER YES,DEFLT
        READ (*,10) YES
 10     FORMAT (A1)
        IANSW=0
        IF (YES .EQ. DEFLT) IANSW=1
        IF (YES .EQ. ' ')   IANSW=1
        RETURN
        END
C
C-------Subroutine to set up the rotation matrices CF - RF,SF ---------------
C
        SUBROUTINE ROTMAT (BEAMD,IHKL,GAMMA,RM,M)
        DIMENSION BEAMD(3),IHKL(3),RM(3,3,M),ROT(3,3)
C
        DO 1 I=1,3
 1      RM(1,I,M) = IHKL(I)
C
        DO 2 I=1,3
 2      RM(3,I,M)=BEAMD(I)
C                                            BD x K.line
        RM(2,1,M)=RM(3,2,M)*RM(1,3,M) - RM(3,3,M)*RM(1,2,M)
        RM(2,2,M)=RM(3,3,M)*RM(1,1,M) - RM(3,1,M)*RM(1,3,M)
        RM(2,3,M)=RM(3,1,M)*RM(1,2,M) - RM(3,2,M)*RM(1,1,M)
C                                           (BD x K.line) x BD
        RM(1,1,M)=RM(2,2,M)*RM(3,3,M) - RM(2,3,M)*RM(3,2,M)
        RM(1,2,M)=RM(2,3,M)*RM(3,1,M) - RM(2,1,M)*RM(3,3,M)
        RM(1,3,M)=RM(2,1,M)*RM(3,2,M) - RM(2,2,M)*RM(3,1,M)
C                                  (normal. each row)
        CALL NORMAT(RM,M)
C                                  (rot. to ref fr.)
        CALL MCBA (0.,0.,GAMMA,ROT,1)
        CALL MMPRD (ROT,1,RM,M,RM,M)
C
        RETURN
        END
C
C
C-------MATRIX MULTIPLACTION------------------------------------
C           A X B = R
        SUBROUTINE MMPRD (A,L,B,M,R,N)
        DIMENSION A(3,3,N),B(3,3,N),R(3,3,N),P(3,3)
        DO 10 J=1,3
        DO 10 I=1,3
        P(I,J)= 0.0
        DO 10 K=1,3
  10    P(I,J) = P(I,J) + A(I,K,L)*B(K,J,M)
        DO 20 I=1,3
        DO 20 J=1,3
  20    R(I,J,N) = P(I,J)
        RETURN
        END
CCC
C-------FORMULATION OF MATRIX REPRESENTING A THREE STEP ROTATION--
C
        SUBROUTINE MCBA (X,Y,Z,CBA,N)
        DIMENSION CBA(3,3,N)
C
        A = SIN(X)
        B = COS(X)
        C = SIN(Y)
        D = COS(Y)
        E = SIN(Z)
        F = COS(Z) 
        CBA(1,1,N)=F*D
        CBA(1,2,N)=F*C*A - E*B
        CBA(1,3,N)=F*C*B + A*E
        CBA(2,1,N)=E*D
        CBA(2,2,N)=A*E*C + F*B
        CBA(2,3,N)=B*E*C - A*F
        CBA(3,1,N)=-C
        CBA(3,2,N)=A*D
        CBA(3,3,N)=B*D
        RETURN
        END
CCC
C-------CALCULATION OF ROT.ANGLES AROUND X, Y, AND Z AXIS-----
C
        SUBROUTINE PROTAX(A11,A21,A31,A32,A33,AL,BE,GA)
C
	Real A11,A21,A31,A32,AL,BE,GA,BE1,BE2
        PI=3.141593
        CF=180/PI
        BE1=ASIN(-A31)
        BE2=PI-ABS(BE1)
        IF (BE1 .LT. 0) BE2=-BE2
        D1=COS(BE1)
        D2=COS(BE2)
        CALL ACANG(A11/D1,A21/D1,GA1)
        CALL ACANG(A11/D2,A21/D2,GA2)
        CALL ACANG(A33/D1,A32/D1,AL1)
        CALL ACANG(A33/D2,A32/D2,AL2)
        TEST1=ABS(AL1)+ABS(BE1)+ABS(GA1)
        TEST2=ABS(AL2)+ABS(BE2)+ABS(GA2)
        IF (TEST1 .LT. TEST2) GOTO 1
        AL=AL2*CF
        BE=BE2*CF
        GA=GA2*CF
        RETURN
 1      AL=AL1*CF
        BE=BE1*CF
        GA=GA1*CF
        RETURN
        END
CCC
C-------DETERMINE ANGLE, SIN AND COS GIVEN-------------------
C
        SUBROUTINE ACANG(C,S,A)
	Real S,C,A
        IF (S .GT. 0.) GOTO 1
        IF (C .GT. 0.) GOTO 2
        A=-ACOS(C)
        RETURN
 2      A=ASIN(S)
        RETURN
 1      A=ACOS(C)
        RETURN
        END
CCC
C-------DETRMINATION OF AXIS OF ROTATION---------------------
C
        SUBROUTINE AXIS (AM,M,AX)
        DIMENSION AM(3,3,M),AX(3)
C
        AX(1) = AM(3,2,M) - AM(2,3,M)
        AX(2) = AM(1,3,M) - AM(3,1,M)
        AX(3) = AM(2,1,M) - AM(1,2,M)
        STA = AX(1)**2 + AX(2)**2 + AX(3)**2
        IF (STA .GT. 1E-6) GOTO 2
        DO 1 I=1,3
 1      AX(I) = SQRT(1+AM(I,I,M))
        AX(2) = SIGN(AX(2),AM(2,1,M))
        AX(3) = SIGN(AX(3),AM(3,1,M))
        STA=3+AM(1,1,M)+AM(2,2,M)+AM(3,3,M)
 2      STA=SQRT(STA)
        DO 3 I=1,3
 3      AX(I)=AX(I)/STA
        RETURN
        END
CCC
C-------CALCULATE MATRIX DM FOR THE 24 CUBIC SYMM. OPERATIONS---
C
        SUBROUTINE SYMMAT(DM)
C
        DIMENSION DM(3,3,24)
        R90=3.141593/2
        DO 80 IXY=1,6
        GO TO (61,62,63,64,65,66),IXY
   61   XA=0.
        YA=0.
        GO TO 67
   62   XA=0.
        YA=R90
        GO TO 67
   63   XA=-R90
        YA=0.
        GO TO 67
   64   XA=0.
        YA=-R90
        GO TO 67
   65   XA=R90
        YA=0.
        GO TO 67
   66   XA=0.
        YA=2.*R90
   67   DO 80 IZ=1,4
        ZI=IZ
        ZA=(ZI-1.)*R90
        N=(IXY-1)*4+IZ
   80   CALL MCBA (XA,YA,ZA,DM,N)
        RETURN
        END
CCC
C-------ROUND REAL VALUE (POS. AND NEG.)
C
        FUNCTION IRND (X)
C
        IF (X .GT. 0.) GOTO 1
        IF (X .LT. 0.) GOTO 2
        IRND=0
        RETURN
 1      IRND=IFIX(X+.5)
        RETURN
 2      IRND=IFIX(X-.5)
        RETURN
        END
CCC
C-------INVERT ROT. MATRIX-----------------------------------------
C
        SUBROUTINE INVERT (A,K)
        DIMENSION A(3,3,K)
C
        C = A(1,2,K)
        A(1,2,K) = A(2,1,K)
        A(2,1,K) = C
        C = A(1,3,K)
        A(1,3,K) = A(3,1,K)
        A(3,1,K) = C
        C = A(2,3,K)
        A(2,3,K) = A(3,2,K)
        A(3,2,K) = C
        RETURN
        END
CCC
C-------NORMALIZE MATRIX  (ROW BY ROW)------------------------
C
        SUBROUTINE NORMAT(A,K)
        DIMENSION A(3,3,K)
C
        DO 1 J=1,3
        STP=SQRT(A(J,1,K)**2+A(J,2,K)**2+A(J,3,K)**2)
        DO 1 I=1,3
  1     A(J,I,K)=A(J,I,K)/STP
        RETURN
        END
C**************************************************************
c
	Subroutine Extrot(infile)
c
	  Implicit  None
	  Real                rr(10,3,3)
	  Double Precision    grainNum,matnum,range(10,2)
	  Integer             numex, i1, i2, j, k
	  Integer             status3, Eof3
	  character*64        answer,infile
	  character*7         progsource
	  Real                a,b,c
	  Logical             numfound,doneNow
c
	  parameter   (Eof3 = -1)
c
	  Common /ExtRots/ rr(10,3,3),range(10,2),numex
c
	  write(9,*)
	  write(9,*) ' enter the number of grain number'
	  write(9,*) ' ranges needing relative rotation'
	  read(9,*) numex
	  Do (i1 = 1,numex)
	    status3 = 0
   8      write(9,*)
	    write(9,10) i1
  10	    format(' enter the lowest and highest grain',
     1          ' numbers of range ',i3,'  (LLL.LLL,HHH.HHH)')
	    read(9,12) (range(i1,j),j=1,2)
  12      format(2(f8.3))
  15 	    write(9,*)
          write(9,20) i1
  20      format(' enter the number assigned to the',
     1           ' relative rotation matrix',/' located',
     2           ' in the input file for range ',i3)
	    read(9,22) matnum
  21	    format(' grain # : ',f8.3)
  22	    format(f8.3)
  23	    format(f8.3,7a)
c
	    doneNow  = .false.
	    numfound = .false.
	    open (3,file = infile,status = 'old')
	    Do While (.not.doneNow)
  1	      Read(3,23,iostat=status3) grainNum,progsource
	      write(9,21) grainNum
	      if (grainNum.eq.matnum) then
	        read(3,*,iostat=status3) a,b
	        Do 24 j = 1,3
  24           read(3,*,iostat=status3) (rr(i1,j,k),k=1,3)
  		   numfound = .true.
	      else
	        read(3,*,iostat=status3) a,b
	        read(3,*,iostat=status3) a,b,c
	        read(3,*,iostat=status3) a,b,c
	        read(3,*,iostat=status3) a,b,c
	      End if
	      if (numfound) then
	        doneNow = .true.
	      else if (status3.ne.Eof3) then
	        doneNow = .false.
	      else if (status3.eq.Eof3) then
	        doneNow = .true.
	      end if
	    Repeat
	    close(3)
c
	    if (.not.numfound) then
	      write(9,*)
		write(9,*) ' didnt find that pattern number'
		goto 8
	    End if
	    write(9,*)
	    write(9,25) i1,(range(i1,j),j=1,2)
  25	    format(' range ',i3,' from grain ',f7.3,' to grain ',
     1            f7.3)
	    write(9,*)
	    write(9,30) (rr(i1,1,k),k=1,3)
	    write(9,31) (rr(i1,2,k),k=1,3)
	    write(9,32) (rr(i1,3,k),k=1,3)
  30      format(' relative rotation matrix:'/3(f11.7,' '))
  31      format(3(f11.7,' '))
  32      format(3(f11.7,' '))
	    answer = 'none'
	    write(9,*)
	    write(9,*) ' are numbers entered correctly ? (n)'
	    read(9,*) answer
	    if (answer.eq.'n') goto 8
	  End do
	  Return
	  End
c
c------------------------------------------------------------------
c
	Subroutine MatMult(m1,ne,m2,k,prod)
c
	  Implicit None
	  Integer      i1,i2,i3,ne,k
	  Real         m1(10,3,3),m2(3,3,2),prod(3,3)
c

	  Do 195 i1 = 1,3
	   Do 195 i2 = 1,3
	    prod(i1,i2) = 0.0
  195	  continue
	  write(9,*) 'ne k ',ne,k
	  Do (i1 = 1,3)
	    write(9,*) (m2(i1,i2,k),i2=1,3)
	  end do
	  
	   Do (i1 = 1,3)
	    Do (i2 = 1,3)
            prod(i1,i2) = m1(ne,i1,1)*m2(1,i2,k)
     1                    + m1(ne,i1,2)*m2(2,i2,k)
     2                    + m1(ne,i1,3)*m2(3,i2,k)
	    end do
	   end do

	write(9,*) ' prod ',(prod(1,i1),i1=1,3)
	write(9,*) '      ',(prod(2,i1),i1=1,3)
	write(9,*) '      ',(prod(3,i1),i1=1,3)
	Return
	End
c
c------------------------------------------------------------------