      PROGRAM INDLSQ
      integer mxl, mxr
      parameter (mxl=500, mxr=2000)
C LEAST SQUARES UNIT CELL REFINEMENT INCLUDING INDEXING IN FORTRAN IV
      DIMENSION ICLASS(10),ICOND(10),JCLASS(27),JCOND(27),XO(3)
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/ADATA/IALF(19),IFI(2,27),IFJ(3,27),IHKL(6,3),JHKL(10,3)
      COMMON/BLKA/IDG(3),FMN(3),AO(3),UO(3),ALO(3),AR(3),A(27),B(6)
      COMMON/BLKB/H(3),U(3),X(3),W(3),AD(3),ALS(3),CSA(3),SNA(3),CR(6)
      COMMON/BLKC/THCMP,D,CONST,SWT,DIFF,JA,DOBS,PWL,ITHT,THTSE,THOBS,
     1 RSDMAX,IRSD,MMN,E(7),SUM
      COMMON/BLKD/TT(6,6),ST(6,6),SERCE(6),SEDCE(6),DD(3),VCME(6,6),
     1 SS(21)
      COMMON/BLKE/JTYP,HH(mxr,3),DCM(mxr),IOP(27),JOP(27),NOPT,DMIN,NR
      COMMON/BLKF/OBSTH(mxl),SW(mxl),WL(mxl),DSP(mxl),HG(mxl,3)
      COMMON/BLKG/IHEAD(17),ICYCLE,IMN,MM,IALF12
      COMMON/BLKI/Z(6)
      COMMON/BLKJ/TOLMN
      COMMON /BLKH/       NCYC,AA(7),DLST(mxl),IOC(mxl),THTMX,TOL,tolp
      DATA ICLASS/3HHKL,3HHHL,3H0KL,3HH0L,3HHK0,3HHH0,3HH00,3H0K0,3H00L,
     *3HHLL/, ICOND/4H H+K,4H H+L,4H K+L,4H HKL,4H-HKL,4H   H,4H   K,4H
     *  L,4H2H+L,4HH+2L/,JCLASS/6*1,4*2,3*3,3*4,3*5,6,7,8,9,4*10/, JCOND
     */1,2,3,0,4,5,6,8,2,9,7,8,3,6,8,2,6,7,1,6,6,7,8,6,8,2,10/
C FIRST CARD OF DATA SET SHOULD BE A TITLE CARD
    5 READ(5,6) ITEST,(IHEAD(L),L=1,17)
    6 FORMAT(18A4)
      IF (ITEST.EQ.IALF(13)) STOP
      IF (ITEST.EQ.IALF(9)) GO TO 15
      WRITE(6,7) ITEST,IALF(9)
    7 FORMAT(1H0,A4,27H IN COLUMNS 1-6 INSTEAD OF ,A4)
   10 READ(5,6) ITEST,(IHEAD(L),L=1,17)
      IF (ITEST.EQ.IALF(13)) STOP
      IF (ITEST.NE.IALF(9)) GO TO 10
   15 WRITE(6,16) (IHEAD(L),L=1,17)
   16 FORMAT(1H1,5X,17A4)
      READ(5,17) ITEST,(AO(L),L=1,3),(IDG(L),FMN(L),L=1,3),ITHT,THTMX,
     * NCYC,TOLMN,TOLMX
   17 FORMAT(A4,2X,3F9.6,3(I3,F5.3),I1,F3.0,I1,2F5.4)
      IF (ITEST.EQ.IALF(13)) STOP
      IF (ITEST.EQ.IALF(10)) GO TO 20
      WRITE(6,7) ITEST,IALF(10)
      GO TO 10
   20 READ(5,21) ITEST,ICRYST,THEMX,PWL,(IOP(L),L=1,27),SEMULT
   21 FORMAT(A4,2X,A4,1X,F3.0,F6.5,27I1,23X,F2.1)
      IF (ITEST.EQ.IALF(13)) STOP
      IF (ITEST.EQ.IALF(11)) GO TO 25
      WRITE(6,7) ITEST,IALF(11)
      GO TO 10
   25 CALL PRMCHK
      IF (ITYP.GE.1.AND.ITYP.LE.7) GO TO 30
      WRITE(6,26) ICRYST
   26 FORMAT(33H0CRYSTAL SYSTEM ON SYSE GIVEN AS ,A4)
      GO TO 10
   30 SUM=0.0
      DO 35 L=1,3
      IF (AO(L).LE.0.0) GO TO 40
      ALO(L)=FLOAT(IDG(L))*60.0+FMN(L)
      IF (ALO(L).LE.0.0.OR.ALO(L).GE.10800.0) GO TO 40
   35 SUM=SUM+ALO(L)
      IF (SUM.LT.21600.0) GO TO 45
   40 WRITE(6,41) (AO(L),L=1,3),(IDG(L),FMN(L),L=1,3)
   41 FORMAT(26H0GIVEN PARAMETERS IN ERROR,3F11.5,2X,3(I6,F7.3))
      GO TO 10
   45 ITHT=MOD(ITHT,2)
      IF (ITHT.NE.0) GO TO 50
      THTMX=0.5*THTMX
      TOLMN=0.5*TOLMN
      TOLMX=0.5*TOLMX
      THEMX=0.5*THEMX
   50 IF (THTMX.LE.0.0) THTMX=20.0
      IF (NCYC.LT.1) NCYC=2
      IF (TOLMN.LE.0.0) TOLMN=0.02
      IF (TOLMX.LE.TOLMN) TOLMX=5.0*TOLMN
      IF (PWL.LE.0.0) PWL=1.5405
      IF (SEMULT .LE. 0.0) SEMULT=2.0
      IALF12=IALF(ITHT+14)
      CALL OPSET
      IMN=((N+3)*N)/2
      DO 55 L=1,3
      AD(L)=AO(L)
      ALS(L)=ALO(L)
      ARAD=0.00029088821*ALO(L)
      CSA(L)=COS(ARAD)
   55 SNA(L)=SIN(ARAD)
      CALL RSCMP
      DO 60 L=1,3
      UO(L)=U(L)
   60 XO(L)=X(L)
      VRO=VR
      VCO=VC
      CALL DCPRN
      CALL RCPRN
      CALL RDHKL
      MM=M
      THEMXO=0.0
      IF(M.GT.0) THEMXO=OBSTH(M)+TOLMX
      THEMX=AMAX1(THEMX,THEMXO)
      IF (THEMX.LE.0.0.OR.THEMX.GT.90.0) THEMX=90.0
      DMIN=0.5*PWL/SIN(  0.017453293*THEMX)
      WRITE(6,76) THTMX,NCYC,TOLMN,TOLMX,THEMX,DMIN,ICRYST
   76 FORMAT(24H01-THETA ANGLES  THTMX =,F6.1,8H  NCYC =,I2,9H  TOLMN =,
     * F7.4,9H  TOLMX =,F7.4,9H  THEMX =,F6.1,8H  DMIN =,F10.6,
     * 2X,A4,7H SYSTEM)
      WRITE(6,77) NOPT
   77 FORMAT(1H0,I2,41H CONDITIONS FOR NON-EXTINCTION CALLED FOR,5X,
     * 18HCLASS CONDITION(S))
      IF (NOPT .LT. 1) GO TO 79
      DO 78 L=1,NOPT
      J=JOP(L)
      I=JCLASS(J)
      K=JCOND(J)
      IF (K .GE. 1 .AND. K .LE. 10) GO TO 782
      WRITE(6,781) IOP(J),IOP(J),IOP(J)
  781 FORMAT(50X,3HHKL,5X,4HH+K=,I1,6HN,H+L=,I1,6HN,K+L=,I1,1HN)
      GO TO 78
  782 WRITE(6,783) ICLASS(I),ICOND(K),IOP(J)
  783 FORMAT(50X,A4,A4,1H=,I1,1HN)
   78 CONTINUE
   79 CALL HKLGEN
      CALL SORTD(NR,DCM,HH,H)
      CALL SKWEZE
      TOL=TOLMX
      ICYCLE=0
   80 CALL SLCR
      IF (MM .LT. N) GO TO 92
      MMN=M-N
      IF (MMN) 90,99,94
   90 IF ((ICYCLE.GE.10).OR.(ICYCLE.GT.NCYC.AND.TOL.LE.TOLMN))GO TO 92
      TOL=0.5*TOL
      IF (TOL .LT. TOLMN) TOL=TOLMN
      NCYC=0
      GO TO 80
   92 WRITE(6,93)
   93 FORMAT(12H0NO SOLUTION)
      GO TO 150
   94 J=N+1
      K=J
      DO 98 I=1,N
      AA(I)=A(K)
      J=J-1
   98 K=K+J
   99 CALL PDSQRT(A,N,B,IERR)
      IF (IERR .EQ. 0) GO TO 104
      WRITE(6,100) IERR
  100 FORMAT(38H0SINGULAR NORMAL EQUATION MATRIX - ROW,I2,
     * 36H LINEARLY DEPENDENT ON PREVIOUS ROWS)
      GO TO 90
  104 CALL COORD
      DO 120 J=1,3
      X(J)=X(J)+CR(J)
      U(J)=U(J)+CR(J+3)
      IF(X(J).GT.0.) GOTO 109
      WRITE(6,108) J
  108 FORMAT(40H0CORRECTED VALUE OF RECIPROCAL CELL EDGE,I2,
     * 27H LESS THAN OR EQUAL TO ZERO)
      GO TO 114
  109 IF(ABS(U(J)).LT.1.) GOTO 120
      WRITE(6,112) J
  112 FORMAT(42H0CORRECTED VALUE OF RECIPROCAL CELL COSINE,I2,
     * 40H EQUALS OR EXCEEDS 1.0 IN ABSOLUTE VALUE)
  114 DO 116 K=1,J
      X(K)=X(K)-CR(K)
  116 U(K)=U(K)-CR(K+3)
      GO TO 90
  120 CONTINUE
      WRITE(6,122)
  122 FORMAT(1H0,24X,1HA,15X,1HB,15X,1HC,13X,5HALPHA,7X,4HBETA,8X,
     * 5HGAMMA,10X,6HVOLUME)
      CALL CMPNUP
      IF (MMN .GT. 0) GO TO 130
      TOL=0.5*TOL
      GO TO 140
  130 DO 134 I=1,N
  134 SUM=SUM-AA(I)*B(I)
      IF (SUM .LT. 0.0) SUM=0.0
      THTSE=SQRT(SUM/FLOAT(MMN))
      CALL RSDPRN
      TOL=SEMULT*THTSE
  140 IF (TOL .LT. TOLMN) TOL=TOLMN
      IF (TOL .GT. TOLMX) TOL=TOLMX
      IF (ICYCLE .LE. NCYC) GO TO 160
      CRMAX=0.0
      DO 146 I=1,N
  146 IF (ABS(B(I)) .GT. CRMAX) CRMAX=ABS( B(I))
      IF (CRMAX .GE. 0.5E-6) GO TO 152
  149 IF (MMN .GT. 0) CALL SESUB
  150 CALL FNLST
      GO TO 5
  152 IF (ICYCLE .LT. 10) GO TO 160
      WRITE(6,154)
  154 FORMAT(28H0NOT CONVERGENT IN 10 CYCLES)
      GO TO 149
  160 CALL RCMPD
      CALL SORTD(NR,DCM,HH,H)
      GO TO 80
      END
      SUBROUTINE PRMCHK
