C.....
C  ..           PROGRAM ABSORB1
C  .. SET UP THE ABSORB.DAT DATA FILE FOR ABSORPTION CORRECTION
C  ..
C  .. GEORGE T. DETITTA, MEDICAL FOUNDATION OF BUFFALO
C  .. FEBRUARY 10, 1982
C.....
      PROGRAM ABSORB1
C-
      CHARACTER TITLE*80, ATIME*8, ADATE*9
      CHARACTER FIL01*80, FIL02*80
      DIMENSION CELL(6)
      CHARACTER*4 WAVE
      CHARACTER*2 ATOMS(10)
      DIMENSION STOIC(10)
      DIMENSION FACES(50,4)
      DIMENSION NGP(3), D1(3),D2(3),D3(3)
      DIMENSION ORMAT(3,3)
      DIMENSION NTRP(50)
      COMMON /INPT01/ TITLE, ATIME, ADATE
      COMMON /INPT02/ FIL01, FIL02, IPATH
      COMMON /INPT03/ CELL
      COMMON /INPT04/ WAVE
      COMMON /LAMBDA/ FLAM
      COMMON /INPT05/ ATOMS
      COMMON /INPT06/ ZCELL, STOIC
      COMMON /INPT07/ FACES
      COMMON /INPT08/ NGP, D1,D2,D3
      COMMON /INPT09/ ORMAT, NDIFF
      COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP
      COMMON /INPT11/ NBPR, A0, A1, A2, A3
      COMMON /INPT12/ NSOL, EMUS, NTRP
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      DIMENSION UBF(3,3), UBI(3,3), UBGF(3,3), UBGI(3,3)
      DIMENSION FS(50,5)
      DIMENSION VR(150,3), MM(150,3)
      PARAMETER (NGMAX=32**3)
      DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3)
      DIMENSION UPTSF(3,3), UPTSI(3,3), UBIG(3,3), DTUB(NGMAX,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
      COMMON /CRYSTL/ TMOLW, RHO, EMUX
      COMMON /MATRIX/ UBF, UBI, UBGF, UBGI
      COMMON /PLANES/ NF, FS
      COMMON /VERTEX/ NV, VR, MM
      COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC
      COMMON /KAPCOR/ EMUK, RAD2, RPT2, UPTSF, UPTSI, UBIG, DTUB
C-
      NIN=IO5
      NOU=IO6
      CALL SETUP
      STOP 'ABSORB1 DONE'
      END
C.....
C  ..           SUBROUTINE BLGEN
C.....
C  .. GIVEN H,K,L PSI AND THE ORIENTATION MATRIX
C  .. CALCULATE 2THETA OMEGA PHI AND CHI FOR THE
C  .. BUSING-LEVY DIFFRACTOMETER
C.....
      SUBROUTINE BLGEN (FH,FK,FL,FPSI, ANGLS)
      DIMENSION UBF(3,3), UBI(3,3), UBGF(3,3), UBGI(3,3)
      COMMON /MATRIX/ UBF, UBI, UBGF, UBGI
      COMMON /LAMBDA/ FLAM
      DIMENSION ANGLS(4)
      DATA CONS/57.2957795/
C-
      TT = 0.
      OM = 0.
      PH = 0.
      CH = 0.
C-    GET CARTESIAN COMPONENTS OF RECIP VECTOR
      XI = FH*UBF(1,1) + FK*UBF(1,2) + FL*UBF(1,3)
      YI = FH*UBF(2,1) + FK*UBF(2,2) + FL*UBF(2,3)
      ZI = FH*UBF(3,1) + FK*UBF(3,2) + FL*UBF(3,3)
C-    CALCULATE DSTAR
      Q = SQRT(XI*XI + YI*YI + ZI*ZI)
C-    TEST FOR THE 000 REFLECTION
      IF (Q .LE. 0.00001) GO TO 999
C-    GET SIN THETA
      SINT = 0.5 * FLAM * Q
      IF (SINT .GE. 1.0) GO TO 999
C-    GET TWO THETA FROM SIN AND COS THETA
      COST = SQRT(1.0 - SINT*SINT)
      TT = 2.0*CONS*ATAN2(SINT,COST)
C-    GET OTHER ANGLES
      IF(XI.EQ.0.0 .AND. YI.EQ.0.0) THEN
      PHI = 0.0
      ELSE
      PHI = ATAN2(YI,XI)
      ENDIF
      CHI = ATAN2(ZI,SQRT(XI*XI+YI*YI))
      PH = AMOD(PHI*CONS+360., 360.)
      CH = AMOD(CHI*CONS+360., 360.)
  999 ANGLS(1) = TT
      ANGLS(2) = OM
      ANGLS(3) = PH
      ANGLS(4) = CH
      RETURN
      END
C.....
C  ..           SUBROUTINE BPRO1
C.....
C  .. IF THE DATA WERE MEASURED WITH A MONOCHROMATOR
C  .. ENTER THE ANALYTICAL FORM OF THE BEAM PROFILE
C  .. EQUATION OF THE FORM
C  ..  F = A0 + A1*DEL +A2*DEL*DEL +A3*DEL*DEL*DEL
C.....
      SUBROUTINE BPRO1
      COMMON /INPT11/ NBPR, A0, A1, A2, A3
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      A0=1.0
      A1=0.0
      A2=0.0
      A3=0.0
      WRITE(IO6,9001)
 9001 FORMAT(1X,'ENTER BEAM PROFILE FLAG:',/,
     *       1X,'NO MONOCHROMATOR......NBPR=0',/,
     *       1X,'SCALE FACTOR..........NBPR=1',/,
     *       1X,'LINEAR FIT............NBPR=2',/,
     *       1X,'QUADRATIC FIT.........NBPR=3',/,
     *       1X,'CUBIC FIT.............NBPR=4'   )
      WRITE(IO6,9002)
 9002 FORMAT(1X,'ENTER NBPR: ',$)
      READ(IO5,9003) NBPR
 9003 FORMAT(I5)
      IF(NBPR.LE.0) GO TO 100
      IF(NBPR.EQ.1) GO TO 200
      IF(NBPR.EQ.2) GO TO 300
      IF(NBPR.EQ.3) GO TO 400
      IF(NBPR.GE.4) GO TO 500
  100 NBPR = 0
      WRITE(IO6,9004)
 9004 FORMAT(1X,'NO MONOCHROMATOR...NO PROFILE')
      RETURN
  200 WRITE(IO6,9005)
 9005 FORMAT(1X,'SCALE FIT')
      WRITE(IO6,9006)
 9006 FORMAT(1X,'ENTER A0:  ',$)
      READ(IO5,9007) A0
 9007 FORMAT(4F10.6)
      RETURN
  300 WRITE(IO6,9008)
 9008 FORMAT(1X,'LINEAR FIT')
      WRITE(IO6,9009)
 9009 FORMAT(1X,'ENTER A0, A1:  ',$)
      READ(IO5,9007) A0,A1
      RETURN
  400 WRITE(IO6,9010)
 9010 FORMAT(1X,'QUADRATIC FIT')
      WRITE(IO6,9011)
 9011 FORMAT(1X,'ENTER A0, A1, A2:  ',$)
      READ(IO5,9007) A0, A1, A2
      RETURN
  500 WRITE(IO6,9012)
 9012 FORMAT(1X,'CUBIC FIT')
      WRITE(IO6,9013)
 9013 FORMAT(1X,'ENTER A0, A1, A2, A3:  ',$)
      READ(IO5,9007) A0, A1, A2, A3
      RETURN
      END
C.....
C  ..           SUBROUTINE C4GEN
C.....
C  .. CALCULATE THE CAD-4 DIFFRACTOMETER ANGLES
C  .. USING THE CODE RIPPED OUT OF THE ROUTINE
C  .. HTOW.  INPUT IS H,K,L AND PSI, OUTPUT IS
C  .. THETA, PHIK, OMGK, AND KAPPA IN DEGREES
C.....
      SUBROUTINE C4GEN (FH,FK,FL,FPSI, ANGLS)
      DIMENSION UBF(3,3), UBI(3,3), UBGF(3,3), UBGI(3,3)
      COMMON /MATRIX/ UBF, UBI, UBGF, UBGI
      COMMON /LAMBDA/ FLAM
      DIMENSION ANGLS(4)
      DATA PI/3.14159263694/
      DATA ALPHA/49.96681/
C
      CON1 = SIND(ALPHA)*SIND(ALPHA)
      CON2 = COSD(ALPHA)
      A89 = 89.
      CON3 = COSD(A89)*COSD(A89)*CON1
C
      DO 10 J=1,4
      ANGLS(J) = 0.0
   10 CONTINUE
C
      FC1 = FH*UBF(1,1) + FK*UBF(1,2) + FL*UBF(1,3)
      FC2 = FH*UBF(2,1) + FK*UBF(2,2) + FL*UBF(2,3)
      FC3 = FH*UBF(3,1) + FK*UBF(3,2) + FL*UBF(3,3)
C
      FPHIB = 180.*ATAN2(-FC1,FC2)/PI
      CO = SQRT(FC1*FC1+FC2*FC2)
      FCHIB = 180.*ATAN2(FC3,CO)/PI
      SI = SQRT(CO*CO+FC3*FC3)*FLAM/2.
      IF(SI.GE.1.)     GO TO 999
      CO = SQRT(1.-SI*SI)
      FTHETA = 180.*ATAN2(SI,CO)/PI
C
      SI = SIND(FPSI)*COSD(FCHIB)
      CO = SIND(FCHIB)
      FOME = 180.*ATAN2(SI,CO)/PI + FTHETA
      SI = -SIND(FPSI)
      CO = SIND(FCHIB)*COSD(FPSI)
      FPHIE = 180.*ATAN2(SI,CO)/PI + FPHIB
      SI = SQRT(SIND(FCHIB)**2.+(SIND(FPSI)*COSD(FCHIB))**2.)
      CO = COSD(FCHIB)*COSD(FPSI)
      FCHIE = 180.*ATAN2(SI,CO)/PI
      IF(FCHIB) 310,330,320
  310 FCHIE = -FCHIE
      FOME = FOME+180.
      FPHIE = FPHIE + 180.
  320 FPHIE = BETW(FPHIE,180.)
  325 IF(FTHETA.LT.0.0)FOME=FOME+180.
      FOME = BETW(FOME,180.)
      GO TO 350
  330 IF(SIND(FPSI)) 310,340,320
  340 FCHIE = 0.
      FPHIE = FPHIB
      FOME = FTHETA
      GO TO 325
  350 CONTINUE
C
  400 FCHIE = BETW(FCHIE,180.)
      SI = SIND(0.5*FCHIE)
      CO = CON1-SI*SI
      IF((CON3-CO).GT.0.)  GO TO 410
      CO = SQRT(CO)
      FKAPPA=360.*ATAN2(SI,CO)/PI
      SI = SI*CON2
      FOMK=180.*ATAN2(SI,CO)/PI
  420 FPHIK = FPHIE-FOMK
      FOMK = FOME-FOMK
      FPHIK = BETW(FPHIK,180.)
      FOMK = BETW(FOMK,180.)
      GO TO 475
  410 FKAPPA = 180.
      FOMK = 90.
      GO TO 420
  475 CONTINUE
C
      ANGLS(1) = FTHETA
      ANGLS(2) = FPHIK
      ANGLS(3) = FOMK
      ANGLS(4) = FKAPPA
 999  RETURN
      END
C.....
C  ..           FUNCTION BETW
C.....
      FUNCTION BETW(A,R)
      BETW=A-2.*R*FLOAT(INT((SIGN(R,A)+A)/(2.*R)))
      END
C.....
C  ..           FUNCTION SIND
C.....
      FUNCTION SIND(A)
      SIND = SIN(A/57.2957795)
      END
C.....
C  ..           FUNCTION COSD
C.....
      FUNCTION COSD(A)
      COSD = COS(A/57.2957795)
      END
C.....
C  ..           SUBROUTINE CELL1
C.....
C  .. INPUT THE UNIT CELL CONSTANTS.
C  .. THE CELL EDGES ARE IN ANGSTROMS AND THE ANGLES
C  .. ARE IN DEGREES.
C.....
      SUBROUTINE CELL1
      DIMENSION CELL(6)
      COMMON /INPT03/ CELL
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      WRITE(IO6,101)
  101 FORMAT(/1X,'ENTER CELL CONSTANTS: ',$)
      READ(IO5,102) CELL
  102 FORMAT(6F10.4)
      RETURN
      END
C.....
C  ..           SUBROUTINE DIFF1
C.....
C  .. INPUT THE DIFFRACTOMETER TYPE
C  ..    DIFFRACTOMETER        NDIFF
C  ..     NICOLET P3/F           1
C  ..     NONIUS CAD-4           2
C  ..     BUSING-LEVY            3
C.....
      SUBROUTINE DIFF1
      DIMENSION ORMAT(3,3)
      COMMON /INPT09/ ORMAT, NDIFF
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
 1000 WRITE(IO6,9001)
 9001 FORMAT(1X,'ENTER THE DIFFRACTOMETER TYPE')
      WRITE(IO6,9002)
 9002 FORMAT(1X,'NICOLET P3/F   NDIFF=1',/,
     *       1X,'NONIUS CAD-4   NDIFF=2',/,
     *       1X,'BUSING-LEVY    NDIFF=3'   )
      WRITE(IO6,9003)
 9003 FORMAT(1X,'ENTER NDIFF: ',$)
      READ(IO5,9004) NDIFF
 9004 FORMAT(I1)
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE EDIT1
C.....
C  .. MAIN EDITING ROUTINE SCHEDULER
C.....
      SUBROUTINE EDIT1
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      CHARACTER*1 TIT, FIL, WAV, STO, PLN, UNC, GRD, LIQ
      CHARACTER*1 ALL, DIF, XXX, KAP, BPF, MAT, RSP
      DATA TIT/'T'/, FIL/'F'/, WAV/'W'/, STO/'S'/
      DATA PLN/'P'/, UNC/'U'/, GRD/'G'/, MAT/'M'/
      DATA ALL/'A'/, DIF/'D'/, XXX/'X'/, KAP/'K'/
      DATA BPF/'B'/, LIQ/'L'/
C-
 1000 WRITE(IO6,9001)
 9001 FORMAT(/1X,'OPTIONS FOR EDITING ARE:',/,
     *        1X,'ALL NEW INPUT.........(A)',/,
     *        1X,'TITLE.................(T)',/,
     *        1X,'FILE NAMES............(F)',/,
     *        1X,'WAVELENGTH............(W)',/,
     *        1X,'STOICHIOMETRY.........(S)',/,
     *        1X,'PLANES................(P)',/,
     *        1X,'UNITCELL..............(U)',/,
     *        1X,'GRIDWORK..............(G)',/,
     *        1X,'DIFFRACTOMETER........(D)',/,
     *        1X,'MATRIX................(M)',/,
     *        1X,'KAPPILARY.............(K)',/,
     *        1X,'BEAM PROFILE..........(B)',/,
     *        1X,'LIQUID TRAPPED........(L)',/,
     *        1X,'EXIT FROM EDITING.....(X)'   )
      WRITE(IO6,9002)
 9002 FORMAT(1X,'ENTER EDIT OPTION: ',$)
      READ(IO5,9003) RSP
 9003 FORMAT(A)
      IF(RSP.EQ.ALL) GO TO 1001
      IF(RSP.EQ.TIT) GO TO 1002
      IF(RSP.EQ.FIL) GO TO 1003
      IF(RSP.EQ.WAV) GO TO 1004
      IF(RSP.EQ.STO) GO TO 1005
      IF(RSP.EQ.PLN) GO TO 1006
      IF(RSP.EQ.UNC) GO TO 1007
      IF(RSP.EQ.GRD) GO TO 1008
      IF(RSP.EQ.DIF) GO TO 1009
      IF(RSP.EQ.MAT) GO TO 1010
      IF(RSP.EQ.KAP) GO TO 1011
      IF(RSP.EQ.BPF) GO TO 1012
      IF(RSP.EQ.LIQ) GO TO 1013
      IF(RSP.EQ.XXX) GO TO 9999
      GO TO 1000
C-
 1001 CALL TITL1
      CALL FILS1
      CALL CELL1
      CALL FACE1
      CALL WAVE1
      CALL EMUX1
      CALL GRID1
      CALL DIFF1
      CALL MATX1
      CALL KAPP1
      CALL BPRO1
      CALL SOLN0
      CALL PAGE2
      CALL TITL2
      CALL FILS2
      CALL BPRO2
      CALL CELL2
      CALL CELL3
      CALL FACE2
      CALL FACE3
      CALL EMUX2
      CALL EMUX3
      CALL GRID2
      CALL GRID3
      CALL MATX2
      CALL MATX3
      CALL KAPP2
      CALL KAPP3
      CALL SOLN4
      GO TO 1000
 1002 CALL TITL1
      GO TO 1000
 1003 CALL FILS1
      GO TO 1000
 1004 CALL WAVE1
      CALL EMUX2
      CALL EMUX3
      GO TO 1000
 1005 CALL EMUX1
      CALL EMUX2
      CALL EMUX3
      GO TO 1000
 1006 CALL FACE1
      CALL FACE2
      CALL FACE3
      CALL GRID2
      CALL GRID3
      CALL KAPP2
      CALL KAPP3
      GO TO 1000
 1007 CALL CELL1
      CALL CELL2
      CALL FACE2
      CALL FACE3
      CALL GRID2
      CALL GRID3
      CALL KAPP2
      CALL KAPP3
      GO TO 1000
 1008 CALL GRID1
      CALL GRID2
      CALL GRID3
      CALL KAPP2
      CALL KAPP3
      GO TO 1000
 1009 CALL DIFF1
      GO TO 1000
 1010 CALL MATX1
      CALL MATX2
      CALL MATX3
      GO TO 1000
 1011 CALL KAPP1
      CALL KAPP2
      CALL KAPP3
      GO TO 1000
 1012 CALL BPRO1
      CALL BPRO2
      GO TO 1000
 1013 CALL SOLN0
      CALL SOLN4
      GO TO 1000
 9999 RETURN
      END
C.....
C  ..           SUBROUTINE EMUX1
C.....
C  .. INPUT THE COMPOSITION OF THE CRYSTAL CHEMICAL UNIT
C  .. AND THE NUMBER OF C.C.U PER UNIT CELL.
C.....
      SUBROUTINE EMUX1
      CHARACTER*2 ATOMS(10)
      COMMON /INPT05/ ATOMS
      DIMENSION STOIC(10)
      COMMON /INPT06/ ZCELL,STOIC
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      CHARACTER*2 SYMBL,BLANK,ATLAB(96)
      DATA BLANK/'  '/
      DATA ATLAB/'H ' ,'HE' ,'LI' ,'BE' ,'B ' ,'C ' ,
     *           'N ' ,'O ' ,'F ' ,'NE' ,'NA' ,'MG' ,
     *           'AL' ,'SI' ,'P ' ,'S ' ,'CL' ,'AR' ,
     *           'K ' ,'CA' ,'SC' ,'TI' ,'V ' ,'CR' ,
     *           'MN' ,'FE' ,'CO' ,'NI' ,'CU' ,'ZN' ,
     *           'GA' ,'GE' ,'AS' ,'SE' ,'BR' ,'KR' ,
     *           'RB' ,'SR' ,'Y ' ,'ZR' ,'NB' ,'MD' ,
     *           'TC' ,'RU' ,'RH' ,'PD' ,'AG' ,'CD' ,
     *           'IN' ,'SN' ,'SB' ,'TE' ,'I ' ,'XE' ,
     *           'CS' ,'BA' ,'LA' ,'CE' ,'PR' ,'ND' ,
     *           'PM' ,'SM' ,'EU' ,'GD' ,'TB' ,'DY' ,
     *           'HO' ,'ER' ,'TM' ,'YB' ,'LU' ,'HF' ,
     *           'TA' ,'W ' ,'RE' ,'OS' ,'IR' ,'PT' ,
     *           'AU' ,'HG' ,'TL' ,'PB' ,'BI' ,'  ' ,
     *           '  ' ,'RA' ,'  ' ,'  ' ,'  ' ,'TH' ,
     *           '  ' ,'U ' ,'  ' ,'PU' ,'  ' ,'  ' / 
C-
      DO 50 I=1,10
      ATOMS(I) = BLANK
      STOIC(I) = 0.0
   50 CONTINUE
      NTYPE = 0
      NMAX = 10
      ZCELL = 0.
C-
      WRITE(IO6,99)
   99 FORMAT(/1X,'ENTER THE NUMBER OF C.C.U PER UNIT CELL: ',$)
      READ(IO5,100) ZCELL
  100 FORMAT(F5.0)
C-
  500 WRITE(IO6,101)
  101 FORMAT(1X,'ENTER CHEMICAL SYMBOL: ',$)
      READ(IO5,102) SYMBL
  102 FORMAT(A)
      IF( SYMBL.EQ.BLANK ) RETURN
      NOK=0
      DO 499 J=1,96
      IF(SYMBL.EQ.ATLAB(J)) NOK=1
  499 CONTINUE
      IF(NOK.NE.1) GO TO 500
      NTYPE = NTYPE + 1
      WRITE(IO6,103)
  103 FORMAT(1X,'ENTER NUMBER OF ATOMS: ',$)
      READ(IO5,104) FATMS
  104 FORMAT(F5.0)
      ATOMS(NTYPE) = SYMBL
      STOIC(NTYPE) = FATMS
      NMAX = NMAX - 1
      IF(NMAX.GT.0) GO TO 500
      RETURN
      END
C.....
C  ..           SUBROUTINE FACE1
C.....
C  .. INPUT THE MILLERIAN SYMBOL AND  PERPENDICULAR DISTANCE
C  .. FOR EACH FACE OF THE CRYSTAL.  PERPENDICULAR DISTANCE
C  .. IN MILLIMETERS PLEASE.
C.....
      SUBROUTINE FACE1
      DIMENSION FACES(50,4)
      COMMON /INPT07/ FACES
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C.....
      WRITE(IO6,9000)
 9000 FORMAT(/1X,'ENTER FACE INFORMATION')
 1000 WRITE(IO6,9001) 
 9001 FORMAT(1X,'ENTER THE FACE NUMBER : ',$)
      READ(IO5,9002) JF
 9002 FORMAT(I2)
      IF (JF.EQ.0) RETURN
      IF (JF.GT.50) RETURN
      WRITE(IO6,9004)
 9004 FORMAT(1X,'ENTER MILLERIAN SYMBOL: ',$)
      READ(IO5,9005) FH, FK, FL
 9005 FORMAT(3F10.0)
      WRITE(IO6,9006)
 9006 FORMAT(1X,'ENTER PERP. DISTANCE(MM): ',$)
      READ(IO5,9007) PD
 9007 FORMAT(F8.5)
      FACES(JF,1) = FH
      FACES(JF,2) = FK
      FACES(JF,3) = FL
      FACES(JF,4) = PD
      GO TO 1000
      END
C.....
C  ..           SUBROUTINE FILS1
C.....
C  .. GET THE INPUT AND OUTPUT FILE NAMES
C.....
      SUBROUTINE FILS1
      CHARACTER FIL01*80, FIL02*80
      COMMON /INPT02/ FIL01, FIL02, IPATH
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C....
      WRITE(IO6,9001)
 9001 FORMAT(1X,'INPUT FILE NAME (DEV:[DIRECTORY]NAME.SCL) : ',$)
      READ(IO5,9002) FIL01
 9002 FORMAT(A)
      WRITE(IO6,9003)
 9003 FORMAT(1X,'OUTPUT FILE NAME (DEV:[DIRECTORY]NAME.ABS) :' ,$)
      READ(IO5,9002) FIL02
    5 WRITE(IO6,9004)
 9004 FORMAT(1X,'WRITE TBAR AND S0 AND S1 COMPONENTS TO OUTPUT FILE?',
     *'  Y OR N.  ',$)
      READ(IO5,9002) ANS
      IF (ANS.EQ.'Y') IPATH=1
      IF (ANS.EQ.'N') IPATH=0
      IF (ANS.NE.'Y'.AND.ANS.NE.'N') GO TO 5
      RETURN
      END
C.....
C  ..           SUBROUTINE GEOM1
C.....
C  .. ASSUMING AN INTACT AND COMPLETE SET OF PRIMARY DATA
C  .. GENERATE A SUMMARY OF CRYSTAL GEOMETRY AND GENERAL INFO
C.....
      SUBROUTINE GEOM1
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      CHARACTER*1 CEL, PLN, VER, LIM, EDG, GRD, MAT
      CHARACTER*1 STO, REC, KAP, ALL, XXX, RSP
      DATA CEL/'C'/, PLN/'P'/, VER/'V'/, LIM/'L'/
      DATA EDG/'E'/, GRD/'G'/, MAT/'M'/, ALL/'A'/
      DATA STO/'S'/, REC/'R'/, XXX/'X'/, KAP/'K'/
C-
  500 WRITE(IO6,9001)
 9001 FORMAT(/1X,'GEOMETRY AND INFORMATION SUMMARY',/,
     *        1X,'OPTIONS INCLUDE: ',/,
     *        1X,'UNIT CELL DATA.................(C)',/,
     *        1X,'CRYSTAL PLANE DATA.............(P)',/,
     *        1X,'VERTEX GENERATION..............(V)',/,
     *        1X,'LIMITS OF CRYSTAL..............(L)',/,
     *        1X,'EDGES OF THE CRYSTAL...........(E)',/,
     *        1X,'GRID SIZE INFORMATION..........(G)',/,
     *        1X,'ORIENTATION MATRIX DATA........(M)',/,
     *        1X,'STOICHIOMETRY DATA.............(S)',/,
     *        1X,'KAPPILARY DATA.................(K)',/,
     *        1X,'ALL OF THE ABOVE...............(A)',/,
     *        1X,'RECAP THE OPTIONS..............(R)',/,
     *        1X,'EXIT FROM GEOM.................(X)'   )
 1000 WRITE(IO6,9002)
 9002 FORMAT(1X,'ENTER OPTION (R TO RECAP): ',$)
      READ(IO5,9003) RSP
 9003 FORMAT(A)
      IF(RSP.EQ.CEL) GO TO 1001
      IF(RSP.EQ.PLN) GO TO 1002
      IF(RSP.EQ.VER) GO TO 1003
      IF(RSP.EQ.LIM) GO TO 1004
      IF(RSP.EQ.EDG) GO TO 1005
      IF(RSP.EQ.GRD) GO TO 1006
      IF(RSP.EQ.MAT) GO TO 1007
      IF(RSP.EQ.STO) GO TO 1008
      IF(RSP.EQ.KAP) GO TO 1009
      IF(RSP.EQ.ALL) GO TO 1010
      IF(RSP.EQ.REC) GO TO 500
      IF(RSP.EQ.XXX) RETURN
C-
 1001 CALL CELL2
      CALL CELL3
      GO TO 1000
C-
 1002 CALL FACE2
      CALL FACE3
      GO TO 1000
C-
 1003 CALL VERT1
      CALL VERT5
      GO TO 1000
C-
 1004 CALL LIMS1
      GO TO 1000
C-
 1005 CALL EDGE1
      GO TO 1000
C-
 1006 CALL GRID2
      CALL GRID3
      GO TO 1000
C-
 1007 CALL MATX2
      CALL MATX3
      GO TO 1000
C-
 1008 CALL EMUX2
      CALL EMUX3
      GO TO 1000
C-
 1009 CALL KAPP2
      CALL KAPP3
      GO TO 1000
C-
 1010 CALL TITL2
      CALL FILS2
      CALL CELL2
      CALL CELL3
      CALL FACE2
      CALL FACE3
      CALL VERT1
      CALL VERT5
      CALL LIMS1
      CALL EDGE1
      CALL MATX2
      CALL MATX3
      CALL EMUX2
      CALL EMUX3
      CALL GRID2
      CALL GRID3
      CALL KAPP2
      CALL KAPP3
      RETURN
      END
C.....
C  ..           SUBROUTINE GRID1
C.....
C  .. INPUT THE RECIPROCAL DIRECTIONS AND NUMBER OF
C  .. GRID POINTS WHICH WILL BE USED FOR THE GAUSSIAN
C  .. SECTIONING OF THE CRYSTAL.  THE THREE RECIPROCAL
C  .. VECTORS MUST NOT BE COPLANAR OR COLINEAR AND IT
C  .. IS ADVANTAGEOUS THAT THEY BE NEARLY ORTHOGONAL.
C  .. THE NUMBER OF GRID POINTS ALONG EACH DIRECTION
C  .. MUST CORRESPOND TO ONE OF THE TABULATED VALUES
C  .. OF GAUSSIAN WEIGHTS AND FRACTIONS FROM THE ROUTINE
C  .. GAUTAB.
C.....
      SUBROUTINE GRID1
      DIMENSION NGP(3), D1(3),D2(3),D3(3)
      COMMON /INPT08/ NGP, D1,D2,D3
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      DIMENSION NOK(19),NP(3), U1(3),U2(3),U3(3)
      DATA NOK/2,3,4,5,6,7,8,9,10,12,16,20,24,32,40,48,64,80,96/
      PARAMETER (NGMAX=32**3)
C-
      WRITE(IO6,9000)
 9000 FORMAT(/1X,'GRID GENERATION DATA')
 1000 WRITE(IO6,9001)
 9001 FORMAT(1X,'ENTER RECIPROCAL VECTOR (HKL) #1: ',$)
      READ(IO5,9002) D1
 9002 FORMAT(3F5.0)
      WRITE(IO6,9003)
 9003 FORMAT(1X,'ENTER RECIPROCAL VECTOR (HKL) #2: ',$)
      READ(IO5,9002) D2
      WRITE(IO6,9004)
 9004 FORMAT(1X,'ENTER RECIPROCAL VECTOR (HKL) #3: ',$)
      READ(IO5,9002) D3
 1001 WRITE(IO6,8005)
 8005 FORMAT(/1X,'ENTER NUMBER OF GRID POINTS ALONG VECTOR #1: ',$)
      READ(IO5,8006) NP(1)
 8006 FORMAT(I2)
      WRITE(IO6,8007)
 8007 FORMAT(1X,'ENTER NUMBER OF GRID POINTS ALONG VECTOR #2: ',$)
      READ(IO5,8006) NP(2)
      WRITE(IO6,8008)
 8008 FORMAT(1X,'ENTER NUMBER OF GRID POINTS ALONG VECTOR #3: ',$)
      READ(IO5,8006) NP(3)
C-    CHECK TO SEE THAT THESE ARE OK
      DO 10 I=1,3
      IF (NP(I) .LE.  2)  NP(I) = 2
      IF (NP(I) .GE. 96)  NP(I) = 96
      DO 5 J=1,18
      K = J + 1
      NBOT = NOK(J)
      NTOP = NOK(K)
      IF (NP(I).GE.NBOT .AND. NP(I).LT.NTOP) NP(I) = NBOT
    5 CONTINUE
   10 CONTINUE
      DO 20 I=1,3
      NGP(I) = NP(I)
   20 CONTINUE
      NGPTS = NGP(1)*NGP(2)*NGP(3)
      IF (NGPTS.LE.NGMAX) RETURN
      WRITE(IO6,8011) NGPTS, NGMAX
 8011 FORMAT(1X,'NUMBER OF POINTS: ',I6,2X,' EXCEEDS MAX OF: ',I6)
      GO TO 1001
      END
C.....
C  ..           SUBROUTINE KAPP1
C.....
C  .. INPUT THE DATA FOR CRYSTALS IN CAPPILARIES
C.....
      SUBROUTINE KAPP1
      COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      NKAP=0
      DKAP=0.0
      TKAP=0.0
      TAU=0.0
      SIG=0.0
      XDSP=0.0
      YDSP=0.0      
      WRITE(IO6,9001)
 9001 FORMAT(/1X,'CAPILLARY TYPE OPTIONS INCLUDE:',/,
     *        1X,'NONE........................(0)',/,
     *        1X,'"GLASKAPILLAREN"............(1)',/,
     *        1X,'"PYREX".....................(2)',/,
     *        1X,'"AR GLASS"..................(3)',/,
     *        1X,'FUSED SILICA................(4)',/,
     *        1X,'LINDEMANN GLASS.............(5)')
      WRITE(IO6,9002)
 9002 FORMAT(1X,'ENTER CAPPILARY TYPE: ',$)
      READ(IO5,9003) NKAP
 9003 FORMAT(I2)
      IF (NKAP.EQ.0) RETURN
      WRITE(IO6,9004)
 9004 FORMAT(1X,'ALL DIMENSIONS AND DISPLACEMENTS IN MM',/,
     *       1X,'ENTER THE DIAMETER OF CAPILLARY: ',$)
      READ(IO5,9005) DKAP
 9005 FORMAT(F5.2)
      WRITE(IO6,9006)
 9006 FORMAT(1X,'ENTER THE CAPPILARY THICKNESS:   ',$)
      READ(IO5,9007) TKAP
 9007 FORMAT(F7.4)
      WRITE(IO6,9008)
 9008 FORMAT(1X,'ENTER THE TUBE DISPLACEMENT ALONG X[PHI]: ',$)
      READ(IO5,9007) XDSP
      WRITE(IO6,9009)
 9009 FORMAT(1X,'ENTER THE TUBE DISPLACEMENT ALONG Y[PHI]: ',$)
      READ(IO5,9007) YDSP
      WRITE(IO6,9010)
 9010 FORMAT(1X,'ENTER TOP ARC ANGLE (SGM): ',$)
      READ(IO5,9007) SGM
      WRITE(IO6,9011)
 9011 FORMAT(1X,'ENTER LOW ARC ANGLE (TAU): ',$)
      READ(IO5,9007) TAU
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE MATX1
C.....
C  .. GET THE DIFFRACTOMETER ORIENTATION MATRIX
C.....
      SUBROUTINE MATX1
      DIMENSION ORMAT(3,3)
      COMMON /INPT09/ ORMAT, NDIFF
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      DO I=1,3
      DO J=1,3
        ORMAT(I,J)=0.0
      END DO
      END DO
      WRITE(IO6,9001)
 9001 FORMAT(/1X,'ENTER THE DIFFRACTOMETER ORIENTATION MATRIX')
 2000 DO 2500 I=1,3
      DO 2200 J=1,3
      WRITE(IO6,9015) I,J
 9015 FORMAT(1X,'ENTER ELEMENT(',I1,',',I1,'): ',$)
      READ(IO5,9016) ELEM
 9016 FORMAT(F10.7)
      ORMAT(I,J) = ELEM
 2200 CONTINUE
 2500 CONTINUE
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE MMSX0
C.....
C  .. CALL THE PROPER MMS-X GENERATING ROUTINES
C.....
      SUBROUTINE MMSX0
      DIMENSION ORMAT(3,3)
      COMMON /INPT09/ ORMAT, NDIFF
      DIMENSION ROT(3,3), SI(3), SD(3)
      COMMON /MMSXS/ SCALE, ROT, SI, SD
      DATA SCALE/10./
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      DIMENSION ANGLS(4), RMTF(3,3)
      CHARACTER*1 STD, BRG, ANS
      DATA STD/'S'/, BRG/'B'/
      DATA CON/57.2957795/
C-
      CALL CELL2
      CALL FACE2
      CALL EMUX2
      CALL MATX2
      WRITE(IO6,9001)
 9001 FORMAT(/1X,'MMS-X OPTIONS: ',/,
     *        1X,'STANDARD VIEW...........(S)',/,
     *        1X,'BRAGG CONDITION VIEW....(B)'   )
      WRITE(IO6,9002)
 9002 FORMAT(1X,'ENTER VIEW (S) OR VIEW (B): ',$)
      READ(IO5,9003) ANS
 9003 FORMAT(A)
      IF (ANS .EQ. STD) GO TO 1000
      IF (ANS .EQ. BRG) GO TO 2000
      RETURN
C-    STANDARD VIEW OF CRYSTAL
 1000 CALL GRID2
      DO 1002 I=1,3
      DO 1001 J=1,3
      ROT(I,J) = 0.0
 1001 CONTINUE
 1002 CONTINUE
      ROT(1,1) = +1.
      ROT(2,2) = +1.
      ROT(3,3) = +1.
      CALL MMSX1
      CALL MMSX3
      RETURN
C-    BRAGG VIEW OF CRYSTAL
 2000 WRITE(IO6,101)
  101 FORMAT(1X,'ENTER H K L (3F): ',$)
      READ(IO5,102) FH, FK, FL
  102 FORMAT(3F7.0)
      WRITE(IO6,103)
  103 FORMAT(1X,'ENTER PSI (1F): ',$)
      READ(IO5,104) PSI
  104 FORMAT(1F10.0)
      IF (NDIFF .EQ. 1)  GO TO 3000
      IF (NDIFF .EQ. 2)  GO TO 4000
      IF (NDIFF .EQ. 3)  GO TO 5000
      RETURN
C-    NICOLET P3/F DIFFRACTOMETER
 3000 CALL P3GEN (FH,FK,FL,PSI, ANGLS)
      WRITE(IO6,105) FH,FK,FL,PSI,ANGLS
  105 FORMAT(1X,'H K L: ',3F8.3,/,
     *       1X,'PSI: ',F8.3,/,
     *       1X,'P3 ANGLES: ',4F8.3 )
      CALL P3RMT (ANGLS, RMTF,ST,CT)
      SI(1) = -CT
      SI(2) = +ST
      SI(3) = 0.0
      SD(1) = +CT
      SD(2) = +ST
      SD(3) = 0.0
      GO TO 6000
C-    NONIUS CAD-4 DIFFRACTOMETER
 4000 CALL C4GEN (FH,FK,FL,PSI, ANGLS)
      WRITE(IO6,106) FH,FK,FL,PSI,ANGLS
  106 FORMAT(1X,'H K L: ',3F8.3,/,
     *       1X,'PSI: ',F8.3,/,
     *       1X,'CAD-4 ANGLES: ',4F8.3 )
      CALL C4RMT (ANGLS, RMTF,ST,CT)
      SI(1) = +CT
      SI(2) = +ST
      SI(3) = 0.0
      SD(1) = -CT
      SD(2) = +ST
      SD(3) = 0.0
      GO TO 6000
C-    BUSING-LEVY DIFFRACTOMETER
 5000 CALL BLGEN (FH,FK,FL,PSI, ANGLS)
      WRITE(IO6,107) FH,FK,FL,PSI,ANGLS
  107 FORMAT(1X,'H K L: ',3F8.3,/,
     *       1X,'PSI: ',F8.3,/,
     *       1X,'BL ANGLES: ',4F8.3 )
      CALL BLRMT (ANGLS, RMTF,ST,CT)
      SI(1) = +ST
      SI(2) = -CT
      SI(3) = 0.0
      SD(1) = +ST
      SD(2) = +CT
      SD(3) = 0.0
      GO TO 6000
C-
 6000 DO 7001 I=1,3
      DO 7002 J=1,3
      ROT(I,J) = RMTF(I,J)
 7002 CONTINUE
 7001 CONTINUE
      CALL MMSX2
      CALL MMSX3
      RETURN
      END
C.....
C  ..           SUBROUTINE MMSX1
C.....
C  .. GENERATE THE VERTICES, CONNECTION TABLE, AND
C  .. THE GAUSSIAN GRID POINTS TO ALLOW THE CRYSTAL
C  .. TO BE SEEN ON THE MMSX DEVICE.
C  .. CRYSTAL AS SEEN IN [PHI] AXIAL SYSTEM
C.....
      SUBROUTINE MMSX1
      CHARACTER TITLE*80, ATIME*8, ADATE*9
      COMMON /INPT01/ TITLE, ATIME, ADATE
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
      DIMENSION UBF(3,3), UBI(3,3), UBGF(3,3), UBGI(3,3)
      COMMON /MATRIX/ UBF, UBI, UBGF, UBGI
      DIMENSION VR(150,3), MM(150,3)
      COMMON /VERTEX/ NV, VR, MM
      PARAMETER (NGMAX=32**3)
      DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3)
      COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC
      DIMENSION ROT(3,3), SI(3), SD(3)
      COMMON /MMSXS/ SCALE, ROT, SI, SD
