C      LATTICE PARAMETER REFINEMENT
C
C
C      THIS PROGRAM IS OF UNCERTAIN ORIGINS.  IT WAS MODIFIED FOR 
C      A DEC-10 COMPUTER AT DUPONT BY K.H.GARDNER.  PROGRAM MADE
C      INTERACTIVE 10/78 BY KHG.  DENSITY OPTION ADDED 1/79
C      BY KHG.
C
	PROGRAM CELLB
C      UNIT CELL REFINEMENT PROGRAM
	REAL    CELN(6)                 ! CELL* PARAMETERS
	REAL    CELL(6)                 ! CELL PARAMETERS
	INTEGER HH(201)                  ! MILLER INDEX
	INTEGER KK(201)                  ! MILLER INDEX
	INTEGER LL(201)                  ! MILLER INDEX
        REAL    LAMBDA                   ! WAVELENGTH
	INTEGER NOT(6)			! KEY FOR FIXED PARAMETERS
        INTEGER NREL                    ! NUMBER OF REFLECTIONS
        REAL    RA                      ! A*
        REAL    RALP                    ! ALPHA*
        REAL    RB                      ! B*
        REAL    RBET                    ! BETA*
        REAL    RC                      ! C*
        REAL    RGAM                    ! GAMMA*
	REAL	SHT(6)			! LQ SHIFT IN CELL PARAMETER
        REAL    SIG(201)                ! WEIGHT/ERROR IN INPUT; ERROR
	CHARACTER*80 TITLE
        REAL    TTH(201)                ! TWO THETA/D IN INPUT; TWO THETA
        REAL    WTG(201)                ! REFLECTION WEIGHT
	DIMENSION A(7,7),SSQ(201),SSC(201),D(7),B(6,6),
	1	NUMTMP(2,3),TEM(6),COSINE(3),SINE(3),CELERR(6)
      EQUIVALENCE (CELERR,SHT)
      COMMON RA,RB,RC,RALP,RBET,RGAM,CELN
      PI=3.1415926536
      NPAR=6
C      OPEN FILE (PREPARED BY PROGRAM CELLA) THAT CONTAINS 
C      REFLECTION INFORMATION.
      OPEN(UNIT=1,FILE='CELL.DAT',TYPE='OLD') 

      OPEN(UNIT=2,FILE='CELL.OUT',TYPE='NEW')
      MM=7
      WRITE(2,74)
74    FORMAT(' REFLECTION FILE TITLE:')
C      READ TITLE FROM REFLECTION FILE
      WRITE(6,74)
      READ(1,72) TITLE
72    FORMAT(A)
      WRITE(6,73) TITLE
73    FORMAT(1H0,A/)
      WRITE(2,73) TITLE
C      READ TITLE FOR RUN FROM INTERACTIVE TERMINAL.
      WRITE(6,75)
75    FORMAT(' TITLE - '$)
      READ (5,72) TITLE
      WRITE (2,73) TITLE
      WRITE(6,76)
76    FORMAT(' ENTER APPROXIMATE UNIT CELL PARAMETERS')
      READ  (5,*) RA,RB,RC,RALP,RBET,RGAM
C
C      REFINEMENT IS PERFORMED IN RECIPROCAL SPACE.  INPUT CAN BE REAL
C	OR RECIPROCAL CELL PARAMETERS.
C      
      WRITE(6,77)
77    FORMAT(' IS UNIT CELL INPUT IN RECIPROCAL UNITS? '$)
	CALL YESNO(IKEY)
C      CHECK TO SEE IF INPUT IS REAL OR RECIPROCAL PARAMETERS
	IF(IKEY.EQ.0) THEN
	      CALL RECIP
	      RA=CELN(1)
	      RB=CELN(2)
	      RC=CELN(3)
	      RALP=CELN(4)
	      RBET=CELN(5)
	      RGAM=CELN(6)
	END IF
      READ(1,79) LAMBDA,IKEY,JKEY
