	SUBROUTINE UPDATE
C
C ROUTINE TO UPDATE STORED FILE
C
	include 'single.inc'
	include 'datacol.inc'
	include 'symmetry.inc'
	include 'alph.inc'
	include 'device.inc'
C
	REWIND 3
	write(3,1)title
1	format(60a1)
	ibad=0
	if(badcrystal)ibad=1
	WRITE(3,8)UB,WAVEL,SCN1,SCN2,MINCT,MAXMIN,TMAX,IBAD,UBO
8	FORMAT(3(3E14.7/)6F12.6,I2,3(/3E14.7))
	DO 110 I=NLST+1,LISTLEN
	  DO 100 J=1,3
	    IHLST(J,I)=0.
	    ANGLST(J,I)=0.0
100	  CONTINUE
	  ANGLST(4,I)=0.0
	  IUSE(I)=0
	  IREF(I)=0
110	CONTINUE
120	WRITE(3,10)NLST,((IHLST(J,I),J=1,3),I=1,LISTLEN)
10	FORMAT(I4,10F7.2/(4X,10F7.2))
	WRITE(3,12)((ANGLST(J,I),J=1,4),I=1,LISTLEN)
12	FORMAT(8F8.3)
	WRITE(3,11)(LSTCOUNT(I),I=1,LISTLEN)
11	FORMAT(8I9)
	WRITE(3,13)(IREF(I),IUSE(I),I=1,LISTLEN)
13	FORMAT(40I2)
	WRITE(3,14)FILeNAMe
14	FORMAT(A80)
	write(3,15)nlat,nstd,stdtime,ihkstd
15	FORMAT(20I4)
	write(3,20)spg
20	format(20a1)
	write(3,22)hklr,iorder
22	format(9i3)
	write(3,23)scan
23	format(10f8.4)
	write(3,24)shell
24	format(10f8.3/10f8.3)
c device constraints
	write(3,26)device_no,fixed_phi,chi_flag,phi_flag
     1  ,dc_psimax,chimn,chimx,phimn,phimx
26	format(4i2,5f10.3)
	REWIND 3
	RETURN
	END
	SUBROUTINE CALCH(ANGL,THPHI)
	include 'single.inc'
C
C CALCULATE VECTOR COMPONENTS FROM ANGLES - NOTE: OMEGA MUST BE DEVIATION
C	FROM BISECTING POSITION
C
	DIMENSION ANGL(4),THPHI(3)
	DATA RAD/57.2958/
C GET SINES AND COSINES OF EULERIAN ANGLES
c	SINO=SIN(OMEGA/RAD)*PARITY(2)
	SINO=SIN(angl(2)/RAD)*PARITY(2)
	SINX=SIN(ANGL(3)/RAD)*PARITY(3)
	SINP=SIN(ANGL(4)/RAD)*PARITY(4)
c	COSO=COS(OMEGA/RAD)
	COSO=COS(angl(2)/RAD)
	COSX=COS(ANGL(3)/RAD)
	COSP=COS(ANGL(4)/RAD)
C CALCULATE SIN THETA
	SINT=SIN(ANGL(1)*0.5/RAD)
C FORM TERMS OF H-PHI
	THPHI(1)=2.0*SINT*(COSO*COSX*COSP-SINO*SINP)/WAVEL
	THPHI(2)=2.0*SINT*(COSO*COSX*SINP+SINO*COSP)/WAVEL
	THPHI(3)=2.0*SINT*(COSO*SINX)/WAVEL
	RETURN
	END
	SUBROUTINE CALCA(THPHI,ANGL,IERR)
C
C CALCULATE ANGLES FROM VECTOR - IERR NON-ZERO IF NOT ACCESSIBLE
C
	DIMENSION THPHI(3),ANGL(4)
	include 'single.inc'
	include 'device.inc'
	data rad/57.2958/
c calculate sin theta
	ierr=0
	dinv=sqrt(thphi(1)**2+thphi(2)**2+thphi(3)**2)
	stheta=dinv*wavel*0.5
	if(dinv.le.1.0e-6.or.stheta.gt.1.0)then
c sin(theta) .le. 0 or sin(theta) .gt. 1.0 is an error
	  ierr=1
	  return
	endif
c calculate two theta
	angl(1)=2.0*rad*asin(stheta)
	if(fixed_phi.eq.0)then