C-
      CHARACTER VL1*1, VL2*2, BLANK*1
      CHARACTER ORG*1, XAX*1, YAX*1, ZAX*1
      DIMENSION VXYZ(150,3), NCON(150,8), NZRO(8)
      DIMENSION RMTX(3,3)
      DIMENSION NPR(150), VCRY(3), VNOR(3), NORG(8)
      DATA NZRO/8*0/
      DATA VL1/'V'/, VL2/'V0'/
      DATA BLANK/' '/, ORG/' '/, XAX/'X'/, YAX/'Y'/, ZAX/'Z'/
C-
      NTOT = NV + 4
      CALL MTM1(ROT,UBGF, RMTX)
C-
      DO 10 I=1,150
      NPR(I) = 0
      DO 5 J=1,8
      NORG(J) = 0
      NCON(I,J) = 0
    5 CONTINUE
   10 CONTINUE
C-
      DO 100 N=1,NV
      DO 50 J=1,3
      VCRY(J) = VR(N,J) - VCOM(J)
   50 CONTINUE
      CALL MTV1 (RMTX,VCRY, VNOR)
      DO 75 J=1,3
      VXYZ(N,J) = SCALE*VNOR(J)
   75 CONTINUE
  100 CONTINUE
C-
      DO 200 I=1,NV-1
      MI1 = MM(I,1)
      MI2 = MM(I,2)
      MI3 = MM(I,3)
      DO 150 J=I+1,NV
      MJ1 = MM(J,1)
      MJ2 = MM(J,2)
      MJ3 = MM(J,3)
      CALL TWOCO (MI1,MI2,MI3, MJ1,MJ2,MJ3, NN)
      IF (NN.LE.0) GO TO 150
      NPR(I) = NPR(I) + 1
      NPR(J) = NPR(J) + 1
      NI = NPR(I)
      NJ = NPR(J)
      IF( NI.GE.9 .OR. NJ.GE.9 )  GO TO 999
      NCON(I,NI) = J
      NCON(J,NJ) = I
  150 CONTINUE
  200 CONTINUE
C-
      MMSX = 9
      OPEN(UNIT=MMSX,FILE='cryst.mms',STATUS='NEW',FORM='FORMATTED')
      WRITE(MMSX,9001) TITLE
      WRITE(MMSX,9002) NTOT
      DO 500 N=1,NV
      IF (N.LE.9)
     *  WRITE(MMSX,9003) VL2,N,(VXYZ(N,J),J=1,3),(NCON(N,J),J=1,8)
      IF (N.GE.10)
     *  WRITE(MMSX,9004) VL1,N,(VXYZ(N,J),J=1,3),(NCON(N,J),J=1,8)
  500 CONTINUE
C-
C     FIND MAXIMA ALONG [PHI] AXIAL SYSTEM VECTORS
C-
      XMAX=0.
      YMAX=0.
      ZMAX=0.
      DO 525 N=1,NV
      X = ABS(VXYZ(N,1))
      Y = ABS(VXYZ(N,2))
      Z = ABS(VXYZ(N,3))
      XMAX = MAX(X,XMAX)
      YMAX = MAX(Y,YMAX)
      ZMAX = MAX(Z,ZMAX)
  525 CONTINUE
C-
C     ORIGIN OF [PHI] AXIAL SYSTEM
C-
      VNOR(1) = 0.0
      VNOR(2) = 0.0
      VNOR(3) = 0.0
      NORG(1) = NV+2
      NORG(2) = NV+3
      NORG(3) = NV+4
      WRITE(MMSX,9005) ORG,VNOR,NORG
C-
C     X-AXIS OF [PHI]
C-
      VNOR(1) = 1.5*XMAX
      VNOR(2) = 0.0
      VNOR(3) = 0.0            
      NORG(1) = NV+1
      NORG(2) = 0
      NORG(3) = 0
      WRITE(MMSX,9005) XAX,VNOR,NORG
C-
C     Y-AXIS OF [PHI]
C
      VNOR(1) = 0.0
      VNOR(2) = 1.5*YMAX
      VNOR(3) = 0.0
      WRITE(MMSX,9005) YAX,VNOR,NORG
C-
C     Z-AXIS OF [PHI]
C
      VNOR(1) = 0.0
      VNOR(2) = 0.0
      VNOR(3) = 1.5*ZMAX
      WRITE(MMSX,9005) ZAX,VNOR,NORG
      WRITE(MMSX,9006) BLANK
      WRITE(MMSX,9006) BLANK
      CLOSE(UNIT=MMSX,STATUS='KEEP')
 9001 FORMAT(A)
 9002 FORMAT(I5)
 9003 FORMAT(A,I1,5X,3F8.3,16X,8I4)
 9004 FORMAT(A,I2,5X,3F8.3,16X,8I4)
 9005 FORMAT(A,2X,5X,3F8.3,16X,8I4)
 9006 FORMAT(A)
C-
  999 RETURN
      END
C.....
C  ..           SUBROUTINE MMSX2
C.....
C  .. GENERATE THE MMS-X FILE NECESSARY TO SEE THE CRYSTAL
C  .. IN THE BRAGG CONDITION OF A PARTICULAR INPUT REFLECTION.
C  .. OUTPUT IS A BRAGG.MMS FILE
C  .. CRYSTAL AS SEEN IN THE [THETA] AXIAL SYSTEM
C.....
      SUBROUTINE MMSX2
      CHARACTER TITLE*80, ATIME*8, ADATE*9
      COMMON /INPT01/ TITLE, ATIME, ADATE
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
      DIMENSION UBF(3,3), UBI(3,3), UBGF(3,3), UBGI(3,3)
      COMMON /MATRIX/ UBF, UBI, UBGF, UBGI
      DIMENSION VR(150,3), MM(150,3)
      COMMON /VERTEX/ NV, VR, MM
      PARAMETER (NGMAX=32**3)
      DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3)
      COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC
      DIMENSION ROT(3,3), SI(3), SD(3)
      COMMON /MMSXS/ SCALE, ROT, SI, SD
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      CHARACTER VL1*1, VL2*2, BLANK*1, COM*1, PLS*1
      CHARACTER SOUR*1, DETC*1
      DIMENSION VXYZ(150,3), NCON(150,8), NZRO(8)
      DIMENSION NSOUR(8), NDETC(8), NORGN(8)
      DIMENSION RMTX(3,3)
      DIMENSION NPR(150), VCRY(3), VNOR(3)
      DIMENSION RMTF(3,3), ANGLS(4)
      DATA NZRO/8*0/
      DATA NORGN/2,3,0,0,0,0,0,0/
      DATA NSOUR/1,0,0,0,0,0,0,0/
      DATA NDETC/1,0,0,0,0,0,0,0/
      DATA VL1/'V'/, VL2/'V0'/, COM/'*'/
      DATA BLANK/' '/, PLS/'+'/
      DATA SOUR/'S'/, DETC/'D'/
      DATA CON/57.2957795/
C-
      NTOT = NV + 3
      CALL MTM1(ROT,UBGF, RMTX)
C-
      DO 500 I=1,150
      NPR(I) = 0
      DO 450 J=1,8
      NCON(I,J) = 0
  450 CONTINUE
  500 CONTINUE
C-
      DO 600 N=1,NV
      DO 550 J=1,3
      VCRY(J) = VR(N,J) - VCOM(J)
  550 CONTINUE
      CALL MTV1 (RMTX,VCRY, VNOR)
      DO 575 J=1,3
      VXYZ(N,J) = SCALE*VNOR(J)
  575 CONTINUE
  600 CONTINUE
C-
      DO 800 I=1,NV-1
      MI1 = MM(I,1)
      MI2 = MM(I,2)
      MI3 = MM(I,3)
      DO 750 J=I+1,NV
      MJ1 = MM(J,1)
      MJ2 = MM(J,2)
      MJ3 = MM(J,3)
      CALL TWOCO (MI1,MI2,MI3, MJ1,MJ2,MJ3, NN)
      IF (NN.LE.0) GO TO 750
      NPR(I) = NPR(I) + 1
      NPR(J) = NPR(J) + 1
      NI = NPR(I)
      NJ = NPR(J)
      IF( NI.GE.9 .OR. NJ.GE.9 )  GO TO 999
      NCON(I,NI) = J + 3
      NCON(J,NJ) = I + 3
  750 CONTINUE
  800 CONTINUE
C-
      MMSX = 9
      OPEN(UNIT=MMSX,FILE='bragg.mms',STATUS='NEW',FORM='FORMATTED')
      WRITE(MMSX,9001) TITLE
      WRITE(MMSX,9002) NTOT
C-
      VCRY(1) = 0.0
      VCRY(2) = 0.0
      VCRY(3) = 0.0
      WRITE(MMSX,9005) COM,VCRY,NORGN
C-
      DO 900 L =1,3
      SI(L) = 1.01*SCALE*SI(L)
      SD(L) = 1.01*SCALE*SD(L)
  900 CONTINUE
      WRITE(MMSX,9005) SOUR,SI,NSOUR
      WRITE(MMSX,9005) DETC,SD,NDETC
C-
      DO 1000 N=1,NV
      IF (N.LE.9)
     *  WRITE(MMSX,9003) VL2,N,(VXYZ(N,J),J=1,3),(NCON(N,J),J=1,8)
      IF (N.GE.10)
     *  WRITE(MMSX,9004) VL1,N,(VXYZ(N,J),J=1,3),(NCON(N,J),J=1,8)
 1000 CONTINUE
C-
      WRITE(MMSX,9006) BLANK
      WRITE(MMSX,9006) BLANK
      CLOSE(UNIT=MMSX,STATUS='KEEP')
 9001 FORMAT(A)
 9002 FORMAT(I5)
 9003 FORMAT(A,I1,5X,3F8.3,16X,8I4)
 9004 FORMAT(A,I2,5X,3F8.3,16X,8I4)
 9005 FORMAT(A,2X,5X,3F8.3,16X,8I4)
 9006 FORMAT(A)
C-
  999 RETURN
      END
C.....
C  ..           SUBROUTINE MMSX3
C.....
C  .. GENERATE LINES REPRESENTING THE LIMITS OF THE CAPILLARY
C  .. SO THAT THESE MAY ALSO BE DISPLAYED ON THE MMSX SCREEN.
C.....
      SUBROUTINE MMSX3
      CHARACTER TITLE*80, ATIME*8, ADATE*9
      COMMON /INPT01/ TITLE, ATIME, ADATE
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
      DIMENSION UBF(3,3), UBI(3,3), UBGF(3,3), UBGI(3,3)
      COMMON /MATRIX/ UBF, UBI, UBGF, UBGI
      DIMENSION VR(150,3), MM(150,3)
      COMMON /VERTEX/ NV, VR, MM
      PARAMETER (NGMAX=32**3)
      DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3)
      COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC
      COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP
      DIMENSION UPTSF(3,3), UPTSI(3,3), UBIG(3,3), DTUB(NGMAX,3)
      COMMON /KAPCOR/ EMUK, RAD2, RPT2, UPTSF, UPTSI, UBIG, DTUB
      DIMENSION ROT(3,3), SI(3), SD(3)
      COMMON /MMSXS/ SCALE, ROT, SI, SD
