      PROGRAM SPHERABS
C***********************************************************************
C*                                                                     *
C* SPHERABS: a FORTRAN77 program for calculating X-ray absorption for  *
C* a spherical single crystal sample.                                  *
C*                                                                     *
C*                           by M.NARDELLI                             *
C*                                                                     *
C* Dipartimento di Chimica Generale ed Inorganica della Universita`    * 
C* degli Studi di Parma, Centro di Studio per la Strutturistica        * 
C* Diffrattometrica del C.N.R., Viale delle Scienze, I-43100 Parma,    * 
C* Italy.                                                              *
C*                                                                     *
C* The program calculates: the absorption correction by interpolating  *
C* data from a table giving the absorption correction factors. The     *
C* table is that for spherical crystals: Int.Tab. vol. C (1992) p. 523;*
C* Tab. 6.3.3.3 for values of muR in the range 0-2.5, and Int.Tab.     *
C* vol.II (1959) p.302; Table 5.3.6 B for muR in the range 2.6-10.0.   *   
C* For cylindrical crystals the table of the correction factors storag-*
C* ged in the data instruction must be changed accordingly.            *
C* The interpolation procedure illustrated by C.W.Dwiggins jr (Acta    *
C* Cryst.(1975).A31,146-148) is used with some modification.           *
C* The program corrects for absorption the observed squared structure  *
C* factors or structure factors in the input table for SHELXL-93 or    *
C* SHELX- 76, respectively (FORMAT(3I4,2F8.2,3X,I1).                   *
C* The muR value, the absorption correction as a function of theta,    *
C* the max., min. and mean values of transmission factors, error       *
C* messages (if any) are quoted in a file named SPHERABS.RIS           *
C* The computer interactively asks :                                   *
C*         the name of the input file (the same as for SHELX76)        *
C*         the name of the output file (input for SHELX76)             *
C*         the name of the compound                                    *
C*         a,b,c,alpha,beta,gamma                                      *
C*         the muR value                                               *
C*         the code for F^2 or F                                       *
C*         the value of lambda (wave-length of the radiation)          *
C* Subroutines used:                                                   *
C*         RCCNST                                                      *
C*         COSTK(X1,X2,X3,Y1,Y2,Y3,AK0,AK1,AK2)                        *
C*         COSTL(T1,T2,T3,T4,Y1,Y2,Y3,Y4,AL0,AL1,AL2,AL3)              *
C***********************************************************************
      CHARACTER*1 TIT(80)
      CHARACTER*14 FILIN,FILOUT
      INTEGER H
      DOUBLE PRECISION ST,AL0,AL1,AL2,AL3,SN(3),CS(3),CR(3),CSR(3)
      DIMENSION C(3),A(3),AF(101,19),RM(101),
     1 TH(19),AS(19),AK(3,19)
      COMMON/GRP1/RD
     1      /GRP2/C,A
     2      /GRP3/SN,CS,CR,CSR
      DATA((AF(I,K),I=1,101),K=1,2)/1.0,1.1609,1.3457,1.5574,1.7994,
     1 2.0755,2.3897,2.7467,3.1511,3.6082,4.1237,4.7035,5.3542,6.082,
     2 6.895,7.801,8.806,9.920,11.151,12.507,13.998,15.632,17.419,
     3 19.369,21.489,23.791,26.3,29.0,31.9,35.0,38.4,42.0,45.8,49.9,
     4 54.3,58.9,63.8,69.0,74.6,80.4,86.5,93.0,99.8,107,114,122,130,
     5 139,148,157,167,177,188,199,210,222,234,247,260,274,288,303,318,
     6 334,350,366,383,401,419,438,457,477,498,519,540,563,585,609,633,
     7 657,683,709,735,762,790,819,848,878,909,940,972,1005,1038,1072,
     8 1107,1143,1180,1217,1255,1294,1333,1.0,1.1609,1.3456,1.5571,
     9 1.7988,2.0743,2.3877,2.7434,3.1461,3.6009,4.1131,4.6886,5.3335,
     1 6.054,6.857,7.750,8.740,9.834,11.040,12.366,13.819,15.408,
     2 17.141,19.025,21.069,23.280,25.6,28.2,30.9,33.9,37.0,40.4,43.9,
     3 47.7,51.7,56.0,60.4,65.1,70.1,75.3,80.7,86.4,92.4,98.6,105,112,
     4 119,126,134,141,150,158,167,176,185,194,204,214,224,235,246,257,
     5 269,280,293,305,318,331,344,357,371,385,400,415,430,445,460,476,
     6 493,509,526,543,560,578,596,614,632,651,670,690,709,729,749,770,
     7 791,812,833,855,877,899,921/
      DATA((AF(I,K),I=1,101),K=3,4)/ 1.0,1.1609,1.3452,1.5561,1.7968,
     8 2.0706,2.3816,2.7336,3.1312,3.5789,4.0815,4.6442,5.2722,5.9710,
     9 6.746,7.604,8.549,9.587,10.725,11.967,13.320,14.788,16.376,
     1 18.089,19.931,21.907,24.0,26.2,28.6,31.2,33.9,36.7,39.7,42.8,
     2 46.0,49.5,53.0,56.8,60.6,64.6,68.8,73.1,77.5,82.1,86.8,91.7,96.7,
     3 102,107,112,118,124,129,135,141,147,153,160,166,173,179,186,193,
     4 200,207,214,221,229,236,244,251,259,267,275,283,291,299,307,315,
     5 324,332,341,349,358,367,375,384,393,402,411,420,429,439,448,457,
     6 467,476,485,495,505,514,1.0,1.1607,1.3447,1.5546,1.7935,2.0647,
     7 2.3715,2.7177,3.1069,3.5431,4.0304,4.5729,5.1747,5.8399,6.573,
     8 7.377,8.256,9.214,10.254,11.380,12.593,13.895,15.290,16.778,
     9 18.361,20.040,21.8,23.7,25.6,27.7,29.9,32.1,34.4,36.8,39.3,41.9,
     1 44.6,47.3,50.2,53.1,56.0,59.1,62.2,65.3,68.6,71.9,75.2,78.6,82.1,
     2 85.6,89.1,92.7,96.4,100,104,108,111,115,119,123,127,131,135,139,
     3 143,148,152,156,160,165,169,173,177,182,186,191,195,200,204,209,
     4 213,218,222,227,232,236,241,246,250,255,260,264,269,274,279,284,
     5 288,293,298,303,308/
      DATA((AF(I,K),I=1,101),K=5,6)/ 1.0,1.1606,1.3439,1.5525,1.7891,
     8 2.0565,2.3578,2.6959,3.0740,3.4952,3.9625,4.4790,5.0476,5.6710,
     9 6.352,7.092,7.894,8.759,9.689,10.685,11.746,12.874,14.067,
     1 15.327,16.652,18.041,19.5,21.0,22.6,24.2,25.9,27.7,29.5,31.3,
     2 33.2,35.2,37.1,39.2,41.2,43.3,45.5,47.6,49.9,52.1,54.4,56.7,59.0,
     3 61.3,63.7,66.1,68.5,71.0,73.4,75.9,78.4,80.9,83.5,86.0,88.6,91.2,
     4 93.8,96.4,99.1,102,104,107,110,112,115,118,121,123,126,129,132,
     5 135,137,140,143,146,149,152,154,157,160,163,166,169,172,175,178,
     6 180,183,186,189,192,195,198,201,204,207,1.0,1.1603,1.3428,1.5497,
     7 1.7833,2.0462,2.3406,2.6691,3.0339,3.4374,3.8816,4.3686,4.9001,
     8 5.4776,6.102,6.775,7.497,8.268,9.088,9.957,10.873,11.837,12.847,
     9 13.902,15.000,16.142,17.3,18.6,19.8,21.1,22.4,23.8,25.2,26.6,
     1 28.1,29.5,31.0,32.6,34.1,35.7,37.2,38.8,40.5,42.1,43.7,45.4,47.1,
     2 48.8,50.5,52.2,53.9,55.6,57.4,59.2,60.9,62.7,64.5,66.3,68.1,69.9,
     3 71.7,73.6,75.4,77.2,79.1,81.0,82.8,84.7,86.6,88.5,90.4,92.2,94.1,
     4 96.0,98.0,99.9,102,104,106,108,110,111,113,115,117,119,121,123,
     5 125,127,129,131,133,135,137,139,141,143,145,147,149/
      DATA((AF(I,K),I=1,101),K=7,8)/ 1.0,1.1600,1.3415,1.5463,1.7765,
     8 2.0340,2.3206,2.6382,2.9882,3.3723,3.7917,4.2474,4.7404,5.2711,
     9 5.8400,6.447,7.092,7.774,8.492,9.246,10.034,10.855,11.708,
     1 12.592,13.504,14.445,15.4,16.4,17.4,18.5,19.5,20.6,21.7,22.8,
     2 23.9,25.1,26.2,27.4,28.6,29.8,31.0,32.2,33.4,34.7,35.9,37.2,38.4,
     3 39.7,41.0,42.2,43.5,44.8,46.1,47.4,48.8,50.1,51.4,52.7,54.1,55.4,
     4 56.8,58.1,59.5,60.8,62.2,63.5,64.9,66.3,67.7,69.0,70.4,71.8,73.2,
     5 74.6,76.0,77.4,78.8,80.2,81.6,83.0,84.4,85.8,87.2,88.6,90.1,91.5,
     6 92.9,94.3,95.8,97.2,98.6,100,101,103,104,106,107,109,110,112,113,
     7 1.0,1.1597,1.3400,1.5426,1.7689,2.0204,2.2984,2.6042,2.9386,
     8 3.3026,3.6966,4.1211,4.5761,5.0617,5.5774,6.123,6.697,7.299,
     9 7.928,8.583,9.262,9.964,10.688,11.433,12.198,12.982,13.8,14.6,
     1 15.4,16.3,17.1,18.0,18.9,19.8,20.7,21.6,22.5,23.4,24.4,25.3,26.3,
     2 27.2,28.2,29.2,30.1,31.1,32.1,33.1,34.1,35.1,36.1,37.1,38.1,39.1,
     3 40.1,41.1,42.2,43.2,44.2,45.2,46.3,47.3,48.4,49.4,50.4,51.5,52.5,
     4 53.6,54.6,55.7,56.7,57.8,58.9,59.9,61.0,62.0,63.1,64.2,65.2,66.3,
     5 67.4,68.4,69.5,70.6,71.7,72.7,73.8,74.9,76.0,77.1,78.2,79.3,80.4,
     6 81.5,82.6,83.7,84.8,85.9,87.0,88.2,89.3/
      DATA((AF(I,K),I=1,101),K=9,10)/ 1.0,1.1593,1.3383,1.5383,1.7604,
     8 2.0056,2.2746,2.5683,2.8869,3.3208,3.6001,3.9945,4.4137,4.8573,
     9 5.3244,5.8143,6.326,6.859,7.411,7.982,8.570,9.175,9.795,10.429,
     1 11.077,11.738,12.4,13.1,13.8,14.5,15.2,15.9,16.6,17.4,18.1,18.9,
     2 19.6,20.4,21.1,21.9,22.7,23.4,24.2,25.0,25.8,26.6,27.4,28.1,28.9,
     3 29.7,30.5,31.3,32.2,33.0,33.8,34.6,35.4,36.2,37.0,37.9,38.7,39.5,
     4 40.3,41.2,42.0,42.8,43.6,44.5,45.3,46.1,47.0,47.8,48.7,49.5,50.3,
     5 51.2,52.0,52.8,53.7,54.5,55.4,56.2,57.1,57.9,58.8,59.6,60.4,61.3,
     6 62.1,63.0,63.9,64.7,65.6,66.4,67.3,68.2,69.1,69.9,70.8,71.7,72.6,
     7 1.0,1.1589,1.3366,1.5339,1.7515,1.9901,2.2500,2.5316,2.8347,
     8 3.1592,3.5048,3.8710,4.2571,4.6625,5.0862,5.5273,5.9849,6.458,
     9 6.946,7.447,7.961,8.486,9.023,9.569,10.125,10.690,11.3,11.8,
     1 12.4,13.0,13.6,14.2,14.8,15.4,16.1,16.7,17.3,17.9,18.6,19.2,19.8,
     2 20.5,21.1,21.8,22.4,23.1,23.7,24.4,25.0,25.7,26.4,27.0,27.7,28.3,
     3 29.0,29.7,30.3,31.0,31.7,32.3,33.0,33.7,34.4,35.0,35.7,36.4,37.1,
     4 37.8,38.4,39.1,39.8,40.5,41.2,41.8,42.5,43.2,43.9,44.6,45.3,45.9,
     5 46.6,47.3,48.0,48.7,49.4,50.1,50.7,51.4,52.1,52.8,53.5,54.2,54.9,
     6 55.6,56.3,57.0,57.7,58.4,59.1,59.8,60.5/
      DATA((AF(I,K),I=1,101),K=11,12)/ 1.0,1.1586,1.3348,1.5293,1.7425,
     8 1.9745,2.2255,2.4952,2.7835,3.0898,3.4135,3.7540,4.1104,4.4819,
     9 4.8676,5.2666,5.6780,6.101,6.535,6.978,7.431,7.893,8.362,8.839,
     1 9.322,9.810,10.3,10.8,11.3,11.8,12.3,12.8,13.4,13.9,14.4,14.9,
     2 15.5,16.0,16.5,17.1,17.6,18.2,18.7,19.3,19.8,20.3,20.9,21.4,22.0,
     3 22.6,23.1,23.7,24.2,24.8,25.3,25.9,26.4,27.0,27.6,28.1,28.7,29.3,
     4 29.8,30.4,31.0,31.5,32.1,32.7,33.2,33.8,34.4,34.9,35.5,36.1,36.6,
     5 37.2,37.8,38.4,38.9,39.5,40.1,40.6,41.2,41.8,42.4,42.9,43.5,44.1,
     6 44.7,45.2,45.8,46.4,46.9,47.5,48.0,48.6,49.1,49.6,50.2,50.7,51.2,
     7 1.0,1.1582,1.3331,1.5248,1.7335,1.9592,2.2015,2.4602,2.7346,
     8 3.0241,3.3280,3.6455,3.9756,4.3175,4.6703,5.0333,5.4057,5.7867,
     9 6.176,6.572,6.975,7.385,7.800,8.220,8.645,9.074,9.51,9.94,
     1 10.4,10.8,11.3,11.7,12.2,12.6,13.1,13.5,14.0,14.5,14.9,15.4,15.8,
     2 16.3,16.8,17.2,17.7,18.2,18.7,19.1,19.6,20.1,20.6,21.0,21.5,22.0,
     3 22.5,22.9,23.4,23.9,24.4,24.9,25.3,25.8,26.3,26.8,27.3,27.8,28.2,
     4 28.7,29.2,29.7,30.2,30.7,31.1,31.6,32.1,32.6,33.1,33.6,34.1,34.6,
     5 35.0,35.5,36.0,36.5,37.0,37.5,38.0,38.5,39.0,39.5,39.9,40.4,40.9,
     6 41.3,41.8,42.2,42.7,43.1,43.5,43.9,44.2/
      DATA((AF(I,K),I=1,101),K=13,14)/ 1.0,1.1579,1.3313,1.5204,1.7249,
     8 1.9445,2.1789,2.4274,2.6892,2.9637,3.2499,3.5470,3.8542,4.1706,
     9 4.4955,4.8281,5.1678,5.5140,5.8662,6.224,6.587,6.955,7.327,7.702,
     1 8.081,8.462,8.85,9.23,9.62,10.0,10.4,10.8,11.2,11.6,12.0,12.4,
     2 12.8,13.2,13.6,14.0,14.4,14.8,15.2,15.7,16.1,16.5,16.9,17.3,17.7,
     3 18.1,18.5,19.0,19.4,19.8,20.2,20.6,21.0,21.5,21.9,22.3,22.7,23.1,
     4 23.6,24.0,24.4,24.8,25.2,25.7,26.1,26.5,26.9,27.3,27.8,28.2,28.6,
     5 29.0,29.4,29.9,30.3,30.7,31.1,31.6,32.0,32.4,32.9,33.3,33.7,34.1,
     6 34.6,35.0,35.4,35.8,36.2,36.6,37.0,37.4,37.7,38.1,38.4,38.7,39.0,
     7 1.0,1.1575,1.3297,1.5162,1.7169,1.9311,2.1583,2.3977,2.6484,
     8 2.9098,3.1807,3.4605,3.7483,4.0432,4.3447,4.6520,4.9647,5.2823,
     9 5.6045,5.9308,6.261,6.595,6.932,7.272,7.614,7.957,8.30,8.64,
     1 8.99,9.34,9.70,10.0,10.4,10.8,11.1,11.5,11.8,12.2,12.6,12.9,13.3,
     2 13.6,14.0,14.4,14.7,15.1,15.5,15.8,16.2,16.6,16.9,17.3,17.7,18.0,
     3 18.4,18.8,19.1,19.5,19.9,20.3,20.6,21.0,21.4,21.7,22.1,22.5,22.9,
     4 23.2,23.6,24.0,24.3,24.7,25.1,25.5,25.8,26.2,26.6,27.0,27.3,27.7,
     5 28.1,28.5,28.8,29.2,29.6,30.0,30.3,30.7,31.1,31.5,31.8,32.2,32.6,
     6 32.9,33.3,33.6,33.9,34.3,34.6,34.9,35.2/
      DATA((AF(I,K),I=1,101),K=15,16)/ 1.0,1.1572,1.3282,1.5126,1.7099,
     8 1.9194,2.1403,2.3719,2.6133,2.8634,3.1216,3.3870,3.6586,3.9360,
     9 4.2183,4.5052,4.7961,5.0907,5.3888,5.6900,5.9942,6.301,6.610,
     1 6.922,7.235,7.548,7.86,8.17,8.49,8.81,9.13,9.45,9.77,10.1,10.4,
     2 10.7,11.1,11.4,11.7,12.0,12.4,12.7,13.0,13.4,13.7,14.0,14.3,14.7,
     3 15.0,15.3,15.7,16.0,16.3,16.7,17.0,17.3,17.7,18.0,18.3,18.7,19.0,
     4 19.3,19.7,20.0,20.3,20.7,21.0,21.3,21.7,22.0,22.3,22.7,23.0,23.4,
     5 23.7,24.0,24.4,24.7,25.0,25.4,25.7,26.0,26.4,26.7,27.1,27.4,27.7,
     6 28.1,28.4,28.7,29.1,29.4,29.7,30.1,30.4,30.7,31.0,31.4,31.7,32.0,
     7 32.3,1.0,1.1570,1.3271,1.5096,1.7041,1.9097,1.1256,2.3508,2.5845,
     8 2.8258,3.0738,3.3276,3.5866,3.8500,4.1174,4.3883,4.6622,4.9390,
     9 5.2184,5.5001,5.7842,6.070,6.358,6.648,6.938,7.229,7.52,7.81,
     1 8.10,8.40,8.70,8.99,9.29,9.59,9.88,10.2,10.5,10.8,11.1,11.4,11.7,
     2 12.0,12.3,12.6,12.9,13.2,13.5,13.8,14.1,14.4,14.7,15.0,15.3,15.6,
     3 15.9,16.2,16.5,16.8,17.2,17.5,17.8,18.1,18.4,18.7,19.0,19.3,19.6,
     4 19.9,20.2,20.5,20.8,21.1,21.4,21.7,22.1,22.4,22.7,23.0,23.3,23.6,
     5 23.9,24.2,24.5,24.8,25.1,25.4,25.7,26.1,26.4,26.7,27.0,27.3,27.6,
     6 27.9,28.2,28.5,28.8,29.1,29.5,29.8,30.1/
      DATA((AF(I,K),I=1,101),K=17,19)/ 1.0,1.1568,1.3262,1.5074,1.6997,
     8 1.9024,2.1145,2.3351,2.5632,2.7979,3.0383,3.2838,3.5334,3.7868,
     9 4.0432,4.3024,4.5641,4.8279,5.0936,5.3613,5.6307,5.9017,6.174,
     1 6.448,6.722,6.996,7.27,7.54,7.82,8.10,8.38,8.66,8.94,9.21,9.49,
     2 9.77,10.1,10.3,10.6,10.9,11.2,11.5,11.7,12.0,12.3,12.6,12.9,13.2,
     3 13.4,13.7,14.0,14.3,14.6,14.9,15.1,15.4,15.7,16.0,16.3,16.6,16.8,
     4 17.1,17.4,17.7,18.0,18.3,18.5,18.8,19.1,19.4,19.7,20.0,20.3,20.5,
     5 20.8,21.1,21.4,21.7,22.0,22.3,22.5,22.8,23.1,23.4,23.7,24.0,24.3,
     6 24.5,24.8,25.1,25.4,25.7,26.0,26.3,26.5,26.8,27.1,27.4,27.7,28.0,
     7 28.3,1.0,1.1567,1.3256,1.5059,1.6970,1.8979,2.1076,2.3253,2.5499,
     8 2.7805,3.0163,3.2566,3.5005,3.7477,3.9974,4.2495,4.5036,4.7595,
     9 5.0170,5.2760,5.5365,5.7982,6.061,6.325,6.589,6.853,7.12,7.38,
     1 7.65,7.92,8.18,8.45,8.72,8.99,9.25,9.52,9.79,10.1,10.3,10.6,10.9,
     2 11.1,11.4,11.7,11.9,12.2,12.5,12.8,13.0,13.3,13.6,13.8,14.1,14.4,
     3 14.6,14.9,15.2,15.5,15.7,16.0,16.3,16.5,16.8,17.1,17.4,17.6,17.9,
     4 18.2,18.4,18.7,19.0,19.3,19.5,19.8,20.1,20.3,20.6,20.9,21.2,21.4,
     5 21.7,22.0,22.2,22.5,22.8,23.1,23.3,23.6,23.9,24.2,24.4,24.7,25.0,
     6 25.2,25.5,25.8,26.1,26.3,26.6,26.9,27.2,1.0,1.1567,1.3254,1.5055,
     7 1.6961,1.8964,2.1063,2.3220,2.5454,2.7747,3.0090,3.2474,3.4894,
     8 3.7344,3.9819,4.2315,4.4830,4.7361,4.9908,5.2468,5.5041,5.7627,
     9 6.022,6.282,6.543,6.803,7.06,7.33,7.59,7.85,8.11,8.38,8.64,8.90,
     1 9.17,9.43,9.69,9.96,10.2,10.5,10.8,11.0,11.3,11.5,11.8,12.1,12.3,
     2 12.6,12.9,13.1,13.4,13.7,13.9,14.2,14.5,14.7,15.0,15.3,15.5,15.8,
     3 16.1,16.3,16.6,16.9,17.1,17.4,17.7,17.9,18.2,18.4,18.7,19.0,19.2,
     4 19.5,19.8,20.0,20.3,20.6,20.8,21.1,21.4,21.6,21.9,22.2,22.4,22.7,
     5 23.0,23.2,23.5,23.8,24.0,24.3,24.6,24.8,25.1,25.4,25.6,25.9,26.2,
     6 26.4,26.7/
      DATA RM/ 0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,
     1 1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,
     2 3.0,3.1,3.2,3.3,3.4,3.5,3.6,3.7,3.8,3.9,4.0,4.1,4.2,4.3,4.4,
     3 4.5,4.6,4.7,4.8,4.9,5.0,5.1,5.2,5.3,5.4,5.5,5.6,5.7,5.8,5.9,6.0,
     4 6.1,6.2,6.3,6.4,6.5,6.6,6.7,6.8,6.9,7.0,7.1,7.2,7.3,7.4,7.5,7.6,
     5 7.7,7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5,8.6,8.7,8.8,8.9,9.0,9.1,9.2,
     6 9.3,9.4,9.5,9.6,9.7,9.8,9.9,10.0/
      DATA TH/.0,.08727,.17453,.26180,.34907,.43633,.52360,.61087,
     1 .69813,.78540,.87266,0.95993,1.04720,1.13446,1.22173,1.30900,
     2 1.39626,1.48353,1.57080/
      RD=0.017453293
      WRITE(*,100)
100   FORMAT(' Key in the input file'/' ?>')
      READ(*,'(A14)')FILIN
      WRITE(*,101)
101   FORMAT(' Key in the output file'/' ?>')
      READ(*,'(A14)')FILOUT
      OPEN(UNIT=13,FILE=FILIN,FORM='FORMATTED')
      REWIND(UNIT=13)
      OPEN(UNIT=15,FILE=FILOUT,FORM='FORMATTED')
      REWIND(UNIT=15)
      OPEN(UNIT=16,FILE='SPHERABS.RIS',FORM='FORMATTED')
      REWIND(UNIT=16)
C-----Read title, parameters
      WRITE(*,102)
102   FORMAT(' Key in the name of the compound'/' ?>')
      READ(*,'(80A1)') TIT
      WRITE(16,'(1X,80A1)') TIT
      WRITE(*,103)
103   FORMAT(' Key in a,b,c,alpha,beta,gamma'/' ?: ')
      READ(*,*) C,A
C-----Calculate reciprocal constants
      CALL RCCNST
C-----Read the muR value and choose the standard muR values for
C-----interpolation
      WRITE(*,104)
104   FORMAT(' Key in the muR value'/' ?: ')
      READ(*,*)X0
      WRITE(*,111)
111   FORMAT(' Key in the code for structure factors'/
     1 '   0  if F^2'/'   1  if F'/' ?:>')
      READ(*,*)KF
      IF(X0.LT.0.1.OR.X0.GT.10) GO TO 1001
      DO 1 I=1,101
      IF(X0.EQ.RM(I+1)) GO TO 2
      IF(X0.GT.RM(I+1)) GO TO 1
      X1=RM(I-1)
      X2=RM(I)
      X3=RM(I+1)
      DO 3 K=1,19
      Y1=AF(I-1,K)
      Y2=AF(I,K)
      Y3=AF(I+1,K)
      CALL COSTK(X1,X2,X3,Y1,Y2,Y3,AK0,AK1,AK2)
3     AS(K)=EXP(AK0+AK1*X0+AK2*X0**2)
      GO TO 9
1     CONTINUE
      WRITE(15,'(''There is something wrong with your data'')')
      STOP
2     DO 6 K=1,19
6     AS(K)=AF(I+1,K)
9     WRITE(16,'('' muR = '',F7.4/'' Theta'',9X,''A*'')')X0
      SAS=0.0
      DO 200 I=1,19
      SAS=SAS+AS(I)
      THD=TH(I)/RD
200   WRITE(16,'(1X,F5.2,5X,F9.4)')THD,AS(I)
      SA=SAS/19
C-----Read wavelength
      WRITE(*,105)
105   FORMAT(' Key in the wavelength'/' ?: ')
      READ(*,*) WL
C      IF(KF.EQ.0)WRITE(16,110)
C110   FORMAT(/18X,'Fo^2',4X,'Fo^2(corr.)'/)
C      IF(KF.EQ.1)WRITE(16,112)
C112   FORMAT(/16X,'Fo',1X,'Fo(corr.)'/)
C-----Read indices and structure factors
      TFMAX=0.0
      TFMIN=10000.0
      NN=0
      SS=0
      TH(1)=0.043633
      AS0=AS(1)
      AS(1)=(AS(1)+AS(2))/2.0
1000  CONTINUE
      READ(13,106)H,K,L,FO,SF,IN
106   FORMAT(3I4,2F8.2,3X,I1)
      IF(H.EQ.0.AND.K.EQ.0.AND.L.EQ.0) GO TO 1002
C-----Calculate sin(theta) and theta
      SNT=(WL/2.)*DSQRT((H*CR(1))**2+(K*CR(2))**2+(L*CR(3))**2+
     1        2.*(K*L*CR(2)*CR(3)*CSR(1)+
     2            L*H*CR(3)*CR(1)*CSR(2)+
     3            H*K*CR(1)*CR(2)*CSR(3)))
      IF(SNT.GT.1.0) WRITE(16,504)SNT
504   FORMAT(' sin(T0)=',F8.5)
      IF(SNT.GT.1.0) GO TO 1000
      T0=ASIN(SNT)
      ST=SNT
C-----Choose the standard theta values for interpolation
      DO 10 I=1,19
      KK=I
      IF(ABS(T0-TH(I)).LE.0.0001) GO TO 15
      IF(T0.GT.TH(2))GO TO 11
      T1=TH(1)
      T2=TH(2)
      T3=TH(3)
      T4=TH(4)
      Y1=AS(1)
      Y2=AS(2)
      Y3=AS(3)
      Y4=AS(4)
      GO TO 14
11    IF(T0.LT.TH(18))GO TO 12
      T1=TH(16)
      T2=TH(17)
      T3=TH(18)
      T4=TH(19)
      Y1=AS(16)
      Y2=AS(17)
      Y3=AS(18)
      Y4=AS(19)
      GO TO 14
12    IF(T0.LT.TH(I).OR.T0.GT.TH(I+1))GO TO 10
      IF(ABS(T0-TH(I)).LE.0.0001)THEN
      KK=I
      GO TO 15
      ENDIF
      IF(ABS(T0-TH(I+1)).LE.0.0001)THEN
      KK=I+1
      GO TO 15
      ENDIF
      IF(T0.GT.TH(2).AND.T0.LT.TH(3))THEN
      T1=TH(2)
      T2=TH(3)
      T3=TH(4)
      T4=TH(5)
      Y1=AS(2)
      Y2=AS(3)
      Y3=AS(4)
      Y4=AS(5)
      GO TO 14
      ENDIF
      T1=TH(I-1)
      T2=TH(I)
      T3=TH(I+1)
      T4=TH(I+2)
      Y1=AS(I-1)
      Y2=AS(I)
      Y3=AS(I+1)
      Y4=AS(I+2)
      GO TO 14
10    CONTINUE
14    CALL COSTL(T1,T2,T3,T4,Y1,Y2,Y3,Y4,AL0,AL1,AL2,AL3)
      ST=ST**2
      AC=(AL0+AL1*ST+AL2*ST**2+AL3*ST**3)
      GO TO 16
15    AC=AS(KK)
16    IF(KF.EQ.1)AC=SQRT(AC)
      IF(TFMAX.LT.1/AC)TFMAX=1/AC
      IF(TFMIN.GT.1/AC)TFMIN=1/AC
      T0GRD=T0/RD
      IF(AC.GT.AS0)WRITE(16,'('' AC is greater than AS(1)'')')
C     WRITE(16,'(''T0='',F6.2,3X,'' AC='',F8.4)')T0GRD,AC
      FCR=FO*AC/AS0
      SFC=SF*AC/AS0
      WRITE(15,107)H,K,L,FCR,SFC,IN
107   FORMAT(3I4,2F8.2,3X,I1)
C     WRITE(16,108)H,K,L,FO,FCR
C108  FORMAT(3I4,2F12.2)
      SS=SS+AC
      NN=NN+1
      GO TO 1000
1002  WRITE(15,106)H,K,L,FO,SF,IN
      SS=NN/SS
C     TFMAX=TFMAX/SS
C     TFMIN=TFMIN/SS
      WRITE(16,114)TFMAX,TFMIN,SS
114   FORMAT(' Transmission factor max. = ',F9.5/
     1 ' Transmission factor min. = ', F9.5/
     2 ' Transmission factor mean = ', F9.5)
      GO TO 2000
1001  WRITE(*,109)
109   FORMAT(' No calculation can be made with your muR value!')
2000  STOP
      END
 
      SUBROUTINE  RCCNST
C-----Calcola le costanti reciproche a*,b*,c*,cos(alpha*),cos(beta*),
C-----cos(gamma*) da quelle dirette.
      DOUBLE PRECISION SN(3),CS(3),CR(3),CSR(3),AD(3),S,VD,AA
      DIMENSION C(3),A(3)
      COMMON/GRP1/RD
     1      /GRP2/C,A
     2      /GRP3/SN,CS,CR,CSR
      S=0D0
      DO 1 I=1,3
      AD(I)=A(I)*RD
      SN(I)=DSIN(AD(I))
      CS(I)=DCOS(AD(I))
1     S=S+AD(I)
      S=S/2D0
      VD=2*C(1)*C(2)*C(3)
      AA=DSIN(S)
      DO 2 I=1,3
2     AA=AA*DSIN(S-AD(I))
      VD=VD*DSQRT(AA)
      CR(1)=C(2)*C(3)*SN(1)/VD
      CR(2)=C(3)*C(1)*SN(2)/VD
      CR(3)=C(1)*C(2)*SN(3)/VD
      CSR(1)=(CS(2)*CS(3)-CS(1))/(SN(2)*SN(3))
      CSR(2)=(CS(3)*CS(1)-CS(2))/(SN(3)*SN(1))
      CSR(3)=(CS(1)*CS(2)-CS(3))/(SN(1)*SN(2))
      RETURN
      END
 
      SUBROUTINE COSTK(X1,X2,X3,Y1,Y2,Y3,AK0,AK1,AK2)
C-----This routine solves a system of three linear equations
      Y1=ALOG(Y1)
      Y2=ALOG(Y2)
      Y3=ALOG(Y3)
      AA=(X2*X3**2-X3*X2**2)
      BB=(X3*X1**2-X1*X3**2)
      CC=(X1*X2**2-X2*X1**2)
      D1=AA+BB+CC
      AK0=(Y1*AA+Y2*BB+Y3*CC)/D1
      AK1=((Y1-AK0)*X2**2-(Y2-AK0)*X1**2)/CC
      AK2=(Y1-AK0-AK1*X1)/X1**2
      RETURN
      END
 
      SUBROUTINE COSTL(T1,T2,T3,T4,Y1,Y2,Y3,Y4,AL0,AL1,AL2,AL3)
C-----This routine solves a system of four linear equations
      COMMON/GRP1/RD
      DOUBLE PRECISION X1,X2,X3,X4,A1,A2,A3,A4,A5,A6,B1,B2,B3,B4,
     1 D,DY1,DY2,DY3,DY4,AL0,AL1,AL2,AL3,DT1,DT2,DT3,DT4
      DY1=Y1
      DY2=Y2
      DY3=Y3
      DY4=Y4
      DT1=T1
      DT2=T2
      DT3=T3
      DT4=T4
      X1=DSIN(DT1)**2
      X2=DSIN(DT2)**2
      X3=DSIN(DT3)**2
      X4=DSIN(DT4)**2
      A1=X3**2*X4**2*(X4-X3)
      A2=X2**2*X4**2*(X4-X2)
      A3=X2**2*X3**2*(X3-X2)
      A4=X1**2*X4**2*(X4-X1)
      A5=X1**2*X3**2*(X3-X1)
      A6=X1**2*X2**2*(X2-X1)
      B1=X2*A1-X3*A2+X4*A3
      B2=X1*A1-X3*A4+X4*A5
      B3=X1*A2-X2*A4+X4*A6
      B4=X1*A3-X2*A5+X3*A6
      D=B1-B2+B3-B4
      AL0=(DY1*B1-DY2*B2+DY3*B3-DY4*B4)/D
      AL1=((DY2-AL0)*A1-(DY3-AL0)*A2+(DY4-AL0)*A3)/B1
      AL2=((DY3-AL0-AL1*X3)*X4**3-(DY4-AL0-AL1*X4)*X3**3)/A1
      AL3=(DY2-AL0-AL1*X2-AL2*X2**2)/X2**3
      RETURN
      END