79    FORMAT(F10.5,2I5)
      IF(IKEY.EQ.1)WRITE(2,80)
      IF(IKEY.EQ.0)WRITE(2,81)
      IF(JKEY.EQ.1)WRITE(2,82)
      IF(JKEY.EQ.0)WRITE(2,83)
      FAM=(2./LAMBDA)**2
      CALL SETUP(NOT,NUMTMP,AXINC,ANINC)
      DO 2    I=1,201
      READ (1,62,END=3   )   HH(I),KK(I),LL(I),TTH(I),SIG(I)
62    FORMAT(3I5,2F10.5,I5)
2     CONTINUE
3     NREL=I-1
      WRITE(6,85) NREL
      WRITE(2,85) NREL
      WRITE(2,66)
	DO 5    I=1,NREL
	      IF(IKEY.EQ.0) THEN  
C		IF INPUT IS D-SPACINGS-CONVERT TO TWO THETA
		      ST=0.5*LAMBDA/TTH(I)
		      CT = SQRT (1.0 - ST*ST)
		      TTH(I) = 180.0 - 360.0*ATAN (CT/ST)/PI
		END IF
	      WTG(I)=1.0/(SIG(I)*FAM*SIN (PI*TTH(I)/180.0))
	      IF(JKEY.EQ.1)WTG(I)=SIG(I)
		SSQ(I)=FAM*(SIN (PI*TTH(I)/360.0))**2
5	CONTINUE
C      CALCULATE REAL SPACE CONSTANTS
	CALL RECIP
	DO 6    I=1,6
		CELL(I)=CELN(I)
6	CONTINUE
	WRITE(2,65) (CELL(I),I=1,6)
7	DO 8    J=1,7
	DO 8    K=J,7
		A(J,K)=0.0
8	CONTINUE
	COSA=COS (PI*RALP/180.0)
	COSB=COS (PI*RBET/180.0)
	COSG=COS (PI*RGAM/180.0)
	SINA=SIN (PI*RALP/180.0)
	SINB=SIN (PI*RBET/180.0)
	SING=SIN (PI*RGAM/180.0)
	S11=RA*RA
	S22=RB*RB
	S33=RC*RC
	S12=2.0*RA*RB*COSG
	S13=2.0*RA*RC*COSB
	S23=2.0*RB*RC*COSA
	DO 14 I=1,NREL
		FH=HH(I)
		FK=KK(I)
		FL=LL(I)
      SSC(I)=FH*FH*S11+FK*FK*S22+FL*FL*S33+FH*FK*S12+FH*FL*S13+FK*FL*S23
      D(7)=(SSQ(I)-SSC(I))*WTG(I)
      D(1)=(2.0*FH*(FH*RA+FK*RB*COSG+FL*RC*COSB))*WTG(I)
      D(2)=(2.0*FK*(FH*RA*COSG+FK*RB+FL*RC*COSA))*WTG(I)
      D(3)=(2.0*FL*(FH*RA*COSB+FK*RB*COSA+FL*RC))*WTG(I)
      D(4)=(-PI*FK*FL*RB*RC*SINA/90.0)*WTG(I)
      D(5)=(-PI*FH*FL*RA*RC*SINB/90.0)*WTG(I)
      D(6)=(-PI*FH*FK*RA*RB*SING/90.0)*WTG(I)
C      CHECK TO SEE IF THERE ARE RELATED PARAMETERS
		DO 10 J=1,2
			DO 9     K=1,3
				IF(NUMTMP(J,K).EQ.0) GO TO 10
				IF(K.EQ.1)IT=NUMTMP(J,1)
				IF(K.EQ.1)GO TO 9    
				D(IT)=D(IT)+D(NUMTMP(J,K))
9			CONTINUE
10		CONTINUE
C      CHECK TO SEE IF THERE ARE FIXED PARAMETERS
		DO 11 J=1,6
			IF(NOT(J).EQ.0)GO TO 11
			D(J)=0.
11		CONTINUE
		DO 13 J=1,7
			DO 12 K=J,7
				A(J,K)=A(J,K)+D(J)*D(K)
12			CONTINUE
13		CONTINUE
14	CONTINUE
	DO 15 I=1,6