C-
      CHARACTER BLANK*1, PLUS*1, MINUS*1
      DIMENSION NZRO(7), TCOR(3), PCOR(3)
      DATA NZRO/7*0/, NLN/72/
      DATA BLANK/' '/, PLUS/'+'/, MINUS/'-'/
      DATA CON/57.29577951/
C-    IF THERE IS NO CAPILLARY RETURN
      IF (NKAP .EQ. 0)   RETURN
C-    GENERATE EVERYTHING YOU NEED FOR TUBE
      CALL KAPP2
      ZMIN = -2.0 * DKAP
      ZMAX = +2.5 * DKAP
C-    OPEN FILE AND START WRITING
      MMSX = 9
      OPEN(UNIT=MMSX,FILE='tube.mms',STATUS='NEW',FORM='FORMATTED')
      WRITE(MMSX,9001) TITLE
      NTOT = 2*NLN
      WRITE(MMSX,9002) NTOT
      RKAP = 0.5 * DKAP
      FLN = NLN
      FLN = 360. / FLN
      DO 500 J=1,NLN
      ANG = (J-1) * FLN
      TCOR(1) = RKAP * COS(ANG/CON)
      TCOR(2) = RKAP * SIN(ANG/CON)
      TCOR(3) = ZMIN
      CALL MTV1 (UPTSF,TCOR, PCOR)
      PCOR(1) = PCOR(1) + XDSP
      PCOR(2) = PCOR(2) + YDSP
      CALL MTV1 (ROT,PCOR, TCOR)
      DO 450 L=1,3
      TCOR(L) = SCALE * TCOR(L)
  450 CONTINUE
      NCON = (J-1)*2 + 2
      WRITE(MMSX,9003) MINUS, TCOR, NCON, NZRO
      TCOR(1) = RKAP * COS(ANG/CON)
      TCOR(2) = RKAP * SIN(ANG/CON)
      TCOR(3) = ZMAX
      CALL MTV1 (UPTSF,TCOR, PCOR)
      PCOR(1) = PCOR(1) + XDSP
      PCOR(2) = PCOR(2) + YDSP
      CALL MTV1 (ROT,PCOR, TCOR)
      DO 475 L=1,3
      TCOR(L) = SCALE * TCOR(L)
  475 CONTINUE
      NCON = (J-1)*2 + 1
      WRITE(MMSX,9003) PLUS, TCOR, NCON, NZRO
  500 CONTINUE
      WRITE(MMSX,9006) BLANK
      WRITE(MMSX,9006) BLANK
      CLOSE(UNIT=MMSX,STATUS='KEEP')
 9001 FORMAT(A)
 9002 FORMAT(I5)
 9003 FORMAT(A,7X,3F8.3,16X,8I4)
 9006 FORMAT(A)
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE P3GEN
C.....
C  .. GIVEN H,K,L PSI AND THE ORIENTATION MATRIX
C  .. CALCULATE 2THETA OMEGA PHI AND CHI FOR THE
C  .. NICOLET P3/F DIFFRACTOMETER
C  .. CODE RIPPED OUT OF P3/F SOURCE CODE; SEE THE
C  .. ROUTINE TWOTH.FOR FOR INSTANCE
C.....
      SUBROUTINE P3GEN (FH,FK,FL,FPSI, ANGLS)
      DIMENSION UBF(3,3), UBI(3,3), UBGF(3,3), UBGI(3,3)
      COMMON /MATRIX/ UBF, UBI, UBGF, UBGI
      COMMON /LAMBDA/ FLAM
      DIMENSION ANGLS(4)
      DATA CONS/57.2957795/
C-
      TT = 0.
      OM = 0.
      PH = 0.
      CH = 0.
C-    GET CARTESIAN COMPONENTS OF RECIP VECTOR
      XI = FH*UBF(1,1) + FK*UBF(1,2) + FL*UBF(1,3)
      YI = FH*UBF(2,1) + FK*UBF(2,2) + FL*UBF(2,3)
      ZI = FH*UBF(3,1) + FK*UBF(3,2) + FL*UBF(3,3)
C-    CALCULATE DSTAR
      Q = SQRT(XI*XI + YI*YI + ZI*ZI)
C-    TEST FOR THE 000 REFLECTION
      IF (Q .LE. 0.00001) GO TO 999
C-    GET SIN THETA
      SINT = 0.5 * FLAM * Q
      IF (SINT .GE. 1.0) GO TO 999
C-    GET TWO THETA FROM SIN AND COS THETA
      COST = SQRT(1.0 - SINT*SINT)
      TT = 2.0*CONS*ATAN2(SINT,COST)
C-    GET OTHER ANGLES
      IF(YI.EQ.0.0) PHI0 =0.0
      IF(YI.NE.0.0)
     *   PHI0 = ATAN2(-XI,YI)
      SINC0 = ZI/Q
      PS = FPSI/CONS
      PHI = PHI0 + ATAN2(-SIN(PS), (COS(PS)*ABS(SINC0)))
      CHI = ATAN2(ZI, (YI*COS(PHI)-XI*SIN(PHI)))
      OMT = ATAN2((XI*COS(PHI)+YI*SIN(PHI)) , (ZI/SIN(CHI)))
      OM = AMOD(OMT*CONS+0.5*TT+360., 360.)
      PH = AMOD(PHI*CONS+360., 360.)
      CH = AMOD(CHI*CONS+360., 360.)
  999 ANGLS(1) = TT
      ANGLS(2) = OM
      ANGLS(3) = PH
      ANGLS(4) = CH
      RETURN
      END
C.....
C  ..           SUBROUTINE RISK1
C.....
C  .. RISK IT AND ACTUALLY CALCULATE A TRANSMISSION
C  .. FACTOR AND ASSOCIATED TBAR
C.....
      SUBROUTINE RISK1
C-
      CHARACTER TITLE*80, ATIME*8, ADATE*9
      CHARACTER FIL01*80, FIL02*80
      DIMENSION CELL(6)
      CHARACTER*4 WAVE
      CHARACTER*2 ATOMS(10)
      DIMENSION STOIC(10)
      DIMENSION FACES(50,4)
      DIMENSION NGP(3), D1(3),D2(3),D3(3)
      DIMENSION ORMAT(3,3)
      DIMENSION NTRP(50)
      COMMON /INPT01/ TITLE, ATIME, ADATE
      COMMON /INPT02/ FIL01, FIL02, IPATH
      COMMON /INPT03/ CELL
      COMMON /INPT04/ WAVE
      COMMON /LAMBDA/ FLAM
      COMMON /INPT05/ ATOMS
      COMMON /INPT06/ ZCELL, STOIC
      COMMON /INPT07/ FACES
      COMMON /INPT08/ NGP, D1,D2,D3
      COMMON /INPT09/ ORMAT, NDIFF
      COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP
      COMMON /INPT11/ NBPR, A0, A1, A2, A3
      COMMON /INPT12/ NSOL, EMUS, NTRP
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      DIMENSION UBF(3,3), UBI(3,3), UBGF(3,3), UBGI(3,3)
      DIMENSION FS(50,5)
      DIMENSION VR(150,3), MM(150,3)
      PARAMETER (NGMAX=32**3)
      DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3)
      DIMENSION UPTSF(3,3), UPTSI(3,3), UBIG(3,3), DTUB(NGMAX,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
      COMMON /CRYSTL/ TMOLW, RHO, EMUX
      COMMON /MATRIX/ UBF, UBI, UBGF, UBGI
      COMMON /PLANES/ NF, FS
      COMMON /VERTEX/ NV, VR, MM
      COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC
      COMMON /KAPCOR/ EMUK, RAD2, RPT2, UPTSF, UPTSI, UBIG, DTUB
C-
      DIMENSION ANGLS(4), SIP(3), SDP(3)
      DIMENSION SIR(3), SDR(3), SIS(3), SDS(3)
      DIMENSION SIT(3), SDT(3)
      DIMENSION PMP(3), PMS(3)
C-
      DOUBLE PRECISION ABSCO1,ABSCO2,ABSCO3,ABSCO4,ABSCOT,TBAR
C-
      WRITE(IO6,9001)
 9001 FORMAT(/1X,'CALCULATE TRANSMISSION AND TBAR',/,
     *        1X,'ENTER MILLER INDICES: ',$)
      READ(IO5,9002) FH,FK,FL
 9002 FORMAT(3F5.0)
      WRITE(IO6,9003) FH,FK,FL
 9003 FORMAT(/1X,'MILLER INDICES: ',3F5.0)
      DPSI=0
      PSIMAX=0
      WRITE(IO6,'(1X,''ENTER DELTA(PSI) (DEFAULT:  10 DEG):  ''$)')
      READ(IO5,'(F10.0)') DPSI
      IF (DPSI.EQ.0) DPSI=10
      IF (DPSI.LT.1) DPSI=1
      WRITE(IO6,'(1X,''ENTER PSI(MAX) (DEFAULT:  360 DEG):  ''$)')
      READ(IO5,'(F10.0)') PSIMAX
      IF (PSIMAX.LE.0.OR.PSIMAX.GT.360) PSIMAX=360
   80 WRITE(IO6,'(1X,''LINE PRINTER COPY OF RESULTS?  Y OR N.  ''$)')
      READ(IO5,'(A)') ANS
      IF (ANS.NE.'Y'.AND.ANS.NE.'N') GO TO 80
      ILP=0
      IF (ANS.EQ.'Y') THEN
        ILP=1
        OPEN (UNIT=NLP,FILE='absorb.lp',STATUS='UNKNOWN')
        DO J=1,1E9
          READ (NLP,'(A1)',END=85) DUMMY
        END DO
   85   WRITE (NLP,'(/1X)')
      END IF
      WRITE(IO6,9004)
 9004 FORMAT(1X,3X,'PSI',18X,'ANGLES',12X,'A CRYSTAL',4X,'A PROFILE',
     *      4X,'A KAPPLRY',4X,'A SOLVENT',4X,'A TOTALS',5X,'TBAR(MM)')
      IF (ILP.EQ.1) WRITE (NLP,9003) FH,FK,FL
      IF (ILP.EQ.1) WRITE (NLP,9004)
      DO 10 J=0,360
      FPSI=DPSI*J
      IF (FPSI.GT.PSIMAX) GO TO 11
      IF(NDIFF.EQ.1)
     * CALL P3GEN(FH,FK,FL,FPSI, ANGLS)
      IF(NDIFF.EQ.2)
     * CALL C4GEN(FH,FK,FL,FPSI, ANGLS)
      IF(NDIFF.EQ.3)
     * CALL BLGEN(FH,FK,FL,FPSI, ANGLS)
      IF (NDIFF.EQ.1)  CALL P3ROT (ANGLS, SIP,SDP,PMP)
      IF (NDIFF.EQ.2)  CALL C4ROT (ANGLS, SIP,SDP,PMP)
      IF (NDIFF.EQ.3)  CALL BLROT (ANGLS, SIP,SDP,PMP)
C-    GET REAL SPACE COMPONENTS OF SI AND SD: SIR AND SDR
      CALL MTV1 (UBGI,SIP, SIR)
      CALL MTV1 (UBGI,SDP, SDR)
C-    GET RECIPROCAL SPACE COMPONENTS OF PM: PMS
      CALL MTV1 (UBI,PMP,  PMS)
      ABSCO1 = 0.
      ABSCO2 = 0.
      ABSCO3 = 0.
      ABSCO4 = 0.
      ABSCOT = 0.
      TBAR   = 0.
      DO 100 N=1,NPTS
      W = WABC(N)
C-    CRYSTAL CONTRIBUTION
      CALL PATH1 (SIR,N, PXI,NPI)
      CALL PATH1 (SDR,N, PXD,NPD)
      PXTL = PXI + PXD
      AXTL = EXP( -EMUX * PXTL )
      TXTL = PXTL * AXTL
C-    CAPILLARY CONTRIBUTION
      AKAP = 1.
      IF(NKAP.EQ.0)  GO TO 50
      CALL MTV1 (UPTSI,SIP, SIT)
      CALL MTV1 (UPTSI,SDP, SDT)
      CALL PATH2 (SIT,N, W1I,W2I,PKI)
      CALL PATH2 (SDT,N, W1D,W2D,PKD)
      PKAP = PKI + PKD
      AKAP = EXP( -EMUK * PKAP )
   50 CONTINUE
C-    BEAM NONHOMOGENIETY CONTRIBUTION
      APRO = 1.
      IF(NBPR.EQ.0) GO TO 60
      CALL BPRO3 (PMS,N, DEL,APRO)
   60 CONTINUE
C-    TRAPPED SOLVENT CONTRIBUTION                    
      ASOL = 1.                                       
      IF(NSOL.EQ.0) GO TO 70                         
      CALL PATH3 (W1I, PSI)
      CALL PATH3 (W1D, PSD)
      PSOL = PSI + PSD
      ASOL = EXP( -EMUS * PSOL )
   70 CONTINUE
C-    TOTAL UP THE VARIOUS CONTRIBUTIONS
      ABSCO1 = ABSCO1 + W * AXTL
      ABSCO2 = ABSCO2 + W * APRO
      ABSCO3 = ABSCO3 + W * AKAP
      ABSCO4 = ABSCO4 + W * ASOL
      ABSCOT = ABSCOT + W * AXTL * APRO * AKAP * ASOL
      TBAR   = TBAR   + W * TXTL
  100 CONTINUE
C-
      AFSQ = FSQ  / ABSCOT
      ASIG = SIG  / ABSCOT
      TBAR = TBAR / ABSCO1
      WRITE(IO6,9005) FPSI, ANGLS,
     *       ABSCO1, ABSCO2, ABSCO3, ABSCO4, ABSCOT, TBAR
 9005 FORMAT(1X,F6.0,3X,4F7.2,6(3X,F10.7) )
      IF (ILP.EQ.1)
     *WRITE(NLP,9005) FPSI, ANGLS,
     *       ABSCO1, ABSCO2, ABSCO3, ABSCO4, ABSCOT, TBAR
   10 CONTINUE
   11 CONTINUE
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE SETUP
C.....
C  .. THE INPUT COMMAND DECODER ROUTINE.  ALL
C  .. OTHER ROUTINES ARE ACCESSED THROUGH IT
C  .. AND ALL MAJOR ROUTINES RETURN IT IT.
C.....
      SUBROUTINE SETUP
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      CHARACTER*1 REA, WRI, EDI, OUT
      CHARACTER*1 GEO, MMS, XXX, RES, TES
      DATA REA/'R'/, WRI/'W'/, EDI/'E'/, OUT/'O'/
      DATA GEO/'G'/, MMS/'M'/, XXX/'X'/, TES/'T'/
C-
      WRITE(IO6,9000)
 9000 FORMAT(1H1,/1X,'PROGRAM ABSORB1...SET UP ABSORB.DAT FILE')
C-
 1000 WRITE(IO6,9001)
 9001 FORMAT(/1X,'SETUP OPTIONS INCLUDE:',/,
     *        1X,'READ EXISTING ABSORB.DAT FILE....(R)',/,
     *        1X,'WRITE WORKING ABSORB.DAT FILE....(W)',/,
     *        1X,'EDIT WORKING ABSORB.DAT FILE.....(E)',/,
     *        1X,'CRYSTAL GEOMETRY CALCULATIONS....(G)',/,
     *        1X,'WRITE AN MMS-X FILE FOR VIEWING..(M)',/,
     *        1X,'OUTPUT A GEOMETRY SUMMARY........(O)',/,
     *        1X,'TEST THE ABSORPTION CALC.........(T)',/,
     *        1X,'EXIT FROM THIS SESSION...........(X)'    )
      WRITE(IO6,9002)
 9002 FORMAT(1X,'ENTER OPTION: ',$)
      READ(IO5,9003) RES
 9003 FORMAT(A)
      IF( RES.EQ.REA ) GO TO 1001
      IF( RES.EQ.WRI ) GO TO 1002
      IF( RES.EQ.EDI ) GO TO 1003
      IF( RES.EQ.GEO ) GO TO 1004
      IF( RES.EQ.MMS ) GO TO 1005
      IF( RES.EQ.OUT ) GO TO 1006
      IF( RES.EQ.TES ) GO TO 1007
      IF( RES.EQ.XXX ) GO TO 1008
      GO TO 1000
 1001 CALL PAGE1
      CALL CELL2
      CALL FACE2
      CALL GRID2
      CALL EMUX2
      CALL MATX2
      CALL KAPP2
      GO TO 1000
 1002 CALL PAGE2
      GO TO 1000
 1003 CALL EDIT1
      GO TO 1000
 1004 CALL GEOM1
      GO TO 1000
 1005 CALL MMSX0
      GO TO 1000
 1006 CALL OUTP1
      GO TO 1000
 1007 CALL RISK1
      GO TO 1000
 1008 CALL PAGE2
      RETURN
      END
C.....
C  ..           SUBROUTINE SOLN0
C.....
C  .. SUB-EDITOR FOR SOLVENT DATA
C.....
      SUBROUTINE SOLN0
      DIMENSION NTRP(50)
      COMMON /INPT12/ NSOL, EMUS, NTRP
      COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      CHARACTER*1 ALL, FLG, MUS, SFS, LST
      CHARACTER*1 XXX, RSP
      DATA ALL/'A'/, FLG/'F'/, MUS/'M'/
      DATA SFS/'S'/, LST/'L'/, XXX/'X'/
C-
 1000 WRITE(IO6,9001)
 9001 FORMAT(/1X,'TRAPPED SOLVENT SUB-EDITOR OPTIONS:',/,
     *        1X,'ALL NEW INPUT.....................(A)',/,
     *        1X,'SET TRAPPED SOLVENT FLAG NSOL.....(F)',/,
     *        1X,'ENTER MU OF SOLVENT...............(M)',/,
     *        1X,'SPECIFY TRAPPING FACES............(S)',/,
     *        1X,'LIST SOLVENT DATA.................(L)',/,
     *        1X,'EXIT TRAPPED SOLVENT SUB-EDITOR...(X)'    )
      WRITE(IO6,9002)
 9002 FORMAT(1X,'ENTER OPTION: ',$)
      READ(IO5,9003) RSP
 9003 FORMAT(A)
      IF(RSP.EQ.ALL)  GO TO 1001
      IF(RSP.EQ.FLG)  GO TO 1002
      IF(RSP.EQ.MUS)  GO TO 1003
      IF(RSP.EQ.SFS)  GO TO 1004
      IF(RSP.EQ.LST)  GO TO 1005
      IF(RSP.EQ.XXX)  RETURN
      GO TO 1000
C-
 1001 CALL SOLN1
      CALL SOLN2
      CALL SOLN3
      CALL SOLN4
      GO TO 1000
C-
 1002 CALL SOLN1
      GO TO 1000
 1003 CALL SOLN2
      GO TO 1000
 1004 CALL SOLN3
      GO TO 1000
 1005 CALL SOLN4
      GO TO 1000
C-
      END
C.....
C  ..           SUBROUTINE SOLN1
C.....
C  .. SPECIFY SOLVENT TRAPPING FLAG NSOL
C  ..      NSOL = 0  NO SOLVENT TRAPPED
C  ..      NSOL = 1  ONE OR MORE FACES TRAP SOLVENT
C.....
      SUBROUTINE SOLN1
      DIMENSION NTRP(50)
      COMMON /INPT12/ NSOL, EMUS, NTRP
      COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      WRITE(IO6,9001)
 9001 FORMAT(/1X,'TRAPPED SOLVENT FLAG:',/,
     *        1X,'NSOL = 0......NO SOLVENT TRAPPED',/,
     *        1X,'NSOL = 1......SOLVENT TRAPPED'      )
      WRITE(IO6,9002)
 9002 FORMAT(1X,'ENTER NSOL: ',$)
      READ(IO5,9003) NSOL
 9003 FORMAT(I5)
      IF(NSOL.NE.0) NSOL = 1
      IF(NSOL.EQ.0) GO TO 1001
      IF(NSOL.EQ.1) GO TO 1002
 1001 WRITE(IO6,9004)
 9004 FORMAT(1X,'NO TRAPPING SOLVENT')
      RETURN
 1002 WRITE(IO6,9005)
 9005 FORMAT(1X,'ONE OR MORE FACES TRAP SOLVENT')
      IF(NKAP.NE.0) RETURN
      WRITE(IO6,9006)
 9006 FORMAT(1X,'!!!!!WARNING!!!!!',/,
     *       1X,'NO CAPILLARY MEANS NO SOLVENT CORRECTION ')
      NSOL = 0
      RETURN
      END
C.....
C  ..           SUBROUTINE SOLN2
C.....
C  .. SPECIFY LINEAR ABSORPTION COEFFICIENT FOR SOLVENT
C  ..      EMUS IN RECIPROCAL MILLIMETERS
C.....
      SUBROUTINE SOLN2
      DIMENSION NTRP(50)
      COMMON /INPT12/ NSOL, EMUS, NTRP
      COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      EMUS=0.0
      IF (NSOL.EQ.0) RETURN
      WRITE(IO6,9001)
 9001 FORMAT(1X,'ENTER EMUS - MU OF SOLVENT(1/MM): ',$)
      READ(IO5,9003) EMUS
 9003 FORMAT(F10.6)
      IF(EMUS.LE.0.0) EMUS = 0.0
      RETURN
      END
C.....
C  ..           SUBROUTINE SOLN3
C.....
C  .. SPECIFY SOLVENT TRAPPING FACES
C.....
      SUBROUTINE SOLN3
      DIMENSION NTRP(50)
      COMMON /INPT12/ NSOL, EMUS, NTRP
      COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP
      DIMENSION FS(50,5)
      COMMON /PLANES/ NF, FS
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      DO I=1,50
        NTRP(I)=0
      END DO
      IF (NSOL.EQ.0.OR.EMUS.EQ.0.0) RETURN
      WRITE(IO6,9001)
 9001 FORMAT(/1X,'SPECIFY FACES TRAPPING SOLVENT ')
 1000 WRITE(IO6,9002)
 9002 FORMAT(1X,'ENTER FACE NUMBER (0 TO EXIT): ',$)
      READ(IO5,9003) MFN
 9003 FORMAT(I5)
      IF(MFN.LE.0) RETURN
      IF(MFN.GT.NF) GO TO 1000
      WRITE(IO6,9004)
 9004 FORMAT(1X,'NOTRAP(0) / TRAP(1): ',$)
      READ(IO5,9003) ITRP
      IF(ITRP.NE.0) ITRP = 1
      NTRP(MFN) = ITRP
      GO TO 1000
      END
C.....
C  ..           SUBROUTINE TITL1
C.....
C  .. GET THE TITLE FOR THIS JOB
C.....
      SUBROUTINE TITL1
      CHARACTER TITLE*80, ATIME*9, ADATE*8
      COMMON /INPT01/ TITLE, ATIME, ADATE
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C.....
      WRITE(IO6,9001)
 9001 FORMAT(/1X,'ENTER TITLE: ',$)
      READ(IO5,9002) TITLE
 9002 FORMAT(A)
      RETURN
      END
C.....
C  ..           SUBROUTINE WAVE1
C.....
C  .. INPUT THE WAVELENGTH SYMBOL (CU,MO,AG) AND THE
C  .. SPECTRAL LINE (KA,KB).  IF THE SYMBOL IS NOT ONE 
C  .. OF THE KNOWN TYPE IT CAN BE PUT IN.
C.....
      SUBROUTINE WAVE1
      CHARACTER*4 WAVE
      COMMON /INPT04/ WAVE
      COMMON /LAMBDA/ FLAM
      COMMON /CRYSTL/ TMOLW, RHO, EMUX
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      CHARACTER*4 SYMTB(6)
      DIMENSION WVLN(6)
      DATA WVLN/1.5418,1.3922,0.7107,0.6323,0.5608,0.4970/
      DATA SYMTB/'CUKA','CUKB','MOKA','MOKB','AGKA','AGKB'/
C-
    1 WRITE(IO6,9001)
 9001 FORMAT(/1X,'RECOGNIZED WAVELENGTH SYMBOLS ARE',/,
     *        1X,'COPPER K-ALPHA........CUKA',/,
     *        1X,'COPPER K-BETA.........CUKB',/,
     *        1X,'MOLYBDENUM K-ALPHA....MOKA',/,
     *        1X,'MOLYBDENUM K-BETA.....MOKB',/,
     *        1X,'SILVER K-ALPHA........AGKA',/,
     *        1X,'SILVER K-BETA.........AGKB',/,
     *        1X,'OTHER X-RAY...........XRAY',/,
     *        1X,'NEUTRON...............NEUT')
      WRITE(IO6,9002)
 9002 FORMAT(1X,'ENTER WAVELENGTH SYMBOL: ',$)
      READ(IO5,9003) WAVE
 9003 FORMAT(A)
      DO 10 J=1,6
      IF( WAVE.NE.SYMTB(J) ) GO TO 10
      FLAM = WVLN(J)
      RETURN
   10 CONTINUE
      WRITE(IO6,9004)
 9004 FORMAT(' ENTER WAVELENGTH (A):  '$)
      READ(IO5,9005) FLAM
 9005 FORMAT(F10.0)
      WRITE(IO6,9006)
 9006 FORMAT(' ENTER LINEAR ABSORPTION COEFFICIENT (MM**-1):  '$)
      READ(IO5,9005) EMUX
      RETURN
      END
C.....
C  ..           BLOCKDATA ABSORB3
C.....
      BLOCKDATA ABSORB3
C-
      CHARACTER TITLE*80, ATIME*8, ADATE*9
      CHARACTER FIL01*80, FIL02*80
      DIMENSION CELL(6)
      CHARACTER*4 WAVE
      CHARACTER*2 ATOMS(10)
      DIMENSION STOIC(10)
      DIMENSION FACES(50,4)
      DIMENSION NGP(3), D1(3),D2(3),D3(3)
      DIMENSION ORMAT(3,3)
      DIMENSION NTRP(50)
      COMMON /INPT01/ TITLE, ATIME, ADATE
      COMMON /INPT02/ FIL01, FIL02, IPATH
      COMMON /INPT03/ CELL
      COMMON /INPT04/ WAVE
      COMMON /LAMBDA/ FLAM
      COMMON /INPT05/ ATOMS
      COMMON /INPT06/ ZCELL, STOIC
      COMMON /INPT07/ FACES
      COMMON /INPT08/ NGP, D1,D2,D3
      COMMON /INPT09/ ORMAT, NDIFF
      COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP
      COMMON /INPT11/ NBPR, A0, A1, A2, A3
      COMMON /INPT12/ NSOL, EMUS, NTRP
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      DATA TITLE/'ABSORB1...VERSION OF 10/II/82'/
      DATA ATIME/'hh:mm:ss'/, ADATE/'dd-mmm-yy'/
      DATA FIL01/'FIL01.DATA'/, FIL02/'FIL02.DATA'/, IPATH/0/
      DATA CELL/1.,1.,1.,90.,90.,90./
      DATA WAVE/'MOKA'/
      DATA FLAM/0.7107/
      DATA ATOMS/10*'  '/
      DATA ZCELL/1./, STOIC/10*0./
      DATA FACES/200*0./
      DATA NGP/3*2/
      DATA D1/1.,0.,0./, D2/0.,1.,0./, D3/0.,0.,1./
      DATA ORMAT/9*0./
      DATA NDIFF/1/
      DATA NKAP/0/, DKAP/1./, TKAP/0.01/, TAU/0./, SGM/0./
      DATA XDSP/0./, YDSP/0./
      DATA NBPR/0/, A0/1./, A1/0./, A2/0./, A3/0./
      DATA NSOL/0/, EMUS/0./, NTRP/50*0/
      DATA IO5/5/, IO6/6/, NLP/20/, NAB/30/, NIN/1/, NOU/1/
      END
C.....
C  ..           SUBROUTINE AMI3
C.....
C  .. INVERT AN ASYMMETRIC 3 X 3 MATRIX
C.....
      SUBROUTINE AMI3(A,B,DET)
      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 .EQ. 0.0) RETURN
      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 10 I=1,3
      DO 15 J=1,3
      B(I,J) = (-1.0)**(I+J)*B(I,J)/DET
  15  CONTINUE
  10  CONTINUE
      RETURN
      END
