CDC   PROGRAM LOUM(INPUT,OUTPUT,TEMP,REFINE,TAPE5=INPUT,TAPE6=OUTPUT,
CDC  1 TAPE8=TEMP,TAPE9=REFINE)
C     PROGRAMME PRINCIPAL
C     VERSION  20 OCTOBRE 1977  LOUM-9
C     MONOCLINIQUE  PARTITION  ESPACE VOLUME   -  DICHOTOMIE
C---MODIFIED (ULCC VERSION ONLY) TO WRITE SOLUTIONS TO REFINE FILE
C---(CHANNEL 9) IN FORMAT SUITABLE FOR REFINEMENT BY FZON(6)
CDC   DIMENSION IDENT(30)
	COMMON D(50),EPSIL(50),Q(50),EPSQ(50),KQ(50),KEPSQ(50),KQT(50),N,V
     1,MH,MH2,MK,MK2,ML,ML2,PAS,PAS2,PAS4,PAS8,PAS16,PAS32,PAS64,CONST,V
     2INF,VSUP
C---USED TO PASS TITLE TO S/R REFSET FOR REFORMATTING (ULCC)
	COMMON /IDTCOM/ IDENT(30)
	READ(5,1)(IDENT(K),K=1,30)
	READ(5,2)N
	READ(5,3)AMAX,BMAX,CMAX,VOLMIN,VOLMAX,BEMIN,BEMAX
	READ(5,4)(D(I),EPSIL(I),I=1,N)
	IF(BEMIN.EQ.0.)BEMIN=90.
	IF(BEMAX.EQ.0.)BEMAX=130.
	WRITE(6,5)(IDENT(K),K=1,30)
	WRITE(6,6)
	WRITE(6,14)(D(I),EPSIL(I),I=1,N)
	WRITE(6,7)AMAX,BMAX,CMAX,BEMIN,BEMAX,VOLMIN,VOLMAX
1	FORMAT(30A2)
2	FORMAT(I3)
3	FORMAT(7F7.2)
4	FORMAT(2F7.4)
5	FORMAT(1H1/1H0,30X,30A2/1H ,30X,60(1H*))
6     FORMAT(1H0/1H0,45X,30HDONNEES FOURNIES A L EXECUTION/1H ,45X,30(1H    
     1*)/1H0,20X,37HDISTANCES RETICULAIRES EXPERIMENTALES,20X,20HERREUR    
     2EXPERIMENTALE)
