CDC    PROGRAM DATAPP(INPUT,OUTPUT,TAPE3,TAPE4,TAPE5=INPUT,PUNCH,
CDC  1TAPE6=OUTPUT,TAPE1=TAPE3,TAPE2,TAPE7)
C4Z DATAPH-DATA PROCESSOR FOR INTEGRATED INTENSITIES(NEUTRON AND X-RAY)
C  FOR LORENTZ, POLARIZATION, EXTINCTION, AND ABSORPTION. SUITABLE FOR
C  DATA OBTAINED AUTOMATICALLY BY DIGITAL STEP SCAN AS WELL AS MANUALLY
C  CDC6600/FORTRAN66///COPPENS,P. , MODIFIED BY W.C.HAMILTON
C  DEPT. OF CHEM. BROOKHAVEN NATIONAL LABORATORY, UPTON,NEW YORK,11973
C  PROCESSES INTEGRATED INTENSITY AND BACKROUNDS. FINDS THE OBSERVED
C  INTENSITY AND CALCULATES THE STANDARD DEVIATION AND WEIGHT TO BE
C  APPLIED. ON OPTION IT WILL APPLY CORRECTIONS FOR, LORENTZ FACTOR,
C  POLARIZATION, ABSORPTION, AND EXTINCTION. INPUT AND OUTPUT MAY BE ON
C  TAPE OR CARDS. CP MEMORY REQUIRED IS APPX. 42000 OCTAL.
C   EXECUTION FIELD LENGTH IS 37000
C     VERSION MARCH 68  USES OMEGA AS ARBITRARY ANGLE
C     DATA PROCESSING PROGRAM VERSION DATAPP FOR PICKER DIFFRACTOMETER
C     CRYSTAL ORIENTATION SPECIFIED BY MATRIX RELATING RECIPROCAL
C     CRYSTAL SYSTEM TO GONIOMETERHEAD SYSTEM
C     INCLUDES OPTIONS FOR ABSORPTION, EXTINCTION AND SPOTSHAPE CORRECTI
	REAL*8 TEMP
	DIMENSION TITLE(20)
	DIMENSION IHSTD(3,30)
	DIMENSION NSERIAL(10)
	DIMENSION ICN(16),XNSO(3,10)
	DIMENSION OBSFI(10),BACFK(10)
	DIMENSION SLO(10),SLB(10),XN(3,10),WI(3,3)
	DIMENSION RCP(6),RECIP(3,3),DIRECT(3,3),DRC(6),SHUFEL(3,3),OM(10)
	DIMENSION X(3,10),OBSI(10),BACK(10),FI(10),EPOA(3000)
	DIMENSION NGAUSS(3),EPOB(1000),DATA(40,4),SCB(10)
	DIMENSION GRID(4,5000),TWOTETA(2)
	DIMENSION NSCOFF(10)
	DIMENSION SUM(40),TT(3),E(3),W(3,3)
	DIMENSION XNOUT(3),XSOUT(3)
	EQUIVALENCE (ICN(1),IABS),(ICN(2),IEXT)
	EQUIVALENCE(ABS,COR)
	EQUIVALENCE(WI,SHUFEL)
	INTEGER H,H1,H2
	DATA NFLAG,NSCALE,JM/1,1,0/
	DATA FACSIG,PI/0.,3.1415926563/
	DATA MDGRID,NDGRID/4,5000/
	DATA BLANKS/4H    /
	 DATA NWALT/0/
	DATA XN/30*0./
	DATA XNSO/30*0./
	DATA RAD/57.2957794636/
	COMMON RCP,RECIP,DIRECT,DRC,SHUFEL,X,OBSI,BACK,FI,EPOA,EPOB
	COMMON ELAMDA,SINSQT,COST,SINT,EMU,NPOINTS,NR
	COMMON     EXT,ABS,DATA,NO,NFACES
	COMMON NGAUSS,VOL,TBAR
	COMMONSUM,W,TT,E  ,GRID
	COMMON /DATDIR/ ICN
	COMMON /METRICS/ AA(3,3),BB(3,3)
	COMMON /ZEROES/ PHIZRO,OMZRO,TWOTHZO,CHIZRO
	COMMON /CN/ IDIST
	COMMON/CRYO/THICK(2),TBO,TB,EMUBER
	DATA EMUBER,THICK(1),THICK(2),ABSBE/.453,.0375,.050,1.0/
31	CONTINUE
	TBAR=0.
	ABS=1.
	JDENT=0
	IDENT=0
	SCLA=1.
	READ 500, TITLE
	PRINT 501, TITLE
	PHIZRO=OMZRO=TWOTHZO=CHIZRO=0.
	LI = 5
	LO = 6
	 CALL RDCL(DRC,RCP,V,ELAMDA,LI,LO)
C
C          PATCH HERE TO USE *RDCL* AND AVOID REPROGRAMMING ENTIRE
C          PROGRAM
C
	 DO 1 I=1,3
	  DO 1 J=1,3
	   DIRECT(I,J)=AA(I,J)
1	   RECIP(I,J)=BB(I,J)
	READ 513,((SHUFEL(I,J),J=1,3),I=1,3)
	 IF (SHUFEL(1,1).EQ.0..AND.SHUFEL(1,2).EQ.0..AND.SHUFEL(1,3).EQ.
	*     0.) CALL UMATRIX(SHUFEL,ELAMDA)
	PRINT5130,((SHUFEL(I,J),J=1,3),I=1,3)
5130	FORMAT(' THE WN MATRIX (BEDNOWITZ-COPPENS DEFINITION) IS'/
	*     3(3F10.6/))
	DO9K=1,3
	DO9L=1,3
9	W(K,L)=SHUFEL(K,L)*RCP(L)
	PRINT 5131, W
5131	FORMAT(' THE U MATRIX (HAMILTON DEFINITION, REFERRED TO NON-NORMAL
	'IZED AXES, = W TRANSPOSE) IS '/3(3F10.6/))
	PRINT 5132
5132	FORMAT(' THE MATRIX USED IN DATAPP FOR ALL CALCULATIONS IS W, WHER
	'E W = U TRANSPOSE (ABOVE)'/)
	CALL AMI3(W,WI,DET)
	READ 503, (ICN(I),I=1,16)
	 IGRID=LOCF(GRID(1,1))
	 DO 10 I=1,16
10	   ICN(I)=ICN(I)+0
	PRINT 504, (ICN(I),I=1,16)
	IDIST=ICN(4)
	 IF(ICN(7) .GE. 0) REWIND 4
	 IF(ICN(8) .NE. 0) REWIND 3
	IF(ICN(7).LE.0)GOTO4
C     IF DATA ARE TO BE CORRECTED FOR ABSORPTION SET UP GAUSSIAN GRID IN
4	IF(ICN(9).LE.0)GOTO40
	PRINT521
	ND=ICN(9)
	DO43I=1,ND
	READ520,H1,H2,H3,OA
520	FORMAT(3F3.0,F9.3) 8
	CALLDSTAR(H1,H2,H3, RCP,DST)
	FA=0.
	CALLANGLES(OA,FA,W,DST,RCP,PHI,CHI,H1,H2,H3)
	DST=DST/2.
	 SINT=DST*ELAMDA
	 TEMP=2.*RAD*ASIN(SINT)
	 ENCODE(8,525,TWOTETA)  TEMP
525	 FORMAT(F8.4)
	 IF (ELAMDA.EQ.0.) TWOTETA(1)=BLANKS;TWOTETA(2) = BLANKS
	PRINT 522, H1,H2,H3,PHI,OA,TWOTETA,DST,CHI
521	FORMAT('-ANGLES CHECK (ALL ANGLES ARE TRUE VALUES, ZEROES MUST BE
	'ADDED)'/'    H   K   L    PHI(T)     OM(T)    TWO THETA(T)  (SINT/
	'L)   CHI(T)'/)
522	FORMAT(1H ,3F4.0,2F10.4,7X,2A4,2F10.4)
43	CONTINUE
40	IF(IABS.EQ.0.AND.IEXT.EQ.0)GOTO44
	CALLGRIT2(GRID,MDGRID,NDGRID)
	Y=0.0
C     SET UP EXPONENTIAL TABLE IF DATA ARE TO BE CORRECTED FOR ABSORPTIO
	DO 41 I=1,3000
	 EPOA(I)=EXP(-Y)
41	Y=Y+0.1
	Y=0.0
	DO 42 I=1,1000
	 EPOB(I)=EXP(-Y)
42	Y=Y+0.0001
	GOTO5
C     MAKE OUTPUT BLANK IF NO COMPONENTS ARE TO BE CALCULATED
44	DO 45I=1,3
	XNOUT(I)=0.0
45	XSOUT(I)=0.0
	TOUT=0.0
5	READ 506, IQ,ELAM,CONE,CTWO,EMU,SC,OM0,FI0,T2M
	 OM0=OM0+0.  ;  FI0=FI0+0.  ;  T2M=T2M+0.  ;  CTWO=CTWO+0.
	 CONE=CONE+0. ;  SC=SC+0.
	 ILIST = 0
	 IF (ELAMDA.EQ.0.) ELAMDA=ELAM
	 IF (PHIZRO.EQ.0.)GOTO8
	  OM0=OMZRO
	  FI0=PHIZRO
8	PRINT 516, ELAMDA,CONE,CTWO,SC,IQ,T2M
	ISHIELD=ICN(12)
516	FORMAT('0LAMDA ='F10.4'   CONE ='F10.4'   CTWO ='F10.4'     SCALE
	'FACTOR = 'F10.4'   NUMBER OF SCALE FACTOR'I5/'  2 THETA MONOCHROMA
	'TOR'F10.2)
	IF(IABS.GT.0.OR.IEXT.GT.0) PRINT 517,EMU,OM0,FI0
517	FORMAT('0MU ='F10.4'   OM ZERO ='F10.4'  FI ZERO = 'F10.4)
	 C2TM2=(COS(T2M/RAD))**2
	IF (ICN(3) .EQ. 0) GO TO 540
	READ 530, START, IDEGREE, NSTD
	READ 531, ((IHSTD(I,J), I = 1, 3), J = 1, NSTD)
	IF (IDEGREE .EQ. 0) PRINT 533
	IF (IDEGREE .EQ. 1) PRINT 534
	PRINT 532, START, NSTD, ((IHSTD(I,J), I = 1, 3), J = 1, NSTD)