C.....
C  ..           SUBROUTINE BLRMT
C.....
C  .. GET THE BUSING-LEVY ROTATION MATRIX FOR THE
C  .. BUSING-LEVY DIFFRACTOMETER ACTA CRYST 22,457-64(1967)
C  .. ANGLES ARE: TWO THETA, OMEGA, PHI, CHI
C.....
      SUBROUTINE BLRMT(ANGLS, RMTF,ST,CT)
      DIMENSION ANGLS(4), RMTF(3,3)
      DATA CONS/57.29577951/
C-
      TTH = ANGLS(1)/CONS
      OME = ANGLS(2)/CONS
      PHI = ANGLS(3)/CONS
      CHI = ANGLS(4)/CONS
      THT = 0.5*TTH
      ST = SIN(THT)
      CT = COS(THT)
      SO = SIN(OME)
      CO = COS(OME)
      SP = SIN(PHI)
      CP = COS(PHI)
      SC = SIN(CHI)
      CC = COS(CHI)
C-    EXPLICIT GENERATION OF ROTATION MATRIX
      RMTF(1,1) = +CO*CC*CP-SO*SP
      RMTF(1,2) = +CO*CC*SP+SO*CP
      RMTF(1,3) = +CO*SC
      RMTF(2,1) = -SO*CC*CP-CO*SP
      RMTF(2,2) = -SO*CC*SP+CO*CP
      RMTF(2,3) = -SO*SC
      RMTF(3,1) = -SC*CP
      RMTF(3,2) = -SC*SP
      RMTF(3,3) = +CC
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE BLROT
C.....
C  .. GET THE [PHI] FRAME COMPONENTS OF THE UNIT
C  .. VECTORS FOR THE BUSING-LEVY DIFFRACTOMETER
C  .. -SI- REVERSE INCIDENT BEAM
C  .. -SD- DIFFRACTED BEAM
C  .. -PM- MONOCHROMATOR VERTICAL PROFILE
C  .. ANGLES ARE: TWO THETA, OMEGA, PHI, CHI
C.....
      SUBROUTINE BLROT(ANGLS, SIP,SDP,PMP)
      DIMENSION ANGLS(4), SIP(3), SDP(3), PMP(3)
      DIMENSION RMTF(3,3), RMTI(3,3), SDT(3), SIT(3), PMT(3)
C-
      CALL BLRMT (ANGLS, RMTF,ST,CT)
C-
C-    TRANSPOSE OF ROTATION MATRIX EQUALS ITS INVERSE
      DO 100 I=1,3
      DO 50  J=1,3
      RMTI(I,J) = RMTF(J,I)
   50 CONTINUE
  100 CONTINUE
C-
      SIT(1) = +ST
      SIT(2) = -CT
      SIT(3) = 0.0
      SDT(1) = +ST
      SDT(2) = +CT
      SDT(3) = 0.0
      PMT(1) = 0.0
      PMT(2) = 0.0
      PMT(3) = 1.0
C-
      CALL MTV1(RMTI,SIT, SIP)
      CALL MTV1(RMTI,SDT, SDP)
      CALL MTV1(RMTI,PMT, PMP)
      RETURN
      END
C.....
C  ..           SUBROUTINE BPRO2
C.....
C  .. WRITE THE ANALYTICAL FORM OF THE BEAM PROFILE
C  .. THE EQUATION IS OF THE FORM
C  ..  F = A0 + A1*DEL +A2*DEL*DEL +A3*DEL*DEL*DEL
C.....
      SUBROUTINE BPRO2
      COMMON /INPT11/ NBPR, A0, A1, A2, A3
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      IF(NBPR.EQ.0) GO TO 100
      IF(NBPR.EQ.1) GO TO 200
      IF(NBPR.EQ.2) GO TO 300
      IF(NBPR.EQ.3) GO TO 400
      IF(NBPR.EQ.4) GO TO 500
  100 WRITE(NOU,9001)
 9001 FORMAT(/1X,'NO MONOCROMATOR / NO BEAM PROFILE')
      RETURN
  200 WRITE(NOU,9002) A0
 9002 FORMAT(/1X,'BEAM: SCALE CORRECTION / A0: ',F8.4)
      RETURN
  300 WRITE(NOU,9003) A0, A1
 9003 FORMAT(/1X,'BEAM: LINEAR CORRECTION / A0,A1: ',2F8.4)
      RETURN
  400 WRITE(NOU,9004) A0, A1, A2
 9004 FORMAT(/1X,'BEAM: QUADRATIC CORRECTION / A0,A1,A2: ',3F8.4)
      RETURN
  500 WRITE(NOU,9005) A0, A1, A2, A3
 9005 FORMAT(/1X,'BEAM: CUBIC CORRECTION /A0,A1,A2,A3: ',4F8.4)
      RETURN
      END
C.....
C  ..          SUBROUTINE BPRO3
C.....
C  .. CALCULATE THE BEAM PROFILE CONTRIBUTION ASSUMING A
C  .. PERFECTLY PARALLEL BEAM AND A NONHOMOGENIETY IN ONE
C  .. DIRECTION ONLY
C  ..        DIFFRACTOMETER         NONHOMOGENIETY DIRECTION
C  ..          NICOLET P3                 HORIZONTAL
C  ..          NONIUS CAD4                 VERTICAL 
C  ..          BUSING-LEVY                 VERTICAL
C.....
C  .. BEAM INTENSITY GIVEN IN ANALYTICAL FORM FOR A LEAST-SQUARES FIT
C  .. OF COUNT RATE VERSUS DISPLACEMENT FROM THE DIFFRACTOMETER CENTER
C  .. COUNT RATE IS IN FRACTIONAL FORM
C  .. DISPLACEMENT IS IN MILLIMETERS
C  .. FIT IS TO AN  EQUATION OF THE FORM
C  ..     BPF = A0 + A1*DEL + A2*DEL**2 + A3*DEL**3
C  ..     BPF IS THE FRACTIONAL COUNT RATE
C  ..     DEL IS THE PROJECTION OF POINT RABC ON NONHOMOGENIETY AXIS
C  ..     A0, A1, A2, A3 ARE FIT TO EXPERIMENTAL POINTS 
C.....
      SUBROUTINE BPRO3 (PM,N, DEL,BPF)
      COMMON /INPT11/ NBPR, A0, A1, A2, A3
      PARAMETER (NGMAX=32**3)
      DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3)
      COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC
      DIMENSION PM(3)
C-
      DEL = PM(1)*RABC(N,1) + PM(2)*RABC(N,2) + PM(3)*RABC(N,3)
      BPF = 1.
      IF(NBPR.EQ.0) RETURN
      BPF = A0
      IF(NBPR.EQ.1) RETURN
      BPF = BPF + A1*DEL
      IF(NBPR.EQ.2) RETURN
      BPF = BPF + A2*DEL*DEL
      IF(NBPR.EQ.3) RETURN
      BPF = BPF + A3*DEL*DEL*DEL
      RETURN
      END
C.....
C  ..           SUBROUTINE C4RMT
C.....
C  .. GET THE BUSING-LEVY ROTATION MATRIX FOR
C  .. THE ENRAF-NONIUS CAD-4 DIFFRACTOMETER
C  .. ANGLES ARE: THETA, PHIK, OMGK, KAPPA
C.....
      SUBROUTINE C4RMT(ANGLS, RMTF,ST,CT)
      DIMENSION ANGLS(4), RMTF(3,3)
      DIMENSION A(3,3), B(3,3), C(3,3), D(3,3), E(3,3)
      DIMENSION DE(3,3), CDE(3,3), BCDE(3,3), ABCDE(3,3)
      CON = 57.2957795
      ALPHA = 49.96681/CON
C-
      THETA = ANGLS(1)/CON
      PHIK  = ANGLS(2)/CON
      OMGK  = ANGLS(3)/CON
      FKAP  = ANGLS(4)/CON
C
C     CORRECTION OF CODE 10 MARCH 1987 (GDT)
C     OMEGA(KAPPA) IS REALLY OMEGA(BUSING-LEVY) PLUS THETA
C
      OMGK  = OMGK - THETA
C-
      SA = SIN(ALPHA)
      CA = COS(ALPHA)
      ST = SIN(THETA)
      CT = COS(THETA)
      SP = SIN(PHIK)
      CP = COS(PHIK)
      SO = SIN(OMGK)
      CO = COS(OMGK)
      SK = SIN(FKAP)
      CK = COS(FKAP)
C-
      A(1,1) = +CO
      A(1,2) = +SO
      A(1,3) = +0.0
      A(2,1) = -SO
      A(2,2) = +CO
      A(2,3) = +0.0
      A(3,1) = +0.0
      A(3,2) = +0.0
      A(3,3) = +1.0
C-
      B(1,1) = +CA
      B(1,2) = +0.0
      B(1,3) = +SA
      B(2,1) = +0.0
      B(2,2) = +1.0
      B(2,3) = +0.0
      B(3,1) = -SA
      B(3,2) = +0.0
      B(3,3) = +CA
C-
      C(1,1) = +CK
      C(1,2) = +SK
      C(1,3) = +0.0
      C(2,1) = -SK
      C(2,2) = +CK
      C(2,3) = +0.0
      C(3,1) = +0.0
      C(3,2) = +0.0
      C(3,3) = +1.0
C-
      D(1,1) = +CA
      D(1,2) = +0.0
      D(1,3) = -SA
      D(2,1) = +0.0
      D(2,2) = +1.0
      D(2,3) = +0.0
      D(3,1) = +SA
      D(3,2) = +0.0
      D(3,3) = +CA
C-
      E(1,1) = +CP
      E(1,2) = +SP
      E(1,3) = +0.0
      E(2,1) = -SP
      E(2,2) = +CP
      E(2,3) = +0.0
      E(3,1) = +0.0
      E(3,2) = +0.0
      E(3,3) = +1.0
C-
      CALL MTM1 (D,E, DE)
      CALL MTM1 (C,DE, CDE)
      CALL MTM1 (B,CDE, BCDE)
      CALL MTM1 (A,BCDE, ABCDE)
      DO 100 I=1,3
      DO  50 J=1,3
      RMTF(I,J) = ABCDE(I,J)
   50 CONTINUE
  100 CONTINUE
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE C4ROT
C.....
C  .. GET THE [PHI] FRAME COMPONENTS OF THE UNIT
C  .. VECTORS FOR THE ENRAF-NONIUS CAD-4 DIFFRACTOMETER
C  .. -SI- REVERSE INCIDENT BEAM
C  .. -SD- DIFFRACTED BEAM
C  .. -PM- MONOCHROMATOR VERTICAL PROFILE
C  .. ANGLES ARE: THETA, PHIK, OMGK, KAPPA
C.....
      SUBROUTINE C4ROT(ANGLS, SIP,SDP,PMP)
      DIMENSION ANGLS(4), SIP(3), SDP(3), PMP(3)
      DIMENSION RMTF(3,3), RMTI(3,3), SIT(3), SDT(3), PMT(3)
C-
      CALL C4RMT (ANGLS, RMTF,ST,CT)
C-    INVERSE OF ROTATION MATRIX EQUALS ITS TRANSPOSE
      DO 100 I=1,3
      DO  50 J=1,3
      RMTI(I,J) = RMTF(J,I)
   50 CONTINUE
  100 CONTINUE
C
      SIT(1) = +CT
      SIT(2) = +ST
      SIT(3) = +0.0
      SDT(1) = -CT 
      SDT(2) = +ST
      SDT(3) = +0.0
C
      PMT(1) = +0.0
      PMT(2) = +0.0
      PMT(3) = +1.0
C-
      CALL MTV1 (RMTI,SIT, SIP)
      CALL MTV1 (RMTI,SDT, SDP)
      CALL MTV1 (RMTI,PMT, PMP)
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE CELL2
C.....
C  .. GIVEN THE UNIT CELL CONSTANTS
C  .. GENERATE THE REAL AND RECIPROCAL CELL CONSTANTS
C  .. THE REAL AND RECIPROCAL METRIC TENSORS
C  .. AND THE VOLUMES.
C.....
      SUBROUTINE CELL2
      DIMENSION CELL(6)
      COMMON /INPT03/ CELL
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
      DATA CON/57.29577951/
C-    WRITE OUT TERMS EXPLICITLY
      A = CELL(1)
      B = CELL(2)
      C = CELL(3)
      ALP = CELL(4)
      BET = CELL(5)
      GAM = CELL(6)
C-    CONVERT ANGLES TO RADIANS
      CALP = COS (ALP/CON)
      SALP = SIN (ALP/CON)
      CBET = COS (BET/CON)
      SBET = SIN (BET/CON)
      CGAM = COS (GAM/CON)
      SGAM = SIN (GAM/CON)
      IF( ALP.EQ.90.) CALP = 0.
      IF( ALP.EQ.90.) SALP = 1.
      IF( BET.EQ.90.) CBET = 0.
      IF( BET.EQ.90.) SBET = 1.
      IF( GAM.EQ.90.) CGAM = 0.
      IF( GAM.EQ.90.) SGAM = 1.
C-    BUILD THE GR MATRIX EXPLICITLY
      GR(1,1) = A*A
      GR(1,2) = A*B*CGAM
      GR(1,3) = A*C*CBET
      GR(2,1) = GR(1,2)
      GR(2,2) = B*B
      GR(2,3) = B*C*CALP
      GR(3,1) = GR(1,3)
      GR(3,2) = GR(2,3)
      GR(3,3) = C*C
C-    INVERT GR TO GET GS AND VOL**2
      CALL AMI3 (GR, GS, V2)
      CVR = SQRT(V2)
      CVS = 1.0/CVR
C-    GET RECIPROCAL CELL EXPLICITLY
      AST = B*C*SALP*CVS
      BST = C*A*SBET*CVS
      CST = A*B*SGAM*CVS
      CALS = (CBET*CGAM - CALP) / (SBET*SGAM)
      SALS = SQRT(1. - CALS*CALS)
      CBES = (CALP*CGAM - CBET) / (SALP*SGAM)
      SBES = SQRT(1. - CBES*CBES)
      CGAS = (CALP*CBET - CGAM) / (SALP*SBET)
      SGAS = SQRT(1. - CGAS*CGAS)
      IF( ABS(CALS).LE.1E-8 ) CALS = 0.
      IF( ABS(CBES).LE.1E-8 ) CBES = 0.
      IF( ABS(CGAS).LE.1E-8 ) CGAS = 0.
      IF( CALS.EQ.0. ) SALS = 1.
      IF( CBES.EQ.0. ) SBES = 1.
      IF( CGAS.EQ.0. ) SGAS = 1.
C-    STUFF RR AND RS
      RR(1) = A
      RR(2) = B
      RR(3) = C
      RR(4) = CALP
      RR(5) = CBET
      RR(6) = CGAM
      RS(1) = AST
      RS(2) = BST
      RS(3) = CST
      RS(4) = CALS
      RS(5) = CBES
      RS(6) = CGAS
C-    BUILD THE B MATRIX EXPLICITLY
      BM(1,1) = AST
      BM(1,2) = BST*CGAS
      BM(1,3) = CST*CBES
      BM(2,1) = 0.0
      BM(2,2) = BST*SGAS
      BM(2,3) = -CST*SBES*CALP
      BM(3,1) = 0.0
      BM(3,2) = 0.0
      BM(3,3) = CST*SBES*SALP
      RETURN
      END
C.....
C  ..           SUBROUTINE CELL3
C.....
C  .. PRINT OR TYPE OUT THE CELL INFORMATION.
C.....
      SUBROUTINE CELL3
      DIMENSION CELL(6)
      COMMON /INPT03/ CELL
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      WRITE(NOU,102) CELL
  102 FORMAT(/1X,'A= ',3X,F7.3,4X,'B=',3X,F7.3,3X,'C=',4X,F7.3,/,
     *   1X,'ALPHA=',F7.3,3X,' BETA=',F7.3,3X,'GAMMA=',F7.3)
      WRITE(NOU,103)
  103 FORMAT(/1X,'REAL METRIC TENSOR')
      WRITE(NOU,104) ( (GR(I,J),J=1,3), I=1,3 )
  104 FORMAT(1X,3E15.7)
      WRITE(NOU,105)
  105 FORMAT(1X,'RECIPROCAL METRIC TENSOR')
      WRITE(NOU,104) ( (GS(I,J),J=1,3), I=1,3 )
      WRITE(NOU,106)
  106 FORMAT(1X,'BUSING-LEVY B MATRIX')
      WRITE(NOU,104) ( (BM(I,J),J=1,3), I=1,3)
      WRITE(NOU,107) CVR
  107 FORMAT(/1X,'REAL CELL VOL(A**3): ',E15.6)
      WRITE(NOU,108) CVS
  108 FORMAT(1X,'RECL. CELL VOL(A**-3):  ',E12.6)
      RETURN
      END
C.....
C  ..           SUBROUTINE LREAL
C.....
C  .. GET THE LENGTH OF A VECTOR GIVEN
C  .. ITS REAL SPACE COMPONENTS
C.....
      SUBROUTINE LREAL(FX,FY,FZ, DREL)
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
C-
      S1 = GR(1,1)*FX + GR(1,2)*FY + GR(1,3)*FZ
      S2 = GR(2,1)*FX + GR(2,2)*FY + GR(2,3)*FZ
      S3 = GR(3,1)*FX + GR(3,2)*FY + GR(3,3)*FZ
      D2 = FX*S1 + FY*S2 + FZ*S3
      IF (D2.LE.0.) D2 = 0.
      DREL = SQRT(D2)
      RETURN
      END
C.....
C  ..           SUBROUTINE LSTAR
C.....
C  .. CALCULATE THE LENGTH OF A RECIPROCAL VECTOR
C.....
      SUBROUTINE LSTAR(FH,FK,FL, DSTR)
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
      R1 = GS(1,1)*FH + GS(1,2)*FK + GS(1,3)*FL
      R2 = GS(2,1)*FH + GS(2,2)*FK + GS(2,3)*FL
      R3 = GS(3,1)*FH + GS(3,2)*FK + GS(3,3)*FL
      D2 = FH*R1 + FK*R2 + FL*R3
      IF (D2.LE.0.)  D2 = 0.
      DSTR = SQRT(D2)
      RETURN
      END
C.....
C  ..           SUBROUTINE EDGE1
C.....
C  .. GENERATE THE EDGES OF A CRYSTAL
C.....
      SUBROUTINE EDGE1
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
      DIMENSION VR(150,3), MM(150,3)
      COMMON /VERTEX/ NV, VR, MM
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      CHARACTER*1 STAR, BLNK, LABL, LABLS(250)
      DIMENSION JS(250), ELENS(250)
      DATA STAR/'*'/, BLNK/' '/
C-
      WRITE(NOU,9001)
 9001 FORMAT(/1X,'DIAGONALS AND EDGES OF CRYSTAL IN MILLIMETERS')
      DO 200 I=1,NV-1
      MI1 = MM(I,1)
      MI2 = MM(I,2)
      MI3 = MM(I,3)
      NS = 0
      DO 150 J=I+1,NV
      NS = NS+1
      MJ1 = MM(J,1)
      MJ2 = MM(J,2)
      MJ3 = MM(J,3)
      CALL TWOCO (MI1,MI2,MI3, MJ1,MJ2,MJ3, NN)
      IF (NN.EQ.1) LABL = STAR
      IF (NN.EQ.0) LABL = BLNK
      DELA = VR(I,1) - VR(J,1)
      DELB = VR(I,2) - VR(J,2)
      DELC = VR(I,3) - VR(J,3)
      CALL LREAL (DELA,DELB,DELC, ELEN)
      ELENS(NS) = ELEN
      LABLS(NS) = LABL
      JS(NS) = J
  150 CONTINUE
      WRITE(NOU,9003) I,(JS(N),ELENS(N),LABLS(N),N=1,NS)
 9003 FORMAT(1X,'VERTEX # ',I2,' TO',/,
     *       5(' #',I2,F7.3,A))
  200 CONTINUE
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE EMUX2
C.....
C  .. GIVEN THE CHEMICAL COMPOSITION OF THE CRYSTAL CHEMICAL UNIT
C  .. AND THE UNIT CELL VOLUME, CALCULATE THE MOLECULAR WEIGHT OF
C  .. THE CRYSTAL CHEMICAL UNIT, THE DENSITY OF THE CRYSTAL, AND
C  .. THE LINEAR ABSORPTION COEFFICIENT OF THE CRYSTAL.  THE VALUES
C  .. OF THE MASS ABSORPTION COEFFICIENTS FOR ELEMENTS 1 TO 96 ARE
C  .. TAKEN FROM THE INTERNATIONAL TABLES VOLUME IV AND ARE TABULATED
C  .. FOR COPPER, MOLYBDENUM AND SILVER RADIATION (BOTH K-ALPHA AND
C  .. K-BETA).  THE VALUES FOR THE ATOMIC WEIGHTS ARE FROM THE SARGENT
C  .. AND WELCH SCIENTIFIC COMPANY "PERIODIC TABLE OF THE ELEMENTS"
C  .. COPYRIGHT 1968 (WITH THE ATOMIC WEIGHT OF LITHIUM CORRECTED).
C  .. VALUES OF THE UNIT CELL VOLUME COME FROM SUBROUTINE CELL2.
C.....
      SUBROUTINE EMUX2
      CHARACTER*4 WAVE
      COMMON /INPT04/ WAVE
      COMMON /LAMBDA/ FLAM
      CHARACTER*2 ATOMS(10)
      COMMON /INPT05/ ATOMS
      DIMENSION   STOIC(10)
      COMMON /INPT06/ ZCELL, STOIC
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
C-
      COMMON /CRYSTL/ TMOLW, RHO, EMUX
      CHARACTER*2 SYMBL,BLANK,ATLAB(96)
      CHARACTER*4 CUKA,CUKB, MOKA,MOKB, AGKA,AGKB
      DIMENSION ECUKA(96), ECUKB(96)
      DIMENSION EMOKA(96), EMOKB(96)
      DIMENSION EAGKA(96), EAGKB(96)
      DIMENSION ATWHT(96), EWORK(96)