C IMPOSES SYMMETRY CONDITIONS ON DIRECT CELL ELEMENTS,SETS ITYP AND
C JTYP ACCORDING TO CRYSTAL SYSTEM
      integer mxl, mxr
      parameter (mxl=500, mxr=2000)
      COMMON/ADATA/IALF(19),IFI(2,27),IFJ(3,27),IHKL(6,3),JHKL(10,3)
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/BLKA/IDG(3),FMN(3),AO(3),UO(3),ALO(3),AR(3),A(27),B(6)
      COMMON/BLKE/JTYP,HH(mxr,3),DCM(mxr),IOP(27),JOP(27),NOPT,DMIN,NR
      ITYP=0
      IF (ICRYST.NE.IALF(1)) GO TO 13
      JTYP=1
      ITYP=1
      N=1
      AO(3)=AO(1)
    4 AO(2)=AO(1)
    5 DO 10 L=1,3
      IDG(L)=90
   10 FMN(L)=0.0
      RETURN
   13 IF (ICRYST.NE.IALF(7)) GO TO 21
      JTYP=2
   15 ITYP=2
      N=2
      AO(2)=AO(1)
      DO 18 L=1,3
      IDG(L)=90
   18 FMN(L)=0.0
      IDG(3)=120
      RETURN
   21 IF (ICRYST.NE.IALF(8)) GO TO 24
      JTYP=8
      GO TO 15
   24 IF (ICRYST.NE.IALF(3)) GO TO 30
      JTYP=3
      ITYP=3
      N=3
      GO TO 5
   30 IF(ICRYST.NE.IALF(4)) GO TO 38
      JTYP=4
      ITYP=4
      N=4
      DO 36 L=1,3,2
      IDG(L)=90
   36 FMN(L)=0.0
      RETURN
   38 IF (ICRYST.NE.IALF(2)) GO TO 42
      JTYP=5
      ITYP=5
      N=2
      GO TO 4
   42 IF (ICRYST.NE.IALF(5)) GO TO 51
      JTYP=7
      ITYP=7
      N=2
      DO 49 L=2,3
      AO(L)=AO(1)
      IDG(L)=IDG(1)
   49 FMN(L)=FMN(1)
      RETURN
   51 IF (ICRYST.NE.IALF(6)) RETURN
      JTYP=6
      ITYP=6
      N=6
      END
      SUBROUTINE OPSET
      integer mxl, mxr
      parameter (mxl=500, mxr=2000)
C DOES BOOKKEEPING FOR SETTING UP EXTINCTION RELATIONS TO USE
      COMMON/BLKE/JTYP,HH(mxr,3),DCM(mxr),IOP(27),JOP(27),NOPT,DMIN,NR
      NOPT=0
      IF (JTYP.EQ.8) IOP(6)=3
      DO 9 L=1,27
      IF (IOP(L).LT.2.OR.IOP(L).GT.6.OR.IOP(L).EQ.5) GO TO 9
      NOPT=NOPT+1
      JOP(NOPT)=L
    9 CONTINUE
      END
      SUBROUTINE RSCMP
C COMPUTES RECIPROCAL CELL ELEMENTS FROM DIRECT CELL ELEMENTS
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/BLKA/IDG(3),FMN(3),AO(3),UO(3),ALO(3),AR(3),A(27),B(6)
      COMMON/BLKB/H(3),U(3),X(3),W(3),AD(3),ALS(3),CSA(3),SNA(3),CR(6)
      DO 10 I=1,3
      J=MOD(I,3)+1
      K=MOD(J,3)+1
      U(I)=(CSA(J)*CSA(K)-CSA(I))/(SNA(J)*SNA(K))
      W(I)=SQRT(1.0-U(I)**2)
      X(J)=1.0/(AD(J)*W(I)*SNA(K))
   10 AR(I)=3437.7468*(1.5707963-ATAN2(U(I),W(I)))
      VR=X(1)*X(2)*X(3)*SNA(1)*W(2)*W(3)
      VC=1.0/VR
      END
      SUBROUTINE DCPRN
C PRINTS OUT DIRECT CELL ELEMENTS
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/BLKA/IDG(3),FMN(3),AO(3),UO(3),ALO(3),AR(3),A(27),B(6)
      COMMON/BLKB/H(3),U(3),X(3),W(3),AD(3),ALS(3),CSA(3),SNA(3),CR(6)
      DO 5  L=1,3
    5 CALL MTDM(ALS(L),IDG(L),FMN(L))
      WRITE(6,7) (AD(L),L=1,3),(IDG(L),FMN(L),L=1,3),VC
    7 FORMAT(16H     DIRECT CELL,1X,3E16.8,1X,3(I5,F7.3),E18.8)
      END
      SUBROUTINE RCPRN
C PRINTS OUT RECIPROCAL CELL ELEMENTS
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/BLKA/IDG(3),FMN(3),AO(3),UO(3),ALO(3),AR(3),A(27),B(6)
      COMMON/BLKB/H(3),U(3),X(3),W(3),AD(3),ALS(3),CSA(3),SNA(3),CR(6)
      DO 6 I=1,3
      AR(I)=3437.7468*(1.5707963-ATAN2(U(I),SQRT(1.0-U(I)**2)))
    6 CALL MTDM(AR(I),IDG(I),FMN(I))
      WRITE(6,8) (X(I),I=1,3),(IDG(I),FMN(I),I=1,3),VR
    8 FORMAT(16H RECIPROCAL CELL,1X,3E16.8,1X,3(I5,F7.3),E18.8)
      END
      SUBROUTINE RDHKL
      integer mxl, mxr
      parameter (mxl=500, mxr=2000)