C	SET DIAGONAL TERM TO ZER0 FOR FIXED PARAMETERS
		IF(NOT(I).EQ.0)GO TO 15
		A(I,I)=1.
15	CONTINUE
C      ZERO TRIANGLE OF MATRIX
	DO 17 I=2,NPAR
		L=I-1
		DO 16 J=1,L
			A(I,J)=0.0
16		CONTINUE
17	CONTINUE
	A(1,1)=1.0/SQRT (A(1,1))
	DO 18 I=2,NPAR
		A(1,I)=A(1,I)*A(1,1)
18	CONTINUE
	DO 24 I=2,NPAR
		L=I-1
		DO 19 J=1,L
			A(I,I)=A(I,I)-A(J,I)*A(J,I)
19		CONTINUE
		A(I,I)=1.0/SQRT (A(I,I))
		LLL=I+1
		DO 21 J=1,L
			DO 20 K=J,L
				A(I,J)=A(I,J)-A(K,I)*A(K,J)
20			CONTINUE
			A(I,J)=A(I,J)*A(I,I)
21		CONTINUE
		IF(I.EQ.NPAR)GO TO 24
		DO 23 J=LLL,NPAR
			DO 22 K=1,L
				A(I,J)=A(I,J)-A(K,I)*A(K,J)
22			CONTINUE
			A(I,J)=A(I,J)*A(I,I)
23		CONTINUE
24	CONTINUE
	DO 25 J=1,NPAR
		A(NPAR,J)=A(NPAR,J)*A(NPAR,NPAR)
25	CONTINUE
	DO 28 I=2,NPAR
		MI=NPAR-I+1
		ATEMP=A(MI,MI)
		DO 27 J=1,MI
			L=I-1
			DO 26 K=1,L
				MK=NPAR-K+1
				A(MI,J)=A(MI,J)-A(MI,MK)*A(MK,J)
26			CONTINUE
			A(MI,J)=A(MI,J)*ATEMP
27		CONTINUE
28	CONTINUE
	DO 30 I=2,NPAR
		L=I-1
		DO 29 J=1,L
			A(J,I)=A(I,J)
29		CONTINUE
30	CONTINUE
C
C      MATRIX INVERSION FINISHED--NEXT CALCULATED SHIFTS
	DO 32 I=1,NPAR
		SHT(I)=0.0
		DO 31 J=1,NPAR
			SHT(I)=SHT(I)+A(I,J)*A(J,MM)
31		CONTINUE
32	CONTINUE
C      APPLY EQUAL SHIFTS TO RELATED PARAMETERS
	DO 34 I=1,2
		DO 33 K=1,3
			IF(NUMTMP(I,K).EQ.0)GO TO 34
			IF(K.EQ.1)IT=NUMTMP(I,1)
			IF(K.EQ.1)GO TO 33
			SHT(NUMTMP(I,K))=SHT(IT)
33		CONTINUE
34	CONTINUE
C315   WRITE ( 2,64) ((A(I,J),J=1,7),I=1,7)
C       WRITE (2,64) (SHT(I), I=1,6)
C      ADD SHIFTS TO RECIP.PARAMETERS THAN CONVERT TO REAL SPACE
	RA=RA+SHT(1)
	RB=RB+SHT(2)
	RC=RC+SHT(3)
	RALP=RALP+SHT(4)
	RBET=RBET+SHT(5)
	RGAM=RGAM+SHT(6)
	CALL RECIP
C      CALCULATE REAL SPACE SHIFTS
	DO 35 J=1,6
		SHT(J)=CELN(J)-CELL(J)
		CELL(J)=CELN(J)
35	CONTINUE
       WRITE (2,65) (CELL(J),J=1,6)
C      CHECK SIZE OF SHIFT.  DO ANOTHER CYCLE?
	DO 36 J=1,3
		IF(ABS(SHT(J)).GT.AXINC)GO TO 7   
36	CONTINUE
	DO 37 J=4,6
		IF(ABS(SHT(J)).GT.ANINC)GO TO 7   
37    CONTINUE
	ICOUNT=0
	DO 38 I=1,6
		IF(NOT(I).EQ.1)ICOUNT=ICOUNT+1