C-
      DATA CUKA/'CUKA'/, CUKB/'CUKB'/
      DATA MOKA/'MOKA'/, MOKB/'MOKB'/
      DATA AGKA/'AGKA'/, AGKB/'AGKB'/
      DATA ATLAB/'H ' ,'HE' ,'LI' ,'BE' ,'B ' ,'C ' ,
     *           'N ' ,'O ' ,'F ' ,'NE' ,'NA' ,'MG' ,
     *           'AL' ,'SI' ,'P ' ,'S ' ,'CL' ,'AR' ,
     *           'K ' ,'CA' ,'SC' ,'TI' ,'V ' ,'CR' ,
     *           'MN' ,'FE' ,'CO' ,'NI' ,'CU' ,'ZN' ,
     *           'GA' ,'GE' ,'AS' ,'SE' ,'BR' ,'KR' ,
     *           'RB' ,'SR' ,'Y ' ,'ZR' ,'NB' ,'MO' ,
     *           'TC' ,'RU' ,'RH' ,'PD' ,'AG' ,'CD' ,
     *           'IN' ,'SN' ,'SB' ,'TE' ,'I ' ,'XE' ,
     *           'CS' ,'BA' ,'LA' ,'CE' ,'PR' ,'ND' ,
     *           'PM' ,'SM' ,'EU' ,'GD' ,'TB' ,'DY' ,
     *           'HO' ,'ER' ,'TM' ,'YB' ,'LU' ,'HF' ,
     *           'TA' ,'W ' ,'RE' ,'OS' ,'IR' ,'PT' ,
     *           'AU' ,'HG' ,'TL' ,'PB' ,'BI' ,'  ' ,
     *           '  ' ,'RA' ,'  ' ,'  ' ,'  ' ,'TH' ,
     *           '  ' ,'U ' ,'  ' ,'PU' ,'  ' ,'  ' / 
      DATA AVOG /1.660565/
      DATA ECUKA/  0.3912,  0.2835,  0.4770,  1.0070,  2.1420,  4.2190,
     *             7.1420, 11.0300, 15.9500, 22.1300, 30.3000, 40.8800,
     *            50.2300, 65.3200, 77.2800, 92.5300,109.2000,119.5000,
     *           148.4000,171.4000,186.0000,202.4000,222.6000,252.3000,
     *           272.5000,304.4000,338.6000, 48.8300, 51.5400, 59.5100,
     *            62.1300, 67.9200, 75.6500, 82.8900, 90.2900, 97.0200,
     *           106.3000,115.3000,127.1000,136.8000,148.8000,158.3000,
     *           167.7000,180.8000,194.1000,205.0000,218.1000,229.3000,
     *           242.1000,253.3000,266.5000,273.4000,291.7000,309.8000,
     *           325.4000,336.1000,353.5000,378.8000,402.2000,417.9000,
     *           441.1000,453.5000,417.9000,426.7000,321.9000,336.6000,
     *           128.4000,134.3000,140.2000,144.7000,152.0000,157.7000,
     *           161.5000,170.5000,178.3000,183.8000,192.2000,198.2000,
     *           207.8000,216.2000,222.2000,232.1000,242.9000,  0.0000,
     *             0.0000,263.7000,  0.0000,  0.0000,  0.0000,306.8000,
     *             0.0000,305.7000,  0.0000,352.9000,  0.0000,  0.0000/
      DATA ECUKB/  0.3882,  0.2623,  0.3939,  0.7742,  1.5900,  3.0930,
     *             5.2150,  8.0620, 11.6600, 16.2400, 22.2300, 30.0800,
     *            37.1400, 48.3700, 57.4400, 68.9000, 81.7900, 89.3400,
     *           111.7000,129.0000,140.8000,153.2000,168.0000,191.1000,
     *           206.7000,233.6000,258.7000,282.8000, 38.7400, 45.3000,
     *            46.6500, 51.4400, 57.0100, 62.3200, 68.0700, 73.2200,
     *            80.1600, 86.7700, 96.1900,103.3000,112.3000,119.7000,
     *           126.9000,137.0000,147.1000,155.6000,165.8000,174.0000,
     *           183.3000,193.1000,203.6000,208.4000,221.9000,235.9000,
     *           247.5000,256.0000,270.8000,289.6000,306.8000,319.8000,
     *           336.4000,346.6000,369.0000,377.2000,399.9000,360.2000,
     *           272.4000,291.7000,288.5000,110.9000,116.5000,121.0000,
     *           123.9000,131.5000,137.3000,141.3000,147.6000,151.2000,
     *           160.6000,166.6000,171.1000,178.6000,187.5000,  0.0000,
     *             0.0000,203.0000,  0.0000,  0.0000,  0.0000,236.6000,
     *             0.0000,236.2000,  0.0000,271.2000,  0.0000,  0.0000/
      DATA EMOKA/  0.3727,  0.2019,  0.1968,  0.2451,  0.3451,  0.5348,
     *             0.7898,  1.1470,  1.5840,  2.2090,  2.9390,  3.9790,
     *             5.0430,  6.5330,  7.8700,  9.6250, 11.6400, 12.6200,
     *            16.2000, 19.0000, 21.0400, 23.2500, 25.2400, 29.2500,
     *            31.8600, 37.7400, 41.0200, 47.2400, 49.3400, 55.4600,
     *            56.9000, 60.4700, 65.9700, 68.8200, 74.6800, 79.1000,
     *            83.0000, 88.0400, 97.5600, 16.1000, 16.9600, 18.4400,
     *            19.7800, 21.3300, 23.0500, 24.4200, 26.3800, 27.7300,
     *            29.1300, 31.1800, 33.0100, 33.9200, 36.3300, 38.3100,
     *            40.4400, 42.3700, 45.3400, 48.5600, 50.7800, 53.2800,
     *            55.5200, 57.9600, 61.1800, 62.7900, 66.7700, 68.8900,
     *            72.1400, 75.6100, 78.9800, 80.2300, 84.1800, 86.3300,
     *            89.5100, 95.7600, 98.7400,100.2000,103.4000,108.6000,
     *           111.3000,114.7000,119.4000,122.8000,125.9000,  0.0000,
     *             0.0000,117.2000,  0.0000,  0.0000,  0.0000, 99.4600,
     *             0.0000, 96.6700,  0.0000, 48.8400,  0.0000,  0.0000/
      DATA EMOKB/  0.3699,  0.1972,  0.1866,  0.2216,  0.2928,  0.4285,
     *             0.6054,  0.8545,  1.1540,  1.5970,  2.0980,  2.8250,
     *             3.5850,  4.6240,  5.5690,  6.8350,  8.2610,  8.9490,
     *            11.5100, 13.5600, 15.0000, 16.6500, 18.0700, 20.9900,
     *            22.8900, 27.2100, 29.5100, 34.1800, 35.7700, 40.2600,
     *            41.6900, 44.2600, 48.5700, 51.2000, 55.5600, 58.6400,
     *            62.0700, 65.5900, 72.5700, 75.2000, 81.2200, 13.2900,
     *            14.3000, 15.4000, 16.6500, 17.6300, 19.1000, 20.1300,
     *            21.1800, 22.6200, 23.9100, 24.6700, 26.5300, 27.8600,
     *            29.5100, 31.0000, 33.1000, 35.5400, 37.0900, 38.8800,
     *            40.5200, 42.4000, 44.7400, 45.9500, 48.8800, 50.3800,
     *            52.7600, 55.0700, 57.9400, 59.2200, 62.0400, 64.1500,
     *            66.0700, 70.5700, 72.4700, 74.1300, 77.2000, 80.2300,
     *            82.3300, 85.3000, 88.2500, 90.5500, 93.5000,  0.0000,
     *             0.0000,100.7000,  0.0000,  0.0000,  0.0000, 73.3400,
     *             0.0000, 72.6300,  0.0000, 78.9900,  0.0000,  0.0000/
      DATA EAGKA/  0.3663,  0.1931,  0.1788,  0.2047,  0.2558,  0.3537,
     *             0.4756,  0.6479,  0.8507,  1.1630,  1.5010,  2.0040,
     *             2.5400,  3.2540,  3.9110,  4.8160,  5.8030,  6.2800,
     *             8.0800,  9.5730, 10.5400, 11.7600, 12.7800, 14.8800,
     *            16.2300, 19.3100, 20.9200, 24.3200, 25.5200, 28.7200,
     *            29.9300, 31.7800, 35.0000, 37.2300, 40.4100, 42.5300,
     *            45.3300, 47.8900, 52.9600, 55.2400, 59.6300, 62.5300,
     *            64.9900, 11.0000, 11.9100, 12.5800, 13.6700, 14.4500,
     *            15.2500, 16.1900, 17.0500, 17.7200, 19.1800, 20.0100,
     *            21.3100, 22.4500, 23.8500, 25.7200, 26.8000, 28.0200,
     *            29.2500, 30.6900, 32.4000, 33.3000, 35.4200, 36.5200,
     *            38.3100, 39.8300, 42.1200, 43.2100, 45.2700, 46.9300,
     *            48.2500, 51.4000, 52.7800, 54.2400, 56.7500, 58.5800,
     *            60.2100, 62.6300, 64.5600, 66.1400, 68.7300,  0.0000,
     *             0.0000, 74.1400,  0.0000,  0.0000,  0.0000, 86.2900,
     *             0.0000, 86.9200,  0.0000, 58.7600,  0.0000,  0.0000/
      DATA EAGKB/  0.3626,  0.1894,  0.1726,  0.1920,  0.2294,  0.3011,
     *             0.3851,  0.5041,  0.6401,  0.8605,  1.0880,  1.4360,
     *             1.8120,  2.3010,  2.7550,  3.4010,  4.0790,  4.4090,
     *             5.6670,  6.7520,  7.4020,  8.2980,  9.0240, 10.5200,
     *            11.4900, 13.6700, 14.7900, 17.2500, 18.1500, 20.4200,
     *            21.4100, 22.7300, 25.1100, 26.9400, 29.2400, 30.7100,
     *            32.9400, 34.8100, 38.4900, 40.4000, 43.6200, 45.7300,
     *            47.7300, 50.3400, 52.0300, 57.0300,  9.7800, 10.3700,
     *            10.9700, 11.5700, 12.1400, 12.7000, 13.8500, 14.3600,
     *            15.3700, 16.2500, 17.1500, 18.5900, 19.3400, 20.1600,
     *            21.0800, 22.1800, 23.4300, 24.1000, 25.6300, 26.4400,
     *            27.7900, 28.7700, 30.5800, 31.4600, 32.9800, 34.2400,
     *            35.1700, 37.3500, 38.3800, 39.6000, 74.5900, 42.6800,
     *            43.9300, 45.8700, 47.1500, 48.2200, 50.4100,  0.0000,
     *             0.0000, 54.4400,  0.0000,  0.0000,  0.0000, 62.8800,
     *             0.0000, 64.2400,  0.0000, 69.7500,  0.0000,  0.0000/
      DATA ATWHT/  1.0080,  4.0026,  6.9410,  9.0122, 10.8110, 12.0112,
     *            14.0067, 15.9994, 18.9984, 20.1830, 22.9898, 24.3120,
     *            26.9815, 28.0860, 30.9738, 32.0640, 35.4530, 39.9480,
     *            39.1020, 40.0800, 44.9560, 47.9000, 50.9420, 51.9960,
     *            54.9380, 55.8430, 58.9330, 58.7100, 63.5400, 65.3700,
     *            69.7200, 72.5900, 74.9220, 78.9600, 79.9090, 83.8000,
     *            85.4700, 87.6200, 88.9050, 91.2200, 92.9060, 95.9400,
     *            98.0000,101.0700,102.9050,106.4000,107.8700,112.4000,
     *           114.8200,118.6900,121.7500,127.6000,126.9040,131.3000,
     *           132.9050,137.3400,138.9100,140.1200,140.9070,144.2400,
     *           147.0000,150.3500,151.9600,157.2500,158.9240,162.5000,
     *           164.9300,167.2600,168.9340,173.0400,174.9700,178.4900,
     *           180.9480,183.8500,186.2000,190.2000,192.2000,195.0900,
     *           196.9670,200.5900,204.3700,207.1900,208.9800,000.0000,
     *           000.0000,222.0000,000.0000,000.0000,000.0000,232.0380,
     *           000.0000,238.0300,000.0000,242.0000,000.0000,000.0000/
      DATA WCUKA/1.5418/, WCUKB/1.3922/
      DATA WMOKA/0.7107/, WMOKB/0.6323/
      DATA WAGKA/0.5608/, WAGKB/0.4970/
      DATA EWORK/96*0.0/
C.....
C  .. NO CALCULATIONS IF NO ATOMS SPECIFIED
C.....
      IF (STOIC(1).LE.0.0)   RETURN
C.....
C-    GET THE MASS-ABSORPTION COEFFICIENTS AND WAVELENGTHS APPROPRIATE
C-    TO THE RADIATION USED
      IF (WAVE.EQ.CUKA) GO TO 100
      IF (WAVE.EQ.CUKB) GO TO 110
      IF (WAVE.EQ.MOKA) GO TO 120
      IF (WAVE.EQ.MOKB) GO TO 130
      IF (WAVE.EQ.AGKA) GO TO 140
      IF (WAVE.EQ.AGKB) GO TO 150
      GO TO 160
  100 FLAM = WCUKA
      DO 101 J=1,96
      EWORK(J) = ECUKA(J)
  101 CONTINUE
      GO TO 200
  110 FLAM = WCUKB
      DO 111 J=1,96
      EWORK(J) = ECUKB(J)
  111 CONTINUE
      GO TO 200
  120 FLAM = WMOKA
      DO 121 J=1,96
      EWORK(J) = EMOKA(J)
  121 CONTINUE
      GO TO 200
  130 FLAM = WMOKB
      DO 131 J=1,96
      EWORK(J) = EMOKB(J)
  131 CONTINUE
      GO TO 200
  140 FLAM = WAGKA
      DO 141 J=1,96
      EWORK(J) = EAGKA(J)
  141 CONTINUE
      GO TO 200
  150 FLAM = WAGKB
      DO 151 J=1,96
      EWORK(J) = EAGKB(J)
  151 CONTINUE
      GO TO 200
  160 DO 161 J=1,96
      EWORK(J) = 0.0
  161 CONTINUE
  200 CONTINUE
C-    ACCUMULATE THE SUMS NECESSARY
      SUM1 = 0.
      SUM2 = 0.
      DO 500 N=1,10
      SYMBL = ATOMS(N)
      FRACT = STOIC(N)
      DO 450 J=1,96
      IF( SYMBL.NE.ATLAB(J) )  GO TO 450
      ATW = ATWHT(J)
      FMU = EWORK(J)
  450 CONTINUE
      FAT = FRACT*ATW
      SUM1 = SUM1 + FAT
      SUM2 = SUM2 + FAT*FMU
  500 CONTINUE
C-    GET THE VARIOUS DERIVED DATA
      TMOLW = SUM1
      RHO = AVOG*ZCELL*SUM1/CVR
      IF( SUM2.EQ.0. ) RETURN
      EMUX = RHO*SUM2/SUM1
C-    CONVERT EMUX TO 1/MM
      EMUX = EMUX/10.
      RETURN
      END
C.....
C  ..           SUBROUTINE EMUX3
C.....
C  .. WRITE OUT THE LINEAR ABSORPTION COEFFICIENT AND
C  .. RELATED PERTINENT DATA
C.....
      SUBROUTINE EMUX3
      CHARACTER*4 WAVE
      COMMON /INPT04/ WAVE
      COMMON /LAMBDA/ FLAM
      CHARACTER*2 ATOMS(10)
      COMMON /INPT05/ ATOMS
      DIMENSION STOIC(10)
      COMMON /INPT06/ ZCELL,STOIC
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      COMMON /CRYSTL/ TMOLW, RHO, EMUX
C-
      NTYPE = 0
      DO 10 J=1,10
      IF(STOIC(J).NE.0.0) NTYPE=NTYPE+1
   10 CONTINUE
      WRITE(NOU,200) ZCELL
  200 FORMAT(/1X,'NUMBER OF C.C.U./UNIT CELL: ',F5.1)
      WRITE(NOU,201)
  201 FORMAT(1X,'CONTENTS OF C.C.U.')
      DO 700 J=1,NTYPE
      WRITE(NOU,202) ATOMS(J), STOIC(J)
  202 FORMAT(1X,A,3X,F5.1)
  700 CONTINUE
      WRITE(NOU,203) TMOLW
  203 FORMAT(1X,'MOLECULAR WEIGHT: ',F7.2)
      WRITE(NOU,204) RHO
  204 FORMAT(1X,'RHO (MG/MM**3): ',F7.4)
      WRITE(NOU,205) WAVE
  205 FORMAT(1X,'WAVE SYMBOL: ',2X,A)
      WRITE(NOU,206) FLAM
  206 FORMAT(1X,'LAMBDA (A): ',F6.4)
      WRITE(NOU,207) EMUX
  207 FORMAT(1X,'CRYSTAL MU (1/MM): ',F8.5)
      RETURN
      END
C.....
C  ..           SUBROUTINE FACE2
C.....
C  .. GIVEN THE MILLERIAN SYMBOLS OF THE CRYSTAL PLANES
C  .. GENERATE THE UNIT VECTORS IN THE [STAR] FRAME
C  .. WHICH ARE PERPENDICULAR TO THE PLANES.
C.....
C  ..       FS(N,1)...A* COMPONENT OF UNIT VECTOR
C  ..       FS(N,2)...B*    "      "   "     "
C  ..       FS(N,3)...C*    "      "   "     "
C  ..       FS(N,4)...PERPENDICULAR DISTANCE IN MM.
C  ..       FS(N,5)...DIST TO FACE FOR SOLVENT CALC
C.....
      SUBROUTINE FACE2
      DIMENSION FACES(50,4)
      COMMON /INPT07/ FACES
      DIMENSION FS(50,5)
      COMMON /PLANES/ NF, FS
C-
      NF = 0
      J = 1
   10 FH = FACES(J,1)
      FK = FACES(J,2)
      FL = FACES(J,3)
      PD = FACES(J,4)
      CALL LSTAR(FH,FK,FL, DSTR)
      IF( DSTR.LE.0. )  RETURN
C-
      FS(J,1) = FH/DSTR
      FS(J,2) = FK/DSTR
      FS(J,3) = FL/DSTR
      FS(J,4) = PD
      FS(J,5) = 0.0
      NF = NF + 1
      J  =  J + 1
      IF( J .LE. 50 ) GO TO 10
      RETURN
      END
C.....
C  ..           SUBROUTINE FACE3
C.....
C  .. PRINT OR TYPE OUT THE CRYSTAL FACE INFORMATION
C.....
      SUBROUTINE FACE3
      DIMENSION FACES(50,4)
      COMMON /INPT07/ FACES
      DIMENSION FS(50,5)
      COMMON /PLANES/ NF, FS
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      WRITE(NOU,9001)
 9001 FORMAT(/1X,'CRYSTAL FACE DATA')
      WRITE(NOU,9002)
 9002 FORMAT(1X,'  #',7X,'MILLER INDICES',5X,
     *'  LENGTH',5X,'[STAR] COMPONENTS OF UNIT VECTOR')
      DO 10 I=1,NF
      WRITE(NOU,9003) I,(FACES(I,J),J=1,4),(FS(I,J),J=1,3)
 9003 FORMAT(1X,I3,3F8.3,2X,F8.5,1X,3E13.4)
   10 CONTINUE
      RETURN
      END
C.....
C  ..           SUBROUTINE FILS2
C.....
C  .. OUTPUT THE NAMES OF THE INPUT AND OUTPUT FILES
C.....
      SUBROUTINE FILS2
      CHARACTER FIL01*80, FIL02*80
      COMMON /INPT02/ FIL01, FIL02, IPATH
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C....
      WRITE(NOU,9001) FIL01
      WRITE(NOU,9002) FIL02
      WRITE(NOU,9003) IPATH
 9001 FORMAT(/1X,'INPUT FILE NAME:   ',A)
 9002 FORMAT( 1X,'OUTPUT FILE NAME:  ',A)
 9003 FORMAT(/1X,
     &'IPATH = ',I1//1X,
     &'IPATH = 0  DO NOT'/1X,
     &'      = 1  DO'/1X,
     &'           WRITE TBAR AND S0 AND S1 COMPONENTS TO OUTPUT FILE')
      RETURN
      END
C.....
C  ..           SUBROUTINE FLIPS
C.....
C  .. GENERATE THE DUMMY FACE VECTOR TO BE
C  .. USED IN THE SECTIONING OF THE CRYSTAL.
C  .. IF THE DISTANCE IS NEGATIVE THEN FLIP
C  .. THE SIGNS OF THE INDICES OF THE PLANE.
C  .. THAT IS, AN HKL PLANE WITH A NEGATIVE
C  .. DISTANCE IS ACTUALLY AN -H-K-L PLANE
C  .. WITH A POSITIVE DISTANCE.
C.....
      SUBROUTINE FLIPS (DVEC,EPSL, FVEC)
      DIMENSION DVEC(3), FVEC(4)
C-
      FVEC(1) = DVEC(1)
      FVEC(2) = DVEC(2)
      FVEC(3) = DVEC(3)
      FVEC(4) = EPSL
      IF (EPSL.GE.0.0) RETURN
      DO 10 N=1,4
      FVEC(N) = -FVEC(N)
   10 CONTINUE
      RETURN
      END
C.....
C  ..          SUBROUTINE GAUTAB
C.....
C  .. GIVEN NGP, THE GAUSSIAN GRID NUMBER, TRANSFER THE
C  .. GAUSSIAN FRACTIONS AND GAUSSIAN WEIGHTS TO WORKING
C  .. ARRAYS PASSED BACK TO CALLING ROUTINE.  THE VALUES
C  .. IN THESE TABLES ARE CALCULATED FROM THE TABLES OF
C  .. ABRAMOWITZ AND STEGUN.  THE TABLES WERE PUNCHED BY
C  .. THAT STAR OF THE STAGE, OUR OWN CATHERINE DEVINE,
C  .. AND CONVERTED TO THE PRESENT FORM BY GDT.
C.....
      SUBROUTINE GAUTAB (NGP, FGAUSS,WGAUSS)
      DIMENSION FGAUSS(96), WGAUSS(96)
      DIMENSION NALLOW(19)
      DIMENSION F02( 2), W02( 2)
      DIMENSION F03( 3), W03( 3)
      DIMENSION F04( 4), W04( 4)
      DIMENSION F05( 5), W05( 5)
      DIMENSION F06( 6), W06( 6)
      DIMENSION F07( 7), W07( 7)
      DIMENSION F08( 8), W08( 8)
      DIMENSION F09( 9), W09( 9)
      DIMENSION F10(10), W10(10)
      DIMENSION F12(12), W12(12)
      DIMENSION F16(16), W16(16)
      DIMENSION F20(20), W20(20)
      DIMENSION F24(24), W24(24)
      DIMENSION F32(32), W32(32)
      DIMENSION F40(40), W40(40)
      DIMENSION F48(48), W48(48)
      DIMENSION F64(64), W64(64)
      DIMENSION F80(80), W80(80)
      DIMENSION F96(96), W96(96)
C     NGAUSS =  2
      DATA F02/.2113248,.7886752/
      DATA W02/.5000000,.5000000/
C     NGAUSS =  3
      DATA F03/.1127017,.5000000,.8872983/
      DATA W03/.2777778,.4444444,.2777778/
C     NGAUSS =  4
      DATA F04/.0694318,.3300095,.6699905,.9305682/
      DATA W04/.1739274,.3260726,.3260726,.1739274/
C     NGAUSS =  5
      DATA F05/.0469101,.2307653,.5000000,.7692347,.9530899/
      DATA W05/.1184634,.2393143,.2844445,.2393143,.1184634/
C     NGAUSS =  6
      DATA F06/.0337653,.1693953,.3806904,.6193096,.8306047,.9662347/
      DATA W06/.0856623,.1803808,.2339569,.2339569,.1803808,.0856623/
C     NGAUSS =  7
      DATA F07/.0254461,.1292344,.2970774,.5000000,.7029226,.8707656, 
     *         .9745539/
      DATA W07/.0647425,.1398527,.1909150,.2089796,.1909150,.1398527, 
     *         .0647425/
C     NGAUSS =  8
      DATA F08/.0198551,.1016667,.2372338,.4082827,.5917173,.7627662, 
     *         .8983333,.9801450/
      DATA W08/.0506142,.1111905,.1568533,.1813419,.1813419,.1568533, 
     *         .1111905,.0506142/
C     NGAUSS =  9
      DATA F09/.0159199,.0819845,.1933143,.3378733,.5000000,.6621267, 
     *         .8066857,.9180155,.9840801/
      DATA W09/.0406372,.0903241,.1303053,.1561735,.1651197,.1561735, 
     *         .1303053,.0903241,.0406372/
C     NGAUSS = 10
      DATA F10/.0130467,.0674683,.1602952,.2833023,.4255629,.5744371, 
     *         .7166977,.8397048,.9325317,.9869533/
      DATA W10/.0333356,.0747257,.1095432,.1346333,.1477621,.1477621, 
     *         .1346333,.1095432,.0747257,.0333356/
C     NGAUSS = 12
      DATA F12/.0092197,.0479413,.1150486,.2063410,.3160843,.4373833, 
     *         .5626167,.6839157,.7936590,.8849514,.9520587,.9907803/
      DATA W12/.0235876,.0534697,.0800392,.1015837,.1167462,.1245735, 
     *         .1245735,.1167462,.1015837,.0800392,.0534697,.0235876/