C READS HKL CARDS AND STORES OBSERVED DATA
      DIMENSION H(3),IH(3)
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/ADATA/IALF(19),IFI(2,27),IFJ(3,27),IHKL(6,3),JHKL(10,3)
      COMMON/BLKC/THCMP,D,CONST,SWT,DIFF,JA,DOBS,PWL,ITHT,THTSE,THOBS,
     * RSDMAX,IRSD,MMN,E(7),SUM
      COMMON/BLKE/JTYP,HH(mxr,3),DCM(mxr),IOP(27),JOP(27),NOPT,DMIN,NR
      COMMON/BLKG/IHEAD(17),ICYCLE,IMN,MM,IALF12
      COMMON/BLKF/OBSTH(mxl),SW(mxl),WL(mxl),DSP(mxl),HG(mxl,3)
      M=0
    9 READ(5,10)ITEST,(IH(L),L=1,3),D,WLC,OBS,WT
   10 FORMAT(A4,2X,3I5,F9.6,F8.6,2F9.5)
      IF (ITEST.EQ.IALF(12)) GO TO 64
      IF (WLC.GT.0.0) PWL=WLC
      IF (WT.LE.0.0) WT=1.0
      IF (OBS.GT.0.0) GO TO 20
      IF (D.GT.0.0) GO TO 25
   16 WRITE(6,17)(IH(L),L=1,3),D,PWL,OBS,WT
   17 FORMAT(18H0CARD WITH INDICES,3I5,5H  D =,F10.6,6H  WL =,F9.6,
     * 13H  OBS ANGLE =,F9.4,8H  WGHT =,F10.5,15H  NOT PROCESSED)
      GO TO 9
   20 THOBS=OBS
      IF (ITHT.EQ.0) THOBS=0.5*THOBS
      IF (THOBS.GE.90.0) GO TO 16
      DOBS=0.5*PWL/SIN(0.017453293*THOBS)
      GO TO 29
   25 CNST=0.5*PWL
      THOBS=THETA(CNST,D)
      IF (THOBS.LE.0.0.OR.THOBS.GE.90.0) GO TO 16
      DOBS=D
   29 IF (NOPT.LT.1) GO TO 45
      IF (IH(1).EQ.0.AND.IH(2).EQ.0.AND.IH(3).EQ.0)    GO TO 45
C CHECKING GIVEN FIXED REFLECTION FOR CONSISTENCY WITH CONDITIONS
C FOR NON-EXTINCTION AS GIVEN ON SYSEXT CARD
      CALL CHKEXT(IH(1),IH(2),IH(3),JEXT)
      IF (JEXT.EQ.0) GO TO 45
      WRITE(6,40) (IH(L),L=1,3)
   40 FORMAT(20H0REFLECTION WITH H =,I5,5H  K =,I5,5H  L =,I5,71H  NOT P
     *ROCESSED - SHOULD NOT EXIST ACCORDING TO SYSEXT CARD INFORMATION)
      GO TO 9
   45 M=M+1
      IF (M.LE.mxl) GO TO 52
c     WRITE(6,50)
c  50 FORMAT(50H0OVER 200 OBSERVATIONS GIVEN - ONLY FIRST 200 USED)
      write(6, '(a, i3, a, i3, a)')
     &  '0OVER ', mxl, ' OBSERVATIONS GIVEN - ONLY FIRST ',
     &  mxl, ' USED'
      M=mxl
      GO TO 64
   52 OBSTH(M)=THOBS
      DSP(M)=DOBS
      WL(M)=PWL
      SW(M)=SQRT(WT)
      DO 60 L=1,3
   60 HG(M,L)=IH(L)
      GO TO 9
C SORT GIVEN REFLECTIONS SO AS TO BE IN DECREASING VALUE FOR D-SPACING
   64 I=M
   65 I=I-1
      IF (I.LT.1) RETURN
      DO 80 J=1,I
      IF (DSP(J).GE.DSP(J+1)) GO TO 80
      DOBS=DSP(J)
      DSP(J)=DSP(J+1)
      DSP(J+1)=DOBS
      THOBS=OBSTH(J)
      OBSTH(J)=OBSTH(J+1)
      OBSTH(J+1)=THOBS
      WLC=WL(J)
      WL(J)=WL(J+1)
      WL(J+1)=WLC
      SWT=SW(J)
      SW(J)=SW(J+1)
      SW(J+1)=SWT
      DO 75 K=1,3
      H(K)=HG(J,K)
      HG(J,K)=HG(J+1,K)
   75 HG(J+1,K)=H(K)
   80 CONTINUE
      GO TO 65
      END
      SUBROUTINE SLCR
      integer mxl, mxr
      parameter (mxl=500, mxr=2000)
C SELECTS REFLECTIONS FOR CURRENT CYCLE AND ATTEMPTS TO INDEX
C REFLECTIONS NOT ALREADY INDEXED
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/ADATA/IALF(19),IFI(2,27),IFJ(3,27),IHKL(6,3),JHKL(10,3)
      COMMON/BLKA/IDG(3),FMN(3),AO(3),UO(3),ALO(3),AR(3),A(27),B(6)
      COMMON/BLKB/H(3),U(3),X(3),W(3),AD(3),ALS(3),CSA(3),SNA(3),CR(6)
      COMMON/BLKC/THCMP,D,CONST,SWT,DIFF,JA,DOBS,PWL,ITHT,THTSE,THOBS,
     * RSDMAX,IRSD,MMN,E(7),SUM
      COMMON/BLKE/JTYP,HH(mxr,3),DCM(mxr),IOP(27),JOP(27),NOPT,DMIN,NR
      COMMON/BLKF/OBSTH(mxl),SW(mxl),WL(mxl),DSP(mxl),HG(mxl,3)
      COMMON/BLKG/IHEAD(17),ICYCLE,IMN,MM,IALF12
      COMMON/BLKH/ NCYC,AA(7),DLST(mxl),IOC(mxl),THTMX,TOL,tolp
      COMMON/BLKI/Z(6)
      COMMON/BLKJ/TOLMN
      ICYCLE=ICYCLE+1
      TOLP=TOL
      IF (ITHT .EQ. 0) TOLP=TOLP+TOLP
      WRITE(6,2) ICYCLE,(IHEAD(L),L=1,17),IALF12,TOLP,ICRYST
    2 FORMAT(6H1CYCLE,I3,4X,17A4,4X,A2,17HTHETA TOLERANCE =,F9.5,4X,A4)
      M=0
      IRSD=0
      SUM=0.0
      RSDMAX=0.0
      WRITE(6,8) IALF12,IALF12,IALF12
    8 FORMAT(4H0  N,4X,1HH,4X,1HK,4X,1HL,5X,6HD CALC,8X,5HD OBS,7X,
     * 6HLAMBDA,5X,A2,10HTHETA CALC,4X,A2,9HTHETA OBS,4X,A2,
     * 10HTHETA DIFF,7X,6HWEIGHT)
      IF (MM .LT. 1) RETURN
      DO 10 L=1,IMN
   10 A(L)=0.0
      IHP=1
      DO 100 IGP=1,MM
      JA=IALF(15)
      SWT=SW(IGP)
      THOBS=OBSTH(IGP)
      ISUM=0
      DO 20 L=1,3
      H(L)=HG(IGP,L)
   20 IF (H(L) .NE. 0.0) ISUM=ISUM+1
      PWL=WL(IGP)
      DOBS=DSP(IGP)
      DMNDF=ABS(DCM(IHP)-DOBS)
      JHP=IHP+1
   28 IF (JHP .GT. NR) GO TO 32
      IF (ABS(DCM(JHP)-DOBS) .GT. DMNDF) GO TO 32
      DMNDF=ABS(DCM(JHP)-DOBS)
      JHP=JHP+1
      GO TO 28
   32 IHP=JHP-1
      IOC(IGP)=IHP
      IF (ISUM .NE. 0) GO TO 70
      D=DCM(IHP)
      DLST(IGP)=D
      CONST=0.5*PWL
      IF (D .LE. CONST) GO TO 99
      THCMP=THETA(CONST,D)
      IF (ABS((THOBS-THCMP)*SWT) .GT. TOL) GO TO 99
      J=IGP-1
      ASSIGN 50 TO JADRS
   40 IF (J .LT. 1 .OR. J .GT. MM) GO TO 48
      IF (ABS(OBSTH(J)-THOBS)*SWT .LT. TOLMN) GO TO 99
      CNST=0.5*WL(J)
      IF (D .LE. CNST) GO TO 48
      THTRY=THETA(CNST,D)
      IF (ABS((OBSTH(J)-THTRY)*SW(J)).LE. TOL) GO TO 99
   48 GO TO JADRS,(50,56)
   50 J=IGP+1
      ASSIGN 56 TO JADRS
      GO TO 40
   56 K=IHP-1
      ASSIGN 66 TO KADRS
   60 IF (K .LT. 1 .OR. K .GT. NR) GO TO 64
      DB=DCM(K)
      IF (DB .LE. CONST) GO TO 64
      THTRY=THETA(CONST,DB)
      IF (ABS(THTRY-THCMP)*SWT .LT. TOLMN) GO TO 99
      IF (ABS(THTRY-THOBS)*SWT .LE. TOL) GO TO 99
   64 GO TO KADRS,(66,68)
   66 K=IHP+1
      ASSIGN 68 TO KADRS
      GO TO 60
   68 IF (THOBS .GT. THTMX .AND. ICYCLE .LE. NCYC) GO TO 100
      DO 69 L=1,3
   69 H(L)=HH(IHP,L)
      GO TO 90
   70 CONST=0.5*PWL
      CALL QCOMP
      IF (Q .LE. 0.0) GO TO 99
      D=1.0/SQRT(Q)
      DLST(IGP)=D
      IF (CONST .GE. D) GO TO 99
      THCMP=THETA(CONST,D)
      IF  (ABS((THCMP-THOBS)*SWT) .GT. TOL) JA=IALF(16)
   90 DIFF=THCMP-THOBS
      IF (JA .EQ. IALF(16)) GO TO 98
      CALL OBSEQ
      CALL HKLPR
      CALL NORME
      GO TO 100
   98 CALL HKLPR
   99 IOC(IGP)=-IOC(IGP)
  100 CONTINUE
      END
      SUBROUTINE QCOMP