38	CONTINUE
	FNMN=NREL-7+ICOUNT
	WRITE(6,66)
	WRITE(6,65) (CELL(J),J=1,6),A(7,7)
C
C      CALCULATE STANDARD DEVIATIONS.
	RA=1001.0*RA/1000.0
	CALL RECIP
	RA=RA*1000.0/1001.0
	DO 39 I=1,6
		B(1,I)=(CELN(I)-CELL(I))*1000.0/RA
39	CONTINUE
	RB=RB*1001.0/1000.0
	CALL RECIP
	RB=RB*1000.0/1001.0
	DO 40 I=1,6
		B(2,I)=(CELN(I)-CELL(I))*1000.0/RB
40	CONTINUE
	RC=RC*1001.0/1000.0
	CALL RECIP
	RC=RC*1000.0/1001.0
	DO 41 I=1,6
		B(3,I)=(CELN(I)-CELL(I))*1000.0/RC
41	CONTINUE
	RALP=RALP+0.01
	CALL RECIP
	RALP=RALP-0.01
	DO 42 I=1,6
		B(4,I)=(CELN(I)-CELL(I))*100.0
42	CONTINUE
	RBET=RBET+0.01
	CALL RECIP
	RBET=RBET-0.01
	DO 43 I=1,6
		B(5,I)=(CELN(I)-CELL(I))*100.0
43	CONTINUE
	RGAM=RGAM+0.01
	CALL RECIP
	RGAM=RGAM-0.01
	DO 44 I=1,6
		B(6,I)=(CELN(I)-CELL(I))*100.0
44	CONTINUE
	DO 47 I=1,2
		DO 46 K=1,3
			IF(NUMTMP(I,K).EQ.0)GO TO 47
			IF(K.EQ.1)IT=NUMTMP(I,1)
			IF(K.EQ.1)GO TO 46
			DO 45 J=1,6
				B(IT,J)=B(IT,J)+B(NUMTMP(I,K),J)
45			CONTINUE
46		CONTINUE
47	CONTINUE
	DO 49 I=1,6
		IF(NOT(I).EQ.0)GO TO 49
		DO 48 J=1,6
			B(I,J)=0.0
48		CONTINUE
49	CONTINUE
50	FNMN=A(7,7)/FNMN
	DO 53 I=1,6
		SHT(I)=0.0
		DO 52 J=1,6
			TEM(J)=0.0
			DO 51 K=1,6
				TEM(J)=TEM(J)+A(J,K)*B(K,I)
51			CONTINUE
			SHT(I)=SHT(I)+TEM(J)*B(J,I)
52		CONTINUE
		SHT(I)=SQRT (SHT(I)*FNMN)
53	CONTINUE
	DO 55 I=1,2
		DO 54 K=1,3
			IF(NUMTMP(I,K).EQ.0)GO TO 55
			IF(K.EQ.1)IT=NUMTMP(I,1)
			IF(K.EQ.1)GO TO 54
			SHT(NUMTMP(I,K))=SHT(IT)
54		CONTINUE
55	CONTINUE
	WRITE (2,68) (SHT(I),I=1,6)
	WRITE(6,68) (SHT(I),I=1,6)
	WRITE (2,71)
	COSA=COS (PI*CELL(4)/180.0)
	COSB=COS (PI*CELL(5)/180.0)
	COSG=COS (PI*CELL(6)/180.0)
	S=SQRT(1.0-COSA*COSA-COSB*COSB-COSG*COSG+2.*COSA*COSB*COSG)
	VOL=CELL(1)*CELL(2)*CELL(3)*S
	WRITE (2,65) RA,RB,RC,RALP,RBET,RGAM,VOL
	WRITE (2,69)
	WRITE (2,70)
	DO 56 I=1,NREL
		TTC=SSC(I)/FAM
		COSQ=    1.0-TTC
		COT=SQRT (COSQ/TTC)
		TTC=180.0-360.0*ATAN (COT)/PI
		DEL=TTH(I)-TTC
	WRITE (2,63)       HH(I),KK(I),LL(I),SSQ(I),SSC(I)
	1      ,TTH(I),TTC,DEL,WTG(I)