C     NGAUSS = 16
      DATA F16/.0052995,.0277125,.0671844,.1222978,.1910619,.2709916, 
     *         .3591982,.4524938,.5475063,.6408018,.7290084,.8089381, 
     *         .8777022,.9328156,.9722875,.9947005/
      DATA W16/.0135763,.0311268,.0475793,.0623145,.0747980,.0845783, 
     *         .0913017,.0947253,.0947253,.0913017,.0845783,.0747980, 
     *         .0623145,.0475793,.0311268,.0135763/
C     NGAUSS = 20
      DATA F20/.0034357,.0180140,.0438828,.0804415,.1268341,.1819731, 
     *         .2445665,.3131469,.3861071,.4617367,.5382633,.6138930, 
     *         .6868531,.7554335,.8180269,.8731660,.9195585,.9561172, 
     *         .9819860,.9965643/
      DATA W20/.0088070,.0203007,.0313360,.0416383,.0509650,.0590972, 
     *         .0658443,.0710481,.0745865,.0763767,.0763767,.0745865, 
     *         .0710481,.0658443,.0590972,.0509650,.0416383,.0313360, 
     *         .0203007,.0088070/
C     NGAUSS = 24
      DATA F24/.0024064,.0126357,.0308627,.0567923,.0899990,.1299379, 
     *         .1759531,.2272892,.2831033,.3424787,.4044406,.4679716, 
     *         .5320284,.5955595,.6575214,.7168968,.7727108,.8240469, 
     *         .8700621,.9100010,.9432077,.9691373,.9873644,.9975936/
      DATA W24/.0061706,.0142657,.0221387,.0296493,.0366732,.0430951, 
     *         .0488093,.0537221,.0577528,.0608352,.0629188,.0639691, 
     *         .0639691,.0629188,.0608352,.0577528,.0537221,.0488093, 
     *         .0430951,.0366732,.0296493,.0221387,.0142657,.0061706/
C     NGAUSS = 32
      DATA F32/.0013680,.0071943,.0176188,.0325469,.0518394,.0753162, 
     *         .1027581,.1339090,.1684778,.2061421,.2465501,.2893244, 
     *         .3340657,.3803563,.4277640,.4758461,.5241538,.5722360, 
     *         .6196437,.6659343,.7106757,.7534500,.7938579,.8315222, 
     *         .8660911,.8972420,.9246838,.9481606,.9674531,.9823812, 
     *         .9928058,.9986320/
      DATA W32/.0035093,.0081372,.0126961,.0171370,.0214180,.0254990, 
     *         .0293420,.0329111,.0361729,.0390969,.0416560,.0438260, 
     *         .0455869,.0469222,.0478193,.0482700,.0482700,.0478193, 
     *         .0469222,.0455869,.0438260,.0416560,.0390969,.0361729, 
     *         .0329111,.0293420,.0254990,.0214180,.0171370,.0126961, 
     *         .0081372,.0035093/
C     NGAUSS = 40
      DATA F40/.0008811,.0046369,.0113701,.0210416,.0335936,.0489506, 
     *         .0670202,.0876939,.1108471,.1363408,.1640216,.1937231, 
     *         .2252665,.2584621,.2931104,.3290030,.3659239,.4036512, 
     *         .4419580,.4806138,.5193862,.5580420,.5963488,.6340761, 
     *         .6709970,.7068896,.7415379,.7747335,.8062770,.8359784, 
     *         .8636592,.8891529,.9123061,.9329798,.9510494,.9664064, 
     *         .9789584,.9886299,.9953631,.9991189/
      DATA W40/.0022607,.0052491,.0082106,.0111229,.0139685,.0167301, 
     *         .0193911,.0219354,.0243479,.0266139,.0287199,.0306531, 
     *         .0324020,.0339560,.0353058,.0364433,.0373616,.0380552, 
     *         .0385199,.0387530,.0387530,.0385199,.0380552,.0373616, 
     *         .0364433,.0353058,.0339560,.0324020,.0306531,.0287199, 
     *         .0266139,.0243479,.0219354,.0193911,.0167301,.0139685, 
     *         .0111229,.0082106,.0052491,.0022607/
C     NGAUSS = 48
      DATA F48/.0006145,.0032349,.0079377,.0147042,.0235062,.0343066, 
     *         .0470605,.0617140,.0782059,.0964669,.1164205,.1379830, 
     *         .1610638,.1855663,.2113877,.2384195,.2665485,.2956568, 
     *         .3256221,.3563188,.3876181,.4193888,.4514976,.4838099, 
     *         .5161901,.5485023,.5806112,.6123819,.6436812,.6743780, 
     *         .7043433,.7334515,.7615805,.7886124,.8144338,.8389362, 
     *         .8620170,.8835795,.9035331,.9217942,.9382860,.9529396, 
     *         .9656934,.9764938,.9852958,.9920623,.9967651,.9993855/
      DATA W48/.0015766,.0036638,.0057386,.0077897,.0098081,.0117854, 
     *         .0137132,.0155836,.0173886,.0191207,.0207726,.0223373, 
     *         .0238084,.0251795,.0264451,.0275998,.0286386,.0295574, 
     *         .0303522,.0310197,.0315571,.0319621,.0322331,.0323688, 
     *         .0323688,.0322331,.0319621,.0315571,.0310197,.0303522, 
     *         .0295574,.0286386,.0275998,.0264451,.0251795,.0238084, 
     *         .0223373,.0207726,.0191207,.0173886,.0155836,.0137132, 
     *         .0117854,.0098081,.0077897,.0057386,.0036638,.0015766/
C     NGAUSS = 64
      DATA F64/.0003475,.0018300,.0044933,.0083318,.0133366,.0194956, 
     *         .0267943,.0352154,.0447389,.0553423,.0670003,.0796853, 
     *         .0933673,.1080138,.1235901,.1400591,.1573819,.1755172, 
     *         .1944223,.2140522,.2343602,.2552985,.2768170,.2988649, 
     *         .3213899,.3443386,.3676564,.3912882,.4151778,.4392686, 
     *         .4635035,.4878249,.5121751,.5364966,.5607314,.5848222, 
     *         .6087118,.6323436,.6556615,.6786101,.7011351,.7231830, 
     *         .7447016,.7656398,.7859478,.8055777,.8244828,.8426182, 
     *         .8599409,.8764099,.8919863,.9066327,.9203147,.9329997, 
     *         .9446577,.9552611,.9647846,.9732057,.9805045,.9866634, 
     *         .9916682,.9955067,.9981701,.9996525/
      DATA W64/.0008916,.0020735,.0032523,.0044234,.0055840,.0067315, 
     *         .0078630,.0089759,.0100687,.0111351,.0121763,.0131888, 
     *         .0141699,.0151174,.0160289,.0169026,.0177361,.0185276, 
     *         .0192751,.0199769,.0206313,.0212367,.0217918,.0222953, 
     *         .0227458,.0231424,.0234841,.0237701,.0239997,.0241724, 
     *         .0242878,.0243455,.0243455,.0242878,.0241724,.0239997, 
     *         .0237701,.0234841,.0231424,.0227458,.0222953,.0217918, 
     *         .0212367,.0206313,.0199769,.0192751,.0185276,.0177361, 
     *         .0169026,.0160289,.0151174,.0141699,.0131888,.0121763, 
     *         .0111351,.0100687,.0089759,.0078630,.0067315,.0055840, 
     *         .0044234,.0032523,.0020735,.0008916/
C     NGAUSS = 80
      DATA F80/.0002231,.0011750,.0028862,.0053543,.0085757,.0125454, 
     *         .0172575,.0227046,.0288786,.0357700,.0433685,.0516622, 
     *         .0606387,.0702843,.0805842,.0915230,.1030836,.1152488, 
     *         .1279998,.1413174,.1551812,.1695700,.1844621,.1998347, 
     *         .2156643,.2319270,.2485979,.2656517,.2830623,.3008033, 
     *         .3188476,.3371678,.3557360,.3745238,.3935027,.4126438, 
     *         .4319180,.4512958,.4707478,.4902443,.5097557,.5292522, 
     *         .5487042,.5680820,.5873562,.6064972,.6254762,.6442640, 
     *         .6628322,.6811524,.6991967,.7169377,.7343483,.7514021, 
     *         .7680730,.7843357,.8001653,.8155379,.8304300,.8448188, 
     *         .8586826,.8720002,.8847512,.8969164,.9084771,.9194158, 
     *         .9297157,.9393613,.9483379,.9566315,.9642300,.9711214, 
     *         .9772954,.9827425,.9874546,.9914243,.9946457,.9971138, 
     *         .9988250,.9997769/
      DATA W80/.0005725,.0013317,.0020902,.0028455,.0035964,.0043420, 
     *         .0050809,.0058120,.0065344,.0072467,.0079481,.0086374, 
     *         .0093134,.0099753,.0106220,.0112526,.0118660,.0124612, 
     *         .0130376,.0135791,.0141299,.0146442,.0151361,.0156051, 
     *         .0160503,.0164710,.0168666,.0172365,.0175803,.0178972, 
     *         .0181869,.0184489,.0186829,.0188882,.0190648,.0192125, 
     *         .0193309,.0194199,.0194792,.0195089,.0195089,.0194792, 
     *         .0194199,.0193309,.0192125,.0190648,.0188882,.0186829, 
     *         .0184489,.0181869,.0178972,.0175803,.0172365,.0168666, 
     *         .0164710,.0160503,.0156051,.0151361,.0146442,.0141299, 
     *         .0135791,.0130376,.0124612,.0118660,.0112526,.0106220, 
     *         .0099753,.0093134,.0086374,.0079481,.0072467,.0065344, 
     *         .0058120,.0050809,.0043420,.0035964,.0028455,.0020902, 
     *         .0013317,.0005725/
C     NGAUSS = 96
      DATA F96/.0001552,.0008178,.0020091,.0037281,.0059730,.0087413, 
     *         .0120304,.0158366,.0201558,.0249836,.0303149,.0361437, 
     *         .0424643,.0492697,.0580527,.0643058,.0725205,.0811883, 
     *         .0902998,.0998456,.1098155,.1201988,.1309847,.1421616, 
     *         .1537178,.1656409,.1779183,.1905371,.2034838,.2167448, 
     *         .2303059,.2441529,.2582710,.2726453,.2872605,.3021012, 
     *         .3171516,.3323958,.3478176,.3634006,.3791284,.3949844, 
     *         .4109516,.4270132,.4431520,.4593512,.4755935,.4918616, 
     *         .5081384,.5244065,.5406488,.5568479,.5729868,.5890484, 
     *         .6050156,.6208716,.6365994,.6521825,.6676043,.6828485, 
     *         .6978988,.7127395,.7273547,.7417290,.7558471,.7696941, 
     *         .7832552,.7965162,.8094629,.8220817,.8343592,.8462822, 
     *         .8578384,.8690153,.8798012,.8901845,.9001544,.9097002, 
     *         .9188117,.9274795,.9356943,.9419473,.9507303,.9575357, 
     *         .9638563,.9696851,.9750164,.9798442,.9841634,.9879696, 
     *         .9912587,.9940271,.9962720,.9979909,.9991822,.9998448/
      DATA W96/.0003984,.0009270,.0014553,.0019823,.0025071,.0030293, 
     *         .0035483,.0040634,.0045743,.0050804,.0055811,.0060758,
     *         .0065641,.0070455,.0075193,.0079853,.0084428,.0088913, 
     *         .0093303,.0097595,.0101784,.0105864,.0109833,.0113686, 
     *         .0117417,.0121024,.0124503,.0127850,.0131061,.0134135, 
     *         .0137065,.0139850,.0142487,.0144973,.0147306,.0149482, 
     *         .0151500,.0153357,.0155051,.0156582,.0157946,.0159144, 
     *         .0160173,.0161031,.0161719,.0162236,.0162580,.0162753, 
     *         .0162753,.0162580,.0162236,.0161719,.0161031,.0160173, 
     *         .0159144,.0157946,.0156582,.0155051,.0153357,.0151500, 
     *         .0149482,.0147306,.0144973,.0142487,.0139850,.0137065, 
     *         .0134135,.0131061,.0127850,.0124503,.0121024,.0117417, 
     *         .0113686,.0109833,.0105864,.0101784,.0097595,.0093303, 
     *         .0088913,.0084428,.0079853,.0075193,.0070455,.0065641, 
     *         .0060758,.0055811,.0050804,.0045743,.0040634,.0035483, 
     *         .0030293,.0025071,.0019823,.0014553,.0009270,.0003984/
      DATA NALLOW/2,3,4,5,6,7,8,9,10,12,16,20,24,32,40,48,64,80,96/
C-    ERROR CONDITIONS REQUIRING RETURN
      DO 200 J=1,96
      FGAUSS(J) = 0.
      WGAUSS(J) = 0.
  200 CONTINUE
      IF (NGP.LT.2 .OR. NGP.GT.96) RETURN
      NOK = 0
      DO 210 J=1,19
      IF( NGP .EQ. NALLOW(J) ) NOK = 1
  210 CONTINUE
      IF( NOK .NE. 1 ) RETURN
C-    FIND THE PROPER GAUSS POINT TABLE
      DO 220 J=1,19
      IF( NGP .EQ. NALLOW(J) ) NGPNT = J
  220 CONTINUE
      GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19)NGPNT
C-
    1 DO 1000 J=1,2
      FGAUSS(J) = F02(J)
      WGAUSS(J) = W02(J)
 1000 CONTINUE
      RETURN
C-
    2 DO 2000 J=1,3
      FGAUSS(J) = F03(J)
      WGAUSS(J) = W03(J)
 2000 CONTINUE
      RETURN
C-
    3 DO 3000 J=1,4
      FGAUSS(J) = F04(J)
      WGAUSS(J) = W04(J)
 3000 CONTINUE
      RETURN
C-
    4 DO 4000 J=1,5
      FGAUSS(J) = F05(J)
      WGAUSS(J) = W05(J)
 4000 CONTINUE
      RETURN
C-
    5 DO 5000 J=1,6
      FGAUSS(J) = F06(J)
      WGAUSS(J) = W06(J)
 5000 CONTINUE
      RETURN
C-
    6 DO 6000 J=1,7
      FGAUSS(J) = F07(J)
      WGAUSS(J) = W07(J)
 6000 CONTINUE
      RETURN
C-
    7 DO 7000 J=1,8
      FGAUSS(J) = F08(J)
      WGAUSS(J) = W08(J)
 7000 CONTINUE
      RETURN
C-
    8 DO 8000 J=1,9
      FGAUSS(J) = F09(J)
      WGAUSS(J) = W09(J)
 8000 CONTINUE
      RETURN
C-
    9 DO 9000 J=1,10
      FGAUSS(J) = F10(J)
      WGAUSS(J) = W10(J)
 9000 CONTINUE
      RETURN
C-
   10 DO 10000 J=1,12
      FGAUSS(J) = F12(J)
      WGAUSS(J) = W12(J)
10000 CONTINUE
      RETURN
C-
   11 DO 11000 J=1,16
      FGAUSS(J) = F16(J)
      WGAUSS(J) = W16(J)
11000 CONTINUE
      RETURN
C-
   12 DO 12000 J=1,20
      FGAUSS(J) = F20(J)
      WGAUSS(J) = W20(J)
12000 CONTINUE
      RETURN
C-
   13 DO 13000 J=1,24
      FGAUSS(J) = F24(J)
      WGAUSS(J) = W24(J)
13000 CONTINUE
      RETURN
C-
   14 DO 14000 J=1,32
      FGAUSS(J) = F32(J)
      WGAUSS(J) = W32(J)
14000 CONTINUE
      RETURN
C-
   15 DO 15000 J=1,40
      FGAUSS(J) = F40(J)
      WGAUSS(J) = W40(J)
15000 CONTINUE
      RETURN
C-
   16 DO 16000 J=1,48
      FGAUSS(J) = F48(J)
      WGAUSS(J) = W48(J)
16000 CONTINUE
      RETURN
C-
   17 DO 17000 J=1,64
      FGAUSS(J) = F64(J)
      WGAUSS(J) = W64(J)
17000 CONTINUE
      RETURN
C-
   18 DO 18000 J=1,80
      FGAUSS(J) = F80(J)
      WGAUSS(J) = W80(J)
18000 CONTINUE
      RETURN
C-
   19 DO 19000 J=1,96
      FGAUSS(J) = F96(J)
      WGAUSS(J) = W96(J)
19000 CONTINUE
      RETURN
C-
      END
C.....
C  ..           SUBROUTINE GRID2
C.....
C  .. GENERATE THE GAUSSIAN GRID POINTS AND WEIGHTS.
C  .. THE POINTS ARE SPECIFIED BY THE COMPONENTS OF
C  .. VECTORS IN THE [REAL] FRAME AND THE WEIGHTS
C  .. ARE PROPORTIONAL TO THE VOLUME ELEMENTS SUR-
C  .. ROUNDING THE POINTS.
C  .. CALCULATE COORDINATES OF POINTS W/R/T CENTER OF MASS
C  ..       PABC(I,3) W/R/T ABSORPTION ORIGIN OA
C  ..       RABC(I,3) W/R/T CENTER OF MASS
C.....
      SUBROUTINE GRID2
      DIMENSION NGP(3), D1(3),D2(3),D3(3)
      COMMON /INPT08/ NGP, D1,D2,D3
      DIMENSION FS(50,5)
      COMMON /PLANES/ NF, FS
      PARAMETER (NGMAX=32**3)
      DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3)
      COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC
      DIMENSION F1(4), F2(4), F3(4), POR(3)
      DIMENSION U1(3), U2(3), U3(3)
      DIMENSION GF1(96),GW1(96), GF2(96),GW2(96), GF3(96),GW3(96)
C-
      NPTS = 0
      VXTL = 0.0
      DO 10 I=1,3
      VCOM(I) = 0.
   10 CONTINUE
      DO 30 I=1,NGMAX
      WABC(I) = 0.
      DO 20 J=1,3
      PABC(I,J) = 0.
      RABC(I,J) = 0.
   20 CONTINUE
   30 CONTINUE
C-
      IF (NF .LE. 4)  RETURN
      NP1 = NGP(1)
      NP2 = NGP(2)
      NP3 = NGP(3)
      CALL GAUTAB (NP1, GF1,GW1)
      CALL GAUTAB (NP2, GF2,GW2)
      CALL GAUTAB (NP3, GF3,GW3)
      CALL VSTAR (VSTR)
      CALL UNIVS (D1, U1)
      CALL UNIVS (D2, U2)
      CALL UNIVS (D3, U3)
C-
      CALL VERT1
      CALL MINMAX (U1, D1MIN,D1MAX)
      D1DEL = D1MAX - D1MIN
C-
      DO 1000 N1 = 1,NP1
      EPS1 = D1MIN + GF1(N1)*D1DEL
      PIE1 = GW1(N1)*D1DEL
      CALL FLIPS (U1,EPS1, F1)
      CALL VERT2 (F1)
      CALL MINMAX (U2, D2MIN,D2MAX)
      D2DEL = D2MAX - D2MIN
C-
      DO 900 N2 = 1,NP2
      EPS2 = D2MIN + GF2(N2)*D2DEL
      PIE2 = GW2(N2)*D2DEL
      CALL FLIPS (U2,EPS2, F2)
      CALL VERT3 (F1,F2)
      CALL MINMAX (U3, D3MIN,D3MAX)
      D3DEL = D3MAX - D3MIN
C-
      DO 800 N3 = 1,NP3
      EPS3 = D3MIN + GF3(N3)*D3DEL
      PIE3 = GW3(N3)*D3DEL
      CALL FLIPS (U3,EPS3, F3)
      CALL VERT4 (F1,F2,F3, POR)
      VOL = PIE1*PIE2*PIE3*VSTR
      NPTS = NPTS + 1
      PABC(NPTS,1) = POR(1)
      PABC(NPTS,2) = POR(2)
      PABC(NPTS,3) = POR(3)
      WABC(NPTS)   = VOL
      VXTL = VXTL + VOL
      DO 700 I = 1,3
      VCOM(I) = VCOM(I) + VOL*POR(I)
  700 CONTINUE
  800 CONTINUE
  900 CONTINUE
 1000 CONTINUE
C-
      DO 1100 I =1,3
      VCOM(I) = VCOM(I)/VXTL
 1100 CONTINUE
      DO 1200 I =1,NPTS
      WABC(I) = WABC(I)/VXTL
 1200 CONTINUE
      DO 1300 I =1,NPTS
      DO 1300 J =1,3
      RABC(I,J) = PABC(I,J) - VCOM(J)
 1300 CONTINUE
C-
      CALL VERT1
      RETURN
      END
C.....
C  ..           SUBROUTINE GRID3
C.....
C  .. OUTPUT INFORMATION RELEVANT TO GRID SECTIONING
C.....
      SUBROUTINE GRID3
      DIMENSION NGP(3), D1(3),D2(3),D3(3)
      COMMON /INPT08/ NGP, D1,D2,D3
      PARAMETER (NGMAX=32**3)
      DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3)
      COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC
      COMMON /CRYSTL/ TMOLW, RHO, EMUX
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      DIMENSION  U1(3),U2(3),U3(3)
      DATA CONS/57.2957795/
C-
      CALL UNIVS (D1, U1)
      CALL UNIVS (D2, U2)
      CALL UNIVS (D3, U3)
      CALL VERT1
      CALL MINMAX (U1, DMIN1,DMAX1)
      CALL MINMAX (U2, DMIN2,DMAX2)
      CALL MINMAX (U3, DMIN3,DMAX3)
      CALL SDS1 (U1,U2, DOT12)
      CALL SDS1 (U1,U3, DOT13)
      CALL SDS1 (U2,U3, DOT23)
      DOT12 = CONS*ACOS(DOT12)
      DOT13 = CONS*ACOS(DOT13)
      DOT23 = CONS*ACOS(DOT23)
      CMASS = VXTL*RHO
C-
      WRITE(NOU,8001)
 8001 FORMAT(/1X,'GRID SECTIONING INDICES AND EXTREMES')
      WRITE(NOU,8002)
 8002 FORMAT(1X,5X,'INDICES',8X,'DMIN',5X,'DMAX',3X,'#PTS')
      WRITE(NOU,8003) D1,DMIN1,DMAX1,NGP(1)
      WRITE(NOU,8003) D2,DMIN2,DMAX2,NGP(2)
      WRITE(NOU,8003) D3,DMIN3,DMAX3,NGP(3)
 8003 FORMAT(1X,3F5.0,4X,F6.3,3X,F6.3,I5)
      WRITE(NOU,8004) DOT12,DOT13,DOT23
 8004 FORMAT(/1X,'ANGLE D1^D2: ',F6.2,/,
     *        1X,'ANGLE D1^D3: ',F6.2,/,
     *        1X,'ANGLE D2^D3: ',F6.2   )
      WRITE(NOU,8005) VXTL
 8005 FORMAT(/1X,'CRYSTAL VOLUME(MM**3): ',E12.5)
      WRITE(NOU,8006) CMASS
 8006 FORMAT(1X,'CRYSTAL MASS(MG):      ',E12.5)
      WRITE(NOU,8007) VCOM
 8007 FORMAT(1X,'C.O.M. [REAL] VECTOR: ',3E12.4)
      RETURN
      END
C.....
C  ..           SUBROUTINE INNIE
C.....
C  .. ANSWERS THE AGE-OLD QUESTION:  IS THE
C  .. GENERATED POINT WITHIN THE CRYSTAL?
C  .. NIN = +1  YES IT IS
C  .. NIN = -1  NO IT IS NOT
C.....
      SUBROUTINE INNIE (POR, NIN)
      DIMENSION FS(50,5)
      COMMON /PLANES/ NF, FS
      DIMENSION POR(3)
      DATA EPSIL/0.0001/
C-
      NIN = +1
      DO 100 M=1,NF
      DOT = POR(1)*FS(M,1)
     *     +POR(2)*FS(M,2)
     *     +POR(3)*FS(M,3)
     *     -EPSIL
      IF( DOT.GE.FS(M,4) ) NIN = -1
  100 CONTINUE
      RETURN
      END
C.....
C  ..           SUBROUTINE KAPP2
C.....
C  .. CALCULATE THE U(PHI,TAU,SGM) MATRIX AND CONVERT
C  .. THE POINTS PABC INTO THEIR [TUBE] COMPONENTS
C.....
      SUBROUTINE KAPP2
      CHARACTER*4 WAVE
      COMMON /INPT04/ WAVE
      COMMON /LAMBDA/ FLAM
      DIMENSION ORMAT(3,3)
      COMMON /INPT09/ ORMAT, NDIFF
      COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP
      PARAMETER (NGMAX=32**3)
      DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3)
      COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
      DIMENSION UBF(3,3), UBI(3,3), UBGF(3,3), UBGI(3,3)
      COMMON /MATRIX/ UBF, UBI, UBGF, UBGI
C-
      DIMENSION UPTSF(3,3), UPTSI(3,3), UBIG(3,3), DTUB(NGMAX,3)
      COMMON /KAPCOR/ EMUK, RAD2, RPT2, UPTSF, UPTSI, UBIG, DTUB