c calculate peak position in bisecting mode
	  angl(2)=0.0
	  angl(3)=atan2(thphi(3),sqrt(thphi(1)**2+thphi(2)**2))*rad
     1		*parity(3)
	  angl(4)=rad*atan2(thphi(2),thphi(1))*parity(4)
	else
c NOT bisecting mode - calculate in fixed phi mode
	  angl(2)=asin(thphi(2)/dinv)*rad*parity(2)
	  angl(3)=atan2(thphi(3),thphi(1))*rad*parity(3)
	  angl(4)=0.0
	endif
c check software limits on motor positions
	call softlimits(angl,0,ierr)
	RETURN
	END
	function sgn(x)
c
	if(x.gt.1.0e-4)then
c x .gt. 0.0
		sgn=1.0
	else
c x .le. 0.0
		sgn=-1.0
	endif
	return
	end
	SUBROUTINE CALPHI(TH,THPHI)
	DIMENSION TH(3),THPHI(3)
	include 'single.inc'
	DO 110 I=1,3
	  THPHI(I)=0.0
	  DO 100 J=1,3
	    THPHI(I)=THPHI(I)+UB(I,J)*TH(J)
100	  CONTINUE
110	CONTINUE
	RETURN
	END
	SUBROUTINE OMTOB(iprnt)
C
C CALCULATE B FROM ORIENTATION MATRIX
C
	include 'single.inc'
	DIMENSION TEMP1(3,3),TEMP2(3,3),temp(3,3)
	DATA RAD/57.2958/
C
C GET UB TRANSPOSE, calculate G=UB(T)*UB and invert
C
	call mattsp(ub,temp1)
	CALL MATMUL(TEMP1,UB,TEMP)
	CALL MINV(TEMP,TEMP2,DET)
	IF(DET.LE.0.0)GO TO 130
	cell(7)=1.0/sqrt(det)
	DO 120 I=1,3
	  CELL(I)=SQRT(TEMP2(I,I))
120	CONTINUE
	CELL(4)=RAD*ACOS(TEMP2(2,3)/(CELL(2)*CELL(3)))
	CELL(5)=RAD*ACOS(TEMP2(1,3)/(CELL(1)*CELL(3)))
	CELL(6)=RAD*ACOS(TEMP2(1,2)/(CELL(1)*CELL(2)))
	if(iprnt.eq.0)TYPE 2,CELL
2	FORMAT('0Cell Constants: ',3F8.4,3F8.3,f9.3)
	CALL CALCB(CELL,B)
	RETURN
130	TYPE 4
4	FORMAT(' Invalid Orientation Matrix - Use COM')
	type *,'UB:',ub
	type *,'B(rec):',temp
	type *,'B(dir):',temp2
	RETURN
	END
	SUBROUTINE PUSHOM
C ROUTINE TO PUSH CURRENT ORIENTATION MATRIX FROM 'UB' TO 'UBO'
	include 'single.inc'
C SAVE CURRENT OM IN 'UBO'
	DO 90 I=1,3
	  DO 80 J=1,3
	    UBO(J,I)=UB(J,I)
80	  CONTINUE
90	CONTINUE
	RETURN
	END
	SUBROUTINE POPOM
C ROUTINE TO POP ORIENTATION MATRIX FROM 'UBO' TO 'UB'
	include 'single.inc'
C RESTORE OLD OM FROM 'UBO'
	DO 90 I=1,3
	  DO 80 J=1,3
	    UB(J,I)=UBO(J,I)
80	  CONTINUE
90	CONTINUE
	RETURN
	END
	FUNCTION ANINT(X)
C
C FUNCTION TO GET NEAREST INTEGRAL VALUE FOR ARGUMENT X
C
	TEMP=X+0.5*SGN(X)
	ANINT=AINT(X)
	RETURN
	END
	FUNCTION ASIN(X)
	ASIN=ATAN2(X,SQRT(1.0-X*X))
	RETURN
	END
	FUNCTION ACOS(X)
	ACOS=ATAN2(SQRT(1.0-X*X),X)
	RETURN
	END
	SUBROUTINE MATMUL(A,B,C)
C
C ROUTINE TO PERFORM MATRIX MULTIPLICATION C=A*B
C
	DIMENSION A(3,3),B(3,3),C(3,3)
C
	DO 120 I=1,3
	  DO 110 J=1,3
	    C(J,I)=0.0
	    DO 100 K=1,3
	      C(J,I)=C(J,I)+A(J,K)*B(K,I)
100	    CONTINUE
110	  CONTINUE
120	CONTINUE
	RETURN
	END
	SUBROUTINE ORTHOG(V1,V2,X)
