C        COPYRIGHT 1985 UNIVERSITY OF GLASGOW                           0001.000
C        GX PACKAGE APRIL 1985                                          0002.000
      PROGRAM NABKWM                                                    0003.000
C.....GAUSSIAN ABSORPTION CORRECTION - K W MUIR FEB 1982                0004.000
C.....BASED ON COPPENS, CRYSTALLOGRAPHIC COMPUTING, ED. AHMED           0005.000
C.....INCORPORATES IDEAS FROM NORTHWESTERN AGNOSTC                      0006.000
C.....AUTOMATIC GRID SELECTION FROM SHELX                               0007.000
C.....PACKED DATA FILE HANDLING ADDED AND DRAWING REMOVED 20FEB84       0008.000
C     MODEL FILE STORAGE                                                0009.000
      CHARACTER COMPID(72)                                              0010.000
      CHARACTER*8 LABEL(200),FTYPE(13)                                  0011.000
      INTEGER CN(200),CONTNT(13)                                        0012.000
      COMMON /MODELC/COMPID,LABEL,FTYPE                                 0013.000
      COMMON /MODEL/ NATOM,NR,NT,NTYPE,NFTYPE(200),P(3,200),            0014.000
     1CN,SOF(200),UIJ(6,200),EP(3,200),ESOF(200),EUIJ(6,200),           0015.000
     2WAVEL,CELL(6),ECELL(6),RR(24,3,4),T(3,4),ICENT,                   0016.000
     3SFAC(14,13),CONTNT,SCALE(4),ETA,NQ,EXTNCT                         0017.000
