	subroutine lsq(iprnt)
c
c Routine to do least-squares refinement of List Data
c
	implicit none
	include 'single.inc'
	real*4 phi(3),h6(3,3),h7(3,3),z2(3,3),z3(3,3)
     1,da(6),db(6),temp(3,3),tempa(4,listlen),det,sum
     2,tthc,p1,vol,sigcell
	integer*4 tempi(3,listlen),i,j,nlsq,k,iprnt,iz
	common /lsqcom/sigcell(6)
c
c Zero Arrays
c
	do 110 i=1,3
	  do 100 j=1,3
	    h6(j,i)=0.0
	    h7(j,i)=0.0
100	  continue
110	continue
c
c transfer angles and indices from list to temp arrays for lsq
c
	nlsq=0
	do 50 i=1,nlst
		if(iuse(i) .eq. 0)goto50
	        nlsq=nlsq+1
		do j=1,3
		  tempa(j,nlsq)=anglst(j,i)
		  tempi(j,nlsq)=INT(ihlst(j,i)+SIGN(0.5,IHLST(J,I)))
		enddo
		tempa(4,nlsq)=anglst(4,i)
50	continue
	if(nlsq.lt.3)then
	    type 1,nlsq
1	    format(' Only ',i3,' Reflection(s) with USE Flag Set,'
     1		,' 3 are Needed for LSQ')
 	    return
	endif
c
c form matrices H(phi) and H of rank 3 x N and do matrix products
c
c H6 = H(phi)*H;  H7 = H(T)*H
c
	do 140 i=1,nlsq
	  call calch(tempa(1,i),phi)
	  do 130 j=1,3
	    do 120 k=1,3
	      h6(k,j)=h6(k,j)+phi(j)*tempi(k,i)
	      h7(k,j)=h7(k,j)+tempi(j,i)*tempi(k,i)
120	    continue
130	  continue
140	continue
	call minv(h7,h7,det)
	if(abs(det).gt.1.0e-4)go to 150
	  type 2
2	  format('0Singular Matrix')
	  return
150	call mattsp(h6,h6)
	call matmul(h6,h7,ub)
	if(nlsq.lt.4)go to 230
	if(iprnt.eq.0)type 4
4	format('0  H  K  L    TTHO    TTHC         Vector Deltas'/)
	do 190 i=1,3
	  h6(i,1)=0.0
190	continue
	do 220 iz=1,nlsq
	  call calch(tempa(1,iz),z2)
	  sum=0.0
	  do 210 i=1,3
	    z3(i,1)=0.0
	    do 200 j=1,3
	      z3(i,1)=z3(i,1)+ub(i,j)*tempi(j,iz)
200	    continue
	    sum=sum+z3(i,1)**2
	    z3(i,1)=z2(i,1)-z3(i,1)
	    h6(i,1)=h6(i,1)+z3(i,1)**2
210	  continue
	  tthc=asin(sqrt(sum)*wavel*0.5)*114.5916
	  if(iprnt.eq.0)then
	    type 6,(tempi(j,iz),j=1,3),tempa(1,iz),tthc,(z3(j,1),j=1,3)
	  endif