530	FORMAT (F10.5, 2I10)
531	FORMAT (4(5X, 3I5))
533	FORMAT ('0LINEAR, TWO-POINT INTERPOLATION FOR SCALING ACCORDING TO
	1 THE VARIATION OF THE STANDARD REFERENCE REFLECTIONS')
534	FORMAT ('0CUBIC POLYNOMIAL, FOUR-POINT INTERPOLATION FOR SCALING A
	1CCORDING TO THE VARIATION OF THE STANDARD REFERENCE REFLECTIONS')
532	FORMAT ('0START = ' F10.5, ' STARTING VALUE OF THE FACTOR FOR SCAL
	1ING ACCORDING TO THE VARIATION OF THE STANDARD REFERENCE REFLECTIO
	2S'/ '0NSTD  = ' I10, ' KINDS OF STANDARD REFERENCE REFLECTIONS'/
	3'0   H    K    L'// (3I5/))
	ICOUNT = 0
540	PRINT 5081
	NLINES = 0
51	NR=0
6	J=NR+1
	IF(ICN(3).EQ.0)     GO TO 100
200	IF(NFLAG.EQ.0) GO TO 202
300	CONTINUE
650	FORMAT(6I3)
100	NFLAG=0
202	CONTINUE
	CALL IMPD (X(1,J), X(2,J), X(3,J), OBSI(J), BACK(J), SCB(J),
	1OM(J), FI(J), SLO(J), SLB(J), OBSFI(J), BACFK(J), IDENT,
	2NSCOFF(J), ICN, NBACK, NSERIAL(J))
	ICOUNT = ICOUNT + 1
	H = X(1,J)
	IF(H.GE.99000) X(1,J)=X(1,J)-100000.
	K = X(2,J)
	L = X(3,J)
C    ICN(11)=1 FOR LAST DIGIT INCLUDED
	IF(ICN(11).NE.0) GO TO 2001
	 OBSI(J)=OBSI(J)+.45
	 BACK(J)=BACK(J)+(NBACK*0.45)
	 OBSFI(J)=OBSFI(J)+0.45
	 BACFK(J)=BACFK(J)+(NBACK*0.45)
C
C          FACTOR OF .45 ACCOUNTS STATISTICALLY FOR LOST INTEGER IN
C          FACTOR OF NBACK*.45 BECAUSE BACFK AND BACK ARE SUM OF NBACK
C          BACKGROUNDS,EACH IN ERROR OF .45
C          LOWEST DECADE OF PICKER INSTRUMENT
C
2001	JM=JM+1
C
C     ICN(11) = 2 FOR STEP SCAN DATA FROM PROGRAM PROFILE
C
	IF(ICN(11).NE.2) GOTO4000
	ISIG=OBSI(J)-BACK(J)
	M=BACK(J)
	GOTO4001
4000	CONTINUE
	ISIG  =    OBSI(J) - BACK(J)*SLO(J)/SLB(J)
	M =(SLO(J)/SLB(J))*BACK(J)
4001	CONTINUE
	IF(ICN(3)-1) 601,602,603
603	IF (H .LT. 99000) GO TO 601
	IDIR = 1
	CALL SCALEB (H, K, L, ISIG, M, NSTD, IHSTD, IDIR, ICOUNT, FACTOR,
	1 ILIST)
	ILIST = 1
	GO TO 601
602	CALL SCALED(H,K,L,ISIG,M,NSCALE,NREJ,JM,SCLA,NLINES)
	SCB(J) = SCB(J)*SCLA
	IF(NREJ.GT.0) GO TO 202
601	IF(IDENT-1)7,70,70
7	NR=NR+1
	IF(10-NR) 82,82,6
70	IF(NR-1)270,82,82
82	CONTINUE
	DO27NO=1,NR
C     SCALE FACTORS MADE EQUAL TO ONE IF NOT SPECIFIED ON CARDS
	IF(SCB(NO))831,831,832
831	SCB(NO)=1.
832	SCA=SCB(NO)*SC
C      SUBTRACT OMEGA ZERO AND FI ZERO
	OM(NO)=OM(NO)-OM0
	FI(NO)=FI(NO)-FI0
	CALL DSTAR(X(1,NO),X(2,NO),X(3,NO),RCP,DST)
	SINOL=DST/2.
	SINT=ELAMDA*SINOL
	SINSQT=SINT**2
	 COST=SQRT(1.-SINSQT)
C     CALL ABSORPTION SUBROUTINE FOR ABSORPTION OR EXTINCTION CALCULATIO
	IF(IABS.EQ.0.AND.IEXT.EQ.0) GOTO19
	CALLANGLES(OM(NO),FI(NO),W,DST,RCP,PHI,CHI,X(1,NO),X(2,NO),X(3,NO)
	1)
	ABS=0.
	CALL DIRCOS(GRID,MDGRID,NDGRID,APHI,XN,XNSO)
	IF(ICN(12).NE.0) ABSBE= EXP(-(TB+TBO)*EMUBER)
	ABS=ABS*ABSBE
	TOUT = TBAR
	DO 523 I = 1,3
	XNOUT(I) = XN(I,NO)
523	XSOUT(I) = XNSO(I,NO)
19	CONTINUE
	 TETA=180./PI*ASIN(SINT)
C     CALCULATE LP FACTOR
21	ELP = 4.*COST*SINT
C     FOR X-RAY DATA APPLY THE POLARIZATION FACTOR
	IF(ICN(5).EQ.0)  ELP=ELP/(C2TM2+(1.-2.*SINSQT)**2)
C
C     ICN(11) = 2 FOR STEP SCAN DATA FROM PROGRAM PROFILE
C
	IF(ICN(11).EQ.2) GOTO240
	IF(SLO(NO).EQ.0) GO TO 232
231	SACK = BACK(NO)*SLO(NO)/SLB(NO)
	TACK =     SACK*SLO(NO)/SLB(NO)
	GO TO 233
232	SACK = BACK(NO)
	TACK = SACK
C     CALCULATE FOBS AND WEIGHT
233	EMINB=OBSI(NO)-SACK
C*******
	IF(ICN(11).EQ.0) GO TO 2002
	SIGSQ=(OBSI(NO)+TACK)+(EMINB*CONE)**2
	GO TO 2003
C
C     ENTER HERE FOR STEP SCAN DATA FROM PROGRAM PROFILE
C
240	SACK=BACK(NO)
	EMINB=OBSI(NO)-SACK
	SIGSQ=SLO(NO)**2
	SIGSQ=SIGSQ+(EMINB*CONE)**2
	GOTO2003
2002	CONTINUE
C*******
C     FACTOR 10 HERE BECAUSE LAST DECADE OF SCALER NOT ON TAPE
	SIGSQ= .1*(OBSFI(NO)+.825+(BACFK(NO)+NBACK*0.825)*(SLO(NO)/SLB(NO)
	1)**2 )  +  (EMINB*CONE)**2
C
C          ABOVE FACTOR OF .825 ACCOUNTS FOR INCREASE IN SIGMA FOR
C          F**2 WHEN LOWEST DECADE LOST, AS IN PICKER INSTRUMENT.
C          FACTOR OF NBACK*0.825 IS BECAUSE BACK AND BACFK ARE SUMS OF
C          NBACK BACKGROUNDS, DEFINED IN IMP ROUTINE
C          SEE WRITE UP
C
2003	SIGMASQ=SIGSQ+(CTWO*EMINB*(1./COR-1.))**2
	 SIGMA=SQRT(SIGMASQ)
	SIGFSQ=SCA*SIGMA*ELP/COR
25	FSQ=SCA*EMINB*ELP/COR
	IF(FSQ.EQ.0.0) GO TO 22
	AFSQ=SIGN(FSQ,1.)
	F=SIGN(SQRT(AFSQ),FSQ)
26	AF = SQRT(AFSQ)
	FSIG = SIGFSQ/(2.0*AF)
	GO TO 23
22	F=0.0
	FSIG=0.0
23	SCALEN=NSCOFF(NO)
	IF(SCALEN.EQ.1.)SCALEN=IQ
	 NWALT=NWALT+1
	WRITE (7)   (X(I,NO), I = 1, 3), OBSI(NO), SACK, SCA, ABS,
	1FSQ, SIGFSQ, F, FSIG, SINOL, TETA, OM(NO), ELP, SCALEN,
	2(XNOUT(I), I = 1, 3), TOUT, NSERIAL(NO), ABSBE,
	3(XSOUT(I), I = 1, 3)
27	CONTINUE
C
C     DETERMINE WETHER JOB IS FINISHED OR ADDITIONAL CARDS ARE TO BE
C     READ
C
270	IF (IDENT .LT. 1) GO TO 51
29	IF (1 .LE. IDENT .AND. IDENT .LE. 3) GO TO 30
	GO TO 300
30	CONTINUE
	ENDFILE 7
	REWIND 7
	IF(ICN(3).LE.1) GO TO 776
	IDIR = 2
	CALL SCALEB (H, K, L, ISIG, M, NSTD, IHSTD, IDIR, ICOUNT, FACTOR,
	1 ILIST)
	FLAST = START
	SUN = 0.0
	NUM = 0
776	NO = 1
777	READ(7,END=701) (X(I,NO), I = 1,3), OBSI(NO), SACK, SCA,ABS,
	1FSQ, SIGFSQ, F, FSIG, SINOL, TETA, OM(NO), ELP, SCALEN,
	2(XNOUT(I), I = 1, 3), TOUT, NSERIAL(NO), ABSBE,
	3(XSOUT(I), I = 1, 3)
702	IF(ICN(3).LE.1) GO TO 703
	IDIR = 3 + IDEGREE
	CALL SCALEB (H, K, L, ISIG, M, NSTD, IHSTD, IDIR, ICOUNT, FACTOR,
	1 ILIST)
	FACTOR = FACTOR*START
	SCA = FACTOR
	FSQ = FACTOR*FSQ
	SIGFSQ = FACTOR*SIGFSQ
	F = SQRT(FACTOR)*F
	FSIG = SQRT(FACTOR)*FSIG
	EMINB = (OBSI(NO) - SACK) * FACTOR
	IF(EMINB.LT.0.0) EMINB = - EMINB
	EMINS = SQRT(EMINB)
703	PRINT  510, (X(I,NO), I = 1, 3), OBSI(NO), SACK, SCA, ABS,
	1FSQ, SIGFSQ, F, FSIG, SINOL, TETA, OM(NO), ELP, SCALEN,
	2(XNOUT(I), I = 1, 3), TOUT,              ABSBE,
	3(XSOUT(I), I = 1, 3), NSERIAL(NO)
	NLINES = NLINES + 2
	IF (NLINES .LT. 51) GO TO 710
	PRINT 5081
	NLINES = 0
710	IF(ICN(3).LE.1) GO TO 713
	SUN = SUN + FACTOR
	NUM = NUM + 1
	IF (NUM .NE. 500) GO TO 711
	FLAST = SUN/NUM
	SUN = 0.0
	NUM = 0
711	IF (0.95 .LT. (FACTOR/FLAST) .AND.
	1              (FACTOR/FLAST) .LT. 1.05) GO TO 713
712	PRINT 1510
1510	FORMAT (1H+, 25X, 1H*)
713	IF (ICN(7) .GE. 0) WRITE (4) JDENT, (X(I,NO), I = 1, 3),
	1FSQ,SIGFSQ,SCALEN,(XNOUT(I),I=1,3),(XSOUT(I),I=1,3),
	2TOUT,NSERIAL(NO)
	IJK1 = X(1,NO)
	IJK2 = X(2,NO)
	IJK3 = X(3,NO)
	IF(ICN(6).EQ.2.OR.ICN(6).EQ.3) WRITE(2,511) IJK1,IJK2,IJK3,
	1 FSQ,SIGFSQ,EMINB,EMINS,SCALEN,(XNOUT(I),I=1,3),(XSOUT(I),I=1,3),
	2 TOUT,NSERIAL(NO)
511	FORMAT(1X3I3,F10.2,F8.2,F8.0,F8.2,F3.0,7F6.4,I7)
	IF(ICN(6).EQ.1.OR.ICN(6).EQ.3) PUNCH 512,   IJK1,IJK2,IJK3,
	1 FSQ,SIGFSQ,EMINB,EMINS,SCALEN,NSERIAL(NO)
512	FORMAT(1X3I3,4F10.2,2I10)
	GO TO 777
701	REWIND 7
	IJK1 = 998
	IF(IDENT.EQ.1) IT=1; GO TO 620
	IF(IDENT.EQ.3)IT=2; GO TO 620
	IJK1 = 999
	IF(ICN(6).LE.1) GO TO 32
	WRITE(2,511) IJK1
	ENDFILE2
	REWIND2
	GO TO 32
620	IF(ICN(6).EQ.2.OR.ICN(6).EQ.3) WRITE(2,511) IJK1
	IF(ICN(6).EQ.1.OR.ICN(6).EQ.3) PUNCH 512,IJK1
	GO TO (5,31),IT
900	FORMAT(I1)
32	 IF(ICN(7) .GE. 0) END FILE 4
	IF(ICN(6).EQ.1.OR.ICN(6).EQ.3) PUNCH 512,IJK1
	 IF(ICN(7) .GE. 0) REWIND4
	 IF(ICN(8) .NE. 0) REWIND 3
700	 PRINT 1000,NWALT
1000	FORMAT('0NUMBER OF REFLECTIONS AT THIS POINT IS'I6)
	PRINT 515
	CALL EXIT
500	FORMAT (20A4)
501	FORMAT (1H1, 20A4)
502	FORMAT(6F9.6)
503	FORMAT(16I5)
504	FORMAT(' ICN (CONTROL INTEGERS) = '16I5)
506	FORMAT(I5,7F10.5,F5.2)
510	FORMAT (3F4.0, F8.0, F6.0, 2F8.4, F10.2, F8.2, F8.2, F6.2, F7.4,
	1F6.2, F6.2, F6.3, F3.0, 4F6.4,       /, 34X, F8.4, 43X, 8X,
	29X, 3F6.4, I6)
513	FORMAT(9F8.5)
5081	FORMAT ('1  H   K   L   (I+B)     B  SCALED     ABS    FSQ   SIG(F
	1SQ)    F   SIG(F) SINT/L THETA OMEGA    LP IQ COMPS. OF NORMALS  T
	2BAR'        / 36X, 'BE ABS', 47X, 6X,     30X,'NO.'/)
515	FORMAT(1H0,30X,40HT H I S  I S  T H E  E N D.              )
	END
	SUBROUTINE AMI3(A,B,DET)
C     INVERTS ASYMMETRIC 3X3 MATRIX
	DIMENSION A(3,3),B(3,3)
	B(1,1)=A(2,2)*A(3,3)-A(3,2)*A(2,3)
	B(2,1)=A(2,1)*A(3,3)-A(3,1)*A(2,3)
	B(3,1)=A(2,1)*A(3,2)-A(3,1)*A(2,2)
	DET=A(1,1)*B(1,1)-A(1,2)*B(2,1)+A(1,3)*B(3,1)
	IF(DET)1,5,1
1	B(1,2)=A(1,2)*A(3,3)-A(3,2)*A(1,3)
	B(2,2)=A(1,1)*A(3,3)-A(3,1)*A(1,3)
	B(3,2)=A(1,1)*A(3,2)-A(3,1)*A(1,2)
	B(1,3)=A(1,2)*A(2,3)-A(2,2)*A(1,3)
	B(2,3)=A(1,1)*A(2,3)-A(2,1)*A(1,3)
	B(3,3)=A(1,1)*A(2,2)-A(2,1)*A(1,2)
	DO 2 I=1,3
	DO 2 J=1,3
2	B(I,J)=(-1.)**(I+J)*B(I,J)/DET
5	RETURN
	END
	SUBROUTINEDSTAR(H,K,L,RCP,DST)
	DIMENSIONRCP(6)
	REALH,K,L
	DST=(H*RCP(1))**2+(K*RCP(2))**2+(L*RCP(3))**2+2.*(H*K*RCP(1)*
	1RCP(2)*RCP(6)+H*L*RCP(1)*RCP(3)*RCP(5)+K*L*RCP(2)*RCP(3)*RCP(4))
	DST=SQRT(DST)
	RETURN
	END
	SUBROUTINEANGLES(OMT,FI,W,DST ,RCP,PHI,CHI,H1,H2,H3)
C
C     CALCULATES DIFFRACTOMETER SETTING ANGLES GIVEN OMEGA
C     CHECK FIRST FOR IMPOSSIBLE OMEGA VALUES
	DIMENSION X(3),RCP(6),W(3,3),DXYZ(3)
	DATA RAD/57.2957795/
	COMMON/ANGDIR/COSCHI,SINCHI,COSOMT,SINOMT,COSPHI,SINPHI
	X(1)=H1;X(2)=H2;X(3)=H3
	DO2I=1,3
	DXYZ(I)=0.
	DO1J=1,3
1	DXYZ(I)=DXYZ(I)+W(I,J)*X(J)
2	DXYZ(I)=DXYZ(I)/DST
	 COSOMT=COS(OMT/RAD)
	 SINOMT=SIN(OMT/RAD)
C
C     CHI LIMITED BETWEEN + AND - 90
C
	SINCHI=-DXYZ(3)/COSOMT
	IF(SINCHI.GT.1.0.OR.SINCHI.LT.-1.0)SINCHI=SIGN(1.0,SINCHI)
	 COSCHI=SQRT(1.-SINCHI**2)
	 CHI=ASIN(SINCHI)*RAD
	Q=SINOMT;R=COSOMT*COSCHI
	IF(Q.EQ.0.AND.R.EQ.0.)GOTO3
	COSPHI = (-Q*DXYZ(1) + R*DXYZ(2))/(Q**2 + R**2)
	SINPHI = -(DXYZ(1)-Q*COSPHI)/R
C
C          ADD NORMALIZATION TO MINIMIZE ERRORS NEAR CHI=90
C
	T=COSPHI**2 + SINPHI**2
	T=1./SQRT(T)
	COSPHI=COSPHI*T  ;  SINPHI=SINPHI*T
	 PHI=ACOS(COSPHI)*RAD
	 PHI=SIGN(PHI,SINPHI)
	RETURN
C     REFLECTIONS AT CHI =90
3	PHI=FI
	 COSPHI=COS(FI/RAD)
	 SINPHI=SIN(FI/RAD)
	RETURN
	END
	SUBROUTINE RDCL (X,Y,V,FLAME,LI,LO)
C
C          X(6) = REAL CELL CONSTANTS
C          Y(6) = RECIPROCAL CELL CONSTANTS
C          V = REAL CELL VOLUME
C          FLAME = WAVELENGTH, ANGSTROMS
C          UPON RETURN, AA(3,3) IS REAL METRIC, BB(3,3) IS RECIPROCAL
C               METRIC
C          LI = INPUT TAPE NUMBER;  LO = OUTPUT TAPE NUMBER
C          AFTER C. K. JOHNSON, ORTEP
C
	DIMENSION A(9),AA(3,3),B(9),BB(3,3),X(6),Y(6)
	COMMON /METRICS/ AA,BB
	READ(LI,1) (A(I),I=1,6),FLAME
1	FORMAT(7F9.6)
	TEMP=ABS(A(4))-1.
	DO 5 J=1,3
	 IF(TEMP) 3,2,2
C
C**** CELL ANGLES IN DEGREES*****
C
2	 A(J+6)=A(J+3)
	 A(J+3)=COS(A(J+6)*0.01745329252)
	 GOTO5
C
C****CELL ANGLE COSINES*****
C
3	 A(J+6)=ARCCOS(A(J+3))
5	 AA(J,J)=A(J)**2
	AA(2,1)=AA(1,2)=A(1)*A(2)*A(6)
	AA(3,1)=AA(1,3)=A(1)*A(3)*A(5)
	AA(3,2)=AA(2,3)=A(2)*A(3)*A(4)
C
C          CALL TO INVERT 3X3 MATRIX****MAY BE CHANGED FOR PARTICULAR
C          PROGRAM
C
	CALL SMI3(AA,BB,V)
C
C
C****CALC RECIPROCAL VALUES FROM BB*****
C
	DO 6 J=1,3
6	 B(J)=SQRT(BB(J,J))
	B(6)=BB(1,2)/(B(1)*B(2))
	B(5)=BB(1,3)/(B(1)*B(3))
	B(4)=BB(2,3)/(B(2)*B(3))
	DO 7 J=1,3
7	 B(J+6)=ARCCOS(B(J+3))
C
C****CHECK TO SEE IF INPUT WAS REAL OR RECIPROCAL
C
	IF(A(1)-1.) 8,10,10
8	V=1./V
	DO 9 J=1,9
	 TEMP=AA(J,1)
	 AA(J,1)=BB(J,1)
	 BB(J,1)=TEMP
	TEMP=A(J)
	 A(J)=B(J)
9	 B(J)=TEMP
10	V=SQRT(V)
	WRITE (LO,11)
11	FORMAT('0DIRECT CELL PARAMETERS'/1H ,15X,'A'14X'B'14X'C'13X'ALPHA'
	'     10X'BETA'11X'GAMMA')
	WRITE(LO,12) (A(I),I=1,3),(A(I),I=7,9),(A(I),I=4,6)
12	FORMAT(1H ,10X,F9.6,2F15.6,3F15.3/1H ,48X'COSINE'F12.8,2F15.8)
	WRITE(LO,13)
13	FORMAT('0RECIPROCAL CELL PARAMETERS'/1H ,15X,2HA*,13X,2HB*,13X,
	*     2HC*,12X,6HALPHA*,9X,5HBETA*,10X,6HGAMMA*)
	WRITE(LO,12) (B(I),I=1,3),(B(I),I=7,9),(B(I),I=4,6)
	DO 14 J=1,6
	 X(J)=A(J)
14	 Y(J)=B(J)
	RETURN
	END
	FUNCTION ARCCOS(X)
C
C****ARCCOS(X) IN DEGREES*****
C
	IF(1.-ABS(X))1,2,2
1	X=SIGN(1.,X)
2	IF(X) 3,4,5
3	ARCCOS=180. + ATAN(SQRT(1.-X*X)/X)*57.29577951
	GOTO6
4	ARCCOS=90.
	GOTO6
5	ARCCOS=ATAN(SQRT(1.-X*X)/X)*57.29577951
6	RETURN
	END
	SUBROUTINE SMI3(A,AI,D)
C
C          SYMMETRICAL (3,3) MATRIX INVERTER
C          A IS MATRIX, AI IS INVERSE, D IS VOLUME (=DETERMINANT OF A)
C
	DIMENSION A(3,3),AI(3,3)
	DO 1 I=1,3
	 K=1+MOD(I+1,3)
	 M=1+MOD(I+3,3)
	 DO 1 J=1,3
	  L=1+MOD(J+1,3)
	  N=1+MOD(J+3,3)
1	  AI(I,J)=A(L,K)*A(N,M)-A(L,M)*A(N,K)
	D=0.0
	DO 2 I=1,3
2	 D=D+A(I,1)*AI(1,I)
	DO 3 I=1,3
	 DO 3 J=1,3
3	  AI(I,J)=AI(I,J)/D
	RETURN
	END
	SUBROUTINE UMATRIX (W,ELAMDA)
C
C          ROUTINE TO CALCULATE ORIENTATION MATRIX (BEDNOWITZ-COPPENS
C          DEFINITION) FROM THREE REFERENCE REFLECTIONS AND ANGLE ZEROES
C          REFERENCE REFLECTIONS AND ANGLE ZEROES
C                       A*(TRANSPOSE) = A(G) * W
C
C          FORMULAS ARE CONSISTENT FOR A DIFFRACTOMETER WITH
C          FOLLOWING PARITIES:
C                  OMEGA POSITIVE CLOCKWISE
C                  PHI   POSITIVE CLOCKWISE
C                  CHI   POSITIVE COUNTERCLOCKWISE
C
	DIMENSION PHI(3),OM(3),CHI(3),TWOTH(3),H(3,3),HH(3,3),W(3,3)
	DIMENSION X(3,3),RCP(6),DST(3)
	COMMON /ZEROES/ PHIZRO,OMZRO,TWOTHZO,CHIZRO
	COMMON RCP
	DATA RADIAN /57.2957794636/
	LI = 105
	LO = 108
	READ (LI,1) PHIZRO,OMZRO,TWOTHZO,CHIZRO
1	FORMAT(4F9.4)
	WRITE(LO,2) PHIZRO,OMZRO,TWOTHZO,CHIZRO
2	FORMAT('0REFERENCE ANGLE ZEROES'/'       PHI ZERO = 'F10.4
	'      /'     OMEGA ZERO = 'F10.4/' TWO THETA ZERO = 'F10.4/
	'     '       CHI ZERO = 'F10.4/)
	READ (LI,3) (((H(I,J),J=1,3),PHI(I),OM(I),TWOTH(I),CHI(I)),I=1,3)
3	FORMAT(3F9.0,4F9.4)
	WRITE (LO,4)(((H(I,J),J=1,3),PHI(I),OM(I),TWOTH(I),CHI(I)),I=1,3)
4	FORMAT(/' WN MATRIX DETERMINED FROM THREE REFERENCE REFLECTIONS,
	'(OBSERVED VALUES)'
	'   /'   H   K   L     PHI        OM    TWO THETA    CHI'
	'     /3(3F4.0,4F10.4/))
	DO 5 I=1,3
	 PHI(I)=(PHI(I)-PHIZRO)/RADIAN
	 OM(I)=(OM(I)-OMZRO)/RADIAN
	 TWOTH(I)=(TWOTH(I)-TWOTHZO)/RADIAN
5	 CHI(I)=(CHI(I)-CHIZRO)/RADIAN
	DO 8 I=1,3
	 CALL DSTAR(H(I,1),H(I,2),H(I,3),RCP,DST(I))
	 TTHETA=ASIN(ELAMDA*.5*DST(I))*RADIAN*2.
	 WRITE (LO,6) I,DST(I),TTHETA
6	FORMAT  ('0VECTOR('I1')    DSTAR = 'F10.6'     TWO THETA (CALC) ='
	*  F10.4  )
	 X(I,1)= -DST(I)*(COS(PHI(I))*SIN(OM(I)) + SIN(PHI(I))*
	*     COS(CHI(I))*COS(OM(I)))  + .0
	 X(I,2)=DST(I)*(-SIN(PHI(I))*SIN(OM(I))+COS(PHI(I))*
	*     COS(CHI(I))*COS(OM(I)) )  + .0
	 X(I,3)=-DST(I)*SIN(CHI(I))*COS(OM(I)) + .0
	 WRITE (LO,7) (X(I,J),J=1,3)
7	FORMAT(1H ,12X,'( X(1) = 'F10.6')'/1H ,12X'( X(2) = 'F10.6')'/
	'     1H ,12X'( X(3) = 'F10.6')'/)
8	CONTINUE
	 CALL AMI3(H,HH,DET)
	WRITE (LO,9) ((HH(I,J),J=1,3),I=1,3)
9	FORMAT ('0HH MATRIX (INVERSE OF H TRIPLES) IS '/3(3F10.6/))
C
C          HERE CALCULATE W TRANSPOSE (W(J,I)) RATHER THAN W ( W(I,J) )
C          BEAUSE LATTER WOULD INVOLVE TRANSPOSING BOTH X AND HH FIRST
C
	DO 12 I=1,3
	 DO 12 J=1,3
	  W(J,I)=0.
	  DO 11 K=1,3
11	   W(J,I)=W(J,I)+HH(I,K)*X(K,J)
12	  W(J,I)=W(J,I)/RCP(I)
C
	RETURN
	END
	SUBROUTINE GRIT2(GRID,MDGRID,NDGRID)
	DIMENSION GRID(MDGRID,NDGRID)
	REAL*8 GAUSS(128,2),ERICH(256)
C     CONSTRUCTS GAUSSIAN GRID WITHIN CRYSTAL. GRID AXES PARALLEL TO XTA
	DIMENSION RCP(6),RECIP(3,3),DIRECT(3,3),DRC(6),SHUFEL(3,3)
	DIMENSION X(3,10),OBSI(10),BACK(10),FI(10),EPOA(3000)
	DIMENSION NGAUSS(3),EPOB(1000),DATA(40,4),SCB(10)
	DIMENSION SUM(40),TT(3),E(3),W(3,3)
	DIMENSION NADRES(3),ABOUND(3,3),DIST(3),XYZ(3),XX(50)
	DIMENSION BBOUND(2,2),DI(2),TEMP(3)
	DIMENSION CBOUND(3,3),DBOUND(2,2)
	DIMENSION ABCD(40,4)
	DIMENSION GU(3,3),XNZ(3)
	COMMON RCP,RECIP,DIRECT,DRC,SHUFEL,X,OBSI,BACK,FI,EPOA,EPOB
	COMMON ELAMDA,SINSQT,COST,SINT,EMU,NPOINTS,NR
	COMMON     EXT,ABS,DATA,NO,NFACES
	COMMON NGAUSS,VOL,TBAR
	COMMONSUM,W,TT,E
	COMMON /CN/ IDIST
	COMMON /DATDIR/ ICN(16)
	COMMON /DST/ XOZ(3,100)
	EQUIVALENCE (ICN(1),IABS),(ICN(2),IEXT)
	DATA IA,IB,IC/1,2,3/
	 EQUIVALENCE(GAUSS,ERICH)
C       FOLLOWING ARE THE TABLES FOR THE GAUSSIAN INTEGRATION
	 DATA (ERICH(I),I=1,64)  /
	1                 .2113248660,.7886751350,.0694318440,.3300094780,
	1                 .6699905220,.9305681560,.0337652430,.1693953070,
	1                 .3806904070,.6193095930,.8306046930,.9662347570,
	1                 .0198550720,.1016667610,.2372337950,.4082826790,
	1                 .5917173210,.7627662050,.8983332390,.9801449280,
	1                 .0130467360,.0674683170,.1602952160,.2833023030,
	1    .4255628310,.5744371700,.7166976970,.8397047840,
	1                 .9325316830,.9869532640,.0092196830,.0479413720,
	1                 .1150486630,.2063410230,.3160842510,.4373332960,
	1                 .5626667040,.6839157500,.7936589770,.8849513370,
	1                 .9520586280,.9907803170,.0068580960,.0357825580,
	1                 .0863993420,.1563535480,.2423756820,.3404438160,
	1                 .4459725260,.5540274740,.6595561840,.7576243180,
	1                 .8436464520,.9136006580,.9642174420,.9931419040,
	1                 .0052995330,.0277124880,.0671843990,.1222977960,
	1                 .1910618780,.2709916110,.3591982250,.4524937450/
	 DATA (ERICH(I),I=65,128) /
	1                 .5475062550,.6408017750,.7290083890,.8089381220,
	1                 .8777022040,.9328156010,.9722875120,.9947004680,
	1                 .9975936100,.9873642780,.9691372760,.9432077640,
	1                 .9100009930,.8700620960,.8240468260,.7727107360,
	1                 .7168967540,.6575213400,.5955594340,.5320284470,
	1                 .4679715530,.4044405660,.3424786600,.2831032460,
	1                 .2272892640,.1759531740,.1299379040,.0899990070,
	1                 .0567922360,.0308627240,.0126357220,.0024063900,
	1                 .9986319310,.9928057560,.9823811280,.9674530380,
	1                 .9481605780,.9246838070,.8972418980,.8660910600,
	1                 .8315221340,.7938578790,.7534499550,.7106756380,
	1                 .6659343010,.6196436810,.5722359810,.5241538330,
	1                 .4758461670,.4277640190,.3803563190,.3340656990,
	1                 .2893243620,.2465500450,.2061421210,.1684778670,
	1                 .1339089400,.1027581020,.0753161930,.0518394220,
	1                 .0325469620,.0176188720,.0071942440,.0013680690/
	 DATA (ERICH(I),I=129,192) /
	1                 .5000000000,.5000000000,.1739274230,.3260725780,
	1                 .3260725780,.1739274230,.0856622460,.1803807870,
	1                 .2339569670,.2339569670,.1803807870,.0856622460,
	1                 .0506142680,.1111905170,.1568533230,.1813418920,
	1                 .1813418920,.1568533230,.1111905170,.0506142680,
	1                 .0333356720,.0747256750,.1095431810,.1346333600,
	1                 .1477621120,.1477621120,.1346333600,.1095431810,
	1                 .0747256750,.0333356720,.0235876680,.0534696630,
	1                 .0800391640,.1015837130,.1167462680,.1245735230,
	1                 .1245735230,.1167462680,.1015837130,.0800391640,
	1                 .0534696630,.0235876680,.0175597300,.0400790440,
	1                 .0607592850,.0786015840,.0927691990,.1025992320,
	1                 .1076319270,.1076319270,.1025992320,.0927691990,
	1                 .0786015840,.0607592850,.0400790440,.0175597300,
	1                 .0135762300,.0311267620,.0475792560,.0623144860,
	1                 .0747979940,.0845782600,.0913017080,.0947253050/
	 DATA (ERICH(I),I=193,256) /
	1                 .0947253050,.0913017080,.0845782600,.0747979940,
	1                 .0623144860,.0475792560,.0311267620,.0135762300,
	1                 .0061706150,.0142656940,.0221387190,.0296492920,
	1                 .0366732410,.0430950810,.0488093260,.0537221350,
	1                 .0577528340,.0608352360,.0629187280,.0639690980,
	1                 .0639690980,.0629187280,.0608352360,.0577528340,
	1                 .0537221350,.0488093260,.0430950810,.0366732410,
	1                 .0296492920,.0221387190,.0142656940,.0061706150,
	1                 .0035093050,.0081371970,.0126960330,.0171369320,
	1                 .0214179490,.0254990300,.0293420470,.0329111110,
	1                 .0361728970,.0390969480,.0416559620,.0438260470,
	1                 .0455869390,.0469221990,.0478193600,.0482700440,
	1                 .0482700440,.0478193600,.0469222000,.0455869390,
	1                 .0438260470,.0416559620,.0390969480,.0361728970,
	1                 .0329111110,.0293420470,.0254990300,.0214179490,
	1                 .0171369320,.0126960330,.0081371970,.0035093050/
	READ102,(NGAUSS(I),I=1,3),NFACES
102	FORMAT(4I5)
	IGRID1=1
	IGRID=4*NGAUSS(1)*NGAUSS(2)*NGAUSS(3) + 1
	IGRID1=LOCF(GRID(1,IGRID1))
	IGRID=LOCF(GRID(4,IGRID))
	IT = IGRID - IGRID1
	PRINT 103,     IGRID1,IGRID,IT
103	FORMAT('-NECESSARY FIELD LENGTH FOR GRID IS FROM 'M6' TO 'M6,
	1/' OR ' ,I6,' DECIMAL WORDS OF MEMORY')
	CALL MM (W,DIRECT,GU)
	IAE = IABS
	IF (IABS.EQ.0) IAE = IEXT
	CALL FANGLE (IAE,NFACES,W,DATA,ABCD)
C     DETERMINE LENGTHS OF RECIPROCAL VECTORS PERPENDICULAR TO BOUNDARY
	DO 50 I=1,NFACES
	DO 51 J=1,3
	TEMP(J)=0.
	DO 51 K=1,3
51	TEMP(J)=TEMP(J)+DATA(I,K)*RECIP(K,J)
	SUM(I)=0.
	DO 53 L=1,3
53	SUM(I)=SUM(I)+TEMP(L)*DATA(I,L)
	 SUM(I)=SQRT(SUM(I))
C     MULTIPLY DISTANCE FROM ORIGIN TO FACE BY LENGTH OF RECIPROCAL VECT
C     NORMAL TO FACE,INSTEAD OF NORMALIZING THIS VECTOR TO UNITY BY DIVI
C     H,K AND L BY LENGTH
	DATA(I,4)=DATA(I,4)*SUM(I)
50	CONTINUE
C     FIND FIRST ADDRESS OF GAUSSIAN CONSTANTS TO BE USED
	DO 1 I=1,3
	K=NGAUSS(I)/2
	IF(K-8)54,54,55
54	NADRES(I)=(K-1)*K
	GO TO 1
55	NADRES(I)=72+(K-12)*6
1	CONTINUE
	MMM=0
	NPOINTS=0
	VOL=0.
C     SCAN OVER ALL POSSIBLE COMBINATIONS OF THREE FACES TO DETERMINE CO
C     TES OF CORNERS OF CRYSTAL
2	N1=NFACES-1
	N2=NFACES-2
	PRINT 110
110	FORMAT('-POINTS OF INTERSECTION OF TRANSFORMED PLANES'//10X'CRYSTA
	'L SYSTEM'68X'ORTHOGONAL SYSTEM (CM)'/'       X         Y         Z
	''35X'FACE NO.       VERTEX' 8X'X' 9X'Y' 9X'Z'/)
	DO 3 K=1,N2
	DO 40 I=1,3
40	ABOUND(1,I)=DATA(K,I)
	K1=K+1
	DO 3 L=K1,N1
	DO 41 I=1,3
41	ABOUND(2,I)=DATA(L,I)
	L1=L+1
	DO 3 M=L1,NFACES
	DO4 I=1,3
4	ABOUND(3,I)=DATA(M,I)
	CALL AMI3(ABOUND,CBOUND,DET)
	IF(DET)5,3,5
5	DIST(1)=DATA(K,4)
	DIST(2)=DATA(L,4)
	DIST(3)=DATA(M,4)
	CALL MV(CBOUND,DIST,XYZ)
C     CHECK WHETHER POINT IS ON CRYSTAL SURFACE, IF OUTSIDE CRYSTAL REJE
C     XYZ ARE FRACTIONAL COORDINATES OF THE POINT OF INTERSECTION OF 3 P
C     SINCE THE DISTANCES FROM THE ORIGIN TO EACH OF THE FACES ARE GIVEN
C     THE COORDINATES XYZ ARE FRACTIONS OF AXES A,B,C OF LENGTH A(CM),B(
	CALL INSIDE(XYZ,NOUT,NFACES,DATA)
	IF(NOUT)3,6,3
6	MMM=MMM+1
	CALL MV(GU,XYZ,XNZ)
	PRINT 111,  XYZ,K,L,M,MMM,XNZ
111	FORMAT (3F10.5,5X ' INTERSECTION OF PLANES '3I4,5X,I5,5X,3F10.6)
	XX(MMM)=XYZ(IA)
	DO 37 IJ=1,3
37	 XOZ(IJ,MMM)= XYZ(IJ)
3	CONTINUE
	IF (IDIST.EQ.0)GOTO9
	CALL DISTA(DIRECT)
9	CALL MINMAX(XX,MMM,XMIN,XMAX)
C     DETERMINE LOWER AND UPPER LIMITS IN X
	NIA =NGAUSS(IA)
10	DO11NPA=1,NIA
C     XYZ(IA) ARE FRACTIONAL COORDINATES OF GRID POINTS ALONG A AXIS
	IF(NIA-1)21,21,22
22	XYZ(IA)=XMIN+(XMAX-XMIN)*GAUSS((NADRES(IA)+NPA),1)
C     R1 IS THE GAUSSIAN WEIGHT CORRESPONDING TO XYZ(IA)
	R1= GAUSS((NADRES(IA)+NPA),2)*(XMAX-XMIN)
21	MMM=0
C     SCAN OVER ALL POSSIBLE COMBINATIONS OF TWO FACES TO DETERMINE COOR
C     OF EDGES OF CRYSTAL FOR A GIVEN X
	DO 12 K=1,N1
	BBOUND(1,1)=DATA(K,IB)
	BBOUND(1,2)=DATA(K,IC)
	DI(1)= DATA(K,4)-DATA(K,IA)*XYZ(IA)
	K1=K+1
	DO 12 L=K1,NFACES
	BBOUND(2,1)=DATA(L,IB)
	BBOUND(2,2)=DATA(L,IC)
	CALL AMI2(BBOUND,DBOUND,DET)
	IF(DET)13,12,13
13	DI(2)= DATA(L,4)-DATA(L,IA)*XYZ(IA)
	XYZ(IB)=DBOUND(1,1)*DI(1)+DBOUND(1,2)*DI(2)
	XYZ(IC)=DBOUND(2,1)*DI(1)+DBOUND(2,2)*DI(2)
	CALL INSIDE(XYZ,NOUT,NFACES,DATA)
	IF(NOUT)12,14,12
14	MMM=MMM+1
	XX(MMM)=XYZ(IB)
12	CONTINUE
C     DETERMINE LOWER AND UPPER LIMITS IN Y FOR GIVEN X
	CALL MINMAX(XX,MMM,YMIN,YMAX)
	NIB=NGAUSS(IB)
	DO11 NPB=1,NIB
	XYZ(IB)=YMIN+(YMAX-YMIN)*GAUSS((NADRES(IB)+NPB),1)
	R2=GAUSS((NADRES(IB)+NPB),2)*R1*(YMAX-YMIN)
	MMM=0
C     SCAN OVER ALL FACES TO DETERMINE UPPER AND LOWER LIMITS OF Z FOR G
	DO 15 I=1,NFACES
	IF(DATA(I,IC))16,15,16
16	XYZ(IC)=(DATA(I,4)-XYZ(IA)*DATA(I,IA)-XYZ(IB)*DATA(I,IB))/DATA(I,I
	1C)
	CALL INSIDE(XYZ,NOUT,NFACES,DATA)
	IF(NOUT)15,17,15
17	MMM=MMM+1
	XX(MMM)=XYZ(IC)
15	CONTINUE
	CALL MINMAX(XX,MMM,ZMIN,ZMAX)
	NIC=NGAUSS(IC)
	DELTAZ=ZMAX-ZMIN
	DO11 NPC=1,NIC
	NPOINTS=NPOINTS+1
C     STORE FRACTIONAL COORDINATES OF SAMPLING POINTS AND CORRESPONDING
	GRID(IC,NPOINTS)=ZMIN+DELTAZ*GAUSS((NADRES(IC)+NPC),1)
	GRID(IB,NPOINTS)=XYZ(IB)
	GRID(IA,NPOINTS)=XYZ(IA)
	GRID(4,NPOINTS)=GAUSS((NADRES(IC)+NPC),2)*R2*DELTAZ
C     VOL IS EQUAL TO THE VOLUME OF THE CRYSTAL AS A FRACTION OF THE VOL
C     A PARALLELOPIPED WHOSE EDGES ARE PARALLEL TO A,B,C AND OF LENGTH A
C     B(CMS),C(CMS)
	VOL=VOL+GRID(4,NPOINTS)
11	CONTINUE
	PRINT106,(NGAUSS(I),I=1,3),NPOINTS
106	FORMAT('-N(1) = 'I2'  N(2) = 'I2'  N(3) = 'I2'          NUMBER OF
	'SAMPLING POINTS 'I6/)
	 SINIA=SQRT(1.-DRC(IA+3)**2)
	CLES=DRC(IB)*DRC(IC)*SINIA/RCP(IA)
	VOLUME=VOL*CLES*1000.
	PRINT 104,VOLUME
104	FORMAT('0VOLUME OF CRYSTAL IS 'E13.5,6H MM**3)
	RETURN
	END
	SUBROUTINE MM(X,Y,Z)
C
C          MULTIPLY TWO MATIRCES
C          Z = X * Y
C
	DIMENSION X(3,3),Y(3,3),Z(3,3)
	DO 1 I=1,3
	 DO 1 K=1,3
	  Z(I,K)=0.
	  DO 1 J=1,3
1	   Z(I,K)=Z(I,K)+X(I,J)*Y(J,K)
	RETURN
	END
	SUBROUTINE VM(X,Y,Z)
C
C          VECTOR = VECTOR TIMES MATRIX
C          Z(3) = X(3) * Y(3,3)
C
	DIMENSION Z(3),X(3),Y(3,3)
	DO 1 J=1,3
	 Z(J)=0.
	 DO 1 I=1,3
1	  Z(J)=Z(J)+X(I)*Y(I,J)
	RETURN
	END
	SUBROUTINE FANGLE(IABS,NFACES,W,DATA,ABCD)
C
C          READS   H,K,L,D(CM) FOR     IABS = 1
C                  PHI,OMEGA,CHI FOR   IABS = 2
C          OR      A,B,C,D(CM) FOR     IABS = 3
C          AND CALCULATES  H,K,L AND A,B,C, PRINTS BOTH, INCLUDING D,
C          AND SETTING ANGLES FOR REFERENCE.
C
C          PROGRAM TOMPAS REQUIRES A,B,C,D, REFERRED TO -X,-Y,Z SYSTEM
C          PROGRAM DATAPP REQUIRES H,K,L,D, REFERRED TO X,Y,Z SYSTEM
C                    ( X,Y,Z, SYSTEM IS THAT DEFINED AT BNL, SEE
C                    WRITE-UP)
C
	DIMENSION T(3,3),WI(3,3),T1(3),TI(3,3),DATA(40,4),ABCD(40,4)
	DIMENSION WTI(3,3),W(3,3),TEMP(3,3),WT(3,3),T2(3)
	COMMON RCP(6)
	COMMON /ZEROES/ PHIZRO,OMZRO,TWOTHZO,CHIZRO
	DATA RAD/57.2957794636/
	DATA TI /-1.,0.,0.,0.,-1.,0.,0.,0.,1./
	EQUIVALENCE (TI,T)
	LI = 105
	LO = 108
	DO 1 I=1,3
	 DO 1 J=1,3
1	  TEMP(I,J)=W(J,I)
	CALL MM(TEMP,TI,WT)
	GOTO(2,11,21),IABS
C
C          HERE READ H,K,L,D(CM) IN BNL (X,Y,Z) SYSTEM
C          CALCULATE P = (A,B,C) AS H/(D*) . W(TRANS) . TI  (EQN. 4.4B,
C               EDMONDS DATAPP NOTES)
C          WHERE TI = T = T(TRANS), AND TRANS = TRANSPOSE
C          W(TRANS) IS TRANSPOSE OF W USED IN MAIN PROGRAM DATAPP
C
2	READ(LI,3) ((DATA(I,J),J=1,4),I=1,NFACES)
3	FORMAT(4F10.6)
	WRITE (LO,4)
4	FORMAT( '-CRYSTAL BOUNDARY PLANES:  H,K,L REFER TO X,Y,Z (BNL) SYS
	'TEM, AND A,B,C REFER TO -X,-Y,Z (SEATTLE) SYSTEM'//'0FACE NO.   H
	'         K          L       DIST(CM)    DSTAR         A         B
	'        C           PHI(O)   OMEGA(O)   CHI(O)'/)
	DO 9 I=1,NFACES
	 CALL DSTAR(DATA(I,1),DATA(I,2),DATA(I,3),RCP,DST)
	 DO 5 J=1,3
5	  TEMP(J,1)=DATA(I,J)/DST
	 CALL VM(TEMP(1,1),WT,TEMP(1,2))
	 SUM=0.
	 DO 6 J=1,3
	  SUM=SUM+TEMP(J,2)**2
6	  ABCD(I,J)=TEMP(J,2)
	 SUM=SQRT(SUM)
	 DO 7 J=1,3
7	  ABCD(I,J)=ABCD(I,J)/SUM
	 ABCD(I,4)=10.*DATA(I,4)
C
	 OMEGA=FI=0.
	 CALL ANGLES (OMEGA,FI,W,DST,RCP,PHI,CHI,DATA(I,1),DATA(I,2),
	*      DATA(I,3) )
	 PHI=PHI+PHIZRO
	 OMEGA=OMEGA+OMZRO
	 CHI=CHI+CHIZRO
C
	 WRITE (LO,8)I,(DATA(I,J),J=1,4),DST,(ABCD(I,J),J=1,3),PHI,OMEGA,
	*     CHI
8	 FORMAT(I5,3(1X,F10.6),2X,1H*,F8.4,2X,F9.6,2X,1H*,1X,3(F10.7),
	*   2X,1H*,1X,2F9.4,F10.4)
9	CONTINUE
	RETURN
C
C          HERE READ OBSERVED SETTING ANGLES TO BRING /FACE/ INTO
C          DIFFRACTING CONDITION, I.E., NORMAL IS ALONG THETA-TWO THETA
C          BISECTOR.
C
C          FAKE H,K,L
C
C          TEMP USED HERE BOTH AS VECTOR(1X3) AND AS MATRIX (3X3)
C          WTI HERE IS W TRANSPOSE INVERSE
C
11	CONTINUE
C
	CALL AMI3 (TEMP,WTI,DET)
	WRITE(LO,13)
13	FORMAT('-CRYSTAL BOUNDARY PARAMETERS FROM OBSERVED SETTING ANGLES(
	'SEE ANGLE ZEROES ABOVE)'//'0FACE NO.   H          K          L
	'   DIST(CM)    DSTAR         A         B         C           PHI(O
	')   OMEGA(O)   CHI(O)'/)
	DO 20 I=1,NFACES
	 READ (LI,3) PHI,OMEGA,CHI,DATA(I,4)
C
C          CALCAULATE COMPONENTS OF X
C
	 XPHI=(PHI-PHIZRO)/RAD
	 XOMEGA=(OMEGA-OMZRO)/RAD
	 XCHI=(CHI-CHIZRO)/RAD
C
C          FOLLOWING VALUES OF T1 ARE NORMALIZED (= H*(1/D*) * WT )
C
	 T1(1) = -COS(XPHI)*SIN(XOMEGA) - SIN(XPHI)*COS(XCHI)*COS(XOMEGA)
	 T1(2) = -SIN(XPHI)*SIN(XOMEGA) + COS(XPHI)*COS(XCHI)*COS(XOMEGA)
	 T1(3) = -SIN(XCHI)*COS(XOMEGA)
C
	 SUM=0.
	 DO 15 J=1,3
15	  SUM=SUM + T1(J)**2
	 DST=SQRT(SUM)
C
C          HERE T1(3) ARE IDENTICALLY -A,-B,C, AND ARE TRANSFORMED TO
C          A,B,C BY TI
C
	 CALL VM(T1,TI,TEMP(1,1))
C
	 CALL VM(T1,WTI,TEMP(1,2))
C
C          (TEMP(J,2),J=1,3) NOW ARE ( H*(1/D*) )
C
C          HERE, FAKE H,K,L BY PICKING LOWEST NON-ZERO VALUE OF H*(1/D*)
C          (=TEMP(J,2),J=1,3) AND USING THAT AS LOWEST POSSIBLE VALUE OF
C          AN INDIVIDUAL H(I)*(1/D*).     ASSUMES THAT NO POSSIBLE
C          VALUE OF H/(D*) CAN BE LESS THAN .1
C
	 DO 16 J=1,3
	  TEMP(J,3)=ABS(TEMP(J,2))
	  IF (TEMP(J,3).GT..1) GOTO16
	  TEMP(J,3)=100000.
	  TEMP(J,2)=0.
16	 CONTINUE
	 T1MIN=AMIN1(TEMP(1,3),TEMP(2,3),TEMP(3,3) )
	 DO 17 J=1,3
17	 TEMP(J,3)=TEMP(J,2)/T1MIN
C
C          NORMALIZE (TEMP(J,1),J=1,3) AND STORE IN (ABCD(I,J),J=1,3)
C          STORE (TEMP(J,2),J=1,3) (=H*(1/D*) ) IN (DATA(I,J),J=1,3)
C
	 SUM=0.
	 DO 18 J=1,3
	 DATA(I,J)=TEMP(J,2)
	 ABCD(I,J)=TEMP(J,1)
18	  SUM=SUM+ABCD(I,J)**2
	 SUM=SQRT(SUM)
	 DO 19 J=1,3
19	  ABCD(I,J)=ABCD(I,J)/SUM
C
	 ABCD(I,4)=10.*DATA(I,4)
	 CALL DSTAR(TEMP(1,3),TEMP(2,3),TEMP(3,3),RCP,DST)
C
C          WRITE FAKE H,K,L VALUES HERE.  NORMALIZED VALUES PASSED TO
C          PROGRAMS IN ARRAY (DATA(I,J) )
C
	 WRITE (LO,8) I,(TEMP(J,3),J=1,3),DATA(I,4),DST,(ABCD(I,J),J=1,3),
	*     PHI,OMEGA,CHI
C
20	CONTINUE
	RETURN
C
C          READ A,B,C,D(MM) REFERRED TO -X,-Y,Z (SEATTLE) AND CALCULATE
C          H*(1/D*) REFERRED TO X,Y,Z (BNL)
C          WTI HERE IS W*TI INVERSE
C
21	CALL AMI3(WT,WTI,DET)
	READ(LI,3) ((ABCD(I,J),J=1,4),I=1,NFACES)
	WRITE (LO,22)
22	FORMAT('-CRYSTAL BOUNDARY PARAMETERS FROM EQUATIONS OF PLANE, '
	'         'AX + BY + CZ -D = 0, REFERRED TO -X,-Y,Z (SEATTLE) SYSTE
	'M'//'0FACE NO.   H          K          L       DIST(CM)    DSTAR
	'       A         B         C           PHI(O)   OMEGA(O)   CHI(O)'
	'/)
	DO 27 I=1,NFACES
	 DO 23 J=1,3
23	  TEMP(J,1)=ABCD(I,J)
	 CALL VM(TEMP(1,1),WTI,TEMP(1,2))
	 DO 24 J=1,3
24	  DATA(I,J)=TEMP(J,2)
	 DATA(I,4)=0.1*ABCD(I,4)
C
C          (DATA(I,J),J=1,3) NOW CONTAIN H*(1/D*)
C          ASSUME NO H*(1/D*) CAN BE LESS THAN .1
C
	 DO 25 J=1,3                           q
	  TEMP(J,1)=ABS(TEMP(J,2))
	  IF (TEMP(J,1).GT..1) GOTO25
	  TEMP(J,1)=100000.
	  TEMP(J,2)=0.
25	 CONTINUE
C
	 DMIN=AMIN1(TEMP(1,1),TEMP(2,1),TEMP(3,1))
	 DO 26 J=1,3
26	  TEMP(J,1)=DATA(I,J)/DMIN
C
C          LOWEST POSSIBLE VALUES OF (H,K,L) NOW IN TEMP(J,1)
C          GET TRUE DSTAR
C
	 CALL DSTAR (TEMP(1,1),TEMP(2,1),TEMP(3,1),RCP,DST)
	 OMT=0.  ;  FI=0.
	 CALL ANGLES(OMT,FI,W,DST,RCP,PHI,CHI,TEMP(1,1),TEMP(2,1),
	*     TEMP(3,1) )
	 PHI=PHI+PHIZRO
	 OMT=OMT+OMZRO
	 CHI=CHI+CHIZRO
C
	 WRITE (LO,8) I,(TEMP(J,1),J=1,3),DATA(I,4),DST,(ABCD(I,J),J=1,3),
	*     PHI,OMT,CHI
C
C          DOES NOT MATTER IF DATA(I,J)J=1,3 ARE H,K,L OR H/D*, K/D*,
C               L/D*, AS ROUTINE *GRIT2* *NORMALIZES* ANYWAY.  IF USE
C               H/D*, K/D*, L/D*, THEN NORMALIZATION LOOP IN *GRIT2* IS
C               A DO NOTHING LOOP
C               HERE TEMP(J,1) ARE LOWEST POSSIBLE (H,K,L) VALUES
C               UNNORMALIZED.   (DATA(I,J),J=1,3) ARE NORMALIZED
C
27	CONTINUE
	RETURN
	END
	SUBROUTINE MV(X,Y,Z)
C     SUBROUTINE MV MULTIPLIES VECTOR Y WITH MATRIX X. Z(I)=X(I,J)*Y(J)
	DIMENSION X(3,3),Y(3),Z(3)
	DO 1 I=1,3
	Z(I)=0.
	DO 1 J=1,3
1	Z(I)=Z(I)+X(I,J)*Y(J)
	RETURN
	END
	SUBROUTINE INSIDE(XYZ,NOUT,NFACES,DATA)
C     CHECKS WHETHER POINT XYZ IS ON SURFACE OF CRYSTAL , SETS NOUT=1 IF
	DIMENSION XYZ(3),DATA(40,4)
	NOUT=0
	DO1 J=1,NFACES
	CHECK=0.
	DO2 I=1,3
2	CHECK=CHECK+DATA(J,I)*XYZ(I)
	TROUBLE=DATA(J,4)-CHECK+.0000000001
	IF(TROUBLE)3,1,1
3	NOUT=1
	GOTO 4
1	CONTINUE
4	RETURN
	END
	SUBROUTINE DISTA(DIRECT)
	DIMENSION DIRECT(3,3),DIR(3),REC(3)
	COMMON /DST/A(3,100)
	PRINT 6
2	READ 10,KA,KB
	IF (KA.EQ.KB)RETURN
	BRA=0.
	DO 4 I=1,3
4	 DIR(I)=A(I,KA)-A(I,KB)
	CALL VM(DIR,DIRECT,REC)
	DO 5 J=1,3
5	 BRA=BRA+DIR(J)*REC(J)
	BRA=SQRT(BRA)
	PRINT 8,KA,KB,BRA
	GOTO2
6	FORMAT(//25X'DISTANCES BETWEEN POINTS OF INTERSECTION(VERTICES)'/)
8	FORMAT (26X'DISTANCE BETWEEN 'I4' AND 'I4'  = 'F10.6,1X'CM')
10	FORMAT(2I5)
	END
	SUBROUTINE MINMAX(A,LA,AMIN,AMAX)
	DIMENSION A(50)
	AMIN=A(1)
	AMAX=A(1)
	DO1 I=2,LA
	 AMIN=AMIN1(AMIN,A(I))
1	 AMAX=AMAX1(AMAX,A(I))
	RETURN
	END
	SUBROUTINE AMI2(A,B,DET)
C     INVERTS ASYMMETRIC 2X2 MATRIX
	DIMENSION A(2,2),B(2,2)
	DET=A(1,1)*A(2,2)-A(1,2)*A(2,1)
	IF(DET)1,5,1
1	B(1,1)= A(2,2)/DET
	B(1,2)=-A(1,2)/DET
	B(2,1)=-A(2,1)/DET
	B(2,2)= A(1,1)/DET
5	RETURN
	END
	SUBROUTINE SCALED (NH,K,L,ISIG,M,NSCALE,NREJ,JM,TODAY,NLINES)
C     VERSION SEPT69 FOR PICKER DIFFRACTOMETER
C     CHANGES SCALE IF WEIGHTED AVERAGE OF STANDARD REFLECTIONS CHANGES
C    1BY MORE THAN ERR
C     SUBROUTINE REQUIRES SCLA=1 AT BEGINNING OF MAIN PROGRAM
	DIMENSION RATIO(3),RSIG(3),RLAST(3)
	DATA ERR,FIRST,IREF/0.02,0.,0/
	 DATA FSC,SCAL/2*1./
	DATA N/3/
	DATA FLIM1,FLIM2/50.,0.02/
	DATA ERRL1,ERRL2/.0001,20000./
	NREJ=0
	IF(NH.GE.99000) GO TO 11
	IREF=0
	RETURN
11	NH=NH-100000
	NREJ=1
	IREF=IREF+1
	RSIG(IREF)=ISIG
	IF(IREF.LT.N) RETURN
	IF(FIRST.EQ.0.) GO TO 13
	PUM=0.; DEN=0.
	DO 12 I=1,N
	RATIO(I)=RSIG(I)/RLAST(I)
	IF(RATIO(I).LT.FLIM2.OR.RATIO(I).GT.FLIM1) GO TO 12
	PUM=PUM+RSIG(I)
	DEN=DEN+RLAST(I)
	IF(RATIO(I).LT.ERRL1.OR.RATIO(I).GT.ERRL2) GOTO15
12	CONTINUE
	IF(DEN.EQ.0.) RETURN
	RALL=PUM/DEN
	 SCAL= SCAL*RALL
	 TEST=SCAL/FSC
	 IF(TEST.LT.1.-ERR.OR.TEST.GT.1.+ERR) FSC=SCAL
	TODAY=1./FSC
	 PRINT 500,JM,RATIO(1),RATIO(2),RATIO(3),RALL,SCAL,FSC
	NLINES = NLINES + 3
500	FORMAT (/, ' DATUM NO.', I5, 5X, 'RATIO 1', F6.3, 5X, 'RATIO 2',
	1F6.3, 5X, 'RATIO 3', F6.3, 5X, 'WEIGHTED RATIO', F6.3, 5X, 'SCALE'
	2, F6.3, 5X, 'SCALE APPLIED', F6.3, /)
	IREF=0
13	FIRST=1.
	DO 14 I=1,3
14	RLAST(I)=RSIG(I)
	RETURN
15	PRINT501,I,RATIO(I)
501	FORMAT(' STANDARDS SUPPRESSED 'I5,F10.5)
	RETURN
	END
	SUBROUTINE SCALEB (JH, JK, JL, IRNET, IBKG, NSTD, IHSTD, IDIR,
	1ICOUNT, FACTOR, ILIST)
C
C     ROBERT H. BLESSING
C     DEPARTMENT OF CHEMISTRY
C     STATE UNIVERSITY OF NEW YORK AT BUFFALO
C     OCTOBER 1973
C
C     SCALING ACCORDING TO THE VARIATION OF THE STANDARD REFERENCE
C     REFLECTIONS USING THE INTERPOLATION SCHEME EMPLOYED IN THE DATA
C     REDUCTION PROGRAM, PIKOT, BY EDWARD L. MC GANDY, DEPARTMENT OF
C     CRYSTALLOGRAPHY, UNIVERSITY OF PITTSBURGH, JULY 1969.
C
	DIMENSION IH(3), IHSTD(3,30), KIND(600), MARK(600), R(600),
	1SIGR(600), R1(3), U(200), V(200)
	IF (NSTD .EQ. 0) RETURN
	GO TO (1, 2, 3, 4), IDIR
C
C     IDENTIFY AND STORE STD. REF. REFLNS.
C
1	IH(1) = JH - 100000
	IH(2) = JK
	IH(3) = JL
	IF(ILIST.EQ.0) LIST = 0
	DO 10 J = 1, NSTD
	K = J
	N = 0
	DO 11 I = 1, 3
	IF (IH(I) .EQ. IHSTD(I,J)) N = N + 1
11	CONTINUE
	IF (N .EQ. 3) GO TO 100
10	CONTINUE
	RETURN
100	LIST = LIST + 1
	KIND(LIST) = K
	MARK(LIST) = ICOUNT
	R(LIST) = IRNET
	SIGR(LIST) = SQRT(IRNET + 2.0*IBKG)
	RETURN
C
C     SET UP SCALING FUNCTION
C
C     FIND FIRST STD. REF. REFLN. OF EACH KIND
C
2	DO 725 J = 1, NSTD
	DO 715 I = 1, LIST
	IF (J .EQ. KIND(I)) GO TO 720
715	CONTINUE
720	R1(J) = R(I)
725	CONTINUE
C
C     CONVERT STD. REF. REFLNS. TO INVERSE CORRECTION FACTORS
C     INITIALIZED TO UNITY                     q
C
	DO 730 I = 1, LIST
	J = KIND(I)
	R(I) = R(I)/R1(J)
730	SIGR(I) = SIGR(I)/R1(J)
C
C     AVERAGE INV. CORR. FACTORS OVER EACH LOCAL GROUP AND FILE AVERAGED
C     VALUES
C
	M = 0
	RNUM = 0.0
	qRDEN = 0.0
	I = 3
	DO 750 J = 1, LIST
	M = M + 1
	RNUM = RNUM + (R(J)/(SIGR(J)**2))
	RDEN = RDEN + (1.0/(SIGR(J)**2))
	IF (J .EQ. LIST) GO TO 751
	IF ((MARK(J + 1) - MARK(J)) .LE. 4) GO TO 750
751	U(I) = MARK(J) - M + 1
	V(I) = RNUM/RDEN
	I = I + 1
	M = 0
	RNUM = 0.0
	RDEN = 0.0
750	CONTINUE
	N = I + 1
C
C     SET INITIAL AND TERMINAL VALUES IN INV. CORR. FACTOR FILE
C
	U(1) = U(3) - 20.0
	U(2) = U(3) - 10.0
	V(1) = V(3)
	V(2) = V(3)
	U(N - 1) = U(N - 2) + 10.0
	U(N) = U(N - 2) + 20.0
	V(N - 1) = V(N - 2)
	V(N) = V(N - 2)
	ICOUNT = 0
	X2 = 0.0
	X3 = 0.0
	I = 1
	RETURN
C
C     INTERPOLATE INV. CORR. FACTOR FILE
C
C     ENTER HERE FOR LINEAR, TWO-POINT INTERPOLATION
C
3	ICOUNT = ICOUNT + 1
	XC = ICOUNT
	IF (XC .LE. X2) GO TO 770
	I = I + 1
	X1 = U(I)
	X2 = U(I + 1)
	Y1 = V(I)
	Y2 = V(I + 1)
770	YC = Y1 + ((Y2 - Y1)/(X2 - X1))*(XC - X1)
	FACTOR = 1.0/YC
	RETURN
C
C     ENTER HERE FOR CUBIC POLYNOMIAL, FOUR-POINT INTERPOLATION
C
4	ICOUNT = ICOUNT + 1
	XC = ICOUNT
	IF (XC .LT. X3) GO TO 765
	X1 = U(I)
	X2 = U(I + 1)
	X3 = U(I + 2)
	X4 = U(I + 3)
	Y1 = V(I)
	Y2 = V(I + 1)
	Y3 = V(I + 2)
	Y4 = V(I + 3)
	I = I + 1
765	YC = 0.5*(Y1*(XC - X2)*(XC - X3)/((X1 - X2)*(X1 - X3))
	1        + Y2*(XC - X1)*(XC - X3)/((X2 - X1)*(X2 - X3))
	2        + Y3*(XC - X1)*(XC - X2)/((X3 - X1)*(X3 - X2))
	3        + Y2*(XC - X3)*(XC - X4)/((X2 - X3)*(X2 - X4))
	4        + Y3*(XC - X2)*(XC - X4)/((X3 - X2)*(X3 - X4))
	5        + Y4*(XC - X2)*(XC - X3)/((X4 - X2)*(X4 - X3)))
	FACTOR = 1.0/YC
	RETURN
	END
	SUBROUTINE DIRCOS(GRID,MDGRID,NDGRID,APHI,XN,XNSO)
C     DETERMINES COMPONENTS OF INCIDENT AND DIFFRACTED BEAMS AND COMPONE
C     OF NORMAL TO DIFFRACTION PLANE AND NORMAL TO INCIDENT BEAM IN
C     DIFFRACTION PLANE
C     SUBSIDIARY SUBROUTINE LENGTH2 FOR CALCULATION OF ABSORPTION
C     SOD AND SD ARE THE COMPONENTS OF THE INCIDENT AND DIFFRACTED RAYS
C     DIRECT SPACE
C     SO AND S ARE THE COMPONENTS OF THE INCIDENT AND DIFFRACTED RAYS IN
C     CAL SP1CE
	DIMENSION GRID(MDGRID,NDGRID)
	DIMENSION XN(3,10),XNSO(3,10),SOD(3),SD(3),TN(3)
	DIMENSION NGAUSS(3),EPOB(1000),DATA(40,4),SCB(10)
	DIMENSION X(3,10),OBSI(10),BACK(10),FI(10),EPOA(3000)
	DIMENSION RCP(6),RECIP(3,3),DIRECT(3,3),DRC(6),SHUFEL(3,3)
	DIMENSION T(3),SUM(40),TT(3),E(3),W(3,3),RT(3,3)
	DIMENSIONS(3),SO(3),R(3,3),WI(3,3)
	EQUIVALENCE(WI,SHUFEL)
	COMMON RCP,RECIP,DIRECT,DRC,SHUFEL,X,OBSI,BACK,FI,EPOA,EPOB
	COMMON/CRYO/THICK(2),TBO,TB,EMUBER
	COMMON ELAMDA,SINSQT,COST,SINT,EMU,NPOINTS,NR
	COMMON     EXT,ABS,DATA,NO,NFACES
	COMMON NGAUSS,VOL,TBAR
	COMMONSUM,W,TT,E
	COMMON /DATDIR/ ICN(16)
	EQUIVALENCE (ICN(1),IABS)
	COMMON/ANGDIR/COSCHI,SINCHI,COSOMT,SINOMT,COSPHI,SINPHI
	IF(IABS.EQ.0)EMU=0.01
	R(1,1)=COSPHI*COSOMT-SINPHI*COSCHI*SINOMT
	R(1,2) = -COSPHI*SINOMT - SINPHI*COSCHI*COSOMT
	R(2,1) =  SINPHI*COSOMT + COSPHI*COSCHI*SINOMT
	R(2,2)=-SINPHI*SINOMT+COSPHI*COSCHI*COSOMT
	R(3,1) = -SINCHI*SINOMT
	R(3,2)=-SINCHI*COSOMT
	R(1,3)=0. ; R(2,3)=0. ; R(3,3)=0.
	DO2I=1,3
	DO2J=1,3
	RT(I,J)=0.
	DO2K=1,3
2	RT(I,J)=RT(I,J)+WI(I,K)*R(K,J)
	T(1)=-COST;T(2)=-SINT       ;T(3)=0.
	TT(1)=-COST;TT(2)=-T(2);TT(3)=0.
	DO3I=1,3
	SO(I)=0.;S(I)=0.
	DO3J=1,3
	SO(I)=RT(I,J)*T(J)+SO(I)
3	S(I)=RT(I,J)*TT(J)+S(I)
	DO7 I=1,3
	SD(I)=0.  ; SOD(I)=0.
	DO7 J=1,3
	SOD(I)=SOD(I)+RECIP(I,J)*SO(J)
7	SD(I)=SD(I)+RECIP(I,J)*S(J)
	CALL NORMAL(SD,SOD,DRC,RCP,RECIP,XN,NO)
	IF(ICN(12).EQ.0) GOTO70
	TBO=0.  ; TB=0.
	DO4 I=1,3
	TBO=TBO+SOD(I)*WI(I,3)
4	TB=TB+SD(I)*WI(I,3)
	ISHIELD=ICN(12)
	TBO=THICK(ISHIELD)/ SQRT(1.-TBO**2)
	TB=THICK(ISHIELD)/ SQRT(1.-TB**2)
70	DO8 I=1,3
8	TN(I)=XN(I,NO)/DRC(I)
	CALL NORMAL(TN,SOD,DRC,RCP,RECIP,XNSO,NO)
	CALLLENGTH2(SD,SOD,GRID,MDGRID,NDGRID)
	RETURN
	END
	SUBROUTINENORMAL(SD,SOD,DRC,RCP,RECIP,XN,NO)
C     SUBROUTINE CALCULATES COMPONENTS OF NORMAL TO DIFFRACTION PLANE
C     ALONG A,B,C.
	DIMENSIONDRCS(3)
	DIMENSIONSD(3),SOD(3),RCP(6),RECIP(3,3),DRC(6),XN(3,10),YN(3)
	DATAIS/0/
	IF(IS.NE.0)GOTO10
	DO11I=1,3
11	DRCS(I)=SQRT(1.-DRC(I+3)**2)
	IS=1
10	I=1;J=2
1	K=6-I-J
C     YN ARE COMPONENTS OF NORMAL ALONG RECIPROCAL AXES
	YN(K)=(SD(I)*SOD(J)-SD(J)*SOD(I))*DRC(I)*DRC(J)*DRCS(K)/RCP(K)
	IF(I.EQ.1)GOTO2
	IF(I.EQ.3)GOTO3
	IF(I.EQ.2)GOTO4
2	I=3;J=1;GOTO1
3	I=2;J=3;GOTO1
4	DO5I=1,3
C     XN ARE COMPONENTS OF NORMAL ALONG DIRECT AXES
	XN(I,NO)=0.
	DO5J=1,3
5	XN(I,NO)=XN(I,NO)+RECIP(I,J)*YN(J)
	FNORM=0.
	DO6 I=1,3
6	FNORM=FNORM+XN(I,NO)*YN(I)
	FNORM=SQRT(FNORM)
	DO7 I=1,3
C     CONVERT TO COMPONENTS ALONG UNIT AXIS AND NORMALIZE
7	XN(I,NO)=XN(I,NO)/FNORM*DRC(I)
	RETURN
	END
	SUBROUTINE LENGTH2(SD,SOD,GRID,MDGRID,NDGRID)
	DIMENSION GRID(MDGRID,NDGRID)
C     CALCULATES PATH LENGTH OF RAY IN CRYSTAL FOR EACH OF THE SAMPLING
C     OF THE GAUSSIAN GRID AND THE ABSORPTION CORRECTION
	DIMENSION RCP(6),RECIP(3,3),DIRECT(3,3),DRC(6),SHUFEL(3,3)
	DIMENSION X(3,10),OBSI(10),BACK(10),FI(10),EPOA(3000)
	DIMENSION NGAUSS(3),EPOB(1000),DATA(40,4),SCB(10)
	DIMENSION SD(3),SOD(3),G(50),H(50),COSIN(50),COSOUT(50)
	DIMENSION NIN(26),NOUT(26)
	DIMENSION SUM(40),TT(3),E(3),W(3,3)
	COMMON RCP,RECIP,DIRECT,DRC,SHUFEL,X,OBSI,BACK,FI,EPOA,EPOB
	COMMON ELAMDA,SINSQT,COST,SINT,EMU,NPOINTS,NR
	COMMON     EXT,ABS,DATA,NO,NFACES
	COMMON NGAUSS,VOL,TBAR
	COMMONSUM,W,TT,E
	MO=1
	LO=1
	TTABS=0.
C     CALCULATE COSINES OF ANGLES OF INCIDENT AND DIFFRACTED RAYS WITH A
	DO 310 K=1,NFACES
	COSIN(MO)=0.
	COSOUT(LO)=0.
	DO 311 I=1,3
	COSIN(MO)=COSIN(MO)+SOD(I)*DATA(K,I)
311	COSOUT(LO)=COSOUT(LO)+SD(I)*DATA(K,I)
C     REJECT FACES FOR WHICH ANGLE BETWEEN INCIDENT BEAM AND FACENORMAL
	IF(COSIN(MO))312,313,313
312	NIN(MO)=K
	MO=MO+1
C     REJECT FACE IF ANGLE BETWEEN DIFFRACTED RAY AND FACE NORMAL IS OBT
313	IF(COSOUT(LO))310,310,314
314	NOUT(LO)=K
	LO=LO+1
310	CONTINUE
	LO=LO-1
	MO=MO-1
C     CALCULATE PATH LENGTH AND ABSORPTION FOR EACH OF THE SAMPLING POIN
	DO 12 LL=1,NPOINTS
	DO 9 K=1,MO
	RDOTN=0.
	NMO=NIN(K)
	DO 4 I=1,3
4	RDOTN=RDOTN+GRID(I,LL)*DATA(NMO,I)
	RCOS=DATA(NMO,4)-RDOTN
	 G(K)=-RCOS/COSIN(K)
9	CONTINUE
	DO 91 K=1,LO
	RDOTN=0.
	NLO=NOUT(K)
	DO 41 I=1,3
41	RDOTN=RDOTN+GRID(I,LL)*DATA(NLO,I)
	RCOS=DATA(NLO,4)-RDOTN
	H(K)=RCOS/COSOUT(K)
91	CONTINUE
	GK=G(1)
	DO 10 K=1,MO
C     CHOOSE MINIMUM PATH AS TRUE PATH LENGTH
10	 GK=AMIN1(GK,G(K))
	HK=H(1)
	DO 11 K=1,LO
C     USING TABLE LOOK-UP CALCULATE EXPONENTIAL FUNCTION
11	 HK=AMIN1(HK,H(K))
	GKHK=GK+HK
	EUR=GKHK+EMU*10. + 0.0005
	EUR=(GK+HK)*EMU*10.+0.0005
	MUR=EUR
	DUR=MUR
	LUR=(EUR-DUR)*1000.
	A=EPOA(MUR+1)
	B=EPOB(LUR+1)
C     MULTIPLY ABSORPTION FOR EACH GRID POINT BY ITS GAUSSIAN WEIGHT
	AB=A*B*GRID(4,LL)
	TAB=GKHK*AB
	TTABS=TTABS+TAB
C     SUM ABSORPTION OVER ALL SAMPLING POINTS
	ABS=ABS+AB
12	CONTINUE
	TBAR=TTABS/ABS
	ABS=ABS/VOL
	RETURN
	END
	SUBROUTINE IMPD(FH,FK,FL,S,B,SCB,OM,PHI,SLO,SLB,OBSFI,BACFK,
	&IDENT,NSCALE,ICN,NBACK,NSER)
	DATA N/1/
	DATA NN /0/
	READ(105,1)IH,FH,IK,FK,IL,FL,B1,S,B2,SLO
1	FORMAT(3(I1,F2.0),3F7.0,F5.1)
	NN = NN + 1
	NSER = NN
	IF(IH.EQ.9)GO TO 90
	IF(IH.EQ.8)N=N+1;IDENT=1;GO TO 80
	IF(IH.EQ.1)FH=-FH
	IF(IK.EQ.1)FK=-FK
	IF(IL.EQ.1)FL=-FL
	IF(FH.EQ.0.0.AND.FK.EQ.4.0.AND.FL.EQ.0.0)FH=FH+100000
	IF(FH.EQ.0.0.AND.FK.EQ.0.0.AND.FL.EQ.4.0)FH=FH+100000
	S = S * 10 + 5
	B = B1*10+5.0 +B2*10+5.0
	SCB=1.0
	OBSFI=S
	BACFK = B
	IDENT=0
	NSCALE=N
	NBACK=2
	SLB=40.0
	OM = 0.0
	PHI = 0.0
	IF(N.EQ.3)SLB=20.0
80	RETURN
90	IDENT=2
	RETURN
	END
