CHARACTER*40 INFILE,FILENAME
	LOGICAL EFLG
	real*8 M1(3,4),M2(3,4),ncs(3,4),cell(3,4),rcell(3,4),pi
c	REAL*8 P1(3),P2(3),T(3)
	character*40 string
	INTEGER*4 MM,NN,I,J

	pi = 4*atan(1.)
	write(6,*)'Enter cell constants a,b,c, alpha, beta, gamma:'
	read(5,*)a,b,c, alpha, beta, gamma
	alpha=alpha*pi/180.
	beta = beta*pi/180.
	gamma=gamma*pi/180.
	
	cell(1,1)=a
	cell(1,2)=b*cos(gamma)
	cell(2,2)=b*sin(gamma)
	cell(1,3)=c*cos(beta)
	cell(2,3)=c*(cos(alpha)-cos(beta)*cos(gamma))/sin(gamma)
	cell(3,3)=sqrt(c*c-(cell(3,1)*cell(3,1)+cell(3,2)*cell(3,2)))


	do 20 i=1,3
	do 20 j=1,3
20	rcell(i,j)=cell(i,j)
	mm=3
	call matinv(rcell,mm,j)

c 	read operator into m1
5	WRITE(6,*) 'ENTER FILENAME rotation-trans op (.o) on orthog Angst:'
	READ(5,1200) INFILE
1200	FORMAT(A)
	IF (INFILE(:5).EQ.'     ') STOP
	  INQUIRE(FILE=INFILE, NAME=FILENAME,EXIST=EFLG)
	  IF (.NOT.EFLG) GOTO 5
        OPEN (UNIT=2,FILE=INFILE,STATUS='OLD',READONLY)
c	READ(2,*) '.SPACE_GROUP_operators    r   12 (3f12.7)'
	READ(2,'(A40)') string
	do 1150 j=1,4
1150	read(2,*)(ncs(I,J),I=1,3)
1151	FORMAT(3f14.7)
	CLOSE(UNIT=2)
	MM=3
	NN=4
	CALL MATPRN (ncs,MM,NN)

c	multiply the three matrices together:   ncs(i,j) is noncryst sym operator; j=4 => translation.
	do 30 i=1,3
	do 30 j=1,3
	m1(i,j)=0
	do 30 k=1,3
30	m1(i,j)=m1(i,j)+ncs(i,k)*cell(k,j)

	do 40 i=1,3
	do 40 j=1,3
	m2(i,j)=0
	do 40 k=1,3
40	m2(i,j)=m2(i,j)+rcell(i,k)*m1(k,j)

	do 50 i=1,3
	m2(i,4)=0
	do 50 k=1,3
50	m2(i,4)=m2(i,4)+rcell(i,k)*ncs(k,4)

	WRITE(6,*)' NCS OPERATOR (fractional):'
	WRITE(6,*)' ROTATE MATRIX'
	DO 1100 I=1,3
1100	WRITE(6,9) (M2(I,J),J=1,3)
	WRITE(6,*)' TRANSLATE '
    	WRITE(6,9) (M2(I,4),I=1,3)
9	FORMAT(' ',3F14.6)
	write(6,*)' this operator will be saved as an O datablock RT.O'
	WRITE(6,*)'  Note this operator is for fractional coordinates'
	WRITE(6,*) 'ENTER NAME FOR output of resulting rot-trans op (.o):'
	READ(5,1200) filename
	open (unit=3,name=FILENAME,status='unknown')
	WRITE(3,*) '.SPACE_GROUP_operators    r   12 (3f14.7)'
	do 1250 j=1,4
1250	write(3,1151)(M2(I,J),I=1,3)
	close (unit=3)

	open (unit=3,name='frac2angs.o',status='unknown')
	WRITE(3,*) '.SPACE_GROUP_operators    r   12 (3f14.7)'
	do 1260 j=1,4
1260	write(3,1151)(cell(I,J),I=1,3)
	close (unit=3)

	open (unit=3,name='angs2frac.o',status='unknown')
	WRITE(3,*) '.SPACE_GROUP_operators    r   12 (3f14.7)'
	do 1270 j=1,4
1270	write(3,1151)(rcell(I,J),I=1,3)
	close (unit=3)

	END

	SUBROUTINE MATPRN(ARRAY,M,N)
	REAL*8 ARRAY(3,4)
	WRITE(6,*)'M, N=',M,N
	DO 10 I=1,M
10	WRITE (6,50) (ARRAY(I,J),J=1,N)
50	format(8f10.5)
  	RETURN
        END


	SUBROUTINE MATINV (ARRAY, NORDER, DET)
	DOUBLE PRECISION ARRAY, AMAX, SAVE
2	DIMENSION ARRAY(3,4), IK (10), JK (10)
10	DET = 1
11 	DO 100  K=1, NORDER
C
C	FIND LARGEST ELEMENT ARRAY(I,J) IN REST OF MATRIX
C
	AMAX=0
21 	DO 30 I=K, NORDER
     	DO 30 J=K, NORDER
23 	IF (DABS (AMAX) - DABS (ARRAY(I,J) ) ) 24, 24, 30
24 	AMAX=ARRAY(I,J)
     	IK(K) =I
     	JK(K) =J
30 	CONTINUE
C
C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K,K)
C
31 	IF (AMAX) 41, 32, 41
32 	DET=0.
     	GO TO 140
41 	I=IK (K)
     	IF (I-K) 21, 51, 43
43 	DO 50 J=1, NORDER
     	SAVE= ARRAY(K,J)
      	ARRAY(K,J)= ARRAY(I,J)
50 	ARRAY(I,J)= -SAVE
51 	J=JK (K)
     	IF (J-K) 21, 61, 53
53 	DO 60 I=1, NORDER
     	SAVE= ARRAY(I,K)
     	ARRAY(I,K)= ARRAY(I,J)
60 	ARRAY(I,J)= -SAVE
C
C	ACCUMULATE ELEMENTS OF INVERSE MATRIX
C
61 	DO 70 I=1, NORDER
     	IF(I-K) 63, 70, 63
63 	ARRAY(I,K) = -ARRAY(I,K)  / AMAX
70 	CONTINUE
71 	DO 80 I=1, NORDER
     	DO 80 J=1, NORDER
     	IF (I-K) 74, 80, 74
74 	IF (J-K) 75, 80, 75
75 	ARRAY(I,J) = ARRAY(I,J) + ARRAY(I,K) *ARRAY(K,J)
80 	CONTINUE
81 	DO 90 J=1, NORDER
     	IF (J-K) 83, 90, 83
83 	ARRAY(K,J) = ARRAY(K,J) / AMAX
90 	CONTINUE
     	ARRAY(K,K) = 1. / AMAX
100 	DET=DET *AMAX
C
C 	RESTORE ORDERING OF MATRIX
C 
101 	DO 130 L=1, NORDER
       	K= NORDER  - L +1
       	J= IK (K)
       	IF (J-K) 111, 111, 105
105 	DO 110  I=1, NORDER
       	SAVE= ARRAY(I,K)
        ARRAY(I,K)= -ARRAY(I,J)
110 	ARRAY(I,J)= SAVE
111 	I= JK(K)
       	IF (I-K) 130, 130, 113
113 	DO 120  J=1, NORDER
       	SAVE= ARRAY(K,J)
       	ARRAY(K,J) = -ARRAY(I,J)
120 	ARRAY(I,J)= SAVE
130 	CONTINUE
140  	RETURN
	END