C                                                                       0018.000
      COMMON EA(3000),EB(1000),EC(128),ED(128),EE(2000,4)               0019.000
      COMMON/ACELL/A(6),B(6),G(3,3),H(3,3),U,V,ANG(3)                   0020.000
      COMMON/AFACE/F(50,4),NIN(50),NOUT(50),COSI(50),COSD(50),NF        0021.000
      COMMON/ADATA/AMU,AFAC,NG(3),NSW,NPR,NA(16),IDRAW                  0022.000
      COMMON/AVERT/VR(100,3),IVR(100,10),NVR,VXTL,JVR(100),NIJ(100,2),  0023.000
     1EDGE(100),NE                                                      0024.000
      COMMON/ADEVS/LIN,LA,LB,LC,LP,LT,LG                                0025.000
      DIMENSION EF(64),EG(64),R(3,3),RI(3,3),S(2,2),SI(2,2)             0026.000
      DIMENSION X(3),Y(3),XX(3),YY(3),IIN(3)                            0027.000
      DATA EH/1.0E5/,EI/1.0E-5/,IRET/0/,NFACES/0/                       0028.000
      DATA EF/                                                          0029.000
     A.2113248660, .0694318440, .3300094780, .0337652430, .1693953070,  0030.000
     B.3806904070, .0198550720, .1016667610, .2372337950, .4082826790,  0031.000
     C.0130467360, .0674683170, .1602952160, .2833023030, .4255628310,  0032.000
     D.0092196830, .0479413720, .1150486630, .2063410230, .3160842510,  0033.000
     E.4373332960, .0068580960, .0357825580, .0863993420, .1563535480,  0034.000
     F.2423756820, .3404438160, .4459725260, .0052995330, .0277124880,  0035.000
     G.0671843990, .1222977960, .1910618780, .2709916110, .3591982250,  0036.000
     H.4524937450, .0024063900, .0126357220, .0308627240, .0567922360,  0037.000
     I.0899990070, .1299379040, .1759531740, .2272892640, .2831032460,  0038.000
     J.3424786600, .4044405660, .4679715530, .0013680690, .0071942440,  0039.000
     K.0176188720, .0325469620, .0518394220, .0753161930, .1027581020,  0040.000
     L.1339089400, .1684778670, .2061421210, .2465500450, .2893243620,  0041.000
     M.3340656990, .3803563190, .4277640190, .4758461670/               0042.000
      DATA EG/                                                          0043.000
     A.5000000000, .1739274230, .3260725780, .0856622460, .1803807870,  0044.000
     B.2339569670, .0506142680, .1111905170, .1568533230, .1813418920,  0045.000
     C.0333356720, .0747256750, .1095431810, .1346333600, .1477621120,  0046.000
     D.0235876680, .0534696630, .0800391640, .1015837130, .1167462680,  0047.000
     E.1245735230, .0175597300, .0400790440, .0607592850, .0786015840,  0048.000
     F.0927691990, .1025992320, .1076319270, .0135762300, .0311267620,  0049.000
     G.0475792560, .0623144860, .0747979940, .0845782600, .0913017080,  0050.000
     H.0947253050, .0061706150, .0142656940, .0221387190, .0296492920,  0051.000
     I.0366732410, .0430950810, .0488093260, .0537221350, .0577528340,  0052.000
     J.0608352360, .0629187280, .0639690980, .0035093050, .0081371970,  0053.000
     K.0126960330, .0171369320, .0214179490, .0254990300, .0293420470,  0054.000
     L.0329111110, .0361728970, .0390969480, .0416559620, .0438260470,  0055.000
     M.0455869390, .0469222000, .0478193600, .0482700440/               0056.000
 1000 FORMAT(///,3X,'H',3X,'K',3X,'L',6X,'FOBS',4X,'SIG',5X,'A',5X,     0057.000
     1'TBAR',4X,'IERR',5X,' THETA VALUES ',/)                           0058.000
 1001 FORMAT(3I4,F10.2,F7.2,2F8.5,I5,3F10.3)                            0059.000
 1002 FORMAT(//,' NO OF REFLECTIONS CORRECTED = ',I6,                   0060.000
     1//,' RANGE OF TRANSMISSION FACTORS ON FSQUARED = ',2F8.5,         0061.000
     2//,' SCALE FACTOR FOR FOBS = ',F8.5)                              0062.000
 1003 FORMAT(//,' ERROR STOP - FACES+VERTICES-EDGES.NE.2 OR',/14X,      0063.000
     1'FACES.GT.100  OR A VERTEX IS COMMON TO MORE THAN 10 FACES',//)   0064.000
 1004 FORMAT(//,' ERROR STOP - EXTEND TABLE OF EXP(-X)',//)             0065.000
 1005 FORMAT(/,' VERTEX COORDINATES IN CM ALONG A, B, AND C',           0066.000
     16X,' INTERSECTING PLANES',/)                                      0067.000
 1006 FORMAT(I5,2X,3F12.5,5X,10(X,I2))                                  0068.000
 1007 FORMAT(/,I4,' EDGES DEFINED BY VERTEX PAIRS. LENGTHS ARE IN CMS', 0069.000
     114(/7(X,I2,'-',I2,F5.3)))                                         0070.000
 1008 FORMAT(//,' ABSORPTION CORRECTION COMPLETED'/)                    0071.000
C                                                                       0072.000
C.....READ MODEL FILE, USERS REQUIREMENTS, AND SET UP TABLES            0073.000
  999 CONTINUE                                                          0074.000
      RTOD=45.0/ATAN(1.0)                                               0075.000
      CALL GABIN(IRET,NFACES)                                           0076.000
      DO 1 I=1,3000                                                     0077.000
    1 EA(I)=EXP(0.1*REAL(1-I))                                          0078.000
      DO 2 I=1,1000                                                     0079.000
    2 EB(I)=EXP(0.0001*REAL(1-I))                                       0080.000
      M=0                                                               0081.000
      DO 5 I=1,16                                                       0082.000
      IF(I.LE.8)GOTO 3                                                  0083.000
      IF(I.EQ.12)GOTO 3                                                 0084.000
      IF(I.NE.16)GOTO 5                                                 0085.000
    3 N=2*M                                                             0086.000
      NA(I)=N                                                           0087.000
      DO 4 J=1,I                                                        0088.000
      K=J+N                                                             0089.000
      L=2*I+N+1-J                                                       0090.000
      EC(K)=EF(J+M)                                                     0091.000
      ED(K)=EG(J+M)                                                     0092.000
      EC(L)=1.0 - EC(K)                                                 0093.000
    4 ED(L)=ED(K)                                                       0094.000
      M=M+I                                                             0095.000
    5 CONTINUE                                                          0096.000
C.....UNIT VECTORS ALONG FACE NORMALS (RECIPROCAL COMPONENTS)           0097.000
      DO 7 I=1,NF                                                       0098.000
      Q=0.                                                              0099.000
      DO 6 J=1,3                                                        0100.000
      DO 6 K=1,3                                                        0101.000
    6 Q=Q+H(J,K)*F(I,J)*F(I,K)                                          0102.000
      Q=1.0/SQRT(Q)                                                     0103.000
      DO 7 J=1,3                                                        0104.000
    7 F(I,J)=F(I,J)*Q                                                   0105.000
C.....CRYSTAL VERTICES (REAL COMPONENTS) AND LIMITS ALONG A             0106.000
      WRITE(LT,1005)                                                    0107.000
      NVR=0                                                             0108.000
      NE=0                                                              0109.000
      XMAX=-EH                                                          0110.000
      XMIN=-XMAX                                                        0111.000
      NF1=NF-1                                                          0112.000
      NF2=NF-2                                                          0113.000
C.....PLANE I                                                           0114.000
      DO 15 I=1,NF2                                                     0115.000
      DO 8 L=1,3                                                        0116.000
    8 R(1,L)=F(I,L)                                                     0117.000
      X(1)=F(I,4)                                                       0118.000
      J1=I+1                                                            0119.000
C.....PLANE J                                                           0120.000
      DO 14 J=J1,NF1                                                    0121.000
      DO 9 L=1,3                                                        0122.000
    9 R(2,L)=F(J,L)                                                     0123.000
      X(2)=F(J,4)                                                       0124.000
C.....PLANE K                                                           0125.000
      J2=J+1                                                            0126.000
      DO 13 K=J2,NF                                                     0127.000
      DO 10 L=1,3                                                       0128.000
   10 R(3,L)=F(K,L)                                                     0129.000
      CALL AMI3(R,RI,DET)                                               0130.000
      IF(ABS(DET).LT.EI)GOTO 13                                         0131.000
      X(3)=F(K,4)                                                       0132.000
      CALL MV(RI,X,Y)                                                   0133.000
C.....REJECT INTERSECTION IF OUTSIDE CRYSTAL                            0134.000
      M=NVR+1                                                           0135.000
      IF(M.GT.100)GOTO 36                                               0136.000
      N=0                                                               0137.000
      DO 11 L=1,NF                                                      0138.000
      Q=F(L,4)+EI-F(L,1)*Y(1)-F(L,2)*Y(2)-F(L,3)*Y(3)                   0139.000
      IF(Q.LT.0.)GOTO 13                                                0140.000
      IF(ABS(Q-EI).GT.EI)GOTO 11                                        0141.000
      N=N+1                                                             0142.000
      IF(N.GT.10)GOTO 36                                                0143.000
      IVR(M,N)=L                                                        0144.000
   11 CONTINUE                                                          0145.000
      JVR(M)=N                                                          0146.000
      DO 12 L=1,3                                                       0147.000
      ANG(L)=A(L)*Y(L)                                                  0148.000
   12 VR(M,L)=Y(L)                                                      0149.000
      CALL EDGES(M)                                                     0150.000
      IF(M.LT.0)GOTO 13                                                 0151.000
      NVR=NVR+1                                                         0152.000
      WRITE(LT,1006)NVR,ANG,(IVR(NVR,L),L=1,N)                          0153.000
      XMAX=AMAX1(XMAX,Y(1))                                             0154.000
      XMIN=AMIN1(XMIN,Y(1))                                             0155.000
   13 CONTINUE                                                          0156.000
   14 CONTINUE                                                          0157.000
   15 CONTINUE                                                          0158.000
      WRITE  (LT,1007)NE,((NIJ(I,J),J=1,2),EDGE(I),I=1,NE)              0159.000
      IF(NF+NVR-NE.NE.2)GOTO 36                                         0160.000
C.....SET UP THE REQUIRED GRID                                          0161.000
      VXTL=0.                                                           0162.000
      NP=0                                                              0163.000
      CALL LIMITS(XMAX,XMIN,MX,NX,1)                                    0164.000
C.....LOOP 1 - CHOOSE X GRID COORDINATES AND WEIGHTS                    0165.000
      DO 24 I=1,MX                                                      0166.000
      X(1)=XMIN+(XMAX-XMIN)*EC(NX+I)                                    0167.000
      RX=(XMAX-XMIN)*ED(NX+I)                                           0168.000
      YMAX=-EH                                                          0169.000
      YMIN=-YMAX                                                        0170.000
C.....FIND LIMITS ON Y FOR GIVEN X                                      0171.000
      DO 19 J=1,NF1                                                     0172.000
      DO 16 L=1,2                                                       0173.000
   16 S(1,L)=F(J,L+1)                                                   0174.000
      J1=J+1                                                            0175.000
      DO 19 K=J1,NF                                                     0176.000
      DO 17 L=1,2                                                       0177.000
   17 S(2,L)=F(K,L+1)                                                   0178.000
      CALL AMI2(S,SI,DET)                                               0179.000
      IF(ABS(DET).LT.EI)GOTO 19                                         0180.000
      Y(1)=F(J,4)-F(J,1)*X(1)                                           0181.000
      Y(2)=F(K,4)-F(K,1)*X(1)                                           0182.000
      X(2)=SI(1,1)*Y(1)+SI(1,2)*Y(2)                                    0183.000
      X(3)=SI(2,1)*Y(1)+SI(2,2)*Y(2)                                    0184.000
      DO 18 L=1,NF                                                      0185.000
      Q=F(L,4)+EI-F(L,1)*X(1)-F(L,2)*X(2)-F(L,3)*X(3)                   0186.000
      IF(Q.LT.0.)GOTO 19                                                0187.000
   18 CONTINUE                                                          0188.000
      YMAX=AMAX1(YMAX,X(2))                                             0189.000
      YMIN=AMIN1(YMIN,X(2))                                             0190.000
   19 CONTINUE                                                          0191.000
      CALL LIMITS(YMAX,YMIN,MY,NY,2)                                    0192.000
C.....LOOP 2 - CHOOSE Y GRID COORDINATES AND WEIGHTS                    0193.000
      DO 23 J=1,MY                                                      0194.000
      X(2)=YMIN+(YMAX-YMIN)*EC(NY+J)                                    0195.000
      RY=(YMAX-YMIN)*ED(NY+J)                                           0196.000
      ZMAX=-EH                                                          0197.000
      ZMIN=-ZMAX                                                        0198.000
C.....FIND LIMITS OF Z FOR GIVEN X AND Y                                0199.000
      DO 21 K=1,NF                                                      0200.000
      IF(ABS(F(K,3)).LT.EI)GOTO 21                                      0201.000
      X(3)=(F(K,4)-F(K,1)*X(1)-F(K,2)*X(2))/F(K,3)                      0202.000
      DO 20 L=1,NF                                                      0203.000
      Q=F(L,4)+EI-F(L,1)*X(1)-F(L,2)*X(2)-F(L,3)*X(3)                   0204.000
      IF(Q.LT.0.)GOTO 21                                                0205.000
   20 CONTINUE                                                          0206.000
      ZMAX=AMAX1(ZMAX,X(3))                                             0207.000
      ZMIN=AMIN1(ZMIN,X(3))                                             0208.000
   21 CALL LIMITS(ZMAX,ZMIN,MZ,NZ,3)                                    0209.000
C.....GET Z FOR GIVEN X,Y. STORE GRID FRACTIONAL COORDINATES AND WEIGHT 0210.000
      DO 22 K=1,MZ                                                      0211.000
      NP=NP+1                                                           0212.000
      IF(NP.GT.2000)GOTO 400                                            0213.000
      EE(NP,1)=X(1)                                                     0214.000
      EE(NP,2)=X(2)                                                     0215.000
      EE(NP,3)=ZMIN+(ZMAX-ZMIN)*EC(NZ+K)                                0216.000
      EE(NP,4)=RX*RY*(ZMAX-ZMIN)*ED(NZ+K)                               0217.000
   22 VXTL=VXTL+EE(NP,4)                                                0218.000
   23 CONTINUE                                                          0219.000
   24 CONTINUE                                                          0220.000
      VOL=VXTL*U*1000.                                                  0221.000
      WRITE(LT,25)NP,VOL                                                0222.000
   25 FORMAT(/I6,' GRID POINTS GIVE A CRYSTAL VOLUME OF ',F12.6,' MM**3 0223.000
     1',/)                                                              0224.000
C.....LOOP THROUGH REFLECTIONS                                          0225.000
      IF(NPR.GT.0)WRITE(LT,1000)                                        0226.000
      AMAX=0.                                                           0227.000
      AMIN=EH                                                           0228.000
      AMUTEN=AMU*10.                                                    0229.000
      ABAR=0.                                                           0230.000
      NREF=0                                                            0231.000
      NERR=0                                                            0232.000
      REWIND(LA)                                                        0233.000
      REWIND(LB)                                                        0234.000
   26 READ(LA,END=35)JCC,JBB,FF,SIGI,(X(I),Y(I),I=1,3)                  0235.000
      JAA=JCC                                                           0236.000
C.....X,Y - DIRECTION COSINES OF REVERSED INCIDENT AND DIFFRACTED       0237.000
C.....RAYS ON RECIPROCAL AXES. MULTIPLY BY RECIPROCAL AXES TO GET       0238.000
C.....UNIT VECTOR COMPONENTS ON REAL AXES XX AND YY.                    0239.000
      ISET=JAA/16777216                                                 0240.000
      JAA=JAA-ISET*16777216                                             0241.000
      IIN(1)=JAA/65536                                                  0242.000
      JAA=JAA-IIN(1)*65536                                              0243.000
      IIN(1)=IIN(1)-128                                                 0244.000
      IIN(2)=JAA/256                                                    0245.000
      IIN(3)=JAA-256*IIN(2)-128                                         0246.000
      IIN(2)=IIN(2)-128                                                 0247.000
      MX=JBB/65536                                                      0248.000
      MNUM=JBB-MX*65536                                                 0249.000
      DO 27 I=1,3                                                       0250.000
      XX(I)=B(I)*X(I)                                                   0251.000
   27 YY(I)=B(I)*Y(I)                                                   0252.000
C.....NOW ELIMINATE FACE IF ANGLE BETWEEN BEAM AND NORMAL GE 90 DEGREES 0253.000
      MX=0                                                              0254.000
      MY=0                                                              0255.000
      DO 30 I=1,NF                                                      0256.000
      CI=0.                                                             0257.000
      CD=0.                                                             0258.000
      DO 28 J=1,3                                                       0259.000
      CI=CI+XX(J)*F(I,J)                                                0260.000
   28 CD=CD+YY(J)*F(I,J)                                                0261.000
      IF(CI.LE.0.)GOTO 29                                               0262.000
      MX=MX+1                                                           0263.000
      NIN(MX)=I                                                         0264.000
      COSI(MX)=CI                                                       0265.000
   29 IF(CD.LE.0.)GOTO 30                                               0266.000
      MY=MY+1                                                           0267.000
      NOUT(MY)=I                                                        0268.000
      COSD(MY)=CD                                                       0269.000
   30 CONTINUE                                                          0270.000
      TAB=0.                                                            0271.000
      ASTAR=0.                                                          0272.000
C.....LOOP OVER POINTS FOR THIS PLANE                                   0273.000
      DO 33 I=1,NP                                                      0274.000
      GK=EH                                                             0275.000
      DO 31 J=1,MX                                                      0276.000
      NX=NIN(J)                                                         0277.000
      Q=(F(NX,4)-F(NX,1)*EE(I,1)-F(NX,2)*EE(I,2)-F(NX,3)*EE(I,3))/      0278.000
     1COSI(J)                                                           0279.000
   31 GK=AMIN1(GK,Q)                                                    0280.000
      HK=EH                                                             0281.000
      DO 32 J=1,MY                                                      0282.000
      NY=NOUT(J)                                                        0283.000
      Q=(F(NY,4)-F(NY,1)*EE(I,1)-F(NY,2)*EE(I,2)-F(NY,3)*EE(I,3))/      0284.000
     1COSD(J)                                                           0285.000
   32 HK=AMIN1(HK,Q)                                                    0286.000
      EX=(GK+HK)*AMUTEN+0.0005                                          0287.000
      M=EX                                                              0288.000
      IF(M.GT.2999)GOTO 37                                              0289.000
      DX=M                                                              0290.000
      N=(EX-DX)*1000.                                                   0291.000
      AB=EA(M+1)*EB(N+1)*EE(I,4)                                        0292.000
      TAB=TAB+AB*(GK+HK)                                                0293.000
   33 ASTAR=ASTAR+AB                                                    0294.000
      TBAR=TAB/ASTAR                                                    0295.000
      ASTAR=ASTAR/VXTL                                                  0296.000
      FF=FF/ASTAR                                                       0297.000
      SIGI=SIGI/ASTAR                                                   0298.000
      JBB=TBAR*10000.+.5                                                0299.000
      JBB=JBB*65536+MNUM                                                0300.000
      WRITE(LB)JCC,JBB,FF,SIGI                                          0301.000
      NREF=NREF+1                                                       0302.000
      AMAX=AMAX1(AMAX,ASTAR)                                            0303.000
      AMIN=AMIN1(AMIN,ASTAR)                                            0304.000
      ABAR=ABAR+ASTAR                                                   0305.000
C.....TEST AND PRINT. CHECK THETA THREE WAYS                            0306.000
      Q=0.                                                              0307.000
      TH1=0.                                                            0308.000
      TH2=0.                                                            0309.000
      TH3=0.                                                            0310.000
      DO 34 I=1,3                                                       0311.000
      X(I)=G(I,1)*XX(1)+G(I,2)*XX(2)+G(I,3)*XX(3)                       0312.000
      Y(I)=G(I,1)*YY(1)+G(I,2)*YY(2)+G(I,3)*YY(3)                       0313.000
C.....X AND Y ARE COMPONENTS ON RECIPROCAL AXES OF XX AND YY            0314.000
      TH1=TH1+XX(I)*Y(I)                                                0315.000
      TH2=TH2+XX(I)*REAL(IIN(I))                                        0316.000
      TH3=TH3+YY(I)*REAL(IIN(I))                                        0317.000
      DO 34 J=1,3                                                       0318.000
   34 Q=Q+H(I,J)*REAL(IIN(I)*IIN(J))                                    0319.000
      Q=1.0/SQRT(Q)                                                     0320.000
      TH1=.5*RTOD*ACOS(-TH1)                                            0321.000
      TH2=RTOD*ASIN(TH2*Q)                                              0322.000
      TH3=RTOD*ASIN(TH3*Q)                                              0323.000
      IERR=0                                                            0324.000
      IF(ABS(TH1-TH2).GT..2 .OR. ABS(TH1-TH3).GT..2 .OR. ABS(TH2-TH3).GT0325.000
     1..2)IERR=1                                                        0326.000
      IF(IERR.EQ.1)NERR=NERR+1                                          0327.000
      FF=SQRT(ABS(FF))                                                  0328.000
      SIGI=0.5*SIGI/FF                                                  0329.000
      IF(NREF.LE.NPR.OR.IERR.EQ.1)                                      0330.000
     1WRITE(LT,1001)IIN,FF,SIGI,ASTAR,TBAR,IERR,TH1,TH2,TH3             0331.000
      IF(NERR.GT.100)GOTO 35                                            0332.000
      IF(IRET.EQ.1.AND.NREF.EQ.NPR)GOTO 35                              0333.000
      GOTO 26                                                           0334.000
   35 ABAR=SQRT(ABAR/REAL(NREF))                                        0335.000
      WRITE(LT,1002)NREF,AMIN,AMAX,ABAR                                 0336.000
      IF(IRET.EQ.0)GOTO 38                                              0337.000
      IRET=0                                                            0338.000
      NFACES=NF                                                         0339.000
      GOTO 999                                                          0340.000
   36 WRITE(LT,1003)                                                    0341.000
      STOP                                                              0342.000
   37 WRITE(LT,1004)                                                    0343.000
      STOP                                                              0344.000
   38 IF(NERR.GT.100)GOTO 40                                            0345.000
C     DO 39 I=1,NQ                                                      0346.000
C  39 SCALE(NQ)=SCALE(NQ)/ABAR                                          0347.000
C     CALL MDLOUT(LIN)                                                  0348.000
      WRITE(LT,1008)                                                    0349.000
   40 STOP                                                              0350.000
  400 WRITE(LT,2000)                                                    0351.000
 2000 FORMAT(' 2000 GRID POINTS IS CURRENT MAXIMUM ')                   0352.000
      STOP                                                              0353.000
      END                                                               0354.000
C-----------------------------------------------------------------------0355.000
C     DETECT DUPLICATE VERTICES AND FIND EDGES                          0356.000
      SUBROUTINE EDGES(M)                                               0357.000
      COMMON/ACELL/A(6),B(6),G(3,3),H(3,3),U,V,ANG(3)                   0358.000
      COMMON/AVERT/VR(100,3),IVR(100,10),NVR,VXTL,JVR(100),NIJ(100,2),  0359.000
     1EDGE(100),NE                                                      0360.000
      COMMON/ADEVS/LIN,LA,LB,LC,LP,LT,LG                                0361.000
 1000 FORMAT(' ERROR STOP - TOO MANY EDGES ',//)                        0362.000
      ME=NE                                                             0363.000
      JM=JVR(M)                                                         0364.000
      DO 3 I=1,NVR                                                      0365.000
      JI=JVR(I)                                                         0366.000
      J=0                                                               0367.000
      DO 1 K=1,JM                                                       0368.000
      DO 1 L=1,JI                                                       0369.000
      IF(IVR(I,L).EQ.IVR(M,K))J=J+1                                     0370.000
    1 CONTINUE                                                          0371.000
      IF(J.EQ.JM.AND.JM.EQ.JI)GOTO 4                                    0372.000
      IF(J.NE.2)GOTO 3                                                  0373.000
      NE=NE+1                                                           0374.000
      IF(NE.GT.100)GOTO 5                                               0375.000
      NIJ(NE,1)=I                                                       0376.000
      NIJ(NE,2)=M                                                       0377.000
      EDGE(NE)=0.                                                       0378.000
      DO 2 K=1,3                                                        0379.000
      DO 2 L=1,3                                                        0380.000
    2 EDGE(NE)=EDGE(NE)+G(K,L)*(VR(I,K)-VR(M,K))*(VR(I,L)-VR(M,L))      0381.000
      EDGE(NE)=SQRT(EDGE(NE) )                                          0382.000
    3 CONTINUE                                                          0383.000
      GOTO 6                                                            0384.000
    4 M=-M                                                              0385.000
      NE=ME                                                             0386.000
      GOTO 6                                                            0387.000
    5 WRITE(LT,1000)                                                    0388.000
      STOP                                                              0389.000
    6 RETURN                                                            0390.000
      END                                                               0391.000
C-----------------------------------------------------------------------0392.000
C     CHOOSE GRID SIZE                                                  0393.000
      SUBROUTINE LIMITS(XMAX,XMIN,MX,NX,I)                              0394.000
      COMMON/ACELL/A(6),B(6),G(3,3),H(3,3),U,V                          0395.000
      COMMON/ADATA/AMU,AFAC,NG(3),NSW,NPR,NA(16)                        0396.000
      MX=MIN0(8,MAX0(1,INT(.5+ 2.5*AMU*A(I)*(XMAX-XMIN)/AFAC)))*2       0397.000
      IF(NSW.NE.0)MX=NG(I)                                              0398.000
      NX=NA(MX/2)                                                       0399.000
      RETURN                                                            0400.000
      END                                                               0401.000
C-----------------------------------------------------------------------0402.000
C     READ MODEL FILE AND USERS REQUIREMENTS                            0403.000
      SUBROUTINE GABIN(IRET,NFACES)                                     0404.000
C     MODEL FILE AND FREE FORMAT READING STORAGE                        0405.000
      CHARACTER COMPID(72),LINE(148)                                    0406.000
      CHARACTER*8 LABEL(200),FTYPE(13)                                  0407.000
      INTEGER CN(200),CONTNT(13),FCH(20),FWD(20)                        0408.000
      COMMON /MODELC/COMPID,LABEL,FTYPE                                 0409.000
      COMMON /MODEL/ NATOM,NR,NT,NTYPE,NFTYPE(200),P(3,200),            0410.000
     1CN,SOF(200),UIJ(6,200),EP(3,200),ESOF(200),EUIJ(6,200),           0411.000
     2WAVEL,CELL(6),ECELL(6),R(24,3,4),T(3,4),ICENT,                    0412.000
     3SFAC(14,13),CONTNT,SCALE(4),ETA,NQ,EXTNCT                         0413.000
C                                                                       0414.000
      COMMON/ACELL/A(6),B(6),G(3,3),H(3,3),U,V,ANG(3)                   0415.000
      COMMON/AFACE/F(50,4),NIN(50),NOUT(50),COSI(50),COSD(50),NF        0416.000
      COMMON/ADATA/AMU,AFAC,NG(3),NSW,NPR,NA(16),IDRAW                  0417.000
      COMMON/ADEVS/LIN,LA,LB,LC,LP,LT,LG                                0418.000
      DIMENSION C(20)                                                   0419.000
      CHARACTER*8 TAG                                                   0420.000
      DATA LINE/148*' '/                                                0421.000
 1000 FORMAT(//10X,'PLEASE WAIT. I AM READING THE MODEL FILE'//)        0422.000
 1001 FORMAT(///20X,'<<<MENU FOR GAUSSIAN ABSORPTION CORRECTION>>>'//   0423.000
     110X,' 1. FACE   H K L D(MMS)',/                                   0424.000
     210X,' 2. MU     LINEAR ABSORPTION COEFFICIENT(CM-1)',/            0425.000
     310X,' 3. GRID   <PERCENT ERROR> OR <NX NY NZ>',/                  0426.000
     410X,' 4. LIST   <N> LIST FIRST N REFLECTIONS',/                   0427.000
     524X,' NEGATIVE N FORCES REINITIALISATION AFTER N REFLECTIONS',/   0428.000
     610X,' 5. <CR>   START CALCULATION',//)                            0429.000
 1002 FORMAT( 10X,'TYPE AN INSTRUCTION')                                0430.000
 1003 FORMAT(72A1)                                                      0431.000
 1004 FORMAT( 10X,'INPUT ERROR - CHECK AND RETYPE LAST LINE')           0432.000
 1005 FORMAT(1H1'GAUSSIAN ABSORPTION CORRECTION FOR',/,1X,72A1,         0433.000
     1//,' LINEAR ABSORPTION COEFFICIENT = ',F6.2,' CM-1',              0434.000
     2/ ,' LIST FIRST ',I6,' REFLECTIONS',                              0435.000
     3//,' FACE  H   K   L    D(CMS)'/)                                 0436.000
 1006 FORMAT(X,I3,X,3F4.0,F9.5)                                         0437.000
 1007 FORMAT(/,' FLEXIBLE GRID BASED ON ',F5.2,' PERCENT ERROR'/)       0438.000
 1008 FORMAT(/,' FIXED GRID OF ',I2,' X ',I2,' X ',I2,' ALONG A, B, AND 0439.000
     1C'/)                                                              0440.000
 1009 FORMAT(/,10X,' CELL DIMENSIONS ',/,' A = ',F10.3,' ALPHA = ',     0441.000
     1F10.2,/,' B = ',F10.3,' BETA  = ',F10.2,/,' C = ',F10.3,          0442.000
     2' GAMMA = ',F10.2,/,' CELL VOLUME = ',F10.3,' CUBIC ANGSTROMS ',/)0443.000
      LIN=1                                                             0444.000
      LA=2                                                              0445.000
      LB=4                                                              0446.000
      LC=5                                                              0447.000
      LP=6                                                              0448.000
      LG=7                                                              0449.000
      LT=8                                                              0450.000
      WRITE(LP,1000)                                                    0451.000
      CALL MDLIN(LIN)                                                   0452.000
      CALL AB(CELL,B,G,H,U,V)                                           0453.000
      CALL GA(G,A)                                                      0454.000
      AMU=0.                                                            0455.000
      AFAC=1.0                                                          0456.000
      NSW=0                                                             0457.000
      NPR=20                                                            0458.000
      NF=0                                                              0459.000
      IDRAW=0                                                           0460.000
      DO 1 I=1,3                                                        0461.000
    1 NG(I)=0                                                           0462.000
      DO 2 I=1,NTYPE                                                    0463.000
    2 AMU=AMU+V*SFAC(12,I)*REAL(CONTNT(I))                              0464.000
    3 WRITE(LP,1001)                                                    0465.000
    4 WRITE(LP,1002)                                                    0466.000
    5 READ(LC,1003)(LINE(I),I=2,73)                                     0467.000
      CALL PARSE(LINE,NFS,FCH,FWD)                                      0468.000
      IF(NFS.EQ.0)GOTO 14                                               0469.000
      CALL AFORMT(LINE,FCH(1),FWD(1),TAG)                               0470.000
      IF(TAG.EQ.'MENU')GOTO 3                                           0471.000
      IF(NFS.EQ.1)GOTO 13                                               0472.000
      NFS=NFS-1                                                         0473.000
      DO 6 I=1,NFS                                                      0474.000
      CALL FFORMT(LINE,FCH(I+1),FWD(I+1),C(I))                          0475.000
    6 CONTINUE                                                          0476.000
      IF(TAG.NE.'FACE')GOTO 8                                           0477.000
      IF(NFS.NE.4.OR.C(4).LT.0..OR.NF.GE.50)GOTO 13                     0478.000
      NF=NF+1                                                           0479.000
      DO 7 I=1,4                                                        0480.000
    7 F(NF,I)=C(I)                                                      0481.000
      F(NF,4)=0.1*F(NF,4)                                               0482.000
      GOTO 5                                                            0483.000
    8 IF(TAG.NE.'MU')GOTO 9                                             0484.000
      AMU=C(1)                                                          0485.000
      GOTO 5                                                            0486.000
    9 IF(TAG.NE.'LIST')GOTO 10                                          0487.000
      IF(C(1).LT.-1.)IRET=1                                             0488.000
      NPR=NINT(ABS(C(1)))                                               0489.000
      GOTO 5                                                            0490.000
   10 IF(TAG.NE.'GRID')GOTO 13                                          0491.000
      IF(NFS.NE.1)GOTO 11                                               0492.000
      AFAC=AMAX1(C(1),0.01)                                             0493.000
      GOTO 5                                                            0494.000
   11 IF(NFS.NE.3)GOTO 13                                               0495.000
      NSW=1                                                             0496.000
      DO 12 I=1,3                                                       0497.000
      J=NINT(C(I))/2                                                    0498.000
      IF(J.LE.8)GOTO 12                                                 0499.000
      IF(J.EQ.12)GOTO 12                                                0500.000
      IF(J.NE.16)GOTO 13                                                0501.000
   12 NG(I)=2*J                                                         0502.000
      IF(NG(1)*NG(2)*NG(3).LE.10240)GOTO 5                              0503.000
   13 WRITE(LP,1004)                                                    0504.000
      GOTO 5                                                            0505.000
   14 WRITE(LT,1005)COMPID,AMU,NPR                                      0506.000
      IF(NF.EQ.0)NF=NFACES                                              0507.000
      DO 15 I=1,NF                                                      0508.000
   15 WRITE(LT,1006)I,(F(I,J),J=1,4)                                    0509.000
      IF(NSW.EQ.0)WRITE(LT,1007)AFAC                                    0510.000
      IF(NSW.NE.0)WRITE(LT,1008)NG                                      0511.000
      FAC=45.0/ATAN(1.0)                                                0512.000
      DO 16 I=1,3                                                       0513.000
   16 ANG(I)=FAC*ACOS(A(I+3))                                           0514.000
      WRITE(LT,1009)(A(I),ANG(I),I=1,3),U                               0515.000
   17 RETURN                                                            0516.000
      END                                                               0517.000
C-----------------------------------------------------------------------0518.000
C     A - REAL CELL GENERATES B - RECIPROCAL CELL, G AND H THE          0519.000
C     REAL AND RECIPROCAL METRIC TENSORS, AND U AND V REAL AND          0520.000
C     RECIPROCAL CELL VOLUMES.                                          0521.000
      SUBROUTINE AB(A,B,G,H,U,V)                                        0522.000
      DIMENSION A(6),B(6),G(3,3),H(3,3)                                 0523.000
      CALL AG(A,G)                                                      0524.000
      CALL AMI3(G,H,U)                                                  0525.000
      CALL GA(H,B)                                                      0526.000
      U=SQRT(U)                                                         0527.000
      V=1.0/U                                                           0528.000
      RETURN                                                            0529.000
      END                                                               0530.000
C-----------------------------------------------------------------------0531.000
C     METRIC FROM CELL                                                  0532.000
      SUBROUTINE AG(A,G)                                                0533.000
      DIMENSION A(6),G(3,3)                                             0534.000
      DO 1 I=1,3                                                        0535.000
      DO 1 J=1,3                                                        0536.000
      G(I,J)=A(I)*A(J)                                                  0537.000
      IF(I.EQ.J)GOTO 1                                                  0538.000
      G(I,J)=G(I,J)*A(9-I-J)                                            0539.000
    1 CONTINUE                                                          0540.000
      RETURN                                                            0541.000
      END                                                               0542.000
C-----------------------------------------------------------------------0543.000
C     CELL FROM METRIC                                                  0544.000
      SUBROUTINE GA(G,A)                                                0545.000
      DIMENSION A(6),G(3,3)                                             0546.000
      DO 1 I=1,3                                                        0547.000
    1 A(I)=SQRT (G(I,I))                                                0548.000
      A(4)=G(2,3)/(A(2)*A(3))                                           0549.000
      A(5)=G(1,3)/(A(1)*A(3))                                           0550.000
      A(6)=G(1,2)/(A(1)*A(2))                                           0551.000
      RETURN                                                            0552.000
      END                                                               0553.000
C-----------------------------------------------------------------------0554.000
C     INVERT UNSYMMETRICAL 2X2 MATRIX                                   0555.000
      SUBROUTINE AMI2(A,B,DET)                                          0556.000
      DIMENSION A(2,2),B(2,2)                                           0557.000
      DET=A(1,1)*A(2,2)-A(1,2)*A(2,1)                                   0558.000
      IF(DET)1,2,1                                                      0559.000
    1 B(1,1)=A(2,2)/DET                                                 0560.000
      B(1,2)=-A(1,2)/DET                                                0561.000
      B(2,1)=-A(2,1)/DET                                                0562.000
      B(2,2)=A(1,1)/DET                                                 0563.000
    2 RETURN                                                            0564.000
      END                                                               0565.000
C-----------------------------------------------------------------------0566.000
C     INVERT UNSYMMETRICAL 3X3 MATRIX                                   0567.000
      SUBROUTINE AMI3(A,B,DET)                                          0568.000
      DIMENSION A(3,3),B(3,3)                                           0569.000
      B(1,1)=A(2,2)*A(3,3)-A(3,2)*A(2,3)                                0570.000
      B(2,1)=A(2,1)*A(3,3)-A(3,1)*A(2,3)                                0571.000
      B(3,1)=A(2,1)*A(3,2)-A(3,1)*A(2,2)                                0572.000
      DET=A(1,1)*B(1,1)-A(1,2)*B(2,1)+A(1,3)*B(3,1)                     0573.000
      IF(DET)1,3,1                                                      0574.000
    1 B(1,2)=A(1,2)*A(3,3)-A(3,2)*A(1,3)                                0575.000
      B(2,2)=A(1,1)*A(3,3)-A(3,1)*A(1,3)                                0576.000
      B(3,2)=A(1,1)*A(3,2)-A(3,1)*A(1,2)                                0577.000
      B(1,3)=A(1,2)*A(2,3)-A(2,2)*A(1,3)                                0578.000
      B(2,3)=A(1,1)*A(2,3)-A(2,1)*A(1,3)                                0579.000
      B(3,3)=A(1,1)*A(2,2)-A(2,1)*A(1,2)                                0580.000
      DO 2 I=1,3                                                        0581.000
      DO 2 J=1,3                                                        0582.000
    2 B(I,J)=(-1.)**(I+J)*B(I,J)/DET                                    0583.000
    3 RETURN                                                            0584.000
      END                                                               0585.000
C-----------------------------------------------------------------------0586.000
C     MATRIX*VECTOR Z(3)=X(3,3)*Y(3)                                    0587.000
      SUBROUTINE MV(X,Y,Z)                                              0588.000
      DIMENSION X(3,3),Y(3),Z(3)                                        0589.000
      DO 1 I=1,3                                                        0590.000
      Z(I)=0.                                                           0591.000
      DO 1 J=1,3                                                        0592.000
    1 Z(I)=Z(I)+X(I,J)*Y(J)                                             0593.000
      RETURN                                                            0594.000
      END                                                               0595.000
C-----------------------------------------------------------------------0596.000
C     A(6) B(6) REAL AND RECIPROCAL CELL                                0597.000
C     C ORTHOGONALISES, D IS ITS INVERSE:  Y = CX AND X = DY            0598.000
      SUBROUTINE OM(A,B,C,D)                                            0599.000
      DIMENSION A(6),B(6),C(3,3),D(3,3)                                 0600.000
      DO 1 I=1,3                                                        0601.000
      DO 1 J=1,3                                                        0602.000
      C(I,J)=0.                                                         0603.000
    1 D(I,J)=0.                                                         0604.000
      SA= SQRT(ABS(1.0 - A(4)*A(4)))                                    0605.000
      SB= SQRT(ABS(1.0 - A(5)*A(5)))                                    0606.000
      SGS=SQRT(ABS(1.0 - B(6)*B(6)))                                    0607.000
      C(1,1)= A(1)*SB*SGS                                               0608.000
      C(2,1)=-A(1)*SB*B(6)                                              0609.000
      C(3,1)= A(1)*A(5)                                                 0610.000
      C(2,2)= A(2)*SA                                                   0611.000
      C(3,2)= A(2)*A(4)                                                 0612.000
      C(3,3)= A(3)                                                      0613.000
      DO 2 I=1,3                                                        0614.000
    2 D(I,I)=1.0/C(I,I)                                                 0615.000
      D(2,1)=-C(2,1)*D(1,1)*D(2,2)                                      0616.000
      D(3,2)=-C(3,2)*D(2,2)*D(3,3)                                      0617.000
      D(3,1)=-D(3,3)*(C(3,1)*D(1,1)+C(3,2)*D(2,1))                      0618.000
      RETURN                                                            0619.000
      END                                                               0620.000
C-----------------------------------------------------------------------0621.000
C     MATRIX MULTIPLICATION C(3,3) = A(3,3) * B(3,3)                    0622.000
      SUBROUTINE MM(A,B,C)                                              0623.000
      DIMENSION A(3,3),B(3,3),C(3,3)                                    0624.000
      DO 1 I=1,3                                                        0625.000
      DO 1 J=1,3                                                        0626.000
      C(I,J)=0.0                                                        0627.000
      DO 2 K=1,3                                                        0628.000
    2 C(I,J)=C(I,J)+A(I,K)*B(K,J)                                       0629.000
    1 CONTINUE                                                          0630.000
      RETURN                                                            0631.000
      END                                                               0632.000