C-
      DIMENSION UPHI1(3,3), UPHI2(3,3), UTS(3,3)
      DIMENSION VKPH(3), VKTB(3), VNTB(3), P(3), D(3)
      DIMENSION ETAB(3,5)
      DATA CON/57.2957795/
      CHARACTER*4 CUKA,CUKB, MOKA,MOKB, AGKA,AGKB
      DATA CUKA/'CUKA'/, CUKB/'CUKB'/
      DATA MOKA/'MOKA'/, MOKB/'MOKB'/
      DATA AGKA/'AGKA'/, AGKB/'AGKB'/
C-
C-    RETURN IF NO CAPILLARY
      IF(NKAP.EQ.0) RETURN
C-
C-    DIMENSIONS OF THE CAPPILARY
      RAD2 = (0.5*DKAP)**2.
      RPT2 = (0.5*DKAP + TKAP)**2.
C-
C-    LOOK UP THE CAPILLARY LINEAR ABSORPTION COEFFICIENT
      IF(WAVE.EQ.CUKA .OR. WAVE.EQ.CUKB) NRAD = 1
      IF(WAVE.EQ.MOKA .OR. WAVE.EQ.MOKB) NRAD = 2
      IF(WAVE.EQ.AGKA .OR. WAVE.EQ.AGKB) NRAD = 3
C-    TABLE OF EMUK CALCULATED BY HAND
C-    ETAB(NRADIATION,NKAPPILARY) IN RECIPROCAL MM
C     NRAD = 1  CU KA    NKAP = 1  "GLASKAPILLAREN"
C            2  MO KA           2  "PYREX 7740"
C            3  AG KA           3  "AR GLASS"
C                               4  FUSED SILICA
C                               5  LINDEMANN GLASS
      ETAB(1,1) = 10.870
      ETAB(1,2) =  7.600
      ETAB(1,3) = 13.298
      ETAB(1,4) =  8.021
      ETAB(1,5) =  1.484
      ETAB(2,1) =  1.058
      ETAB(2,2) =  0.733
      ETAB(2,3) =  1.360
      ETAB(2,4) =  0.807
      ETAB(2,5) =  0.161
      ETAB(3,1) =  0.538
      ETAB(3,2) =  0.375
      ETAB(3,3) =  0.700
      ETAB(3,4) =  0.411
      ETAB(3,5) =  0.094
      EMUK = ETAB(NRAD,NKAP)
C-    NICOLET P3/F UPHI MATRIX
      UPHI1(1,1) = +1.
      UPHI1(1,2) = +0.
      UPHI1(1,3) = +0.
      UPHI1(2,1) = +0.
      UPHI1(2,2) = +1.
      UPHI1(2,3) = +0.
      UPHI1(3,1) = +0.
      UPHI1(3,2) = +0.
      UPHI1(3,3) = +1.
C-    NONIUS CAD-4 AND BUSING-LEVY UPHI MATRIX
      UPHI2(1,1) = +0.
      UPHI2(1,2) = -1.
      UPHI2(1,3) = +0.
      UPHI2(2,1) = +1.
      UPHI2(2,2) = +0.
      UPHI2(2,3) = +0.
      UPHI2(3,1) = +0.
      UPHI2(3,2) = +0.
      UPHI2(3,3) = +1.
C-    BUILD THE UTS MATRIX EXPLICITLY
      CTAU = COS (TAU/CON)
      STAU = SIN (TAU/CON)
      CSGM = COS (SGM/CON)
      SSGM = SIN (SGM/CON)
      UTS(1,1) = +CTAU
      UTS(1,2) = +STAU*SSGM
      UTS(1,3) = -STAU*CSGM
      UTS(2,1) = +0.
      UTS(2,2) = +CSGM
      UTS(2,3) = +SSGM
      UTS(3,1) = +STAU
      UTS(3,2) = -CTAU*SSGM
      UTS(3,3) = +CTAU*CSGM
C-    GET UPTSF AND ITS INVERSE UPTSI
      IF(NDIFF.EQ.1) CALL MTM1 (UPHI1,UTS, UPTSF)
      IF(NDIFF.EQ.2) CALL MTM1 (UPHI2,UTS, UPTSF)
      IF(NDIFF.EQ.3) CALL MTM1 (UPHI2,UTS, UPTSF)
      CALL AMI3 (UPTSF,UPTSI,UDET)
      IF(UDET.EQ.0.0)  RETURN
C-    GET UBIG
      CALL MTM1 (UPTSI,UBGF, UBIG)
C-    CONVERT KAPILLARY DISPLACEMENT FROM [PHI] TO [TUBE]
      VKPH(1) = XDSP
      VKPH(2) = YDSP
      VKPH(3) = 0.0
      CALL MTV1 (UPTSI,VKPH, VKTB)
C-    CONVERT CENTER-OF-MASS VECTOR FROM [REAL] TO [TUBE]
      CALL MTV1 (UBIG,VCOM, VNTB)
C-    FORM THE VECTOR { -N -K }
      VNK1 = -VNTB(1) -VKTB(1)
      VNK2 = -VNTB(2) -VKTB(2)
      VNK3 = -VNTB(3) -VKTB(3)
C-    CONVERT THE PABC VECTORS FROM [REAL] TO [TUBE]
C-    AND FORM THE DTUB VECTORS  P - N - K
      DO 1000 N = 1,NPTS
      P(1) = PABC(N,1)
      P(2) = PABC(N,2)
      P(3) = PABC(N,3)
      CALL MTV1 (UBIG,P, D)
      DTUB(N,1) = D(1) + VNK1
      DTUB(N,2) = D(2) + VNK2
      DTUB(N,3) = D(3) + VNK3
 1000 CONTINUE
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE KAPP3
C.....
C  .. CALCULATE AND OUTPUT THE [TUBE] COORDINATES OF THE
C  .. CRYSTAL VERTICES.  CHECK TO SEE THAT ALL THE VERTICES 
C  .. ARE WITHIN THE TUBE RADIUS.
C.....
      SUBROUTINE KAPP3
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
      DIMENSION VR(150,3), MM(150,3)
      COMMON /VERTEX/ NV, VR, MM
      PARAMETER (NGMAX=32**3)
      DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3)
      COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC
      DIMENSION UPTSF(3,3), UPTSI(3,3), UBIG(3,3), DTUB(NGMAX,3)
      COMMON /KAPCOR/ EMUK, RAD2, RPT2, UPTSF, UPTSI, UBIG, DTUB
C-
      DIMENSION VKPH(3), VKTB(3), VNTB(3), P(3), D(3)
      DIMENSION VTUB(150,3)
C-    RETURN IF NO KAPILLARY
      IF(NKAP.EQ.0) RETURN
C-    CONVERT KAPILLARY DISPLACEMENT FROM [PHI] TO [TUBE]
      VKPH(1) = XDSP
      VKPH(2) = YDSP
      VKPH(3) = 0.0
      CALL MTV1 (UPTSI,VKPH, VKTB)
C-    CONVERT CENTER-OF-MASS VECTOR FROM [REAL] TO [TUBE]
      CALL MTV1 (UBIG,VCOM, VNTB)
C-    FORM THE VECTOR { -N -K }
      VNK1 = -VNTB(1) -VKTB(1)
      VNK2 = -VNTB(2) -VKTB(2)
      VNK3 = -VNTB(3) -VKTB(3)
C-    CONVERT THE VERTEX VECTORS FROM [REAL] TO [TUBE]
      DO 1000 N = 1,NV
      P(1) = VR(N,1)
      P(2) = VR(N,2)
      P(3) = VR(N,3)
      CALL MTV1 (UBIG,P, D)
      VTUB(N,1) = D(1) + VNK1
      VTUB(N,2) = D(2) + VNK2
      VTUB(N,3) = D(3) + VNK3
 1000 CONTINUE
C-
      IF(NKAP.EQ.1) WRITE(NOU,9101)
      IF(NKAP.EQ.2) WRITE(NOU,9102)
      IF(NKAP.EQ.3) WRITE(NOU,9103)
      IF(NKAP.EQ.4) WRITE(NOU,9104)
      IF(NKAP.EQ.5) WRITE(NOU,9105)
 9101 FORMAT(/1X,'CRYSTAL IN A "GLASKAPILLAREN" CAPILLARY')
 9102 FORMAT(/1X,'CRYSTAL IN A "PYREX 7740" CAPILLARY')
 9103 FORMAT(/1X,'CRYSTAL IN A "AR GLASS" CAPILLARY')
 9104 FORMAT(/1X,'CRYSTAL IN A FUSED SILICA CAPILLARY')
 9105 FORMAT(/1X,'CRYSTAL IN A LINDEMANN GLASS CAPILLARY')
      WRITE(NOU,9004) EMUK
 9004 FORMAT(1X,'CAPILLARY LINEAR ABSORPTION COEFF(MM): ',F6.3)
      WRITE(NOU,9005) DKAP, TKAP
 9005 FORMAT(1X,'CAPILLARY DIAMETER: ',F6.3,/,
     *       1X,'CAPILLARY THICKNESS:',F6.3   )
      WRITE(NOU,9006) SGM, TAU
 9006 FORMAT(1X,'TOP ARC(SGM): ',F7.3,/,
     *       1X,'LOW ARC(TAU): ',F7.3   )
      WRITE(NOU,9007) XDSP, YDSP
 9007 FORMAT(1X,'X[PHI] DISPLACEMENT(MM): ',F8.5,/,
     *       1X,'Y[PHI] DISPLACEMENT(MM): ',F8.5   )
C-    WRITE OUT THE DISTANCES FROM THE [TUBE] ORIGIN
C-    TO THE VERTICES IN THE CAPILLARY X,Y PLANE
      WRITE(NOU,9008)
 9008 FORMAT(/1X,'DISTANCES OF VERTICES FROM CAPILLARY WALL IN MM')
      WRITE(NOU,9009)
 9009 FORMAT(1X,'NVERT',4X,'DIST',20X,'[TUBE] COMPONENTS')
      NWARN = 0
      RAD1 = SQRT(RAD2)
      DO 2000 N=1,NV
      DIST = SQRT(VTUB(N,1)**2.+VTUB(N,2)**2.)
      IF(DIST.GE.RAD1) NWARN = 1
      DIST = RAD1 - DIST
      WRITE(NOU,9010) N, DIST, (VTUB(N,I),I=1,3)
 9010 FORMAT(1X,I4,2X,F9.6,5X,3E15.6)
 2000 CONTINUE
C-
      IF(NWARN.NE.0) WRITE(NOU,9011)
 9011 FORMAT(/1X,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!',/,
     *        1X,'ONE OR MORE VERTICES CALCULATED NOT IN TUBE',/,
     *        1X,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'   )
      RETURN
      END
C.....
C  ..           SUBROUTINE LIMS1
C.....
C  .. CALCULATE AND OUTPUT THE LIMITS OF THE
C  .. CRYSTAL IN VARIOUS PRINCIPAL DIRECTIONS
C.....
      SUBROUTINE LIMS1
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      DIMENSION G1(3),U1(3)
      DIMENSION GG(13,3)
C-
      GG( 1,1) = +1.
      GG( 1,2) = +0.
      GG( 1,3) = +0.
C-
      GG( 2,1) = +0.
      GG( 2,2) = +1.
      GG( 2,3) = +0.
C-
      GG( 3,1) = +0.
      GG( 3,2) = +0.
      GG( 3,3) = +1.
C-
      GG( 4,1) = +1.
      GG( 4,2) = +1.
      GG( 4,3) = +0.
C-
      GG( 5,1) = +1.
      GG( 5,2) = +0.
      GG( 5,3) = +1.
C-
      GG( 6,1) = +0.
      GG( 6,2) = +1.
      GG( 6,3) = +1.
C-
      GG( 7,1) = +1.
      GG( 7,2) = -1.
      GG( 7,3) = +0.
C-
      GG( 8,1) = +1.
      GG( 8,2) = +0.
      GG( 8,3) = -1.
C-
      GG( 9,1) = +0.
      GG( 9,2) = +1.
      GG( 9,3) = -1.
C-
      GG(10,1) = +1.
      GG(10,2) = +1.
      GG(10,3) = +1.
C-
      GG(11,1) = +1.
      GG(11,2) = +1.
      GG(11,3) = -1.
C-
      GG(12,1) = +1.
      GG(12,2) = -1.
      GG(12,3) = +1.
C-
      GG(13,1) = +1.
      GG(13,2) = -1.
      GG(13,3) = -1.
C-
      WRITE(NOU,9001)
 9001 FORMAT(/1X,'LIMITS OF CRYSTAL IN MILLIMETERS')
      WRITE(NOU,9002)
 9002 FORMAT(1X,2X,'DIRECTION',12X,'DMIN',9X,'DMAX',8X,'DELT')
      DO 100 I=1,13
      G1(1) = GG(I,1)
      G1(2) = GG(I,2)
      G1(3) = GG(I,3)
      CALL UNIVS (G1, U1)
      CALL MINMAX (U1, DMIN, DMAX)
      DDEL = ABS(DMAX-DMIN)
      WRITE(NOU,9003) G1,DMIN,DMAX,DDEL
 9003 FORMAT(1X,3F4.0,5X,3F12.6)
  100 CONTINUE
      RETURN
      END
C.....
C  ..           SUBROUTINE MATX2
C.....
C  .. CALCULATE THE UBF UBI UBGF AND UBGI MATRICES
C.....
      SUBROUTINE MATX2
      DIMENSION ORMAT(3,3)
      COMMON /INPT09/ ORMAT, NDIFF
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
      DIMENSION UBF(3,3), UBI(3,3), UBGF(3,3), UBGI(3,3)
      COMMON /MATRIX/ UBF, UBI, UBGF, UBGI
C-
      DO 1200 I=1,3
      DO 1100 J=1,3
      UBF(I,J) = ORMAT(I,J)
 1100 CONTINUE
 1200 CONTINUE
      CALL AMI3 (UBF, UBI,UDET)
      CALL MTM1 (UBF,GR, UBGF)
      CALL AMI3 (UBGF, UBGI,UDET)
      RETURN
      END
C.....
C  ..           SUBROUTINE MATX3
C  .. OUTPUT THE DATA ON THE ORIENTATION MATRIX AND DIFFRACTOMETER TYPE
C.....
      SUBROUTINE MATX3
      DIMENSION ORMAT(3,3)
      COMMON /INPT09/ ORMAT, NDIFF
      DIMENSION UBF(3,3), UBI(3,3), UBGF(3,3), UBGI(3,3)
      COMMON /MATRIX/ UBF, UBI, UBGF, UBGI
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      IF(NDIFF.EQ.1)  GO TO 100
      IF(NDIFF.EQ.2)  GO TO 200
      IF(NDIFF.EQ.3)  GO TO 300
  100 WRITE(NOU,8001)
      GO TO 500
  200 WRITE(NOU,8002)
      GO TO 500
  300 WRITE(NOU,8003)
      GO TO 500
 8001 FORMAT(/1X,'NICOLET P3/F DIFFRACTOMETER')
 8002 FORMAT(/1X,'NONIUS CAD-4 DIFFRACTOMETER')
 8003 FORMAT(/1X,'BUSING-LEVY DIFFRACTOMETER')
  500 CONTINUE
C-
      WRITE(NOU,9001)
 9001 FORMAT(1X,'ORIENTATION MATRIX UBF(I,J)')
      DO 1000 I=1,3
      WRITE(NOU,9002) (UBF(I,J),J=1,3)
 1000 CONTINUE
 9002 FORMAT(1X,3E15.6)
      RETURN
      END
C.....
C  ..           SUBROUTINE MINMAX
C.....
C  .. GIVEN A LIST OF COMPONENTS OF VERTEX VECTORS
C  .. IN THE [REAL] FRAME SEARCH FOR THE MINIMUM
C  .. AND MAXIMUM OF THE DOT-PRODUCT U*-DOT-V
C  .. GIVEN UNIT VECTOR U* IN THE [STAR] FRAME
C.....
      SUBROUTINE MINMAX(UVECT, DMIN,DMAX)
      DIMENSION VR(150,3), MM(150,3)
      COMMON /VERTEX/ NV, VR, MM
      DIMENSION UVECT(3)
C-
      U1 = UVECT(1)
      U2 = UVECT(2)
      U3 = UVECT(3)
      DMIN = +1E10
      DMAX = -1E10
      DO 10 J=1,NV
      V1 = VR(J,1)
      V2 = VR(J,2)
      V3 = VR(J,3)
      DOT = U1*V1 + U2*V2 + U3*V3
      IF (DOT.LE.DMIN) DMIN = DOT
      IF (DOT.GE.DMAX) DMAX = DOT
   10 CONTINUE
      RETURN
      END
C.....
C  ..           SUBROUTINE MTM1
C.....
C  .. MULTIPLY MATRIX B BY MATRIX A
C  .. TO GET MATRIX C   C = A B
C.....
      SUBROUTINE MTM1 (A,B, C)
      DIMENSION A(3,3), B(3,3), C(3,3)
      DO 10 I=1,3
      DO 5  J=1,3
      C(I,J) = 0.
    5 CONTINUE
   10 CONTINUE
C-
      DO 30 I=1,3
      DO 25 J=1,3
      DO 20 N=1,3
      C(I,J) = C(I,J) + A(I,N)*B(N,J)
   20 CONTINUE
   25 CONTINUE
   30 CONTINUE
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE MTV1
C.....
C  .. MULTIPLY VECTOR -R- BY MATRIX -G- TO GET 
C  .. VECTOR -S-
C.....
      SUBROUTINE MTV1(G, R, S)
      DIMENSION G(3,3), R(3), S(3)
      S(1) = G(1,1)*R(1) + G(1,2)*R(2) + G(1,3)*R(3)
      S(2) = G(2,1)*R(1) + G(2,2)*R(2) + G(2,3)*R(3)
      S(3) = G(3,1)*R(1) + G(3,2)*R(2) + G(3,3)*R(3)
      RETURN
      END
C.....
C  ..           SUBROUTINE OUTP1
C.....
C  .. PRINT OUT A SUMMARY OF INPUT AND PRIMARY OUTPUT DATA
C.....
      SUBROUTINE OUTP1
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      OPEN(UNIT=NLP,FILE='absorb.lp',STATUS='UNKNOWN')
      DO I=1,1E9
        READ (NLP,'(A1)',END=1) DUMMY
      END DO
 1    WRITE (NLP,'(/1X)')
      I=NOU
      NOU=NLP
      CALL TITL2
      CALL FILS2
      CALL BPRO2
      CALL CELL2
      CALL CELL3
      CALL FACE2
      CALL FACE3
      CALL VERT1
      CALL VERT5
      CALL LIMS1
      CALL EDGE1
      CALL MATX2
      CALL MATX3
      CALL EMUX2
      CALL EMUX3
      CALL GRID2
      CALL GRID3
      CALL KAPP2
      CALL KAPP3
      CALL SOLN4
      CLOSE(UNIT=NLP,STATUS='KEEP')
      NOU=I
      RETURN
      END
C.....
C  ..           SUBROUTINE P3RMT
C.....
C  .. GET THE BUSING-LEVY ROTATION MATRIX FOR THE
C  .. NICOLET P3/F DIFFRACTOMETER
C  .. ANGLES ARE: TWO THETA, OMEGA, PHI, CHI
C.....
      SUBROUTINE P3RMT(ANGLS, RMTF,ST,CT)
      DIMENSION ANGLS(4), RMTF(3,3)
      DATA CONS/57.29577951/
C-
      TTH = ANGLS(1)/CONS
      OME = ANGLS(2)/CONS
      PHI = ANGLS(3)/CONS
      CHI = ANGLS(4)/CONS
      THT = 0.5*TTH
      OMT = OME - THT
      ST = SIN(THT)
      CT = COS(THT)
      SO = SIN(OMT)
      CO = COS(OMT)
      SP = SIN(PHI)
      CP = COS(PHI)
      SC = SIN(CHI)
      CC = COS(CHI)
C-    EXPLICIT GENERATION OF ROTATION MATRIX
      RMTF(1,1) = +CO*CP + SO*CC*SP
      RMTF(1,2) = +CO*SP - SO*CC*CP
      RMTF(1,3) = -SO*SC
      RMTF(2,1) = +SO*CP - CO*CC*SP
      RMTF(2,2) = +SO*SP + CO*CC*CP
      RMTF(2,3) = +CO*SC
      RMTF(3,1) = +SC*SP
      RMTF(3,2) = -SC*CP
      RMTF(3,3) = +CC
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE P3ROT
C.....
C  .. GET THE [PHI] FRAME COMPONENTS OF THE UNIT
C  .. VECTORS FOR THE NICOLET P3/F DIFFRACTOMETER
C  .. -SI- REVERSE INCIDENT BEAM
C  .. -SD- DIFFRACTED BEAM
C  .. -PM- MONOCHROMATOR HORIZONTAL PROFILE
C  .. ANGLES ARE: TWO THETA, OMEGA, PHI, CHI
C.....
      SUBROUTINE P3ROT(ANGLS, SIP,SDP,PMP)
      DIMENSION ANGLS(4), SIP(3), SDP(3), PMP(3)
      DIMENSION RMTF(3,3), RMTI(3,3), SDT(3), SIT(3), PMT(3)
C-
      CALL P3RMT (ANGLS, RMTF,ST,CT)
C-
C-    TRANSPOSE OF ROTATION MATRIX EQUALS ITS INVERSE
      DO 100 I=1,3
      DO 50  J=1,3
      RMTI(I,J) = RMTF(J,I)
   50 CONTINUE
  100 CONTINUE
C-
      SIT(1) = -CT
      SIT(2) = +ST
      SIT(3) = 0.0
      SDT(1) = +CT
      SDT(2) = +ST
      SDT(3) = 0.0
      PMT(1) = +CT
      PMT(2) = +ST
      PMT(3) = +0.0
C-
      CALL MTV1(RMTI,SIT, SIP)
      CALL MTV1(RMTI,SDT, SDP)
      CALL MTV1(RMTI,PMT, PMP)
      RETURN
      END
C.....
C  ..           SUBROUTINE PAGE1
C.....
C  .. READ IN FROM FILE NAB THE ARCHIVAL QUALITY DATA
C.....
      SUBROUTINE PAGE1
      CHARACTER TITLE*80, ATIME*8, ADATE*9
      CHARACTER FIL01*80, FIL02*80
      DIMENSION CELL(6)
      CHARACTER*4 WAVE
      CHARACTER*2 ATOMS(10)
      DIMENSION STOIC(10)
      DIMENSION FACES(50,4)
      DIMENSION NGP(3), D1(3),D2(3),D3(3)
      DIMENSION ORMAT(3,3)
      DIMENSION NTRP(50)
      COMMON /INPT01/ TITLE, ATIME, ADATE
      COMMON /INPT02/ FIL01, FIL02, IPATH
      COMMON /INPT03/ CELL
      COMMON /INPT04/ WAVE
      COMMON /LAMBDA/ FLAM
      COMMON /INPT05/ ATOMS
      COMMON /INPT06/ ZCELL, STOIC
      COMMON /INPT07/ FACES
      COMMON /INPT08/ NGP, D1,D2,D3
      COMMON /INPT09/ ORMAT, NDIFF
      COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP
      COMMON /INPT11/ NBPR, A0, A1, A2, A3
      COMMON /INPT12/ NSOL, EMUS, NTRP
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      COMMON /CRYSTL/ TMOLW, RHO, EMUX
C-
      OPEN(UNIT=NAB,FILE='absorb.dat',STATUS='OLD',ERR=999)
C-
      READ(NAB,7001) TITLE
      READ(NAB,7001) ATIME
      READ(NAB,7001) ADATE
      READ(NAB,7001) FIL01
      READ(NAB,7001) FIL02
      READ(NAB,7010) IPATH
      READ(NAB,7002) CELL
      READ(NAB,7001) WAVE
      READ(NAB,7002) FLAM
      READ(NAB,7003) ATOMS
      READ(NAB,7004) ZCELL
      READ(NAB,7005) STOIC
      READ(NAB,7002) EMUX
      READ(NAB,7006) ( (FACES(I,J),J=1,4), I=1,50)
      READ(NAB,7007) NGP
      READ(NAB,7008) D1,D2,D3
      READ(NAB,7009) ( (ORMAT(I,J),J=1,3),  I=1,3)
      READ(NAB,7010) NDIFF
      READ(NAB,7011) NKAP, DKAP, TKAP
      READ(NAB,7012) TAU, SGM, XDSP, YDSP
      READ(NAB,7013) NBPR, A0, A1, A2, A3
      READ(NAB,7014) NSOL, EMUS
      READ(NAB,7015) NTRP
      CLOSE(UNIT=NAB,STATUS='KEEP')
      RETURN
C-
  999 WRITE(IO6,8001)
      RETURN
C-
 7001 FORMAT(A)
 7002 FORMAT(6F10.5)
 7003 FORMAT(10A)
 7004 FORMAT(F5.2)
 7005 FORMAT(10F7.2)
 7006 FORMAT(3F10.3,F10.5,/)
 7007 FORMAT(3I5)
 7008 FORMAT(3F6.2,/)
 7009 FORMAT(3E15.7)
 7010 FORMAT(I5)
 7011 FORMAT(I5,F7.4,F7.4)
 7012 FORMAT(4F10.5)
 7013 FORMAT(I5,4F10.6)
 7014 FORMAT(I5,F10.6)
 7015 FORMAT(50I1)
C-
8001  FORMAT(1X,'ERROR OPENING THE ABSORB.DAT FILE')
      END
C.....
C  ..           SUBROUTINE PAGE2
C.....
C  .. WRITE OUT ON FILE NAB THE ARCHIVAL QUALITY DATA
C.....
      SUBROUTINE PAGE2
      CHARACTER TITLE*80, ATIME*8, ADATE*9
      CHARACTER FIL01*80, FIL02*80
      DIMENSION CELL(6)
      CHARACTER*4 WAVE
      CHARACTER*2 ATOMS(10)
      DIMENSION STOIC(10)
      DIMENSION FACES(50,4)
      DIMENSION NGP(3), D1(3),D2(3),D3(3)
      DIMENSION ORMAT(3,3)
      DIMENSION NTRP(50)
      COMMON /INPT01/ TITLE, ATIME, ADATE
      COMMON /INPT02/ FIL01, FIL02, IPATH
      COMMON /INPT03/ CELL
      COMMON /INPT04/ WAVE
      COMMON /LAMBDA/ FLAM
      COMMON /INPT05/ ATOMS
      COMMON /INPT06/ ZCELL, STOIC
      COMMON /INPT07/ FACES
      COMMON /INPT08/ NGP, D1,D2,D3
      COMMON /INPT09/ ORMAT, NDIFF
      COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP
      COMMON /INPT11/ NBPR, A0, A1, A2, A3
      COMMON /INPT12/ NSOL, EMUS, NTRP
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      COMMON /CRYSTL/ TMOLW, RHO, EMUX
      OPEN(UNIT=NAB,FILE='absorb.dat.new',STATUS='NEW')
      WRITE(NAB,7001) TITLE
      WRITE(NAB,7001) ATIME
      WRITE(NAB,7001) ADATE
      WRITE(NAB,7001) FIL01
      WRITE(NAB,7001) FIL02
      WRITE(NAB,7010) IPATH
      WRITE(NAB,7002) CELL
      WRITE(NAB,7001) WAVE
      WRITE(NAB,7002) FLAM
      WRITE(NAB,7003) ATOMS
      WRITE(NAB,7004) ZCELL
      WRITE(NAB,7005) STOIC
      WRITE(NAB,7002) EMUX
      WRITE(NAB,7006) ( (FACES(I,J),J=1,4), I=1,50)
      WRITE(NAB,7007) NGP
      WRITE(NAB,7008) D1,D2,D3
      WRITE(NAB,7009) ( (ORMAT(I,J),J=1,3),  I=1,3)
      WRITE(NAB,7010) NDIFF
      WRITE(NAB,7011) NKAP, DKAP, TKAP
      WRITE(NAB,7012) TAU, SGM, XDSP, YDSP
      WRITE(NAB,7013) NBPR, A0, A1, A2, A3
      WRITE(NAB,7014) NSOL, EMUS
      WRITE(NAB,7015) NTRP
      CLOSE(UNIT=NAB,STATUS='KEEP')
 7001 FORMAT(A)
 7002 FORMAT(6F10.5)
 7003 FORMAT(10A)
 7004 FORMAT(F5.2)
 7005 FORMAT(10F7.2)
 7006 FORMAT(3F10.3,F10.5,/)
 7007 FORMAT(3I5)
 7008 FORMAT(3F6.2,/)
 7009 FORMAT(3E15.7)
 7010 FORMAT(I5)
 7011 FORMAT(I5,F7.4,F7.4)
 7012 FORMAT(4F10.5)
 7013 FORMAT(I5,4F10.6)
 7014 FORMAT(I5,F10.6)
 7015 FORMAT(50I1)
      RETURN
      END
C.....
C  ..           SUBROUTINE PATH1
C.....
C  .. FIND THE ENTRY/EXIT FACE OF THE BEAM
C  .. WITH [REAL] SPACE COMPONENTS -SR- FOR
C  .. THE POINT WITH [REAL] COMPONENTS -P-
C  .. AND CALCULATE THE PATH LENGTH -PL-
C  .. IN MILLIMETERS
C.....
      SUBROUTINE PATH1 (SR,JPT, PL,NX)
      DIMENSION FS(50,5)
      COMMON /PLANES/ NF, FS
      PARAMETER (NGMAX=32**3)
      DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3)
      COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC
      DIMENSION SR(3)