C COMPUTES RECIPROCAL OF SQUARE OF DISTANCE FROM ORIGIN FOR PLANE HKL
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/BLKI/Z(6)
      COMMON/BLKB/H(3),U(3),X(3),W(3),AD(3),ALS(3),CSA(3),SNA(3),CR(6)
      Q= (Z(1)*H(1)+Z(4)*H(2)+Z(6)*H(3))*H(1)
     1  +(Z(4)*H(1)+Z(2)*H(2)+Z(5)*H(3))*H(2)
     2  +(Z(6)*H(1)+Z(5)*H(2)+Z(3)*H(3))*H(3)
      END
      SUBROUTINE FNLST
      integer mxl, mxr
      parameter (mxl=500, mxr=2000)
C LISTS THEORETICAL REFLECTIONS IN DESCENDING ORDER FOR D-SPACING
C REFLECTIONS WITH FIXED MILLER INDICES INDICATED BY ***
      DIMENSION IH(3),JH(3)
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/ADATA/IALF(19),IFI(2,27),IFJ(3,27),IHKL(6,3),JHKL(10,3)
      COMMON/BLKC/THCMP,D,CONST,SWT,DIFF,JA,DOBS,PWL,ITHT,THTSE,THOBS,
     * RSDMAX,IRSD,MMN,E(7),SUM
      COMMON/BLKE/JTYP,HH(mxr,3),DCM(mxr),IOP(27),JOP(27),NOPT,DMIN,NR
      COMMON/BLKF/OBSTH(mxl),SW(mxl),WL(mxl),DSP(mxl),HG(mxl,3)
      COMMON/BLKG/IHEAD(17),ICYCLE,IMN,MM,IALF12
      COMMON/BLKH/ NCYC,AA(7),DLST(mxl),IOC(mxl),THTMX,TOL,tolp
      WRITE(6,2) (IHEAD(L),L=1,17),ICRYST
    2 FORMAT(1H1,17A4,35H HKL LISTING - *** REFERS TO FIXED,,
     * 13HR TO REJECTS ,A4)
      WL1=PWL
      I2=0
      IF (MM .LT. 1) GO TO 10
      WL1=WL(1)
      I2=IABS(IOC(1))
      D1=DSP(1)
   10  WL2=WL1
      IF (NR .LT. 1) GO TO 14
      D1=DCM(1)
   14 D2=D1
      IGP=1
      IHP=1
      WRITE(6,18) IALF12,IALF12,IALF12
   18 FORMAT(4H0  N,4X,1HH,4X,1HK,4X,1HL,5X,6HD CALC,8X,5HD OBS,7X,
     * 6HLAMBDA,5X,A2,10HTHETA CALC,4X,A2,9HTHETA OBS,4X,A2,
     * 10HTHETA DIFF,8X,6HWEIGHT)
   20 IF (IHP .LE. NR) GO TO 22
      IF (IGP .GT. MM) RETURN
      GO TO 40
   22 DO 24 L=1,3
   24 IH(L)=IFIX(HH(IHP,L)+SIGN(0.5,HH(IHP,L)))
      D=DCM(IHP)
      IF (IHP .EQ. I2) GO TO 40
      WLC=WL1
      IF (ABS(D-D2) .LT. ABS(D1-D)) WLC=WL2
      THCMP=THETA(0.5*WLC,D)
      IF (ITHT .EQ. 0) THCMP=THCMP+THCMP
      WRITE(6,9930) IHP,(IH(L),L=1,3),ialf(18),D,WLC,THCMP
 9930 FORMAT(1H ,I3,3I5,A2,F11.6,F25.6,F16.5)
   30 FORMAT(1H ,I3,3I5,F13.6,F25.6,F16.5)
   32 IHP=IHP+1
      GO TO 20
   40 DO 42 L=1,3
   42 JH(L)=IFIX(HG(IGP,L)+SIGN(0.5,HG(IGP,L)))
      DF=DLST(IGP)
      WT=SW(IGP)**2
      DOBS=DSP(IGP)
      THOBS=OBSTH(IGP)
      WLC=WL(IGP)
      THCLC=THETA(0.5*WLC,DF)
      DIFF=THCLC-THOBS
      IF (ITHT .NE. 0) GO TO 54
      THCLC=THCLC+THCLC
      THOBS=THOBS+THOBS
      DIFF=DIFF+DIFF
   54 JA=ialf(17)
      IF (IOC(IGP) .LT. 0) JA=IALF(16)
      if (ioc(igp) .lt. 0 .and. abs(diff) .le. tolp) ja = ialf(19)
      IF (JH(1).NE.0.OR.JH(2).NE.0.OR.JH(3).NE.0) GO TO 70
      WRITE(6,58) IHP,(IH(L),L=1,3),JA,DF,DOBS,WLC,THCLC,THOBS,DIFF,WT
   58 FORMAT(1H ,I3,3I5,A2,F11.6,F13.6,F12.6,F16.5,F15.4,F16.5,F14.5)
   60 IGP=IGP+1
      D1=D2
      WL1=WL2
      IF(IGP.GT.MM) GOTO 68
      IF (IABS(IOC(IGP)) .EQ. I2) GO TO 40
      I2=IABS(IOC(IGP))
      D2=DCM(I2)
      WL2=WL(IGP)
      GO TO 32
   68 I2=0
      GO TO 32
   70 THCMP=THETA(0.5*WLC,D)
      IF (ITHT .EQ. 0) THCMP=THCMP+THCMP
      IF (IGP .LT. 2) GO TO 75
      IF (IABS(IOC(IGP-1)) .NE. I2) GO TO 75
   73 WRITE(6,74) (JH(L),L=1,3),JA,DF,DOBS,WLC,THCLC,THOBS,DIFF,WT
   74 FORMAT(4H ***,3I5,A2,F11.6,F13.6,F12.6,F16.5,F15.4,F16.5,F14.5)
      GO TO 60
   75 IF(DLST(IGP).GT.D) GOTO 80
      WRITE(6,30) IHP,(IH(L),L=1,3),D,WLC,THCMP
      GO TO 73
   80 WRITE(6,74) (JH(L),L=1,3),JA,DF,DOBS,WLC,THCLC,THOBS,DIFF,WT
      WRITE(6,30) IHP,(IH(L),L=1,3),D,WLC,THCMP
      GO TO 60
      END
      SUBROUTINE PDSQRT(A,N,B,IERR)