7     FORMAT(//1H0,55X,10(1H*)////1H0,10X,11HA MAXIMUM =,F6.2/1H0,10X,11
     1HB MAXIMUM =,F6.2/1H0,10X,11HC MAXIMUM =,F6.2/1H0,10X,14HBETA MINI
     2MUM =,F6.2/1H0,10X,14HBETA MAXIMUM =,F6.2/1H0,10X,16HVOLUME MINIMU
     3M =,F9.2/1H0,10X,16HVOLUME MAXIMUM =,F9.2///1H0,55X,10(1H*)/)
8	FORMAT(1H1/1H0,50X,15HSYSTEME CUBIQUE/1H ,50X,15(1H*)///)
9	FORMAT(////1H0,50X,32HPAS D INDEXATION DANS CE DOMAINE///)
10	FORMAT(1H1/1H0,50X,18HSYSTEME TETRAGONAL/1H ,50X,18(1H*)///)
11	FORMAT(1H1/1H0,50X,17HSYSTEME HEXAGONAL/1H ,50X,17(1H*)///)
12	FORMAT(1H1/1H0,50X,22HSYSTEME ORTHORHOMBIQUE/1H ,50X,22(1H*)///)
C---REFINEMENT INITIALISATION (ULCC)
	CALL REFSET
13	FORMAT(1H0,50X,3HFIN/1H ,50X,3(1H*))
14	FORMAT(1H ,33X,F10.4,36X,F10.4)
15	FORMAT(1H ,5X,22HSYSTEME MONOCLINIQUE :)
16	FORMAT(3F10.4)
70    FORMAT(1H ,10X,44HDOMAINE VOLUME EXAMINE(TOLERANCE SUPERIEURE)/1H
     1,10X,44(1H*)/1H0,20X,15HVOLUME MINIMUM=,F9.2,4HA**3,10X,15HVOLUME
     2MAXIMUM=,F9.2,4HA**3//)
	CONST=1.0E08
	DO 100 I=1,N
	PR=D(I)*D(I)
	Q(I)=1./PR
	KQ(I)=CONST*Q(I)
	EPSQ(I)=(2.*EPSIL(I))/(PR*D(I))
100	KEPSQ(I)=CONST*EPSQ(I)
	PAS=0.5
	PAS2=PAS/2.
	PAS4=PAS/4.
	PAS8=PAS/8.
	PAS16=PAS/16.
	PAS32=PAS/32.
	PAS64=PAS/64.
	PI=3.1415926/180.
	BEMIN=BEMIN*PI
	BEMAX=BEMAX*PI
	AMIN=D(1)-EPSIL(1)
	ZB=D(1)/D(2)
	ZZ=ABS(2.-ZB)
	IF (ZZ-5.0E-04)19,18,18
19	BMIN=D(3)-EPSIL(3)
	GO TO 17
18	BMIN=D(2)-EPSIL(2)
17	AMAX=AMAX-PAS2
	BMAX=BMAX-PAS2
	CMAX=CMAX-PAS2
	CMIN=2.5
	NA=(AMAX-BMIN)/PAS+1.
	NB=(BMAX-CMIN)/PAS+1.
	NC=(CMAX-CMIN)/PAS+1.
	V=VOLMAX
	VSUP=VOLMIN
	PASVOL=400.
	NVOL=(VOLMAX-VOLMIN)/PASVOL
	NVOL=NVOL+1
	DO 1000 KVOL=1,NVOL
	VV=V
	VINF=VSUP
	VSUP=VINF+PASVOL
	IF(VSUP.GT.VOLMAX)VSUP=VOLMAX
	IF(VINF.GE.VOLMAX)GO TO 1100
	WRITE(6,70)VINF,VSUP
	  VORT=V
51	 WRITE(6,15)
	CALL MONOC(BMIN,NA,NB,CMIN,NC,BEMAX,BEMIN)
	DIFVOL=ABS(V-VORT)
	IF(DIFVOL-0.02) 60,60,61
60	WRITE(6,9)
61	DIFVOL=ABS(V-VV)
	IF(DIFVOL-0.02)1000,1000,1100
1000	CONTINUE
1100	WRITE(6,13)
C---COPY SAVED SOLUTIONS TO REFINE FILE (ULCC)
	CALL REFSAV
C---READ SOLUTIONS, ETC. FROM REFINE FILE AND USE MODIFIED VERSION
C--- OF DUTCH ZONE-INDEXING PROGRAM FOR NORMALISATION, REFINEMENT
C--- AND EVALUATION (ULCC)
	CALL FZONRF
	STOP
	END
	
	SUBROUTINE MONOC(AMIN,NA,NB,CMIN,NC,BEMAX,BEMIN)
	DIMENSION QCAL(50,12),XH(50,12),XK(50,12),XL(50,12)
	COMMON D(50),EPSIL(50),Q(50),EPSQ(50),KQ(50),KEPSQ(50),KQT(50),N,V
     1,MH,MH2,MK,MK2,ML,ML2,PAS,PAS2,PAS4,PAS8,PAS16,PAS32,PAS64,CONST,V
     2INF,VSUP
C---QRCELL AND QCELL ARE FOR PASSING SOLUTIONS TO S/R REFINC (ULCC)
	DIMENSION QRCELL(6),QCELL(6)
	EQUIVALENCE (QAR,QRCELL(1)),(QBR,QRCELL(2)),(QCR,QRCELL(3)),
     1 (QDR,QRCELL(4)),(QER,QRCELL(5)),(QFR,QRCELL(6))
	DATA QDR,QFR/0.,0./
1     FORMAT(1H0,5X,2HA=,F8.5,5X,2HB=,F8.5,5X,2HC=,F8.5,5X,5HBETA=,F8.3,
     15X,33HVOLUME DE LA MAILLE ELEMENTAIRE =,F10.3/1H0,7X,14HD EXPERIME
     2NTAL,12X,1HH,12X,1HK,12X,1HL,10X,9HD CALCULE,12X,5HDE-DC,12X,7HEPS
     3ILON/1H0)
2     FORMAT(1H ,10X,F8.5,4X,3(10X,F3.0),9X,F8.5,11X,F8.5,9X,F8.5,10X,F8
     1.5)
3	FORMAT(1H0,50X,9(1H-)/)
4     FORMAT(1H ,20X,23HDOMAINE ANGLE EXAMINE :,5X,13HBETA MINIMUM=,F10.
     12,6HDEGRES,5X,13HBETA MAXIMUM=,F10.2,6HDEGRES)
	AA=100.
	BB=100.
	CC=100.
	BEBE=3.
	BMIN=CMIN
	PI=3.1415926
	PIDEG=PI/180.
	PASBE=PIDEG*2.5
	NBETA=((BEMAX-BEMIN)/PASBE)+1.
	PAB2=PASBE/2.
	PAB4=PASBE/4.
	PAB8=PASBE/8.
	PAB16=PASBE/16.
	PAB32=PASBE/32.
	PAB64=PASBE/64.
	DO 1000 IBET1=1,NBETA
	BETW=IBET1-1
	BEINF1=BEMIN+(BETW*PASBE)
	BESUP1=BEINF1+PASBE
	BBETM=BEINF1/PIDEG
	BBETP=BESUP1/PIDEG
	WRITE(6,4)BBETM,BBETP
	SM1=1./(SIN(BEINF1))**2
	SP1=1./(SIN(BESUP1))**2
	CSM1=COS(BEINF1)*SM1
	CSP1=COS(BESUP1)*SP1
	DO 1100 IA=1,NA
	AW=IA-1
	AMOI1=AMIN+(AW*PAS)
	APLU1=AMOI1+PAS
	APMO=APLU1*SIN(BEINF1)
	DO 21 I=1,10
	PP=I
	DA=APMO/PP
	IF(D(5)-DA)21,22,22
21	CONTINUE
22	MH=I
	DO 55 J=I,20
	PP=J
	DA=APMO/PP
	IF(D(N)-DA)55,51,51
55	CONTINUE
51	MH2=J
	KRAP1=CONST*SP1/(AMOI1*AMOI1)
	KRAM1=CONST*SM1/(APLU1*APLU1)
	DO 1200 IB=1,NB
	BW=IB-1
	BMOI1=BMIN+(BW*PAS)
	BPLU1=BMOI1+PAS
	IF((APLU1.LT.D(1)).AND.(BPLU1.LT.D(1)))GO TO 1200
	DO 23 I=1,10
	PP=I
	DB=BPLU1/PP
	IF(D(5)-DB)23,24,24
23	CONTINUE
24	MK=I
	DO 85 J=I,20
	PP=J
	DB=BPLU1/PP
	IF(D(N)-DB)85,61,61
85	CONTINUE
61	MK2=J
	KRBM1=CONST/(BPLU1*BPLU1)
	KRBP1=CONST/(BMOI1*BMOI1)
	DO 1300 IC=1,NC
	CW=IC-1
	CMOI1=CMIN+(CW*PAS)
	CPLU1=CMOI1+PAS
	IF(CMOI1.GT.AMOI1)GO TO 1200
	VTEST=AMOI1*CMOI1*BMOI1*SIN(BESUP1)
	IF(V-VTEST)1200,802,802
802	IF(VTEST.LT.VINF.OR.VTEST.GT.VSUP)GO TO 1300
	CPMO=CPLU1*SIN(BEINF1)
	DO 25 I=1,10
	PP=I
	DC=CPMO/PP
	IF(D(5)-DC)25,26,26
25	CONTINUE
26	ML=I
	DO 75 J=I,20
	PP=J
	DC=CPMO/PP
	IF(D(N)-DC)75,71,71
75	CONTINUE
71	ML2=J
	KRCM=CONST*SM1/(CPLU1*CPLU1)
	KRCP=CONST*SP1/(CMOI1*CMOI1)
	KRACM=2.*CONST*CSM1/(APLU1*CPLU1)
	KRACP=2.*CONST*CSP1/(AMOI1*CMOI1)
	CALL MODHKL(KRAM1,KRAP1,KRCM,KRCP,KRACM,KRACP,KRBM1,KRBP1)
	DO 120 I=6,N
	IF(KQT(I).LT.0)GO TO 1300
120	CONTINUE
	DO 2000 IBET2=1,2
	BETW=IBET2-1
	BEINF2=BEINF1+(BETW*PAB2)
	BESUP2=BEINF2+PAB2
	SM2=1./(SIN(BEINF2))**2
	SP2=1./(SIN(BESUP2))**2
	CSM2=COS(BEINF2)*SM2
	CSP2=COS(BESUP2)*SP2
	DO 2100 IA2=1,2
	AW=IA2-1
	AMOI2=AMOI1+(AW*PAS2)
	APLU2=AMOI2+PAS2
	KRAM2=CONST*SM2/(APLU2*APLU2)
	KRAP2=CONST*SP2/(AMOI2*AMOI2)
	DO 2200 IB2=1,2
	BW=IB2-1
	BMOI2=BMOI1+(BW*PAS2)
	BPLU2=BMOI2+PAS2
	KRBM2=CONST/(BPLU2*BPLU2)
	KRBP2=CONST/(BMOI2*BMOI2)
	DO 2300 IC2=1,2
	CW=IC2-1
	CMOI2=CMOI1+(CW*PAS2)
	CPLU2=CMOI2+PAS2
	ATEST=ABS(AMOI2-AA)
	BTEST=ABS(BMOI2-BB)
	CTEST=ABS(CMOI2-CC)
	BETEST=ABS(BEINF2-BEBE)
	IF(ATEST.LT.0.30.AND.BTEST.LT.0.30.AND.CTEST.LT.0.30.AND.BETEST.LT
     1.0.035)GO TO 2300
	KRCM=CONST*SM2/(CPLU2*CPLU2)
	KRCP=CONST*SP2/(CMOI2*CMOI2)
	KRACM=2.*CONST*CSM2/(APLU2*CPLU2)
	KRACP=2.*CONST*CSP2/(AMOI2*CMOI2)
	CALL MODHKL(KRAM2,KRAP2,KRCM,KRCP,KRACM,KRACP,KRBM2,KRBP2)
	DO 122 I=6,N
	IF(KQT(I).LT.0)GO TO 2300
122	CONTINUE
	DO 3000 IBET3=1,2
	BETW=IBET3-1
	BEINF3=BEINF2+(BETW*PAB4)
	BESUP3=BEINF3+PAB4
	SM3=1./(SIN(BEINF3))**2
	SP3=1./(SIN(BESUP3))**2
	CSM3=COS(BEINF3)*SM3
	CSP3=COS(BESUP3)*SP3
	DO 3100 IA3=1,2
	AW=IA3-1
	AMOI3=AMOI2+(AW*PAS4)
	APLU3=AMOI3+PAS4
	KRAM3=CONST*SM3/(APLU3*APLU3)
	KRAP3=CONST*SP3/(AMOI3*AMOI3)
	DO 3200 IB3=1,2
	BW=IB3-1
	BMOI3=BMOI2+(BW*PAS4)
	BPLU3=BMOI3+PAS4
	KRBM3=CONST/(BPLU3*BPLU3)
	KRBP3=CONST/(BMOI3*BMOI3)
	DO 3300 IC3=1,2
	CW=IC3-1
	CMOI3=CMOI2+(CW*PAS4)
	CPLU3=CMOI3+PAS4
	ATEST=ABS(AMOI3-AA)
	BTEST=ABS(BMOI3-BB)
	CTEST=ABS(CMOI3-CC)
	BETEST=ABS(BEINF3-BEBE)
	IF(ATEST.LT.0.30.AND.BTEST.LT.0.30.AND.CTEST.LT.0.30.AND.BETEST.LT
     1.0.035)GO TO 3300
	KRCM=CONST*SM3/(CPLU3*CPLU3)
	KRCP=CONST*SP3/(CMOI3*CMOI3)
	KRACM=2.*CONST*CSM3/(APLU3*CPLU3)
	KRACP=2.*CONST*CSP3/(AMOI3*CMOI3)
	CALL MODHKL(KRAM3,KRAP3,KRCM,KRCP,KRACM,KRACP,KRBM3,KRBP3)
	DO 124 I=6,N
	IF(KQT(I).LT.0)GO TO 3300
124	CONTINUE
	DO 4000 IBET4=1,2
	BETW=IBET4-1
	BEINF4=BEINF3+(BETW*PAB8)
	BESUP4=BEINF4+PAB8
	SM4=1./(SIN(BEINF4))**2
	SP4=1./(SIN(BESUP4))**2
	CSM4=COS(BEINF4)*SM4
	CSP4=COS(BESUP4)*SP4
	DO 4100 IA4=1,2
	AW=IA4-1
	AMOI4=AMOI3+(AW*PAS8)
	APLU4=AMOI4+PAS8
	KRAM4=CONST*SM4/(APLU4*APLU4)
	KRAP4=CONST*SP4/(AMOI4*AMOI4)
	DO 4200 IB4=1,2
	BW=IB4-1
	BMOI4=BMOI3+(BW*PAS8)
	BPLU4=BMOI4+PAS8
	KRBM4=CONST/(BPLU4*BPLU4)
	KRBP4=CONST/(BMOI4*BMOI4)
	DO 4300 IC4=1,2
	CW=IC4-1
	CMOI4=CMOI3+(CW*PAS8)
	CPLU4=CMOI4+PAS8
	ATEST=ABS(AMOI4-AA)
	BTEST=ABS(BMOI4-BB)
	CTEST=ABS(CMOI4-CC)
	BETEST=ABS(BEINF4-BEBE)
	IF(ATEST.LT.0.30.AND.BTEST.LT.0.30.AND.CTEST.LT.0.30.AND.BETEST.LT
     1.0.035)GO TO 4300
	KRCM=CONST*SM4/(CPLU4*CPLU4)
	KRCP=CONST*SP4/(CMOI4*CMOI4)
	KRACM=2.*CONST*CSM4/(APLU4*CPLU4)
	KRACP=2.*CONST*CSP4/(AMOI4*CMOI4)
	CALL MODHKL(KRAM4,KRAP4,KRCM,KRCP,KRACM,KRACP,KRBM4,KRBP4)
	DO 125 I=6,N
	IF(KQT(I).LT.0)GO TO 4300
125	CONTINUE
	DO 5000 IBET5=1,2
	BETW=IBET5-1
	BEINF5=BEINF4+(BETW*PAB16)
	BESUP5=BEINF5+PAB16
	SM5=1./(SIN(BEINF5))**2
	SP5=1./(SIN(BESUP5))**2
	CSM5=COS(BEINF5)*SM5
	CSP5=COS(BESUP5)*SP5
	DO 5100 IA5=1,2
	AW=IA5-1
	AMOI5=AMOI4+(AW*PAS16)
	APLU5=AMOI5+PAS16
	KRAM5=CONST*SM5/(APLU5*APLU5)
	KRAP5=CONST*SP5/(AMOI5*AMOI5)
	DO 5200 IB5=1,2
	BW=IB5-1
	BMOI5=BMOI4+(BW*PAS16)
	BPLU5=BMOI5+PAS16
	KRBM5=CONST/(BPLU5*BPLU5)
	KRBP5=CONST/(BMOI5*BMOI5)
	DO 5300 IC5=1,2
	CW=IC5-1
	CMOI5=CMOI4+(CW*PAS16)
	CPLU5=CMOI5+PAS16
	ATEST=ABS(AMOI5-AA)
	BTEST=ABS(BMOI5-BB)
	CTEST=ABS(CMOI5-CC)
	BETEST=ABS(BEINF5-BEBE)
	IF(ATEST.LT.0.30.AND.BTEST.LT.0.30.AND.CTEST.LT.0.30.AND.BETEST.LT
     1.0.035)GO TO 5300
	KRCM=CONST*SM5/(CPLU5*CPLU5)
	KRCP=CONST*SP5/(CMOI5*CMOI5)
	KRACM=2.*CONST*CSM5/(APLU5*CPLU5)
	KRACP=2.*CONST*CSP5/(AMOI5*CMOI5)
	CALL MODHKL(KRAM5,KRAP5,KRCM,KRCP,KRACM,KRACP,KRBM5,KRBP5)
c	DO 126 I=6,N
	I = 1
126  	IF(KQT(I).LT.0)GO TO 5300
C  126 CONTINUE
	I = I + 1
	IF(I .LE. N)GO TO 126
C      DO 6000 IBET6 =1,2
	IBET6 = 1
6001	BETW=IBET6-1
	BEINF6=BEINF5+(BETW*PAB32)
	BESUP6=BEINF6+PAB32
	SM6=1./(SIN(BEINF6))**2
	SP6=1./(SIN(BESUP6))**2
	CSM6=COS(BEINF6)*SM6
	CSP6=COS(BESUP6)*SP6
C      DO 6100 IA6=1,2
	IA6 = 1
6101	AW=IA6-1
	AMOI6=AMOI5+(AW*PAS32)
	APLU6=AMOI6+PAS32
	KRAM6=CONST*SM6/(APLU6*APLU6)
	KRAP6=CONST*SP6/(AMOI6*AMOI6)
C      DO 6200 IB6=1,2
	IB6 = 1
6201	BW=IB6-1
	BMOI6=BMOI5+(BW*PAS32)
	BPLU6=BMOI6+PAS32
	KRBM6=CONST/(BPLU6*BPLU6)
	KRBP6=CONST/(BMOI6*BMOI6)
C      DO 6300 IC6=1,2
	IC6 = 1
6301	CW=IC6-1
	CMOI6=CMOI5+(CW*PAS32)
	CPLU6=CMOI6+PAS32
	ATEST=ABS(AMOI6-AA)
	BTEST=ABS(BMOI6-BB)
	CTEST=ABS(CMOI6-CC)
	BETEST=ABS(BEINF6-BEBE)
	IF(ATEST.LT.0.30.AND.BTEST.LT.0.30.AND.CTEST.LT.0.30.AND.BETEST.LT
     1.0.035)GO TO 6300
	KRCM=CONST*SM6/(CPLU6*CPLU6)
	KRCP=CONST*SP6/(CMOI6*CMOI6)
	KRACM=2.*CONST*CSM6/(APLU6*CPLU6)
	KRACP=2.*CONST*CSP6/(AMOI6*CMOI6)
	CALL MODHKL(KRAM6,KRAP6,KRCM,KRCP,KRACM,KRACP,KRBM6,KRBP6)
C      DO 127 I=6,N
	I = 6
127	IF(KQT(I).LT.0)GO TO 6300
C  127 CONTINUE
	I = I + 1
	IF(I .LE. N)GO TO 127
C      DO 7000 IBET7=1,2
	IBET7 = 1
7001	BETW=IBET7-1
	BEINF7=BEINF6+(BETW*PAB64)
	BESUP7=BEINF7+PAB64
	SM7=1./(SIN(BEINF7))**2
	SP7=1./(SIN(BESUP7))**2
	CSM7=COS(BEINF7)*SM7
	CSP7=COS(BESUP7)*SP7
C      DO 7100 IA7=1,2
	IA7 = 1
7101	AW=IA7-1
	AMOI7=AMOI6+(AW*PAS64)
	APLU7=AMOI7+PAS64
	RAM7=SM7/APLU7**2
	RAP7=SP7/AMOI7**2
	KRAM7=CONST*SM7/(APLU7*APLU7)
	KRAP7=CONST*SP7/(AMOI7*AMOI7)
C      DO 7200 IB7=1,2
	IB7 = 1
7201	BW=IB7-1
	BMOI7=BMOI6+(BW*PAS64)
	BPLU7=BMOI7+PAS64
	RBM7=1./BPLU7**2
	RBP7=1./BMOI7**2
	KRBM7=CONST/(BPLU7*BPLU7)
	KRBP7=CONST/(BMOI7*BMOI7)
C      DO 7300 IC7=1,2
	IC7 = 1
7301	CW=IC7-1
	CMOI7=CMOI6+(CW*PAS64)
	CPLU7=CMOI7+PAS64
	ATEST=ABS(AMOI7-AA)
	BTEST=ABS(BMOI7-BB)
	CTEST=ABS(CMOI7-CC)
	BETEST=ABS(BEINF7-BEBE)
	IF(ATEST.LT.0.30.AND.BTEST.LT.0.30.AND.CTEST.LT.0.30.AND.BETEST.LT
     1.0.035)GO TO 7300
	RCM=SM7/CPLU7**2
	RCP=SP7/CMOI7**2
	RACM=2.*CSM7/(APLU7*CPLU7)
	RACP=2.*CSP7/(AMOI7*CMOI7)
	KRCM=CONST*SM7/(CPLU7*CPLU7)
	KRCP=CONST*SP7/(CMOI7*CMOI7)
	KRACM=2.*CONST*CSM7/(APLU7*CPLU7)
	KRACP=2.*CONST*CSP7/(AMOI7*CMOI7)
	CALL MODHKL(KRAM7,KRAP7,KRCM,KRCP,KRACM,KRACP,KRBM7,KRBP7)
C      DO 128 I=6,N
	I = 6
128	IF(KQT(I).LT.0)GO TO 7300
C  128 CONTINUE
	I = I + 1
	IF(I .LE. N)GO TO 128
	AA=AMOI7+(PAS64/2.)
	BB=BMOI7+(PAS64/2.)
	CC=CMOI7+(PAS64/2.)
	BEBE=BEINF7+(PAB64/2.)
	SINBE=SIN(BEBE)
	RSINB=1./SINBE
	AST=RSINB/AA
	CST=RSINB/CC
	QAR=AST*AST
	QBR=1./(BB*BB)
	QCR=CST*CST
	QER=-2.*AST*CST*SQRT(1.-SINBE*SINBE)
C      DO 15 I=1,N
C      DO 16 J=1,12
C   16 QCAL(I,J)=-1.
C   15 CONTINUE
	I = 1
16	J = 1
15	QCAL(I,J) = -1
	J = J + 1
	IF(J .LE. 12)GO TO 15
	I = I + 1
	IF(I .LE. N)GO TO 16
	NML2=2*ML2+1
C      DO 80 M=1,MH2
	M = 1
8001	MXH=M-1
	CARH=MXH*MXH
	DM1=CARH*RAM7
	DP1=CARH*RAP7
C      DO 90 K=1,MK2
	K = 1
9001	    MXK=K-1
	CARK=MXK*MXK
	DM2=(CARK*RBM7)+DM1
	DP2=(CARK*RBP7)+DP1
C      DO 100 L=1,NML2
	L = 1
10001	 MXL=-ML2-1+L
	CARL=MXL*MXL
	PROHL=MXH*MXL
	DM3=(CARL*RCM)+DM2
	DP3=(CARL*RCP)+DP2
	IF(MXL)101,102,102
101	QQM=DM3-(PROHL*RACP)
	IF((QQM-EPSQ(N)).GT.Q(N))GO TO 100
	QQP=DP3-(PROHL*RACM)
	GO TO 103
102	QQM=DM3-(PROHL*RACM)
	IF((QQM-EPSQ(N)).GT.Q(N))GO TO 90
	QQP=DP3-(PROHL*RACP)
103	IF((QQP+EPSQ(1)).LT.Q(1))GO TO 100
C      DO 119 I=1,N
	I = 1
1191	IF((QQM-EPSQ(I)).GT.Q(I))GO TO 119
	IF((QQP+EPSQ(I)).LT.Q(I))GO TO 119
	J=1
111	IF(QCAL(I,J))112,113,113
113	J=J+1
	IF(J-10)111,111,119
112	QCAL(I,J)=CARH*QAR+CARK*QBR+CARL*QCR+PROHL*QER
	XH(I,J)=MXH
	XK(I,J)=MXK
	XL(I,J)=MXL
C  119 CONTINUE
119	I = I + 1
	IF(I .LE. N)GO TO 1191
C  100 CONTINUE
100	L = L + 1
	IF(L .LE. NML2)GO TO 10001
C   90 CONTINUE
90	K = K + 1
	IF(K .LE. MK2)GO TO 9001
C   80 CONTINUE
80	M = M + 1
	IF(M .LE. MH2)GO TO 8001
	ANG7=BEBE/PIDEG
	V=AA*BB*CC*SIN(BEBE)
C---STORE SOLUTION FOR SUBSEQUENT COPYING TO REFINE FILE (ULCC)
C      DO 190 I=1,6
	I = 1
190	QCELL(I)=10000.*QRCELL(I)
	I = I + 1
	IF(I .LE. 6)GO TO 190
	CALL REFINC(QCELL)
	WRITE(6,1)AA,BB,CC,ANG7,V
C      DO 210 I=1,N
	I = 1
2101	J=1
201	IF(QCAL(I,J))210,202,202
202	DCAL=1./SQRT(QCAL(I,J))
	DIF=DCAL-D(I)
	WRITE(6,2)D(I),XH(I,J),XK(I,J),XL(I,J),DCAL,DIF,EPSIL(I)
	J=J+1
	GO TO 201
C  210 CONTINUE
210	I = I + 1
	IF(I .LE. N)GO TO 2101
	 WRITE(6,3)
C 7300 CONTINUE
7300	IC7 = IC7 + 1
	IF(IC7 .LE. 2)GO TO 7301
C 7200 CONTINUE
7200	IB7 = IB7 + 1
	IF(IB7 .LE. 2)GO TO 7200
C 7100 CONTINUE
7100	IA7 = IA7 + 1
	IF(IA7 .LE. 2)GO TO 7101
C 7000 CONTINUE
7000	IBET7 = IBET7 + 1
	IF(IBET7 .LE. 2)GO TO 7001
C 6300 CONTINUE
6300	IC6 = IC6 + 1
	IF(IC6 .LE. 2)GO TO 6301
C 6200 CONTINUE
6200	IB6 = IB6 + 1
	IF(IB6 .LE. 2)GO TO 6201
C 6100 CONTINUE
6100	IA6 = IA6 + 1
	IF(IA6 .LE. 2)GO TO 6101
C 6000 CONTINUE
6000	IBET6 = IBET6 + 1
	IF(IBET6 .LE. 2)GO TO 6001
5300	CONTINUE
5200	CONTINUE
5100	CONTINUE
5000	CONTINUE
4300	CONTINUE
4200	CONTINUE
4100	CONTINUE
4000	CONTINUE
3300	CONTINUE
3200	CONTINUE
3100	CONTINUE
3000	CONTINUE
2300	CONTINUE
2200	CONTINUE
2100	CONTINUE
2000	CONTINUE
1300	CONTINUE
1200	CONTINUE
1100	CONTINUE
1000	CONTINUE
	RETURN
	END
	
	SUBROUTINE MODHKL(KRAM,KRAP,KRCM,KRCP,KRACM,KRACP,KRBM,KRBP)
	COMMON D(50),EPSIL(50),Q(50),EPSQ(50),KQ(50),KEPSQ(50),KQT(50),N,V
     1,MH,MH2,MK,MK2,ML,ML2,PAS,PAS2,PAS4,PAS8,PAS16,PAS32,PAS64,CONST,V
     2INF,VSUP
	DO 15 I=1,N
	KQT(I)=-1
15	CONTINUE
	NML=2*ML+1
	DO 10 M=1,MH
	MXH=M-1
	MCARH=MXH*MXH
	MM1=MCARH*KRAM
	MP1=MCARH*KRAP
	DO 20 K=1,MK
	MXK=K-1
	MCARK=MXK*MXK
	MM2=(MCARK*KRBM)+MM1
	MP2=(MCARK*KRBP)+MP1
	DO 30 L=1,NML
	MXL=-ML+L-1
	MCARL=MXL*MXL
	MPROHL=MXH*MXL
	MM3=(MCARL*KRCM)+MM2
	MP3=(MCARL*KRCP)+MP2
	IF(MXL)50,51,51
50	KQQM=MM3-(MPROHL*KRACP)
	IF((KQQM-KEPSQ(5)).GT.KQ(5))GO TO 30
	KQQP=MP3-(MPROHL*KRACM)
	GO TO 52
51	KQQM=MM3-(MPROHL*KRACM)
	IF((KQQM-KEPSQ(5)).GT.KQ(5))GO TO 20
	KQQP=MP3-(MPROHL*KRACP)
52	IF((KQQP+KEPSQ(1)).LT.KQ(1))GO TO 30
	DO 40 I=1,5
	KBORIN=KQQM-KEPSQ(I)
	IF(KBORIN.GT.KQ(I))GO TO 40
	KBORSU=KQQP+KEPSQ(I)
	IF(KBORSU.GE.KQ(I))KQT(I)=1
40	CONTINUE
30	CONTINUE
20	CONTINUE
10	CONTINUE
	DO 70 I=1,5
	IF(KQT(I).LT.0)GO TO 3000
70	CONTINUE
	NML2=2*ML2+1
	DO 80 M=1,MH2
	MXH=M-1
	MCARH=MXH*MXH
	MM1=MCARH*KRAM
	MP1=MCARH*KRAP
	DO 90 K=1,MK2
	MXK=K-1
	MCARK=MXK*MXK
	MM2=(MCARK*KRBM)+MM1
	MP2=(MCARK*KRBP)+MP1
	DO 100 L=1,NML2
	MXL=-ML2-1+L
	MCARL=MXL*MXL
	MPROHL=MXH*MXL
	MM3=(MCARL*KRCM)+MM2
	MP3=(MCARL*KRCP)+MP2
	IF(MXL)101,102,102
101	KQQM=MM3-(MPROHL*KRACP)
	IF((KQQM-KEPSQ(N)).GT.KQ(N))GO TO 100
	KQQP=MP3-(MPROHL*KRACM)
	GO TO 103
102	KQQM=MM3-(MPROHL*KRACM)
	IF((KQQM-KEPSQ(N)).GT.KQ(N))GO TO 90
	KQQP=MP3-(MPROHL*KRACP)
103	IF((KQQP+KEPSQ(6)).LT.KQ(6))GO TO 100
	DO 110 I=6,N
	KBORIN=KQQM-KEPSQ(I)
	IF(KBORIN.GT.KQ(I))GO TO 110
	KBORSU=KQQP+KEPSQ(I)
	IF(KBORSU.GE.KQ(I))KQT(I)=1
110	CONTINUE
100	CONTINUE
90	CONTINUE
80	CONTINUE
3000	RETURN
	END
	
	SUBROUTINE REFSET
C---INITIALISES TEMP STORAGE FOR REFINEMENT
	COMMON D(50),EPSIL(50),Q(50),EPSQ(50),KQ(50),KEPSQ(50),KQT(50),N,V
     1,MH,MH2,MK,MK2,ML,ML2,PAS,PAS2,PAS4,PAS8,PAS16,PAS32,PAS64,CONST,V
     2INF,VSUP
	COMMON /REFCOM/ TITLE(20),ALAMBD,AMOLWT,DOBS,NSOL,NOBS,ITEMCH,
     1 IREFCH
C---USED TO PASS TITLE TO S/R REFSET FOR REFORMATTING (ULCC)
	COMMON /IDTCOM/ IDENT(30)
	ITEMCH=38
	IREFCH=39
        OPEN(UNIT=ITEMCH, STATUS='SCRATCH', DISP='DELETE')
        OPEN(UNIT=IREFCH, STATUS='SCRATCH', DISP='DELETE')
	NSOL=0
	NOBS=N
	REWIND ITEMCH
	WRITE(ITEMCH,1) (IDENT(I),I=1,30)
1	FORMAT(30A2,20X)
	REWIND ITEMCH
	READ(ITEMCH,2) (TITLE(I),I=1,20)
2	FORMAT(20A4)
	ALAMBD=0.
	AMOLWT=0.
	DOBS=0.
	REWIND ITEMCH
	RETURN
	END
	
	SUBROUTINE REFINC(QCELL)
C---INCREMENTS NSOL AND WRITES QCELL TO FILE ITEMCH FOR TEMP STORAGE
C---PRIOR TO RE-READING, OUTPUT AND REFINEMENT (QCELL IS IN QU.)
	DIMENSION QCELL(6)
	COMMON /REFCOM/ TITLE(20),ALAMBD,AMOLWT,DOBS,NSOL,NOBS,ITEMCH,
     1 IREFCH
	NSOL=NSOL+1
	WRITE(ITEMCH,1) (QCELL(I),I=1,6)
1	FORMAT(6F10.5)
	RETURN
	END
	
	SUBROUTINE REFSAV
C---CALLED WHEN SOLUTION-FINDING IS COMPLETE. ALL NSOL SOLUTIONS FOUND
C---ARE READ BACK FROM TEMPORARY STORAGE CHANNEL ITEMCH AND WRITTEN TO
C---CHANNEL IREFCH IN THE FORM OF CARD IMAGES THAT MAKE UP A DATA
C---DECK FOR FZON(6), WITH THE NSOL SOLUTIONS AS SPECIFIED LATTICES
C---TO BE NORMALISED, REFINED AND EVALUATED
	COMMON D(50),EPSIL(50),Q(50),EPSQ(50),KQ(50),KEPSQ(50),KQT(50),N,V
     1,MH,MH2,MK,MK2,ML,ML2,PAS,PAS2,PAS4,PAS8,PAS16,PAS32,PAS64,CONST,V
     2INF,VSUP
	COMMON /REFCOM/ TITLE(20),ALAMBD,AMOLWT,DOBS,NSOL,NOBS,ITEMCH,
     1 IREFCH
	DIMENSION QCELL(6)
	ENDFILE ITEMCH
	REWIND ITEMCH
	REWIND IREFCH
	NR=NSOL
	IF(NR.GT.99) NR=99
	WRITE(IREFCH,100) (TITLE(I),I=1,20),ALAMBD,NR,AMOLWT,DOBS
100	FORMAT(20A4/
     1 1H9,19X,F10.6,3X,1H1,6X,I4,3X,1H1,2X,F10.3,F10.4,10X)
	WRITE(IREFCH,150) (D(I),I=1,NOBS)
150	FORMAT(F10.6,70X)
	WRITE(IREFCH,160)
160	FORMAT(1H0,79X)
	DO 200 J=1,NSOL
	READ(ITEMCH,180) (QCELL(I),I=1,6)
	WRITE(IREFCH,180) (QCELL(I),I=1,6)
180	FORMAT(6F10.5,20X)
200	CONTINUE
	WRITE(IREFCH,250)
250	FORMAT(1H0,79X/3HEND,77X)
	ENDFILE IREFCH
	REWIND IREFCH
	RETURN
	END
	
	SUBROUTINE FZONRF
C  ---------------------------------------------------------------------
C---THIS IS A CUT-DOWN VERSION OF JAN VISSER'S ZONE-INDEXING PROGRAM,
C--- DESIGNED TO BE CALLED AS A SUBROUTINE FOR NORMALISING, REFINING,
C--- EVALUATING AND DISPLAYING SPECIFIED SOLUTIONS. IT EXPECTS DATA
C--- IN THE FORM OF A DECK OF CARD IMAGES IN THE =REFINE= FILE
C--- (CHANNEL 9). AS THIS CAN ALSO BE PRINTED TO GIVE A LISTING OF THE
C--- INPUT DATA, SUBROUTINE ECHO HAS BEEN MADE A DUMMY IN THIS VERSION.
C-- *** BASED ON VERSION 6
C-- *** INCLUDING UNIVERSITY OF SURREY MODS UP TO FZON6-7A AND 8
C-- *** DATED OCTOBER 1977
C-- *** ROBIN SHIRLEY, CHEMICAL PHYSICS DEPARTMENT
C  ---------------------------------------------------------------------
C
C             THIS PROGRAM WAS WRITTEN BY
C                        J.W. VISSER
C             OF THE TECHNISCH PHYSISCHE DIENST TNO-TH
C             P.O.BOX 155 DELFT NETHERLANDS.
C              IT CONTAINS MANY IDEAS AND SUGGESTIONS FROM
C                       PROF.DR. P.M. DE WOLFF
C             OF THE TECHNICAL UNIVERSITY (T.H) DELFT.
C
C
C             A DESCRIPTION OF THE PROGRAM WAS PUBLISHED IN
C             JOURNAL OF APPLIED CRYSTALLOGRAPHY.(1969) 2,89.
C
C
C             THIS WORK WAS LARGELY SPONSORED BY THE
C             JOINT COMMITTEE ON POWDER  DIFFRACTION STANDARDS
C             NOW AT  1601 PARK LANE, SWARTHMORE, PA. 19081
C             AND THE CENTRAL ORGANIZATION FOR  TNO, THE HAGUE,
C             NETHERLANDS.
C
C
C-----------------------------------------------------------------------
C     THE INPUT OF THIS PROGRAM SHOULD CONSIST OF0
C     1. A TITLE CARD,CONTAINING UP TO 80 CHARACTERS
C     2. A PARAMETER CARD,WHICH MAY BE  COMPLETELY BLANK.
C     3. A NUMBER OF LINES CARDS, CONTAINING THE POSITIONS OF THE LINES
C        THESE CAN BE ENTERED AS D-,Q- OR TWOTHETA VALUES. ( 8 F10.5 )
C        YOU MAY LEAVE ANY NUMBER OF POSITIONS PER CARD  BLANK.
C        THIS MEANS THAT YOU CAN ALSO ENTER ONLY ONE LINE PER CARD,
C        IN ANY ONE OF THE POSITIONS OF 8F10.5
C     4. A BLANK CARD, TO FINISH THE INPUT OF ONE PROBLEM.
C
C        REPEAT THE ABOVE  AS MANY TIMES AS PROBLEMS,   -- THEN --
C     5. A CARD CONTAINING THE WORD  --  END  -- IN THE FOUR FIRST COLUM
C        EITHER RIGHT- OR LEFT ADJUSTED,  TO FINISH THE INPUT.
C
C  THE PROGRAM WILL TAKE ITS OWN (DEFAULT) PARAMETERS.
C  FOR THE INITIATED THE SUBROUTINE READER GIVES FULL DETAILS ABOUT
C  ALL THE  POSIBILITIES OF THE PROGRAM.
C
C  THE INPUT OF THE PROGRAM IS ON UNIT 5,
C  THE OUTPUT IS ON UNIT  6.
C  ON UNIT 4 YOU MAY GET TEST OUTPUT, THESE ARE THE CARD IMAGES
C  OF THE INPUT CARDS.
C
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C                         GENERAL FLOW OF SUBROUTINES
C
C-----------------------------------------------------------------------
C
C
C        *******            *******           **********
C        *     *            *     * <-------> *  ECHO  *
C        *     *            *  R  *           **********
C        *     *            *     *
C        *     * <--------> *  E  *           **********
C        *     *            *     * <-------> * MANUAL *
C        *     *            *  A  *           **********
C        *     *            *     *
C        *     *            *  D  *           **********
C        *     *            *     * <-------> * INSTRS *
C        *     *            *  E  *           **********
C        *     *            *     *
C        *     *            *  R  *           **********
C        *     *            *     * <-------> * SRTRON *
C        *     *            *******           **********
C        *     *
C        *  M  *            *********
C        *     *  <-------> * NORM3 *
C        *     *            *********
C        *     *
C        *     *            *******           **********
C        *  A  *            *     *  <------> * WISSEL *
C        *     *            *  R  *           **********
C        *     *            *     *           *********
C        *     *            *  E  *  <------> * WRISP *
C        *     *            *     *           *********
C        *     *  <-------> *  F  *           *********
C        *     *            *     *  <------> * ARRAY *
C        *  I  *            *  P  *           *********
C        *     *            *     *           ********
C        *     *            *  A  *  <------> * MINV *
C        *     *            *     *           ********
C        *     *            *  R  *           *********
C        *     *            *     *  <------> * NORM3 *
C        *     *            *     *           *********
C        *     *            *     *           **********
C        *     *            *     *  <------> * LISTM  *
C        *     *            *******           **********
C        *     *
C        *     *
C        *     *            **********
C        *  N  *  <-------> * UNTCEL *
C        *     *            **********
C        *     *
C        *     *
C        *     *            **********
C        *     *  <-------> * SYMTST *
C        *     *            **********
C        *     *
C        *     *
C        *******
C
C-----------------------------------------------------------------------
C
	DIMENSION STORE(200,6 ),BRAVTP(6),ISYST(3,3)
	INTEGER  H,H1
C-----------------------------------------------------------------------
C
C        ****************************
C        * ATTENTION   H = INTEGER  *
C        ****************************
C
C-----------------------------------------------------------------------
	LOGICAL SELF,ENL,TEST,LIST
	LOGICAL RECTAN
CDC---LEVEL STATEMENTS FOR CDC 7600 FOLLOW---REMOVE FOR OTHER MACHINES--
CDC   LEVEL 2, LINCO,NQ2,NSOLMX,ROSTER
C---COMMON BLOCK /TEMRES/, LENGTH 11613 (OCTAL), 5003 (DECIMAL)
	COMMON /LINES /INP,IOUT,JOB,MAX,NAME(20),INTSY(40),OBS(40),CEL(10)
	COMMON /PARMS / TOL2,TOL3,WAVEL,FMINEQ,FMINIM,SELF,ENL,TEST,LIST,
     *                NSYST(3),OBSMT(40),OBSPT(40)
	COMMON /TEMRES/ LINCO,NQ2,NSOLMX,ROSTER(500,10)
	COMMON /REFP/ IREF
	DATA  BRAVTP/1HA,1HB,1HC,1HI,1HF,1HP/
	DATA ISYST(1,1)/4HORTH/,ISYST(1,2)/4HORHO/,ISYST(1,3)/4HMBIC/,
     1     ISYST(2,1)/4HMONO/,ISYST(2,2)/4HCLIN/,ISYST(2,3)/4HIC  /,
     2     ISYST(3,1)/4HTRIC/,ISYST(3,2)/4HLINI/,ISYST(3,3)/4HIC   /
C
	JOB=0
1	RAD=57.29578
	INP=39
C---INPUT WILL BE TAKEN FROM THE =REFINE= FILE
	IOUT=6
	NRMAX=200
	NSOLMX=500
	NREFMX=10
	NORDER=4
	DO 10 I=1,100
	DO 5 J=1,10
5	ROSTER(I,J)=0.0
10	CONTINUE
	DO 30 I=1,200
	DO 20 J=1,6
20	STORE(I,J)=0.0
30	CONTINUE
C
	CALL       READER(NZ1,NZ2,NR,WMOL,DOBS)
	IF (NR .EQ. -1) GOTO 8000
C
	N4=NR
	NR=0
C  ---------------------------------------------------------------------
C--THE LATTICES FOUND WILL BE PUT INTO A STANDARD SHAPE.(NORMALISED)
C  ---------------------------------------------------------------------
	IF (ENL) WRITE (IOUT,9936)
	DO  1990 J=1,N4
	IF (ROSTER(J,7) .LT. 1.0) GO TO  1990
	DO 1950  I=1,6
1950	CEL(I)=ROSTER(J,I)
	CALL NORM3
	DO 1955 I=1,6
1955	ROSTER(J,I)=CEL(I)
	DO 1960 K=1,3
	IF (ROSTER(J,K) .LT. FMINIM) GO TO 1990
1960	CONTINUE
	NR=NR+1
	DO 1965 K=1,8
1965	ROSTER(NR,K)=ROSTER(J,K)
	IF(ENL) WRITE(IOUT,9937) NR,(ROSTER(NR,K),K=1,8)
1990	CONTINUE
C-----------------------------------------------------------------------
C  THE LATTICES ARE NOW NORMALIZED,AND THE PATHOLOGICAL CASES HAVE
C  BEEN DROPPED. THE NEXT STEP IS TO TRY AND REFINE THE LATTICES THAT
C  ARE LEFT.
C  ---------------------------------------------------------------------
C  ---------------------------------------------------------------------
C  FIRST WE ELIMINATE 'EQUAL' LATTICES .
	NR1=NR-1
	DO 2050 H=1,NR1
	H1=H+1
	DO 2040 K=H1,NR
	U=0.0
	DO 2030 L=1,6
2030	U=U+ABS( ROSTER(H,L)- ROSTER(K,L))
	IF (U.GT. FMINEQ) GO TO 2040
	IF ( ROSTER(H,7).GE.  ROSTER(K,7)) GO TO 2035
	 ROSTER(H,7)=0.0
	GO TO 2040
2035	 ROSTER(K,7)=0.0
2040	CONTINUE
2050	CONTINUE
C-----------------------------------------------------------------------
C  A FIRST ROUND OF LEAST SQUARES (IN REFPAR) WILL NOW BE RUN.
C  ---------------------------------------------------------------------
2055	NSOL=0
	NREFCY=3
	DO 2100 I=1,NR
	IF ( ROSTER(I,7) .LT. 1.0) GO TO 2100
	DO 2060 K=1,10
2060	CEL(K)= ROSTER(I,K)
	CEL(7)=0.0
	CEL(10)=CEL(8)
	NSOL=NSOL+1
	DO 2065 IREF=1,NREFCY
	CALL REFPAR(NSOL)
2065	CONTINUE
	IREF=0
	DO 2070 K=1,10
2070	ROSTER(NSOL,K)=CEL(K)
2100	CONTINUE
	NR=NSOL
C  ---------------------------------------------------------------------
C--A DOUBLE SORT WILL NOW BE CARRIED OUT, FIRST ON THE NUMBER OF
C  LINES INDEXED(POS. 7)   AND SECONDLY ON THE FIGURE OF MERIT (POS. 8)
C--OF THE LATTICES.
C  ---------------------------------------------------------------------
	WRITE(IOUT,9950) NAME
	WRITE(IOUT,9951)
	JJ=100
	DO 3050 L=1,NR
	K=0
	AA=0.
	DO 3040 H=L,NR
	M= ROSTER(H,7)+0.5
	B= ROSTER(H,8)
	IF (M-K) 3040,3020,3010
3010	K=M
	I=H
	AA=B
	GO TO 3040
3020	IF(B.LE.AA) GO TO 3040
	AA=B
	I=H
3040	CONTINUE
	IF (K.GE.JJ) GO TO 3045
	WRITE(IOUT,9952)
	JJ=K
3045	DO 3048 N=1,10
	A= ROSTER(I,N)
	 ROSTER(I,N)= ROSTER(L,N)
	 ROSTER(L,N)=A
	CEL(N)=ROSTER(L,N)
3048	CONTINUE
	CALL UNTCEL(WMOL,DOBS,VOL,0)
	KKK=INT(ROSTER(L,9)+0.5)
	WRITE(IOUT,9953) (ROSTER(L,N),N=1,8),BRAVTP(KKK),
     * ROSTER(L,10), VOL
3050	CONTINUE
C-----------------------------------------------------------------------
C--END OF SORT AND PRINT.
C-----------------------------------------------------------------------
	WRITE(IOUT,9963)
C-----------------------------------------------------------------------
C--WE NOW PRINT OUT THE BEST FOUR SOLUTIONS WITH AN INDEXED LIST
C  OF THE FIRST TWENTY LINES.  IF YOU HAVE SET  'ENL' THEN YOU ALSO
C--GET  THE Q-VALUES PRINTED OUT AS A  LATTICE.
C-----------------------------------------------------------------------
	WRITE(IOUT,9964) NREFMX
	N=0
	SELF=.TRUE.
	DO 3052 I=1,3
3052	IF (NSYST(I) .EQ.1) N=I
	IF (N .NE. 0) SELF=.FALSE.
	IF (SELF) GO TO 3058
	MIN=-1
	N=0
	M=0
	DO 3054 K=1,3
	IF (NSYST(K) .EQ. 1) WRITE(IOUT,9965) (ISYST(K,I),I=1,3)
	IF (NSYST(K) .EQ. 1) N=N+1
	IF (NSYST(K) .EQ.(-1)) WRITE(IOUT,9966) (ISYST(K,I),I=1,3)
	IF (NSYST(K) .EQ.(-1) ) M=M+1
3054	CONTINUE
	IF (N .GT. 0) MIN=1
	IF (M .GT. 0) MIN=0
	DO 3056 M=1,MAX
	DIF=OBS(M)-OBSMT(M)
	OBSMT(M)=OBSMT(M)-DIF
3056	OBSPT(M)=OBSPT(M)+DIF
C  THE TOLERANCES HAVE BEEN DOUBLED WHEN YOU WANT TO  FORCE A CERTAIN
C  CRYSTAL SYSTEM.
3058	NSOL=0
	DO 3100 I=1,NR
	IF(NSOL.EQ.NREFMX) GOTO 3105
	IF (SELF) GO TO 3090
	NORDE=1
	DO 3080 J=4,6
3080	IF (ABS(ROSTER(I,J)) .GT. FMINEQ) NORDE=NORDE+1
	IF (NORDE .EQ. 4) NORDE=3
	IF (NSYST(NORDE) .LT. MIN) GO TO 3100
3090	NSOL=NSOL+1
	DO 3060 J=1,10
3060	CEL(J)= ROSTER(I,J)
	CALL REFPAR(NSOL)
	DO 3070 J=1,10
3070	 ROSTER(I,J)=CEL(J)
3100	CONTINUE
C  AFTER THIS REFINEMENT CYCLE, A NEW DOUBLE SORT IS MADE
3105	NSOL=NREFMX
	IF(NR.LT.NREFMX) NSOL=NR
	DO 3250 L=1,NSOL
	K=0
	AA=0.
	DO 3240 H=L,NR
	M=ROSTER(H,7)+0.5
	B=ROSTER(H,8)
	IF (M-K) 3240,3220,3210
3210	K=M
	I=H
	AA=B
	GO TO 3240
3220	IF(B.LE.AA) GO TO 3240
	AA=B
	I=H
3240	CONTINUE
	DO 3248 N=1,10
	A=ROSTER(I,N)
	ROSTER(I,N)=ROSTER(L,N)
	ROSTER(L,N)=A
3248	CONTINUE
3250	CONTINUE
	WRITE(IOUT,9954) NAME,NREFMX
	DO 3110 I=1,NSOL
	KKK=INT(ROSTER(I,9)+0.5)
3110	WRITE(IOUT,9955) ( ROSTER(I,N),N=1,8),  ROSTER(I,10),BRAVTP(KKK)
	WRITE(IOUT,9956)
	DO 3120 I=1,NSOL
	DO 3115 J=1,10
3115	CEL(J)= ROSTER(I,J)
	CALL UNTCEL(WMOL,DOBS,VOL,1)
3120	CONTINUE
C-----------------------------------------------------------------------
C--PRINTS THE VALUES OF THE PARAMETERS OF THE DIRECT LATTICE.
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
	DO 3150  I=1,NSOL
	DO 3140 J=1,10
3140	CEL(J)=ROSTER(I,J)
	CALL SYMTST
3150	CONTINUE
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C  THE BEST SOLUTIONS ARE TESTED ON HIGHER SYMMETRIES,THAT MAY BE
C  CONCEALED IN THE Q-VALUES.
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
	WRITE(IOUT,9957)
	WRITE(IOUT,9959)
	WRITE(IOUT,9958) NAME
C  ---------------------------------------------------------------------
	GO TO 1
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C--- FINAL RETURN
8000	CONTINUE
	RETURN
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
9936  FORMAT(
     160H1        LATTICES REDUCED AND PUT INTO STANDARD FORM        //
     260H SERIAL   A       B       C        D       E       F       H,
     320HITS   PEDIGREE        )
9937	FORMAT(1H ,I5,6F8.1, 2F8.0)
9950	FORMAT(1H1,20A4)
9951  FORMAT(1H0,14X,47HTRIAL LATTICES AFTER INDEXING AND LS REFINEMENT
     *//14X,52HONLY THE FIRST 20 LINES ARE CONSIDERED FOR INDEXING  //
     *59X,24HLINES  FIGURE    BRAVAIS/
     *95H      A        B        C        D         E        F    INDEXE
     *D  OF MERIT  LATTICE  PEDIGREE   ,6HVOLUME )
9952	FORMAT(1H )
9953	FORMAT(1H ,6F9.1,F9.0,F9.1,8X,A1,F9.0,F8.1 )
9954	FORMAT(1H1,20A4 /// 30X,
     1 4HTHE ,I2,24H MOST PROBABLE SOLUTIONS//// 58X,
     234HLINES  FIGURE OF           BRAVAIS/
     360H    Q(A)     Q(B)     Q(C)      Q(D)     Q(E)     Q(F)   IND,
     433HEXED   MERIT   PEDIGREE  LATTICE        )
9955	FORMAT(1H ,6F9.2,F9.0,F9.1,F9.0,8X,A1)
9956  FORMAT(1H0,20X,38HTHE DIRECT CONSTANTS OF THESE LATTICES//
     *60H     A        B        C       ALFA     BETA     GAMMA      ,
     * 6HVOLUME)
9957	FORMAT(1H0///40X,24H ****   ATTENTION   ****//
     160H 1. ANY 'SOLUTION' WITH A FIGURE OF MERIT LESS THAN 4.0 IS W,
     260HORTHLESS AND SHOULD BE REJECTED RIGHT AWAY.                 //
     360H 2. ANY 'SOLUTION'THAT LEAVES MORE THAN TWO VERY WEAK LINES ,
     460HUNEXPLAINED  IS WORTHLESS.   HOWEVER, IF THE FIGURE OF MERIT/
     560H    IS GREATER THAN 10.0 IT MIGHT BE WORTHWHILE TO EXAMINE  ,
     660HTHE INPUT DATA  AND THE 'SOLUTION' MORE CLOSELY.            //
     760H 3. EVEN 'SOLUTIONS' THAT INDEX ALL LINES WITH A FIGURE OF M,
     860HERIT GREATER THAN 10.0 SHOULD NOT BE ACCEPTED UNCRITICALLY  /
     933H    WITHOUT FURTHER INVESTIGATION//)
9958	FORMAT(1H0,20A4)
9959  FORMAT(71H0*******************************************************
     1***************/2H *,68X,1H*/2H *,  69H  THIS PROGRAM GIVES A NUMB
     2ER OF SUGGESTIONS THAT MAY HELP YOU TO   */ 71H *  FIND THE UNIT C
     3ELL THAT EXPLAINS YOUR EXPERIMENTAL DATA. HOWEVER  */2H *,68X,
     4 1H*/ 2H *,5X,58HWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
     4WWWWWWWWW,         6H     */8H *     W,56X,7HW     */ 8H *     W, 
     5 63H   IT IS ENTIRELY THE USERS RESPONSIBILITY TO DECIDE    W     
     6* /71H *     W   WHETHER ANY OF THE SUGGESTED UNIT CELLS IS THE   
     7    W     *  / 71H *     W   CORRECT CELL.                        
     8                W     * / 8H *     W, 56X, 7HW     * / 7H *      ,
     9 58(1HW), 6H     * / 2H *,68X,1H*/1H ,70(1H*) )
9963	FORMAT(7H0NOTES0//
     160H  1. 'LINES INDEXED' REFERS TO THE NUMBER OF LINES FROM THE ,
     260HFIRST 20 OBSERVED LINES THAT HAVE BEEN ACCEPTED AS INDEXED  /
     360H     BY THIS LATTICE, UP TO A MAXIMUM OF 20, WHICH WOULD REP,
     460HRESENT COMPLETE SUCCESS.                                    //
     560H  2. 'FIGURE OF MERIT' IS DEFINED AS QOBS(20)/(2*MEANDIFF*NC,
     660HALC) WHERE DIFFERENCES ARE TAKEN FOR INDEXED LINES ONLY     /
     760H     AND 'NCALC' IS THE NUMBER OF LINES GENERATED UP TO Q=1.,
     860H003*Q(20),REGARDLESS OF WHETHER OR NOT THEY CONTRIBUTE TO   /
     960H     THE INDEXINGS, BUT EXCLUDING SYSTEMATIC ABSENCES DUE TO,
     *27H THE BRAVAIS LATTICE TYPE. /
     A60H     IT IS CLOSELY RELATED TO  M20 ,DE WOLFF(J.APPL.CRYST.19,
     B60H68,1972), ALLOWING FOR LESS THAN 20 SUCCESSFUL INDEXINGS.   )
9964  FORMAT(5H0THE ,I2,49H BEST SOLUTIONS ARE NOW PRINTED MORE EXTENSIV
     1ELY.)
9965	FORMAT( 1H0,3A4, 23H LATTICES ARE PRINTED.  )
9966	FORMAT( 1H0,3A4,27H LATTICES ARE NOT PRINTED.  )
	END
	
	SUBROUTINE ARRAY(MODE,I,J,N,M,S,D)
C
C  COPIED FROM  THE I.B.M.  SCIENTIFIC SUBROUTINE PACKAGE.
C  TO BE USED IN COMBINATION WITH THE SUBROUTINE  MINV (MATRIX INVERSION
C
	DIMENSION S(36),D(36)
	NI=N-I
	IF(MODE .GT. 1) GO TO 20
	IJ=I*J+1
	NM=N*J+1
	DO 10 K=1,J
	NM=NM-NI
	DO 10 L=1,I
	IJ=IJ-1
	NM=NM-1
	D(NM)=S(IJ)
10	CONTINUE
	GO TO 40
20	IJ=0
	NM=0
	DO 30 K=1,J
	DO 25 L=1,I
	IJ=IJ+1
	NM=NM+1
	S(IJ)=D(NM)
25	CONTINUE
	NM=NM+NI
30	CONTINUE
40	RETURN
	END
	
	SUBROUTINE ECHO
C-----------------------------------------------------------------------
C  THE PURPOSE OF THIS SUBROUTINE IS TO PRINT THE CARD IMAGES OF THE INP
C     CARDS ON A SEPARATE PAGE
C  THIS IS VERY USEFUL FOR TEST PURPOSES AND FOR DOCUMENTATION.
C
C  ON SOME MACHINES IT IS NOT POSSIBLE AT ALL TO  BACKSPACE  THE
C  INPUT FILE.  IN THIS CASE YOU COULD READ THE INPUT DATA FROM
C  THE INPUT FILE(IN ECHO), WRITE THEM AGAIN ON A TEMPORARY FILE
C  AND READ THEM FROM THIS TEMPORARY FILE IN SUBROUTINE  READER.
C  SINCE THIS IS RATHER  'MACHINE-DEPENDENT', I LEAVE IT TO YOUR
C  OWN IMAGINATION AND INGENUITY..............
C-----------------------------------------------------------------------
C  YOU CAN ALSO GET RID OF THE WHOLE  PROBLEM BY REPLACING THE
C  SUBROUTINE BY A  DUMMY SUBROUTINE, AS FOLLOWS
C
C     SUBROUTINE  ECHO
C     RETURN
C     END
C
C-----------------------------------------------------------------------
	RETURN
	END
	
	SUBROUTINE  LISTM(NCALC)
	DIMENSION INSTOR(100,3),QSTOR(100)
	LOGICAL SELF,ENL,TEST,LIST
	COMMON /LINES /INP,IOUT,JOB,MAX,NAME(20),INTSY(40),OBS(40),CEL(10)
	COMMON /PARMS / TOL2,TOL3,WAVEL,FMINEQ,FMINIM,SELF,ENL,TEST,LIST,
     *                NSYST(3),OBSMT(40),OBSPT(40)
	COMMON /LINSTR/ INSTOR,QSTOR
	SWAVEL=WAVEL/200.
	NOLIM=20
	IF(MAX.LE.20) NOLIM=MAX-1
C--SORT CALCULATED VALUES.
	DO 100 I=1,NCALC
	Q=1.E+6
	DO  60 J=I,NCALC
	IF (QSTOR(J).GT. Q) GO TO 60
	Q=QSTOR(J)
	N=J
60	CONTINUE
	U=QSTOR(I)
	QSTOR(I)=QSTOR(N)
	QSTOR(N)=U
	DO  80 K=1,3
	INX  =INSTOR(I,K)
	INSTOR(I,K)=INSTOR(N,K)
80	INSTOR(N,K)= INX
100	CONTINUE
	NO=1
	NC=1
	UL=0.
	NWR=21
	IF(MAX.LT.21) NWR=MAX
	WRITE(IOUT,9903) NWR
	WRITE(IOUT,9904)
C--PRINT QCALC AND QOBS IN ORDER OF TWOTHETA.
110	CONTINUE
	IF (NC .GT. NCALC) GO TO 200
	IF (NO .GT. NWR) GO TO 115
	IF (QSTOR(NC) .GT. OBS(NO)) GO TO 200
115	CONTINUE
	W=SQRT(QSTOR(NC))
	U=SWAVEL*W
	W=100./W
	U=114.59156*ATAN(U/SQRT(1.0-U*U))
	IF (U-UL .GT. 0.05) WRITE(IOUT,9905)
	WRITE(IOUT,9901) U,W,(INSTOR(NC,K),K=1,3),QSTOR(NC)
	UL=U
	NC=NC+1
	IF (NC .LE. NCALC ) GO TO 110
	IF (NO .GE. NWR) GO TO 300
200	W=SQRT(OBS(NO))
	U=SWAVEL*W
	W=100./W
	U=114.59156*ATAN(U/SQRT(1.0-U*U))
	IF (U-UL .GT. 0.05) WRITE(IOUT,9905)
	WRITE(IOUT,9902) U,W, OBS(NO)
	UL=U
	NO=NO+1
	IF(NO.LE.NOLIM) GOTO 110
	IF (NC .LE. NCALC) GO TO 110
300	RETURN
9901	FORMAT(1H ,F10.2,F10.3,2X,3I4,F8.1)
9902	FORMAT(1H ,F10.2,F10.3,14H     *OBSERVED ,F8.1)
9903  FORMAT(53H0LIST OF ALL OBSERVED AND CALCULATED LINES UP TO THE ,  
     1 I2,3H-TH/74H IF THE TWOTHETAS DIFFER MORE THAN 0.05 DEGREES, A BL 
     2ANK LINE IS INSERTED.)
9904	FORMAT(45H0  TWOTHETA       D       H   K   L     Q          )
9905	FORMAT(3X)
	END
	
	SUBROUTINE MINV(A,N,D,L,M)
C
C  MATRIX INVERSION. COPIED FROM THE I.B.M. SCIENTIFIC SUBROUTINE PACKAG
C
	DIMENSION A(36),L(6),M(6)
	D=1.0
	NK=-N
	DO 80 K=1,N
	NK=NK+N
	L(K)=K
	M(K)=K
	KK=NK+K
	BIGA=A(KK)
	DO 20 J=K,N
	IZ=N*(J-1)
	DO 20 I=K,N
	IJ=IZ+I
10	IF(ABS(BIGA).GE.ABS(A(IJ))) GO TO 20
	BIGA=A(IJ)
	L(K)=I
	M(K)=J
20	CONTINUE
	J=L(K)
	IF(J .LE. K) GO TO 35
	KI=K-N
	DO 30 I=1,N
	KI=KI+N
	HOLD=-A(KI)
	JI=KI-K+J
	A(KI)=A(JI)
	A(JI)=HOLD
30	CONTINUE
35	I=M(K)
	IF(I .LE. K) GO TO 45
	JP=N*(I-1)
	DO 40 J=1,N
	JK=NK+J
	JI=JP+J
	HOLD=-A(JK)
	A(JK)=A(JI)
	A(JI)=HOLD
40	CONTINUE
45	IF (ABS(BIGA) .GT. 0.0001) GO TO 48
	D=0.0
	RETURN
48	DO 55 I=1,N
	IF(I .EQ. K) GO TO 55
	IK=NK+I
	A(IK)=A(IK)/(-BIGA)
55	CONTINUE
	DO 65 I=1,N
	IK=NK+I
	HOLD=A(IK)
	IJ=I-N
	DO 65 J=1,N
	IJ=IJ+N
	IF(I .EQ. K) GO TO 65
	IF(J .EQ. K) GO TO 65
	KJ=IJ-I+K
	A(IJ)=HOLD*A(KJ)+A(IJ)
65	CONTINUE
	KJ=K-N
	DO 75 J=1,N
	KJ=KJ+N
	IF(J .EQ. K) GO TO 75
	A(KJ)=A(KJ)/BIGA
75	CONTINUE
	D=D*BIGA
	A(KK)=1.0/BIGA
80	CONTINUE
	K=N
100	K=(K-1)
	IF(K .LE. 0) GO TO 150
	I=L(K)
	IF(I .LE. K) GO TO 120
	JQ=N*(K-1)
	JR=N*(I-1)
	DO 110 J=1,N
	JK=JQ+J
	HOLD=A(JK)
	JI=JR+J
	A(JK)=-A(JI)
	A(JI)=HOLD
110	CONTINUE
120	J=M(K)
	IF(J .LE. K) GO TO 100
	KI=K-N
	DO 130 I=1,N
	KI=KI+N
	HOLD=A(KI)
	JI=KI-K+J
	A(KI)=-A(JI)
	A(JI)=HOLD
130	CONTINUE
	GO TO 100
150	RETURN
	END
	
	SUBROUTINE NORM3
C-----------------------------------------------------------------------
C  THIS SUBROUTINE 'NORMALIZES THE LATTICES AND PUTS THEM INTO A
C  STANDARD SHAPE
C-----------------------------------------------------------------------
C
	DIMENSION  OLD(6),          IZON(3,6),QVAL(5,5,5)
	INTEGER H,H1
	LOGICAL SELF,ENL,TEST,LIST,CHANGE
	COMMON /LINES /INP,IOUT,JOB,MAX,NAME(20),INTSY(40),OBS(40),CEL(10)
	COMMON /PARMS / TOL2,TOL3,WAVEL,FMINEQ,FMINIM,SELF,ENL,TEST,LIST,
     *                NSYST(3),OBSMT(40),OBSPT(40)
	EQUIVALENCE (CEL(1),A),(CEL(2),B),(CEL(3),C),(CEL(4),D),
     *(CEL(5),E),(CEL(6),F),(FMINEQ,EQMIN)
	DO 80 I=1,6
80	OLD(I)=CEL(I)
C  INITIALIZING *OLD*
	GO TO 100
C
C-----------------------------------------------------------------------
C  THE NEXT PART IS THE OLD PROCEDURE 'MAZEN'
C  MAZEN FINDS THE  VALUES OF D,E AND F IN A LATTICE, WHEN THE COORDINAT
C  OF A,B,AND C HAVE BEEN GIVEN.
10	I1=IZON(H,1)+IZON(K,1)-3
	I2=IZON(H,2)+IZON(K,2)-3
	I3=IZON(H,3)+IZON(K,3)-3
	X=QVAL(I1,I2,I3)
	I1=IZON(H,1)-IZON(K,1)+3
	I2=IZON(H,2)-IZON(K,2)+3
	I3=IZON(H,3)-IZON(K,3)+3
	Y=QVAL(I1,I2,I3)
	Z=(X-Y)/2.0
	GO TO (601,602,603),ISW
C-----------------------------------------------------------------------
C
C
C-----------------------------------------------------------------------
C  THE FOLLOWING IS THE PROCEDURE 'CHANGE'
C  IF ANY CONSTANT HAS BEEN  CHANGED,THE PROGRAM REPEATS UNTIL NO
C  MORE CHANGES ARE DETECTED.
20	U=0.0
	DO 30 I=1,6
	U=U+ ABS(CEL(I)   -OLD(I))
	OLD(I)=CEL(I)
30	CONTINUE
	CHANGE=.FALSE.
	IF (U.GT.0.1 ) CHANGE=.TRUE.
	GO TO (310,610),ISW
C-----------------------------------------------------------------------
C
100	IF (TEST) WRITE(IOUT,9930) A,B,C,D,E,F
	QMAX=10000.0
C--CHECK ON  CERTAIN RELATIONS BETWEEN A,B,C,D,E,F  WHICH MEAN THAT
C--THE LATTICE HAS HIGHER SYMMETRY.
200	DO 300 II=1,3
	X=A
	A=B
	B=C
	C=X
	X=D
	D=E
	E=F
	F=X
	IF (B .LT. FMINIM) GO TO 1000
	IF (D.GT. 0.0)  GO TO 205
	D=-D
	E=-E
205	IF (E.GT. 0.0)  GO TO 210
	E=-E
	F=-F
210	IF (ABS(A-C).LT.EQMIN  .AND.   E.GT. EQMIN) GO TO 220
	GO TO 230
220	A=(A+C+E)/4.
	C=A-E/2.
	F=(D+F)/2.
	D=ABS(D-F)
	E=0
230	IF (ABS(E-A).GT.EQMIN ) GO TO 240
	A=A/4.
	C=C-A
	E=0.0
	F=F/2.
	D=ABS(D-F)
240	IF (ABS(E-C).GT.EQMIN) GO TO 250
	C=C/4.0
	A=A-C
	E=0.0
	D=D/2.
	F=ABS(D-F)
250	IF (ABS(D) .LE. EQMIN) D=0.0
	IF (ABS(E) .LE. EQMIN) E=0.0
	IF (ABS(F) .LE. EQMIN) F=0.0
300	CONTINUE
	IF (TEST) WRITE(IOUT,9931) A,B,C,D,E,F
	ISW=1
	GO TO 20
C--TEST ON CHANGES.
310	IF (CHANGE) GO TO 200
C--CALCULATE  THE Q-VALUES NEAR THE ORIGIN
	DO 400 L1=3,4
	L=L1-3
	DO 390 K1=2,4
	K=K1-3
	DO 380 H1=2, 4
	H=H1-3
	IF (L  .EQ.0 .AND. K .LT.1) GO TO 330
	X=A*FLOAT(H*H)+B*FLOAT(K*K)+C*FLOAT(L*L)+D*FLOAT(K*L)+E*FLOAT(H*L)
     1  +F*FLOAT(H*K)
	IF (X .GT.(-0.1))GO TO 320
	A=0.0
	GO  TO  1000
320	QVAL (H1,K1,L1)=X
	GO TO 380
330	QVAL(H1,K1,L1)=QMAX
380	CONTINUE
390	CONTINUE
400	CONTINUE
	QVAL(4,3,3) =A
C--NOW FIND THE SMALLEST Q-VALUES.
	N1=0
410	N1=N1+1
	U=QMAX-10.0
	DO 500 L=3,4
	DO 490 K=2,4
	DO 480 H=2,4
	IF (QVAL(H,K,L).GE.U) GO TO 480
	U=QVAL(H,K,L)
	IZON(N1,3)=L
	IZON(N1,1)=H
	IZON(N1,4)=H
	IZON(N1,2)=K
	IZON(N1,5)=K
480	CONTINUE
490	CONTINUE
500	CONTINUE
	I1=IZON(N1,1)
	I2=IZON(N1,2)
	I3=IZON(N1,3)
	QVAL(I1,I2,I3)=QMAX
	IF (N1.LT.3) GO TO 410
C--CHECK WHETHER THEY ARE IN ONE PLANE.
	II=0
	DO  510 I=1,3
	K=IZON(1,I+1)-3
	L=IZON(2,I+2)-3
	M=IZON(2,I+1)-3
	N=IZON(1,I+2)-3
	J=K*L-M*N
	II=II+J*(IZON(3,I)-3)
510	CONTINUE
	IF (II .NE. 0) GO TO 520
	N1=2
	GO TO  410
C--CALCULATE THE RECIPROCAL LATTICE AGAIN., TO DETERMINE D,E, AND F.
520	DO 580 L1=1,5
	L=L1-3
	DO 570 K1=1,5
	K=K1-3
	DO 560 H1=1,5
	H=H1-3
	Q=A*FLOAT(H*H)+B*FLOAT(K*K)+C*FLOAT(L*L)+D*FLOAT(K*L)+E*FLOAT(H*L)
     1  +F*FLOAT(H*K)
	IF (Q .GT.(-0.1)) GO TO 550
	A=0.0
	GO TO 1000
C  IF ANY Q-VALUE IS NEGATIVE,THE LATTICE SHOULD BE DISCARDED
550	QVAL(H1,K1,L1)=Q
560	CONTINUE
570	CONTINUE
580	CONTINUE
	I11=IZON(1,1)
	I12=IZON(1,2)
	I13=IZON(1,3)
	I21=IZON(2,1)
	I22=IZON(2,2)
	I23=IZON(2,3)
	I31=IZON(3,1)
	I32=IZON(3,2)
	I33=IZON(3,3)
	B=QVAL(I11,I12,I13)
	A=QVAL(I21,I22,I23)
	C=QVAL(I31,I32,I33)
C  I KNOW, IKNOW,    BUT THIS PIDGIN-FORTRAN WORKS ON ALL MACHINES....
C--NOW FIND D,E AND F.
	ISW=1
	H=1
	K=3
	GO TO  10
601	D=Z
	ISW=2
	K=2
	GO TO 10
602	F=Z
	ISW=3
	H=3
	GO TO 10
603	E=Z
	IF (TEST) WRITE(IOUT,9934) (CEL(I),I=1,6)
	ISW=2
	GO TO 20
C--TEST ON CHANGES.
610	IF (CHANGE) GO TO 200
C--FROM HERE TO LABEL 670, H COUNTS THE NUMBER OF TERMS AMONG D,E,F
C--THAT ARE EQUAL TO ZERO, WHILE K POINTS OUT WHICH ONE IF THERE
C--IS ONLY ONE.
	H=0
	IF (ABS(D).GT.EQMIN) GO TO 620
	D=0.0
	GO TO  621
620	H=H+1
	K=1
621	IF (ABS(E).GT.EQMIN) GO TO 630
	E=0.0
	GO TO 631
630	H=H+1
	K=2
631	IF (ABS(F).GT.EQMIN) GO TO 640
	F=0.0
	GO TO 641
640	H=H+1
	K=3
641	IF (H.NE. 1)  GO  TO  670
	IF (K.NE. 1)  GO  TO  650
C--REARRANGE INTO STANDARD MONOCLINIC SETTING.
	U=A
	A=B
	B=U
	E=ABS(D)
	D=0.0
	F=0.0
650	IF (K .NE.3)  GO TO 660
C--REARRANGE INTO STANDARD MONOCLINIC SETTING.
	U=A
	A=B
	B=C
	C=U
	E=ABS(F)
	D=0.0
	F=0.0
660	IF (E .LE.A) GO TO 670
	C=A+C-E
	E=A+A-E
C--ADJUST THE SIGNS OF TRICLINIC COSINE TERMS SO THAT ONLY QF CAN BE -VE
670	IF (D .GT.(-0.1))   GO TO 675
	D=-D
	E=-E
675	IF (E .GT. (-0.1) ) GO TO 680
	E=-E
	F=-F
680	IF (H.EQ.3 .AND. ABS(A+B+F-D-E) .LT. EQMIN ) GO TO 690
C--THIS IS A MONOCLINIC LATTICE IN DISGUISE.
	GO TO 700
690	E=(A-B)*0.5
	A=(A+B+F)/4.
	B=C-A
	C=A-F/2.
	D=0.0
	F=0.0
C  THIS WAS THE TRANSFORMATION OF CERTAIN MONOCLINIC CENTERED LATTICES.
700	CONTINUE
	IF (TEST) WRITE(IOUT,9935) (CEL(I),I=1,6)
1000	RETURN
9930	FORMAT(11H LATTICE IS, 6F8.1)
9931	FORMAT(11H REDUCED   , 6F8.1)
9934	FORMAT(11H REDUCED2  ,6F8.1)
9935	FORMAT(11H REDUCED3  , 6F8.1)
	END
	
	SUBROUTINE READER(NZ1,NZ2,NR,WMOL,DOBS)
C  ---------------------------------------------------------------------
C----THE FIRST CARD MUST CONTAIN THE NAME OF THE PROBLEM (TITLE).
C    THE NEXT CARD IS USUALLY COMPLETELY BLANK. THE PROGRAM
C    THEN TAKES THE DEFAULT VALUES.. IF YOU WANT TO SPECIFY YOUR
C    OWN PARAMETERS, YOU SHOULD DO THIS AS FOLLOWS0
C  COLUMN  FORMAT  NAME   DEFAULT   MEANING
C    1      A1     MAN      0       PRINTS INSTRUCTIONS,UNLESS MAN=9
C    2      A1     INSTR    0       PRINTS INSTRUCTIONS FOR SPECIAL
C                                   FEATURES IF INSTR=1
C    3      A1     INTENS   0       READS INTENSITIES BESIDES LINE
C                                   POSITIONS IF INTENS=1
C   5- 6    I2     NSYST(1) 0       ORTHORHOMBIC SOLUTION. +1=YES, -1=NO
C   7- 8    I2     NSYST(2) 0       MONOCLINIC SYSTEM.   0 MEANS  INDIFF
C   9-10    I2     NSYST(3) 0       TRICLINIC SYSTEM
C  11-15   F5.2    TOL2    3.0      TOLERANCE ON 2-DIMENSIONAL SEARCH
C  16-20   F5.2    TOL3    4.5      TOLERANCE ON 3-DIMENSIONAL SEARCH
C  21-30   F10.5   WAVEL   1.5405   WAVELENGTH
C  31-32    I2     LINCO    0       NUMBER OF GIVEN LINE COMBINATIONS
C                                   FOR THE 2-DIMENSIONAL SEARCH
C  33-34    I2     K (ENL)  0       ENL0 OUTPUT OF INTERMEDIATE RESULTS(
C  35-36    I2     NQ1      3       NQ1 AND NQ2 DETERMINE THE NUMBERS OF
C  37-38    I2     NQ2      6       LINES TO BE USED IN COMBINATIONS
C                                   FOR FINDING ZONES .
C  39-40    I2     NZ1      6       NZ1 AND NZ2 DO THE SAME FOR
C  41-42    I2     NZ2      6       FINDING LATTICES FROM ZONES.
C  43-44    I2     NR       0       REFINEMENT AND EVALUATION RUN ONLY
C                                   IF NR>0. THE DECK OF COMPLETE
C                                   LATTICES TO BE REFINED (6F10.2)
C                                   WILL FOLLOW THE LINES DECK, ALSO
C                                   TERMINATED BY A CARD WITH ZERO OR
C                                   BLANKS IN COLS 1-10.
C  45-46    I2     INP      5       NUMBER OF THE UNIT ON WHICH TO READ
C                                   THE LINE POSITIONS.
C  47-48    I2     LST      0       LIST OF CALCULATED LINES IF LST=1
C  49-50    I2     J (TEST) 0       TEST OUTPUT IF J>0
C  51-60   F10.5   WMOL    0.0      MOLECULAR WEIGHT.
C  61-70   F10.5   DOBS    0.0      OBSERVED DENSITY.
C  71-80   F10.5   TOLG    3.0      TOLERANCE ON MATCH BETWEEN CALCULATE
C                                   AND OBSERVED Q-VALUES.
C  ---------------------------------------------------------------------
C  ---------------------------------------------------------------------
	DIMENSION QIN(200),IIN(200),INTY(8),ISYST(3,3)
	LOGICAL SELF,ENL,TEST,LIST
CDC---LEVEL STATEMENTS FOR CDC 7600 FOLLOW---REMOVE FOR OTHER MACHINES--
CDC   LEVEL 2, LINCO,NQ2,NSOLMX,ROSTER
C---COMMON BLOCK /TEMRES/, LENGTH 11613 (OCTAL), 5003 (DECIMAL)
	COMMON /LINES /INP,IOUT,JOB,MAX,NAME(20),INTSY(40),OBS(40),CEL(10)
	COMMON /PARMS / TOL2,TOL3,WAVEL,FMINEQ,FMINIM,SELF,ENL,TEST,LIST,
     *                NSYST(3),OBSMT(40),OBSPT(40)
	COMMON /TEMRES/ LINCO,NQ2,NSOLMX,ROSTER(500,10)
	DATA ISYST(1,1)/4HORTH/,ISYST(1,2)/4HORHO/,ISYST(1,3)/4HMBIC/,
     1     ISYST(2,1)/4HMONO/,ISYST(2,2)/4HCLIN/,ISYST(2,3)/4HIC  /,
     2     ISYST(3,1)/4HTRIC/,ISYST(3,2)/4HLINI/,ISYST(3,3)/4HIC   /
	DATA NEND1/4HEND /,NEND2/4H END/
	DATA NINE/1H9/,IONE/1H1/,IBLA/3H   /
C  ---------------------------------------------------------------------
C  ---------------------------------------------------------------------
C  AFTER THE PARAMETER CARD WE MUST ENTER THE LINES,  EITHER AS
C  TWOTHETA, D- OR Q- VALUES.
C  THE PROGRAM WILL FIND OUT WHAT KIND OF VALUES IT HAS.
C  THE NUMBER OF LINES MUST BE AT LEAST 20  AND MUST NOT EXCEED 200.
C  A MAXIMUM OF 40 LINES IS ACTUALLY USED, THE REST IS HARMLESS.
C  ONLY THE FIRST SIX LINES SHOULD BE IN THEIR CORRECT ORDER, THE REST
C  CAN BE ENTERED IN A RANDOM WAY.  EVEN LINES THAT SHOULD BE AMONG
C  THE FIRST SIX CAN BE ENTERED LATER.
C  THE LINES CARDS SHOULD BE TERMINATED BY A BLANK CARD.
C  THE PROGRAM TESTS ON THIS  BLANK CARD AND THEN COUNTS THE NUMBER
C  OF REFLECTIONS.
C  ---------------------------------------------------------------------
9990	CONTINUE
C  ---------------------------------------------------------------------
C--IF YOU SET LINCO  UNEQUAL  ZERO YOU MUST NOW ADD LINCO CARDS.
C  THERE ARE SEVERAL POSSIBILITIES0
C  A. YOU ENTER TWO Q-VALUES (2F10.5) PER CARD.
C    IN THIS CASE THE PROGRAM TREATS THEM AS ONE OF HIS OWN
C    LINE COMBINATIONS AND TRIES TO FIND THE ANGLE IN FINZON  ETC.
C  B. YOU ENTER THREE CONSTANTS (3F10.5).
C     THE PROGRAM RECOGNIZES THIS AS A COMPLETE ZONE, BUT WILL
C     DETERMINE THE QUALITY OF THIS ZONE. THIS MAY LEAD TO AN
C     UNEXPECTED DOUBLING OF AXES.
C     NOTE0 A RECTANGULARZONE MUST BE ENTERED AS  A,B, 0.1  FOR READING
C     IN CORRECTLY.  THE PROGRAM WILL VERY SOON  REPLACE 0.1 BY 0.0
C  C. YOU ENTER 4 VALUES  (4F10.5).  THE FOURTH VALUE IS THE
C     QUALITY OF THE ZONE.  IF YOU SET THIS HIGH ENOUGH, THE ZONE
C     WILL CERTAINLY BE USED IN LATTICE-FINDING.
C     ENTERING A SINGLE DIGIT IN COLUMN 31 ASSURES YOU OF A HIGH QUALITY
C  ---------------------------------------------------------------------
C-----------------------------------------------------------------------
C
	IOUT=38
C---OUTPUT FROM =READER= IS DISCARDED BY DIVERTING IT TO CHANNEL 8
	CEL(10)=0.
	JOB=JOB+1
	CALL ECHO
	READ(INP,9901)  NAME
C     READ(INP,9901,END=9999)  NAME           (IBM ONLY)
C  FOR   IBM MACHINES THE TRICK WITH THE  END  CARD IS NOT NECESSARY.
	IF( NAME(1) .EQ.NEND1)GOTO 9999
	IF( NAME(1) .EQ.NEND2)GOTO 9999
	RAD=57.29578
	DO 6 I=1,8
6	INTY(I)=IBLA
	DO 7 I=1,200
7	IIN(I)=IBLA
	DO 8 I=1,40
8	INTSY(I)=IBLA
	TOLG=3.
	ENL=.TRUE.
	TEST=.TRUE.
	LIST=.FALSE.
	SELF=.FALSE.
	CALL ECHO
	READ(INP,9902) MAN,INSTR,INTENS,NSYST,TOL2,TOL3,WAVEL,LINCO,K,
     * NQ1,NQ2,NZ1,NZ2,NR,INP,LST,J,WMOL,DOBS,TOLG
C--PARAMETER CARD
	IF(NR.GT.0) GOTO 9
	IOUT=6
	WRITE(IOUT,9000)
9000	FORMAT(27H0NO SOLUTIONS TO BE REFINED/)
	GOTO 9999
9	CONTINUE
	WRITE(IOUT,9920) NAME
	WRITE(IOUT,9950)
	IF (TOL2.LT. 0.1) TOL2=3.0
	IF (TOL3.LT. 0.1) TOL3=4.5
	IF (WAVEL.LT.0.1) WAVEL=1.54050
	IF ( K   .EQ. 0) ENL =.FALSE.
	IF (NQ1 .EQ. 0) NQ1=3
	IF (NQ2 .EQ. 0) NQ2=6
	IF (NZ1 .EQ. 0) NZ1=6
	IF (NZ2 .EQ. 0) NZ2=6
	IF (INP .EQ. 0) INP=39
	IF (LST .EQ. 1) LIST=.TRUE.
	IF ( J   .EQ. 0) TEST=.FALSE.
	IF (TEST) SELF=.TRUE.
	IF (TOLG .LT. 0.1) TOLG=3.0
	L=0
10	CONTINUE
	CALL ECHO
	IF (INTENS .EQ. IONE) GO TO 12
	READ(INP,9903)  (CEL(J),J=1,8)
	GO TO 14
12	READ(INP,9912)  (  CEL(J),INTY(J)  ,J=1,8)
14	U=0.0
	DO 15 I=1,8
15	U=U+ABS(CEL(I) )
C  IT IS POSSIBLE TO GET RID OF A LINE BY SIMPLY PUNCHING A MINUS  SIGN
C  IN FRONT OF THE VALUE(PROVIDED THERE WAS ENOUGH SPACE IN THE FORMAT).
C  THIS TRICK NOW ALSO WORKS WITH ONLY ONE LINE PER CARD.
	IF (U .LT. 0.1) GO TO 30
	DO  20 I=1,8
	IF (CEL(I) .LT. 0.1) GOTO 20
	L=L+1
	QIN(L)  =CEL(I)
	IIN(L)=INTY(I)
20	CONTINUE
	GO TO 10
C-----------------------------------------------------------------------
C--END OF READING LINES CARDS
C-----------------------------------------------------------------------
30	L=L+1
	DO 40 I=1,8
	MAX=L-I
	IF (QIN(MAX) .GT. 0.1) GO TO 90
40	CONTINUE
	GO TO 90
90	CONTINUE
	WRITE(IOUT, 9918)  MAX
	WRITE(IOUT, 9915) ( QIN(M),IIN(M) ,M=1,MAX)
C-----------------------------------------------------------------------
C--WE NOW KNOW THE NUMBER OF LINES. NEXT WE TRY TO FIND OUT,
C--WHAT KIND OF VALUES WE HAVE  AND PRINT THEM OUT.
C-----------------------------------------------------------------------
	IF (MAX .GE. 40) MAX=40
	DO 150 I=1,8
150	CEL(I)=QIN(I)
	DO 160 M=1,5
	IF(CEL(M+1) .LT. CEL(M)) GO TO 180
160	CONTINUE
C-----------------------------------------------------------------------
C--IF INPUT IS IN DESCENDING ORDER, WE HAVE  D-VALUES
C-----------------------------------------------------------------------
	DO 170 M=1,MAX
	OBS(M)=QIN(M)
	INTSY(M)=IIN(M)
170	CONTINUE
	ISTYG=1
	CALL SRTRON(ISTYG)
	IF(OBS(MAX) .GE. 180.0) GO TO 290
C-----------------------------------------------------------------------
C--IF THE LAST INPUT VALUE IS GREATER THAN 180 WE DO NOT HAVE TWOTHETAS.
C-----------------------------------------------------------------------
	A=WAVEL/200.0
	B=1.0/(A*A)
	C=1.0/(2.0*RAD)
	WRITE(IOUT,9906)
	WRITE(IOUT,9909)
	DO 280 M=1,MAX
	W=OBS(M)
	D=SIN(C*W)
	OBS(M)=B*D*D
280	CONTINUE
	WRITE(IOUT,9911)  WAVEL
	GO TO 300
290	CONTINUE
	WRITE(IOUT,9921)  MAX
	WRITE(IOUT,9906)
	WRITE(IOUT,9908)
	GO TO 302
180	DO 200 M=1,MAX
	OBS(M)=QIN(M)
	INTSY(M)=IIN(M)
200	CONTINUE
	ISTYG=0
	CALL SRTRON(ISTYG)
	WRITE(IOUT,9906)
	WRITE(IOUT,9907)
	DO 210 M=1,MAX
	OBS(M)=10000.0/(OBS(M)*OBS(M))
210	CONTINUE
300	WRITE(IOUT,9921)  MAX
	WRITE(IOUT,9910)
302	WRITE(IOUT,9916)  ( OBS(M),INTSY(M) ,M=1,MAX)
	IF (MAX .GE. 20 ) GO  TO 305
	WRITE(IOUT, 9930) MAX
	GO TO 9990
305	CONTINUE
C-----------------------------------------------------------------------
C--THE ARRAY OF Q-VALUES HAS NOW BEEN ESTABLISHED.
C-----------------------------------------------------------------------
	ST=1.0E-3*TOLG
	IF (TEST) WRITE(IOUT,9936)
	DO 306 M=1,MAX
	W=OBS(M)
	T=TOLG
	IF (W .GT. 1000.) T=ST*W
	OBSMT(M)=W-T
	OBSPT(M)=W+T
	IF (TEST) WRITE(IOUT,9937) W,OBSMT(M),OBSPT(M)
306	CONTINUE
C--THE RANGE FOR MATCHING Q-VALUES HAS BEEN SET
	IF (NQ1 .LE. NQ2) GO TO 308
	I=NQ2
	NQ2=NQ1
	NQ1=I
308	IF (NZ1 .LE. NZ2) GO TO 309
	I=NZ2
	NZ2=NZ1
	NZ1=I
309	FMINIM=1.
C--- FMINIM HAS BEEN CHANGED FROM 0.1*OBS(1) IN ORDER TO INHIBIT
C--- THE TEST ON MINIMUM QA,QB,QC WHEN INCORPORATED INTO S/R FZONRF
	FMINEQ=OBS(20)*5.0E-4
	IF (FMINEQ .LT. 0.8)  FMINEQ=0.8
	A=4.0 *FMINEQ
	C=2.0*FMINEQ
	DO 310 M=2,MAX
	B=OBS(M)-OBS(M-1)
	IF  (B .GT. A) GO TO 310
	A=B
	IF ( B .GE. C) GO TO 310
	I=M-1
	D=RAD*WAVEL*B/( 200.0*SQRT(OBS(M)))
	WRITE(IOUT,9914) I,B,D
310	CONTINUE
	FMINEQ=0.25*A
	IF (TEST) ENL=.TRUE.
	IF (.NOT. ENL) WRITE(IOUT,9927)
	IF (ENL .AND. (.NOT. TEST)  ) WRITE(IOUT,9923)
	IF (TEST) WRITE(IOUT, 9924)
	WRITE (IOUT,9925) FMINIM
	WRITE(IOUT,9926) FMINEQ
	WRITE(IOUT,9922) TOL2,TOL3
	IF (ABS(TOLG-3.0) .GT. 0.01) WRITE(IOUT,9928) TOLG
	IF (NR.EQ.0) GOTO 315
C-----------------------------------------------------------------------
C---READ THE COMPLETE LATTICES (REFINEMENT RUN ONLY)
	WRITE(IOUT,9941)
	I=0
	IOUT=6
C---REINSTATE IOUT IN ORDER TO PRINT THE LIST OF SOLUTIONS AS INPUT
	WRITE(IOUT,9312)
9312	FORMAT(48H1** LIST OF CELLS TO BE REFINED AND EVALUATED **//)
311	CALL ECHO
	READ(INP,9940) (CEL(J),J=1,6)
	IF (CEL(1).EQ.0.) GOTO 314
	I=I+1
	IF (I.GT.NSOLMX) GOTO 313
	DO 312 J=1,6
312	ROSTER(I,J)=CEL(J)
	ROSTER(I,7)=99
	ROSTER(I,8)=I
313	IF (I.EQ.NSOLMX+1) WRITE(IOUT,9943) NSOLMX
	WRITE(IOUT, 9942) I,(CEL(J),J=1,6)
	GOTO 311
314	NR=I
	IF (NR.GT.NSOLMX) NR=NSOLMX
	WRITE(IOUT,9944) NR
315	CONTINUE
	IOUT=8
C---SUPPRESS OUTPUT AGAIN
500	CONTINUE
	DO 510 N=1,3
510	IF (NSYST(N) .EQ. 1) WRITE(IOUT,9952) (ISYST(N,I),I=1,3)
	WRITE(IOUT,9934) WMOL
	WRITE(IOUT,9935) DOBS
	IOUT=6
C---HAVE NOW FINISHED DISCARDING OUTPUT
	RETURN
9999	WRITE(IOUT,9998)
	NR=-1
	RETURN
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C  THE FOLLOWING CARDS ARE USEFUL FOR TEST PURPOSES0
C     WRITE(IOUT,9945) TOL2,TOL3,TOLG,WAVEL,LINCO,SELF,ENL,TEST,
C    * NQ1,NQ2,NZ1,NZ2 ,NR,INP,WMOL,DOBS
C9945 FORMAT(53H   TOL2  TOL3  TOLG  WAVEL   LINCO   SELF   ENL  TEST,
C    *33H NQ1 NQ2 NZ1 NZ2 NR INP WMOL DOBS/
C    *  3F6.1, F8.5, I8, 3L6,    6I4 ,2F8.5)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
9901	FORMAT(20A4)
9902	FORMAT(3A1,1X,3I2,2F5.2,F10.5,4I2,2I1,I4,3I2,3F10.5)
9903	FORMAT( 8F10.5 )
9906	FORMAT(23H0THE INPUT CONSISTED OF   )
9907	FORMAT(1H+, 23X,10H D-VALUES.)
9908	FORMAT(1H+, 23X,10H Q-VALUES.)
9909	FORMAT(1H+, 23X,12H TWO THETAS. )
9910	FORMAT(45H0THE INPUT LEADS TO THE FOLLOWING Q-VALUES   )
9911	FORMAT(20H0THE WAVELENGTH IS  , F8.4)
9912	FORMAT( 8(F7.3,A3) )
9914	FORMAT(40H0XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX      /
     * 41H THE DIFFERENCE IN Q-VALUES BETWEEN LINE ,  I3,
     * 22H AND THE NEXT IS ONLY  , F6.2 /
     * 14H THIS IS ONLY  ,F5.2,  18H DEGREES TWO THETA     /
     * 40H            PLEASE  CHECK  YOUR  INPUT     /
     * 56H THE PROGRAM WILL CONTINUE, BUT MAY NOT GIVE THE CORRECT,
     * 59H ANSWERS SINCE THE LIMIT FOR  'EQUAL' VALUES CAN BE TOO LOW /
     * 40H XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX      )
9915	FORMAT( 10(1X,F8.3,1X,A3)  )
9916	FORMAT( 10(1X,F7.1,2X,A3)  )
9917	FORMAT(15H INPUT ZONE NR., I3, 4F8.1)
9918	FORMAT(37H0THE INPUT CONSISTED OF THE FOLLOWING,I4,7H LINES   )
9919	FORMAT(3F10.2, F10.1)
9920  FORMAT(50H1THIS PROGRAM TRIES TO FIND THE UNITCELL FROM THE ,
     *18H POWDER PATTERN OF ///20X,20A4///26H WITH A METHOD PROPOSED BY,
     *61H  DE WOLFF AND PROGRAMMED BY VISSER (J.APPL.CRYST.2,89(1969))/
     *40H            FORTRAN  VERSION    6         )
9921	FORMAT(41H0THE NUMBER OF LINES THAT WILL BE USED IS ,I4)
9922	FORMAT(43H0(COL 11-15)THE TOLERANCE ON THE TRIAL F IS, F4.1//
     *44H (COL 16-20) THE TOLERANCE ON THE TRIAL D IS , F4.1)
9923	FORMAT(38H0INTERMEDIATE RESULTS ARE PRINTED OUT./
     * 58H PUNCH A 1 IN COLUMN 50 OF THE PARAMETER CARD IF YOU WANT  ,
     * 18H THE TEST OUTPUT.   )
9924	FORMAT(39H0TESTRUN, ALL TEST RESULTS ARE PRINTED.)
9925	FORMAT(44H0THE LOWER LIMIT FOR THE VALUES OF A,B,C IS ,  F8.2)
9926  FORMAT(66H0VALUES WILL BE CONSIDERED TO BE EQUAL WHEN THEY DIFFER
     * LESS THAN , F4.1 )
9927  FORMAT(55H0ONLY A MINIMUM OF INTERMEDIATE RESULTS WILL BE PRINTED/
     * 59H PUNCH A 1 IN COLUMN 34 OF THE PARAMETER CARD IF YOU WANT  ,
     * 14H MORE OUTPUT.  )
9928	FORMAT(48H ***********************************************/
     1       40H * THE OVERALL TOLERANCE IS CHANGED INTO,F6.2,2H */
     2       48H *********************************************** )
9930	FORMAT(32H0THE NUMBER OF USEFUL LINES  IS ,I4 /
     *60H SINCE THIS IS LESS THAN THE MINIMUM NUMBER OF LINES THAT TH ,
     *60HE PROGRAM NEEDS (20),YOU CANNOT EXPECT RELIABLE  RESULTS.   //
     *24H THIS PROBLEM IS SKIPPED  )
9931	FORMAT(21H0THE PROGRAM WILL USE ,I3,24H GIVEN LINE COMBINATIONS)
9932	FORMAT(35H0THE PROGRAM WILL COMBINE THE FIRST ,I3,
     * 21H LINES WITH THE FIRST  ,I3)
9933	FORMAT(42H0THE PROGRAM WILL TRY TO COMBINE THE FIRST  ,I3,
     * 21H ZONES WITH THE FIRST ,I3)
9934	FORMAT(30H0THE GIVEN MOLECULAR WEIGHT IS,  F10.3)
9935	FORMAT(30H0THE GIVEN OBSERVED DENSITY IS,  F10.4)
9936  FORMAT(63H0THE Q-VALUES WITH THEIR LIMITS OF ERROR,Q +/- TOLERANCE
     1 (TOLG)/ 32H       Q        Q-TOL     Q+TOL    )
9937	FORMAT(1H ,3F10.1)
9940	FORMAT( 8F10.2 )
9941  FORMAT(57H0INPUT OF COMPLETE LATTICES FOLLOWS (REFINEMENT RUN ONLY
     1)/)
9942	FORMAT(18H INPUT LATTICE NR., I3, 8F10.2 )
9943	FORMAT(1H0,40H** EXCEEDS STORAGE ALLOCATED, ONLY FIRST,I4,
     1 21H WILL BE PROCESSED **/)
9944	FORMAT(30H0LATTICE INPUT FINISHED: FIRST,I4,18H WILL BE PROCESSED)
9950  FORMAT(65H0*******************************************************
     1*********/  2H *,62X,1H*/ 65H *   ARE YOUR INPUT DATA REALLY THE B
     2EST THAT CAN BE OBTAINED   */2H *,62X,1H*/   65H *****************
     3***********************************************    ///
     480H FINDING THE UNIT CELL  DEPENDS FOR 95 PERCENT ON THE QUALITY O
     5F THE INPUT DATA./118H A RANDOM ERROR OF 0.03 DEGREES TWOTHETA CAN
     6 USUALLY BE TOLERATED, BUT A SYSTEMATIC (ZERO POINT) ERROR OF 0.02
     7 DEGREES/77H IS PROBABLY DISASTROUS.  CHECK YOUR INPUT LINES AGAIN
     8ST THEIR HIGHER ORDERS.//  83H DO NOT USE A DEBIJE-SCHERRER CAMERA
     9  UNLESS THE UNIT CELL TO BE EXPECTED IS SMALL./76H WOULD YOU LIKE
     O TO SOLVE A JIGSAW PUZZLE WHEN HALF THE PIECES ARE MISSING---)
9952  FORMAT(15H0LATTICES WITH , 3A4, 35H SYMMETRY WILL GET AN EXTRA CHA
     *NCE  )
9998	FORMAT(16H END OF PROGRAM   )
	END
	
	SUBROUTINE REFPAR(NSOL)
C-----------------------------------------------------------------------
C--REFPAR REFINES  THE PARAMETERS, PRINTS OUT THE LATTICE IF ENL=.TRUE.
C--AND PRINTS OUT A LIST OF   TWOTHETA,D,Q,H,K,L.
C-----------------------------------------------------------------------
	DIMENSION HULP(6),VAR(6),FYN(6),WE(10),WAR(6),QQSTOR(100),
     1 MATR1(6),MATR2(6),REF(6,6)
	DIMENSION KSTOR(20,5),LSTOR(20,5),QSTOR(20,5),INSTOR(100,3)
	INTEGER HSTOR(20,5)
	INTEGER TELOB(25),BR(21,5),P,  H,H1,HS,HMIN,HMAX,HMAX1,
     1 COUNT,ORDE
	LOGICAL SELF,ENL,TEST,LIST
	LOGICAL HEVEN,KEVEN,LEVEN,REPEAT,MIS,DRUK,PRIN,NOPRIN,LIST2
	LOGICAL FIRST4,IGOOD,MGOOD,ROUND2
	REAL M20MIN
	COMMON /LINES /INP,IOUT,JOB,MAX,NAME(20),INTSY(40),OBS(40),CEL(10)
	COMMON /PARMS / TOL2,TOL3,WAVEL,FMINEQ,FMINIM,SELF,ENL,TEST,LIST,
     *                NSYST(3),OBSMT(40),OBSPT(40)
	COMMON /LINSTR/ INSTOR,QQSTOR
	COMMON /REFP/ IREF
	EQUIVALENCE (INSTOR(1,1),HSTOR(1,1)),(INSTOR(1,2),KSTOR(1,1)),
     1  (INSTOR(1,3),LSTOR(1,1)),(QQSTOR(1),QSTOR(1,1) )
C-----------------------------------------------------------------------
	SWAVEL=WAVEL/200.0
	QMAX=OBS(20)*(1.0+0.002*TOL3)
	LIST2=.FALSE.
	ROUND2=IREF.EQ.0
	FIRST4=NSOL.LE.4
	I20MIN=15
	M20MIN=7.0
	IGOOD=CEL(7).GE.I20MIN
	MGOOD=CEL(8).GE.M20MIN
	IBRAV=6
	DRUK=.FALSE.
	DRUK= ROUND2 .AND. (FIRST4 .OR. MGOOD)
	IF (CEL(7) .LT. 2.0) GO TO 90
C--THE FIRST TIME THAT REFPAR IS CALLED, CEL(7)=0.
	IBRAV= CEL(9)
90	CONTINUE
	REPEAT=.TRUE.
	IVORDE=7
	PRIN= DRUK .AND. ENL .AND. MGOOD
	NOPRIN=.NOT. PRIN
	GO TO 100
95	IF ( (.NOT. DRUK) .AND. REPEAT)  CALL NORM3
100	NIET=0
	ISPAT=1
	ORDE=3
	DO 110 J=4,6
	IF (ABS(CEL(J)) .GT. 0.5) GO TO 105
	CEL(J)=0.
	GO TO 110
105	ORDE=ORDE+1
110	CONTINUE
	IF (IVORDE .EQ. ORDE) REPEAT=.NOT. REPEAT
	IVORDE=ORDE
	IF (ORDE .NE. 5) GO TO 200
	IF (ABS(CEL(6)) .LT. 0.1) GO TO 180
	IF (ABS(CEL(4)) .LT. 0.1) GO TO 160
	CALL WISSEL(2,3)
	GO  TO 180
160	CALL WISSEL(1,3)
180	IF (ABS(CEL(5)).LT. ABS(CEL(4))) CALL WISSEL(1,2)
200	CONTINUE
	DO 211 I=1,21
	DO 210 J=1,5
210	BR(I,J)=0
	TELOB(I)=0
211	CONTINUE
220	IF (DRUK) WRITE(IOUT,9916)  NAME,NSOL
	IF (PRIN) WRITE(IOUT,9926) (CEL(J),J=1,10)
	IF (TEST .AND. NOPRIN) WRITE(IOUT,9926) (CEL(J),J=1,10)
C
C                 START GENERATING  Q-VALUES
C
	NN=0
	ASQ=CEL(1)
	BSQ=CEL(2)
	CSQ=CEL(3)
	ALFA=CEL(4)
	BETA=CEL(5)
	GAMMA=CEL(6)
	HMAX=INT(SQRT(QMAX/ASQ))+1
	KMAX=INT(SQRT(QMAX/BSQ))+1
	LMAX=INT(SQRT(QMAX/CSQ))+1
	N20=0
	KEVEN=.FALSE.
	LMIN=-LMAX
	IF(ORDE .EQ. 3) LMIN=0
	HMIN=0
	IF(ORDE .GE. 5) HMIN=-HMAX
	IF (PRIN) WRITE(IOUT,9917)
	KMAX1=KMAX+1
C  K-LOOP STARTS HERE
	DO 1050 K1=1,KMAX1
	K=K1-1
	RK=K
	IF(NOPRIN) GO TO 520
	WRITE(IOUT,9918)  K
	ISPAT=1
520	C1=RK*RK*BSQ
	C2=RK*ALFA
	C3=RK*GAMMA
	KEVEN=.NOT. KEVEN
	HS=HMIN
	IF(K .EQ.  0) HS=0
	HEVEN=.TRUE.
	IF(2*(HS/2) .EQ. HS) HEVEN=.FALSE.
	HMAX1=HMAX-HS+1
C  H-LOOP STARTS HERE
	DO 1040 H1=1,HMAX1
	H=H1+HS-1
	RH=H
	C4=C1+RH*(C3+RH*ASQ)
	C5=C2+RH*BETA
	IF(NOPRIN) GO TO 530
	WRITE(IOUT,9923)
	ISPAT=1
	IF((K+IABS(H)) .EQ. 0 .AND. ORDE .GT. 4) ISPAT=ISPAT+LMAX
	IF(ORDE .EQ. 4 . AND. H .EQ. 0) ISPAT=ISPAT+LMAX
530	HEVEN=.NOT. HEVEN
	LB=LMIN
	IF(H .NE. 0) GO TO 540
	IF(K .EQ. 0.OR. ORDE .EQ. 4) LB=0
540	LEVEN=.TRUE.
	IF(2*(LB/2) .EQ. LB) LEVEN=.FALSE.
	LMAX1=LMAX+1-LB
C  L-LOOP STARTS HERE
	DO 1030 L1=1,LMAX1
	L=L1-1+LB
	LEVEN=.NOT. LEVEN
	Q=C4+FLOAT(L)*(C5+FLOAT(L)*CSQ)
	IF (Q .GT. 1.0) GO TO 630
	IF (PRIN) CALL WRISP(ISPAT,0. ,0)
	GO TO 1020
630	IF(Q .GE. QMAX) GO TO 1020
	IF (.NOT. DRUK)      GO TO 700
	GO TO (640,650,660,670,680,700,950),IBRAV
640	IF(KEVEN.AND.LEVEN.OR..NOT.KEVEN.AND..NOT.LEVEN) GO TO 700
	GO TO 1020
650	IF(HEVEN.AND.LEVEN.OR..NOT.HEVEN.AND..NOT.LEVEN) GO TO 700
	GO TO 1020
660	IF(HEVEN.AND.KEVEN.OR..NOT.HEVEN.AND..NOT.KEVEN) GO TO 700
	GO TO 1020
670	IF (2*((H+K+L)/2) .EQ. H+K+L) GO TO 700
C--I-LATTICE
	GO TO 1020
680	IF(HEVEN.AND.KEVEN.OR..NOT.HEVEN.AND..NOT.KEVEN) GO TO 690
	GO TO 1020
690	IF(KEVEN.AND.LEVEN.OR..NOT.KEVEN.AND..NOT.LEVEN) GO TO 700
	GO TO 1020
700	N20=N20+1
	IF (.NOT. LIST2) GO TO 720
	NN=NN+1
	INSTOR(NN,1)=H
	INSTOR(NN,2)=K
	INSTOR(NN,3)=L
	QQSTOR(NN)=Q
	GO TO 1030
720	CONTINUE
	IF (PRIN) CALL WRISP(ISPAT,Q,0)
	MIS=.TRUE.
	DO 800 J=1,20
	IF (OBSPT(J) .LT. Q ) GO TO 800
	IF (OBSMT(J) .GT. Q ) GO TO 950
	P=J
	IF (NOPRIN) GO TO 750
	CALL WRISP(ISPAT,Q,1)
	MIS=.FALSE.
	GO TO 780
750	CONTINUE
	IF(KEVEN.AND.LEVEN.OR..NOT.KEVEN.AND..NOT.LEVEN)
     * BR(P,1)=BR(P,1)+1
	IF(HEVEN.AND.LEVEN.OR..NOT.HEVEN.AND..NOT.LEVEN)
     * BR(P,2)=BR(P,2)+1
	IF(HEVEN.AND.KEVEN.OR..NOT.HEVEN.AND..NOT.KEVEN)
     * BR(P,3)=BR(P,3)+1
	IF (2*((H+K+L)/2) .EQ. H+K+L) BR(P,4)=BR(P,4)+1
	IF ((HEVEN.AND.KEVEN.OR..NOT.HEVEN.AND..NOT.KEVEN)  .AND.
     * (KEVEN.AND.LEVEN.OR..NOT.KEVEN.AND..NOT.LEVEN))BR(P,5)=BR(P,5)+1
780	I=TELOB(P)+1
	IF (I .GT.  5) GO TO 800
	TELOB(P)=I
	HSTOR(P,I)=H
	KSTOR(P,I)=K
	LSTOR(P,I)=L
	QSTOR(P,I)=Q
800	CONTINUE
950	CONTINUE
1020	IF (PRIN) ISPAT=ISPAT+1
C  END OF K,H,L-LOOP ON 1050,1040,1030 RESPECTIVELY
1030	CONTINUE
	IF (PRIN) CALL WRISP(ISPAT,Q,2)
1040	CONTINUE
1050	CONTINUE
C
C             END GENERATING Q
C
	IF (.NOT. LIST2) GO TO 1100
	CALL LISTM(NN)
	 WRITE(IOUT,9941)  (CEL(J),J=1,6),N20,NAME,NSOL
	GO TO  1360
1100	DO 1110 H=1,6
	DO 1105 K=1,6
1105	REF(H,K)=0.0
1110	CONTINUE
	DO 1120 H=1,6
	FYN(H)=0.0
	WAR(H)=0.0
1120	CONTINUE
	T=0.0
	RMS=0.0
	SUMEPS=1.E-3
C  WHEN THE NUMBER OF LINES INDEXED  EQUALS THE NUMBER OF PARAMETERS,
C  SUMEPS WILL KEEP ITS STARTING VALUE.  THE VALUE ZERO CAUSES ACCIDENTS
	COUNT=0
	IF (DRUK)  WRITE(IOUT,9910)
	DO 1180 P=1,20
	NR=TELOB(P)
	QOBS=OBS(P)
	IF (.NOT. DRUK) GO TO 1122
	W=SQRT(QOBS)
	U=SWAVEL*W
	W=100.0/W
	U=114.59156*ATAN(U/SQRT(1.0-U*U))
	WRITE(IOUT,9911) U,W,QOBS,INTSY(P)
1122	CONTINUE
	IF (NR .GT. 0) GO TO 1124
	IF (DRUK) WRITE(IOUT,9912)
	NIET=NIET+1
	GO TO 1180
1124	CONTINUE
	EPS=1000.0
	SUMW=0.0
	DO 1130 N=1,NR
	D=ABS(QOBS-QSTOR(P,N))
	IF(D .LT. EPS)EPS=D
	IF (D .GT. 9.0) D=9.0
C--NECESSARY TO PREVENT UNDERFLOW ON  PDP-10 (AND UNIVAC?)
	W=EXP(-D*D)
	WE(N)=W
	SUMW=SUMW+W
1130	CONTINUE
	SUMEPS=SUMEPS+EPS
	DO 1170 N=1,NR
	H=HSTOR(P,N)
	K=KSTOR(P,N)
	L=LSTOR(P,N)
	Q=QSTOR(P,N)
	IF (.NOT. DRUK) GO TO 1042
	W=SQRT(Q)
	U=SWAVEL*W
	W=100.0/W
	U=114.59156*ATAN(U/SQRT(1.0-U*U))
	WRITE(IOUT,9914) U,W,Q,H,K,L
1042	CONTINUE
	DIF=QOBS-Q
	COUNT=COUNT+1
	W=WE(N)/SUMW
	RMS=RMS+W*DIF*DIF
	T=T+W
	VAR(1)=H*H
	VAR(2)=K*K
	VAR(3)=L*L
	VAR(4)=H*L
	VAR(5)=K*L
	VAR(6)=H*K
	DO 1150 J=1,ORDE
	HULP(J)=W*VAR(J)
	WAR(J)=WAR(J)+W*DIF*VAR(J)
1150	CONTINUE
	DO 1160 J=1,ORDE
	DO 1155 I=1,ORDE
1155	REF(I,J)=REF(I,J)+VAR(I)*HULP(J)
1160	CONTINUE
1170	CONTINUE
1180	CONTINUE
	CALL ARRAY(2,ORDE,ORDE,6,6,REF,REF)
	CALL MINV(REF,ORDE,DETER,MATR1,MATR2)
	CALL ARRAY(1,ORDE,ORDE,6,6,REF,REF)
	IF (ABS(DETER) .GE. 0.001) GO TO 1185
	IF (TEST) WRITE(IOUT,9931)
	GO TO 1210
1185	IF(T .GT. FLOAT(ORDE+1)) RMS=RMS/(T-FLOAT(ORDE))
	DO 1195 I=1,ORDE
	DO 1190 J=1,ORDE
1190	FYN(I)=FYN(I)+REF(I,J)*WAR(J)
1195	CONTINUE
	DO 1200 N=1,3
1200	CEL(N)=CEL(N)+FYN(N)
	IF (ORDE .GT. 3) CEL(5)=CEL(5)+FYN(4)
	IF (ORDE .GT. 4) CEL(4)=CEL(4)+FYN(5)
	IF (ORDE .GT. 5) CEL(6)=CEL(6)+FYN(6)
1210	NIET=20-NIET
	IF (DRUK) GO TO 1350
	IF (REPEAT) GO TO 95
	DO 1320 J=1,5
	N=0
	DO 1310 M=1,20
1310	IF (BR(M,J) .GT. 0 ) N=N+1
	BR(21,J)=N
1320	CONTINUE
	IBRAV=6
	DO 1330 J=1,5
	IF ( BR(21,J) .GE. NIET) IBRAV=J
1330	CONTINUE
1350	CEL(7)=NIET
	CEL(9)=IBRAV
	CEL(8)=FLOAT(NIET)*OBS(20)/(FLOAT(2*N20)*SUMEPS)
	IF (PRIN) WRITE(IOUT,9940) (CEL(J),J=7,10)
	IF (DRUK) WRITE(IOUT,9941)  (CEL(J),J=1,6),N20,NAME,NSOL
	LIST2=LIST .AND. DRUK
	LIST2=LIST2.AND.FIRST4.AND.MGOOD
	PRIN=.FALSE.
	NOPRIN=.TRUE.
	IF (LIST2 .AND. (N20 .LE. 100) ) GO TO 220
	IF (LIST2) WRITE(IOUT,9942) N20
1360	RETURN
9910	FORMAT(41H0TWOTHETA     D       Q       H    K    L)
9911	FORMAT(1H0, F8.2, F8.3, F9.1, 18H    *OBSERVED**    ,A3)
9912	FORMAT(1H+,45X ,  12H NOT INDEXED)
9914	FORMAT(1H ,F8.2,F8.3,F9.1, 3I5)
9916	FORMAT(1H1,20A4,14H   SOLUTION NR,I3)
9917	FORMAT(36H0LATTICE PRINTED OUT FOR INSPECTION //
     156H DISPLAYED IN LAYERS UP --K-- ,    --L-- ACROSS THE PAGE,
     222H --H-- DOWN THE PAGE.   //)
9918	FORMAT( 8H ***  H ,I2,8H  L  ***   // )
9923	FORMAT(1H0)
9926	FORMAT(1H ,6F9.2, F9.0,F9.1, 2F9.0 )
9931	FORMAT(1H ,14X,19HMATRIX IS SINGULAR )
9940	FORMAT(1H0//14H LINES INDEXED,F5.0,16H FIGURE OF MERIT, F8.1,
     *13H LATTICE TYPE, F6.0, 9H PEDIGREE,F8.0)
9941  FORMAT(1H0//46H RECIPROCAL LATTICE CONSTANTS AFTER REFINEMENT/
     * 6F9.2, 5H N20=,I4 ///1H ,20A4,14H   SOLUTION NR, I3)
9942	FORMAT(35H0THE NUMBER OF CALCULATED LINES IS  ,I4,
     *50H. NO LISTING IS PRINTED WHEN N20 IS MORE THAN 100   /
     * 30H SEE END OF SUBROUTINE REFPAR  )
	END
	
	SUBROUTINE SRTRON(ISTYG)
C
C
C     THIS SUBPROGRAM SORTS THE TWOTHETA-,D- OR Q-VALUES (POSSIBLY
C     TOGETHER WITH INTENSITIES).THE TWOTHETA- OR Q-VALUES ARE PUT IN
C     INCREASING ORDER,THE D-VALUES IN DECREASING ORDER.
C
C     ISTYG=1 - DATA ARE TWOTHETA- OR Q-VALUES (TO PUT IN INCREASING
C               ORDER).
C     ISTYG=0 - DATA ARE D-VALUES (TO PUT IN DECREASING ORDER).
	COMMON /LINES /INP,IOUT,JOB,MAX,NAME(20),INTSY(40),OBS(40),CEL(10)
C
	FAKTOR=1.0
	IF(ISTYG.EQ.0)FAKTOR=-1.0
20	IEIND=MAX-1
C
	DO 70 I=1,IEIND
	IF(FAKTOR*(OBS(I+1)-OBS(I)))30,70,70
30	OBSTR=OBS(I+1)
	IINT=INTSY(I+1)
	DO 40 J=1,I
C
	IF(FAKTOR*(OBSTR-OBS(J)))50,40,40
40	CONTINUE
C
50	MEIND=I+1-J
C
	DO 60 M=1,MEIND
	K=I+1-M
	OBS(K+1)=OBS(K)
	INTSY(K+1)=INTSY(K)
60	CONTINUE
	OBS(J)=OBSTR
	INTSY(J)=IINT
70	CONTINUE
	DO 80 I=1,IEIND
	IF (ABS(OBS(I+1)-OBS(I) ) .GT. 0.0001) GO TO 80
	WRITE(IOUT,99) OBS(I)
	OBS(I+1)=OBS(MAX)
	INTSY(I+1)=INTSY(MAX)
	MAX=MAX-1
	GO TO 20
80	CONTINUE
C
	RETURN
99	FORMAT(41H0XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX   /
     * 17H THE INPUT VALUE  ,F10.5, 15H APPEARS TWICE   /
     * 27H           CHECK YOUR INPUT    /
     * 41H XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX   )
	END
	
	SUBROUTINE  SYMTST
C  SUBROUTINE SYMTEST TRIES TO FIND OUT WHETHER THE LATTICES FOUND
C  ARE LATTICES OF A HIGHER SYMMETRY IN DISGUISE.
	DIMENSION NAM(3,2)
	INTEGER  H
	LOGICAL SELF,ENL,TEST,LIST
	COMMON /LINES /INP,IOUT,JOB,MAX,NAME(20),INTSY(40),OBS(40),CEL(10)
	COMMON /PARMS / TOL2,TOL3,WAVEL,FMINEQ,FMINIM,SELF,ENL,TEST,LIST,
     *                NSYST(3),OBSMT(40),OBSPT(40)
	EQUIVALENCE(CEL(1),A),(CEL(2),B),(CEL(3),C),(CEL(5),E)
	DATA NAM/4HIS P,4HROBA,4HBLY ,4H   M,4HIGHT,4H BE /
	NPASS=0
	DIFMIN=FMINEQ
	IF (CEL(4)+ABS(CEL(6)) .GT. FMINEQ) GO TO 3010
5	DIFMIN=3.0*DIFMIN
	NPASS=NPASS+1
	IF ( E .LT. FMINEQ) GO TO 1000
C  THE NEXT PART   FINDS OUT WHETHER A MONOCLINIC LATTICE IS A DISGUISED
C  RHOMBOHEDRIC LATTICE.
	AA=A
	IF (ABS(E-2.*B) .GT. DIFMIN) GO TO 20
	IF (ABS(C-3.*B) .GT. DIFMIN) GO TO 10
	GO TO 200
10	IF (ABS(A-3.*B) .GT. DIFMIN) GO TO 20
	AA=C
	GO TO 200
20	IF (ABS(A+B-C) .GT. DIFMIN) GO TO 30
	IF (ABS(2.*ABS(B-A)-E).GT.DIFMIN) GO TO 30
	GO TO 200
30	IF (ABS(B+C-4.*A).GT. DIFMIN) GO TO 40
	IF (ABS(E-ABS(B-C)) .GT. DIFMIN) GO TO 40
	GO TO 200
40	IF (ABS(2.*C-3.*E).GT.DIFMIN) GO TO 50
	IF (ABS(2.*B+E-6.*A).GT. DIFMIN*FLOAT(NPASS) ) GO TO 50
	GO TO 200
50	IF (ABS(2.*A-3.*E) .GT. DIFMIN) GO TO 60
	IF (ABS(2.*B+E-6.*C).GT. DIFMIN*FLOAT(NPASS) ) GO TO 60
	AA=C
	GO TO 200
60	IF (NPASS .EQ. 1) GO TO 5
	GO TO 3010
200	X=B/3.0
	Z=AA-X
	IF (Z .LT. 1.) RETURN
	WRITE(IOUT,900)  (CEL(J),J=1,6),(NAM(J,NPASS),J=1,3)
	A=200./SQRT(B)
	C=100./SQRT(Z)
	WRITE(IOUT,901) X,Z, A, C
	X=4.*X
	A=0.5*A
	IF (CEL(9).LT. 6.) WRITE(IOUT,902) X,A
	RETURN
1000	DO 1200 I=1,2
	J=3-I
	G=CEL(J)
	H=INT(A/G+0.5 )
	K=INT(B/G+0.5 )
	L=INT(C/G+0.5 )
	IH=0
	IK=0
	IL=0
	IF (ABS(A-G*FLOAT(H)) .LT. DIFMIN) IH=1
	IF (ABS(B-G*FLOAT(K) ) .LT. DIFMIN) IK=1
	IF (ABS(C-G*FLOAT(L) ) .LT. DIFMIN) IL=1
	II=IH+IK+IL -2
	IF (II) 1200,1020,1010
1010	WRITE(IOUT,900)  (CEL(J),J=1,6),(NAM(J,NPASS),J=1,3)
	X=100./SQRT(G)
	WRITE(IOUT, 911) G,X
	RETURN
1020	IF (IK*IL .EQ. 0) GO TO 1050
	N=L/K
	X=B
	Z=A
	GO TO 1500
1050	IF (IH*IL .EQ. 0) GO TO 1100
	N=L/H
	X=A
	Z=B
	GO TO 1500
1100	N=H/K
	X=B
	Z=C
	GO TO 1500
1200	CONTINUE
C     IF (NPASS .EQ. 1) GO TO 5
	GO TO 3010
1500	WRITE(IOUT,900)  (CEL(J),J=1,6),(NAM(J,NPASS),J=1,3)
	A=100./SQRT(X)
	C=100./SQRT(Z)
	IF (N-3) 1025,1030,1040
1025	WRITE(IOUT,912) X,Z, A,C
	RETURN
1030	A=2.0*A/SQRT(3.0)
	WRITE(IOUT,913) X,Z,A,C
	RETURN
1040	WRITE(IOUT,914) N
	RETURN
3010	AA=OBS(MAX)
	Z=0.
	DIFMIN=FMINEQ
	DO 3020 J=1,6
	X=ABS(CEL(J))
	IF (X .LT. FMINEQ) GO TO 3020
	IF (X .LT. AA) AA=X
C  FINDS THE SMALLEST NON-ZERO VALUE OF THE CONSTANT, AA.
	IF (X .GT. Z ) Z=X
C  AND THE LARGEST VALUE,  Z.
3020	CONTINUE
	DO 3050 I=1,3
	X=AA/FLOAT(I)
	IF (X .LT. 0.5*FMINIM) RETURN
	L=INT(Z/X +0.5)
	X=Z/FLOAT(L)
C  REFINES THE VALUE OF THE CONSTANT BY USING THE LARGEST MULTIPLE.
	N=0
	DO 3030 J=1,6
	G=ABS(CEL(J))
	L=INT(G/X+0.5)
	IF (ABS(G-X*FLOAT(L)) .LE. DIFMIN) N=N+1
3030	CONTINUE
	IF (N .LT. 6) GO TO 3050
	GO TO 3100
3050	CONTINUE
	RETURN
3100	WRITE(IOUT,900) (CEL(J),J=1,6),(NAM(J,2),J=1,3)
	Z=100./SQRT(X)
	WRITE(IOUT,911) X,Z
	RETURN
900   FORMAT(52H0THE LATTICE WITH THE RECIPROCAL LATTICE CONSTANTS0 ,
     * 6F9.2 /1H ,3A4 )
901	FORMAT(1H+, 12X, 24H RHOMBOHEDRAL, WITH THE   ,
     *45H (HEXAGONAL) RECIPROCAL LATTICE CONSTANTS0 X=,F9.2,4H  Z=,F9.2/
     * 29H AND THE DIRECT CONSTANTS  A=, F7.4,4H  C=,F7.4 )
902   FORMAT(21H PROBABLY X SHOULD BE, F10.2,8H  AND A0  ,F8.4)
911   FORMAT(1H+,12X,44H CUBIC WITH A RECIPROCAL LATTICE CONSTANT OF  ,
     *        F10.2,33H AND A DIRECT LATTICE CONSTANT OF , F10.4)
912   FORMAT(1H+,12X,42H TETRAGONAL WITH THE RECIPROCAL CONSTANTS  ,
     * 2F10.2,29H AND THE DIRECT CONSTANTS  A=, F7.4,4H  C=,F7.4 )
913   FORMAT(1H+,12X,41H HEXAGONAL WITH THE RECIPROCAL CONSTANTS   ,
     * 2F10.2,29H AND THE DIRECT CONSTANTS  A=, F7.4,4H  C=,F7.4 )
914   FORMAT(1H+,12X,43H A FREAK. THERE IS A RATIO IN THE CONSTANTS   ,
     * 32H OF THE RECIPROCAL  LATTICE OF0  ,I3)
	END
	
	SUBROUTINE UNTCEL(WMOL,DOBS,VOLUM,ISW)
C-----------------------------------------------------------------------
C  THIS SUBROUTINE CALCULATES THE DIRECT CONSTANTS FROM THE CONSTANTS
C  OF THE RECIPROCAL LATTICE.
C  IF THE MOLECULAR WEIGHT AND THE DENSITY HAVE BEEN GIVEN,
C  IT ALSO CALCULATES THE X-RAY DENSITY,DX,  AND THE NUMBER OF MOLECULES
C  PER UNIT CELL,  N.
C-----------------------------------------------------------------------
	DIMENSION          REC(3,6),REEL(3,6),WAAR(6)
	COMMON /LINES /INP,IOUT,JOB,MAX,NAME(20),INTSY(40),OBS(40),CEL(10)
C
	DO 1100 J=1,6
1100	WAAR(J)= CEL(J)
	DO 1220 J=1,3
	REC(1,J)=SQRT(WAAR(J))
	REC(1,J+3)=REC(1,J)
1220	CONTINUE
	DO 1230 J=1,3
	W=WAAR(J+3)/(REC(1,J+1)*REC(1,J+2)*2.0)
	REC(2,J+3)=W
	REC(2,J)=W
	REC(3,J+3)=SQRT(1.0-W*W)
	REC(3,J)=REC(3,J+3)
1230	CONTINUE
	DO 1240 J=1,3
	W=(REC(2,J+1)*REC(2,J+2)-REC(2,J))/(REC(3,J+1)*REC(3,J+2))
	REEL(2,J+3)=W
	REEL(2,J)=W
	REEL(3,J)=SQRT(1.0-W*W)
	REEL(3,J+3)=REEL(3,J)
1240	CONTINUE
	FNNN=REEL(3,1)*REEL(3,2)*REC(3,3)
	DO 1250 J=1,3
	CEL(J)=REEL(3,J)/(FNNN*REC(1,J))*100.0
1250	CONTINUE
	DO 1260 J=4,6
	CEL(J)=90.0-57.29578*ATAN(REEL(2,J)/REEL(3,J))
1260	CONTINUE
	VOLUM=FNNN*CEL(1)*CEL(2)*CEL(3)
	IF (ISW .EQ. 0) RETURN
	WRITE(IOUT,9930) (CEL(J),J=1,6) ,VOLUM
	IF (WMOL .LT. 0.1) RETURN
	PARTD=WMOL/(0.6025*VOLUM)
	IF (DOBS .LT. 0.1) GO TO 1270
	N=IFIX(1.1*DOBS/PARTD)
	DX=PARTD*FLOAT(N)
	WRITE(IOUT,9931)  N,DX,DOBS
	RETURN
1270	WRITE(IOUT,9932) PARTD
	RETURN
9930	FORMAT(1H ,6F9.3,  F12.2)
9931	FORMAT(1H+,68X, 5H   Z=,I2,6H   DX=, F7.4, 8H   DOBS=, F7.4 )
9932	FORMAT(1H+,68X,11H   DCALC=Z*, F8.4)
	END
	
	SUBROUTINE WISSEL(I,J)
	COMMON /LINES /INP,IOUT,JOB,MAX,NAME(20),INTSY(40),OBS(40),CEL(10)
	A=CEL(I)
	CEL(I)=CEL(J)
	 CEL(J)=A
	A=CEL(I+3)
	CEL(I+3)=CEL(J+3)
	CEL(J+3)=A
	RETURN
	END
	
	SUBROUTINE WRISP(NPOS,VAR,ICODE)
C  THE FUNCTION OF  WRISP  IS TO WRITE THE Q-VALUES OF THE CALCULATED
C  LINES IN A NEAT SCHEME, RECOGNISABLE FOR THOSE WHO  KNOW HOW TO
C  HANDLE THESE SCHEMES.
C     IT WRITES  VAR  (THE Q-VALUE) IN THE FORMAT F8.1, PRECEDED BY
C  (NPOS-1)*9 SPACES.   IF THE CALCULATED VALUE OF Q CORRESPONDS TO
C  AN OBSERVED VALUE, AN X IS ADDED LATER.
C
C  THE PARAMETERS OF WRISP ARE0
C  1. NPOS     = POSITION OF Q-VALUE
C  2. VAR      = VALUE OF Q
C  3. ICODE    = CODE TO TELL WRISP WHAT  TO DO
C  ICODE=0  PUT VAR IN POSITION
C       =1  PUT X IN POSITION.
C       =2  PRINT THE LINE WITH Q-VALUES.
C  THIS PARTICULAR FORM OF WRISP TAKES A LITTLE BIT OF COMPUTING TIME
C  BUT IT IS KIND TO THE PRINTER.  IT SHOULD BE USED WHEN YOU HAVE A
C  RELATIVELY FAST COMPUTER AND A SLOW PRINTER (E.G. A SLOW TERMINAL
C  ON A FAST COMPUTER).
	DIMENSION IFIGRS(10),ITXT(133)
	COMMON /LINES /INP,IOUT,JOB,MAX,NAME(20),INTSY(40),OBS(40),CEL(10)
	DATA IBLANK/1H /,IX/1HX/,IPNT/1H./,IFIGRS(1)/1H0/,IFIGRS(2)/1H1/,
     * IFIGRS(3)/1H2/,IFIGRS(4)/1H3/,IFIGRS(5)/1H4/,IFIGRS(6)/1H5/,
     2 IFIGRS(7)/1H6/,IFIGRS(8)/1H7/,IFIGRS(9)/1H8/,IFIGRS(10)/1H9/
	DATA NUSED/0/,ISTAR/1H*/
	IF (NUSED .GT. 0) GO TO 20
	GO TO 400
20	NUSED=NUSED+1
	IF (ICODE-1) 100,200,300
C
100	IF (NPOS .LE. 14) GO TO 150
C  ONLY 132 PRINT POSITIONS AVAILABLE,  15*9=135!
	WRITE(IOUT,9990)ITXT
	NPOS=1
	DO 110 K=1,133
110	ITXT(K)=IBLANK
150	N=9*(NPOS-1)
	IF (VAR .GT. 0.1) GO TO 156
	DO 152 I=3,9
152	ITXT(N+I)=ISTAR
	ITXT(N+6)=IFIGRS(1)
	GO TO 180
156	NR=IFIX(VAR*10.0+0.5)
	IREST=NR/10
	NFIG=NR-10*IREST
	ITXT(N+8)=IFIGRS(NFIG+1)
	ITXT(N+7)=IPNT
	NR=IREST
	DO 160 I=1,6
	INR=7-I
	IREST=NR/10
	NFIG=NR-10*IREST
	NR=IREST
	ITXT(N+INR)=IFIGRS(NFIG+1)
	IF (NR .EQ. 0) GO TO 180
160	CONTINUE
180	RETURN
C
200	N=9*NPOS
	ITXT(N)=IX
	RETURN
C
300	WRITE(IOUT,9990)ITXT
	NPOS=1
C
400	ITXT(1)=IFIGRS(1)
	DO 410 K=2,133
410	ITXT(K)=IBLANK
	IF (NUSED .EQ. 0) GO TO 20
	RETURN
9990	FORMAT(133A1)
	END
	