C-
      PL = 1E12
      NX = 0
      S1 = SR(1)
      S2 = SR(2)
      S3 = SR(3)
      P1 = PABC(JPT,1)
      P2 = PABC(JPT,2)
      P3 = PABC(JPT,3)
C-
      DO 1000 N=1,NF
      PN = 0.
      U1 = FS(N,1)
      U2 = FS(N,2)
      U3 = FS(N,3)
      DV = FS(N,4)
      SU = U1*S1 + U2*S2 + U3*S3
      IF( SU.LE.0.0 )  GO TO 999
      PU = U1*P1 + U2*P2 + U3*P3
      PN = ( DV - PU ) / SU
      IF( PN.GE.PL ) GO TO 999
      PL = PN
      NX = N
  999 FS(N,5) = PN
 1000 CONTINUE
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE PATH2
C.....
C  .. CALCULATE THE PATH LENGTH THROUGH THE CAPILLARY
C  .. WINN IS DISTANCE IN MM TO INNER WALL
C  .. WOUT IS DISTANCE IN MM TO OUTER WALL
C  .. PCAP IS DISTANCE IN MM THROUGH CAPILLARY
C.....
      SUBROUTINE PATH2 (S,JPT, WINN,WOUT,PCAP)
      PARAMETER (NGMAX=32**3)
      DIMENSION UPTSF(3,3), UPTSI(3,3), UBIG(3,3), DTUB(NGMAX,3)
      COMMON /KAPCOR/ EMUK, RAD2, RPT2, UPTSF, UPTSI, UBIG, DTUB
      DIMENSION S(3)
      S1 = S(1)
      S2 = S(2)
      D1 = DTUB(JPT,1)
      D2 = DTUB(JPT,2)
      PT1 = D1*S1 + D2*S2
      PT2 = S1*S1 + S2*S2
      PT3 = D1*D1 + D2*D2
      WINN = (-PT1 + SQRT(PT1*PT1 - PT2*(PT3-RAD2))) / PT2
      WOUT = (-PT1 + SQRT(PT1*PT1 - PT2*(PT3-RPT2))) / PT2
      PCAP = WOUT - WINN
      RETURN
      END
C.....
C  ..           SUBROUTINE PATH3
C.....
C  .. CALCULATE PATH LENGTH THROUGH SOLVENT
C.....
      SUBROUTINE PATH3 (WINN, PSOL)
      DIMENSION NTRP(50)
      COMMON /INPT12/ NSOL, EMUS, NTRP
      DIMENSION FS(50,5)
      COMMON /PLANES/ NF, FS
C-
      PSOL = 0.0
      DO 1000 N=1,NF
C-    DOES PLANE N TRAP SOLVENT?
      IF(NTRP(N).EQ.0) GO TO 1000
C-    IS RAY ACUTE TO FACE NORMAL?
      F = FS(N,5)
      IF(F.LE.0.) GO TO 1000
C-    DOES RAY INTERSECT FACE WITHIN CAPILLARY?
      IF(F.GE.WINN) GO TO 1000
C-    RAY GOES THROUGH SOLVENT
      PF = WINN - F
C-    MAXIMUM DISTANCE IS CORRECT
      IF( PF .GE. PSOL ) PSOL = PF
 1000 CONTINUE
      RETURN
      END
C.....
C  ..           SUBROUTINE RTOS1
C.....
C  .. CONVERT VECTOR COMPONENTS IN THE [REAL] FRAME
C  .. TO COMPONENTS IN THE [STAR] FRAME.
C.....
      SUBROUTINE RTOS1(R,S)
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
      DIMENSION R(3), S(3)
      S(1) = GR(1,1)*R(1) + GR(1,2)*R(2) + GR(1,3)*R(3)
      S(2) = GR(2,1)*R(1) + GR(2,2)*R(2) + GR(2,3)*R(3)
      S(3) = GR(3,1)*R(1) + GR(3,2)*R(2) + GR(3,3)*R(3)
      RETURN
      END
C.....
C  ..           SUBROUTINE SDS1
C.....
C  .. GIVEN TWO VECTORS -S1- AND -S2- GET THEIR DOT PRODUCT
C.....
      SUBROUTINE SDS1 (S1, S2, DOT)
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
      DIMENSION S1(3), S2(3)
C-
      X1 = GS(1,1)*S1(1) + GS(1,2)*S1(2) + GS(1,3)*S1(3)
      X2 = GS(2,1)*S1(1) + GS(2,2)*S1(2) + GS(2,3)*S1(3)
      X3 = GS(3,1)*S1(1) + GS(3,2)*S1(2) + GS(3,3)*S1(3)
C-
      DOT = S2(1)*X1 + S2(2)*X2 + S2(3)*X3
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE SOLN4
C.....
C  .. OUTPUT THE SOLVENT INFORMATION
C.....
      SUBROUTINE SOLN4
      DIMENSION NTRP(50)
      COMMON /INPT12/ NSOL, EMUS, NTRP
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      DIMENSION FS(50,5)
      COMMON /PLANES/ NF, FS
C-
      IF (NSOL .NE. 0) GO TO 1001
      WRITE(NOU,9001)
 9001 FORMAT(/1X,'NO TRAPPED SOLVENT')
      RETURN
 1001 WRITE(NOU,9002)
 9002 FORMAT(/1X,'SOLVENT TRAPPED BY ONE OR MORE FACES')
      WRITE(NOU,9003) EMUS
 9003 FORMAT(1X,'MU OF SOLVENT(1/MM): ',F7.3)
      WRITE(NOU,9004)
 9004 FORMAT(1X,'STATUS: NOT TRAPPING(0); TRAPPING(1)')
      WRITE(NOU,9005)
 9005 FORMAT(1X,'FACE#',3X,'HKL',12X,'STATUS')
      DO 1002 J=1,NF
      WRITE(NOU,9006) J, (FS(J,I),I=1,3), NTRP(J)
 9006 FORMAT(1X,I5,3F5.0,5X,I1)
 1002 CONTINUE
      RETURN
      END
C.....
C  ..           SUBROUTINE TITL2
C.....
C  .. OUTPUT THE TITLE DATE AND PROGRAM HEADER
C.....
      SUBROUTINE TITL2
      CHARACTER TITLE*80, ATIME*8, ADATE*9
      COMMON /INPT01/ TITLE, ATIME, ADATE
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
      CALL TIME(ATIME)
      CALL DATE(ADATE)
      WRITE(NOU,9001)  ATIME, ADATE
 9001 FORMAT(/1X,'PROGRAM ABSORB1...VERSION 1',5X,A,3X,A)
      WRITE(NOU,9002) TITLE
 9002 FORMAT(1X,A)
      RETURN
      END
C.....
C  ..           SUBROUTINE TWOCO
C.....
C  .. "OH.  I GET IT.    NOBODY CARES!"
C  ..          THE FOONEWAL DEWECTOR, 1982
C  .. FIND OUT IF TWO SETS OF THREE NUMBERS
C  .. HAVE ANY TWO IN COMMON.  A ROUTINE TO
C  .. DISTINGUISH TRUE EDGES FROM GENERAL
C  .. DIAGONALS.
C.....
C  .. NTWO = 0  NO MATCH
C  .. NTWO = 1  A MATCH
C.....
      SUBROUTINE TWOCO (I1,I2,I3, J1,J2,J3, NTWO)
      NTWO = 0
C-
      IF( I1.EQ.J1 .AND. I2.EQ.J2 ) NTWO = 1
      IF( I1.EQ.J1 .AND. I2.EQ.J3 ) NTWO = 1
      IF( I1.EQ.J1 .AND. I3.EQ.J2 ) NTWO = 1
      IF( I1.EQ.J1 .AND. I3.EQ.J3 ) NTWO = 1
C-
      IF( I1.EQ.J2 .AND. I2.EQ.J1 ) NTWO = 1
      IF( I1.EQ.J2 .AND. I2.EQ.J3 ) NTWO = 1
      IF( I1.EQ.J2 .AND. I3.EQ.J1 ) NTWO = 1
      IF( I1.EQ.J2 .AND. I3.EQ.J3 ) NTWO = 1
C-
      IF( I1.EQ.J3 .AND. I2.EQ.J1 ) NTWO = 1
      IF( I1.EQ.J3 .AND. I2.EQ.J2 ) NTWO = 1
      IF( I1.EQ.J3 .AND. I3.EQ.J1 ) NTWO = 1
      IF( I1.EQ.J3 .AND. I3.EQ.J2 ) NTWO = 1
C-
      IF( I2.EQ.J1 .AND. I1.EQ.J2 ) NTWO = 1
      IF( I2.EQ.J1 .AND. I1.EQ.J3 ) NTWO = 1
      IF( I2.EQ.J1 .AND. I3.EQ.J2 ) NTWO = 1
      IF( I2.EQ.J1 .AND. I3.EQ.J3 ) NTWO = 1
C-
      IF( I2.EQ.J2 .AND. I1.EQ.J1 ) NTWO = 1
      IF( I2.EQ.J2 .AND. I1.EQ.J3 ) NTWO = 1
      IF( I2.EQ.J2 .AND. I3.EQ.J1 ) NTWO = 1
      IF( I2.EQ.J2 .AND. I3.EQ.J3 ) NTWO = 1
C-
      IF( I2.EQ.J3 .AND. I1.EQ.J1 ) NTWO = 1
      IF( I2.EQ.J3 .AND. I1.EQ.J2 ) NTWO = 1 
      IF( I2.EQ.J3 .AND. I3.EQ.J1 ) NTWO = 1
      IF( I2.EQ.J3 .AND. I3.EQ.J2 ) NTWO = 1
C-
      IF( I3.EQ.J1 .AND. I1.EQ.J2 ) NTWO = 1
      IF( I3.EQ.J1 .AND. I1.EQ.J3 ) NTWO = 1
      IF( I3.EQ.J1 .AND. I2.EQ.J2 ) NTWO = 1
      IF( I3.EQ.J1 .AND. I2.EQ.J3 ) NTWO = 1
C-
      IF( I3.EQ.J2 .AND. I1.EQ.J1 ) NTWO = 1
      IF( I3.EQ.J2 .AND. I1.EQ.J3 ) NTWO = 1
      IF( I3.EQ.J2 .AND. I2.EQ.J1 ) NTWO = 1
      IF( I3.EQ.J2 .AND. I2.EQ.J3 ) NTWO = 1
C-
      IF( I3.EQ.J3 .AND. I1.EQ.J1 ) NTWO = 1
      IF( I3.EQ.J3 .AND. I1.EQ.J2 ) NTWO = 1
      IF( I3.EQ.J3 .AND. I2.EQ.J1 ) NTWO = 1
      IF( I3.EQ.J3 .AND. I2.EQ.J2 ) NTWO = 1
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE UNIVR
C.....
C  .. GIVEN THE COMPONENTS OF A REAL VECTOR DVEC
C  .. GET THE COMPONENTS OF THE UNIT VECTOR
C  .. UVEC PARALLEL TO DVEC
C.....
      SUBROUTINE UNIVR (DVEC, UVEC)
      DIMENSION DVEC(3), UVEC(3)
C-
      UVEC(1) = 0.
      UVEC(2) = 0.
      UVEC(3) = 0.
      D1 = DVEC(1)
      D2 = DVEC(2)
      D3 = DVEC(3)
      CALL LREAL (D1,D2,D3, DREL)
      IF (DREL.LE.0.0) RETURN
      UVEC(1) = D1/DREL
      UVEC(2) = D2/DREL
      UVEC(3) = D3/DREL
      RETURN
      END
C.....
C  ..           SUBROUTINE UNIVS
C.....
C  .. GIVEN THE COMPONENTS OF A RECIPROCAL VECTOR
C  .. DVEC, GET THE COMPONENTS OF THE UNIT VECTOR
C  .. UVEC PARALLEL TO DVEC
C.....
      SUBROUTINE UNIVS (DVEC, UVEC)
      DIMENSION DVEC(3), UVEC(3)
C-
      UVEC(1) = 0.
      UVEC(2) = 0.
      UVEC(3) = 0.
      D1 = DVEC(1)
      D2 = DVEC(2)
      D3 = DVEC(3)
      CALL LSTAR (D1,D2,D3, DSTR)
      IF (DSTR.LE.0.0) RETURN
      UVEC(1) = D1/DSTR
      UVEC(2) = D2/DSTR
      UVEC(3) = D3/DSTR
      RETURN
      END
C.....
C  ..           SUBROUTINE VERT1
C.....
C  .. FROM A LIST OF UNIT VECTORS PERPENDICULAR TO THE
C  .. CRYSTAL FACES AND THEIR PERPENDICULAR DISTANCES
C  .. GENERATE THE  [REAL] COMPONENTS OF THE VERTEX
C  .. VECTORS.
C.....
      SUBROUTINE VERT1
      DIMENSION FS(50,5)
      COMMON /PLANES/ NF, FS
      DIMENSION VR(150,3), MM(150,3)
      COMMON /VERTEX/ NV, VR, MM
      DIMENSION FF(3,3),FI(3,3), DT(3),POR(3)
C-
      NV = 0
      IF (NF .LE. 3)  RETURN
      DO 30 I=1,NF-2
      FF(1,1) = FS(I,1)
      FF(1,2) = FS(I,2)
      FF(1,3) = FS(I,3)
      DT(1)   = FS(I,4)
C-
      DO 20 J=I+1,NF-1
      FF(2,1) = FS(J,1)
      FF(2,2) = FS(J,2)
      FF(2,3) = FS(J,3)
      DT(2)   = FS(J,4)
C-
      DO 10 K=J+1,NF
      FF(3,1) = FS(K,1)
      FF(3,2) = FS(K,2)
      FF(3,3) = FS(K,3)
      DT(3)   = FS(K,4)
C-
      CALL AMI3 (FF, FI, DET)
      IF (DET.EQ.0.)  GO TO 10
      CALL MTV1 (FI, DT, POR)
      CALL INNIE (POR, NIN)
      IF (NIN.EQ.-1)  GO TO 10
C-
      NV = NV + 1
      VR(NV,1) = POR(1)
      VR(NV,2) = POR(2)
      VR(NV,3) = POR(3)
C-
      MM(NV,1) = I
      MM(NV,2) = J
      MM(NV,3) = K
C-
   10 CONTINUE
   20 CONTINUE
   30 CONTINUE
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE VERT2
C.....
C  .. GIVEN A PARTICULAR CRYSTAL PLANE -F1- GENERATE
C  .. THE COMPONENTS IN THE [REAL] FRAMES
C  .. OF ALL THE VERTICES FORMED BY -F1-.
C.....
      SUBROUTINE VERT2(F1)
      DIMENSION FS(50,5)
      COMMON /PLANES/ NF, FS
      DIMENSION VR(150,3), MM(150,3)
      COMMON /VERTEX/ NV, VR, MM
      DIMENSION FF(3,3),FI(3,3), DT(3),POR(3)
      DIMENSION F1(4)
C-
      NV = 0
      IF (NF .LE. 2)  RETURN
      FF(1,1) = F1(1)
      FF(1,2) = F1(2)
      FF(1,3) = F1(3)
      DT(1)   = F1(4)
C-
      DO 20 J=1,NF-1
      FF(2,1) = FS(J,1)
      FF(2,2) = FS(J,2)
      FF(2,3) = FS(J,3)
      DT(2)   = FS(J,4)
C-
      DO 10 K=J+1,NF
      FF(3,1) = FS(K,1)
      FF(3,2) = FS(K,2)
      FF(3,3) = FS(K,3)
      DT(3)   = FS(K,4)
C-
      CALL AMI3 (FF, FI, DET)
      IF (DET.EQ.0.)  GO TO 10
      CALL MTV1 (FI, DT, POR)
      CALL INNIE (POR, NIN)
      IF (NIN.EQ.-1)  GO TO 10
C-
      NV = NV + 1
      VR(NV,1) = POR(1)
      VR(NV,2) = POR(2)
      VR(NV,3) = POR(3)
C-
   10 CONTINUE
   20 CONTINUE
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE VERT3
C.....
C  .. GIVEN TWO PARTICULAR CRYSTAL PLANES -F1- AND -F2-
C  .. GENERATE THE [REAL] COMPONENTS OF ALL
C  .. THE VERTICES FORMED WITH -F1- AND -F2-.
C.....
      SUBROUTINE VERT3(F1,F2)
      DIMENSION FS(50,5)
      COMMON /PLANES/ NF, FS
      DIMENSION VR(150,3), MM(150,3)
      COMMON /VERTEX/ NV, VR, MM
      DIMENSION FF(3,3),FI(3,3), DT(3),POR(3)
      DIMENSION F1(4), F2(4)
C-
      NV = 0
      IF (NF .LE. 1)  RETURN
      FF(1,1) = F1(1)
      FF(1,2) = F1(2)
      FF(1,3) = F1(3)
      DT(1)   = F1(4)
C-
      FF(2,1) = F2(1)
      FF(2,2) = F2(2)
      FF(2,3) = F2(3)
      DT(2)   = F2(4)
C-
      DO 10 K=1,NF
      FF(3,1) = FS(K,1)
      FF(3,2) = FS(K,2)
      FF(3,3) = FS(K,3)
      DT(3)   = FS(K,4)
C-
      CALL AMI3 (FF, FI, DET)
      IF (DET.EQ.0.)     GO TO 10
      CALL MTV1 (FI, DT, POR)
      CALL INNIE (POR, NIN)
      IF (NIN.EQ.-1)   GO TO 10
C-
      NV = NV + 1
      VR(NV,1) = POR(1)
      VR(NV,2) = POR(2)
      VR(NV,3) = POR(3)
C-
   10 CONTINUE
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE VERT4
C.....
C  .. GIVEN THREE PARTICULAR FACE VECTORS -F1- -F2- -F3-
C  .. WITH COMPONENTS IN THE [STAR] FRAME GET THEIR POINT
C  .. OF INTERSECTION (VERTEX) IN THE [REAL] FRAMES.
C  .. DO NOT CHECK IF POINT IS IN CRYSTAL SINCE
C  .. IT MUST BE BY THE METHOD OF GENERATION.
C.....
      SUBROUTINE VERT4(F1,F2,F3, POR)
      DIMENSION FF(3,3),FI(3,3), DT(3),POR(3)
      DIMENSION F1(4), F2(4), F3(4)
C-
      FF(1,1) = F1(1)
      FF(1,2) = F1(2)
      FF(1,3) = F1(3)
      DT(1)   = F1(4)
C-
      FF(2,1) = F2(1)
      FF(2,2) = F2(2)
      FF(2,3) = F2(3)
      DT(2)   = F2(4)
C-
      FF(3,1) = F3(1)
      FF(3,2) = F3(2)
      FF(3,3) = F3(3)
      DT(3)   = F3(4)
C-
      CALL AMI3 (FF, FI, DET)
      IF (DET.EQ.0.)     RETURN
      CALL MTV1 (FI, DT, POR)
C-
      RETURN
      END
C.....
C  ..           SUBROUTINE VERT5
C.....
C  .. GIVEN A SET OF VERTEX VECTORS AND THEIR CONNECTION
C  .. TABLE, PRINT OR TYPE OUT THE LIST OF VECTORS.
C.....
      SUBROUTINE VERT5
      DIMENSION VR(150,3), MM(150,3)
      COMMON /VERTEX/ NV, VR, MM
      COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU
C-
      WRITE(NOU,9001)
 9001 FORMAT(/1X,'VERTEX VECTORS')
      WRITE(NOU,9002)
 9002 FORMAT(1X,'  #','    FACES   ',3X,
     * '     [REAL] COMPONENTS OF VERTICES')
      DO 10 I=1,NV
      WRITE(NOU,9003) I,(MM(I,J),J=1,3),(VR(I,J),J=1,3)
 9003 FORMAT(1X,I3,3I4,3X,3E13.4)
   10 CONTINUE
      RETURN
      END
C.....
C  ..           SUBROUTINE VSTAR
C.....
C  .. GET THE SCALE FACTOR WHICH YIELDS THE VOLUME
C  .. OF THE GRID POINT WHEN MULTIPLIED BY THE THREE
C  .. GAUSSIAN FRACTIONAL LENGTHS EPS1, EPS2, EPS3
C  .. SCALE = D*1 X D*2 X D*3 / DET X V*
C.....
      SUBROUTINE VSTAR (VSTR)
      DIMENSION NGP(3), D1(3),D2(3),D3(3)
      COMMON /INPT08/ NGP, D1,D2,D3
      DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3)
      COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM
      DIMENSION FHKL(3,3), GHKL(3,3)
C-
      FH = D1(1)
      FK = D1(2)
      FL = D1(3)
      CALL LSTAR (FH,FK,FL, D1STR)
      FH = D2(1)
      FK = D2(2)
      FL = D2(3)
      CALL LSTAR (FH,FK,FL, D2STR)
      FH = D3(1)
      FK = D3(2)
      FL = D3(3)
      CALL LSTAR (FH,FK,FL, D3STR)
      FHKL(1,1) = D1(1)
      FHKL(1,2) = D1(2)
      FHKL(1,3) = D1(3)
      FHKL(2,1) = D2(1)
      FHKL(2,2) = D2(2)
      FHKL(2,3) = D2(3)
      FHKL(3,1) = D3(1)
      FHKL(3,2) = D3(2)
      FHKL(3,3) = D3(3)
      CALL AMI3 (FHKL, GHKL,VDET)
      VSTR = (D1STR*D2STR*D3STR) / (VDET*CVS)
      RETURN
      END