C SOLUTION OF N LINEAR ALGEBRAIC EQUATIONS BY SQUARE ROOT METHOD.
C UPPER TRIANGLE OF GIVEN POSITIVE DEFINITE SYSTEM STORED SO THAT
C ELEMENT IN ROW I AND COLUMN J IS A((2*N+4-2)*(I-1)/2+J-2+U), WITH
C COLUMN N+U CONTAINING RIGHT HAND SIDE EQUATION CONSTANTS.
C B() IS SOLUTION VECTOR. IERR=-1 IF N NOT A POSITIVE INTEGER OR IF N
C EXCEEDS PROGRAMMED LIMIT, IERR=0 IF SOLUTIONS CARRIED TO COMPLETION,
C ELSE IERR=ROW FOR WHICH FIRST NON-+OSITIVE DIAGONAL PIVOT DETECTED.
      DIMENSION A(27),B(6)
      IF (N.GE.1.AND.N.LE.6) GO TO 1
      IERR=-1
      RETURN
    1 I2=N+1
      DO 35 I=1,N
      J=N
      K=I
      SUM=0.0
      IL=I-1
      IF (IL.LT.1) GO TO 15
      DO 12 L=1,IL
      FAC=A(K)
      SUM=SUM-FAC*FAC
      K=K+J
   12 J=J-1
   15 SUM=SUM+A(K)
      IF (SUM.GT.0.0) GO TO 20
      IERR=I
      RETURN
   20 FAC=1.0/SQRT(SUM)
      A(K)=FAC
      I1=I+1
      DO 35 J=I1,I2
      SUM=0.0
      K=N
      IJ=J
      II=I
      IF (IL.LT.1) GO TO 30
      DO 25 L=1,IL
      AI=A(II)
      SUM=SUM-AI*A(IJ)
      II=II+K
      IJ=IJ+K
   25 K=K-1
   30 SUM=SUM+A(IJ)
   35 A(IJ)=FAC*SUM
      J=((N+3)*N)/2
      I=N
   42 IF(I.LT.1) GOTO 998
      SUM=A(J)
      J=J-1
      K=I+1
      L=N
   48 IF(L.LT.K) GOTO 55
      SUM=SUM-A(J)*B(L)
      J=J-1
      L=L-1
      GO TO 48
   55 B(I)=A(J)*SUM
      J=J-1
      I=I-1
      GO TO 42
  998 IERR=0
      END
      SUBROUTINE RCMPD
      integer mxl, mxr
      parameter (mxl=500, mxr=2000)
C COMPUTES D SPACING FOR GENERATED REFLECTIONS
      COMMON/BLKI/Z(6)
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/BLKB/H(3),U(3),X(3),W(3),AD(3),ALS(3),CSA(3),SNA(3),CR(6)
      COMMON/BLKE/JTYP,HH(mxr,3),DCM(mxr),IOP(27),JOP(27),NOPT,DMIN,NR
      DO 5 I=1,3
      J=MOD(I,3)+1
      K=MOD(J,3)+1
      Z(I)=X(I)**2
    5 Z(I+3)=X(I)*X(J)*U(K)
      DO 20 IC=1,NR
      DO 10 J=1,3
   10 H(J)=HH(IC,J)
      CALL QCOMP
      IF (Q .GT. 0.0) GO TO 16
      DCM(IC)=0.0
      GO TO 20
   16 DCM(IC)=1.0/SQRT(Q)
   20 CONTINUE
      END
      SUBROUTINE COORD
C COMPUTES CORRECTIONS TO RECIPROCAL CELL ELEMENTS
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/BLKA/IDG(3),FMN(3),AO(3),UO(3),ALO(3),AR(3),A(27),B(6)
      COMMON/BLKB/H(3),U(3),X(3),W(3),AD(3),ALS(3),CSA(3),SNA(3),CR(6)
      GO TO (10,20,30,40,50,60,70),ITYP
   10 DO 15 J=1,3
      CR(J)=B(1)
   15 CR(J+3)=0.0
      RETURN
   20 DO 25 J=1,2
   25 CR(J)=B(1)
      CR(3)=B(2)
      DO 28 J=4,6
   28 CR(J)=0.0
      RETURN
   30 DO 35 J=1,3
      CR(J)=B(J)
   35 CR(J+3)=0.0
      RETURN
   40 DO 45 J=1,3
   45 CR(J)=B(J)
      CR(4)=0.0
      CR(5)=B(4)
      CR(6)=0.0
      RETURN
   50 DO 52 J=1,2
   52 CR(J)=B(1)
      CR(3)=B(2)
      DO 55 J=4,6
   55 CR(J)=0.0
      RETURN
   60 DO 65 J=1,6
   65 CR(J)=B(J)
      RETURN
   70 DO 75 J=1,3
      CR(J)=B(1)
   75 CR(J+3)=B(2)
      END
      SUBROUTINE CMPNUP
C COMPUTES DIRECT CELL PARAMETERS FROM CORRECTED RECIPROCAL CELL
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/BLKA/IDG(3),FMN(3),AO(3),UO(3),ALO(3),AR(3),A(27),B(6)
      COMMON/BLKB/H(3),U(3),X(3),W(3),AD(3),ALS(3),CSA(3),SNA(3),CR(6)
      DO 9 J=1,3
      UO(J)=AR(J)
      W(J)=SQRT(1.0-U(J)**2)
      AR(J)=3437.7468*(1.5707963-ATAN2(U(J),W(J)))
      CALL MTDM(AR(J),IDG(J),FMN(J))
    9 UO(J)=AR(J)-UO(J)
      VRO=VR
      VR=SQRT(1.0-U(1)**2-U(2)**2-U(3)**2+2.0*U(1)*U(2)*U(3))*X(1)*X(2)
     1*X(3)
      VRD=VR-VRO
      WRITE(6,12) (X(I),I=1,3),(IDG(I),FMN(I),I=1,3),VR
   12 FORMAT(16H RECIPROCAL CELL,1X,3E16.8,1X,3(I5,F7.3),E18.8)
      WRITE(6,14) (CR(J),J=1,3),(UO(J),J=1,3),VRD
   14 FORMAT(16H R C CORRECTIONS,1X,3E16.8,1X,3F12.3,E18.8)
      DO 20 I=1,3
      J=MOD(I,3)+1
      K=MOD(J,3)+1
      AO(I)=AD(I)
      AD(I)=X(J)*X(K)*W(I)/VR
      ALO(I)=ALS(I)
      IF (ITYP .GT. 5) GO TO 18
      IF (I .EQ. 2 .AND. ITYP .EQ. 4) GO TO 18
      IF (I .EQ. 3 .AND. ITYP .EQ. 2) GO TO 18
      CSA(I)=0.0
      GO TO 19
   18 CSA(I)=(U(J)*U(K)-U(I))/(W(J)*W(K))
   19 ALS(I)=3437.7468*(1.5707963-ATAN2(CSA(I),SQRT(1.0-CSA(I)**2)))
   20 CALL MTDM(ALS(I),IDG(I),FMN(I))
      VCO=VC
      VC=1.0/VR
      WRITE(6,24) (AD(I),I=1,3),(IDG(I),FMN(I),I=1,3),VC
   24 FORMAT(16H     DIRECT CELL,1X,3E16.8,1X,3(I5,F7.3),E18.8)
      DO 28 L=1,3
      FMN(L)=AD(L)-AO(L)
   28 UO(L)=ALS(L)-ALO(L)
      VCD=VC-VCO
      WRITE(6,30) (FMN(L),L=1,3),(UO(L),L=1,3),VCD
   30 FORMAT(16H D C CORRECTIONS,1X,3E16.8,1X,3F12.3,E18.8)
      END
      SUBROUTINE RSDPRN
