C        VERSION FOR GX PACKAGE APRIL 1985                              0001.000
C     **************                                                    0002.000
      PROGRAM NDIFKWM                                                   0003.000
C     **************                                                    0004.000
C                                                                       0005.000
C   EMPIRICAL ABSORPTION CORRECTION PROGRAM                             0006.000
C                                                                       0007.000
C   A low order term Fourier series is used to model an absorption surfa0008.000
C   for differences between the observed structure factors and those    0009.000
C   obtained from a structure factor calculation (after isotropic least 0010.000
C   squares refinement).   Spherical polar angles are used to define the0011.000
C   incident and diffracted beam path directions.   This reduces the num0012.000
C   of angles required to four - viz. PHI(p), MU(p), PHI(s) and MU(s).  0013.000
C   A theta-dependent correction is evaluated to allow for diffracted   0014.000
C   beams with different path lengths occurring at the same polar angles0015.000
C   Again, a low order term Fourier series is used; in Bragg angle, THET0016.000
C   This version is general for any 4-circle diffractometer data collect0017.000
C   geometry: however, routine ANGLES may require extra code for the    0018.000
C   calculation of the Eulerian angles from the Miller indices using the0019.000
C   orientation matrix.                                                 0020.000
C                                                                       0021.000
C   Author - Nigel Walker, Chemistry Dept., Queen Mary College, London E0022.000
C   Reference - WALKER & STUART (1983) Acta Cryst. A39, 158-166.        0023.000
C               WALKER & STUART (1984) Acta Cryst. (in preparation).    0024.000
C   GLASGOW GX DIFABS MARK II 18-JUNE-84 KW MUIR                        0025.000
C                                                                       0026.000
C   DIMENSIONING OF ARRAYS IN COMMON /NORMAL/ AND  SUBROUTINE SOLVEN    0027.000
C   ----------------------------------------------------------------    0028.000
C   If NCO=NACOEF+NBCOEF, then arrays ANORM & VEC should have minimum   0029.000
C   dimensions (NCO,NCO); ATB & TERMS (NCO); ROWJ (2*NCO); and WORK (5*N0030.000
C                                                                       0031.000
C   VERSION 8.08                                                        0032.000
C   THIS VERSION USES NON-LINEAR LEAST-SQUARES WITH UNIT WEIGHTS        0033.000
C   INPUTS SHELX 'LIST 5' DATA AND OUTPUTS 'HKLF 3' TYPE DATA           0034.000
C   GENERAL ROUTINE FOR CALCULATION OF SPHERICAL POLAR ANGLES           0035.000
C                                                                       0036.000
      LOGICAL TCOR                                                      0037.000
      CHARACTER*9 TODAY                                                 0038.000
      CHARACTER*8 NOW                                                   0039.000
      REAL DUM(4),X(3),MU,MUMIN,MUMAX                                   0040.000
      INTEGER NP(16)                                                    0041.000
CGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG0042.000
      CHARACTER COMPID(72),LINE(148)                                    0043.000
      CHARACTER*8 LABEL(200),FTYPE(13)                                  0044.000
      INTEGER CN(200),CONTNT(13),FCH(20),FWD(20)                        0045.000
      COMMON /MODELC/COMPID,LABEL,FTYPE                                 0046.000
      COMMON /MODEL/ NATOM,NR,NT,NTYPE,NFTYPE(200),P(3,200),            0047.000
     1CN,SOF(200),UIJ(6,200),EP(3,200),ESOF(200),EUIJ(6,200),           0048.000
     2WAVE,CELL(6),ECELL(6),R(24,3,4),T(3,4),ICENT,                     0049.000
     3SFAC(14,13),CONTNT,SCALE(4),ETA,NQ,EXTNCT                         0050.000
CGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG0051.000
      COMMON /NORMAL/ ANORM(44,44),ATB(44),ROWJ(88),TERMS(44)           0052.000
      COMMON /SOLVE/ NSIZE,FILTER,FILSTA,FILDEC,FILEND                  0053.000
      COMMON /REFL/ FC,FO,FM,PHI(2),MU(2),SINTH,COSTH                   0054.000
      COMMON /DRVR/ LI,LD,LO,LF,LP,NCY,NSHAPE,NACOEF,NBCOEF,SFMU,TCOR   0055.000
      COMMON /SRES/ NUMREF,SUMRES,SUMFS,SUMRSQ,SUMFSQ                   0056.000
      COMMON/ORIENT/TWOTHZ,OMZ,PHIZ,PLUS(3),IPOS,IDIFRT,ORIENT(3,3),DST 0057.000
      COMMON /CELLP/ TA,TB,TC,ZA,ZB,ZC                                  0058.000
      DATA STHMAX,MUMIN,MUMAX /0.,100.,-100./                           0059.000
C                                                                       0060.000
C---------------------------------------------------------------------  0061.000
C   PDP11/RSX & VAX/VMS SPECIFIC CALLS                                  0062.000
C     (simply 'comment-out' on other machines)                          0063.000
C                                                                       0064.000
C     CALL ERRSET(63,.TRUE.,.FALSE.,.FALSE.,.FALSE.,31)                 0065.000
C     CALL DATE(TODAY)                                                  0066.000
C     CALL TIME(NOW)                                                    0067.000
C---------------------------------------------------------------------  0068.000
C                                                                       0069.000
C   SET UP I/O UNIT NUMBERS                                             0070.000
      LI=1                                                              0071.000
      LD=2                                                              0072.000
      LO=3                                                              0073.000
      LF=4                                                              0074.000
      LP=8                                                              0075.000
C                                                                       0076.000
C   OUTPUT HEADER.                                                      0077.000
      WRITE(LP,1200)                                                    0078.000
1200  FORMAT(' DIFABS (V08.8) VERSION USING NON-LINEAR LEAST-SQUARES'/  0079.000
     1' FOR ENRAF-NONIUS CAD4 DATA COLLECTED IN BISECTING GEOMETRY')    0080.000
C                                                                       0081.000
C     *********                                                         0082.000
      NACOEF=34                                                         0083.000
      NBCOEF=10                                                         0084.000
C     *********                                                         0085.000
C                                                                       0086.000
C                                                                       0087.000
C     **************                                                    0088.000
      NSIZE=44                                                          0089.000
      FILTER=0.00005                                                    0090.000
      FILSTA=0.0001                                                     0091.000
      FILDEC=0.00005                                                    0092.000
      FILEND=0.0                                                        0093.000
      NCYCLE=5                                                          0094.000
C     **************                                                    0095.000
C                                                                       0096.000
C                                                                       0097.000
C   READ INPUT PARAMETERS                                               0098.000
CGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG0099.000
      LP=6                                                              0100.000
      CALL START(ORIENT,IS,LP)                                          0101.000
      SFACT=1.0                                                         0102.000
      LP=8                                                              0103.000
      IF(IS.EQ.0)WRITE(LP,1201)                                         0104.000
 1201 FORMAT(' THE CORRECTION WILL NOT DEPEND ON THETA '//)             0105.000
CGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG0106.000
      TCOR=.FALSE.                                                      0107.000
      IF(IS.EQ.1)TCOR=.TRUE.                                            0108.000
C       READ(LI,*)NCYCLE,FILTER,FILSTA,FILDEC,FILEND                    0109.000
      IF(.NOT.TCOR)NCYCLE=NCYCLE-1                                      0110.000
C                                                                       0111.000
C   GET CELL PARAMETERS                                                 0112.000
      CALL ORTOCELL                                                     0113.000
C                                                                       0114.000
C   SET UP DIFFRACTOMETER CHARACTERISTICS                               0115.000
C   the following is for Enraf-Nonius CAD4                              0116.000
      IDIFRT=0                                                          0117.000
C     IDIFRT=1   Syntex P-1                                             0118.000
C     IDIFRT=2   Syntex P21                                             0119.000
C     IDIFRT=3   Picker FACS-I                                          0120.000
C     IDIFRT=4   Philips PW1100                                         0121.000
C                                                                       0122.000
C   SET UP GEOMETRY CHARACTERISTICS                                     0123.000
C   the following is for bisecting geometry                             0124.000
      IPOS=0                                                            0125.000
C     IPOS=1   data collected in parallel position                      0126.000
C     IPOS=2     "      "      " a general   "                          0127.000
C                                                                       0128.000
C   SET ANGLE ZEROES AND DIRECTION SENSES.                              0129.000
C                                                                       0130.000
      TWOTHZ=0.0                                                        0131.000
      OMZ=0.0                                                           0132.000
      PHIZ=0.0                                                          0133.000
      IF(IDIFRT.EQ.4) GO TO 110                                         0134.000
      PLUS(1)=1.0                                                       0135.000
      PLUS(2)=-1.0                                                      0136.000
      PLUS(3)=-1.0                                                      0137.000
      GO TO 120                                                         0138.000
110   PLUS(1)=1.0                                                       0139.000
      PLUS(2)=1.0                                                       0140.000
      PLUS(3)=1.0                                                       0141.000
120   CONTINUE                                                          0142.000
C                                                                       0143.000
C   INITIALIZE RESIDUES                                                 0144.000
      NF=-1                                                             0145.000
      CALL RES(NF)                                                      0146.000
C                                                                       0147.000
C   SET COUNTERS                                                        0148.000
      NREF=0                                                            0149.000
      NTOT=0                                                            0150.000
      NCUTO=0                                                           0151.000
      NCUTC=0                                                           0152.000
      NLOW=0                                                            0153.000
      DO 15 J=1,16                                                      0154.000
15    NP(J)=0                                                           0155.000
C                                                                       0156.000
C   SET UP COEFFICIENT STARTING VALUES (ZERO)                           0157.000
C   AND INITIALIZE NORM EQUATIONS                                       0158.000
      NCO=NACOEF+NBCOEF                                                 0159.000
      DO 20 I=1,NCO                                                     0160.000
      ATB(I)=0.0                                                        0161.000
      TERMS(I)=0.0                                                      0162.000
      DO 25 J=1,I                                                       0163.000
      ANORM(I,J)=0.0                                                    0164.000
25    CONTINUE                                                          0165.000
20    CONTINUE                                                          0166.000
C                                                                       0167.000
C                                                                       0168.000
C     ************                                                      0169.000
      SFMU=WAVE*7.                                                      0170.000
C     ************                                                      0171.000
C                                                                       0172.000
C   NOW GET TO WORK...                                                  0173.000
C                                                                       0174.000
      NCY=1                                                             0175.000
      NSHAPE=NACOEF                                                     0176.000
      WRITE(LP,1556)COMPID                                              0177.000
 1556 FORMAT(1H0,72A1,                                                  0178.000
     1///19X,' *** EMPIRICAL ABSORPTION CORRECTION ***',//,             0179.000
     1'   H   K   L  SINT/L   CHIE   PHIE  OMEGA  PHI 1  PHI 2   MU 1   0180.000
     2MU 2     FOBS    SIGF   NREF',/)                                  0181.000
C                                                                       0182.000
C   DO LOOP THRU REFLECTION DATA.                                       0183.000
C                                                                       0184.000
1000  READ (LD,END=290)JAA,JBB,FF,SIGI,AC,BC,SS                         0185.000
      CALL KRAKJA(JAA,ISET,IH,IK,IL )                                   0186.000
      IUN=0                                                             0187.000
      FOBS=SQRT(FF)                                                     0188.000
      SIGF=0.5*SIGI/FOBS                                                0189.000
      FCALC=SQRT(AC*AC+BC*BC)                                           0190.000
      IF(IH.EQ.0.AND.IK.EQ.0.AND.IL.EQ.0)GOTO 290                       0191.000
      NREF=NREF+1                                                       0192.000
      IUSE=0                                                            0193.000
C   IF UNOBSERVED, WRITE DUMMY RECORD                                   0194.000
      IF(FOBS.LE.0.0) GO TO 250                                         0195.000
      NTOT=NTOT+1                                                       0196.000
C                                                                       0197.000
C                                                                       0198.000
C   CALCULATE SETTING ANGLES FROM ORIENTATION MATRIX.                   0199.000
      STHLSQ=0.25*(IH*IH*TA+IK*IK*TB+IL*IL*TC+IL*IH*ZB                  0200.000
     ++IK*IL*ZA+IH*IK*ZC)                                               0201.000
      SINL=SQRT(STHLSQ)                                                 0202.000
      DST=2.0*SINL                                                      0203.000
      X(1)=IH                                                           0204.000
      X(2)=IK                                                           0205.000
      X(3)=IL                                                           0206.000
      CALL ANGLES (OMEGA,PHIE,CHIE,X,IERR)                              0207.000
      IF(IERR.NE.1) GOTO 147                                            0208.000
      WRITE (LP,1220) IH,IK,IL                                          0209.000
1220  FORMAT(' *** ERROR IN SETTING ANGLES FOR REFLECTION',3I4,' ***')  0210.000
      GOTO 250                                                          0211.000
C                                                                       0212.000
147   SINTH=WAVE*SINL                                                   0213.000
      COSTH=SQRT(1. - SINTH*SINTH)                                      0214.000
C   MAX THETA                                                           0215.000
      IF(SINTH.LE.STHMAX) GOTO 148                                      0216.000
      STHMAX=SINTH                                                      0217.000
      CTHMAX=COSTH                                                      0218.000
148   CONTINUE                                                          0219.000
C                                                                       0220.000
C   CALCULATE SPHERICAL POLAR COORDS FROM 4-CIRCLE DIFFRACT. EULERIAN AN0221.000
C   calculate primary & secondary beam angles                           0222.000
C   for bisecting geometry (not Picker)                                 0223.000
      IF(IPOS.EQ.0.AND.IDIFRT.NE.3) CALL BISECT (CHIE,PHIE)             0224.000
C   or for the general case                                             0225.000
      IF(IPOS.NE.0.OR.IDIFRT.EQ.3) CALL GEN4C (OMEGA,CHIE,PHIE)         0226.000
C                                                                       0227.000
C   MIN AND MAX MU                                                      0228.000
      IF(MU(1).LT.MUMIN) MUMIN=MU(1)                                    0229.000
      IF(MU(1).GT.MUMAX) MUMAX=MU(1)                                    0230.000
      IF(IPOS.EQ.0)GOTO 149                                             0231.000
      IF(MU(2).LT.MUMIN) MUMIN=MU(2)                                    0232.000
      IF(MU(2).GT.MUMAX) MUMAX=MU(2)                                    0233.000
C                                                                       0234.000
C   TAKE THE SQUARE ROOT OF MU * SFMU                                   0235.000
149   IF(MU(1).GT.0.) MU(1)=SQRT(MU(1))*SFMU                            0236.000
      IF(MU(1).LT.0.) MU(1)=-SQRT(ABS(MU(1)))*SFMU                      0237.000
      IF(IPOS.EQ.0)GOTO 150                                             0238.000
      IF(MU(2).GT.0.) MU(2)=SQRT(MU(2))*SFMU                            0239.000
      IF(MU(2).LT.0.) MU(2)=-SQRT(ABS(MU(2)))*SFMU                      0240.000
      GOTO 155                                                          0241.000
C   FOR BISECTING GEOMETRY                                              0242.000
150   MU(2)=MU(1)                                                       0243.000
C                                                                       0244.000
C   COUNT REFLS. IN EACH OCTANT                                         0245.000
155   DO 160 J=1,2                                                      0246.000
      JJ=(J-1)*4                                                        0247.000
      M=0                                                               0248.000
      IF(MU(J).LT.0.)M=8                                                0249.000
      N=PHI(J)/90. + 1.                                                 0250.000
      IOCT=N+JJ+M                                                       0251.000
      NP(IOCT)=NP(IOCT) + 1                                             0252.000
160   CONTINUE                                                          0253.000
C                                                                       0254.000
C                                                                       0255.000
      FO=FOBS*SFACT                                                     0256.000
      FC=FCALC                                                          0257.000
      SIGF=SIGF*SFACT                                                   0258.000
C                                                                       0259.000
C   TEST FOR INCLUSION IN COEFFICIENT DETERMINATION                     0260.000
C                                                                       0261.000
C   DON'T ALLOW REFLS .LT. 3*SIGF                                       0262.000
      FCUT=FO/SIGF                                                      0263.000
      IF(FCUT.GE.3.) GOTO 200                                           0264.000
      NCUTO=NCUTO + 1                                                   0265.000
      GOTO 250                                                          0266.000
C                                                                       0267.000
C   DON'T ALLOW REFLS WITH FCALC .LT. 3*SIGF(OBS)                       0268.000
200   FCUT=FC/SIGF                                                      0269.000
      IF(FCUT.GE.3.) GOTO 202                                           0270.000
      NCUTC=NCUTC + 1                                                   0271.000
      GOTO 250                                                          0272.000
C                                                                       0273.000
C   DON'T ALLOW REFLS IF SCALED FOBS .LT. 0.5*FCALC (POSSIBLE EXTINCTION0274.000
202   TWOBS=FO*2.                                                       0275.000
      IF(FC.LT.TWOBS) GOTO 205                                          0276.000
      NLOW=NLOW + 1                                                     0277.000
      GOTO 250                                                          0278.000
C                                                                       0279.000
C   STORE PHI,MU AND THETA FOR LATER USE                                0280.000
205   IUSE=1                                                            0281.000
250   WRITE(LF)IUSE,FC,FO,PHI,MU,SINTH,COSTH                            0282.000
C                                                                       0283.000
C                                                                       0284.000
C   SET UP OBSERVATIONAL AND NORMAL EQUATIONS                           0285.000
      IF(IUSE.EQ.1) CALL ABNORM                                         0286.000
      IF(NREF.LE.20.OR.MOD(NREF,100).EQ.0)                              0287.000
     1WRITE(LP,1500)IH,IK,IL,SINL,CHIE,PHIE,OMEGA,PHI,MU,FOBS,SIGF,NREF 0288.000
1500      FORMAT(3I4,F8.4,7F7.1,F9.2,F8.1,I6)                           0289.000
      GOTO 1000                                                         0290.000
C                                                                       0291.000
290   CONTINUE                                                          0292.000
C                                                                       0293.000
C                                                                       0294.000
C   REPORT REFLECTION SUMMARY                                           0295.000
      NZERO=NREF-NTOT                                                   0296.000
      WRITE(LP,1600)NREF,NZERO,NCUTO,NCUTC,NLOW                         0297.000
1600  FORMAT(/I6,' Reflections input',                                  0298.000
     +       /I6,' Rejected with Fobs=0',                               0299.000
     +       /I6,' Rejected with Fobs/SigFobs .LT. 3',                  0300.000
     +       /I6,' Rejected with Fcalc/SigFobs .LT. 3',                 0301.000
     +       /I6,' Rejected with Fcalc .GT. 2*Fobs')                    0302.000
C                                                                       0303.000
      CALL ARCTAN (THMAX,CTHMAX,STHMAX)                                 0304.000
C                                                                       0305.000
      WRITE(LP,1700)THMAX,MUMIN,MUMAX,NP                                0306.000
1700  FORMAT(/' Maximum THETA =',F6.2,' degrees',                       0307.000
     +       /' Minimum MU =',F6.2,' degrees',                          0308.000
     +       /' Maximum MU =',F6.2,' degrees',                          0309.000
     + //' NUMBER OF REFLECTION PATHS IN EACH OCTANT',                  0310.000
     + /12X,'Primary',20X,'Secondary'/2(5X,'0-90 -180 -270 -360',4X)    0311.000
     + /' +MU',4I5,8X,4I5 / ' -MU',4I5,8X,4I5)                          0312.000
C                                                                       0313.000
C   REPORT RESIDUALS BEFORE CORRECTION                                  0314.000
      WRITE(LP,1800)                                                    0315.000
1800  FORMAT(//' RESIDUALS BEFORE CORRECTION')                          0316.000
      NF=1                                                              0317.000
      CALL RES(NF)                                                      0318.000
C                                                                       0319.000
C   SOLVE EQUATIONS                                                     0320.000
      CALL SOLVEN                                                       0321.000
C                                                                       0322.000
C                                                                       0323.000
C   LOOP ROUND DOING CYCLES                                             0324.000
C                                                                       0325.000
      IF(TCOR)NSHAPE=NACOEF+NBCOEF                                      0326.000
C                                                                       0327.000
      DO 300 NCY=2,NCYCLE                                               0328.000
C                                                                       0329.000
C   INITIALIZE RESIDUALS AND NORMAL EQUATIONS                           0330.000
      NF=-1                                                             0331.000
      CALL RES(NF)                                                      0332.000
      DO 310 I=1,NSHAPE                                                 0333.000
      ATB(I)=0.0                                                        0334.000
      DO 320 J=1,I                                                      0335.000
      ANORM(I,J)=0.0                                                    0336.000
320   CONTINUE                                                          0337.000
310   CONTINUE                                                          0338.000
C                                                                       0339.000
      REWIND LF                                                         0340.000
      DO 400 NRECRD=1,NREF                                              0341.000
      READ(LF)IUSE,FC,FO,PHI,MU,SINTH,COSTH                             0342.000
      IF(IUSE.EQ.1)CALL ABNORM                                          0343.000
400   CONTINUE                                                          0344.000
C                                                                       0345.000
      WRITE(LP,1900)NCY-1                                               0346.000
1900  FORMAT(//' RESIDUALS AFTER CYCLE',I2)                             0347.000
      NF=1                                                              0348.000
      CALL RES(NF)                                                      0349.000
C                                                                       0350.000
      CALL SOLVEN                                                       0351.000
C                                                                       0352.000
300   CONTINUE                                                          0353.000
C                                                                       0354.000
C                                                                       0355.000
C   APPLY CORRECTION TO ALL DATA                                        0356.000
      CALL ABSPLY(NREF)                                                 0357.000
C                                                                       0358.000
      WRITE(LP,210) (TERMS(J),J=1,NSHAPE)                               0359.000
210   FORMAT(/' ***VALUES FOR COEFFICIENTS***'//,                       0360.000
     +  2(6F10.4/)/ , 2(4F10.4/)/ , 2(7F10.4/)/  , 2(5F10.4/)/ )        0361.000
C                                                                       0362.000
C   PRINT EXAMPLE OF ABSORPTION SURFACE                                 0363.000
      CALL ABSPRI(MUMIN,MUMAX,THMAX)                                    0364.000
C                                                                       0365.000
C                                                                       0366.000
C                                                                       0367.000
      WRITE(LP,1100)                                                    0368.000
 1100 FORMAT(//' DIFABS COMPLETED - SET MODEL FILE SCALE TO UNITY'//)   0369.000
      STOP                                                              0370.000
      END                                                               0371.000
C                                                                       0372.000
C                                                                       0373.000
C                                                                       0374.000
CSUBROUTINESSUBROUTINESSUBROUTINESSUBROUTINESSUBROUTINESSUBROUTINESSUBRO0375.000
CSUBROUTINESSUBROUTINESSUBROUTINESSUBROUTINESSUBROUTINESSUBROUTINESSUBRO0376.000
CSUBROUTINESSUBROUTINESSUBROUTINESSUBROUTINESSUBROUTINESSUBROUTINESSUBRO0377.000
CSUBROUTINESSUBROUTINESSUBROUTINESSUBROUTINESSUBROUTINESSUBROUTINESSUBRO0378.000
C                                                                       0379.000
C                                                                       0380.000
C                                                                       0381.000
C                                                                       0382.000
C                                                                       0383.000
      SUBROUTINE ABNORM                                                 0384.000
C     *****************                                                 0385.000
C                                                                       0386.000
C                                                                       0387.000
C  ROUTINE TO DERIVE ABSORPTION COEFFICIENTS.                           0388.000
C                                                                       0389.000
C   To be modelled:-    Fc = k Fo . A(p,s)                              0390.000
C   where:                                                              0391.000
C         A(p,s) = A(p).A(s).A(p-s)                                     0392.000
C   using model:                                                        0393.000
C         A(p,s) = SUM ?Pn,m ?sin(nop+mup) + Qn,m cos(nop+mup)]?        0394.000
C                 *SUM ?Pn,m ?sin(nos+mus) + Qn,m cos(nos+mus)]?        0395.000
C                 *SUM ?Rn sin(nO) + Sn cos(nO)?                        0396.000
C   where:                                                              0397.000
C     o=phi, u=mu and O=theta                                           0398.000
C     p=primary and s=secondary, p-s=2O, and n and m are integers.      0399.000
C     k = scale factor from least squares refinement.                   0400.000
C                                                                       0401.000
C   Coefficients for zero order terms are set to unity and taken out    0402.000
C   of summation, hence:                                                0403.000
C                   Fc = k Fo . ?1 + A'(p)?.?1 + A'(s)?.?1 + A'(O)?     0404.000
C                                                                       0405.000
C   Terms Pn,m and Qn,m used are:-                                      0406.000
C    a) n=1,2,3,4,5,6 and m=0  b) n=0 and m=1,2,3,4                     0407.000
C    c) n=1,2,3,4 and m=1, and n=m=2,3,4    ie. a total of 34           0408.000
C   and terms Rn and Sn used are:-                                      0409.000
C      n=1,2,3,4,5    ie. a total of 10                                 0410.000
C                                                                       0411.000
C   All coefficients have starting values of zero.                      0412.000
C                                                                       0413.000
C                                                                       0414.000
      LOGICAL TCOR                                                      0415.000
      REAL MU                                                           0416.000
      COMMON /NORMAL/ ANORM(44,44),ATB(44),ROWJ(88),TERMS(44)           0417.000
      COMMON /REFL/ FC,FO,FM,PHI(2),MU(2),SINTH,COSTH                   0418.000
      COMMON /DRVR/ LI,LD,LO,LF,LP,NCY,NSHAPE,NACOEF,NBCOEF,SFMU,TCOR   0419.000
C                                                                       0420.000
C                                                                       0421.000
C   GET TRIG TERMS                                                      0422.000
      CALL ABSTOR                                                       0423.000
C                                                                       0424.000
C   APPLY CORRECTION FACTOR FROM PREVIOUS CYCLE TO FOBS                 0425.000
      IF(NCY.EQ.1) GOTO 5                                               0426.000
      CALL AFACT(PS1,PS2,PS3,ABSFAC)                                    0427.000
      FM=FO*ABSFAC                                                      0428.000
      GOTO 6                                                            0429.000
C                                                                       0430.000
5     FM=FO                                                             0431.000
      PS1=1.0                                                           0432.000
      PS2=1.0                                                           0433.000
      PS3=1.0                                                           0434.000
C                                                                       0435.000
6     NF=0                                                              0436.000
      CALL RES(NF)                                                      0437.000
C                                                                       0438.000
      DEL=FC-FM                                                         0439.000
C                                                                       0440.000
C   CALCULATE DERIVATIVES                                               0441.000
C                                                                       0442.000
      DCONST=FO*PS3                                                     0443.000
      DO 10 J=1,NACOEF                                                  0444.000
      JJ=J+NACOEF                                                       0445.000
      ROWJ(J)=DCONST *(ROWJ(J)*PS2 + PS1*ROWJ(JJ))                      0446.000
10    CONTINUE                                                          0447.000
C                                                                       0448.000
      IF(.NOT.TCOR.OR.NCY.EQ.1) GOTO 30                                 0449.000
C                                                                       0450.000
      NSTART=NACOEF+1                                                   0451.000
      DCONST=FO*PS1*PS2                                                 0452.000
      DO 20 J=NSTART,NSHAPE                                             0453.000
      JJ=J+NACOEF                                                       0454.000
      ROWJ(J)=DCONST*ROWJ(JJ)                                           0455.000
20    CONTINUE                                                          0456.000
C                                                                       0457.000
C                                                                       0458.000
C   NOTE... SHAPE DERIVS. ARE IN ROWJ (1 thro' NSHAPE)                  0459.000
C                                                                       0460.000
C   ADD INTO NORMAL EQUATIONS                                           0461.000
C                                                                       0462.000
C                                                                       0463.000
30    DO 40 K=1,NSHAPE                                                  0464.000
      DO 50 L=1,K                                                       0465.000
      ANORM(K,L)=ANORM(K,L) + ROWJ(K)*ROWJ(L)                           0466.000
50    CONTINUE                                                          0467.000
C                                                                       0468.000
C   R.H.S.                                                              0469.000
      ATB(K)=ATB(K) + ROWJ(K)*DEL                                       0470.000
40    CONTINUE                                                          0471.000
C                                                                       0472.000
C                                                                       0473.000
      RETURN                                                            0474.000
      END                                                               0475.000
C                                                                       0476.000
C                                                                       0477.000
C                                                                       0478.000
C                                                                       0479.000
C                                                                       0480.000
      SUBROUTINE BISECT (CHIE,PHIE)                                     0481.000
C     *****************************                                     0482.000
C                                                                       0483.000
C   CALCULATE SPHERICAL POLAR COORDINATES PHI AND MU                    0484.000
C   FROM EULERIAN DIFFRACTOMETER ANGLES PHIE, CHIE AND THETA            0485.000
C   FOR A CRYSTAL IN THE BISECTING POSITION.                            0486.000
C   ALL ANGLES IN DEGREES.                                              0487.000
C                                                                       0488.000
      REAL MU                                                           0489.000
      COMMON /REFL/ FC,FO,FM,PHI(2),MU(2),SINTH,COSTH                   0490.000
C                                                                       0491.000
      DATA DTOR/0.0174533/                                              0492.000
      DATA RAD/57.2957795/                                              0493.000
C                                                                       0494.000
      RCHIE=CHIE*DTOR                                                   0495.000
      SCHI= SIN(RCHIE)                                                  0496.000
      CCHI= COS(RCHIE)                                                  0497.000
C                                                                       0498.000
C...CALCULATE MU                                                        0499.000
      SMU= SINTH * SCHI                                                 0500.000
      CMU= SQRT(1. - SMU*SMU)                                           0501.000
      CALL ARCTAN(XMU,CMU,SMU)                                          0502.000
C                                                                       0503.000
C                                                                       0504.000
C...CALCULATE DELTA PHI                                                 0505.000
      TANTH= SINTH/COSTH                                                0506.000
      DPHI= ATAN(TANTH * CCHI)*RAD                                      0507.000
C                                                                       0508.000
C                                                                       0509.000
C   CALCULATE PHI AND MU FOR PRIMARY AND SECONDARY BEAM PATHS           0510.000
C   ONLY ASSUMPTION IS A(p)=A(-s)                                       0511.000
C   Huber & Kopfmann (1969) Acta Cryst. A25, 143-152                    0512.000
C                                                                       0513.000
C                                                                       0514.000
C                                                                       0515.000
      PHI(1)= PHIE + DPHI                                               0516.000
      PHI(2)= PHIE + 180. - DPHI                                        0517.000
      CALL NORM(PHI(1))                                                 0518.000
      CALL NORM(PHI(2))                                                 0519.000
C                                                                       0520.000
      MU(1)= XMU                                                        0521.000
      MU(2)= XMU                                                        0522.000
C                                                                       0523.000
      RETURN                                                            0524.000
      END                                                               0525.000
C                                                                       0526.000
C                                                                       0527.000
C                                                                       0528.000
C                                                                       0529.000
C                                                                       0530.000
      SUBROUTINE GEN4C (OMEGA,CHIE,PHIE)                                0531.000
C     **********************************                                0532.000
C                                                                       0533.000
C   GENERAL ROUTINE TO CALCULATE SPHERICAL POLAR COORDINATES            0534.000
C   PHI AND MU FROM EULERIAN DIFFRACTOMETER ANGLES -                    0535.000
C   ALL ANGLES IN DEGREES.                                              0536.000
C   NOTE: OMEGA is defined such that the CHI axis makes an angle        0537.000
C         of THETA + OMEGA with the primary beam direction,             0538.000
C         ie. OMEGA=0 in the bisecting position.                        0539.000
C                                                                       0540.000
      REAL MU,MUP,MUS                                                   0541.000
      COMMON /REFL/ FC,FO,FM,PHI(2),MU(2),SINTH,COSTH                   0542.000
C                                                                       0543.000
      DATA DTOR/0.0174533/                                              0544.000
      DATA RAD/57.2957795/                                              0545.000
C                                                                       0546.000
      THETA= ATAN(SINTH/COSTH)                                          0547.000
      ROMEGA= OMEGA*DTOR                                                0548.000
      OMP= THETA + ROMEGA                                               0549.000
      OMS= THETA - ROMEGA                                               0550.000
      RCHIE=CHIE*DTOR                                                   0551.000
      SCHI= SIN(RCHIE)                                                  0552.000
      CCHI= COS(RCHIE)                                                  0553.000
C                                                                       0554.000
C...CALCULATE MU                                                        0555.000
      SMU= SIN(OMP) * SCHI                                              0556.000
      CMU= SQRT(1. - SMU*SMU)                                           0557.000
      CALL ARCTAN(MUP,CMU,SMU)                                          0558.000
      SMU= SIN(OMS) * SCHI                                              0559.000
      CMU= SQRT(1. - SMU*SMU)                                           0560.000
      CALL ARCTAN(MUS,CMU,SMU)                                          0561.000
C                                                                       0562.000
C                                                                       0563.000
C...CALCULATE DELTA PHI                                                 0564.000
      DPHIP= ATAN(CCHI * TAN(OMP))*RAD                                  0565.000
      DPHIS= ATAN(CCHI * TAN(OMS))*RAD                                  0566.000
C                                                                       0567.000
C                                                                       0568.000
C   CALCULATE PHI AND MU FOR PRIMARY AND SECONDARY BEAM PATHS           0569.000
C   ONLY ASSUMPTION IS A(p)=A(-s)                                       0570.000
C   Huber & Kopfmann (1969) Acta Cryst. A25, 143-152                    0571.000
C                                                                       0572.000
C                                                                       0573.000
C                                                                       0574.000
      PHI(1)= PHIE + DPHIP                                              0575.000
      PHI(2)= PHIE + 180. - DPHIS                                       0576.000
      CALL NORM(PHI(1))                                                 0577.000
      CALL NORM(PHI(2))                                                 0578.000
C                                                                       0579.000
      MU(1)= MUP                                                        0580.000
      MU(2)= MUS                                                        0581.000
C                                                                       0582.000
      RETURN                                                            0583.000
      END                                                               0584.000
C                                                                       0585.000
C                                                                       0586.000
C                                                                       0587.000
C                                                                       0588.000
C                                                                       0589.000
C                                                                       0590.000
      SUBROUTINE ABSTOR                                                 0591.000
C     *****************                                                 0592.000
C                                                                       0593.000
C   CALCULATES SIN & COS OF POLAR ANGLES                                0594.000
C                                                                       0595.000
C   WORK OUT TERMS USING EXPANSIONS FOR TRIG FUNCTIONS.                 0596.000
C                                                                       0597.000
C   NOTE: if bisecting geometry, MU expansions only performed once      0598.000
C         since MU(s)=MU(p).                                            0599.000
C                                                                       0600.000
      LOGICAL TCOR                                                      0601.000
      REAL MU                                                           0602.000
      COMMON /NORMAL/ ANORM(44,44),ATB(44),ROWJ(88),TERMS(44)           0603.000
      COMMON /REFL/ FC,FO,FM,PHI(2),MU(2),ST,CT                         0604.000
      COMMON/ORIENT/TWOTHZ,OMZ,PHIZ,PLUS(3),IPOS,IDIFRT,ORIENT(3,3),DST 0605.000
      COMMON /DRVR/ LI,LD,LO,LF,LP,NCY,NSHAPE,NACOEF,NBCOEF,SFMU,TCOR   0606.000
C                                                                       0607.000
      DATA DTOR/0.0174533/                                              0608.000
C                                                                       0609.000
C                                                                       0610.000
      SPP=SIN(PHI(1)*DTOR)                                              0611.000
      CPP=COS(PHI(1)*DTOR)                                              0612.000
C                                                                       0613.000
      SPS=SIN(PHI(2)*DTOR)                                              0614.000
      CPS=COS(PHI(2)*DTOR)                                              0615.000
C                                                                       0616.000
      SMP=SIN(MU(1)*DTOR)                                               0617.000
      CMP=COS(MU(1)*DTOR)                                               0618.000
C                                                                       0619.000
      IF(IPOS.EQ.0.AND.IDIFRT.NE.3)GOTO 5                               0620.000
      SMS=SIN(MU(2)*DTOR)                                               0621.000
      CMS=COS(MU(2)*DTOR)                                               0622.000
C                                                                       0623.000
C                                                                       0624.000
5     SPP2=CPP*SPP                                                      0625.000
      SPP2=SPP2+SPP2                                                    0626.000
      CPP2=CPP*CPP-SPP*SPP                                              0627.000
C                                                                       0628.000
      SPP3=CPP*SPP2+CPP2*SPP                                            0629.000
      CPP3=CPP*CPP2-SPP*SPP2                                            0630.000
C                                                                       0631.000
      SPP4=CPP2*SPP2                                                    0632.000
      SPP4=SPP4+SPP4                                                    0633.000
      CPP4=CPP2*CPP2-SPP2*SPP2                                          0634.000
C                                                                       0635.000
      SPP5=CPP2*SPP3+CPP3*SPP2                                          0636.000
      CPP5=CPP2*CPP3-SPP2*SPP3                                          0637.000
C                                                                       0638.000
      SPP6=CPP3*SPP3                                                    0639.000
      SPP6=SPP6+SPP6                                                    0640.000
      CPP6=CPP3*CPP3-SPP3*SPP3                                          0641.000
C                                                                       0642.000
      SPS2=CPS*SPS                                                      0643.000
      SPS2=SPS2+SPS2                                                    0644.000
      CPS2=CPS*CPS-SPS*SPS                                              0645.000
C                                                                       0646.000
      SPS3=CPS*SPS2+CPS2*SPS                                            0647.000
      CPS3=CPS*CPS2-SPS*SPS2                                            0648.000
C                                                                       0649.000
      SPS4=CPS2*SPS2                                                    0650.000
      SPS4=SPS4+SPS4                                                    0651.000
      CPS4=CPS2*CPS2-SPS2*SPS2                                          0652.000
C                                                                       0653.000
      SPS5=CPS2*SPS3+CPS3*SPS2                                          0654.000
      CPS5=CPS2*CPS3-SPS2*SPS3                                          0655.000
C                                                                       0656.000
      SPS6=CPS3*SPS3                                                    0657.000
      SPS6=SPS6+SPS6                                                    0658.000
      CPS6=CPS3*CPS3-SPS3*SPS3                                          0659.000
C                                                                       0660.000
      SMP2=CMP*SMP                                                      0661.000
      SMP2=SMP2+SMP2                                                    0662.000
      CMP2=CMP*CMP-SMP*SMP                                              0663.000
C                                                                       0664.000
      SMP3=CMP2*SMP+CMP*SMP2                                            0665.000
      CMP3=CMP2*CMP-SMP2*SMP                                            0666.000
C                                                                       0667.000
      SMP4=CMP2*SMP2                                                    0668.000
      SMP4=SMP4+SMP4                                                    0669.000
      CMP4=CMP2*CMP2-SMP2*SMP2                                          0670.000
C                                                                       0671.000
      IF(IPOS.EQ.0.AND.IDIFRT.NE.3)GOTO 10                              0672.000
      SMS2=CMS*SMS                                                      0673.000
      SMS2=SMS2+SMS2                                                    0674.000
      CMS2=CMS*CMS-SMS*SMS                                              0675.000
C                                                                       0676.000
      SMS3=CMS2*SMS+CMS*SMS2                                            0677.000
      CMS3=CMS2*CMS-SMS2*SMS                                            0678.000
C                                                                       0679.000
      SMS4=CMS2*SMS2                                                    0680.000
      SMS4=SMS4+SMS4                                                    0681.000
      CMS4=CMS2*CMS2-SMS2*SMS2                                          0682.000
      GOTO 15                                                           0683.000
C                                                                       0684.000
C   FOR BISECTING GEOMETRY                                              0685.000
10    SMS=SMP                                                           0686.000
      CMS=CMP                                                           0687.000
      SMS2=SMP2                                                         0688.000
      CMS2=CMP2                                                         0689.000
      SMS3=SMP3                                                         0690.000
      CMS3=CMP3                                                         0691.000
      SMS4=SMP4                                                         0692.000
      CMS4=CMP4                                                         0693.000
C                                                                       0694.000
15    SXP11=CPP*SMP+SPP*CMP                                             0695.000
      CXP11=CPP*CMP-SPP*SMP                                             0696.000
      SXS11=CPS*SMS+SPS*CMS                                             0697.000
      CXS11=CPS*CMS-SPS*SMS                                             0698.000
C                                                                       0699.000
      SXP22=CPP2*SMP2+SPP2*CMP2                                         0700.000
      CXP22=CPP2*CMP2-SPP2*SMP2                                         0701.000
      SXS22=CPS2*SMS2+SPS2*CMS2                                         0702.000
      CXS22=CPS2*CMS2-SPS2*SMS2                                         0703.000
C                                                                       0704.000
      SXP33=CPP3*SMP3+SPP3*CMP3                                         0705.000
      CXP33=CPP3*CMP3-SPP3*SMP3                                         0706.000
      SXS33=CPS3*SMS3+SPS3*CMS3                                         0707.000
      CXS33=CPS3*CMS3-SPS3*SMS3                                         0708.000
C                                                                       0709.000
      SXP44=CPP4*SMP4+SPP4*CMP4                                         0710.000
      CXP44=CPP4*CMP4-SPP4*SMP4                                         0711.000
      SXS44=CPS4*SMS4+SPS4*CMS4                                         0712.000
      CXS44=CPS4*CMS4-SPS4*SMS4                                         0713.000
C                                                                       0714.000
      SXP21=CPP2*SMP+SPP2*CMP                                           0715.000
      CXP21=CPP2*CMP-SPP2*SMP                                           0716.000
      SXS21=CPS2*SMS+SPS2*CMS                                           0717.000
      CXS21=CPS2*CMS-SPS2*SMS                                           0718.000
C                                                                       0719.000
      SXP31=CPP3*SMP+SPP3*CMP                                           0720.000
      CXP31=CPP3*CMP-SPP3*SMP                                           0721.000
      SXS31=CPS3*SMS+SPS3*CMS                                           0722.000
      CXS31=CPS3*CMS-SPS3*SMS                                           0723.000
C                                                                       0724.000
      SXP41=CPP4*SMP+SPP4*CMP                                           0725.000
      CXP41=CPP4*CMP-SPP4*SMP                                           0726.000
      SXS41=CPS4*SMS+SPS4*CMS                                           0727.000
      CXS41=CPS4*CMS-SPS4*SMS                                           0728.000
C                                                                       0729.000
C                                                                       0730.000
      ROWJ(1)=SPP                                                       0731.000
      ROWJ(2)=SPP2                                                      0732.000
      ROWJ(3)=SPP3                                                      0733.000
      ROWJ(4)=SPP4                                                      0734.000
      ROWJ(5)=SPP5                                                      0735.000
      ROWJ(6)=SPP6                                                      0736.000
C                                                                       0737.000
      ROWJ(7)=CPP                                                       0738.000
      ROWJ(8)=CPP2                                                      0739.000
      ROWJ(9)=CPP3                                                      0740.000
      ROWJ(10)=CPP4                                                     0741.000
      ROWJ(11)=CPP5                                                     0742.000
      ROWJ(12)=CPP6                                                     0743.000
C                                                                       0744.000
      ROWJ(13)=SMP                                                      0745.000
      ROWJ(14)=SMP2                                                     0746.000
      ROWJ(15)=SMP3                                                     0747.000
      ROWJ(16)=SMP4                                                     0748.000
C                                                                       0749.000
      ROWJ(17)=CMP                                                      0750.000
      ROWJ(18)=CMP2                                                     0751.000
      ROWJ(19)=CMP3                                                     0752.000
      ROWJ(20)=CMP4                                                     0753.000
C                                                                       0754.000
      ROWJ(21)=SXP11                                                    0755.000
      ROWJ(22)=SXP21                                                    0756.000
      ROWJ(23)=SXP31                                                    0757.000
      ROWJ(24)=SXP41                                                    0758.000
      ROWJ(25)=SXP22                                                    0759.000
      ROWJ(26)=SXP33                                                    0760.000
      ROWJ(27)=SXP44                                                    0761.000
C                                                                       0762.000
      ROWJ(28)=CXP11                                                    0763.000
      ROWJ(29)=CXP21                                                    0764.000
      ROWJ(30)=CXP31                                                    0765.000
      ROWJ(31)=CXP41                                                    0766.000
      ROWJ(32)=CXP22                                                    0767.000
      ROWJ(33)=CXP33                                                    0768.000
      ROWJ(34)=CXP44                                                    0769.000
C                                                                       0770.000
      ROWJ(35)=SPS                                                      0771.000
      ROWJ(36)=SPS2                                                     0772.000
      ROWJ(37)=SPS3                                                     0773.000
      ROWJ(38)=SPS4                                                     0774.000
      ROWJ(39)=SPS5                                                     0775.000
      ROWJ(40)=SPS6                                                     0776.000
C                                                                       0777.000
      ROWJ(41)=CPS                                                      0778.000
      ROWJ(42)=CPS2                                                     0779.000
      ROWJ(43)=CPS3                                                     0780.000
      ROWJ(44)=CPS4                                                     0781.000
      ROWJ(45)=CPS5                                                     0782.000
      ROWJ(46)=CPS6                                                     0783.000
C                                                                       0784.000
      ROWJ(47)=SMS                                                      0785.000
      ROWJ(48)=SMS2                                                     0786.000
      ROWJ(49)=SMS3                                                     0787.000
      ROWJ(50)=SMS4                                                     0788.000
C                                                                       0789.000
      ROWJ(51)=CMS                                                      0790.000
      ROWJ(52)=CMS2                                                     0791.000
      ROWJ(53)=CMS3                                                     0792.000
      ROWJ(54)=CMS4                                                     0793.000
C                                                                       0794.000
      ROWJ(55)=SXS11                                                    0795.000
      ROWJ(56)=SXS21                                                    0796.000
      ROWJ(57)=SXS31                                                    0797.000
      ROWJ(58)=SXS41                                                    0798.000
      ROWJ(59)=SXS22                                                    0799.000
      ROWJ(60)=SXS33                                                    0800.000
      ROWJ(61)=SXS44                                                    0801.000
C                                                                       0802.000
      ROWJ(62)=CXS11                                                    0803.000
      ROWJ(63)=CXS21                                                    0804.000
      ROWJ(64)=CXS31                                                    0805.000
      ROWJ(65)=CXS41                                                    0806.000
      ROWJ(66)=CXS22                                                    0807.000
      ROWJ(67)=CXS33                                                    0808.000
      ROWJ(68)=CXS44                                                    0809.000
C                                                                       0810.000
      IF(.NOT.TCOR.OR.NCY.EQ.1) GOTO 20                                 0811.000
C                                                                       0812.000
      ST2=CT*ST                                                         0813.000
      ST2=ST2+ST2                                                       0814.000
      CT2=CT*CT-ST*ST                                                   0815.000
C                                                                       0816.000
      ST3=CT*ST2+CT2*ST                                                 0817.000
      CT3=CT*CT2-ST*ST2                                                 0818.000
C                                                                       0819.000
      ST4=CT2*ST2                                                       0820.000
      ST4=ST4+ST4                                                       0821.000
      CT4=CT2*CT2-ST2*ST2                                               0822.000
C                                                                       0823.000
      ST5=CT2*ST3+CT3*ST2                                               0824.000
      CT5=CT2*CT3-ST2*ST3                                               0825.000
C                                                                       0826.000
C                                                                       0827.000
      ROWJ(69)=ST                                                       0828.000
      ROWJ(70)=ST2                                                      0829.000
      ROWJ(71)=ST3                                                      0830.000
      ROWJ(72)=ST4                                                      0831.000
      ROWJ(73)=ST5                                                      0832.000
C                                                                       0833.000
      ROWJ(74)=CT                                                       0834.000
      ROWJ(75)=CT2                                                      0835.000
      ROWJ(76)=CT3                                                      0836.000
      ROWJ(77)=CT4                                                      0837.000
      ROWJ(78)=CT5                                                      0838.000
C                                                                       0839.000
20    RETURN                                                            0840.000
      END                                                               0841.000
C                                                                       0842.000
C                                                                       0843.000
C                                                                       0844.000
C                                                                       0845.000
C                                                                       0846.000
      SUBROUTINE AFACT (PS1,PS2,PS3,ABSFAC)                             0847.000
C     *************************************                             0848.000
C                                                                       0849.000
C   ROUTINE TO CALCULATE ABSORPTION FACTOR                              0850.000
C                                                                       0851.000
      LOGICAL TCOR                                                      0852.000
      COMMON /NORMAL/ ANORM(44,44),ATB(44),ROWJ(88),TERMS(44)           0853.000
      COMMON /DRVR/ LI,LD,LO,LF,LP,NCY,NSHAPE,NACOEF,NBCOEF,SFMU,TCOR   0854.000
C                                                                       0855.000
C                                                                       0856.000
      PS1=1.0                                                           0857.000
      PS2=1.0                                                           0858.000
      PS3=1.0                                                           0859.000
C                                                                       0860.000
      DO 10 J=1,NACOEF                                                  0861.000
      JJ=J+NACOEF                                                       0862.000
      PS1=PS1 + TERMS(J)*ROWJ(J)                                        0863.000
      PS2=PS2 + TERMS(J)*ROWJ(JJ)                                       0864.000
10    CONTINUE                                                          0865.000
C                                                                       0866.000
      IF(.NOT.TCOR.OR.NCY.LE.2) GOTO 30                                 0867.000
C                                                                       0868.000
      NSTART=NACOEF+1                                                   0869.000
      DO 20 J=NSTART,NSHAPE                                             0870.000
      JJ=J+NACOEF                                                       0871.000
      PS3=PS3 + TERMS(J)*ROWJ(JJ)                                       0872.000
20    CONTINUE                                                          0873.000
C                                                                       0874.000
30    ABSFAC = PS1*PS2*PS3                                              0875.000
C                                                                       0876.000
C                                                                       0877.000
      RETURN                                                            0878.000
      END                                                               0879.000
C                                                                       0880.000
C                                                                       0881.000
C                                                                       0882.000
C                                                                       0883.000
C                                                                       0884.000
C                                                                       0885.000
C                                                                       0886.000
C                                                                       0887.000
C                                                                       0888.000
      SUBROUTINE RES (NF)                                               0889.000
C     *******************                                               0890.000
C                                                                       0891.000
C  ROUTINE TO ACCUMULATE RESIDUAL INFORMATION                           0892.000
C                                                                       0893.000
C                                                                       0894.000
      LOGICAL TCOR                                                      0895.000
      REAL MU                                                           0896.000
      COMMON /REFL/ FC,FO,FM,PHI(2),MU(2),SINTH,COSTH                   0897.000
      COMMON /DRVR/ LI,LD,LO,LF,LP,NCY,NSHAPE,NACOEF,NBCOEF,SFMU,TCOR   0898.000
      COMMON /SRES/ NUMREF,SUMRES,SUMFS,SUMRSQ,SUMFSQ                   0899.000
C                                                                       0900.000
      IF(NF.NE.-1) GOTO 100                                             0901.000
C                                                                       0902.000
C  NF=-1      INITIALIZE                                                0903.000
      NUMREF=0                                                          0904.000
      SUMRES=0.0                                                        0905.000
      SUMFS=0.0                                                         0906.000
      SUMRSQ=0.0                                                        0907.000
      SUMFSQ=0.0                                                        0908.000
      RETURN                                                            0909.000
C                                                                       0910.000
100   IF(NF.NE.0) GOTO 200                                              0911.000
C                                                                       0912.000
C  NF=0       ADD IN THIS REFLECTION                                    0913.000
      NUMREF=NUMREF+1                                                   0914.000
      SUMFS=SUMFS + FC                                                  0915.000
      SUMFSQ=SUMFSQ + FC*FC                                             0916.000
      RESID=ABS( FC-FM )                                                0917.000
      SUMRES=SUMRES + RESID                                             0918.000
      SUMRSQ=SUMRSQ + RESID*RESID                                       0919.000
      RETURN                                                            0920.000
C                                                                       0921.000
C  NF=1       END      REPORT RESULTS                                   0922.000
200   CONTINUE                                                          0923.000
      RM=SUMRES/SUMFS                                                   0924.000
      RMS=SQRT(SUMRSQ/SUMFSQ)                                           0925.000
      WRITE(LP,210) NUMREF, RM , RMS                                    0926.000
210   FORMAT(/' **ANALYSIS OF RESIDUALS**'/' For the',I6,               0927.000
     +' reflections'/' the relative residual was ',F6.4,                0928.000
     +'   and the relative rms residual was ',F6.4/)                    0929.000
      RETURN                                                            0930.000
      END                                                               0931.000
C                                                                       0932.000
C                                                                       0933.000
C                                                                       0934.000
C                                                                       0935.000
C                                                                       0936.000
      SUBROUTINE ABSPRI (MUMIN,MUMAX,THMAX)                             0937.000
C     *************************************                             0938.000
C                                                                       0939.000
C                                                                       0940.000
C   PRINT ABSORPTION SURFACE AND THETA CURVE                            0941.000
C                                                                       0942.000
C                                                                       0943.000
      LOGICAL TCOR                                                      0944.000
      REAL MU,MUMIN,MUMAX,A(18)                                         0945.000
      COMMON /NORMAL/ ANORM(44,44),ATB(44),ROWJ(88),TERMS(44)           0946.000
      COMMON /REFL/ FC,FO,FM,PHI(2),MU(2),SINTH,COSTH                   0947.000
      COMMON /DRVR/ LI,LD,LO,LF,LP,NCY,NSHAPE,NACOEF,NBCOEF,SFMU,TCOR   0948.000
C                                                                       0949.000
C                                                                       0950.000
      DATA DTOR/0.0174533/                                              0951.000
C                                                                       0952.000
C   PRINT PHI AT INTERVALS OF 10 DEGREES                                0953.000
C   AND MU AT INTERVALS OF 5 DEGREES                                    0954.000
      MAX=MUMAX/5. + 0.8                                                0955.000
      MIN=MUMIN/5. - 0.8                                                0956.000
      MRANGE=MAX-MIN + 1                                                0957.000
      MUS=MAX*5                                                         0958.000
      MUF=MIN*5                                                         0959.000
C                                                                       0960.000
C                                                                       0961.000
      DO 100 J=1,2                                                      0962.000
      NPS=(J-1)*180                                                     0963.000
      NPF=NPS+170                                                       0964.000
      WRITE(LP,50)NPS,NPF,MUF,MUS,(JJ,JJ=NPS,NPF,10)                    0965.000
50    FORMAT(//' ABSORPTION SURFACE FOR PHI =',I4,' --->',I4,           0966.000
     +' WITH MU RANGE =',I3,' --->',I3, / 6X,18(I3,2X) )                0967.000
C                                                                       0968.000
      DO 90 I=1,MRANGE                                                  0969.000
      II=(MAX - (I-1))*5                                                0970.000
      RMU=II                                                            0971.000
      IF(II.GT.0) RMU=SQRT(RMU)*SFMU                                    0972.000
      IF(II.LT.0) RMU=-SQRT(ABS(RMU))*SFMU                              0973.000
      MU(1)=RMU                                                         0974.000
      MU(2)=RMU                                                         0975.000
C                                                                       0976.000
      DO 80 K=1,18                                                      0977.000
      KK=(K-1)*10 + NPS                                                 0978.000
      PHI(1)=KK                                                         0979.000
      PHI(2)=KK                                                         0980.000
C                                                                       0981.000
      CALL ABSTOR                                                       0982.000
      CALL AFACT (PS1,PS2,PS3,ABSFAC)                                   0983.000
      A(K)=PS1*PS2                                                      0984.000
80    CONTINUE                                                          0985.000
C                                                                       0986.000
      WRITE(LP,200)II,A                                                 0987.000
200   FORMAT(I4,1X,18F5.2)                                              0988.000
90    CONTINUE                                                          0989.000
C                                                                       0990.000
100   CONTINUE                                                          0991.000
C                                                                       0992.000
      IF(.NOT.TCOR) GOTO 500                                            0993.000
C                                                                       0994.000
C   PRINT THETA DEPENDENT CURVE AT 5 DEGREE INTERVALS                   0995.000
      MAX=THMAX/5. + 0.8                                                0996.000
      MRANGE=MAX+1                                                      0997.000
      MEND=MAX*5                                                        0998.000
C                                                                       0999.000
      DO 350 I=1,MRANGE                                                 1000.000
      THET=(I-1)*5.                                                     1001.000
      SINTH=SIN(THET*DTOR)                                              1002.000
      COSTH=COS(THET*DTOR)                                              1003.000
      CALL ABSTOR                                                       1004.000
      CALL AFACT(PS1,PS2,PS3,ABSFAC)                                    1005.000
      A(I)=PS3                                                          1006.000
C                                                                       1007.000
350   CONTINUE                                                          1008.000
C                                                                       1009.000
      WRITE(LP,400) MEND,(I,I=0,MEND,5)                                 1010.000
400   FORMAT(/' Curve for THETA = 0 --->',I3,/3X,18(I7))                1011.000
      WRITE(LP,410) (A(I),I=1,MRANGE)                                   1012.000
410   FORMAT(5X,18F7.3)                                                 1013.000
C                                                                       1014.000
C                                                                       1015.000
500   RETURN                                                            1016.000
      END                                                               1017.000
C                                                                       1018.000
C                                                                       1019.000
C                                                                       1020.000
C                                                                       1021.000
C                                                                       1022.000
      SUBROUTINE ABSPLY (NREF)                                          1023.000
C     ************************                                          1024.000
C                                                                       1025.000
C                                                                       1026.000
C   APPLY ABSORPTION CORRECTION TO ALL DATA                             1027.000
C                                                                       1028.000
C                                                                       1029.000
      LOGICAL TCOR                                                      1030.000
      REAL MU,DUM(4)                                                    1031.000
      COMMON /NORMAL/ ANORM(44,44),ATB(44),ROWJ(88),TERMS(44)           1032.000
      COMMON /REFL/ FC,FO,FM,PHI(2),MU(2),SINTH,COSTH                   1033.000
      COMMON /DRVR/ LI,LD,LO,LF,LP,NCY,NSHAPE,NACOEF,NBCOEF,SFMU,TCOR   1034.000
      DATA AMIN,AMAX /100.,-100./                                       1035.000
C                                                                       1036.000
      NTOT=0                                                            1037.000
      ATOT=0.                                                           1038.000
      NSET=1                                                            1039.000
      SFMUSQ=SFMU*SFMU                                                  1040.000
C                                                                       1041.000
C   INITIALIZE RESIDUALS                                                1042.000
      NF=-1                                                             1043.000
      CALL RES(NF)                                                      1044.000
C                                                                       1045.000
C                                                                       1046.000
      NF=0                                                              1047.000
      REWIND LD                                                         1048.000
      REWIND LF                                                         1049.000
      WRITE(LP,1001)                                                    1050.000
 1001 FORMAT(//'   H  K  L   FCALC   FCORR  SIGCOR    FOBS    SIGF  NREF1051.000
     1',/)                                                              1052.000
      DO 100 NRECRD=1,NREF                                              1053.000
      READ(LD)JAA,JBB,FF,SIGI,AC,BC,SS                                  1054.000
      CALL KRAKJA(JAA,ISET,IH,IK,IL)                                    1055.000
      IUN=0                                                             1056.000
      FOBS=SQRT(FF)                                                     1057.000
      SIGF=0.5*SIGI/FOBS                                                1058.000
      FCALC=SQRT(AC*AC+BC*BC)                                           1059.000
      READ(LF) IUSED,FC,FO,PHI,MU,SINTH,COSTH                           1060.000
      NTOT=NTOT + 1                                                     1061.000
C                                                                       1062.000
      CALL ABSTOR                                                       1063.000
      CALL AFACT(PS1,PS2,PS3,A)                                         1064.000
      IF(A.GT.AMIN) GOTO 2                                              1065.000
      AMIN=A                                                            1066.000
      APHINP=PHI(1)                                                     1067.000
      APHINS=PHI(2)                                                     1068.000
      AMUNP=MU(1)*ABS(MU(1))/SFMUSQ                                     1069.000
      AMUNS=MU(2)*ABS(MU(2))/SFMUSQ                                     1070.000
2     IF(A.LT.AMAX) GOTO 4                                              1071.000
      AMAX=A                                                            1072.000
      APHIXP=PHI(1)                                                     1073.000
      APHIXS=PHI(2)                                                     1074.000
      AMUXP=MU(1)*ABS(MU(1))/SFMUSQ                                     1075.000
      AMUXS=MU(2)*ABS(MU(2))/SFMUSQ                                     1076.000
4     ATOT=ATOT + A                                                     1077.000
C                                                                       1078.000
C                                                                       1079.000
      IF(IUSED.NE.1) GOTO 10                                            1080.000
      FM=FO*A                                                           1081.000
      CALL RES(NF)                                                      1082.000
C                                                                       1083.000
10    XOBS=FOBS*A                                                       1084.000
      XIGF=SIGF*A                                                       1085.000
      FF=XOBS*XOBS                                                      1086.000
      SIGI=2.0*XOBS*XIGF                                                1087.000
      WRITE (LO)JAA,JBB,FF,SIGI                                         1088.000
      IF(NRECRD.LE.20.OR.MOD(NRECRD,100).EQ.0)WRITE(LP,201)IH,IK,IL,FCAL1089.000
     1C,XOBS,XIGF,FOBS,SIGF,NRECRD                                      1090.000
  201 FORMAT(1H ,3I3,5F8.1,I6)                                          1091.000
C30   WRITE(LO,35) IH,IK,IL,FOBS,SIGF,NSET                              1092.000
C35   FORMAT(3I4,2F8.2,I4)                                              1093.000
C                                                                       1094.000
100   CONTINUE                                                          1095.000
C                                                                       1096.000
C   END FILE WITH A BLANK LINE FOR SHELX                                1097.000
C     WRITE(LO,40)                                                      1098.000
C40   FORMAT(32X)                                                       1099.000
C                                                                       1100.000
C                                                                       1101.000
C   REPORT RESIDUALS AFTER CORRECTION                                   1102.000
      WRITE(LP,300)                                                     1103.000
300   FORMAT(//' RESIDUALS AFTER CORRECTION')                           1104.000
      NF=1                                                              1105.000
      CALL RES(NF)                                                      1106.000
C                                                                       1107.000
C                                                                       1108.000
C   REPORT CORRECTION SUMMARY                                           1109.000
      AVER=ATOT/NTOT                                                    1110.000
      WRITE(LP,200)NTOT,AMIN,APHINP,AMUNP,APHINS,AMUNS,AMAX,            1111.000
     +APHIXP,AMUXP,APHIXS,AMUXS,AVER                                    1112.000
200   FORMAT(//I6,' Reflections corrected'                              1113.000
     +       /' Minimum absorption correction =',F6.3,                  1114.000
     +  5X,'at  PHI(P)=',F5.1,' MU(P)=',F5.1,                           1115.000
     +  '   PHI(S)=',F5.1,' MU(S)=',F5.1,                               1116.000
     +       /' Maximum absorption correction =',F6.3,                  1117.000
     +  5X,'at  PHI(P)=',F5.1,' MU(P)=',F5.1,                           1118.000
     +  '   PHI(S)=',F5.1,' MU(S)=',F5.1,                               1119.000
     +       /' Average absorption correction =',F6.3)                  1120.000
C                                                                       1121.000
C                                                                       1122.000
      RETURN                                                            1123.000
      END                                                               1124.000
C                                                                       1125.000
C                                                                       1126.000
C                                                                       1127.000
C                                                                       1128.000
C                                                                       1129.000
C                                                                       1130.000
C                                                                       1131.000
      SUBROUTINE SOLVEN                                                 1132.000
C     *****************                                                 1133.000
C                                                                       1134.000
C                                                                       1135.000
      LOGICAL TCOR                                                      1136.000
      REAL VEC(44,44),EIG(44),UM1(44),WORK(220),RDASH(44),SDASH(44)     1137.000
      COMMON /NORMAL/ ANORM(44,44),ATB(44),ROWJ(88),TERMS(44)           1138.000
      COMMON /SOLVE/ NSIZE,FILTER,FILSTA,FILDEC,FILEND                  1139.000
      COMMON /DRVR/ LI,LD,LO,LF,LP,NCY,NSHAPE,NACOEF,NBCOEF,SFMU,TCOR   1140.000
C   FOLLOWING EQUIVALENCING TO SAVE SOME SPACE                          1141.000
      EQUIVALENCE (WORK(1),RDASH(1)),(WORK(45),SDASH(1))                1142.000
      EQUIVALENCE (ROWJ(1),EIG(1)),(ROWJ(45),UM1(1))                    1143.000
C                                                                       1144.000
C Pre-scale normal matrix.                                              1145.000
C Invert normal matrix using HARWELL s/r for eigenvalue inversion.      1146.000
C                                                                       1147.000
C Get inverse diagonal scaling matrix UM1 (U**-1) from normal matrix.   1148.000
C The elements are SQRT(diagonal element of normal matrix)              1149.000
C                                                                       1150.000
C                                                                       1151.000
      AINFL=FILSTA - (NCY-1)*FILDEC                                     1152.000
      IF(AINFL.LT.FILEND)AINFL=FILEND                                   1153.000
C      WRITE(LP,1030)FILTER,AINFL                                       1154.000
C1030  FORMAT(/' Eigenvalue filtering; filter value=',F10.6,            1155.000
C     +',  inflation constant=',F10.6)                                  1156.000
C                                                                       1157.000
      DO 10 I=1,NSHAPE                                                  1158.000
      UM1(I)=1.0                                                        1159.000
      IF(ANORM(I,I).NE.0.0) UM1(I)=1.0/SQRT(ANORM(I,I))                 1160.000
10    CONTINUE                                                          1161.000
C                                                                       1162.000
C Scale matrix and vector(ATB) by UM1                                   1163.000
C                                                                       1164.000
      DO 20 J=1,NSHAPE                                                  1165.000
      UM1J=UM1(J)                                                       1166.000
      ATB(J)=ATB(J)*UM1J                                                1167.000
      DO 30 K=1,J                                                       1168.000
      ANORM(J,K)=ANORM(J,K)*UM1(K)*UM1J                                 1169.000
30    CONTINUE                                                          1170.000
20    CONTINUE                                                          1171.000
C                                                                       1172.000
C Solve normal equations, EA06C returns eigenvectors and eigenvalues    1173.000
C                                                                       1174.000
      CALL EA06C(ANORM,EIG,VEC,NSHAPE,NSIZE,NSIZE,WORK)                 1175.000
C                                                                       1176.000
C Write eigenvalues                                                     1177.000
C                                                                       1178.000
C      WRITE(LP,1000)                                                   1179.000
C1000  FORMAT(/,' Eigenvalues')                                         1180.000
C      WRITE(LP,1010)(EIG(I),I=1,NSHAPE)                                1181.000
C1010  FORMAT(10F8.4)                                                   1182.000
C                                                                       1183.000
C To diagonalize, premultiply by column eigenvectors                    1184.000
C  VEC is matrix of eignvectors.                                        1185.000
C                                                                       1186.000
      DO 40 I=1,NSHAPE                                                  1187.000
      SUM=0.0                                                           1188.000
      DO 50 J=1,NSHAPE                                                  1189.000
      SUM=SUM+VEC(J,I)*ATB(J)                                           1190.000
50    CONTINUE                                                          1191.000
      RDASH(I)=SUM                                                      1192.000
40    CONTINUE                                                          1193.000
C                                                                       1194.000
C Solve diagonalization with inflated eigenvalue weighting & filtering  1195.000
C                                                                       1196.000
      DO 60 I=1,NSHAPE                                                  1197.000
      SDASH(I)=0.0                                                      1198.000
      IF(EIG(I).GT.FILTER)SDASH(I)=RDASH(I)/(EIG(I)+AINFL)              1199.000
60    CONTINUE                                                          1200.000
C                                                                       1201.000
C Rotate shifts back and rescale using UM1                              1202.000
C                                                                       1203.000
      DO 70 J=1,NSHAPE                                                  1204.000
      SUM=0.0                                                           1205.000
      DO 80 I=1,NSHAPE                                                  1206.000
      SUM=SUM+VEC(J,I)*SDASH(I)                                         1207.000
80    CONTINUE                                                          1208.000
      ATB(J)=SUM*UM1(J)                                                 1209.000
70    CONTINUE                                                          1210.000
C                                                                       1211.000
C                                                                       1212.000
C Write coefficient shifts                                              1213.000
C                                                                       1214.000
      WRITE(LP,1040)                                                    1215.000
1040  FORMAT(/,' Coefficient shifts')                                   1216.000
      WRITE(LP,1050)(ATB(I),I=1,NSHAPE)                                 1217.000
1050  FORMAT(1X,10F8.5)                                                 1218.000
C                                                                       1219.000
C Apply shifts to original defining variables                           1220.000
      DO 90 J=1,NSHAPE                                                  1221.000
      TERMS(J)=TERMS(J)+ATB(J)                                          1222.000
90    CONTINUE                                                          1223.000
C                                                                       1224.000
C Write coefficients of absorption factor expansion                     1225.000
C                                                                       1226.000
C      WRITE(LP,1020)                                                   1227.000
C1020  FORMAT(/,' New coefficients')                                    1228.000
C      WRITE(LP,*)(TERMS(I),I=1,NSHAPE)                                 1229.000
      RETURN                                                            1230.000
      END                                                               1231.000
C                                                                       1232.000
C                                                                       1233.000
C                                                                       1234.000
C                                                                       1235.000
C                                                                       1236.000
C                                                                       1237.000
C                                                                       1238.000
C                                                                       1239.000
      SUBROUTINE ORTOCELL                                               1240.000
C     ***************                                                   1241.000
C                                                                       1242.000
C   GENERAL SUBROUTINE TO CALCULATE CELL PARAMETERS FROM THE            1243.000
C   ORIENTATION MATRIX                                                  1244.000
C                                                                       1245.000
      LOGICAL TCOR                                                      1246.000
      COMMON /DRVR/ LI,LD,LO,LF,LP,NCY,NSHAPE,NACOEF,NBCOEF,SFMU,TCOR   1247.000
      COMMON/ORIENT/TWOTHZ,OMZ,PHIZ,PLUS(3),IPOS,IDIFRT,ORIENT(3,3),DST 1248.000
      COMMON /CELLP/ TA,TB,TC,ZA,ZB,ZC                                  1249.000
      REAL DIRECT(3,3),W(3),X(3),Y(3),Z(3)                              1250.000
      DATA RAD/57.2957795/                                              1251.000
C                                                                       1252.000
      CALL AMI3(ORIENT,DIRECT,DET)                                      1253.000
C                                                                       1254.000
      DO 11 I=1,3                                                       1255.000
      W(I)=0.                                                           1256.000
      Y(I)=0.                                                           1257.000
      DO 12 J=1,3                                                       1258.000
      W(I)=W(I)+ORIENT(J,I)**2                                          1259.000
      Y(I)=Y(I)+DIRECT(I,J)**2                                          1260.000
12    CONTINUE                                                          1261.000
      Y(I)=SQRT(Y(I))                                                   1262.000
11    CONTINUE                                                          1263.000
C                                                                       1264.000
      DO 13 I=1,3                                                       1265.000
      J=I+1                                                             1266.000
      K=J+1                                                             1267.000
      IF(J.GT.3)J=J-3                                                   1268.000
      IF(K.GT.3)K=K-3                                                   1269.000
      X(I)=0.                                                           1270.000
      Z(I)=0.                                                           1271.000
      DO 14 L=1,3                                                       1272.000
      X(I)=X(I) + ORIENT(L,J)*ORIENT(L,K)                               1273.000
      Z(I)=Z(I) + DIRECT(J,L)*DIRECT(K,L)                               1274.000
14    CONTINUE                                                          1275.000
      X(I)=2*X(I)                                                       1276.000
      COSZ=Z(I)/Y(J)/Y(K)                                               1277.000
      SINZ=SQRT(1. - COSZ*COSZ)                                         1278.000
      Z(I)=ATAN(SINZ/COSZ)*RAD                                          1279.000
      IF(Z(I).LT.0.)Z(I)=Z(I)+180.                                      1280.000
13    CONTINUE                                                          1281.000
C                                                                       1282.000
      TA=W(1)                                                           1283.000
      TB=W(2)                                                           1284.000
      TC=W(3)                                                           1285.000
      ZA=X(1)                                                           1286.000
      ZB=X(2)                                                           1287.000
      ZC=X(3)                                                           1288.000
C                                                                       1289.000
      WRITE(LP,100)((ORIENT(I,J),J=1,3),I=1,3),Y,Z                      1290.000
100   FORMAT(/' Orientation Matrix:'/,3(5X,3F10.6/),/                   1291.000
     +' Unit Cell Parameters from Orientation Matrix:'/,                1292.000
     +                          '  a=',F7.3,'  b=',F7.3,'  c=',F7.3,    1293.000
     +'   alpha=',F7.2,'  beta=',F7.2,'  gamma=',F7.2/)                 1294.000
C                                                                       1295.000
      RETURN                                                            1296.000
      END                                                               1297.000
C                                                                       1298.000
C                                                                       1299.000
C                                                                       1300.000
C                                                                       1301.000
C                                                                       1302.000
C                                                                       1303.000
      SUBROUTINE AMI3 (A,B,DET)                                         1304.000
C     *************************                                         1305.000
C                                                                       1306.000
C.... INVERTS ASYMMETRIC 3X3 MATRIX                                     1307.000
C                                                                       1308.000
      REAL A(3,3),B(3,3)                                                1309.000
C                                                                       1310.000
      B(1,1)=A(2,2)*A(3,3)-A(3,2)*A(2,3)                                1311.000
      B(2,1)=A(2,1)*A(3,3)-A(3,1)*A(2,3)                                1312.000
      B(3,1)=A(2,1)*A(3,2)-A(3,1)*A(2,2)                                1313.000
      DET=A(1,1)*B(1,1)-A(1,2)*B(2,1)+A(1,3)*B(3,1)                     1314.000
      IF (DET) 12,16,12                                                 1315.000
12    B(1,2)=A(1,2)*A(3,3)-A(3,2)*A(1,3)                                1316.000
      B(2,2)=A(1,1)*A(3,3)-A(3,1)*A(1,3)                                1317.000
      B(3,2)=A(1,1)*A(3,2)-A(3,1)*A(1,2)                                1318.000
      B(1,3)=A(1,2)*A(2,3)-A(2,2)*A(1,3)                                1319.000
      B(2,3)=A(1,1)*A(2,3)-A(2,1)*A(1,3)                                1320.000
      B(3,3)=A(1,1)*A(2,2)-A(2,1)*A(1,2)                                1321.000
      DO 14 I=1,3                                                       1322.000
      DO 14 J=1,3                                                       1323.000
      B(I,J)=(-1.)**(I+J)*B(I,J)/DET                                    1324.000
14    CONTINUE                                                          1325.000
C                                                                       1326.000
16    RETURN                                                            1327.000
      END                                                               1328.000
C                                                                       1329.000
C                                                                       1330.000
C                                                                       1331.000
C                                                                       1332.000
C                                                                       1333.000
C                                                                       1334.000
      SUBROUTINE ANGLES (OMEGA,PHI,CHI,X,IERR)                          1335.000
C     ****************************************                          1336.000
C                                                                       1337.000
C.... CALCULATES DIFFRACTOMETER SETTING ANGLES GIVEN OMEGA              1338.000
C                                                                       1339.000
      COMMON/ORIENT/TWOTHZ,OMZ,PHIZ,PLUS(3),IPOS,IDIFRT,ORIENT(3,3),DST 1340.000
      DIMENSION X(3),E(3)                                               1341.000
      DATA RAD/57.2957795/                                              1342.000
      DATA PI/3.1415926563/                                             1343.000
C                                                                       1344.000
      DO 14 I=1,3                                                       1345.000
      E(I)=0.                                                           1346.000
      DO 12 J=1,3                                                       1347.000
   12 E(I)=E(I)+ORIENT(I,J)*X(J)                                        1348.000
   14 E(I)=E(I)/DST                                                     1349.000
      IERR=0                                                            1350.000
C                                                                       1351.000
      IF(IPOS.NE.1) GO TO 36                                            1352.000
C                                                                       1353.000
C     CHI=90, DATA COLLECTED IN PARALLEL POSITION.                      1354.000
C                                                                       1355.000
      CHI=90.                                                           1356.000
      IF(ABS(E(3)).GT.1.) GO TO 64                                      1357.000
   16 OMEGA=ACOS(ABS(E(3)))                                             1358.000
      IF(-E(3).GE.0.) GO TO 20                                          1359.000
   18 OMEGA=PI-OMEGA                                                    1360.000
   20 CPHI=E(1)/SIN(OMEGA)                                              1361.000
      SPHI=-E(2)/SIN(OMEGA)                                             1362.000
      IF (CPHI.GT.1..OR.SPHI.GT.1.) GO TO 64                            1363.000
      CALL ARCTAN (PHI,CPHI,SPHI)                                       1364.000
   22 IF (E(1).LE.0.) GO TO 26                                          1365.000
   24 IF(E(2)) 30,30,28                                                 1366.000
   26 IF (E(2))34,34,32                                                 1367.000
   28 PHI=360. - ACOS(CPHI)*RAD                                         1368.000
      GO TO 35                                                          1369.000
   30 PHI=ACOS(CPHI)*RAD                                                1370.000
      GO TO 35                                                          1371.000
C*****THIS MAY BE IN ERROR BECAUSE OF THE WAY ARCTAN WORKS:             1372.000
   32 PHI= 180. + PHI                                                   1373.000
      GO TO 35                                                          1374.000
   34 PHI=180. - PHI                                                    1375.000
   35 OMEGA=OMEGA*RAD                                                   1376.000
      GO TO 62                                                          1377.000
C                                                                       1378.000
   36 IF(IDIFRT.NE.3) GO TO 80                                          1379.000
C                                                                       1380.000
C    DATA COLLECTED IN BISECTING POSITION.  (PICKER GEOMETRY).          1381.000
C                                                                       1382.000
      C=ABS(COS(OMEGA/RAD))                                             1383.000
      IF (C-0.000001) 38,38,42                                          1384.000
   38 IF(ABS(E(3))-0.000001) 40,40,64                                   1385.000
   40 CHI=0.0                                                           1386.000
      CPHI=E(1)                                                         1387.000
      SPHI=-E(2)                                                        1388.000
      GO TO 58                                                          1389.000
   42 SCHI=-E(3)/C                                                      1390.000
      IF(ABS(SCHI).GT.1.0) SCHI=SIGN(1.0,SCHI)                          1391.000
      A=SIN(OMEGA/RAD)                                                  1392.000
      CCHI=SQRT(1.0-SCHI**2)                                            1393.000
      CALL ARCTAN (CHI,CCHI,SCHI)                                       1394.000
      B=C*CCHI                                                          1395.000
      DENO=A**2+B**2                                                    1396.000
      IF (DENO-0.000001) 46,46,48                                       1397.000
   46 PHI=0.0                                                           1398.000
      GO TO 62                                                          1399.000
   48 CPHI=(A*E(1)+B*E(2))/DENO                                         1400.000
      IF(ABS(CPHI).GT.1.0) CPHI=SIGN(1.0,CPHI)                          1401.000
      IF (B-0.0000001) 52,52,54                                         1402.000
   52 SPHI=-E(2)/A                                                      1403.000
      GO TO 56                                                          1404.000
   54 SPHI=(E(1)-A*CPHI)/B                                              1405.000
   56 IF(ABS(SPHI).GT.1.0) SPHI=SIGN(1.0,SPHI)                          1406.000
   58 CALL ARCTAN (PHI,CPHI,SPHI)                                       1407.000
   62 OMEGA=OMZ+OMEGA/PLUS(1)                                           1408.000
      GO TO 90                                                          1409.000
   64 IERR=1                                                            1410.000
      RETURN                                                            1411.000
C                                                                       1412.000
C     DATA COLLECTED IN BISECTING POSITION.  (NONIUS & SYNTEX GEOMETRY) 1413.000
C                                                                       1414.000
   80 OMEGA=0.0                                                         1415.000
      C=1.0                                                             1416.000
      A=0.0                                                             1417.000
      SCHI=-E(3)                                                        1418.000
      IF(ABS(SCHI).GT.1.0) SCHI=SIGN(1.0,SCHI)                          1419.000
C                                                                       1420.000
C     COS(CHI) ALWAYS POSITIVE IF COLLECTING IN LOWER                   1421.000
C     HEMISPHERE  (-90<CHI<+90).                                        1422.000
C                                                                       1423.000
      CCHI=SQRT(1.0-SCHI**2)                                            1424.000
      CALL ARCTAN (CHI,CCHI,SCHI)                                       1425.000
      IF(CCHI-0.0001) 46,46,84                                          1426.000
   84 CPHI=E(2)/CCHI                                                    1427.000
      IF(ABS(CPHI).GT.1.0) CPHI=SIGN(1.0,CPHI)                          1428.000
      SPHI=E(1)/CCHI                                                    1429.000
      IF(ABS(SPHI).GT.1.0) SPHI=SIGN(1.0,SPHI)                          1430.000
      CALL ARCTAN (PHI,CPHI,SPHI)                                       1431.000
   90 PHI=PHIZ+PHI/PLUS(2)                                              1432.000
      CHI=CHI/PLUS(3)                                                   1433.000
      RETURN                                                            1434.000
      END                                                               1435.000
C                                                                       1436.000
C                                                                       1437.000
C                                                                       1438.000
C                                                                       1439.000
C                                                                       1440.000
      SUBROUTINE ARCTAN (A,COSA,SINA)                                   1441.000
C     *******************************                                   1442.000
C                                                                       1443.000
      DATA RAD/57.2957795/                                              1444.000
C                                                                       1445.000
C   CALCULATES ARC TAN IN DEGREES.                                      1446.000
C                                                                       1447.000
      IF(COSA) 10,20,10                                                 1448.000
10    A=ATAN(SINA/COSA)*RAD                                             1449.000
      IF(COSA) 30,30,40                                                 1450.000
20    A=90.0                                                            1451.000
      IF(SINA) 30,30,40                                                 1452.000
30    A=A+180.0                                                         1453.000
40    RETURN                                                            1454.000
      END                                                               1455.000
C                                                                       1456.000
C                                                                       1457.000
C                                                                       1458.000
C                                                                       1459.000
      SUBROUTINE NORM (A)                                               1460.000
C     *******************                                               1461.000
C                                                                       1462.000
      DO 14 I=1,10                                                      1463.000
      IF (A) 12,16,16                                                   1464.000
12    A=A+360.                                                          1465.000
14    CONTINUE                                                          1466.000
16    A=AMOD(A,360.0)                                                   1467.000
      RETURN                                                            1468.000
      END                                                               1469.000
C                                                                       1470.000
C                                                                       1471.000
C                                                                       1472.000
C                                                                       1473.000
C                                                                       1474.000
C                                                                       1475.000
C                                                                       1476.000
C                                                                       1477.000
C                                                                       1478.000
      SUBROUTINE EA06C(A,VALUE,VECTOR,M,IA,IV,W)                        1479.000
C     ******************************************                        1480.000
C                                                                       1481.000
C**  18/03/70 LAST LIBRARY UPDATE                                       1482.000
C                                                                       1483.000
C  Given a real mXm symmetric matrix A = ?aij? this routine finds all it1484.000
C  eigenvalues (lambda)i i=1,2,.....,m and eigenvectors xj i=1,2,...,m. 1485.000
C  i.e. finds the non-trivial solutions of Ax=(lambda)x                 1486.000
C                                                                       1487.000
C  Argument list                                                        1488.000
C  -------------                                                        1489.000
C                                                                       1490.000
C  SUBROUTINE EA06C(A,VALUE,X,M,IA,IX,W)                                1491.000
C                                                                       1492.000
C  Arguments set by the user -                                          1493.000
C                                                                       1494.000
C  A  is a REAL 2-dimensional array with first dimension IA. The user mu1495.000
C     store the lower triangle of matrix A into the lower triangle of th1496.000
C     array A. i.e. put a(ij) into A(I,J) for i>=j and I>=J.            1497.000
C                                                                       1498.000
C  M  is an INTEGER and should be set to the order m of the matrix      1499.000
C                                                                       1500.000
C  IA is an INTEGER set to the first dimension of the array A, i.e. if t1501.000
C     allocation for the array A was specified by                       1502.000
C                                                                       1503.000
C            DIMENSION A(100,50)                                        1504.000
C                                                                       1505.000
C     then IA would be set to 100                                       1506.000
C                                                                       1507.000
C  IX is an INTEGER set to the first dimension of the 2-dimensional arra1508.000
C                                                                       1509.000
C                                                                       1510.000
C  Arguments set by the routine -                                       1511.000
C                                                                       1512.000
C  A  the space above the diagonal of the array A is used by the routine1513.000
C     work space, the lower triangle set by the user will remain unchang1514.000
C     on return.                                                        1515.000
C                                                                       1516.000
C  VALUE is a REAL array in which the routine puts the eigenvalues (lamb1517.000
C     i=1,2,...,m. These are not necessarily in any order.              1518.000
C                                                                       1519.000
C  X  is a REAL 2-dimensional array, with first dimension IX, containing1520.000
C     eigenvectors. The components of the eigenvector x(i) corresponding1521.000
C     eigenvalue (lambda)i (in VALUE(I)) are placed in X(J,I) J=1,2,...,1522.000
C     The eigenvectors are normalized so that xT(i)x(i)=1 i=1,2,...,m.  1523.000
C                                                                       1524.000
C  W  is a REAL array used by the routine for work space. It must have  1525.000
C     dimension at least 5*M.                                           1526.000
C                                                                       1527.000
C                                    ------------                       1528.000
C                                                                       1529.000
C                                                                       1530.000
      REAL A(IA,M),VALUE(M),VECTOR(IV,M),W(1)                           1531.000
C                                                                       1532.000
      M1=M+1                                                            1533.000
      W(1)=A(1,1)                                                       1534.000
      IF(M-2)30,10,20                                                   1535.000
   10 W(2)=A(2,2)                                                       1536.000
      W(4)=A(2,1)                                                       1537.000
      GO TO 30                                                          1538.000
C                                                                       1539.000
   20 CALL MC04B(A,W,W(M1),M,IA,W(M+M1))                                1540.000
C                                                                       1541.000
   30 CALL EA08C(W,W(M1),VALUE,VECTOR,M,IV,W(M+M1))                     1542.000
C                                                                       1543.000
      IF(M.LE.2)RETURN                                                  1544.000
      DO 70 L=1,M                                                       1545.000
      DO 70 II=3,M                                                      1546.000
      I=M-II+1                                                          1547.000
      IF(W(M1+I))40,70,40                                               1548.000
   40 PP=0.                                                             1549.000
      I1=I+1                                                            1550.000
      DO 50 K=I1,M                                                      1551.000
   50 PP=PP+A(I,K)*VECTOR(K,L)                                          1552.000
      PP=PP/(A(I,I+1)*W(M1+I))                                          1553.000
      DO 60 K=I1,M                                                      1554.000
   60 VECTOR(K,L)=VECTOR(K,L)+PP*A(I,K)                                 1555.000
   70 CONTINUE                                                          1556.000
      RETURN                                                            1557.000
      END                                                               1558.000
C                                                                       1559.000
C                                                                       1560.000
C***  *  *  *  *  *  *  *  *  *  *  *  *  *                             1561.000
C                                                                       1562.000
C                                                                       1563.000
      SUBROUTINE MC04B(A,ALPHA,BETA,M,IA,Q)                             1564.000
C     *************************************                             1565.000
C                                                                       1566.000
      DIMENSION A(IA,1),ALPHA(1),BETA(1),Q(1)                           1567.000
C                                                                       1568.000
      ALPHA(1)=A(1,1)                                                   1569.000
      DO 20 J=2,M                                                       1570.000
      J1=J-1                                                            1571.000
      DO 10 I=1,J1                                                      1572.000
      A(I,J)=A(J,I)                                                     1573.000
   10 CONTINUE                                                          1574.000
      ALPHA(J)=A(J,J)                                                   1575.000
   20 CONTINUE                                                          1576.000
      M1=M-1                                                            1577.000
      M2=M-2                                                            1578.000
      DO 170 I=1,M2                                                     1579.000
      PP=0.0                                                            1580.000
      I1=I+1                                                            1581.000
      DO 30 J=I1,M                                                      1582.000
      PP=PP+A(I,J)**2                                                   1583.000
   30 CONTINUE                                                          1584.000
      PP1=SQRT(PP)                                                      1585.000
      IF(A(I,I+1))50,40,40                                              1586.000
   40 BETA(I+1)=-PP1                                                    1587.000
      GO TO 60                                                          1588.000
   50 BETA(I+1)=PP1                                                     1589.000
   60 IF(PP)170,170,70                                                  1590.000
   70 H=PP-BETA(I+1)*A(I,I+1)                                           1591.000
      A(I,I+1)=A(I,I+1)-BETA(I+1)                                       1592.000
      DO 120 KI=I1,M                                                    1593.000
      QJ=0.0                                                            1594.000
      DO 80 KJ=I1,KI                                                    1595.000
      QJ=QJ+A(KJ,KI)*A(I,KJ)                                            1596.000
   80 CONTINUE                                                          1597.000
      IF(KI-M)90,110,110                                                1598.000
   90 I2=KI+1                                                           1599.000
      DO 100 KJ=I2,M                                                    1600.000
      QJ=QJ+A(KI,KJ)*A(I,KJ)                                            1601.000
  100 CONTINUE                                                          1602.000
  110 Q(KI)=QJ/H                                                        1603.000
  120 CONTINUE                                                          1604.000
      BIGK=0.0                                                          1605.000
      DO 130 KJ=I1,M                                                    1606.000
      BIGK=BIGK+A(I,KJ)*Q(KJ)                                           1607.000
  130 CONTINUE                                                          1608.000
      BIGK=BIGK/(2.0*H)                                                 1609.000
      DO 140 KJ=I1,M                                                    1610.000
      Q(KJ)=Q(KJ)-BIGK*A(I,KJ)                                          1611.000
  140 CONTINUE                                                          1612.000
      DO 160 KI=I1,M                                                    1613.000
      DO 150 KJ=KI,M                                                    1614.000
      A(KI,KJ)=A(KI,KJ)-Q(KI)*A(I,KJ)-Q(KJ)*A(I,KI)                     1615.000
  150 CONTINUE                                                          1616.000
  160 CONTINUE                                                          1617.000
  170 CONTINUE                                                          1618.000
      DO 180 I=2,M                                                      1619.000
      H=ALPHA(I)                                                        1620.000
      ALPHA(I)=A(I,I)                                                   1621.000
      A(I,I)=H                                                          1622.000
  180 CONTINUE                                                          1623.000
      BETA(M)=A(M-1,M)                                                  1624.000
      RETURN                                                            1625.000
      END                                                               1626.000
C                                                                       1627.000
C                                                                       1628.000
C  *  *  *  *  *  *  *  *  *  *  *  *  *  *                             1629.000
C                                                                       1630.000
C                                                                       1631.000
      SUBROUTINE EA08C(A,B,VALUE,VEC,M,IV,W)                            1632.000
C     **************************************                            1633.000
C                                                                       1634.000
C    19/03/70 LAST LIBRARY UPDATE                                       1635.000
C                                                                       1636.000
C  THIS USES QR ITERATION TO FIND THE EIGENVALUES AND EIGENVECTORS      1637.000
C  OF THE SYMMETRIC TRIDIAGONAL MATRIX WHOSE DIAGONAL ELEMENTS ARE      1638.000
C  A(I),I=1,M AND OFF-DIAGONAL ELEMENTS ARE B(I),I=2,M.  THE ARRAY      1639.000
C  W IS USED FOR WORKSPACE AND MUST HAVE DIMENSION AT LEAST 2*M.        1640.000
C  WE TREAT VEC AS IF IT HAD DIMENSIONS (IV,M).                         1641.000
C                                                                       1642.000
C                                                                       1643.000
      REAL A(M),B(M),VALUE(M),VEC(1),W(1)                               1644.000
      DATA EPS/  1E-5/,A34/0./                                          1645.000
C                                                                       1646.000
      CALL EA09C(A,B,W(M+1),M,W)                                        1647.000
C                                                                       1648.000
C  SET VEC TO THE IDENTITY MATRIX.                                      1649.000
      DO 20 I=1,M                                                       1650.000
      VALUE(I)=A(I)                                                     1651.000
      W(I)=B(I)                                                         1652.000
      K=(I-1)*IV+1                                                      1653.000
      L=K+M-1                                                           1654.000
      DO 10 J=K,L                                                       1655.000
   10 VEC(J)=0.                                                         1656.000
   20 VEC(K+I-1)=1.                                                     1657.000
      ITER=0                                                            1658.000
      IF(M.EQ.1)RETURN                                                  1659.000
      N2=M                                                              1660.000
C  EACH QR ITERATION IS PERFORMED OF ROWS AND COLUMNS N1 TO N2          1661.000
   30 DO 40 II=2,N2                                                     1662.000
      N1=2+N2-II                                                        1663.000
      IF(ABS(W(N1)).LE.(ABS(VALUE(N1-1))+ABS(VALUE(N1)))*EPS)GO TO 50   1664.000
   40 CONTINUE                                                          1665.000
      N1=1                                                              1666.000
   50 IF(N2.NE.N1) GO TO 60                                             1667.000
      N2=N2-1                                                           1668.000
      IF(N2-1)90,90,30                                                  1669.000
   60 ROOT=W(M+N2)                                                      1670.000
      ITER=ITER+1                                                       1671.000
      N2M1=N2-1                                                         1672.000
      A22=VALUE(N1)                                                     1673.000
      A12=A22-ROOT                                                      1674.000
      A23=W(N1+1)                                                       1675.000
      A13=A23                                                           1676.000
      DO 80 I=N1,N2M1                                                   1677.000
      A33=VALUE(I+1)                                                    1678.000
      IF(I.NE.N2M1)A34=W  (I+2)                                         1679.000
      S=SIGN(SQRT(A12*A12+A13*A13),A12)                                 1680.000
      SI=A13/S                                                          1681.000
      CO=A12/S                                                          1682.000
      JK=I*IV+1                                                         1683.000
      J1=JK-IV                                                          1684.000
      J2=J1+MIN0(M,I+ITER)-1                                            1685.000
      DO 70 JI=J1,J2                                                    1686.000
      V1=VEC(JI)                                                        1687.000
      V2=VEC(JK)                                                        1688.000
      VEC(JI)=V1*CO+V2*SI                                               1689.000
      VEC(JK)=V2*CO-V1*SI                                               1690.000
   70 JK=JK+1                                                           1691.000
      IF(I.NE.N1)  W(I)=S                                               1692.000
      A11=CO*A22+SI*A23                                                 1693.000
      A12=CO*A23+SI*A33                                                 1694.000
      A13=SI*A34                                                        1695.000
      A21=CO*A23-SI*A22                                                 1696.000
      A22=CO*A33-SI*A23                                                 1697.000
      A23=CO*A34                                                        1698.000
      VALUE(I)=A11*CO+A12*SI                                            1699.000
      A12=-A11*SI+A12*CO                                                1700.000
      W(I+1)=A12                                                        1701.000
   80 A22=A22*CO-A21*SI                                                 1702.000
      VALUE(N2)=A22                                                     1703.000
      GO TO 30                                                          1704.000
C  RAYLEIGH QUOTIENT                                                    1705.000
   90 DO 110 J=1,M                                                      1706.000
      K=(J-1)*IV                                                        1707.000
      XX=VEC(K+1)**2                                                    1708.000
      XAX=XX*A(1)                                                       1709.000
      DO 100 I=2,M                                                      1710.000
      XX=XX+VEC(K+I)**2                                                 1711.000
  100 XAX=XAX+VEC(K+I)*(2.*B(I)*VEC(K+I-1)+A(I)*VEC(K+I))               1712.000
  110 VALUE(J)=XAX/XX                                                   1713.000
      RETURN                                                            1714.000
      END                                                               1715.000
C                                                                       1716.000
C                                                                       1717.000
C  *  *  *  *  *  *  *  *  *  *  *  *  *  *                             1718.000
C                                                                       1719.000
C                                                                       1720.000
      SUBROUTINE EA09C(A,B,VALUE,M,OFF)                                 1721.000
C     *********************************                                 1722.000
C                                                                       1723.000
C    18/03/70 LAST LIBRARY UPDATE                                       1724.000
C                                                                       1725.000
C                                                                       1726.000
      REAL A(M),B(M),VALUE(M),OFF(M)                                    1727.000
      DATA A34/0./,EPS/0.6E-7/                                          1728.000
C                                                                       1729.000
      VALUE(1)=A(1)                                                     1730.000
      IF(M.EQ.1)RETURN                                                  1731.000
      DO 10 I=2,M                                                       1732.000
      VALUE(I)=A(I)                                                     1733.000
   10 OFF(I)=B(I)                                                       1734.000
C  EACH QR ITERATION IS PERFORMED OF ROWS AND COLUMNS N1 TO N2          1735.000
      N2=M                                                              1736.000
   20 IF(N2.LE.1)RETURN                                                 1737.000
      DO 30 II=2,N2                                                     1738.000
      N1=2+N2-II                                                        1739.000
      IF(ABS(OFF(N1)).LE.(ABS(VALUE(N1-1))+ABS(VALUE(N1)))*EPS)GO TO 40 1740.000
   30 CONTINUE                                                          1741.000
      N1=1                                                              1742.000
   40 IF(N2.NE.N1) GO TO 50                                             1743.000
      N2=N2-1                                                           1744.000
      GO TO 20                                                          1745.000
C  ROOT  IS THE EIGENVALUE OF THE BOTTOM 2*2 MATRIX THAT IS NEAREST     1746.000
C  TO THE LAST MATRIX ELEMENT AND IS USED TO ACCELERATE THE             1747.000
C  CONVERGENCE                                                          1748.000
   50 BB=(VALUE(N2)-VALUE(N2-1))*0.5                                    1749.000
      CC=OFF(N2)*OFF(N2)                                                1750.000
      SBB=1.                                                            1751.000
      IF(BB.LT.0.)SBB=-1.                                               1752.000
      ROOT=VALUE(N2)+CC/(BB+SBB*SQRT(BB*BB+CC))                         1753.000
      N2M1=N2-1                                                         1754.000
      A22=VALUE(N1)                                                     1755.000
      A12=A22-ROOT                                                      1756.000
      A23=OFF(N1+1)                                                     1757.000
      A13=A23                                                           1758.000
      DO 60 I=N1,N2M1                                                   1759.000
      A33=VALUE(I+1)                                                    1760.000
      IF(I.NE.N2M1)A34=OFF(I+2)                                         1761.000
      S=SQRT(A12*A12+A13*A13)                                           1762.000
      SI=A13/S                                                          1763.000
      CO=A12/S                                                          1764.000
      IF(I.NE.N1)OFF(I)=S                                               1765.000
      A11=CO*A22+SI*A23                                                 1766.000
      A12=CO*A23+SI*A33                                                 1767.000
      A13=SI*A34                                                        1768.000
      A21=CO*A23-SI*A22                                                 1769.000
      A22=CO*A33-SI*A23                                                 1770.000
      A23=CO*A34                                                        1771.000
      VALUE(I)=A11*CO+A12*SI                                            1772.000
      A12=-A11*SI+A12*CO                                                1773.000
      OFF(I+1)=A12                                                      1774.000
   60 A22=A22*CO-A21*SI                                                 1775.000
      VALUE(N2)=A22                                                     1776.000
      GO TO 20                                                          1777.000
      END                                                               1778.000
      SUBROUTINE START(OR,IS,LP)                                        1779.000
C     MODEL FILE AND FREE FORMAT READING STORAGE                        1780.000
      CHARACTER COMPID(72),LINE(148)                                    1781.000
      CHARACTER*8 LABEL(200),FTYPE(13)                                  1782.000
      INTEGER CN(200),CONTNT(13),FCH(20),FWD(20)                        1783.000
      COMMON /MODELC/COMPID,LABEL,FTYPE                                 1784.000
      COMMON /MODEL/ NATOM,NR,NT,NTYPE,NFTYPE(200),P(3,200),            1785.000
     1CN,SOF(200),UIJ(6,200),EP(3,200),ESOF(200),EUIJ(6,200),           1786.000
     2WAVE,CELL(6),ECELL(6),R(24,3,4),T(3,4),ICENT,                     1787.000
     3SFAC(14,13),CONTNT,SCALE(4),ETA,NQ,EXTNCT                         1788.000
C                                                                       1789.000
      DIMENSION RCP(6),OR(3,3),H(3,3),RR(3,3),RRI(3,3)                  1790.000
      CHARACTER*8 TAG                                                   1791.000
      DATA LINE/148*' '/                                                1792.000
 1000 FORMAT(/8X,' ****** EMPIRICAL ABSORPTION CORRECTION ***** ',      1793.000
     1//10X,' PLEASE WAIT. I AM READING THE MODEL FILE'//)              1794.000
 1001 FORMAT(72A1)                                                      1795.000
 1002 FORMAT(10X,' TYPE ORIENTATION MATRIX. ONE ROW PER LINE',/)        1796.000
 1003 FORMAT(10X,' RETYPE LINE. ONE ROW CONSISTS OF THREE NUMBERS',/)   1797.000
 1004 FORMAT(10X,' ORIENTATION MATRIX INCONSISTENT WITH CELL IN MODEL FI1798.000
     1LE',/)                                                            1799.000
 1005 FORMAT(//10X,' WAVELENGTH = ',F8.5,' ANGSTROMS ',/                1800.000
     110X,' TYPE A NEW VALUE OF THE WAVELENGTH OR <CR> TO START',//)    1801.000
 1006 FORMAT(/10X,' HAVE THE INDICES BEEN TRANSFORMED? ',/)             1802.000
 1007 FORMAT(10X,' THEN TYPE ROWS OF T WHERE HNEW = T.HOLD ',/)         1803.000
 1008 FORMAT(/10X,' INCLUDE THETA DEPENDENCE IN CORRECTION?',/)         1804.000
      IS=-1                                                             1805.000
      LIN=1                                                             1806.000
      LC=5                                                              1807.000
      WRITE(LP,1000)                                                    1808.000
      CALL MDLIN(LIN)                                                   1809.000
      CALL AB(CELL,RCP,OR,H,U,V)                                        1810.000
C READ THE ORIENTATION MATRIX: COLUMNS ARE A*,B*,C*                     1811.000
    5 WRITE(LP,1002)                                                    1812.000
      DO 8 I=1,3                                                        1813.000
    6 READ(LC,1001)(LINE(J),J=2,73)                                     1814.000
      CALL PARSE(LINE,NFS,FCH,FWD)                                      1815.000
      IF(NFS.EQ.3)GOTO 7                                                1816.000
      WRITE(LP,1003)                                                    1817.000
      GOTO 6                                                            1818.000
    7 DO 8 J=1,3                                                        1819.000
      CALL FFORMT(LINE,FCH(J),FWD(J),OR(I,J))                           1820.000
    8 CONTINUE                                                          1821.000
C CHECK FOR INDEX TRANSFORMATIONS                                       1822.000
    9 WRITE(LP,1006)                                                    1823.000
      READ(LC,1001)(LINE(J),J=2,73)                                     1824.000
      CALL PARSE(LINE,NFS,FCH,FWD)                                      1825.000
      IF(NFS.EQ.0)GOTO 9                                                1826.000
      CALL AFORMT(LINE,FCH(1),FWD(1),TAG)                               1827.000
      IF(TAG.EQ.'NO'.OR.TAG.EQ.'N')GOTO 15                              1828.000
      IF(TAG.EQ.'YES'.OR.TAG.EQ.'Y')GOTO 10                             1829.000
      GOTO 9                                                            1830.000
C READ THE TRANSFORMATION MATRIX                                        1831.000
   10 WRITE(LP,1007)                                                    1832.000
      DO 13 I=1,3                                                       1833.000
   11 READ(LC,1001)(LINE(J),J=2,73)                                     1834.000
      CALL PARSE(LINE,NFS,FCH,FWD)                                      1835.000
      IF(NFS.EQ.3)GOTO 12                                               1836.000
      WRITE(LP,1003)                                                    1837.000
      GOTO 11                                                           1838.000
   12 DO 13 J=1,3                                                       1839.000
      CALL FFORMT(LINE,FCH(J),FWD(J),RR(I,J))                           1840.000
   13 CONTINUE                                                          1841.000
C NOW TRANSFORM OR SO THAT IT REFERS TO THE NEW INDEXING                1842.000
      CALL AMI3(RR,RRI,DET)                                             1843.000
      CALL MAM(OR,RRI,RR)                                               1844.000
      DO 14 I=1,3                                                       1845.000
      DO 14 J=1,3                                                       1846.000
   14 OR(I,J)=RR(I,J)                                                   1847.000
C IN ANY EVENT CHECK THAT OR GIVES SAME UNIT CELL AS IN MODEL           1848.000
   15 DO 17 I=1,3                                                       1849.000
      DO 17 J=I,3                                                       1850.000
      Q=(OR(1,I)*OR(1,J)+OR(2,I)*OR(2,J)+OR(3,I)*OR(3,J))/(RCP(I)*      1851.000
     1   RCP(J))                                                        1852.000
      IF(I.EQ.J)Q=ABS(1.-Q)                                             1853.000
      IF(I.NE.J)Q=ABS(Q - RCP(9-I-J))*10.                               1854.000
      IF(Q.LT.01)GOTO 17                                                1855.000
      WRITE(LP,1004)                                                    1856.000
      GOTO 5                                                            1857.000
   17 CONTINUE                                                          1858.000
   18 WRITE(LP,1008)                                                    1859.000
      READ(LC,1001)(LINE(J),J=2,73)                                     1860.000
      CALL PARSE(LINE,NFS,FCH,FWD)                                      1861.000
      IF(NFS.EQ.0)GOTO 18                                               1862.000
      CALL AFORMT(LINE,FCH(1),FWD(1),TAG)                               1863.000
      IF(TAG.EQ.'YES'.OR.TAG.EQ.'Y')IS=1                                1864.000
      IF(TAG.EQ.'NO'.OR.TAG.EQ.'N')IS=0                                 1865.000
      IF(IS.LT.0)GOTO 18                                                1866.000
C                            VERIFY WAVELENGTH AND START                1867.000
      IF(WAVE.LE.01)WAVE=0.71069                                        1868.000
   19 WRITE(LP,1005)WAVE                                                1869.000
      READ(LC,1001)(LINE(J),J=2,73)                                     1870.000
      CALL PARSE(LINE,NFS,FCH,FWD)                                      1871.000
      IF(NFS.EQ.0)GOTO 20                                               1872.000
      CALL FFORMT(LINE,FCH(1),FWD(1),WAVE)                              1873.000
      GOTO 19                                                           1874.000
   20 RETURN                                                            1875.000
      END                                                               1876.000
C-----------------------------------------------------------------------1877.000
C     A - REAL CELL GENERATES B - RECIPROCAL CELL, G AND H THE          1878.000
C     REAL AND RECIPROCAL METRIC TENSORS, AND U AND V REAL AND          1879.000
C     RECIPROCAL CELL VOLUMES.                                          1880.000
      SUBROUTINE AB(A,B,G,H,U,V)                                        1881.000
      DIMENSION A(6),B(6),G(3,3),H(3,3)                                 1882.000
      CALL AG(A,G)                                                      1883.000
      CALL AMI3(G,H,U)                                                  1884.000
      CALL GA(H,B)                                                      1885.000
      U=SQRT(U)                                                         1886.000
      V=1.0/U                                                           1887.000
      RETURN                                                            1888.000
      END                                                               1889.000
C-----------------------------------------------------------------------1890.000
C     METRIC FROM CELL                                                  1891.000
      SUBROUTINE AG(A,G)                                                1892.000
      DIMENSION A(6),G(3,3)                                             1893.000
      DO 1 I=1,3                                                        1894.000
      DO 1 J=1,3                                                        1895.000
      G(I,J)=A(I)*A(J)                                                  1896.000
      IF(I.EQ.J)GOTO 1                                                  1897.000
      G(I,J)=G(I,J)*A(9-I-J)                                            1898.000
    1 CONTINUE                                                          1899.000
      RETURN                                                            1900.000
      END                                                               1901.000
C-----------------------------------------------------------------------1902.000
C     CELL FROM METRIC                                                  1903.000
      SUBROUTINE MAM (A,B,C)                                            1904.000
C     **********************                                            1905.000
C.... MULTIPLIES MATRIX-A-BY MATRIX-B-                                  1906.000
C....   A(3,3)*B(3,3)=C(3,3)                                            1907.000
      DIMENSION A(3,3), B(3,3), C(3,3)                                  1908.000
      DO 12 I=1,3                                                       1909.000
        DO 12 K=1,3                                                     1910.000
        C(I,K)=0.0                                                      1911.000
        DO 12 J=1,3                                                     1912.000
   12   C(I,K)=C(I,K)+A(I,J)*B(J,K)                                     1913.000
      RETURN                                                            1914.000
      END                                                               1915.000
      SUBROUTINE GA(G,A)                                                1916.000
      DIMENSION A(6),G(3,3)                                             1917.000
      DO 1 I=1,3                                                        1918.000
    1 A(I)=SQRT (G(I,I))                                                1919.000
      A(4)=G(2,3)/(A(2)*A(3))                                           1920.000
      A(5)=G(1,3)/(A(1)*A(3))                                           1921.000
      A(6)=G(1,2)/(A(1)*A(2))                                           1922.000
      RETURN                                                            1923.000
      END                                                               1924.000
      SUBROUTINE KRAKJA(JAA,NSET,IH,IK,IL)                              1925.000
      JA=JAA                                                            1926.000
      NSET=JA/16777216                                                  1927.000
      JA=JA-16777216*NSET                                               1928.000
      IH=JA/65536                                                       1929.000
      JA=JA-65536*IH                                                    1930.000
      IH=IH-128                                                         1931.000
      IK=JA/256                                                         1932.000
      IL=JA-256*IK-128                                                  1933.000
      IK=IK-128                                                         1934.000
      RETURN                                                            1935.000
      END                                                               1936.000
C                                                                       1937.000
C                                                                       1938.000
C***** END DIFABS ******************************************************1939.000