C
C ROUTINE TO CREATE ORTHOGONAL MAXTRIX X FROM V1 AND V2
C
	DIMENSION V1(3),V2(3),X(3,3)
C
C COPY V1 TO FIRST COLUMN OF X AND NORMALIZE
C
	DO 100 I=1,3
	  X(I,1)=V1(I)
100	CONTINUE
	CALL UNITL(X(1,1))
C
C CALCULATE CROSS PRODUCT OF V1 AND V2 AND NORMALIZE
C
	CALL VCROSS(V1,V2,X(1,3))
	CALL UNITL(X(1,3))
C
C CALCULATE CROSS OF COL 3 AND COL 1 OF X
C
	CALL VCROSS(X(1,3),X(1,1),X(1,2))
	RETURN
	END
	SUBROUTINE VCROSS(U,V,W)
C
C CALCULATE CROSS PRODUCT - W = U x V
C
	DIMENSION U(3),V(3),W(3)
C
	W(1)=U(2)*V(3)-V(2)*U(3)
	W(2)=V(1)*U(3)-U(1)*V(3)
	W(3)=U(1)*V(2)-V(1)*U(2)
	RETURN
	END
	SUBROUTINE UNITL(V)
C
C REPLACE VECTOR V WITH PARALLEL UNIT VECTOR
C
	DIMENSION V(3)
C
	SUM=SQRT(V(1)**2+V(2)**2+V(3)**2)
	IF(SUM.LT.1.0E-6)SUM=1.0E-6
	DO 100 I=1,3
	  V(I)=V(I)/SUM
100	CONTINUE
	RETURN
	END
	SUBROUTINE MINV(A,B,DET)
C
C INVERT MATRIX A INTO B, A AND B MAY BE SAME MATRIX
C
	DIMENSION A(3,3),B(3,3),C(3,3),work(6)
C
c copy A to B
c
	do i=1,3
	    do j=1,3
		b(j,i)=a(j,i)
	    enddo
	enddo
c Call MXLNEQ to do the work
	call mxlneq(b,3,3,det,jrank,1.0e-6,work,0)
	if(jrank.ne.3)det=0.0
	RETURN
	END
	SUBROUTINE CALCB(CELL,B)
C
C ROUTINE TO FORM MATRIX B FROM CELL PARAMETERS
C
	DIMENSION CELL(6),B(3,3),G(3,3),RCELL(6)
	DATA RAD/57.2958/
C
	DO I=1,2
	  G(I,I)=CELL(I)**2
	  DO J=I+1,3
	    G(I,J)=CELL(I)*CELL(J)*COS(CELL(9-I-J)/RAD)
	    G(J,I)=G(I,J)
	  enddo
	enddo
	G(3,3)=CELL(3)**2
	CALL MINV(G,G,DET)
	if(det.le.0.0)then
	    type 10,det
10	    format(' Error in Cell: Determinant of G = ',f12.3)
	    return
	endif
	DO I=1,3
	  RCELL(I)=SQRT(G(I,I))
	enddo
	RCELL(4)=G(2,3)/(RCELL(2)*RCELL(3))
	RCELL(5)=G(1,3)/(RCELL(1)*RCELL(3))
	RCELL(6)=G(1,2)/(RCELL(1)*RCELL(2))
	SNGM=SQRT(1.0-RCELL(6)**2)
	B(1,1)=RCELL(1)
	B(2,1)=0.0
	B(3,1)=0.0
	B(1,2)=RCELL(2)*RCELL(6)
	B(2,2)=RCELL(2)*SNGM
	B(3,2)=0.0
	B(1,3)=RCELL(3)*RCELL(5)
	TEMP=(RCELL(4)-RCELL(5)*RCELL(6))/SNGM
	B(2,3)=RCELL(3)*TEMP
	B(3,3)=RCELL(3)*SQRT(1.0-RCELL(5)**2-TEMP**2)
	RETURN
	END
	subroutine mattsp(a,b)
c
c routine to transpose 3x3 matrix A and output it in B, A and B may be
c    the same
	dimension a(3,3),b(3,3),c(3,3)
c
	do 110 i=1,3
	  do 100 j=1,3
	    c(j,i)=a(i,j)
100	  continue
110	continue
	do 130 i=1,3
	  do 120 j=1,3
	    b(j,i)=c(j,i)
120	  continue
130	continue
	return
	end
      subroutine psimat(angl,pmat)