C PRINTS OUT MAXIMUM RESIDUAL AND STANDARD ERROR OF UNIT WEIGHT
C OBSERVATION IN TERMS OT THETA OR 2-THETA ACCORDING TO ITHT
      COMMON/BLKC/THCMP,D,CONST,SWT,DIFF,JA,DOBS,PWL,ITHT,THTSE,THOBS,
     1 RSDMAX,IRSD,MMN,E(7),SUM
      RM=RSDMAX
      TSE=THTSE
      IF (ITHT.NE.0) GO TO 9
      RM=RM+RM
      TSE=TSE+TSE
    9 WRITE(6,10) RM,IRSD,TSE,MMN
   10 FORMAT(40H LARGEST RESIDUAL REDUCED TO UNIT WEIGHT,F8.5,3X,
     1 3HOBS,I4,3X,26HSTANDARD ERROR UNIT WT OBS,F8.5,3X,
     2 18HDEGREES OF FREEDOM,I4)
      END
      SUBROUTINE HKLPR
C PRINTS H,K,L,CALCULATED AND OBSERVED D-S AND THETAS
      DIMENSION IH(3)
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/BLKB/H(3),U(3),X(3),W(3),AD(3),ALS(3),CSA(3),SNA(3),CR(6)
      COMMON/BLKC/THCMP,D,CONST,SWT,DIFF,JA,DOBS,PWL,ITHT,THTSE,THOBS,
     1 RSDMAX,IRSD,MMN,E(7),SUM
      WT=SWT*SWT
      TC=THCMP
      OB=THOBS
      TD=DIFF
      IF (ITHT.NE.0) GO TO 12
      TC=TC+TC
      OB=OB+OB
      TD=TD+TD
   12 DO 13 L=1,3
   13 IH(L)=IFIX(H(L)+SIGN(0.5,H(L)))
      WRITE(6,15) M,(IH(L),L=1,3),JA,D,DOBS,PWL,TC,OB,TD,WT
   15 FORMAT(I4,3I5,A2,F11.6,F13.6,F12.6,F16.5,F15.5,F16.5,F14.5)
      END
      SUBROUTINE HKLGEN
C GENERATES MILLER INDICES CORRESPONDING TO PLANES AT DISTANCES FROM
C DMIN ON UP - MAXIMUM OF mxr REFLECTIONS MAY BE GENERATED
      integer mxl, mxr
      parameter (mxl=500, mxr=2000)
      COMMON/ADATA/IALF(19),IFI(2,27),IFJ(3,27),IHKL(6,3),JHKL(10,3)
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/BLKA/IDG(3),FMN(3),AO(3),UO(3),ALO(3),AR(3),A(27),B(6)
      COMMON/BLKB/H(3),U(3),X(3),W(3),AD(3),ALS(3),CSA(3),SNA(3),CR(6)
      COMMON/BLKC/THCMP,D,CONST,SWT,DIFF,JA,DOBS,PWL,ITHT,THTSE,THOBS,
     1 RSDMAX,IRSD,MMN,E(7),SUM
      COMMON/BLKE/JTYP,HH(mxr,3),DCM(mxr),IOP(27),JOP(27),NOPT,DMIN,NR
      COMMON/BLKI/Z(6)
      QMAX=1.0/(DMIN*DMIN)
      NR=0
      DO 10 I=1,3
      J=MOD(I,3)+1
      K=MOD(J,3)+1
      Z(I)=X(I)**2
   10 Z(I+3)=X(I)*X(J)*U(K)
      E(1)=Z(1)
      E(2)=Z(2)
      E(3)=Z(4)
      TD=E(1)*E(2)-E(3)*E(3)
      IL=0
   15 H(3)=IL
      E(4)=(Z(5)+Z(5))*H(3)
      E(5)=(Z(6)+Z(6))*H(3)
      E(6)=Z(3)*H(3)*H(3)-QMAX
      TN=E(3)*E(5)-E(1)*E(4)
      DSCRS=TN*TN-(4.0*E(1)*E(6)-E(5)*E(5))*TD
      IF (DSCRS.LT.0.0) RETURN
      DSCRM=SQRT(DSCRS)
      IK1=IFIX((TN-DSCRM)*0.5/TD)
      IK2=IFIX((TN+DSCRM)*0.5/TD)
      GO TO (25,27,27,27,27,29,29,27),JTYP
   25 IF (IK1.LT.   IL) IK1 =   IL
      GO TO 29
   27 IF (IK1.LT.0) IK1=0
   29 IK=IK2
      IF (IK .LT. IK1) GO TO 58
   31 H(2)=IK
      GT=(E(3)+E(3))*H(2)+E(5)
      HT=(Z(5)*H(3)+Z(2)*H(2))*H(2)+(Z(5)*H(2)+Z(3)*H(3))*H(3)-QMAX
      DEN=Z(1)+Z(1)
      DSCRS=GT*GT-(DEN+DEN)*HT
      IF (DSCRS .LT. 0.0) GO TO 57
      DSCRM=SQRT(DSCRS)
      IH1=IFIX(-(GT+DSCRM)/DEN)
      IH2=IFIX((DSCRM-GT)/DEN)
      GO TO (42,42,40,44,42,44,44,38),JTYP
   38 IF(MOD(IL,3).EQ.0) GO TO 42
   40 IF (IH1.LT.0) IH1=0
      GO TO 44
   42 IF (IH1.LT.IK) IH1=IK
   44 IH=IH2
      IF (IH .LT. IH1) GO TO 57
   46 H(1)=IH
      IF (NOPT.LT.1)GO TO 50
      CALL CHKEXT(IH,IK,IL,JEXT)
      IF(JEXT.NE.0) GO TO 56
   50 CALL QCOMP
      IF (Q.LE.0.0) GO TO 56
      NR=NR+1
      DCM(NR)=1.0/SQRT(Q)
      DO 52 J=1,3
   52 HH(NR,J)=H(J)
      IF (NR.LT.mxr) GO TO 56
      CALL SORTD(NR,DCM,HH,FMN)
      CALL SKWEZE
      IF (NR.LT.mxr) GO TO 56
c     WRITE(6,54) DMIN
c  54 FORMAT(43H0MORE THAN mxr REFLECTIONS EXIST FOR DMIN =,F10.5)
      write(6, '(a, i4, a, f10.5)')
     & '0MORE THAN ', mxr, ' REFLECTIONS EXIST FOR DMIN =', dmin
      stop
c     RETURN
   56 IH=IH-1
      IF(IH.GE.IH1)GO TO 46
   57 IK=IK-1
      IF (IK.GE.IK1) GO TO 31
   58 IL=IL+1
      GO TO 15
      END
      SUBROUTINE SORTD(N,D,HH,H)
C SORTS ARRAY D SO THAT ELEMENTS ARE IN DECREASING ORDER - MILLER
C INDICES IN ARRAY HH MADE TO CORRESPOND
      integer mxl, mxr
      parameter (mxl=500, mxr=2000)
      DIMENSION D (mxr),HH(mxr,3),H(3)
      I=N
    2 I=I-1
      IF (I.LT.1) RETURN
      DO 10 J=1,I
      IF(D(J).GE.D(J+1)) GO TO 10
      FAC=D(J)
      D(J)=D(J+1)
      D(J+1)=FAC
      DO 9 K=1,3
      H(K)=HH(J,K)
      HH(J,K)=HH(J+1,K)
    9 HH(J+1,K)=H(K)
   10 CONTINUE
      GO TO 2
      END
      SUBROUTINE DLQOT
C COMPUTES PARTIAL DERIVATIVES OF RECIPROCAL OF D SQUARED
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/BLKA/IDG(3),FMN(3),AO(3),UO(3),ALO(3),AR(3),A(27),B(6)
      COMMON/BLKB/H(3),U(3),X(3),W(3),AD(3),ALS(3),CSA(3),SNA(3),CR(6)
      DO 15 J=1,3
      K=MOD(J,3)+1
      L=MOD(K,3)+1
      B(J)=FAC*(X(J)*H(J)**2+H(J)*H(L)*X(L)*U(K)+H(J)*H(K)*X(K)*U(L))
   15 B(J+3)=FAC*H(K)*H(L)*X(K)*X(L)
      END
      FUNCTION THETA(C,D)
C COMPUTES ARCSIN(C/D) PROVIDED C.LE.D. ELSE THETA SET TO 90.00
      FAC=D*D-C*C
      IF (FAC.GE.0.0 .AND. (D*D+C*C) .GT. 0.0) GO TO 7
      FAC=90.0
      GO TO 8
    7 FAC=57.29578*ATAN2(C,SQRT(FAC))
    8 THETA=FAC
      END
      SUBROUTINE OBSEQ