6	  format(1x,3i3,2f8.2,2x,3f9.5)
220	continue
230	if(iprnt.eq.0)type 8,((ub(i,j),j=1,3),i=1,3)
8	format('0Orientation Matrix'//(3f12.5))
	call omtob(iprnt)
	call update
	if(nlsq.lt.4)return
	do 250 i=1,3
	  do 240 j=1,3
	    z2(j,i)=sqrt(h7(i,i)/(nlsq-3)*h6(j,1))
240	  continue
250	continue
	do 260 i=1,6
	  db(i)=0.0
260	continue
	p1=0.0
	do 300 i=1,3
	  do 290 j=1,3
	    ub(j,i)=ub(j,i)+z2(j,i)        !add a bit to an element of UB
	    call mattsp(ub,temp)
	    call matmul(temp,ub,z3)        !new value of G* in Z3
	    call minv(z3,z3,det)
c calculate direct space volume and lattice constants
	    vol=1.0/sqrt(det)
	    do 270 k=1,3
	      da(k)=sqrt(z3(k,k))
270	    continue
	    da(4)=acos(z3(3,2)/(da(3)*da(2)))*57.2958
	    da(5)=acos(z3(3,1)/(da(3)*da(1)))*57.2958
	    da(6)=acos(z3(1,2)/(da(1)*da(2)))*57.2958
	    do 280 k=1,6
	      db(k)=db(k)+(da(k)-cell(k))**2
280	    continue
	    p1=p1+(vol-cell(7))**2
c restore UB element
	    ub(j,i)=ub(j,i)-z2(j,i)
290	  continue
300	continue
	do i=1,6
	  db(i)=sqrt(db(i))
	  sigcell(i)=db(i)
	enddo
	do i=4,6
	   sigcell(i)=sin(cell(i)/57.2958)*sigcell(i)/57.2958
	enddo
	p1=sqrt(p1)
	if(iprnt.eq.0)type 10,db,p1
10	format(' ESDs',12x,3f8.4,3f8.3,f9.3)
	return
	end
	subroutine lsvec
c    reads reflection data from reflection list and performs
c    a constrained least-squares refinement of the
c    unit cell parameters, writing full results to the screen.
c
c    coded 6-May-81 by Russell Ralph and Larry Finger
c			Geophysical Laboratory, Washington DC
c
c	method from Shoemaker and Bassi (1970), Acta Cryst. A26, 97.
c
c Modified for Diffractometer Control System by LWF,  6-Feb-92
c
	include 'single.inc'
	dimension a(3,3),yb(3,3),s(3,3),ss(3),vect(3,100),aa1(6),as1(6)
	1 ,bb(3,3),u(3,3),v(3,3),v0(3,3),v0inv(3,3),d(3,3),angl(4)
	2 ,anglc(4),xx(12),hh(3,100),wk3(12),atemp(7,9),
	3 dcov(9,9),derivs(7,9),sigs(7),tmp1(3),tmp2(3)
	integer h(3,100)
	common /clsq/x(3,3),a1(9),a2(6),vol1,vol2,xa(15,16),xb(12)
C
	DATA RAD/57.2958/,FFACT/5.E-4/
C
C  FORMAT STATEMENTS
C
10	format('0Crystal Type: 0-Triclinic, 1-Monoclinic (B-unique)'
	1,', 2-Monoclinic (C-unique),'/'$3-Orthorhombic, 4-Tetragonal,'
	2,' 5-Hexagonal, 6-Cubic? ')
12	FORMAT(//'0Results for Constrained Refinement')
16	FORMAT('0Cell   ',3F10.5,3F8.3,F10.3)
19	FORMAT(' Sigmas ',3F10.5,3F8.3,F10.3)
24	FORMAT('0',T17,'OBS',T25,'CAL',T40,'Vector Deltas'
	1/'   h   k   l',T17,'2TH',T25,'2TH',
	2 t32,'Del X',T42,'Del Y',T52,'Del Z',T64,'Length'/)
26	FORMAT(3I4,2F8.3,3F10.5,2X,F10.5)
C  INITIALIZE COUNTERS AND BEGIN UNCONSTRAINED PROBLEM
	NUM=0
	DO 180 J=1,3
	DO 180 I=1,3
	A(I,J)=0.
	YB(I,J)=0.
180	X(I,J)=0.
	DO 190 J=1,12
	DO 185 I=1,12
185	XA(I,J)=0.
190	XB(J)=0.
c
c transfer angles and indices from list to temp arrays for lsq
c
	num=0
	do i=1,nlst
		if(iuse(i) .ne. 0) then
	          num=num+1
		  do j=1,3
		    h(j,num)=ANINT(ihlst(j,i))
		  enddo
		  call calch(anglst(1,i),vect(1,num))
		endif
	enddo
	if(num.lt.3) then
	  type *,'Too Few Reflections for Least Squares'
	  return
	endif
c
C  Do initial unconstrained lsq:
c
c   form normal equations matrices a, yb, and xa
c
	do k=1,num
	  DO I=1,3
	    DO J=1,3
	      A(J,I)=A(J,I)+H(I,k)*H(J,k)
	      YB(J,I)=YB(J,I)+VECT(I,k)*H(J,k)
	    enddo
	  enddo
	enddo
	DO J=1,3
	    DO I=1,3
		XA(I,J)=A(I,J)
		XA(I+3,J+3)=A(I,J)
		XA(I+6,J+6)=A(I,J)
	    enddo
	enddo
C  INVERT MATRIX A
	CALL MXLNEQ(A,3,3,DET,JRANK,1.0E-6,WK3,0)
	IF(JRANK.NE.3)then
	  type *,'Matrix is Singular'
	  return
	endif
C  SOLVE FOR ORIENTATION MATRIX X
	DO I=1,3
	    DO J=1,3
		DO K=1,3
		    X(I,J)=X(I,J)+A(J,K)*YB(K,I)
		enddo
	    enddo
	    SS(I)=0.0
	enddo
c
c Get Cell Parameters
	call params
C
C  BEGIN CONSTRAINED LEAST-SQUARES
C
c
c Get Crystal system
c
199	type 10
	read(5,'(i5)',end=9999,err=199)icrstp
	if(icrstp.eq.0)return
C
c  construct the true orientation matrix bb corresponding
c  to the unconstrained cell parameters
c
C
c  construct the true orientation matrix bb corresponding
c  to the unconstrained cell parameters
	WRITE(5,12)
c
C  FORM AND INVERT BB AND OBTAIN THE TRUE ORTHOGONAL ROTATION U FROM X
c
	CALL BCALC(BB,A1,A2)
	CALL MINV(BB,BB,DET)
	call matmul(x,bb,u)
c  constrain reciprocal lattice to have proper symmetry
c alpha never refined in cases treated
	a2(4)=0.0
	GO TO (740,750,760,770,780,790),ICRSTP
C MONOCLINIC - B UNIQUE
740	A2(6)=0.0
	GO TO 800
C MONOCLINIC - C UNIQUE
750	A2(5)=0.0
	GO TO 800
C ORTHORHOMBIC
760	A2(5)=0.0
	A2(6)=0.0
	GO TO 800
C TETRAGONAL
770	A2(1)=SQRT(A2(1)*A2(2))
	A2(2)=A2(1)
	A2(5)=0.0
	A2(6)=0.0
	GO TO 800
C HEXAGONAL
780	A2(1)=SQRT(A2(1)*A2(2))
	A2(2)=A2(1)
	A2(5)=0.0
	A2(6)=0.5
	GO TO 800
C CUBIC
790	A2(1)=(A2(1)*A2(2)*A2(3))**0.33333333
	A2(2)=A2(1)
	A2(3)=A2(1)
	A2(5)=0.0
	A2(6)=0.0
C CONSTRUCT CONSTRAINED B
800	CALL RECIP(A2,VOL2,A1,VOL1)
	CALL BCALC(BB,A1,A2)
C  CONSTRUCT CONSTRAINED GUESS V0 FOR ORIENTATION MATRIX
	call matmul(u,bb,v0)
C  INVERT V0
	CALL MINV(V0,V0INV,DET)
	IF(abs(det).le.0.00001)then
	  type *,'Normal Equations Matrix is Singular.'
	  return
	endif
C  CALCULATE "OBSERVED" INDICES HH
	DO I=1,NUM
	  DO J=1,3
	    HH(J,I)=0.
	    DO K=1,3
	      HH(J,I)=HH(J,I)+V0INV(J,K)*VECT(K,I)
	    enddo
	  enddo
	enddo
C  SET UP EXPANDED VECTOR XB
	DO K=1,NUM
	  DO I=1,3
	    XB(I)=XB(I)+H(I,K)*HH(1,K)
	    XB(I+3)=XB(I+3)+H(I,K)*HH(2,K)
	    XB(I+6)=XB(I+6)+H(I,K)*HH(3,K)
	  enddo
	enddo
C TREAT CONSTRAINED PROBLEM
	CALL CONSTR(ICRSTP,IERR,NPARM)
	IF(ierr.ne.0)then
	  type *,'Normal Equations Matrix is Singular.'
	  return
	endif
	DO I=1,3
	  D(1,I)=XB(I)
	  D(2,I)=XB(I+3)
	  D(3,I)=XB(I+6)
	enddo
C  EXTRACT ORIENTATION MATRIX V AND CELL PARAMETERS FROM D
	call matmul(v0,d,v)
	DO I=1,3
	  DO J=1,3
	    X(I,J)=V(I,J)
	  enddo
	enddo
	CALL PARAMS
	WRITE(5,24)
C  BACK CALCULATE "OBSERVED" INDICES AND EVALUATE VARIANCE-
C    -COVARIANCE MATRIX DCOV FOR D
	FACT=0.
	DO K=1,NUM
	    DO I=1,3
		TMP1(I)=0.0
		TMP2(I)=-HH(I,K)
		DO JJ=1,3
		    TMP1(I)=TMP1(I)+X(I,JJ)*H(JJ,K)
		    TMP2(I)=TMP2(I)+D(I,JJ)*H(JJ,K)
		enddo
	    enddo
	    CALL calca(vect(1,k),angl,ierr)
	    CALL calca(tmp1,anglc,ierr)
	    TEMP=0.0
	    DO JJ=1,3
		TMP1(JJ)=VECT(JJ,K)-TMP1(JJ)
		TEMP=TEMP+TMP1(JJ)**2
	    enddo
	    FACT=FACT+TMP2(1)**2+TMP2(2)**2+TMP2(3)**2
	    TEMP=SQRT(TEMP)
	    WRITE(5,26)(H(JJ,K),JJ=1,3),ANGL(1),ANGLC(1),TMP1,TEMP
	enddo
	TEMP=FACT
	FACT=FACT/(3*NUM-NPARM)
	DO J=1,9
	    DO I=1,9
		DCOV(I,J)=FACT*XA(I,J)
	    enddo
	enddo
c
	WRITE(5,16) A1(1),A1(2),A1(3),A1(7),A1(8),A1(9),VOL1
c
C  PERTURB D ONE ELEMENT AT A TIME, PROPAGATE RESULTS TO
C    ORIENTATION MATRIX V, AND ACCUMULATE PARTIAL DERIVATIVES
C    FOR EACH D.L. PARAMETER IN MATRIX DERIVS
	DO 1505 K=1,6
	AA1(K)=A1(K)
1505	CONTINUE
	VOLL1=VOL1
	DO 1690 J=1,3
	DO 1690 I=1,3
	D(I,J)=D(I,J)+FFACT
	DO 1510 JJ=1,3
	DO 1510 II=1,3
	X(II,JJ)=0.
	DO 1510 K=1,3
1510	X(II,JJ)=X(II,JJ)+V0(II,K)*D(K,JJ)
	CALL PARAMS
	DO 1520 K=1,6
1520	AS1(K)=0.
	L=J+3*(I-1)
	DO 1530 K=1,6
1530	DERIVS(K,L)=(A1(K)-AA1(K))/FFACT
	DERIVS(7,L)=(VOL1-VOLL1)/FFACT
	D(I,J)=D(I,J)-FFACT
1690	CONTINUE
C  CALCULATE AND PRINT STANDARD DEVIATIONS
	DO 1700 K=1,7
	DO 1700 I=1,9
	ATEMP(K,I)=0.0
	DO 1700 J=1,9
1700	ATEMP(K,I)=ATEMP(K,I)+DERIVS(K,J)*DCOV(J,I)
	DO 1710 K=1,7
	DO 1710 I=1,7
	DCOV(K,I)=0.0
	DO 1710 J=1,9
1710	DCOV(K,I)=DCOV(K,I)+ATEMP(K,J)*DERIVS(I,J)
C CLEAR NON-REFINED ROWS AND COLUMNS OF DCOV
	DO 1715 I=1,6
	DCOV(I,4)=0.0
	DCOV(4,I)=0.0
	IF(ICRSTP.EQ.2)GO TO 1712
	DCOV(6,I)=0.0
	DCOV(I,6)=0.0
1712	IF(ICRSTP.EQ.1)GO TO 1715
	DCOV(I,5)=0.0
	DCOV(5,I)=0.0
1715	CONTINUE
	DO 1720 K=1,7
	FACT=DCOV(K,K)
	IF(FACT.LE.0.0)FACT=1.0E-24
1720	SIGS(K)=SQRT(FACT)
	DO 1725 J=1,7
	DO 1725 K=1,7
1725	ATEMP(J,K)=DCOV(J,K)/(SIGS(J)*SIGS(K))
	DO 1730 K=4,6
1730	SIGS(K)=RAD*SIGS(K)*SQRT(1.0-AA1(K)**2)
	WRITE(5,19) SIGS
	I=3*NUM
c	WRITE(5,28)I,NPARM,TEMP
c	WRITE(5,20)((DCOV(I,J),I=J,6),J=1,6)
c	WRITE(5,22)((ATEMP(I,J),I=J,6),J=1,6)
C  TERMINATE
9999	return
	END
	SUBROUTINE BCALC(B,A,ASTAR)
C CALCULATES MATRIX B FROM REAL AND RECIPROCAL LATTICES
	DIMENSION B(3,3),A(6),ASTAR(6)
	B(1,1)=ASTAR(1)
	B(1,2)=ASTAR(2)*ASTAR(6)
	B(1,3)=ASTAR(3)*ASTAR(5)
	B(2,1)=0.
	B(2,2)=ASTAR(2)*SQRT(1.0-ASTAR(6)*ASTAR(6))
	B(2,3)=-ASTAR(3)*SQRT(1.0-ASTAR(5)*ASTAR(5))*A(4)
	B(3,1)=0.
	B(3,2)=0.
	B(3,3)=1.0/A(3)
	RETURN
	END
	SUBROUTINE CONSTR(CRYSTP,IERR,NPARM)
C
C SUBROUTINE TO APPLY CONSTRAINTS FOR CRYSTAL OF 'CRYSTP' TO NORMAL
C  EQUATIONS MATRIX XA, SOLVE SYSTEM AND RETURN XA INVERSE IN
C  XA, SOLUTION IN XB. METHOD OF LAGRANGIAN MULTIPLIERS USED
C  FOR CONSTRAINTS. IERR = 0 IF NO ERROR, NON-ZERO MEANS ERROR.
C
	DIMENSION NTEMP(14)
	common /clsq/x(3,3),a1(9),a2(6),vol1,vol2,xa(15,16),xb(12)
	INTEGER CRYSTP,NCNST(6),CNSTYP(5,6),NCOEF(10)
	1,ROWTP(3,10),COEFTP(3,10)
	REAL COEF(15)
C
C CRYSTP = CRYSTAL TYPE: 0-TRICLINIC, 1-MONOCLINIC(B-UNIQUE)
C    2-MONOCLINIC (C-UNIQUE), 3-ORTHORHOMBIC, 4-TETRAGONAL
C    5-HEXAGONAL, 6-CUBIC
C
C NCNST(J) = NUMBER OF CONSTRAINTS FOR CRYSTAL TYPE J
C
C CNSTYP(J,I) = CONSTRAINT TYPE FOR J'TH CONSTRAINT OF
C     I'TH CRYSTAL TYPE
C
C NCOEF(J) = NUMBER OF COEFFICIENTS IN CONSTRAINT EQUATION I
C
C ROWTP(J,I) = ROW (COLUMN) NUMBER IN NORMAL EQUATIONS MATRIX
C     OF J'TH PARAMETER IN I'TH CONSTRAINT EQUATION
C
C COEFTP(J,I) = COEFFICIENT TYPE OF J'TH TERM IN I'TH
C     CONSTRAINT EQUATION
C
	DATA NCNST/2,2,3,4,4,5/
	DATA CNSTYP/9,10,3*0,7,8,3*0,1,2,3,2*0,4,1,2,3,0
	1,6,1,7,8,0,4,5,1,2,3/
	DATA NCOEF/2,2,2,2,2,3,3,3,3,3/
	DATA ROWTP/2,4,0,3,7,0,6,8,0,1,5,0,1,9,0
	1,1,4,5,7,3,6,7,8,6,6,8,2,6,4,2/
	DATA COEFTP/1,2,15,3,4,15,5,6,15,7,8,15,7,8,15
	1,7,7,8,4,3,9,10,6,11,5,6,12,13,2,14/
C
C EVALUATE MULTIPLIER COEFFICIENTS
C
	COEF(1)=A2(1)/A2(2)
	COEF(2)=A2(2)/A2(1)
	COEF(3)=A2(1)/A2(3)
	COEF(4)=A2(3)/A2(1)
	COEF(5)=A2(2)/A2(3)
	COEF(6)=A2(3)/A2(2)
	COEF(7)=1.0
	COEF(8)=-1.0
	COEF(9)=A2(2)*A2(6)/A2(3)
	COEF(10)=-A2(3)*A2(6)/A2(1)
	COEF(11)=A2(2)*(1.0-A2(6)*A2(6))/A2(3)
	COEF(12)=A2(1)*A2(5)/A2(2)
	COEF(13)=-A2(2)*A2(5)/A2(3)
	COEF(14)=A2(1)*(1.0-A2(5)*A2(5))/A2(2)
	COEF(15)=0.0
C
C CLEAR CONSTRAINT PORTION OF MATRIX XA
C
	DO 100 J=10,16
	DO 100 I=1,9
100	XA(I,J)=0.0
	DO 110 J=1,16
	DO 110 I=10,15
110	XA(I,J)=0.0
	IF(CRYSTP.LT.1.OR.CRYSTP.GT.6)GO TO 210
C
C GET NUMBER OF CONSTRAINT EQUATIONS
C
	NEQ=NCNST(CRYSTP)
	NPARM=9-NEQ
C
C LOOP THROUGH CONSTRAINTS
C
	DO 130 I=1,NEQ
C
C GET CONSTRAINT TYPE AND NUMBER OF COEFFICIENTS
C
	ITYP=CNSTYP(I,CRYSTP)
	NUMB=NCOEF(ITYP)
C
C LOOP THROUGH TERMS IN CONSTRAINT EQUATION
C
	DO 120 J=1,NUMB
C
C GET ROW AND COLUMN OF TERM TO MODIFY AND COEFFICIENT OF MODIFICATION
C
	IROW=ROWTP(J,ITYP)
	NCF=COEFTP(J,ITYP)
	XA(IROW,9+I)=COEF(NCF)
	XA(9+I,IROW)=COEF(NCF)
120	CONTINUE
130	CONTINUE
	NRANK=9+NEQ
C
C COPY VECTOR TO AUGMENTED POSITION, INVERT MATRIX AND SOLVE EQUATIONS
C
	DO 140 I=1,9
140	XA(I,NRANK+1)=XB(I)
	CALL MXLNEQ(XA,NRANK,15,DET,JRANK,1.0E-6,NTEMP,1)
	IF(JRANK.NE.NRANK)GO TO 200
	DO 150 I=1,9
150	XB(I)=XA(I,NRANK+1)
	IERR=0
	RETURN
C
C MATRIX SINGULAR
C
200	IERR=1
	RETURN
C
C CRYSTAL TYPE NOT PROPER
C
210	IERR=2
	RETURN
	END
      SUBROUTINE MXLNEQ(A,NRANK,IDA,DET,JRANK,EPS,NTEMP,NCV)
C
C MATRIX INVERSION AND EQUATIONS SOLUTION ROUTINE FROM R. L. HOTCHKISS
C    UNIVERSITY OF MINNESOTA.
C
C ROUTINE USES METHOD OF GAUSS-JORDAN DOUBLE PIVOTING.
C
C A IS MATRIX OF LINEAR SYSTEM WITH COEFFICIENTS AUGMENTED BY VECTOR
C NRANK IS INPUT ORDER OF SYSTEM
C IDA IS FIRST DIMENSION OF A IN CALLING PROGRAM
C DET IS OUTPUT VALUE OF DETERMINANT
C JRANK IS OUTPUT RANK OF SYSTEM
C EPS IS TEST VALUE USED TO DETERMINE SINGULARITY
C NTEMP IS INTEGER ARRAY OF SIZE NRANK
C NCV IS NUMBER OF COLUMNS AUGMENTED. IF NCV=0, INVERSION ONLY.
C
	DIMENSION A(2),NTEMP(2)
C	 INITIALIZE
	NC=NRANK+NCV
	DET=1.0
C	 START MASTER REDUCTION	LOOP
	DO 240 K=1,NRANK
	PIV=0.0
	KSUB=IDA*(K-1)
C	 START SEARCH FOR PIVOT
	DO 110 J=K,NRANK
	JSUB=IDA*(J-1)
	DO 110 I=K,NRANK
	IJ=JSUB+I
C	 CHECK IF ABS(A(I,J)).GT.ABS(PIVOT)
	P=ABS(A(IJ))
	IF(P-PIV)110,110,100
C	 FOUND NEW VALUE FOR TRIAL PIVOT
  100 PIV=P
	PIVOT=A(IJ)
	M=J
	L=I
  110 CONTINUE
C	 CHECK SIZE OF PIVOT
	IF(PIV-EPS)120,130,130
C	 SINGULARITY RETURN
  120 JRANK=K-1
	RETURN
C	 VALID PIVOT FOUND, UPDATE DETERMINANT
  130 DET=DET*PIVOT
C	 SAVE ORIGINAL ROW AND COLUMN OF PIVOT
	NTEMP(K)=256*L+M
	IF(L-K)140,160,140
C	 NEED TO PUT PIVOT ELEMENT INTO	A(K,K),	INTERCHANGE ROWS
  140 DET=-DET
	DO 150 J=1,NC
	JSUB=IDA*(J-1)
	LJ=L+JSUB
	KJ=K+JSUB
	TEMP=A(LJ)
	A(LJ)=A(KJ)
  150 A(KJ)=TEMP
  160 IF(M-K)170,190,170
C	 INTERCHANGE COLUMNS
  170 DET=-DET
	JM=IDA*(M-1)
	JK=KSUB
	DO 180 J=1,NRANK
	JM=JM+1
	JK=JK+1
	TEMP=A(JM)
	A(JM)=A(JK)
  180 A(JK)=TEMP
C	 REDUCE	PIVOT ROW
  190 PIVOT=1.0/PIVOT
	DO 200 J=1,NC
	KJ=IDA*(J-1)+K
  200 A(KJ)=PIVOT*A(KJ)
	KK=K+KSUB
	A(KK)=0.0
C	 REDUCE	NON-PIVOT ROWS
	DO 230 I=1,NRANK
	IK=I+KSUB
	TEMP=A(IK)
	IF(I-K)210,230,210
  210 DO 220 J=1,NC
	IJ=IDA*(J-1)+I
	KJ=IDA*(J-1)+K
  220 A(IJ)=A(IJ)-TEMP*A(KJ)
	A(IK)=-TEMP*PIVOT
  230 CONTINUE
  240 A(KK)=PIVOT
C	 INVERSION COMPLETE, RESTORE ORIGINAL ROW AND COLUMN ORDER
	JRANK=NRANK
	K=NRANK
	DO 300 J=1,NRANK
	KSUB=IDA*(K-1)
	M=NTEMP(K)/256
	L=MOD(NTEMP(K),256)
	IF(M-K)250,270,250
C	 INTERCHANGE ROWS
  250 MSUB=IDA*(M-1)
	DO 260 I=1,NRANK
	IK=KSUB+I
	IM=MSUB+I
	TEMP=A(IK)
	A(IK)=A(IM)
  260 A(IM)=TEMP
  270 IF(L-K)280,300,280
C	 INTERCHANGE COLUMNS
  280 DO 290 I=1,NC
	KI=IDA*(I-1)+K
	LI=IDA*(I-1)+L
	TEMP=A(KI)
	A(KI)=A(LI)
  290 A(LI)=TEMP
  300 K=K-1
C	 JOB COMPLETE, RETURN TO CALLING PROGRAM
	RETURN
	END
	SUBROUTINE PARAMS
C
C  EXTRACTS DIRECT LATTICE PARAMETERS A1 AND VOL1 AND RECIPROCAL
C    LATTICE PARAMETERS A2 AND VOL2 FROM ORIENTATION MATRIX X
C
	DIMENSION XTX(3,3)
	common /clsq/x(3,3),a1(9),a2(6),vol1,vol2,xa(15,16),xb(12)
	DATA RAD/57.2958/
C  WORK RECIPROCAL LATTICE FROM METRIC
	DO 50 I=1,3
	DO 40 J=1,3
	XTX(I,J)=0.
	DO 40 K=1,3
	XTX(I,J)=XTX(I,J)+X(K,I)*X(K,J)
40	CONTINUE
	A2(I)=SQRT(XTX(I,I))
50	CONTINUE
	A2(4)=XTX(2,3)/(A2(2)*A2(3))
	A2(5)=XTX(3,1)/(A2(3)*A2(1))
	A2(6)=XTX(1,2)/(A2(1)*A2(2))
	CALL RECIP(A2,VOL2,A1,VOL1)
	DO 60 I=7,9
60	A1(I)=RAD*ATAN2(SQRT(1.0-A1(I-3)**2),A1(I-3))
	RETURN
	END
	SUBROUTINE RECIP(A2,VOL2,A1,VOL1)
C
C ROUTINE TO ACCEPT CELL PARAMETERS AND COSINES IN A2, CALCULATE VOLUME
C  THEN CALCULATE RECIPROCAL CELL AND VOLUME
	DIMENSION A1(6),A2(6),SINES(3)
C
	VOL2=1.-A2(4)*A2(4)-A2(5)*A2(5)-A2(6)*A2(6)+2.*A2(4)*A2(5)*A2(6)
	VOL2=A2(1)*A2(2)*A2(3)*SQRT(VOL2)
C  WORK DIRECT LATTICE FROM RECIPROCAL LATTICE
	VOL1=1./VOL2
	DO 60 I=1,3
	SINES(I)=SQRT(1.-A2(I+3)*A2(I+3))
60	CONTINUE
	A1(1)=A2(2)*A2(3)*SINES(1)/VOL2
	A1(2)=A2(3)*A2(1)*SINES(2)/VOL2
	A1(3)=A2(1)*A2(2)*SINES(3)/VOL2
	A1(4)=(A2(5)*A2(6)-A2(4))/(SINES(2)*SINES(3))
	A1(5)=(A2(6)*A2(4)-A2(5))/(SINES(3)*SINES(1))
	A1(6)=(A2(4)*A2(5)-A2(6))/(SINES(1)*SINES(2))
	RETURN
	END