56	CONTINUE
C	THIS PATCH IS TO GIVE A DOBS-DCALC TABLE   KHG 7/29/80
	WRITE(2,156)
156	FORMAT(//'    H    K    L   D OBS     D CALC'/)
	DO 157 I=1,NREL
		DD=SQRT(1/SSQ(I))
		DDD=SQRT(1/SSC(I))
		WRITE(2,158) HH(I),KK(I),LL(I),DD,DDD
157	CONTINUE
158	FORMAT(1X,3I5,2F10.4)
	WRITE(2,86)
	WRITE(2,66)
	WRITE(2,65) (CELL(J),J=1,6),A(7,7)
	WRITE(2,68) (SHT(I),I=1,6)
	DO 58 I=1,3
		IF(CELL(I+3).EQ.90.)GO TO 57
		COSINE(I)=COS(PI*CELL(I+3)/180.)
		SINE(I)=SIN(PI*CELL(I+3)/180.)
		GO TO 58
57		COSINE(I)=0.
		SINE(I)=1.
58	CONTINUE
	DO 59 I=4,6
		CELERR(I)=SINE(I-3)*CELERR(I)*PI/180.
59	CONTINUE
	DVDA=CELL(2)*CELL(3)*CELERR(1)
	DVDB=CELL(1)*CELL(3)*CELERR(2)
	DVDC=CELL(1)*CELL(2)*CELERR(3)
	DVDAL=CELERR(4)*(COSINE(2)*COSINE(3)-COSINE(1))
	DVDBE=CELERR(5)*(COSINE(1)*COSINE(3)-COSINE(2))
	DVDGA=CELERR(6)*(COSINE(1)*COSINE(2)-COSINE(3))
	VERROR=S*(DVDA+DVDB+DVDC)+CELL(1)*CELL(2)*CELL(3)/S*(
	1      DVDAL+DVDBE+DVDGA)
	WRITE(6,87) VOL,VERROR
	WRITE(2,87) VOL,VERROR
	WRITE(6,88)
	CALL YESNO(IKEY)
	IF(IKEY.EQ.1)CALL DENSTY(VOL,VERROR)
	CALL EXIT
63    FORMAT(3I5,2F10.5,3F9.4,F9.2)
64    FORMAT(1H 7E12.5)
65    FORMAT(1H0,6F9.4,E12.5//)
66    FORMAT(/1H0,7X,'A',8X,'B',8X,'C',6X,'ALPHA',4X,'BETA',
     1     5X,'GAMMA',2X,'SUM WT SQ DEV')
67    FORMAT(9H SHIFTS =/1H 6F9.4,E12.5)
68    FORMAT(22H STANDARD DEVIATIONS =/1H 6F9.4//)
69    FORMAT(1H ,3X,1HH,4X,1HK,4X,1HL,2X,8H(1/D)**2,2X,
     1 8H(1/D)**2,2X,7H2 THETA,2X,7H2 THETA,2X,7H2 THETA,
     2     3X,6HWEIGHT)
70    FORMAT(20X,3HOBS,7X,4HCALC,6X,3HOBS,5X,4HCALC,5X,4HDIFF)
71    FORMAT(1H0,3X,'A*',7X,'B*',7X,'C*',5X,'ALPHA*',3X,
     1      'BETA*',4X,'GAMMA*',2X,'VOLUME')
80    FORMAT(' THE REFLECTIONS ARE GIVEN AS TWO THETA VALUES')
81    FORMAT(' THE REFLECTIONS ARE GIVEN AS D-SPACINGS')
82    FORMAT(' THE REFLECTIONS HAVE BEEN GIVEN  WEIGHTS')
83    FORMAT(' THE REFLECTIONS HAVE BEEN GIVEN ESTIMATED ERRORS')
84    FORMAT(' LAMBDA = ',F10.5/)
85    FORMAT(/1X,I5,' REFLECTIONS')
86    FORMAT(////72('*')//)
87    FORMAT(' CELL VOLUME = ',F10.3/
     1' ERROR IN CELL VOLUME = ',F10.3/)
88    FORMAT(' DO YOU WANT TO CALCULATE DENSITY?  '$)
      END
C
      SUBROUTINE RECIP
      DIMENSION CELN(6)
      COMMON RA,RB,RC,RALP,RBET,RGAM,CELN
      PI=3.1415926536
      SINA=SIN (PI*RALP/180.0)
      SINB=SIN (PI*RBET/180.0)
      SING=SIN (PI*RGAM/180.0)
      COSA=COS (PI*RALP/180.0)
      COSB=COS (PI*RBET/180.0)
      COSG=COS (PI*RGAM/180.0)
      VOL=1.0-COSA*COSA-COSB*COSB-COSG*COSG+2.0*COSA*COSB*COSG
      VOL=RA*RB*RC*SQRT (VOL)
      CELN(1)=RB*RC*SINA/VOL
      CELN(2)=RA*RC*SINB/VOL
      CELN(3)=RA*RB*SING/VOL
      CT=(COSB*COSG-COSA)/(SINB*SING)
      ST=SQRT (1.0-CT*CT)
      CELN(4)=90.0-180.0*ATAN (CT/ST)/PI
      CT=(COSA*COSG-COSB)/(SINA*SING)
      ST=SQRT (1.0-CT*CT)
      CELN(5)=90.0-180.0*ATAN (CT/ST)/PI
      CT=(COSA*COSB-COSG)/(SINA*SINB)
      ST=SQRT (1.0-CT*CT)
      CELN(6)=90.0-180.0*ATAN (CT/ST)/PI
      RETURN
      END
C
      SUBROUTINE SETUP(NOT,NUMTMP,AXINC,ANINC)
      DIMENSION NOT(6),NUM(6),DA(6),DUM(6),HOLD(6),NUMTMP(2,3)
      DATA DA/'A','B','C','AL','BE','GA'/
      DATA BLANK/' '/, BLNK/'  '/
1     DO 2    I=1,6
      NOT(I)=0
2     NUM(I)=0
      WRITE(6,17)
      CALL YESNO(IKEY)
      IF(IKEY.EQ.0)GO TO 4   
      WRITE(6,18)
      READ(5,19) (HOLD(I),I=1,3)
      WRITE(6,20)
      READ(5,21) (HOLD(I),I=4,6)
      DO 3    I=1,6
      DO 3    J=1,6
3     IF(HOLD(I).EQ.DA(J))NOT(J)=1
4     WRITE(6,22)
      CALL YESNO(IKEY)
      IF(IKEY.EQ.0)GO TO 6   
      WRITE(6,23)
      READ(5,19) (HOLD(I),I=1,3)
      WRITE(6,24)
      READ(5,21) (HOLD(I),I=4,6)
      DO 5    I=1,6
      DO 5    J=1,6
5     IF(HOLD(I).EQ.DA(J))NUM(J)=1
6     DO 10 K=1,2
      MIN=1+(K-1)*3
      MAX=MIN+2
      DO 7    I=MIN,MAX
      IF(NUM(I).EQ.1)GO TO 8   
7     CONTINUE
      GO TO 10
8     NEWMIN=I+1
      DO 9    J=NEWMIN,MAX
      IF(NUM(J).EQ.1)NOT(J)=1
9     CONTINUE
10    CONTINUE
      AXINC=.0005
      ANINC=.005
C     WRITE(6,25)
C     READ(5,26) AXIS,ANGLE
C     IF(AXIS.NE.0)AXINC=AXIS
C     IF(ANGLE.NE.0.)ANINC=ANGLE
      WRITE(6,27)
      CALL YESNO(IKY)
      IF(IKY.NE.1)GO TO 1   
      DO 11 I=1,3
      DUM(I)=BLANK
11    IF(NUM(I).NE.0)DUM(I)=DA(I)
      WRITE(2,28) (DUM(I),I=1,3)
      DO 12 I=4,6
      DUM(I)=BLNK
12    IF(NUM(I).NE.0)DUM(I)=DA(I)
      WRITE(2,29) (DUM(I),I=4,6)
      DO 13 I=1,3
      DUM(I)=BLANK
13    IF(NOT(I).NE.0)DUM(I)=DA(I)
      DO 14 I=4,6
      DUM(I)=BLNK
14    IF(NOT(I).NE.0)DUM(I)=DA(I)
      WRITE(2,30) DUM
      DO 16 I=1,2
      M=1
      DO 15 J=1,3
      NUMTMP(I,J)=0
      L=J+(I-1)*3
      IF(NUM(L).NE.1)GO TO 15
      NUMTMP(I,M)=L
      M=M+1
15    CONTINUE
16    CONTINUE
      RETURN
17    FORMAT(' ARE THERE PARAMETER(S) TO BE HELD AT FIXED VALUE(S) ? '$)
18    FORMAT(' LENGTH(S) (A,B,C) '$)
19    FORMAT(3(1A1,1X))
20    FORMAT(' ANGLE(S) (AL,BE,GA) '$)
21    FORMAT(3(1A2,1X))
22    FORMAT(' ARE THERE PARAMETERS TO BE RELATED? '$)
23    FORMAT(' RELATED LENGTH PARAMETERS (A,B,C) '$)
24    FORMAT(' RELATED ANGLE PARAMETERS (AL,BE,GA) '$)
25    FORMAT(' ENTER THE SIZE OF AXIS AND ANGLE SHIFTS'/
     1' TO WHICH THE REFINEMENT IS TO BE CARRIED.  TO USE'/
     2' THE DEFAULT VALUES (0.0005A & 0.005DEG) HIT CR')
26    FORMAT(2F) 
27    FORMAT(' OK ? '$)
28    FORMAT(' RELATED LENGTHS - ',3(1A1,1X))
29    FORMAT(' RELATED ANGLES - ',3(1A2,1X))
30    FORMAT(' PARAMETERS TO BE HELD AT FIXED VALUES'/
     1     3(1A1,2X),5X,3(1A2,2X))
      END
C
C
C
      SUBROUTINE DENSTY(VOL,DLTVOL)
      INTEGER FORMLA(60),BLANK,LABEL(2,10),FUNC(8),FUNCT
      DIMENSION TYPE(10),NTYPE(10),ITYPE(10)
      DIMENSION ELEM1(49),ELEM2(49),ELEMNT(98)
      DIMENSION WT1(49),WT2(49),WT(98)
      EQUIVALENCE (ELEM1,ELEMNT(1)),(ELEM2,ELEMNT(50))
      EQUIVALENCE (WT1,WT(1)),(WT2,WT(50))
      DATA BLANK,BLANK2/1H ,2H  /
      data iodine/' j'/
      DATA ELEM1/'H ','HE','LI','BE','B ','C ','N ','O ','F ',
     1  'NE','NA','MG','AL','SI','P ','S ','CL','AR','K ',
     2  'CA','SC','TI','V ','CR','MN','FE','CO','NI','CU',
     3  'ZN','GA','GE','AS','SE','BR','KR','RB','SR','Y ',
     4  'ZR','NB','MO','TC','RU','RH','PD','AG','CD','IN'/
      DATA ELEM2/
     5  'SN','SB','TE','I ','XE','CS','BA','LA','CE','PR',
     6  'ND','PM','SM','EU','GD','TB','DY','HO','ER','TM',
     7  'YB','LU','HF','TA','W ','RE','OS','IR','PT','AU',
     8  'HG','TL','PB','BI','PO','AT','RN','FR','RA','AC',
     9  'TH','PA','U ','NP','PU','AM','CM','BK','CF'/
      DATA WT1/1.00797,4.0026,6.939,9.0122,10.811,12.01115,
     114.0067,15.9994,18.9984,20.183,22.9898,24.312,26.9815,
     228.086,30.9738,32.064,35.453,39.948,39.102,40.08,44.96,
     347.90,50.942,51.996,54.938,55.847,58.933,58.71,63.54,
     465.37,69.72,72.59,74.922,78.96,79.909,83.80,85.47,87.62
     5,88.905,91.22,92.906,95.94,98.,101.07,102.905,106.4,107
     6.870,112.40,114.82/
      DATA WT2/ 118.69,121.75,127.60,126.904,131.30,
     7132.905,137.34,138.91,140.12,140.907,
     8144.24,147.,150.35,151.96,157.25,158.924,162.50,164.930
     9,167.26,168.934,173.04,174.97,178.49,180.948,183.85,186
     1.2,190.2,192.2,195.09,196.967,200.59,204.37,207.19,208.
     1980,210.,210.,222.,223.,226.,227.,232.038,231.,238.03,
     2237.,242.,243.,247.,247.,249./
1     I=0
      WRITE(6,9   )
      WRITE(6,10  )
2     II=I+1
      WRITE(6,11) II
      READ(5,12) ATNAME,NOAT
      IF(ATNAME.EQ.BLANK2)GO TO 3   
      I=I+1
      NTYPE(I)=NOAT
      TYPE(I)=ATNAME
      GO TO 2   
3     NTBL=I
      WRITE(6,13)( (TYPE(I),NTYPE(I)),I=1,NTBL)
      WRITE(6,14)
      CALL YESNO(IKEY)
      IF(IKEY.NE.1)GO TO 1   
      WRITE(2,13) ((TYPE(I),NTYPE(I)),I=1,NTBL)
C
C      CHECK TO SEE IF ELEMENT SYMBOLS ARE CORRECT
C
C
      DO 6    I=1,NTBL
	      DO 4    J=1,98
		      IF(TYPE(I).EQ.ELEMNT(J))GO TO 5   
4	     CONTINUE
	      WRITE(6,20) TYPE(I)
	      GO TO 1   
5	     ITYPE(I)=J
6     CONTINUE
	      WEIGHT=0
      DO 7    I=1,NTBL
	      IAN=(ITYPE(I))
	      WEIGHT=WEIGHT+NTYPE(I)*WT(IAN)
7     CONTINUE
      WRITE(6,21) WEIGHT
      WRITE(2,21) WEIGHT
      WRITE(6,22)
      READ(5,*) Z
      IF(Z.LE.1)Z=1
      WRITE(2,15) Z
      DENST=1.6604*WEIGHT*Z/VOL
      DLTDEN=1.6604*WEIGHT*Z/(VOL-DLTVOL) - DENST
      WRITE(6,24) DENST,DLTDEN
      WRITE(2,24) DENST,DLTDEN
      RETURN
9     FORMAT('  INPUT MOLECULAR COMPOSITION '//)
10    FORMAT('  ATOM TYPE,NUMBER ')
11    FORMAT(I5,1X$)
12    FORMAT(A2,I)
13    FORMAT('  EMPIRICAL FORMULA '//10(4X,A2,I2))
14    FORMAT('  OK?   '$)
15    FORMAT('  NO. OF MOLECULES IN UNIT CELL - ',F5.2)
16    FORMAT(' ENTER EMPERICAL FORMULA. SEPARATE ELEMENTS BY
     1COMMA.'/' E.G. MN2,NI1,O10,H24,C46,N6'/)
17    FORMAT(60A1)
18    FORMAT(' THERE ARE ',I2,' ATOM TYPES'/
     1' THE EMPERICAL FORMULA IS ',10(A2,I3,3X))
19    FORMAT(' OK?   '$)
20    FORMAT(' IS ',A2,' A NEW ELEMENT?  I DO NOT RECOGNIZE IT.'/)
21    FORMAT(' THE MOLECULAR WEIGHT IS ',F10.4/)
22    FORMAT(' NO. OF MOLECULES PER UNIT CELL = '$)
24    FORMAT(' CALCULATED DENSITY = ',F10.3/
     1     8X,'EST. ERROR = ',F10.3)
      END
C
	SUBROUTINE YESNO(IKEY)
C	IF SYM IS 'Y' THAN IKEY IS 1 OTHERWISE IKEY=0
	CHARACTER *1 SYM
	IKEY=0
	READ(5,1) SYM
1	FORMAT(A)
	IF (SYM.EQ.'Y'.OR.SYM.EQ.'y') iKEY=1
	RETURN
	END