C FORMS OBSERVATION EQUATION FOR ACCEPTED REFLECTIONS
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/BLKA/IDG(3),FMN(3),AO(3),UO(3),ALO(3),AR(3),A(27),B(6)
      COMMON/BLKB/H(3),U(3),X(3),W(3),AD(3),ALS(3),CSA(3),SNA(3),CR(6)
      COMMON/BLKC/THCMP,D,CONST,SWT,DIFF,JA,DOBS,PWL,ITHT,THTSE,THOBS,
     1 RSDMAX,IRSD,MMN,E(7),SUM
      M=M+1
      DIV=COS(0.017453293*THCMP)
      IF(DIV.GT.0.) GOTO 8
      M=M-1
      RETURN
    8 FAC=(57.29578*D*CONST*SWT)/DIV
      EQCNS=SWT*DIFF
      CALL DLQOT
      E(N+1)=-EQCNS
      GO TO (10,20,30,40,50,60,70),ITYP
   10 E(1)=B(1)+B(2)+B(3)
      RETURN
   20 E(1)=B(1)+B(2)
      E(2)=B(3)
      RETURN
   30 DO 35 J=1,3
   35 E(J)=B(J)
      RETURN
   40 DO 45 J=1,3
   45 E(J)=B(J)
      E(4)=B(5)
      RETURN
   50 E(1)=B(1)+B(2)
      E(2)=B(3)
      RETURN
   60 DO 65 J=1,6
   65 E(J)=B(J)
      RETURN
   70 DO 75 J=1,2
   75 E(J)=B(3*J-2)+B(3*J-1)+B(3*J)
      END
      SUBROUTINE NORME
C ACCUMULATES NORMAL EQUATIONS AND SUM OF SQUARES
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/BLKA/IDG(3),FMN(3),AO(3),UO(3),ALO(3),AR(3),A(27),B(6)
      COMMON/BLKC/THCMP,D,CONST,SWT,DIFF,JA,DOBS,PWL,ITHT,THTSE,THOBS,
     1 RSDMAX,IRSD,MMN,E(7),SUM
      EQUIVALENCE (RSMX,RSDMAX),(IRSD,MAXM)
      L=0
      NP1=N+1
      DO 20 J=1,N
      DO 20 K=J,NP1
      L=L+1
   20 A(L)=A(L)+E(J)*E(K)
      SUM=SUM+E(NP1)**2
      IF (ABS(E(NP1)).LE.ABS(RSMX)) RETURN
      MAXM=M
      RSMX=E(NP1)
      END
      SUBROUTINE MTDM(AMN,IDGS,FMNS)
C CONVERTS MINUTES TO DEGREES AND MINUTES
      IDGS=AMN/60.0
      FMNS=ABS(AMN-60.0*FLOAT(IDGS))
      IF(FMNS.LT.59.9995) RETURN
      IDGS=IDGS+1
      FMNS=0.0
      END
      SUBROUTINE PDNVRT(A,N,B)
C COMPLETES THE INVERSION OF POSITIVE DEFINITE SYMMETRIC MATRIX OF
C ORDER N,STARTING WITH THE TRIANGULARIZED MATRIX RESULTING FROM
C PDSQRT, THE UPPER TRIANGLE OF THE INVERSE BEING STORED SO THAT ITS
C ITH ROW AND JTH COLUMN ELEMENT IS A((2*N+2-I)*(I-1)/2+J-1)
C WITH 1.LE.I.LE.J.LE.N.
      DIMENSION A(27),B(6)
      IE=0
      NP1=N+1
      NP2=N+2
      INM=((N+3)*N)/2-1
      II=1
      DO 395 I=1,N
      B(I)=A(II)
      J=I+1
  155 IF (J.GT.N) GO TO 200
      SUM=0.0
      IJ=II+J-I
      L=J-1
      DO 170 K=I,L
      SUM=SUM-A(IJ)*B(K)
  170 IJ=IJ+NP1-K
      B(J)=A(IJ)*SUM
      J=J+1
      GO TO 155
  200 J=N
      IJ=INM
  210 SUM=B(J)
      K=N
  220 IF (K.LE.J) GO TO 250
      SUM=SUM-A(IJ)*B(K)
      IJ=IJ-1
      K=K-1
      GO TO 220
  250 B(J)=A(IJ)*SUM
      J=J-1
      IF (J.LT.I) GO TO 300
      IJ=IJ-2
      GO TO 210
  300 DO 315 J=I,N
      IE=IE+1
  315 A(IE)=B(J)
  395 II=II+NP2-I
      END
      SUBROUTINE SESUB
C COMPUTES STANDARD ERROR AND VARIANCE - COVARIANCE MATRIX
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/BLKA/IDG(3),FMN(3),AO(3),UO(3),ALO(3),AR(3),A(27),B(6)
      COMMON/BLKB/H(3),U(3),X(3),W(3),AD(3),ALS(3),CSA(3),SNA(3),CR(6)
      COMMON/BLKC/THCMP,D,CONST,SWT,DIFF,JA,DOBS,PWL,ITHT,THTSE,THOBS,
     1 RSDMAX,IRSD,MMN,E(7),SUM
      COMMON/BLKD/TT(6,6),ST(6,6),SERCE(6),SEDCE(6),DD(3),VCME(6,6),
     1 SS(21)
      DO 2 I=1,21
    2 SS(I)=0.0
      FACTR=THTSE*THTSE
      CALL PDNVRT(A,N,B)
      IEC=((N+1)*N)/2
      DO 8 I=1,IEC
    8 A(I)=A(I)*FACTR
      GO TO (10,20,30,40,50,60,70),ITYP
   10 DO 11 I=1,3
   11 SS(I)=A(1)
      SS(7)=A(1)
      SS(8)=A(1)
      SS(12)=A(1)
      GO TO 80
   20 SS(1)=A(1)
      SS(2)=A(1)
      SS(7)=A(1)
      SS(3)=A(2)
      SS(8)=A(2)
      SS(12)=A(3)
      GO TO 80
   30 DO 32 I=1,3
   32 SS(I)=A(I)
      DO 34 I=7,8
   34 SS(I)=A(I-3)
      SS(12)=A(6)
      GO TO 80
   40 DO 42 I=1,3
   42 SS(I)=A(I)
      SS(5)=A(4)
      DO 44 I=7,8
   44 SS(I)=A(I-2)
      DO 46 I=7,9
   46 SS(2*I-4)=A(I)
      SS(19)=A(10)
      GO TO 80
   50 GO TO 20
   60 DO 62 I=1,21
   62 SS(I)=A(I)
      GO TO 80
   70 DO 72 I=1,3
   72 SS(I)=A(1)
      SS(7)=A(1)
      SS(8)=A(1)
      SS(12)=A(1)
      DO 74 I=4,6
   74 SS(I)=A(2)
      DO 76 I=9,11
   76 SS(I)=A(2)
      DO 78 I=13,15
   78 SS(I)=A(2)
      DO 79 I=16,21
   79 SS(I)=A(3)
   80 II=1
      DO 82 I=1,3
      SERCE(I)=SQRT(SS(II))
   82 II=II+7-I
      DO 84 I=4,6
      SERCE(I)=3437.7468*SQRT(SS(II))/W(I-3)
   84 II=II+7-I
      WRITE(6,86) (SERCE(I),I=1,6)
   86 FORMAT(16H R C STNDRD ERRS,1X,3E16.8,1X,3F12.3)
      DO 88 I=1,3
   88 SNA(I)=SQRT(1.0-CSA(I)**2)
      DO  95 I=1,3
      DO  95 J=1,3
      TT(I,J)=0.0
   95 IF(I.EQ.J) TT(I,J)=-AD(I)/X(I)
      DO 100 I=4,6
      DO 100 J=1,3
  100 TT(I,J)=0.0
      DO 105 I=1,3
      J=MOD(I,3)+1
      K=MOD(J,3)+1
  105 DD(I)=-AD(J)*AD(K)*X(J)*X(K)
      DO 915 I=1,3
      DO 115 J=4,6
      IF (I-J+3) 108,110,108
  108 TT(I,J)=DD(J-3)*AD(I)*CSA(J-3)
      GO TO 115
  110 K=MOD(J,3)+1
      L=MOD(K,3)+1
      TT(I,J)=DD(J-3)*AD(J-3)*CSA(K)*CSA(L)
  115 CONTINUE
  915 CONTINUE
      DO 120 I=1,3
      J=MOD(I,3)+1
      K=MOD(J,3)+1
  120 DD(I)=1.0/(W(J)*W(K))
      DO 925 I=4,6
      DO 125 J=4,6
      IF(I.NE.J) GOTO 123
      TT(I,J)=DD(I-3)/SNA(J-3)
      GO TO 125
  123 IJ12=12-I-J
      TT(I,J)=DD(I-3)*CSA(IJ12)/SNA(J-3)
  125 CONTINUE
  925 CONTINUE
      DO 130 I=1,6
      DO 130 J=1,6
      SUM=0.0
      IJ=I
      DO 128 K=1,6
      SUM=SUM+SS(IJ)*TT(J,K)
      IJ=IJ+1
  128 IF(I.GT.K) IJ=IJ+5-K
  130 ST(I,J)=SUM
      WRITE(6,132)
  132 FORMAT(39H DIRECT CELL VARIANCE-COVARIANCE MATRIX,75X,3HROW)
      DO 140 I=1,6
      DO 138 J=1,6
      SUM=0.0
      DO 136 K=1,6
  136 SUM=SUM+TT(I,K)*ST(K,J)
  138 VCME(I,J)=SUM
  140 WRITE(6,142) (VCME(I,L),L=1,6),I
  142 FORMAT(16X,6E16.8,I4)
      DO 145 I=1,6
  145 SEDCE(I)=SQRT(VCME(I,I))
      DO 148 I=4,6
  148 SEDCE(I)=3437.7468*SEDCE(I)
      FACTR=AD(1)*AD(2)*AD(3)
      FACTR=FACTR*FACTR/VC
      DO 150 I=1,3
      J=MOD(I,3)+1
      K=MOD(J,3)+1
      CR(I)=VC/AD(I)
  150 CR(I+3)=FACTR*SNA(I)*(CSA(I)-CSA(J)*CSA(K))
      FACTR=0.0
      DO 152 I=1,6
      DO 152 J=1,6
  152 FACTR=FACTR+CR(I)*CR(J)*VCME(I,J)
      IF(FACTR.LT.0.) FACTR=0.
      FACTR=SQRT(FACTR)
      WRITE(6,160) (SEDCE(I),I=1,6),FACTR
  160 FORMAT(16H D C STNDRD ERRS,1X,3E16.8,1X,3F12.3,E18.8)
      END
      SUBROUTINE CHKEXT(IH,IK,IL,JEXT)