c
c Routine to calculate matrix to rotate about the PSI direction
c
      real angl(4),pmat(3,3),rad,coso,sino,cosx,sinx,sinp,cosp
      include 'single.inc'
c
c 'ANGL' contains diffractometer angles for PSI = 0
c 'PMAT' will be desired rotation matrix
c
      data rad/57.29578/
c
c Get sines and cosines of Euler angles
c
      coso=cos(angl(2)/rad)
      cosx=cos(angl(3)/rad)
      cosp=cos(angl(4)/rad)
      sino=sin(angl(2)/rad)*parity(2)
      sinx=sin(angl(3)/rad)*parity(3)
      sinp=sin(angl(4)/rad)*parity(4)
c fill in matrix
      pmat(1,1) =  coso*cosx*cosp-sino*sinp
      pmat(1,2) =  coso*cosx*sinp+sino*cosp
      pmat(1,3) =  coso*sinx
      pmat(2,1) = -sino*cosx*cosp-coso*sinp
      pmat(2,2) = -sino*cosx*sinp+coso*cosp
      pmat(2,3) = -sino*sinx
      pmat(3,1) = -sinx*cosp
      pmat(3,2) = -sinx*sinp
      pmat(3,3) =  cosx
      return
      end
      subroutine psiang(pmat,psi,angl,ierr)
c
c routine to calculate diffractometer angles for a given value of PSI
c
c 'PMAT' is rotation matrix from routine PSIMAT
c 'PSI' is desired value of PSI, in degrees
c 'ANGL(1)' contains value of two theta at entry - NOT used in this routine
c 'ANGL(2-4)' is replaced with new Euler angles (omega, chi,phi) 
c         at desired value of PSI
c 'IERR' is non-zero if reflection is inaccessible
c
      real pmat(3,3),psi,angl(4),cospsi,sinpsi,psrot(3,3),rad
     1,psimat(3,3)
      integer i
      include 'single.inc'
      data rad/57.29578/
c
      cospsi=cos(psi/rad)      
      sinpsi=sin(psi/rad)
c
c Apply rotation to pmat - build psimat
c
      psimat(1,1)=1.0
      psimat(2,1)=0.0
      psimat(3,1)=0.0
      psimat(1,2)=0.0
      psimat(2,2)=cospsi
      psimat(3,2)=-sinpsi
      psimat(1,3)=0.0
      psimat(2,3)=sinpsi
      psimat(3,3)=cospsi
c do matrix operation
      call matmul(psimat,pmat,psrot)
c
c Solve for new diffractometer angles
c
c Check for special case of cos(chi) = +/- 1
c
      if(abs(abs(psrot(3,3))-1.0).le.0.000002)then
c Omega and phi are colinear
	angl(2)=0.0
	angl(3)=acos(psrot(3,3))*rad*parity(3)
	angl(4)=atan2(-psrot(1,1),psrot(1,2))*rad*parity(4)
     1	  +90.0*parity(2)
      else
c  Omega
        angl(2)=atan2(-psrot(2,3),psrot(1,3))*rad*parity(2)
c Chi
        temp=sqrt(psrot(3,1)**2+psrot(3,2)**2)
        angl(3)=atan2(temp,psrot(3,3))*rad*parity(3)
c Phi      
        angl(4)=atan2(-psrot(3,2),-psrot(3,1))*rad*parity(4)
      endif
c
c Test angles for in bounds conditions
c
      if(abs(angl(2)).gt.90.0)then
c Omega magnitude greater than 90 - select other set
        angl(2)=angl(2)+180.0
        angl(3)=-angl(3)
        angl(4)=angl(4)+180.0
c Reset all angles to be between -180 and 180
        do i=2,4
          if(angl(i).gt.180.0)angl(i)=angl(i)-360.0
        enddo
      endif
c check if angles can be reached
      call softlimits(angl,0,ierr)
      return
      end
      subroutine face
c
c Routine to calculate equations of faces from h,k,l,d and UB
c
	implicit none
	real th(3),thphi(3),d
	integer nc
c
	write(*,'(a)')' Enter h,k,l, and distance'
	do while (.true.)
	    write(*,'(a)')'$? '
	    read(*,'(q,4f8.0)')nc,th,d
	    if (nc.eq.0) return
	    call calphi(th,thphi)	! Get vector on diffractometer
	    call unitl(thphi)		! Convert to unit length
	    write(*,'(a,4f9.5)')' Coefficients: ',thphi,d
	enddo
	end	