C CHECKS MILLER INDICES IH,IK,IL FOR EXTINCTIONS - JEXT SET TO 1
C IF CURRENT REFLECTION NOT ACCEPTABLE, JEXT SET TO 0 IF IT IS
      integer mxl, mxr
      parameter (mxl=500, mxr=2000)
      COMMON/ADATA/IALF(19),IFI(2,27),IFJ(3,27),IHKL(6,3),JHKL(10,3)
      COMMON/BLKE/JTYP,HH(mxr,3),DCM(mxr),IOP(27),JOP(27),NOPT,DMIN,NR
      JEXT=0
      DO 22 I=1,NOPT
      L=JOP(I)
      IDIV=IOP(L)
      DO 10 J=1,2
      IFIJ=IFI(J,L)
      IF (IFIJ.LT.1.OR.IFIJ.GT.6) GO TO 12
      IANS=IH*IHKL(IFIJ,1)+IK*IHKL(IFIJ,2)+IL*IHKL(IFIJ,3)
      IF (IANS.NE.0) GO TO 22
   10 CONTINUE
   12 DO 20 K=1,3
      IFJK=IFJ(K,L)
      IF (IFJK.LT.1.OR.IFJK.GT.10) GO TO 22
      IARG=JHKL(IFJK,1)*IH+JHKL(IFJK,2)*IK+JHKL(IFJK,3)*IL
      IF (MOD(IARG,IDIV).EQ.0) GO TO 20
      JEXT=1
      RETURN
   20 CONTINUE
   22 CONTINUE
      END
      SUBROUTINE TSTHKL(ITYP,IH,JH,IEQL)
C TESTS FOR EXACT EQUALITY OF D-SPACING FOR TWO REFLECTIONS
      DIMENSION IH(3),JH(3),ITEST(6,7),IHKLC(6,10)
      DATA ITEST/8,5*0, 3,7,4*0, 1,2,3,3*0, 1,2,3,6,2*0, 3,10,4*0, 1,2,3
     *,4,5,6, 8,9,4*0/, IHKLC/2*1,4*0, 2*2,4*0, 2*3,4*0, 1,2,4*0, 2,3,4*
     *0, 3,1,4*0, 3*1,3*2, 2*1,2*2,2*3, 1,2*2,2*3,1, 2*1,2*2,2*0/
      DO 30 I=1,6
      INTRY=ITEST(I,ITYP)
      IF (INTRY.LT.1.OR.INTRY.GT.10) GO TO 32
      ISUM=0
      J=1
    5 K=IHKLC(J,INTRY)
      IF (K.LT.1.OR.K.GT.3) GO TO 20
      L=IHKLC(J+1,INTRY)
      ISUM=ISUM+IH(K)*IH(L)-JH(K)*JH(L)
      J=J+2
      IF (J.LT.6) GO TO 5
   20 IF (ISUM.EQ.0) GO TO 30
      IEQL=1
      RETURN
   30 CONTINUE
   32 IEQL=0
      END
      SUBROUTINE SKWEZE
C ELIMINATES REFLECTIONS WITH IDENTICAL D-SPACINGS
      DIMENSION IH(3),JH(3)
      integer mxl, mxr
      parameter (mxl=500, mxr=2000)
      COMMON ICRYST,ITYP,N,M,VR,VRO,VC,VCO,Q,FAC
      COMMON/BLKE/JTYP,HH(mxr,3),DCM(mxr),IOP(27),JOP(27),NOPT,DMIN,NR
      TOL=1.0E-5
      I=1
   12 J=I+1
      IF (J.GT.NR) RETURN
      IF ((DCM(I)-DCM(J)).GE.TOL) GO TO 40
      DO 15 K=1,3
      IH(K)=IFIX(HH(I,K)+SIGN(0.5,HH(I,K)))
   15 JH(K)=IFIX(HH(J,K)+SIGN(0.5,HH(J,K)))
   19 CALL TSTHKL(ITYP,IH,JH,IEQL)
      IF (IEQL.NE.0) GO TO 27
      J=J+1
      IF (J.GT.NR) GO TO 27
      IF ((DCM(I)-DCM(J)).GE.TOL) GO TO 27
      DO 25 K=1,3
   25 JH(K)=IFIX(HH(J,K)+SIGN(0.5,HH(J,K)))
      GO TO 19
   27 JMI=J-I-1
      J2=NR-JMI
      J1=I+1
      IF (J2 .LT. J1) GO TO 35
      DO 32 K=J1,J2
      K2= K+JMI
      DCM(K)=DCM(K2)
      DO 32 L=1,3
   32 HH(K,L)=HH(K2,L)
   35 NR=J2
   40 I=I+1
      GO TO 12
      END
      BLOCK DATA
      COMMON /ADATA/IALF(19),IFI(2,27),IFJ(3,27),IHKL(6,3),JHKL(10,3)
      DATA IALF/4HCUBI,4HTETR,4HORTH,4HMONO,4HRHOM,4HTRIC,4HHXAG,4HHXAR,
     *4HTITL,4HPARA,4HSYSE,3HEND,4HFINI,2H2-,2H  ,2H R,2H F,2H  ,2H C/,
     *IFI/12*0,1,0,1,0,1,0,1,0,4,0,4,0,4,0,5,0,5,0,5,0,6,0,6,0,6,0,1,6,5
     *,6,4,6,4,5,2,0,2,0,2,0,2,0/,IFJ/1,0,0,2,0,0,3,0,0,1,2,3,4,0,0,5,0,
     *0,6,0,0,8,0,0,2,0,0,9,0,0,7,0,0,8,0,0,3,0,0,6,0,0,8,0,0,2,0,0,6,0,
     *0,7,0,0,1,0,0,6,0,0,6,0,0,7,0,0,8,0,0,6,0,0,8,0,0,2,0,0,10,0,0/,
     *IHKL/1,0,1,1,0,0,-1,1,0,0,1,0,0,-1,-1,0,0,1/, JHKL/1,1,0,1,-1,1,0,
     *0,2,1,1,0,1,1,1,0,1,0,0,0,0,1,1,1,1,0,0,1,1,2/
      END
