Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:07:30

0001 C==========================================================================
0002 C.  Library of programs for the generation of nucleus-nucleus interactions
0003 C.  and the study of nucleus-induced cosmic ray showers
0004 C.
0005 C.  September 2001  changes in FPNI, and SIGMA_INI,
0006 C.                  new SIGMA_PP, SIGMA_PPI, SIGMA_KP  (R. Engel)
0007 C.
0008 C.  may  1996       small bug  corrected by Dieter Heck in NUC_CONF
0009 C.
0010 C.  march 1996      small modification to the superposition code
0011 C.
0012 C.  February 1996   change to FPNI to give an interaction length
0013 C.                   also  at very low energy
0014 C.
0015 C.  Version 1.01  september 1995
0016 C.       (small corrections P.L.)
0017 C.       the random number generator is called as S_RNDM(0)
0018 C.  ------------------------------------------------------
0019 C.  Version 1.00  April 1992
0020 C.
0021 C.  Authors:
0022 C.
0023 C.     J. Engel
0024 C.     T.K Gaisser
0025 C.     P.Lipari
0026 C.     T. Stanev
0027 C.
0028 C.  This set of routines  when used in  the simulation of cosmic ray
0029 C.  showers have only three  "contact points" with the "external world"
0030 C.
0031 C.    (i) SUBROUTINE NUC_NUC_INI
0032 C.        (no  calling arguments)
0033 C.         to be called once during general initialization
0034 C.    (ii) SUBROUTINE HEAVY (IA, E0)
0035 C.         where IA (integer) is the mass number of the projectile
0036 C.         nucleus  and E0 (TeV) is the energy per nucleon
0037 C.         The output (positions of first interaction for the IA
0038 C.         nucleons of the projectile) is  contained in  the common block:
0039 C.           COMMON /C1STNC/ XX0(60),XX(60),YY(60),AX(60),AY(60)
0040 C.         In detail:
0041 C.             XX0(j)   (g cm-2) =  position of interaction
0042 C.             XX(j) (mm)    x-distance from shower axis
0043 C.             YY(j) (mm)    y-distance from shower axis
0044 C.             AX(j) (radiants)  Theta_x with respect to original direction
0045 C.             AY(j) (radiants)  Theta_y with respect to original direction
0046 C.
0047 C.    (iii)  FUNCTION FPNI (E,L)
0048 C.           Interaction length in air.
0049 C.           E (TeV) is the energy of the particle, L is the particle
0050 C.           code (NOTE: "Sibyll" codes are used : L =1-18)
0051 C.           WANRNING : The nucleus-nucleus cross section
0052 C.           tabulated in the program are "matched" to the p-Air
0053 C.           cross section calculated  with this FUNCTION, in other words
0054 C.           they are both calculated with the same input pp cross section
0055 C==========================================================================
0056 
0057       SUBROUTINE NUC_NUC_INI
0058 C...Initialization for the generation of nucleus-nucleus interactions
0059 C.  INPUT : E0 (TeV) Energy per nucleon of the beam nucleus
0060 C........................................................................
0061       SAVE
0062 
0063       CALL NUC_GEOM_INI                       ! nucleus profiles
0064       CALL SIGMA_INI                          ! initialize pp cross sections
0065 
0066       RETURN
0067       END
0068 
0069 
0070       SUBROUTINE FRAGM1 (IA,NW, NF, IAF)
0071 C...Nuclear Fragmentation
0072 C.  total dissolution of nucleus
0073 C..........................................
0074       SAVE
0075 
0076       DIMENSION IAF(60)
0077       NF = IA-NW
0078       DO J=1,NF
0079          IAF(J) = 1
0080       ENDDO
0081       RETURN
0082       END
0083 C->
0084       SUBROUTINE FRAGM2 (IA,NW, NF, IAF)
0085 C...Nuclear Fragmentation
0086 C.  Spectator in one single fragment
0087 C..........................................
0088       SAVE
0089 
0090       DIMENSION IAF(60)
0091       IF (IA-NW .GT. 0)  THEN
0092          NF = 1
0093          IAF(1) = IA-NW
0094       ELSE
0095          NF = 0
0096       ENDIF
0097       RETURN
0098       END
0099 
0100 C====================================================================
0101 C...Code of fragmentation  of spectator nucleons
0102 C.  based on Jon Engel  abrasion-ablation algorithms
0103 C====================================================================
0104 
0105       BLOCK DATA FRAG_DATA
0106       SAVE
0107 
0108 C...Data for the fragmentation of  nucleus  projectiles
0109       COMMON /FRAGMOD/A(10,10,20),AE(10,10,20),ERES(10,10),NFLAGG(10,10)
0110       DATA (NFLAGG(I, 1),I=1,10)  /
0111      +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
0112       DATA (NFLAGG(I, 2),I=1,10)  /
0113      +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
0114       DATA (NFLAGG(I, 3),I=1,10)  /
0115      +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
0116       DATA (NFLAGG(I, 4),I=1,10)  /
0117      +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
0118       DATA (NFLAGG(I, 5),I=1,10)  /
0119      +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
0120       DATA (NFLAGG(I, 6),I=1,10)  /
0121      +    0,  0,  0,  0,  0,  0,  0,  1,  1,  1 /
0122       DATA (NFLAGG(I, 7),I=1,10)  /
0123      +    1,  1,  1,  1,  1,  1,  1,  1,  1,  1 /
0124       DATA (NFLAGG(I, 8),I=1,10)  /
0125      +    1,  1,  1,  1,  1,  1,  1,  1,  1,  1 /
0126       DATA (NFLAGG(I, 9),I=1,10)  /
0127      +    1,  1,  1,  1,  1,  1,  1,  1,  1,  1 /
0128       DATA (NFLAGG(I,10),I=1,10)  /
0129      +    1,  1,  1,  1,  1,  1,  1,  1,  1,  1 /
0130       DATA (A(I, 1, 1),I=1,10)  /
0131      +  .438E-01,.172    ,.283    ,.511    ,.715    ,.920    ,1.19    ,
0132      +  1.37    ,1.65    ,2.14     /
0133       DATA (A(I, 1, 2),I=1,10)  /
0134      +  .147E-01,.249E-01,.439E-01,.592E-01,.776E-01,.886E-01,.108    ,
0135      +  .117    ,.126    ,.128     /
0136       DATA (A(I, 1, 3),I=1,10)  /
0137      +  .216E-02,.627E-02,.834E-02,.108E-01,.144E-01,.152E-01,.196E-01,
0138      +  .200E-01,.210E-01,.224E-01 /
0139       DATA (A(I, 1, 4),I=1,10)  /
0140      +  .593E-01,.653E-01,.116    ,.145    ,.184    ,.204    ,.234    ,
0141      +  .257    ,.271    ,.248     /
0142       DATA (A(I, 1, 5),I=1,10)  /
0143      +  .000E+00,.918E-02,.362E-02,.805E-02,.436E-02,.728E-02,.466E-02,
0144      +  .707E-02,.932E-02,.130E-01 /
0145       DATA (A(I, 1, 6),I=1,10)  /
0146      +  .000E+00,.180E-02,.247E-02,.208E-02,.224E-02,.214E-02,.226E-02,
0147      +  .233E-02,.230E-02,.194E-02 /
0148       DATA (A(I, 1, 7),I=1,10)  /
0149      +  .000E+00,.106E-02,.703E-03,.687E-03,.739E-03,.674E-03,.819E-03,
0150      +  .768E-03,.756E-03,.720E-03 /
0151       DATA (A(I, 1, 8),I=1,10)  /
0152      +  .000E+00,.000E+00,.188E-02,.130E-02,.138E-02,.117E-02,.124E-02,
0153      +  .119E-02,.111E-02,.829E-03 /
0154       DATA (A(I, 1, 9),I=1,10)  /
0155      +  .000E+00,.000E+00,.302E-03,.258E-03,.249E-03,.208E-03,.248E-03,
0156      +  .222E-03,.210E-03,.187E-03 /
0157       DATA (A(I, 1,10),I=1,10)  /
0158      +  .000E+00,.000E+00,.000E+00,.235E-03,.222E-03,.172E-03,.181E-03,
0159      +  .166E-03,.152E-03,.124E-03 /
0160       DATA (A(I, 1,11),I=1,10)  /
0161      +  .000E+00,.000E+00,.000E+00,.238E-03,.179E-03,.145E-03,.156E-03,
0162      +  .138E-03,.129E-03,.111E-03 /
0163       DATA (A(I, 1,12),I=1,10)  /
0164      +  .000E+00,.000E+00,.000E+00,.368E-03,.400E-03,.255E-03,.262E-03,
0165      +  .221E-03,.182E-03,.112E-03 /
0166       DATA (A(I, 1,13),I=1,10)  /
0167      +  .000E+00,.000E+00,.000E+00,.000E+00,.753E-04,.712E-04,.527E-04,
0168      +  .537E-04,.538E-04,.487E-04 /
0169       DATA (A(I, 1,14),I=1,10)  /
0170      +  .000E+00,.000E+00,.000E+00,.000E+00,.103E-03,.589E-04,.578E-04,
0171      +  .468E-04,.385E-04,.269E-04 /
0172       DATA (A(I, 1,15),I=1,10)  /
0173      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.444E-04,.372E-04,
0174      +  .318E-04,.284E-04,.218E-04 /
0175       DATA (A(I, 1,16),I=1,10)  /
0176      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.487E-04,.473E-04,
0177      +  .338E-04,.243E-04,.122E-04 /
0178       DATA (A(I, 1,17),I=1,10)  /
0179      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.121E-04,.117E-04,
0180      +  .932E-05,.792E-05,.583E-05 /
0181       DATA (A(I, 1,18),I=1,10)  /
0182      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.147E-04,
0183      +  .101E-04,.756E-05,.496E-05 /
0184       DATA (A(I, 1,19),I=1,10)  /
0185      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.755E-05,
0186      +  .612E-05,.505E-05,.341E-05 /
0187       DATA (A(I, 1,20),I=1,10)  /
0188      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
0189      +  .630E-05,.444E-05,.282E-05 /
0190       DATA (A(I, 2, 1),I=1,10)  /
0191      +  .269    ,.510    ,.738    ,1.12    ,1.46    ,1.83    ,2.22    ,
0192      +  2.57    ,3.00    ,3.67     /
0193       DATA (A(I, 2, 2),I=1,10)  /
0194      +  .121    ,.133    ,.190    ,.234    ,.293    ,.332    ,.395    ,
0195      +  .431    ,.468    ,.502     /
0196       DATA (A(I, 2, 3),I=1,10)  /
0197      +  .227E-01,.374E-01,.474E-01,.578E-01,.722E-01,.794E-01,.960E-01,
0198      +  .102    ,.110    ,.120     /
0199       DATA (A(I, 2, 4),I=1,10)  /
0200      +  .287    ,.196    ,.270    ,.314    ,.373    ,.408    ,.462    ,
0201      +  .498    ,.529    ,.523     /
0202       DATA (A(I, 2, 5),I=1,10)  /
0203      +  .000E+00,.433E-01,.218E-01,.384E-01,.263E-01,.385E-01,.298E-01,
0204      +  .405E-01,.504E-01,.671E-01 /
0205       DATA (A(I, 2, 6),I=1,10)  /
0206      +  .000E+00,.151E-01,.177E-01,.159E-01,.173E-01,.173E-01,.187E-01,
0207      +  .196E-01,.201E-01,.191E-01 /
0208       DATA (A(I, 2, 7),I=1,10)  /
0209      +  .000E+00,.457E-02,.607E-02,.610E-02,.677E-02,.670E-02,.784E-02,
0210      +  .787E-02,.806E-02,.803E-02 /
0211       DATA (A(I, 2, 8),I=1,10)  /
0212      +  .000E+00,.000E+00,.702E-02,.536E-02,.558E-02,.510E-02,.554E-02,
0213      +  .546E-02,.538E-02,.489E-02 /
0214       DATA (A(I, 2, 9),I=1,10)  /
0215      +  .000E+00,.000E+00,.190E-02,.199E-02,.205E-02,.191E-02,.221E-02,
0216      +  .214E-02,.213E-02,.204E-02 /
0217       DATA (A(I, 2,10),I=1,10)  /
0218      +  .000E+00,.000E+00,.000E+00,.226E-02,.219E-02,.195E-02,.208E-02,
0219      +  .204E-02,.203E-02,.194E-02 /
0220       DATA (A(I, 2,11),I=1,10)  /
0221      +  .000E+00,.000E+00,.000E+00,.213E-02,.195E-02,.175E-02,.191E-02,
0222      +  .183E-02,.179E-02,.166E-02 /
0223       DATA (A(I, 2,12),I=1,10)  /
0224      +  .000E+00,.000E+00,.000E+00,.588E-03,.186E-02,.137E-02,.141E-02,
0225      +  .128E-02,.117E-02,.947E-03 /
0226       DATA (A(I, 2,13),I=1,10)  /
0227      +  .000E+00,.000E+00,.000E+00,.000E+00,.554E-03,.562E-03,.454E-03,
0228      +  .485E-03,.505E-03,.509E-03 /
0229       DATA (A(I, 2,14),I=1,10)  /
0230      +  .000E+00,.000E+00,.000E+00,.000E+00,.490E-03,.533E-03,.531E-03,
0231      +  .476E-03,.437E-03,.369E-03 /
0232       DATA (A(I, 2,15),I=1,10)  /
0233      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.427E-03,.382E-03,
0234      +  .358E-03,.340E-03,.294E-03 /
0235       DATA (A(I, 2,16),I=1,10)  /
0236      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.239E-03,.298E-03,
0237      +  .238E-03,.196E-03,.134E-03 /
0238       DATA (A(I, 2,17),I=1,10)  /
0239      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.299E-04,.893E-04,
0240      +  .796E-04,.744E-04,.683E-04 /
0241       DATA (A(I, 2,18),I=1,10)  /
0242      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.127E-03,
0243      +  .107E-03,.916E-04,.720E-04 /
0244       DATA (A(I, 2,19),I=1,10)  /
0245      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.397E-04,
0246      +  .630E-04,.565E-04,.461E-04 /
0247       DATA (A(I, 2,20),I=1,10)  /
0248      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
0249      +  .511E-04,.459E-04,.402E-04 /
0250       DATA (A(I, 3, 1),I=1,10)  /
0251      +  .708    ,1.02    ,1.41    ,1.91    ,2.42    ,3.00    ,3.53    ,
0252      +  4.09    ,4.71    ,5.57     /
0253       DATA (A(I, 3, 2),I=1,10)  /
0254      +  .397    ,.410    ,.539    ,.648    ,.795    ,.910    ,1.06    ,
0255      +  1.17    ,1.29    ,1.42     /
0256       DATA (A(I, 3, 3),I=1,10)  /
0257      +  .845E-01,.122    ,.157    ,.190    ,.232    ,.262    ,.307    ,
0258      +  .335    ,.366    ,.402     /
0259       DATA (A(I, 3, 4),I=1,10)  /
0260      +  .210    ,.379    ,.450    ,.490    ,.574    ,.636    ,.709    ,
0261      +  .769    ,.820    ,.849     /
0262       DATA (A(I, 3, 5),I=1,10)  /
0263      +  .000E+00,.102    ,.675E-01,.104    ,.858E-01,.115    ,.102    ,
0264      +  .129    ,.154    ,.194     /
0265       DATA (A(I, 3, 6),I=1,10)  /
0266      +  .000E+00,.392E-01,.615E-01,.593E-01,.649E-01,.674E-01,.735E-01,
0267      +  .779E-01,.817E-01,.828E-01 /
0268       DATA (A(I, 3, 7),I=1,10)  /
0269      +  .000E+00,.539E-02,.222E-01,.238E-01,.269E-01,.280E-01,.320E-01,
0270      +  .334E-01,.350E-01,.361E-01 /
0271       DATA (A(I, 3, 8),I=1,10)  /
0272      +  .000E+00,.000E+00,.838E-02,.130E-01,.133E-01,.131E-01,.141E-01,
0273      +  .144E-01,.149E-01,.152E-01 /
0274       DATA (A(I, 3, 9),I=1,10)  /
0275      +  .000E+00,.000E+00,.228E-02,.647E-02,.688E-02,.687E-02,.772E-02,
0276      +  .786E-02,.811E-02,.824E-02 /
0277       DATA (A(I, 3,10),I=1,10)  /
0278      +  .000E+00,.000E+00,.000E+00,.664E-02,.828E-02,.802E-02,.845E-02,
0279      +  .869E-02,.902E-02,.930E-02 /
0280       DATA (A(I, 3,11),I=1,10)  /
0281      +  .000E+00,.000E+00,.000E+00,.338E-02,.735E-02,.710E-02,.767E-02,
0282      +  .767E-02,.776E-02,.756E-02 /
0283       DATA (A(I, 3,12),I=1,10)  /
0284      +  .000E+00,.000E+00,.000E+00,.280E-03,.262E-02,.349E-02,.342E-02,
0285      +  .322E-02,.312E-02,.291E-02 /
0286       DATA (A(I, 3,13),I=1,10)  /
0287      +  .000E+00,.000E+00,.000E+00,.000E+00,.618E-03,.161E-02,.138E-02,
0288      +  .148E-02,.155E-02,.166E-02 /
0289       DATA (A(I, 3,14),I=1,10)  /
0290      +  .000E+00,.000E+00,.000E+00,.000E+00,.313E-03,.128E-02,.161E-02,
0291      +  .150E-02,.144E-02,.134E-02 /
0292       DATA (A(I, 3,15),I=1,10)  /
0293      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.645E-03,.118E-02,
0294      +  .115E-02,.111E-02,.103E-02 /
0295       DATA (A(I, 3,16),I=1,10)  /
0296      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.117E-03,.497E-03,
0297      +  .581E-03,.501E-03,.401E-03 /
0298       DATA (A(I, 3,17),I=1,10)  /
0299      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.115E-04,.997E-04,
0300      +  .202E-03,.203E-03,.206E-03 /
0301       DATA (A(I, 3,18),I=1,10)  /
0302      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.877E-04,
0303      +  .242E-03,.263E-03,.226E-03 /
0304       DATA (A(I, 3,19),I=1,10)  /
0305      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.158E-04,
0306      +  .881E-04,.152E-03,.136E-03 /
0307       DATA (A(I, 3,20),I=1,10)  /
0308      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
0309      +  .358E-04,.997E-04,.117E-03 /
0310       DATA (A(I, 4, 1),I=1,10)  /
0311      +  .945    ,1.29    ,1.40    ,1.98    ,2.73    ,3.17    ,3.77    ,
0312      +  4.29    ,4.78    ,5.54     /
0313       DATA (A(I, 4, 2),I=1,10)  /
0314      +  .581    ,.599    ,.645    ,.839    ,1.10    ,1.25    ,1.47    ,
0315      +  1.64    ,1.78    ,1.99     /
0316       DATA (A(I, 4, 3),I=1,10)  /
0317      +  .127    ,.182    ,.202    ,.264    ,.344    ,.387    ,.455    ,
0318      +  .504    ,.549    ,.611     /
0319       DATA (A(I, 4, 4),I=1,10)  /
0320      +  .183    ,.464    ,.351    ,.444    ,.642    ,.659    ,.772    ,
0321      +  .830    ,.882    ,.930     /
0322       DATA (A(I, 4, 5),I=1,10)  /
0323      +  .000E+00,.122    ,.803E-01,.136    ,.134    ,.173    ,.164    ,
0324      +  .203    ,.239    ,.300     /
0325       DATA (A(I, 4, 6),I=1,10)  /
0326      +  .000E+00,.393E-01,.766E-01,.872E-01,.108    ,.111    ,.123    ,
0327      +  .132    ,.139    ,.145     /
0328       DATA (A(I, 4, 7),I=1,10)  /
0329      +  .000E+00,.416E-02,.289E-01,.360E-01,.454E-01,.477E-01,.549E-01,
0330      +  .583E-01,.618E-01,.654E-01 /
0331       DATA (A(I, 4, 8),I=1,10)  /
0332      +  .000E+00,.000E+00,.761E-02,.157E-01,.214E-01,.205E-01,.233E-01,
0333      +  .241E-01,.255E-01,.271E-01 /
0334       DATA (A(I, 4, 9),I=1,10)  /
0335      +  .000E+00,.000E+00,.238E-02,.803E-02,.123E-01,.123E-01,.140E-01,
0336      +  .145E-01,.153E-01,.160E-01 /
0337       DATA (A(I, 4,10),I=1,10)  /
0338      +  .000E+00,.000E+00,.000E+00,.695E-02,.150E-01,.154E-01,.166E-01,
0339      +  .172E-01,.181E-01,.192E-01 /
0340       DATA (A(I, 4,11),I=1,10)  /
0341      +  .000E+00,.000E+00,.000E+00,.355E-02,.104E-01,.143E-01,.156E-01,
0342      +  .158E-01,.164E-01,.165E-01 /
0343       DATA (A(I, 4,12),I=1,10)  /
0344      +  .000E+00,.000E+00,.000E+00,.112E-03,.276E-02,.568E-02,.736E-02,
0345      +  .684E-02,.691E-02,.661E-02 /
0346       DATA (A(I, 4,13),I=1,10)  /
0347      +  .000E+00,.000E+00,.000E+00,.000E+00,.740E-03,.222E-02,.339E-02,
0348      +  .352E-02,.382E-02,.409E-02 /
0349       DATA (A(I, 4,14),I=1,10)  /
0350      +  .000E+00,.000E+00,.000E+00,.000E+00,.369E-03,.160E-02,.322E-02,
0351      +  .375E-02,.375E-02,.355E-02 /
0352       DATA (A(I, 4,15),I=1,10)  /
0353      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.750E-03,.190E-02,
0354      +  .298E-02,.319E-02,.299E-02 /
0355       DATA (A(I, 4,16),I=1,10)  /
0356      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.260E-03,.673E-03,
0357      +  .117E-02,.156E-02,.126E-02 /
0358       DATA (A(I, 4,17),I=1,10)  /
0359      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.283E-05,.131E-03,
0360      +  .363E-03,.618E-03,.690E-03 /
0361       DATA (A(I, 4,18),I=1,10)  /
0362      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.205E-03,
0363      +  .378E-03,.709E-03,.844E-03 /
0364       DATA (A(I, 4,19),I=1,10)  /
0365      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.654E-05,
0366      +  .150E-03,.341E-03,.527E-03 /
0367       DATA (A(I, 4,20),I=1,10)  /
0368      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
0369      +  .957E-04,.197E-03,.406E-03 /
0370       DATA (A(I, 5, 1),I=1,10)  /
0371      +  1.16    ,1.70    ,2.19    ,2.79    ,3.33    ,3.90    ,4.49    ,
0372      +  5.07    ,5.66    ,6.38     /
0373       DATA (A(I, 5, 2),I=1,10)  /
0374      +  .779    ,.899    ,1.09    ,1.28    ,1.51    ,1.71    ,1.96    ,
0375      +  2.18    ,2.39    ,2.62     /
0376       DATA (A(I, 5, 3),I=1,10)  /
0377      +  .167    ,.263    ,.334    ,.408    ,.482    ,.548    ,.632    ,
0378      +  .700    ,.767    ,.840     /
0379       DATA (A(I, 5, 4),I=1,10)  /
0380      +  .203    ,.565    ,.845    ,.867    ,.906    ,.961    ,1.08    ,
0381      +  1.13    ,1.21    ,1.25     /
0382       DATA (A(I, 5, 5),I=1,10)  /
0383      +  .000E+00,.129    ,.152    ,.237    ,.208    ,.268    ,.258    ,
0384      +  .312    ,.368    ,.450     /
0385       DATA (A(I, 5, 6),I=1,10)  /
0386      +  .000E+00,.460E-01,.126    ,.174    ,.182    ,.188    ,.208    ,
0387      +  .219    ,.233    ,.239     /
0388       DATA (A(I, 5, 7),I=1,10)  /
0389      +  .000E+00,.289E-02,.380E-01,.611E-01,.788E-01,.845E-01,.974E-01,
0390      +  .103    ,.111    ,.117     /
0391       DATA (A(I, 5, 8),I=1,10)  /
0392      +  .000E+00,.000E+00,.137E-01,.223E-01,.374E-01,.436E-01,.488E-01,
0393      +  .488E-01,.524E-01,.547E-01 /
0394       DATA (A(I, 5, 9),I=1,10)  /
0395      +  .000E+00,.000E+00,.162E-02,.114E-01,.198E-01,.263E-01,.315E-01,
0396      +  .323E-01,.348E-01,.364E-01 /
0397       DATA (A(I, 5,10),I=1,10)  /
0398      +  .000E+00,.000E+00,.000E+00,.149E-01,.240E-01,.320E-01,.428E-01,
0399      +  .436E-01,.469E-01,.493E-01 /
0400       DATA (A(I, 5,11),I=1,10)  /
0401      +  .000E+00,.000E+00,.000E+00,.562E-02,.194E-01,.290E-01,.408E-01,
0402      +  .460E-01,.492E-01,.500E-01 /
0403       DATA (A(I, 5,12),I=1,10)  /
0404      +  .000E+00,.000E+00,.000E+00,.476E-04,.106E-01,.134E-01,.191E-01,
0405      +  .227E-01,.264E-01,.253E-01 /
0406       DATA (A(I, 5,13),I=1,10)  /
0407      +  .000E+00,.000E+00,.000E+00,.000E+00,.281E-02,.679E-02,.879E-02,
0408      +  .123E-01,.165E-01,.190E-01 /
0409       DATA (A(I, 5,14),I=1,10)  /
0410      +  .000E+00,.000E+00,.000E+00,.000E+00,.542E-04,.847E-02,.125E-01,
0411      +  .144E-01,.173E-01,.192E-01 /
0412       DATA (A(I, 5,15),I=1,10)  /
0413      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.652E-02,.982E-02,
0414      +  .129E-01,.159E-01,.192E-01 /
0415       DATA (A(I, 5,16),I=1,10)  /
0416      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.109E-03,.688E-02,
0417      +  .751E-02,.845E-02,.905E-02 /
0418       DATA (A(I, 5,17),I=1,10)  /
0419      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.823E-06,.237E-02,
0420      +  .318E-02,.446E-02,.569E-02 /
0421       DATA (A(I, 5,18),I=1,10)  /
0422      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.604E-03,
0423      +  .610E-02,.673E-02,.827E-02 /
0424       DATA (A(I, 5,19),I=1,10)  /
0425      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.716E-06,
0426      +  .412E-02,.519E-02,.617E-02 /
0427       DATA (A(I, 5,20),I=1,10)  /
0428      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
0429      +  .710E-03,.543E-02,.674E-02 /
0430       DATA (A(I, 6, 1),I=1,10)  /
0431      +  1.36    ,2.08    ,2.67    ,3.30    ,3.94    ,4.62    ,5.18    ,
0432      +  3.60    ,3.64    ,3.95     /
0433       DATA (A(I, 6, 2),I=1,10)  /
0434      +  1.07    ,1.33    ,1.58    ,1.82    ,2.10    ,2.44    ,2.74    ,
0435      +  1.78    ,1.73    ,1.80     /
0436       DATA (A(I, 6, 3),I=1,10)  /
0437      +  .158    ,.276    ,.402    ,.506    ,.609    ,.700    ,.802    ,
0438      +  .638    ,.629    ,.658     /
0439       DATA (A(I, 6, 4),I=1,10)  /
0440      +  .308    ,.739    ,1.02    ,1.12    ,1.26    ,1.35    ,1.57    ,
0441      +  1.94    ,1.71    ,1.55     /
0442       DATA (A(I, 6, 5),I=1,10)  /
0443      +  .000E+00,.217    ,.183    ,.324    ,.276    ,.395    ,.393    ,
0444      +  .558    ,.602    ,.681     /
0445       DATA (A(I, 6, 6),I=1,10)  /
0446      +  .000E+00,.658E-01,.251    ,.267    ,.299    ,.326    ,.386    ,
0447      +  .452    ,.475    ,.409     /
0448       DATA (A(I, 6, 7),I=1,10)  /
0449      +  .000E+00,.198E-02,.774E-01,.136    ,.149    ,.164    ,.187    ,
0450      +  .210    ,.238    ,.256     /
0451       DATA (A(I, 6, 8),I=1,10)  /
0452      +  .000E+00,.000E+00,.290E-01,.122    ,.139    ,.128    ,.129    ,
0453      +  .137    ,.147    ,.167     /
0454       DATA (A(I, 6, 9),I=1,10)  /
0455      +  .000E+00,.000E+00,.699E-03,.617E-01,.750E-01,.801E-01,.905E-01,
0456      +  .974E-01,.105    ,.122     /
0457       DATA (A(I, 6,10),I=1,10)  /
0458      +  .000E+00,.000E+00,.000E+00,.310E-01,.112    ,.127    ,.140    ,
0459      +  .143    ,.155    ,.176     /
0460       DATA (A(I, 6,11),I=1,10)  /
0461      +  .000E+00,.000E+00,.000E+00,.277E-02,.889E-01,.143    ,.150    ,
0462      +  .175    ,.184    ,.208     /
0463       DATA (A(I, 6,12),I=1,10)  /
0464      +  .000E+00,.000E+00,.000E+00,.202E-04,.343E-01,.959E-01,.109    ,
0465      +  .115    ,.112    ,.116     /
0466       DATA (A(I, 6,13),I=1,10)  /
0467      +  .000E+00,.000E+00,.000E+00,.000E+00,.186E-02,.435E-01,.512E-01,
0468      +  .744E-01,.856E-01,.103     /
0469       DATA (A(I, 6,14),I=1,10)  /
0470      +  .000E+00,.000E+00,.000E+00,.000E+00,.144E-04,.427E-01,.786E-01,
0471      +  .911E-01,.993E-01,.108     /
0472       DATA (A(I, 6,15),I=1,10)  /
0473      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.466E-02,.518E-01,
0474      +  .848E-01,.109    ,.119     /
0475       DATA (A(I, 6,16),I=1,10)  /
0476      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.655E-05,.330E-01,
0477      +  .586E-01,.617E-01,.594E-01 /
0478       DATA (A(I, 6,17),I=1,10)  /
0479      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.228E-06,.328E-02,
0480      +  .190E-01,.301E-01,.454E-01 /
0481       DATA (A(I, 6,18),I=1,10)  /
0482      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.218E-04,
0483      +  .272E-01,.501E-01,.707E-01 /
0484       DATA (A(I, 6,19),I=1,10)  /
0485      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.146E-06,
0486      +  .441E-02,.378E-01,.556E-01 /
0487       DATA (A(I, 6,20),I=1,10)  /
0488      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
0489      +  .160E-03,.204E-01,.679E-01 /
0490       DATA (A(I, 7, 1),I=1,10)  /
0491      +  .522    ,.862    ,1.14    ,1.40    ,1.70    ,1.94    ,2.26    ,
0492      +  2.48    ,2.72    ,3.95     /
0493       DATA (A(I, 7, 2),I=1,10)  /
0494      +  .314    ,.450    ,.588    ,.692    ,.834    ,.936    ,1.09    ,
0495      +  1.18    ,1.28    ,1.80     /
0496       DATA (A(I, 7, 3),I=1,10)  /
0497      +  .814E-01,.147    ,.189    ,.226    ,.272    ,.302    ,.351    ,
0498      +  .378    ,.406    ,.658     /
0499       DATA (A(I, 7, 4),I=1,10)  /
0500      +  .252    ,.864    ,1.01    ,.851    ,.837    ,.774    ,.763    ,
0501      +  .757    ,.748    ,1.55     /
0502       DATA (A(I, 7, 5),I=1,10)  /
0503      +  .000E+00,.225    ,.180    ,.276    ,.193    ,.240    ,.190    ,
0504      +  .228    ,.259    ,.681     /
0505       DATA (A(I, 7, 6),I=1,10)  /
0506      +  .000E+00,.485E-01,.272    ,.273    ,.253    ,.216    ,.206    ,
0507      +  .197    ,.191    ,.409     /
0508       DATA (A(I, 7, 7),I=1,10)  /
0509      +  .000E+00,.137E-02,.752E-01,.137    ,.152    ,.134    ,.125    ,
0510      +  .119    ,.116    ,.256     /
0511       DATA (A(I, 7, 8),I=1,10)  /
0512      +  .000E+00,.000E+00,.220E-01,.155    ,.175    ,.155    ,.116    ,
0513      +  .977E-01,.858E-01,.167     /
0514       DATA (A(I, 7, 9),I=1,10)  /
0515      +  .000E+00,.000E+00,.326E-03,.695E-01,.881E-01,.106    ,.897E-01,
0516      +  .782E-01,.706E-01,.122     /
0517       DATA (A(I, 7,10),I=1,10)  /
0518      +  .000E+00,.000E+00,.000E+00,.261E-01,.124    ,.131    ,.156    ,
0519      +  .141    ,.121    ,.176     /
0520       DATA (A(I, 7,11),I=1,10)  /
0521      +  .000E+00,.000E+00,.000E+00,.785E-03,.864E-01,.130    ,.170    ,
0522      +  .182    ,.172    ,.208     /
0523       DATA (A(I, 7,12),I=1,10)  /
0524      +  .000E+00,.000E+00,.000E+00,.896E-05,.225E-01,.105    ,.126    ,
0525      +  .126    ,.135    ,.116     /
0526       DATA (A(I, 7,13),I=1,10)  /
0527      +  .000E+00,.000E+00,.000E+00,.000E+00,.542E-03,.427E-01,.553E-01,
0528      +  .744E-01,.980E-01,.103     /
0529       DATA (A(I, 7,14),I=1,10)  /
0530      +  .000E+00,.000E+00,.000E+00,.000E+00,.515E-05,.377E-01,.831E-01,
0531      +  .985E-01,.104    ,.108     /
0532       DATA (A(I, 7,15),I=1,10)  /
0533      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.285E-02,.495E-01,
0534      +  .871E-01,.106    ,.119     /
0535       DATA (A(I, 7,16),I=1,10)  /
0536      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.110E-05,.284E-01,
0537      +  .588E-01,.657E-01,.594E-01 /
0538       DATA (A(I, 7,17),I=1,10)  /
0539      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.722E-07,.176E-02,
0540      +  .170E-01,.305E-01,.454E-01 /
0541       DATA (A(I, 7,18),I=1,10)  /
0542      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.148E-05,
0543      +  .213E-01,.492E-01,.707E-01 /
0544       DATA (A(I, 7,19),I=1,10)  /
0545      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.323E-07,
0546      +  .722E-02,.359E-01,.556E-01 /
0547       DATA (A(I, 7,20),I=1,10)  /
0548      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
0549      +  .461E-05,.155E-01,.679E-01 /
0550       DATA (A(I, 8, 1),I=1,10)  /
0551      +  .630    ,.974    ,1.29    ,1.58    ,1.89    ,2.16    ,2.49    ,
0552      +  2.75    ,3.02    ,3.95     /
0553       DATA (A(I, 8, 2),I=1,10)  /
0554      +  .328    ,.459    ,.613    ,.735    ,.879    ,.994    ,1.15    ,
0555      +  1.27    ,1.38    ,1.80     /
0556       DATA (A(I, 8, 3),I=1,10)  /
0557      +  .748E-01,.121    ,.164    ,.197    ,.235    ,.265    ,.310    ,
0558      +  .339    ,.370    ,.658     /
0559       DATA (A(I, 8, 4),I=1,10)  /
0560      +  .194    ,.211    ,.337    ,.344    ,.339    ,.351    ,.390    ,
0561      +  .419    ,.442    ,1.55     /
0562       DATA (A(I, 8, 5),I=1,10)  /
0563      +  .000E+00,.869E-01,.725E-01,.113    ,.810E-01,.106    ,.951E-01,
0564      +  .120    ,.143    ,.681     /
0565       DATA (A(I, 8, 6),I=1,10)  /
0566      +  .000E+00,.288E-01,.102    ,.922E-01,.857E-01,.845E-01,.932E-01,
0567      +  .983E-01,.102    ,.409     /
0568       DATA (A(I, 8, 7),I=1,10)  /
0569      +  .000E+00,.668E-03,.533E-01,.575E-01,.493E-01,.482E-01,.539E-01,
0570      +  .558E-01,.582E-01,.256     /
0571       DATA (A(I, 8, 8),I=1,10)  /
0572      +  .000E+00,.000E+00,.205E-01,.808E-01,.510E-01,.409E-01,.406E-01,
0573      +  .394E-01,.389E-01,.167     /
0574       DATA (A(I, 8, 9),I=1,10)  /
0575      +  .000E+00,.000E+00,.999E-04,.647E-01,.385E-01,.325E-01,.325E-01,
0576      +  .316E-01,.314E-01,.122     /
0577       DATA (A(I, 8,10),I=1,10)  /
0578      +  .000E+00,.000E+00,.000E+00,.169E-01,.834E-01,.611E-01,.565E-01,
0579      +  .533E-01,.519E-01,.176     /
0580       DATA (A(I, 8,11),I=1,10)  /
0581      +  .000E+00,.000E+00,.000E+00,.107E-03,.769E-01,.922E-01,.805E-01,
0582      +  .745E-01,.711E-01,.208     /
0583       DATA (A(I, 8,12),I=1,10)  /
0584      +  .000E+00,.000E+00,.000E+00,.180E-05,.143E-01,.983E-01,.775E-01,
0585      +  .627E-01,.541E-01,.116     /
0586       DATA (A(I, 8,13),I=1,10)  /
0587      +  .000E+00,.000E+00,.000E+00,.000E+00,.157E-04,.346E-01,.507E-01,
0588      +  .479E-01,.455E-01,.103     /
0589       DATA (A(I, 8,14),I=1,10)  /
0590      +  .000E+00,.000E+00,.000E+00,.000E+00,.752E-06,.248E-01,.721E-01,
0591      +  .728E-01,.611E-01,.108     /
0592       DATA (A(I, 8,15),I=1,10)  /
0593      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.686E-04,.356E-01,
0594      +  .731E-01,.791E-01,.119     /
0595       DATA (A(I, 8,16),I=1,10)  /
0596      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.838E-07,.151E-01,
0597      +  .470E-01,.567E-01,.594E-01 /
0598       DATA (A(I, 8,17),I=1,10)  /
0599      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.759E-08,.400E-04,
0600      +  .193E-01,.313E-01,.454E-01 /
0601       DATA (A(I, 8,18),I=1,10)  /
0602      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.385E-07,
0603      +  .921E-02,.353E-01,.707E-01 /
0604       DATA (A(I, 8,19),I=1,10)  /
0605      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.219E-08,
0606      +  .348E-03,.226E-01,.556E-01 /
0607       DATA (A(I, 8,20),I=1,10)  /
0608      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
0609      +  .212E-07,.149E-01,.679E-01 /
0610       DATA (A(I, 9, 1),I=1,10)  /
0611      +  .736    ,1.13    ,1.49    ,1.82    ,2.20    ,2.49    ,2.86    ,
0612      +  3.17    ,3.49    ,3.95     /
0613       DATA (A(I, 9, 2),I=1,10)  /
0614      +  .339    ,.492    ,.658    ,.789    ,.958    ,1.08    ,1.25    ,
0615      +  1.37    ,1.50    ,1.80     /
0616       DATA (A(I, 9, 3),I=1,10)  /
0617      +  .680E-01,.110    ,.150    ,.180    ,.222    ,.247    ,.289    ,
0618      +  .318    ,.349    ,.658     /
0619       DATA (A(I, 9, 4),I=1,10)  /
0620      +  .110    ,.104    ,.157    ,.156    ,.210    ,.205    ,.246    ,
0621      +  .274    ,.300    ,1.55     /
0622       DATA (A(I, 9, 5),I=1,10)  /
0623      +  .000E+00,.379E-01,.347E-01,.477E-01,.486E-01,.576E-01,.569E-01,
0624      +  .732E-01,.893E-01,.681     /
0625       DATA (A(I, 9, 6),I=1,10)  /
0626      +  .000E+00,.223E-01,.354E-01,.312E-01,.436E-01,.400E-01,.489E-01,
0627      +  .548E-01,.600E-01,.409     /
0628       DATA (A(I, 9, 7),I=1,10)  /
0629      +  .000E+00,.338E-03,.149E-01,.142E-01,.215E-01,.188E-01,.248E-01,
0630      +  .278E-01,.307E-01,.256     /
0631       DATA (A(I, 9, 8),I=1,10)  /
0632      +  .000E+00,.000E+00,.553E-02,.862E-02,.150E-01,.106E-01,.145E-01,
0633      +  .165E-01,.181E-01,.167     /
0634       DATA (A(I, 9, 9),I=1,10)  /
0635      +  .000E+00,.000E+00,.375E-04,.641E-02,.111E-01,.792E-02,.112E-01,
0636      +  .127E-01,.140E-01,.122     /
0637       DATA (A(I, 9,10),I=1,10)  /
0638      +  .000E+00,.000E+00,.000E+00,.112E-01,.200E-01,.127E-01,.176E-01,
0639      +  .200E-01,.220E-01,.176     /
0640       DATA (A(I, 9,11),I=1,10)  /
0641      +  .000E+00,.000E+00,.000E+00,.244E-04,.261E-01,.162E-01,.232E-01,
0642      +  .263E-01,.287E-01,.208     /
0643       DATA (A(I, 9,12),I=1,10)  /
0644      +  .000E+00,.000E+00,.000E+00,.455E-06,.635E-02,.121E-01,.186E-01,
0645      +  .201E-01,.207E-01,.116     /
0646       DATA (A(I, 9,13),I=1,10)  /
0647      +  .000E+00,.000E+00,.000E+00,.000E+00,.146E-05,.922E-02,.116E-01,
0648      +  .145E-01,.165E-01,.103     /
0649       DATA (A(I, 9,14),I=1,10)  /
0650      +  .000E+00,.000E+00,.000E+00,.000E+00,.135E-06,.128E-01,.202E-01,
0651      +  .215E-01,.220E-01,.108     /
0652       DATA (A(I, 9,15),I=1,10)  /
0653      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.237E-05,.229E-01,
0654      +  .259E-01,.271E-01,.119     /
0655       DATA (A(I, 9,16),I=1,10)  /
0656      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.100E-07,.534E-02,
0657      +  .210E-01,.193E-01,.594E-01 /
0658       DATA (A(I, 9,17),I=1,10)  /
0659      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.915E-09,.847E-06,
0660      +  .119E-01,.125E-01,.454E-01 /
0661       DATA (A(I, 9,18),I=1,10)  /
0662      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.298E-08,
0663      +  .101E-01,.242E-01,.707E-01 /
0664       DATA (A(I, 9,19),I=1,10)  /
0665      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.196E-09,
0666      +  .243E-05,.234E-01,.556E-01 /
0667       DATA (A(I, 9,20),I=1,10)  /
0668      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
0669      +  .575E-09,.364E-02,.679E-01 /
0670       DATA (A(I,10, 1),I=1,10)  /
0671      +  .959    ,1.46    ,1.92    ,2.34    ,2.80    ,3.24    ,3.64    ,
0672      +  4.05    ,4.48    ,3.95     /
0673       DATA (A(I,10, 2),I=1,10)  /
0674      +  .343    ,.516    ,.692    ,.836    ,1.01    ,1.16    ,1.31    ,
0675      +  1.46    ,1.61    ,1.80     /
0676       DATA (A(I,10, 3),I=1,10)  /
0677      +  .512E-01,.837E-01,.115    ,.138    ,.169    ,.195    ,.220    ,
0678      +  .245    ,.270    ,.658     /
0679       DATA (A(I,10, 4),I=1,10)  /
0680      +  .274E-01,.361E-01,.510E-01,.562E-01,.703E-01,.828E-01,.877E-01,
0681      +  .996E-01,.111    ,1.55     /
0682       DATA (A(I,10, 5),I=1,10)  /
0683      +  .000E+00,.850E-02,.875E-02,.118E-01,.124E-01,.170E-01,.154E-01,
0684      +  .194E-01,.237E-01,.681     /
0685       DATA (A(I,10, 6),I=1,10)  /
0686      +  .000E+00,.345E-02,.519E-02,.533E-02,.691E-02,.842E-02,.844E-02,
0687      +  .987E-02,.113E-01,.409     /
0688       DATA (A(I,10, 7),I=1,10)  /
0689      +  .000E+00,.722E-04,.130E-02,.135E-02,.189E-02,.240E-02,.235E-02,
0690      +  .281E-02,.331E-02,.256     /
0691       DATA (A(I,10, 8),I=1,10)  /
0692      +  .000E+00,.000E+00,.283E-03,.272E-03,.394E-03,.557E-03,.480E-03,
0693      +  .616E-03,.775E-03,.167     /
0694       DATA (A(I,10, 9),I=1,10)  /
0695      +  .000E+00,.000E+00,.457E-05,.122E-03,.192E-03,.275E-03,.225E-03,
0696      +  .292E-03,.373E-03,.122     /
0697       DATA (A(I,10,10),I=1,10)  /
0698      +  .000E+00,.000E+00,.000E+00,.119E-03,.185E-03,.278E-03,.201E-03,
0699      +  .274E-03,.364E-03,.176     /
0700       DATA (A(I,10,11),I=1,10)  /
0701      +  .000E+00,.000E+00,.000E+00,.140E-05,.129E-03,.200E-03,.137E-03,
0702      +  .188E-03,.252E-03,.208     /
0703       DATA (A(I,10,12),I=1,10)  /
0704      +  .000E+00,.000E+00,.000E+00,.207E-07,.307E-04,.518E-04,.278E-04,
0705      +  .421E-04,.608E-04,.116     /
0706       DATA (A(I,10,13),I=1,10)  /
0707      +  .000E+00,.000E+00,.000E+00,.000E+00,.306E-07,.252E-04,.111E-04,
0708      +  .188E-04,.295E-04,.103     /
0709       DATA (A(I,10,14),I=1,10)  /
0710      +  .000E+00,.000E+00,.000E+00,.000E+00,.321E-08,.220E-04,.104E-04,
0711      +  .162E-04,.243E-04,.108     /
0712       DATA (A(I,10,15),I=1,10)  /
0713      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.770E-08,.632E-05,
0714      +  .105E-04,.162E-04,.119     /
0715       DATA (A(I,10,16),I=1,10)  /
0716      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.117E-09,.199E-05,
0717      +  .321E-05,.492E-05,.594E-01 /
0718       DATA (A(I,10,17),I=1,10)  /
0719      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.888E-11,.323E-09,
0720      +  .106E-05,.192E-05,.454E-01 /
0721       DATA (A(I,10,18),I=1,10)  /
0722      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.174E-10,
0723      +  .131E-05,.218E-05,.707E-01 /
0724       DATA (A(I,10,19),I=1,10)  /
0725      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.994E-12,
0726      +  .233E-09,.104E-05,.556E-01 /
0727       DATA (A(I,10,20),I=1,10)  /
0728      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
0729      +  .144E-11,.724E-06,.679E-01 /
0730       DATA (AE(I, 1, 1),I=1,10)  /
0731      +  7.27    ,6.29    ,7.76    ,6.70    ,8.17    ,7.34    ,8.70    ,
0732      +  8.02    ,7.37    ,6.18     /
0733       DATA (AE(I, 1, 2),I=1,10)  /
0734      +  7.41    ,7.52    ,8.14    ,8.20    ,8.96    ,9.05    ,9.96    ,
0735      +  10.0    ,10.1    ,9.86     /
0736       DATA (AE(I, 1, 3),I=1,10)  /
0737      +  7.72    ,7.69    ,9.17    ,8.99    ,10.6    ,10.5    ,12.1    ,
0738      +  12.1    ,12.0    ,11.5     /
0739       DATA (AE(I, 1, 4),I=1,10)  /
0740      +  7.90    ,8.48    ,9.50    ,9.94    ,10.8    ,11.4    ,12.2    ,
0741      +  12.8    ,13.3    ,13.8     /
0742       DATA (AE(I, 1, 5),I=1,10)  /
0743      +  .000E+00,8.52    ,9.59    ,10.1    ,11.1    ,11.8    ,12.7    ,
0744      +  13.3    ,13.8    ,14.4     /
0745       DATA (AE(I, 1, 6),I=1,10)  /
0746      +  .000E+00,9.00    ,10.7    ,11.7    ,13.2    ,14.2    ,15.6    ,
0747      +  16.5    ,17.3    ,18.0     /
0748       DATA (AE(I, 1, 7),I=1,10)  /
0749      +  .000E+00,9.01    ,11.1    ,11.9    ,14.3    ,15.0    ,17.4    ,
0750      +  18.0    ,18.6    ,18.8     /
0751       DATA (AE(I, 1, 8),I=1,10)  /
0752      +  .000E+00,.000E+00,11.2    ,12.4    ,14.5    ,15.7    ,17.6    ,
0753      +  18.8    ,19.9    ,20.9     /
0754       DATA (AE(I, 1, 9),I=1,10)  /
0755      +  .000E+00,.000E+00,11.4    ,12.7    ,15.5    ,16.6    ,19.3    ,
0756      +  20.2    ,21.1    ,21.7     /
0757       DATA (AE(I, 1,10),I=1,10)  /
0758      +  .000E+00,.000E+00,.000E+00,13.2    ,15.8    ,17.3    ,19.9    ,
0759      +  21.2    ,22.4    ,23.2     /
0760       DATA (AE(I, 1,11),I=1,10)  /
0761      +  .000E+00,.000E+00,.000E+00,13.2    ,16.3    ,17.8    ,20.8    ,
0762      +  22.1    ,23.3    ,24.2     /
0763       DATA (AE(I, 1,12),I=1,10)  /
0764      +  .000E+00,.000E+00,.000E+00,13.4    ,16.2    ,18.2    ,21.0    ,
0765      +  22.8    ,24.4    ,25.9     /
0766       DATA (AE(I, 1,13),I=1,10)  /
0767      +  .000E+00,.000E+00,.000E+00,.000E+00,16.5    ,18.4    ,21.6    ,
0768      +  23.2    ,24.8    ,26.2     /
0769       DATA (AE(I, 1,14),I=1,10)  /
0770      +  .000E+00,.000E+00,.000E+00,.000E+00,16.7    ,19.0    ,22.3    ,
0771      +  24.3    ,26.1    ,27.4     /
0772       DATA (AE(I, 1,15),I=1,10)  /
0773      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,19.1    ,22.8    ,
0774      +  24.7    ,26.6    ,28.2     /
0775       DATA (AE(I, 1,16),I=1,10)  /
0776      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,19.2    ,23.0    ,
0777      +  25.3    ,27.5    ,29.5     /
0778       DATA (AE(I, 1,17),I=1,10)  /
0779      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,19.6    ,23.3    ,
0780      +  25.6    ,27.8    ,29.6     /
0781       DATA (AE(I, 1,18),I=1,10)  /
0782      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,23.6    ,
0783      +  26.2    ,28.5    ,30.4     /
0784       DATA (AE(I, 1,19),I=1,10)  /
0785      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,23.7    ,
0786      +  26.3    ,28.8    ,31.0     /
0787       DATA (AE(I, 1,20),I=1,10)  /
0788      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
0789      +  26.5    ,29.2    ,31.5     /
0790       DATA (AE(I, 2, 1),I=1,10)  /
0791      +  8.74    ,8.16    ,9.25    ,8.45    ,9.46    ,8.90    ,9.83    ,
0792      +  9.38    ,8.96    ,8.15     /
0793       DATA (AE(I, 2, 2),I=1,10)  /
0794      +  8.96    ,9.30    ,9.95    ,10.0    ,10.8    ,10.9    ,11.7    ,
0795      +  11.8    ,11.9    ,11.8     /
0796       DATA (AE(I, 2, 3),I=1,10)  /
0797      +  9.44    ,9.66    ,11.0    ,11.0    ,12.3    ,12.5    ,13.7    ,
0798      +  13.9    ,14.0    ,13.8     /
0799       DATA (AE(I, 2, 4),I=1,10)  /
0800      +  8.86    ,9.81    ,10.8    ,11.2    ,12.0    ,12.6    ,13.4    ,
0801      +  14.0    ,14.5    ,15.1     /
0802       DATA (AE(I, 2, 5),I=1,10)  /
0803      +  .000E+00,10.2    ,11.4    ,12.0    ,12.9    ,13.6    ,14.5    ,
0804      +  15.1    ,15.7    ,16.3     /
0805       DATA (AE(I, 2, 6),I=1,10)  /
0806      +  .000E+00,10.7    ,12.5    ,13.5    ,15.1    ,16.0    ,17.5    ,
0807      +  18.3    ,19.2    ,19.9     /
0808       DATA (AE(I, 2, 7),I=1,10)  /
0809      +  .000E+00,11.5    ,12.9    ,13.9    ,16.1    ,17.0    ,19.1    ,
0810      +  19.8    ,20.6    ,21.0     /
0811       DATA (AE(I, 2, 8),I=1,10)  /
0812      +  .000E+00,.000E+00,12.4    ,13.8    ,15.9    ,17.2    ,19.1    ,
0813      +  20.3    ,21.4    ,22.3     /
0814       DATA (AE(I, 2, 9),I=1,10)  /
0815      +  .000E+00,.000E+00,13.4    ,14.5    ,17.1    ,18.3    ,20.9    ,
0816      +  21.9    ,23.0    ,23.7     /
0817       DATA (AE(I, 2,10),I=1,10)  /
0818      +  .000E+00,.000E+00,.000E+00,14.9    ,17.5    ,19.1    ,21.6    ,
0819      +  22.9    ,24.1    ,25.0     /
0820       DATA (AE(I, 2,11),I=1,10)  /
0821      +  .000E+00,.000E+00,.000E+00,15.0    ,18.0    ,19.6    ,22.4    ,
0822      +  23.8    ,25.2    ,26.2     /
0823       DATA (AE(I, 2,12),I=1,10)  /
0824      +  .000E+00,.000E+00,.000E+00,16.2    ,17.3    ,19.4    ,22.2    ,
0825      +  24.0    ,25.7    ,27.2     /
0826       DATA (AE(I, 2,13),I=1,10)  /
0827      +  .000E+00,.000E+00,.000E+00,.000E+00,17.8    ,19.8    ,22.9    ,
0828      +  24.6    ,26.2    ,27.7     /
0829       DATA (AE(I, 2,14),I=1,10)  /
0830      +  .000E+00,.000E+00,.000E+00,.000E+00,19.1    ,20.4    ,23.7    ,
0831      +  25.7    ,27.6    ,29.1     /
0832       DATA (AE(I, 2,15),I=1,10)  /
0833      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,20.5    ,24.1    ,
0834      +  26.1    ,28.1    ,29.9     /
0835       DATA (AE(I, 2,16),I=1,10)  /
0836      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,20.9    ,23.9    ,
0837      +  26.4    ,28.7    ,30.7     /
0838       DATA (AE(I, 2,17),I=1,10)  /
0839      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,22.4    ,24.2    ,
0840      +  26.7    ,29.0    ,30.9     /
0841       DATA (AE(I, 2,18),I=1,10)  /
0842      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,24.8    ,
0843      +  27.3    ,29.7    ,31.8     /
0844       DATA (AE(I, 2,19),I=1,10)  /
0845      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,26.1    ,
0846      +  27.3    ,29.9    ,32.3     /
0847       DATA (AE(I, 2,20),I=1,10)  /
0848      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
0849      +  27.4    ,30.1    ,32.6     /
0850       DATA (AE(I, 3, 1),I=1,10)  /
0851      +  11.0    ,11.0    ,11.7    ,11.3    ,11.9    ,11.4    ,12.1    ,
0852      +  11.7    ,11.5    ,11.0     /
0853       DATA (AE(I, 3, 2),I=1,10)  /
0854      +  11.2    ,12.0    ,12.7    ,12.9    ,13.6    ,13.7    ,14.4    ,
0855      +  14.6    ,14.7    ,14.6     /
0856       DATA (AE(I, 3, 3),I=1,10)  /
0857      +  12.1    ,12.6    ,13.7    ,13.9    ,15.0    ,15.2    ,16.3    ,
0858      +  16.5    ,16.7    ,16.7     /
0859       DATA (AE(I, 3, 4),I=1,10)  /
0860      +  12.6    ,11.3    ,12.4    ,13.0    ,13.8    ,14.2    ,15.0    ,
0861      +  15.6    ,16.1    ,16.6     /
0862       DATA (AE(I, 3, 5),I=1,10)  /
0863      +  .000E+00,12.6    ,13.7    ,14.4    ,15.3    ,16.0    ,16.8    ,
0864      +  17.5    ,18.1    ,18.6     /
0865       DATA (AE(I, 3, 6),I=1,10)  /
0866      +  .000E+00,14.0    ,14.6    ,15.8    ,17.4    ,18.4    ,19.8    ,
0867      +  20.6    ,21.5    ,22.2     /
0868       DATA (AE(I, 3, 7),I=1,10)  /
0869      +  .000E+00,16.0    ,15.2    ,16.3    ,18.3    ,19.3    ,21.1    ,
0870      +  22.0    ,22.8    ,23.5     /
0871       DATA (AE(I, 3, 8),I=1,10)  /
0872      +  .000E+00,.000E+00,15.6    ,15.1    ,17.2    ,18.6    ,20.6    ,
0873      +  21.8    ,22.9    ,23.8     /
0874       DATA (AE(I, 3, 9),I=1,10)  /
0875      +  .000E+00,.000E+00,17.8    ,16.3    ,18.8    ,20.1    ,22.5    ,
0876      +  23.6    ,24.7    ,25.6     /
0877       DATA (AE(I, 3,10),I=1,10)  /
0878      +  .000E+00,.000E+00,.000E+00,17.5    ,19.0    ,20.7    ,23.1    ,
0879      +  24.5    ,25.8    ,26.8     /
0880       DATA (AE(I, 3,11),I=1,10)  /
0881      +  .000E+00,.000E+00,.000E+00,19.2    ,19.4    ,21.1    ,23.8    ,
0882      +  25.4    ,26.8    ,28.0     /
0883       DATA (AE(I, 3,12),I=1,10)  /
0884      +  .000E+00,.000E+00,.000E+00,20.7    ,19.6    ,19.7    ,22.4    ,
0885      +  24.4    ,26.2    ,27.9     /
0886       DATA (AE(I, 3,13),I=1,10)  /
0887      +  .000E+00,.000E+00,.000E+00,.000E+00,21.6    ,20.4    ,23.2    ,
0888      +  25.1    ,26.9    ,28.5     /
0889       DATA (AE(I, 3,14),I=1,10)  /
0890      +  .000E+00,.000E+00,.000E+00,.000E+00,23.5    ,22.0    ,23.8    ,
0891      +  26.1    ,28.1    ,29.9     /
0892       DATA (AE(I, 3,15),I=1,10)  /
0893      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,23.7    ,24.2    ,
0894      +  26.3    ,28.5    ,30.4     /
0895       DATA (AE(I, 3,16),I=1,10)  /
0896      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,25.4    ,24.8    ,
0897      +  25.6    ,28.1    ,30.5     /
0898       DATA (AE(I, 3,17),I=1,10)  /
0899      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,26.9    ,26.8    ,
0900      +  26.1    ,28.4    ,30.8     /
0901       DATA (AE(I, 3,18),I=1,10)  /
0902      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,28.8    ,
0903      +  27.6    ,29.0    ,31.5     /
0904       DATA (AE(I, 3,19),I=1,10)  /
0905      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,30.5    ,
0906      +  29.2    ,28.9    ,31.5     /
0907       DATA (AE(I, 3,20),I=1,10)  /
0908      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
0909      +  31.0    ,30.0    ,31.7     /
0910       DATA (AE(I, 4, 1),I=1,10)  /
0911      +  13.0    ,13.2    ,14.8    ,14.2    ,14.2    ,14.1    ,14.5    ,
0912      +  14.4    ,14.3    ,14.0     /
0913       DATA (AE(I, 4, 2),I=1,10)  /
0914      +  13.5    ,14.5    ,16.1    ,15.9    ,16.0    ,16.3    ,16.8    ,
0915      +  17.0    ,17.1    ,17.2     /
0916       DATA (AE(I, 4, 3),I=1,10)  /
0917      +  14.9    ,15.3    ,17.2    ,17.1    ,17.5    ,17.8    ,18.6    ,
0918      +  18.9    ,19.1    ,19.3     /
0919       DATA (AE(I, 4, 4),I=1,10)  /
0920      +  15.1    ,13.5    ,16.4    ,16.7    ,16.4    ,17.3    ,17.8    ,
0921      +  18.5    ,19.0    ,19.6     /
0922       DATA (AE(I, 4, 5),I=1,10)  /
0923      +  .000E+00,15.6    ,17.5    ,17.7    ,17.8    ,18.6    ,19.2    ,
0924      +  19.9    ,20.3    ,21.1     /
0925       DATA (AE(I, 4, 6),I=1,10)  /
0926      +  .000E+00,18.0    ,18.4    ,19.2    ,19.8    ,20.9    ,22.0    ,
0927      +  23.1    ,23.6    ,24.7     /
0928       DATA (AE(I, 4, 7),I=1,10)  /
0929      +  .000E+00,27.4    ,19.1    ,19.8    ,20.7    ,21.8    ,23.2    ,
0930      +  24.4    ,24.9    ,25.9     /
0931       DATA (AE(I, 4, 8),I=1,10)  /
0932      +  .000E+00,.000E+00,18.9    ,18.9    ,19.3    ,21.1    ,22.5    ,
0933      +  24.0    ,24.7    ,26.0     /
0934       DATA (AE(I, 4, 9),I=1,10)  /
0935      +  .000E+00,.000E+00,21.1    ,19.7    ,20.7    ,22.3    ,24.0    ,
0936      +  25.6    ,26.3    ,27.7     /
0937       DATA (AE(I, 4,10),I=1,10)  /
0938      +  .000E+00,.000E+00,.000E+00,21.0    ,21.1    ,22.9    ,24.6    ,
0939      +  26.5    ,27.3    ,29.0     /
0940       DATA (AE(I, 4,11),I=1,10)  /
0941      +  .000E+00,.000E+00,.000E+00,21.3    ,22.4    ,23.1    ,25.0    ,
0942      +  27.1    ,27.9    ,29.8     /
0943       DATA (AE(I, 4,12),I=1,10)  /
0944      +  .000E+00,.000E+00,.000E+00,36.6    ,21.5    ,22.2    ,23.1    ,
0945      +  25.6    ,26.8    ,29.1     /
0946       DATA (AE(I, 4,13),I=1,10)  /
0947      +  .000E+00,.000E+00,.000E+00,.000E+00,22.9    ,23.1    ,23.7    ,
0948      +  26.2    ,27.3    ,29.6     /
0949       DATA (AE(I, 4,14),I=1,10)  /
0950      +  .000E+00,.000E+00,.000E+00,.000E+00,30.5    ,23.6    ,25.0    ,
0951      +  26.9    ,28.2    ,30.7     /
0952       DATA (AE(I, 4,15),I=1,10)  /
0953      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,25.4    ,26.2    ,
0954      +  27.2    ,28.3    ,31.0     /
0955       DATA (AE(I, 4,16),I=1,10)  /
0956      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,24.5    ,25.9    ,
0957      +  27.4    ,27.6    ,30.7     /
0958       DATA (AE(I, 4,17),I=1,10)  /
0959      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,43.3    ,28.4    ,
0960      +  27.5    ,27.9    ,30.9     /
0961       DATA (AE(I, 4,18),I=1,10)  /
0962      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,27.2    ,
0963      +  29.1    ,29.0    ,31.4     /
0964       DATA (AE(I, 4,19),I=1,10)  /
0965      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,51.3    ,
0966      +  30.6    ,29.5    ,31.4     /
0967       DATA (AE(I, 4,20),I=1,10)  /
0968      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
0969      +  28.8    ,30.6    ,32.4     /
0970       DATA (AE(I, 5, 1),I=1,10)  /
0971      +  15.0    ,14.9    ,15.5    ,15.4    ,15.9    ,15.8    ,16.2    ,
0972      +  16.2    ,16.1    ,15.9     /
0973       DATA (AE(I, 5, 2),I=1,10)  /
0974      +  15.4    ,16.1    ,17.0    ,17.4    ,18.0    ,18.2    ,18.7    ,
0975      +  18.9    ,19.0    ,19.1     /
0976       DATA (AE(I, 5, 3),I=1,10)  /
0977      +  17.1    ,17.2    ,18.3    ,18.7    ,19.3    ,19.6    ,20.3    ,
0978      +  20.6    ,20.8    ,20.9     /
0979       DATA (AE(I, 5, 4),I=1,10)  /
0980      +  14.7    ,14.8    ,15.0    ,16.0    ,17.0    ,17.7    ,18.1    ,
0981      +  19.0    ,19.4    ,20.0     /
0982       DATA (AE(I, 5, 5),I=1,10)  /
0983      +  .000E+00,16.7    ,17.6    ,18.1    ,18.6    ,19.2    ,19.7    ,
0984      +  20.4    ,20.8    ,21.2     /
0985       DATA (AE(I, 5, 6),I=1,10)  /
0986      +  .000E+00,17.8    ,18.2    ,19.2    ,20.0    ,21.0    ,21.9    ,
0987      +  23.0    ,23.6    ,24.3     /
0988       DATA (AE(I, 5, 7),I=1,10)  /
0989      +  .000E+00,35.2    ,18.9    ,20.3    ,20.6    ,21.5    ,22.6    ,
0990      +  23.7    ,24.2    ,24.7     /
0991       DATA (AE(I, 5, 8),I=1,10)  /
0992      +  .000E+00,.000E+00,16.4    ,18.9    ,18.8    ,19.6    ,20.7    ,
0993      +  22.3    ,23.1    ,23.9     /
0994       DATA (AE(I, 5, 9),I=1,10)  /
0995      +  .000E+00,.000E+00,33.9    ,19.8    ,20.3    ,20.7    ,21.9    ,
0996      +  23.4    ,24.1    ,24.8     /
0997       DATA (AE(I, 5,10),I=1,10)  /
0998      +  .000E+00,.000E+00,.000E+00,18.0    ,20.0    ,21.4    ,22.0    ,
0999      +  23.8    ,24.6    ,25.4     /
1000       DATA (AE(I, 5,11),I=1,10)  /
1001      +  .000E+00,.000E+00,.000E+00,26.4    ,20.4    ,21.2    ,22.3    ,
1002      +  23.8    ,24.7    ,25.5     /
1003       DATA (AE(I, 5,12),I=1,10)  /
1004      +  .000E+00,.000E+00,.000E+00,41.7    ,18.2    ,19.8    ,21.1    ,
1005      +  22.6    ,23.4    ,24.6     /
1006       DATA (AE(I, 5,13),I=1,10)  /
1007      +  .000E+00,.000E+00,.000E+00,.000E+00,22.5    ,20.0    ,21.7    ,
1008      +  22.8    ,23.7    ,24.7     /
1009       DATA (AE(I, 5,14),I=1,10)  /
1010      +  .000E+00,.000E+00,.000E+00,.000E+00,54.1    ,19.9    ,21.9    ,
1011      +  23.2    ,24.3    ,25.3     /
1012       DATA (AE(I, 5,15),I=1,10)  /
1013      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,21.2    ,22.2    ,
1014      +  23.6    ,24.9    ,25.5     /
1015       DATA (AE(I, 5,16),I=1,10)  /
1016      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,44.9    ,21.9    ,
1017      +  23.8    ,25.2    ,25.6     /
1018       DATA (AE(I, 5,17),I=1,10)  /
1019      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,47.8    ,22.7    ,
1020      +  23.8    ,24.9    ,26.3     /
1021       DATA (AE(I, 5,18),I=1,10)  /
1022      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,35.5    ,
1023      +  23.9    ,25.9    ,26.6     /
1024       DATA (AE(I, 5,19),I=1,10)  /
1025      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,64.3    ,
1026      +  24.1    ,25.7    ,27.1     /
1027       DATA (AE(I, 5,20),I=1,10)  /
1028      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
1029      +  34.0    ,25.7    ,27.7     /
1030       DATA (AE(I, 6, 1),I=1,10)  /
1031      +  16.6    ,16.5    ,16.8    ,16.7    ,17.0    ,16.5    ,16.7    ,
1032      +  18.3    ,18.9    ,19.0     /
1033       DATA (AE(I, 6, 2),I=1,10)  /
1034      +  16.2    ,16.6    ,17.2    ,17.4    ,17.9    ,17.4    ,17.7    ,
1035      +  20.7    ,22.0    ,22.6     /
1036       DATA (AE(I, 6, 3),I=1,10)  /
1037      +  18.9    ,18.7    ,18.8    ,18.6    ,18.9    ,18.6    ,18.9    ,
1038      +  21.0    ,22.3    ,22.9     /
1039       DATA (AE(I, 6, 4),I=1,10)  /
1040      +  18.3    ,12.7    ,14.2    ,15.0    ,15.7    ,16.1    ,16.3    ,
1041      +  16.5    ,17.9    ,19.0     /
1042       DATA (AE(I, 6, 5),I=1,10)  /
1043      +  .000E+00,15.7    ,15.1    ,15.3    ,16.5    ,16.4    ,16.4    ,
1044      +  17.0    ,18.3    ,19.4     /
1045       DATA (AE(I, 6, 6),I=1,10)  /
1046      +  .000E+00,22.9    ,14.9    ,15.2    ,16.2    ,16.9    ,17.4    ,
1047      +  18.2    ,19.5    ,21.1     /
1048       DATA (AE(I, 6, 7),I=1,10)  /
1049      +  .000E+00,40.7    ,18.4    ,15.9    ,17.1    ,17.7    ,18.9    ,
1050      +  19.5    ,20.3    ,21.1     /
1051       DATA (AE(I, 6, 8),I=1,10)  /
1052      +  .000E+00,.000E+00,23.3    ,16.2    ,16.3    ,17.3    ,18.7    ,
1053      +  19.5    ,20.3    ,21.1     /
1054       DATA (AE(I, 6, 9),I=1,10)  /
1055      +  .000E+00,.000E+00,49.2    ,19.0    ,19.1    ,19.4    ,20.2    ,
1056      +  20.8    ,21.6    ,22.0     /
1057       DATA (AE(I, 6,10),I=1,10)  /
1058      +  .000E+00,.000E+00,.000E+00,27.2    ,21.2    ,20.8    ,21.4    ,
1059      +  22.3    ,22.8    ,23.3     /
1060       DATA (AE(I, 6,11),I=1,10)  /
1061      +  .000E+00,.000E+00,.000E+00,45.6    ,25.0    ,22.8    ,23.9    ,
1062      +  23.6    ,24.3    ,24.4     /
1063       DATA (AE(I, 6,12),I=1,10)  /
1064      +  .000E+00,.000E+00,.000E+00,45.8    ,29.7    ,25.1    ,25.3    ,
1065      +  25.3    ,26.0    ,26.3     /
1066       DATA (AE(I, 6,13),I=1,10)  /
1067      +  .000E+00,.000E+00,.000E+00,.000E+00,42.7    ,29.0    ,28.0    ,
1068      +  27.0    ,27.2    ,27.6     /
1069       DATA (AE(I, 6,14),I=1,10)  /
1070      +  .000E+00,.000E+00,.000E+00,.000E+00,62.0    ,32.0    ,30.0    ,
1071      +  29.8    ,29.5    ,29.6     /
1072       DATA (AE(I, 6,15),I=1,10)  /
1073      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,44.5    ,34.4    ,
1074      +  32.7    ,31.5    ,31.8     /
1075       DATA (AE(I, 6,16),I=1,10)  /
1076      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,75.6    ,37.1    ,
1077      +  34.6    ,34.4    ,34.4     /
1078       DATA (AE(I, 6,17),I=1,10)  /
1079      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,51.2    ,45.2    ,
1080      +  39.0    ,37.5    ,36.4     /
1081       DATA (AE(I, 6,18),I=1,10)  /
1082      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,74.9    ,
1083      +  42.3    ,39.9    ,38.3     /
1084       DATA (AE(I, 6,19),I=1,10)  /
1085      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,69.5    ,
1086      +  50.7    ,42.3    ,41.4     /
1087       DATA (AE(I, 6,20),I=1,10)  /
1088      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
1089      +  66.3    ,48.0    ,43.4     /
1090       DATA (AE(I, 7, 1),I=1,10)  /
1091      +  27.0    ,25.8    ,26.3    ,26.2    ,26.7    ,26.7    ,27.1    ,
1092      +  27.1    ,27.2    ,19.0     /
1093       DATA (AE(I, 7, 2),I=1,10)  /
1094      +  29.1    ,28.9    ,29.7    ,30.3    ,31.0    ,31.4    ,32.0    ,
1095      +  32.3    ,32.7    ,22.6     /
1096       DATA (AE(I, 7, 3),I=1,10)  /
1097      +  31.6    ,29.7    ,30.9    ,31.4    ,32.5    ,33.1    ,34.0    ,
1098      +  34.6    ,35.1    ,22.9     /
1099       DATA (AE(I, 7, 4),I=1,10)  /
1100      +  27.4    ,19.9    ,20.8    ,22.8    ,24.6    ,26.4    ,28.2    ,
1101      +  29.6    ,30.8    ,19.0     /
1102       DATA (AE(I, 7, 5),I=1,10)  /
1103      +  .000E+00,24.6    ,24.1    ,25.0    ,27.2    ,28.7    ,30.7    ,
1104      +  31.8    ,32.9    ,19.4     /
1105       DATA (AE(I, 7, 6),I=1,10)  /
1106      +  .000E+00,35.6    ,25.2    ,25.6    ,27.9    ,30.4    ,32.7    ,
1107      +  34.6    ,36.3    ,21.1     /
1108       DATA (AE(I, 7, 7),I=1,10)  /
1109      +  .000E+00,45.4    ,30.9    ,28.2    ,29.0    ,31.2    ,34.0    ,
1110      +  35.8    ,37.4    ,21.1     /
1111       DATA (AE(I, 7, 8),I=1,10)  /
1112      +  .000E+00,.000E+00,38.2    ,29.6    ,29.4    ,30.3    ,33.2    ,
1113      +  35.5    ,37.6    ,21.1     /
1114       DATA (AE(I, 7, 9),I=1,10)  /
1115      +  .000E+00,.000E+00,59.3    ,34.5    ,33.7    ,32.9    ,35.4    ,
1116      +  37.6    ,39.6    ,22.0     /
1117       DATA (AE(I, 7,10),I=1,10)  /
1118      +  .000E+00,.000E+00,.000E+00,44.5    ,37.8    ,37.5    ,37.2    ,
1119      +  39.0    ,41.4    ,23.3     /
1120       DATA (AE(I, 7,11),I=1,10)  /
1121      +  .000E+00,.000E+00,.000E+00,67.0    ,43.6    ,42.0    ,40.8    ,
1122      +  41.4    ,43.0    ,24.4     /
1123       DATA (AE(I, 7,12),I=1,10)  /
1124      +  .000E+00,.000E+00,.000E+00,49.9    ,50.9    ,44.6    ,43.9    ,
1125      +  44.2    ,44.2    ,26.3     /
1126       DATA (AE(I, 7,13),I=1,10)  /
1127      +  .000E+00,.000E+00,.000E+00,.000E+00,67.2    ,50.5    ,48.7    ,
1128      +  48.1    ,47.2    ,27.6     /
1129       DATA (AE(I, 7,14),I=1,10)  /
1130      +  .000E+00,.000E+00,.000E+00,.000E+00,68.1    ,55.2    ,52.3    ,
1131      +  51.5    ,51.6    ,29.6     /
1132       DATA (AE(I, 7,15),I=1,10)  /
1133      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,68.7    ,58.6    ,
1134      +  56.5    ,55.7    ,31.8     /
1135       DATA (AE(I, 7,16),I=1,10)  /
1136      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,89.3    ,62.9    ,
1137      +  60.0    ,59.1    ,34.4     /
1138       DATA (AE(I, 7,17),I=1,10)  /
1139      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,56.0    ,72.9    ,
1140      +  66.3    ,64.2    ,36.4     /
1141       DATA (AE(I, 7,18),I=1,10)  /
1142      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,105.    ,
1143      +  71.3    ,68.3    ,38.3     /
1144       DATA (AE(I, 7,19),I=1,10)  /
1145      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,73.4    ,
1146      +  76.8    ,72.4    ,41.4     /
1147       DATA (AE(I, 7,20),I=1,10)  /
1148      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
1149      +  107.    ,79.9    ,43.4     /
1150       DATA (AE(I, 8, 1),I=1,10)  /
1151      +  35.5    ,35.3    ,35.7    ,35.7    ,36.3    ,36.3    ,36.7    ,
1152      +  36.7    ,36.7    ,19.0     /
1153       DATA (AE(I, 8, 2),I=1,10)  /
1154      +  40.6    ,41.4    ,41.9    ,42.3    ,43.2    ,43.5    ,44.0    ,
1155      +  44.3    ,44.5    ,22.6     /
1156       DATA (AE(I, 8, 3),I=1,10)  /
1157      +  45.4    ,45.7    ,46.4    ,47.0    ,48.1    ,48.7    ,49.4    ,
1158      +  49.8    ,50.2    ,22.9     /
1159       DATA (AE(I, 8, 4),I=1,10)  /
1160      +  43.9    ,44.3    ,43.4    ,45.1    ,47.3    ,48.7    ,49.6    ,
1161      +  50.5    ,51.3    ,19.0     /
1162       DATA (AE(I, 8, 5),I=1,10)  /
1163      +  .000E+00,49.3    ,49.6    ,50.5    ,53.2    ,54.2    ,55.4    ,
1164      +  56.1    ,56.8    ,19.4     /
1165       DATA (AE(I, 8, 6),I=1,10)  /
1166      +  .000E+00,59.1    ,53.0    ,55.4    ,58.0    ,60.0    ,61.2    ,
1167      +  62.5    ,63.6    ,21.1     /
1168       DATA (AE(I, 8, 7),I=1,10)  /
1169      +  .000E+00,54.5    ,57.1    ,59.2    ,62.3    ,64.4    ,66.0    ,
1170      +  67.3    ,68.5    ,21.1     /
1171       DATA (AE(I, 8, 8),I=1,10)  /
1172      +  .000E+00,.000E+00,65.9    ,62.1    ,65.1    ,67.6    ,69.4    ,
1173      +  71.1    ,72.6    ,21.1     /
1174       DATA (AE(I, 8, 9),I=1,10)  /
1175      +  .000E+00,.000E+00,72.2    ,67.1    ,70.5    ,73.1    ,75.1    ,
1176      +  76.8    ,78.4    ,22.0     /
1177       DATA (AE(I, 8,10),I=1,10)  /
1178      +  .000E+00,.000E+00,.000E+00,80.1    ,75.0    ,78.0    ,80.0    ,
1179      +  82.1    ,83.9    ,23.3     /
1180       DATA (AE(I, 8,11),I=1,10)  /
1181      +  .000E+00,.000E+00,.000E+00,94.5    ,82.2    ,82.8    ,85.1    ,
1182      +  87.3    ,89.2    ,24.4     /
1183       DATA (AE(I, 8,12),I=1,10)  /
1184      +  .000E+00,.000E+00,.000E+00,56.8    ,92.5    ,87.2    ,89.4    ,
1185      +  91.9    ,94.1    ,26.3     /
1186       DATA (AE(I, 8,13),I=1,10)  /
1187      +  .000E+00,.000E+00,.000E+00,.000E+00,116.    ,96.2    ,94.4    ,
1188      +  97.0    ,99.2    ,27.6     /
1189       DATA (AE(I, 8,14),I=1,10)  /
1190      +  .000E+00,.000E+00,.000E+00,.000E+00,78.1    ,104.    ,102.    ,
1191      +  102.    ,105.    ,29.6     /
1192       DATA (AE(I, 8,15),I=1,10)  /
1193      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,128.    ,111.    ,
1194      +  109.    ,110.    ,31.8     /
1195       DATA (AE(I, 8,16),I=1,10)  /
1196      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,104.    ,118.    ,
1197      +  117.    ,115.    ,34.4     /
1198       DATA (AE(I, 8,17),I=1,10)  /
1199      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,64.4    ,138.    ,
1200      +  124.    ,122.    ,36.4     /
1201       DATA (AE(I, 8,18),I=1,10)  /
1202      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,133.    ,
1203      +  133.    ,132.    ,38.3     /
1204       DATA (AE(I, 8,19),I=1,10)  /
1205      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,83.6    ,
1206      +  146.    ,139.    ,41.4     /
1207       DATA (AE(I, 8,20),I=1,10)  /
1208      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
1209      +  166.    ,147.    ,43.4     /
1210       DATA (AE(I, 9, 1),I=1,10)  /
1211      +  43.3    ,43.2    ,43.6    ,43.8    ,44.1    ,44.3    ,44.7    ,
1212      +  44.8    ,44.8    ,19.0     /
1213       DATA (AE(I, 9, 2),I=1,10)  /
1214      +  50.9    ,51.4    ,52.0    ,52.6    ,53.1    ,53.6    ,54.2    ,
1215      +  54.5    ,54.7    ,22.6     /
1216       DATA (AE(I, 9, 3),I=1,10)  /
1217      +  58.0    ,58.4    ,59.3    ,60.1    ,60.7    ,61.5    ,62.3    ,
1218      +  62.7    ,63.1    ,22.9     /
1219       DATA (AE(I, 9, 4),I=1,10)  /
1220      +  62.0    ,63.9    ,63.7    ,65.7    ,65.5    ,67.5    ,68.2    ,
1221      +  68.9    ,69.7    ,19.0     /
1222       DATA (AE(I, 9, 5),I=1,10)  /
1223      +  .000E+00,72.2    ,72.5    ,74.2    ,74.2    ,76.1    ,77.0    ,
1224      +  77.8    ,78.6    ,19.4     /
1225       DATA (AE(I, 9, 6),I=1,10)  /
1226      +  .000E+00,80.4    ,80.5    ,83.1    ,83.0    ,85.5    ,86.8    ,
1227      +  88.1    ,89.2    ,21.1     /
1228       DATA (AE(I, 9, 7),I=1,10)  /
1229      +  .000E+00,63.4    ,88.5    ,91.3    ,91.1    ,94.0    ,95.8    ,
1230      +  97.3    ,98.6    ,21.1     /
1231       DATA (AE(I, 9, 8),I=1,10)  /
1232      +  .000E+00,.000E+00,98.8    ,98.6    ,97.8    ,102.    ,104.    ,
1233      +  106.    ,108.    ,21.1     /
1234       DATA (AE(I, 9, 9),I=1,10)  /
1235      +  .000E+00,.000E+00,84.1    ,107.    ,107.    ,111.    ,113.    ,
1236      +  116.    ,117.    ,22.0     /
1237       DATA (AE(I, 9,10),I=1,10)  /
1238      +  .000E+00,.000E+00,.000E+00,116.    ,115.    ,119.    ,122.    ,
1239      +  125.    ,127.    ,23.3     /
1240       DATA (AE(I, 9,11),I=1,10)  /
1241      +  .000E+00,.000E+00,.000E+00,111.    ,123.    ,127.    ,131.    ,
1242      +  134.    ,137.    ,24.4     /
1243       DATA (AE(I, 9,12),I=1,10)  /
1244      +  .000E+00,.000E+00,.000E+00,65.6    ,136.    ,135.    ,140.    ,
1245      +  143.    ,146.    ,26.3     /
1246       DATA (AE(I, 9,13),I=1,10)  /
1247      +  .000E+00,.000E+00,.000E+00,.000E+00,146.    ,144.    ,149.    ,
1248      +  152.    ,155.    ,27.6     /
1249       DATA (AE(I, 9,14),I=1,10)  /
1250      +  .000E+00,.000E+00,.000E+00,.000E+00,88.7    ,152.    ,158.    ,
1251      +  162.    ,165.    ,29.6     /
1252       DATA (AE(I, 9,15),I=1,10)  /
1253      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,181.    ,167.    ,
1254      +  171.    ,174.    ,31.8     /
1255       DATA (AE(I, 9,16),I=1,10)  /
1256      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,117.    ,174.    ,
1257      +  180.    ,183.    ,34.4     /
1258       DATA (AE(I, 9,17),I=1,10)  /
1259      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,72.0    ,201.    ,
1260      +  189.    ,192.    ,36.4     /
1261       DATA (AE(I, 9,18),I=1,10)  /
1262      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,151.    ,
1263      +  198.    ,201.    ,38.3     /
1264       DATA (AE(I, 9,19),I=1,10)  /
1265      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,95.2    ,
1266      +  220.    ,210.    ,41.4     /
1267       DATA (AE(I, 9,20),I=1,10)  /
1268      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
1269      +  192.    ,217.    ,43.4     /
1270       DATA (AE(I,10, 1),I=1,10)  /
1271      +  62.1    ,62.1    ,62.6    ,62.9    ,63.3    ,63.3    ,64.0    ,
1272      +  64.0    ,64.0    ,19.0     /
1273       DATA (AE(I,10, 2),I=1,10)  /
1274      +  75.1    ,75.4    ,76.3    ,76.8    ,77.6    ,77.9    ,78.8    ,
1275      +  79.0    ,79.3    ,22.6     /
1276       DATA (AE(I,10, 3),I=1,10)  /
1277      +  87.5    ,88.3    ,89.4    ,90.2    ,91.3    ,91.9    ,93.0    ,
1278      +  93.5    ,93.9    ,22.9     /
1279       DATA (AE(I,10, 4),I=1,10)  /
1280      +  104.    ,104.    ,105.    ,106.    ,107.    ,108.    ,109.    ,
1281      +  110.    ,110.    ,19.0     /
1282       DATA (AE(I,10, 5),I=1,10)  /
1283      +  .000E+00,122.    ,122.    ,123.    ,124.    ,125.    ,126.    ,
1284      +  127.    ,128.    ,19.4     /
1285       DATA (AE(I,10, 6),I=1,10)  /
1286      +  .000E+00,138.    ,139.    ,140.    ,142.    ,143.    ,144.    ,
1287      +  146.    ,147.    ,21.1     /
1288       DATA (AE(I,10, 7),I=1,10)  /
1289      +  .000E+00,85.3    ,158.    ,159.    ,161.    ,162.    ,164.    ,
1290      +  166.    ,167.    ,21.1     /
1291       DATA (AE(I,10, 8),I=1,10)  /
1292      +  .000E+00,.000E+00,176.    ,177.    ,179.    ,181.    ,183.    ,
1293      +  184.    ,186.    ,21.1     /
1294       DATA (AE(I,10, 9),I=1,10)  /
1295      +  .000E+00,.000E+00,114.    ,199.    ,201.    ,202.    ,205.    ,
1296      +  206.    ,207.    ,22.0     /
1297       DATA (AE(I,10,10),I=1,10)  /
1298      +  .000E+00,.000E+00,.000E+00,218.    ,219.    ,220.    ,224.    ,
1299      +  225.    ,226.    ,23.3     /
1300       DATA (AE(I,10,11),I=1,10)  /
1301      +  .000E+00,.000E+00,.000E+00,150.    ,238.    ,238.    ,243.    ,
1302      +  244.    ,245.    ,24.4     /
1303       DATA (AE(I,10,12),I=1,10)  /
1304      +  .000E+00,.000E+00,.000E+00,85.8    ,255.    ,255.    ,261.    ,
1305      +  262.    ,263.    ,26.3     /
1306       DATA (AE(I,10,13),I=1,10)  /
1307      +  .000E+00,.000E+00,.000E+00,.000E+00,195.    ,272.    ,279.    ,
1308      +  279.    ,280.    ,27.6     /
1309       DATA (AE(I,10,14),I=1,10)  /
1310      +  .000E+00,.000E+00,.000E+00,.000E+00,115.    ,290.    ,296.    ,
1311      +  297.    ,298.    ,29.6     /
1312       DATA (AE(I,10,15),I=1,10)  /
1313      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,263.    ,313.    ,
1314      +  314.    ,315.    ,31.8     /
1315       DATA (AE(I,10,16),I=1,10)  /
1316      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,150.    ,330.    ,
1317      +  331.    ,332.    ,34.4     /
1318       DATA (AE(I,10,17),I=1,10)  /
1319      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,90.0    ,319.    ,
1320      +  349.    ,349.    ,36.4     /
1321       DATA (AE(I,10,18),I=1,10)  /
1322      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,196.    ,
1323      +  366.    ,367.    ,38.3     /
1324       DATA (AE(I,10,19),I=1,10)  /
1325      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,122.    ,
1326      +  387.    ,384.    ,41.4     /
1327       DATA (AE(I,10,20),I=1,10)  /
1328      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
1329      +  247.    ,401.    ,43.4     /
1330       DATA (ERES(I, 1),I=1,10)  / 10*0./
1331       DATA (ERES(I, 2),I=1,10)  / 10*0./
1332       DATA (ERES(I, 3),I=1,10)  / 10*0./
1333       DATA (ERES(I, 4),I=1,10)  / 10*0./
1334       DATA (ERES(I, 5),I=1,10)  / 10*0./
1335       DATA (ERES(I, 6),I=1,10)  /
1336      +     0.000,   0.000,   0.000,   0.000,   0.000,   0.000,   0.000,
1337      +     2.780,   2.880,   2.890 /
1338       DATA (ERES(I, 7),I=1,10)  /
1339      +     1.500,   2.460,   2.510,   2.610,   2.700,   2.920,   3.070,
1340      +     3.200,   3.330,   2.890 /
1341       DATA (ERES(I, 8),I=1,10)  /
1342      +     4.470,   4.350,   4.390,   4.550,   4.660,   4.890,   4.980,
1343      +     5.100,   5.220,   2.890 /
1344       DATA (ERES(I, 9),I=1,10)  /
1345      +     7.480,   7.380,   7.370,   7.480,   7.510,   7.630,   7.660,
1346      +     7.750,   7.820,   2.890 /
1347       DATA (ERES(I,10),I=1,10)  /
1348      +    15.270,  15.190,  15.200,  15.370,  15.380,  15.430,  15.540,
1349      +    15.590,  15.630,   2.890 /
1350       END
1351 C->
1352       SUBROUTINE FRAGM (IAT,IAP, NW,B, NF, IAF)
1353 C...Nuclear Fragmentation, Abrasion-ablation model,
1354 C...Based on Jon Engel's routines ABRABL
1355 C...This most recent version adds for all prefragment
1356 C...masses > 10 the model calculation for the fragment
1357 C...mass distribution and the energy carried by the fragment
1358 C...of W. Friedmann
1359 C...The average values are used to implement the model
1360 C...in the montecarlo fashion / TSS, Dec '91
1361 C.
1362 C.  INPUT: IAP = mass of incident nucleus
1363 C.         IAT = mass of target   nucleus
1364 C.         NW = number of wounded nucleons in the beam nucleus
1365 C.         B  = impact parameter in the interaction
1366 C.
1367 C.  OUTPUT : NF = number of fragments  of the spectator nucleus
1368 C.           IAF(1:NF) = mass number of each fragment
1369 C.           PF(3,60) in common block /FRAGMENTS/ contains
1370 C.           the three momentum components (MeV/c) of each
1371 C.           fragment in the projectile frame
1372 C..............................................................
1373       SAVE
1374 
1375       COMMON /FRAGMENTS/ PPP(3,60)
1376       COMMON /FRAGMOD/A(10,10,20),AE(10,10,20),ERES(10,10),NFLAGG(10,10)
1377       DIMENSION IAF(60)
1378       DIMENSION AA(10), EAA(10)
1379       DATA AA/10.,15.,20.,25.,30.,35.,40.,45.,50.,56./
1380       DATA EAA/1.,2.,4.,6.,8.,10.,12.,16.,20.,30/
1381       AP=IAP
1382       AT=IAT
1383       NPF = IAP - NW
1384       IF (NPF .EQ. 0) THEN
1385          NF = 0
1386          RETURN
1387       ENDIF
1388 
1389       EB = ESTAR(AP,AT, B)
1390       EBP = ESTARP (NPF, NW)
1391 C CONTRIBUTION TO E* FROM ENERGY DEPOSITED BY SECONDARIES
1392       EB = EB + EBP
1393 C TOTAL E* IS THE SUM OF THE TWO COMPONENTS
1394 
1395 C.....Prefragment transverse momentum (MeV/nucleon)...
1396             FK = FERMK(AP)
1397 C FERMI MOMENTUM OF THE PROJECTILE NUCLEUS
1398             IF (NW .LT. IAP) THEN
1399             SIG = FK*SQRT(NW*NPF/(AP-1.))/3.162
1400 C GAUSSIAN SIGMA IN ALL THREE DIRECTION
1401             ELSE
1402             SIG = FK/3.162
1403 C THIS IS NOT CORRECT, TOO LARGE !!!!!!!!!!!!!!
1404             ENDIF
1405 ctp060203             PPFX = SIG*GASDEV(0)/NPF
1406 ctp060203             PPFY = SIG*GASDEV(0)/NPF
1407              PPFX = SIG*GASDEV(NW)/NPF
1408              PPFY = SIG*GASDEV(IAP)/NPF
1409 C THREE MOMENTUM COMPONENTS PER NUCLEON FOR THE PREFRAGMENT
1410 
1411 C.............Crude model for small prefragment mass .......
1412             IF (NPF .LT. 10) THEN
1413                  CALL EVAP(NPF, EB, EPS, NNUC, NALP)
1414 C                  EPS IS THE KINETIC ENERGY CARRIED BY THE EVAPORATED NUCLEONS
1415                ETOT = 938. + EPS
1416                  PP = SQRT((ETOT*ETOT - 8.79844E5)/3.)
1417 C                  AVERAGE MOMENTUM OF EVAPORATED NUCLEONS IN EACH DIRECTION
1418                  NUC = NPF - NNUC - 4*NALP
1419                  NF = 0
1420                  IF (NUC .GT. 0) THEN
1421                     NF = NF + 1
1422                     IAF(NF) = NUC
1423                     PPP(1,NF) = NUC*PPFX
1424                     PPP(2,NF) = NUC*PPFY
1425                  ENDIF
1426                  IF (NALP .NE. 0) THEN
1427                  DO I=1,NALP
1428                    NF = NF + 1
1429                     IAF(NF) = 4
1430                    CALL SINCO(S1,C1)
1431                    CALL SINCO(S2,C2)
1432                    PXE = 4.*PP*S1*S2
1433                    PYE = 4.*PP*S1*C2
1434                    PPP(1,NF) = 4.*PPFX + PXE
1435                    PPP(2,NF) = 4.*PPFY + PYE
1436                    PPP(1,1) = PPP(1,1) - PXE
1437                    PPP(2,1) = PPP(2,1) - PYE
1438                  ENDDO
1439                  ENDIF
1440                  IF (NNUC .NE. 0) THEN
1441                  DO I=1,NNUC
1442                     NF = NF + 1
1443                     IAF(NF) = 1
1444                     CALL SINCO(S1,C1)
1445                     CALL SINCO(S2,C2)
1446                     PXE = PP*S1*S2
1447                     PYE = PP*S1*C2
1448                     PPP(1,NF) = 4.*PPFX + PXE
1449                     PPP(2,NF) = 4.*PPFY + PYE
1450                     PPP(1,1) = PPP(1,1) - PXE
1451                     PPP(2,1) = PPP(2,1) - PYE
1452                  ENDDO
1453                  ENDIF
1454                  RETURN
1455             ENDIF
1456 
1457 C.........More refined model calculation .............
1458       JA = NPF/5 -1
1459       IF (JA .LT. 10) THEN
1460       IF ((NPF - AA(JA)) .GT. (AA(JA+1)-NPF)) JA = JA + 1
1461       ENDIF
1462       ARAT = FLOAT(NPF)/AA(JA)
1463       DO J=1,10
1464       IF (EB .LT. EAA(J)) GO TO 29
1465       ENDDO
1466       JE = 10
1467       GO TO 39
1468    29      JE = J
1469    39      IF (JE .GT. 1 .AND. JE .NE. 10) THEN
1470       IF ((EB - EAA(J-1)) .LT. (EAA(J)-EB)) JE = J - 1
1471       ENDIF
1472       ERAT = EB/EAA(JE)
1473         IF (EB .LT. 1.) THEN
1474         ERAT = EB
1475         ENDIF
1476 C INTERPOLATE BETWEEN EB=0. (NOTHING HAPPENS) AND EB = 1. MeV
1477 
1478          IF (JA .EQ. 10 .AND. JE .GT. 6) THEN
1479             WRITE(*,*)' JA=',JA,',   JE=',JE
1480          ENDIF
1481    43    ESUM = 0.
1482       NSUM = 0
1483       JF = 0
1484       DO J=20,1,-1
1485       FR =  A(JA, JE, J)*ARAT*ERAT
1486       N1 = 1 + FR
1487       FR1 = FR/FLOAT(N1)
1488       DO K=1, N1
1489       IF (S_RNDM(k) .LT. FR1) THEN
1490       JF = JF + 1
1491       IAF(JF) = J
1492       NSUM = NSUM + J
1493       EKIN = ERAT*AE(JA,JE, J)
1494          IF (EKIN .GT. 0.) THEN
1495          ESUM = ESUM + EKIN
1496          ETOT = 938.*IAF(JF) + EKIN
1497 ctp         PP = SQRT(2.*(ETOT*ETOT - IAF(JF)**2*8.79844E5)/3.)  !numerical precision problem
1498          PP = SQRT(2.*EKIN*(ETOT + IAF(JF)*938.)/3.)
1499          CALL SINCO(S1,C1)
1500          CALL SINCO(S2,C2)
1501          PPP(1,JF) = PP*S1*S2 + IAF(JF)*PPFX
1502          PPP(2,JF) = PP*S1*C2 + IAF(JF)*PPFY
1503          ENDIF
1504          IF (NSUM .GT. NPF) THEN
1505 C         WRITE(*,*)' WARNING, NSUM=', NSUM,',  NPF=',NPF
1506 C         WRITE(*,*)'  ARAT =', ARAT
1507          GO TO 43
1508         ELSE
1509         IF (NSUM .EQ. NPF) THEN
1510         GO TO 44
1511         ENDIF
1512         ENDIF
1513       ENDIF
1514       ENDDO
1515       ENDDO
1516       IF (NFLAGG(JA,JE) .EQ. 0) THEN
1517 C 'THE RESIDUE' IS A NUCLEAR FRAGMENT
1518       JF = JF + 1
1519       IAF(JF) = NPF - NSUM
1520       F1 = NPF*EB - ESUM
1521       IF (F1 .LT. 0.) F1 = 0.
1522 C GIVE THE REST OF EB TO THE FRAGMENT
1523       EKIN = F1
1524          IF (EKIN .GT. 0.) THEN
1525          ETOT = 938.*IAF(JF) + EKIN
1526 ctp         PP = SQRT(2.*(ETOT*ETOT - IAF(JF)**2*8.79844E5)/3.)  !numerical precision problem
1527          PP = SQRT(2.*EKIN*(ETOT + IAF(JF)*938.)/3.)
1528          CALL SINCO(S1,C1)
1529          CALL SINCO(S2,C2)
1530          PPP(1,JF) = PP*S1*S2 + IAF(JF)*PPFX
1531          PPP(2,JF) = PP*S1*C2 + IAF(JF)*PPFY
1532          ENDIF
1533       ELSE
1534 C 'THE RESIDUE' CONSISTS OF SPECTATOR NUCLEONS
1535       N1 = NPF - NSUM
1536       DO K=1,N1
1537       JF = JF + 1
1538       IAF(JF) = 1
1539       EKIN = ERAT*ERES(JA,JE)
1540          IF (EKIN .GT. 0.) THEN
1541          ETOT = 938.*IAF(JF) + EKIN
1542 ctp         PP = SQRT(2.*(ETOT*ETOT - IAF(JF)**2*8.79844E5)/3.)  !numerical precision problem
1543          PP = SQRT(2.*EKIN*(ETOT + IAF(JF)*938.)/3.)
1544          CALL SINCO(S1,C1)
1545          CALL SINCO(S2,C2)
1546          PPP(1,JF) = PP*S1*S2 + PPFX
1547          PPP(2,JF) = PP*S1*C2 + PPFY
1548          ENDIF
1549       ENDDO
1550       ENDIF
1551   44  NF = JF
1552       RETURN
1553       END
1554 C->
1555       FUNCTION ESTARP (NPF, NW)
1556 C CONTRIBUTION TO E* FROM ENERGY DEPOSITED BY SECONDARIES
1557 C VERY NAIVE VERSION INCORPORATING HUEFFNER'S IDEAS
1558       SAVE
1559 
1560       APF = NPF
1561       F1 = 15.3/APF**0.666666666
1562 C AVERAGE KINETIC ENERGY/NUCLEON IN PREFRAGMENT (MeV)
1563 C PER PATHLENGTH EQUAL TO THE PREFRAGMENT RADIUS
1564       ESTARP = 0.
1565       DO I=1,NW
1566       IF (S_RNDM(i) .GT. 0.5) THEN
1567 ctp060203      F2 = F1*RDIS(0)
1568       F2 = F1*RDIS()
1569       ESTARP = ESTARP + F2
1570       ENDIF
1571       ENDDO
1572 C SAMPLE RANDOMLY PER WOUNDED NUCLEON, x NW
1573       RETURN
1574       END
1575 
1576 ctp060203      function rdis(Idum)
1577       function rdis()
1578       SAVE
1579 
1580       dimension probr(20)
1581       data probr/
1582      *      0.10000, 0.15748, 0.21778, 0.28605, 0.36060,
1583      *      0.43815, 0.51892, 0.60631, 0.70002, 0.79325,
1584      *      0.88863, 0.98686, 1.10129, 1.21202, 1.32932,
1585      *      1.44890, 1.57048, 1.70139, 1.83417, 2.00000/
1586       nr = 20.*S_RNDM(0) + 1
1587       if (nr .eq. 1) then
1588       f1 = 0.
1589       else
1590       f1 = probr(nr-1)
1591       endif
1592       dr = probr(nr) - f1
1593       rdis = f1 + dr*S_RNDM(0)
1594       return
1595       end
1596 
1597 
1598       function estar(ap,at,b)
1599       implicit real*8(a-h,o-z)
1600       SAVE
1601 
1602       real*4 ap,at,b,estar
1603       sigma=4.5  !total n-n cross section in fm**2
1604       rt=.82*at**.3333 !target radius
1605       rp=.82*ap**.3333 !projectile radius
1606       alpha=rt**2/rp**2
1607       beta=b**2/rt**2
1608       f=at*sigma/(3.14159*rt**2)
1609       alf = log(f)
1610       alalf = log(alpha)
1611       gfac=0
1612       gfac1=0
1613       s1=0.
1614       s2=0.
1615       s3=0.
1616       ii=1
1617       do n=0,10 ! This limit may not need to be so high.
1618          if(n.ge.2) then
1619             gfac1=gfac
1620             gfac=gfac+log(float(n))
1621          endif
1622          g0=n*alf -n*beta*alpha/(n+alpha)+alalf
1623          g1=g0-log(alpha+n)-gfac
1624          g2=(n+2)*log(f)-(n+2)*beta*alpha/(n+2+alpha)
1625      >      +log(n+2+alpha+beta*alpha**2)-3*log(n+2+alpha)-gfac
1626          g3=g0-2*log(n+alpha)-gfac1
1627          ii=-ii
1628          s1=s1+ii*exp(g1)
1629          s2=s2+ii*exp(g2)
1630          if(n.ge.1) s3=s3+ii*exp(g3)
1631       enddo
1632 
1633       pb=s1
1634       e1b=197.**2/(2*938.*rp**2*pb) *s2
1635 c      a=b*(s3/pb-1)
1636 c      a=-b*s3/pb
1637 c      e2b=-.5* 938. * (41./(ap**.333))**2 * a**2 /(197.**2)
1638 c      estar=e1b+e2b
1639       estar = e1b
1640       return
1641       end
1642 
1643       subroutine evap(npf,eb,eps,nnuc,nalp)
1644       SAVE
1645 
1646       eps=7.5+sqrt(8*eb)
1647       n=min(npf*int(eb/eps),npf)
1648       nalp=n/5
1649       nnuc=n-4*nalp
1650       return
1651       end
1652 C->
1653       FUNCTION FERMK(A)
1654       SAVE
1655 
1656       DIMENSION AA(6), FK(6)
1657       DATA AA/4., 6., 12., 24., 40., 57./
1658       DATA FK/130.,169.,221.,235.,251.,260./
1659       DO I=2,4
1660       IF (A .LT. AA(I)) GO TO 25
1661       ENDDO
1662       I = 5
1663    25      F11 = AA(I-1)
1664       F12 = AA(I)
1665       F13 = AA(I+1)
1666       F21 = FK(I-1)
1667       F22 = FK(I)
1668       F23 = FK(I+1)
1669       FERMK = QUAD_INT(A,F11,F12,F13, F21,F22,F23)
1670       RETURN
1671       END
1672 
1673 C========================================================================
1674 C. Multiple interaction structure
1675 C========================================================================
1676 
1677       SUBROUTINE INT_NUC (IA, IB, SIG0, SIGEL)
1678 C...Compute with a montecarlo code  the  "multiple interaction structure"
1679 C.  of a nucleus-nucleus interaction
1680 C.
1681 C.  INPUT : IA            = mass of target nucleus
1682 C.          IB            = mass of projectile nucleus
1683 C.          SIG0 (mbarn)  = inelastic pp cross section
1684 C.          SIGEL(mbarn)  = elastic pp cross section
1685 C.
1686 C.  OUTPUT : in common block /CNUCMS/
1687 C.           B = impact parameter (fm)
1688 C.           BMAX = maximum impact parameter for generation
1689 C.           NTRY = number of "trials" before one interaction
1690 C.           NA = number of wounded nucleons in A
1691 C.           NB =    "        "        "     in B
1692 C.           NI = number of nucleon-nucleon inelastic interactions
1693 C.           NAEL = number of elastically scattered nucleons in  A
1694 C.           NBEL =    "         "           "          "    in  B
1695 C.           JJA(J)  [J=1:IA]   = number of inelastic interactions
1696 C.                                of J-th nucleon of nucleus A
1697 C.           JJB(J)  [J=1:IB]   = number of inelastic interactions
1698 C.                                of J-th nucleon of nucleus B
1699 C.           JJAEL(J)  [J=1:IA]   = number of elastic interactions
1700 C.                                of J-th nucleon of nucleus A
1701 C.           JJBEL(J)  [J=1:IB]   = number of elastic interactions
1702 C.                                of J-th nucleon of nucleus B
1703 C.           JJINT(J,K)  [J=1:NB, K=1:NA]  (0 = no interaction)
1704 C.                                         (1 = interaction )
1705 C.                                         between nucleon J of A and K of B
1706 C-----------------------------------------------------------------------------
1707       SAVE
1708 
1709       PARAMETER (IAMAX=56)
1710       COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
1711      +         ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
1712      +         ,JJAEL(IAMAX), JJBEL(IAMAX)
1713       DIMENSION XA(IAMAX), YA(IAMAX), XB(IAMAX), YB(IAMAX)
1714       DATA PI /3.1415926/
1715       SIGT = SIG0 + SIGEL
1716       R2  = 0.1 * SIG0/PI
1717       R2T = 0.1 * SIGT/PI
1718       BMAX = 15.                             ! fm
1719       NTRY = 0
1720       CALL NUC_CONF (IA, XA, YA)
1721       CALL NUC_CONF (IB, XB, YB)
1722       NI = 0
1723       NIEL = 0
1724       DO JA=1,IA
1725          JJA(JA) = 0
1726          JJAEL(JA) = 0
1727       ENDDO
1728       DO JB=1,IB
1729          JJB(JB) = 0
1730          JJBEL(JB) = 0
1731          DO JA=1,IA
1732             JJINT(JB,JA) = 0
1733          ENDDO
1734       ENDDO
1735 1000  B = BMAX*SQRT(S_RNDM(0))
1736       PHI = 2.*PI*S_RNDM(0)
1737       BX = B*COS(PHI)
1738       BY = B*SIN(PHI)
1739       NTRY = NTRY+1
1740       DO JA=1,IA
1741          DO JB=1,IB
1742             S = (XA(JA)-XB(JB)-BX)**2 + (YA(JA)-YB(JB)-BY)**2
1743             IF (S .LT. R2)  THEN
1744                NI = NI + 1
1745                JJA(JA) = JJA(JA)+1
1746                JJB(JB) = JJB(JB)+1
1747                JJINT(JB,JA) = 1
1748             ELSE IF (S .LT. R2T)  THEN
1749                NIEL = NIEL + 1
1750                JJAEL(JA) = JJAEL(JA)+1
1751                JJBEL(JB) = JJBEL(JB)+1
1752             ENDIF
1753          ENDDO
1754       ENDDO
1755       IF (NI + NIEL .EQ. 0)  GOTO 1000
1756       NA = 0
1757       NB = 0
1758       NAEL = 0
1759       NBEL = 0
1760       DO JA=1,IA
1761          IF (JJA(JA) .GT. 0)  THEN
1762             NA = NA + 1
1763          ELSE
1764             IF (JJAEL(JA) .GT. 0)  NAEL = NAEL+1
1765          ENDIF
1766       ENDDO
1767       DO JB=1,IB
1768          IF (JJB(JB) .GT. 0)  THEN
1769             NB = NB + 1
1770          ELSE
1771             IF (JJBEL(JB) .GT. 0)  NBEL = NBEL+1
1772          ENDIF
1773       ENDDO
1774       RETURN
1775       END
1776 
1777        SUBROUTINE NUC_CONF (IA, XX, YY)
1778 C...This routine generates the configuration  of a nucleus
1779 C.  need an initialization call to NUC_GEOM_INI
1780 C.
1781 C.  INPUT  : IA = mass number of the nucleus
1782 C.  OUTPUT : XX(1:IA), YY(1:IA) (fm) = position in impact parameter
1783 C.                                     space of the IA nucleons
1784 C...................................................................
1785       SAVE
1786 
1787       PARAMETER (IAMAX=56)
1788       DIMENSION XX(IAMAX), YY(IAMAX)
1789       PARAMETER (NB=401)
1790       COMMON /CPROFA/ ZMIN, DZ, BBZ(NB,IAMAX)
1791       DATA PI /3.1415926/
1792       DO J=1,IA
1793          Z = S_RNDM(j)
1794          JZ = INT((Z-ZMIN)/DZ)+1
1795 CDH
1796          JZ = MIN(JZ,400)
1797          T = (Z-ZMIN)/DZ - FLOAT(JZ-1)
1798          B = BBZ(JZ,IA)*(1.-T) + BBZ(JZ+1,IA)*T
1799          PHI = 2.*PI*S_RNDM(jz)
1800          XX(J) = B*COS(PHI)
1801          YY(J) = B*SIN(PHI)
1802       ENDDO
1803       RETURN
1804       END
1805 
1806       SUBROUTINE NUC_GEOM_INI
1807 C...Initialize all nucleus profiles
1808       SAVE
1809 
1810       PARAMETER (NB=401)
1811       PARAMETER (IAMAX=56)
1812       COMMON /CPROF/ DB, BMAX, BB(NB), TB(NB), A
1813       COMMON /CPROFA/ ZMIN, DZ, BBZ(NB,IAMAX)
1814       DIMENSION FFB(NB), GGB(NB)
1815       DATA PI /3.1415926/
1816       CALL SHELL_INI
1817       CALL WOOD_SAXON_INI
1818       DO IA= 2,IAMAX
1819            JA = IA
1820          CALL NUC_PROFIL(JA)
1821          DO K=1,NB
1822            FFB(K) = BB(K)*TB(K) * (2.*PI)
1823          ENDDO
1824          GGB(1) = 0.
1825          GGB(NB) = 1.
1826          DO K=2,NB-1
1827            GGB(K) = GGB(K-1) + FFB(K-1)*DB
1828          ENDDO
1829          CALL INVERT_ARRAY(GGB,0.,DB,NB, BBZ(1,IA), ZMIN, DZ)
1830       ENDDO
1831       RETURN
1832       END
1833 
1834       SUBROUTINE NUC_PROFIL (JA)
1835 C...Compute the profile function T(b)
1836 C.  normalised as INT[d2b T(b) = 1]
1837 C.  INPUT : JA = integer mass number of nucleus
1838 C...............................................
1839       SAVE
1840 
1841       PARAMETER (NB=401)
1842       EXTERNAL DENSA
1843       REAL DENSA
1844       COMMON /CC01/  B
1845       COMMON /CCDA/ JJA
1846       COMMON /CPROF/ DB, BMAX, BB(NB), TB(NB), A
1847       BMAX = 7.5
1848       DB = BMAX/FLOAT(NB-1)
1849       JJA = JA
1850       A = JA
1851       DO JB=1,NB
1852         B = DB*FLOAT(JB-1)
1853         BB(JB) = B
1854         IF (JA .LE. 18)  THEN
1855             TB(JB) = PROFNUC (B, JA)
1856          ELSE
1857             TB(JB) = 2.*GAUSS (DENSA,0.,BMAX)
1858          ENDIF
1859       ENDDO
1860       RETURN
1861       END
1862 
1863       SUBROUTINE NUC1_PROFIL (AA)
1864 C...Compute the profile function T(b)
1865 C.  normalised as INT[d2b T(b) = 1]
1866 C.  INPUT : AA = mass number of nucleus
1867 C...............................................
1868       SAVE
1869 
1870       PARAMETER (NB=401)
1871       EXTERNAL DENSA
1872       REAL DENSA
1873       COMMON /CC01/  B
1874       COMMON /CPROF/ DB, BMAX, BB(NB), TB(NB), A
1875       A = AA
1876       IA1 = INT(AA)
1877       IA2 = IA1 + 1
1878       U = AA - FLOAT(IA1)
1879       BMAX = 7.5
1880       DB = BMAX/FLOAT(NB-1)
1881       DO JB=1,NB
1882          B = DB*FLOAT(JB-1)
1883          BB(JB) = B
1884          IF (A .LE. 18.)  THEN
1885              T1 = PROFNUC (B, IA1)
1886              T2 = PROFNUC (B, IA2)
1887           ELSE
1888              JJA = IA1
1889              T1 = 2.*GAUSS (DENSA,0.,BMAX)
1890              JJA = IA2
1891              T2 = 2.*GAUSS (DENSA,0.,BMAX)
1892           ENDIF
1893           TB(JB) = (1.-U)*T1  + U*T2
1894       ENDDO
1895       RETURN
1896       END
1897 
1898 C===========================================================================
1899 C.   Code about nuclear densities
1900 C===========================================================================
1901 
1902       FUNCTION DENS_NUC (R, JA)
1903 C....Nuclear density (normalised to 1)
1904 C.   for a nucleus of mass number JA
1905 C.   INPUT R = radial coordinate  (fm)
1906 C.         JA = integer mass number
1907 C.  OUTPUT (fm**-3)
1908 C--------------------------------------------------------
1909       SAVE
1910 
1911       COMMON /CWOOD/ RR0(19:56), AA0(19:56), CC0(19:56)
1912       IF (JA .GT. 18)  THEN
1913          DENS_NUC = WOOD_SAXON(R,JA)
1914       ELSE IF (JA .NE. 4)  THEN
1915          DENS_NUC = HELIUM(R)
1916       ELSE
1917          DENS_NUC = SHELL(R,JA)
1918       ENDIF
1919       RETURN
1920       END
1921 
1922       FUNCTION WOOD_SAXON (R, JA)
1923 C....Wood-Saxon nuclear density (normalised to 1)
1924 C.   for a nucleus of mass number A.
1925 C.   INPUT R =  (fm)
1926 C.         JA = mass number
1927 C.   OUTPUT (fm**-3)
1928 C------------------------------------------------------
1929       SAVE
1930 
1931       COMMON /CWOOD/ RR0(19:56), AA0(19:56), CC0(19:56)
1932       WOOD_SAXON = CC0(JA)/(1.+EXP((R-RR0(JA))/AA0(JA)))
1933       RETURN
1934       END
1935 
1936       FUNCTION HELIUM (R)
1937 C... Helium density from Barrett and Jackson
1938 C.   INPUT R = r coordinate (fm)
1939 C.   OUTPUT (fm**-3)
1940 C........................................................
1941       SAVE
1942 
1943       DATA R0 /0.964/, CA /0.322/   ! fm
1944       DATA W /0.517/, CC /5.993224E-02/
1945       HELIUM = CC*(1.+W*(R/R0)**2)/(1. + EXP((R-R0)/CA))
1946       RETURN
1947       END
1948 
1949       FUNCTION SHELL (R,JA)
1950 C...Density in the shell model
1951       COMMON /CSHELL/ RR0(18), RR02(18)
1952       SAVE
1953 
1954       DATA PI /3.1415926/
1955       R0 = RR0(JA)
1956       C1 = MIN(1.,4./FLOAT(JA))
1957       CS = 1./(R0**3*PI**(1.5))
1958       CP = 2.*CS/3.
1959       FS = EXP(-(R/R0)**2)
1960       FP = (R/R0)**2 * FS
1961       SHELL = C1*CS*FS + (1.-C1)*CP*FP
1962       RETURN
1963       END
1964 
1965       FUNCTION PROFNUC (B, JA)
1966 C...This function return
1967 C.  the profile T(b) for a nucleus of mass number A
1968 C.  INPUT B = impact parameter (GeV**-1)
1969 C.        JA = integer mass number
1970 C.  OUTPUT  (fm**-2)
1971 C.
1972 C.  The  density of the nucleus is the `shell model density'
1973 C.  the parameter r0 must beinitialized in the common block
1974 C.............................................................
1975       SAVE
1976 
1977       COMMON /CSHELL/ RR0(18), RR02(18)
1978       DATA PI /3.1415926/
1979       B2 = B*B
1980       ARG = B2/RR02(JA)
1981       TS = EXP(-ARG)
1982       TP = TS*(2.*B2+RR02(JA))/(3.*RR02(JA))
1983       CS = MIN(1.,4./FLOAT(JA))
1984       PROFNUC = (CS*TS + (1.-CS)*TP)/(PI*RR02(JA))
1985       RETURN
1986       END
1987 
1988       SUBROUTINE SHELL_INI
1989 C...Initialize the parameter  of the shell model
1990 C.  for the nuclei with    6 < A < 18
1991 C..............................................
1992       SAVE
1993 
1994       COMMON /CSHELL/ RR0(18), RR02(18)
1995       DIMENSION RR(18)
1996 C...Data on Sqrt[<r**2>]  in fermi
1997       DATA RR /0.81,2.095,1.88,1.674, -1.,2.56,2.41,-1.,2.519,2.45
1998      +          ,2.37, 2.460, 2.440, 2.54, 2.58, 2.718, 2.662,2.789 /
1999       DO JA=1,18
2000          A = FLOAT(JA)
2001          RMED = RR(JA)
2002          IF (RMED .LE. 0.)   RMED = 0.5*(RR(JA-1) + RR(JA+1))
2003          C = MAX(1.5,(5./2. - 4./A) )
2004          R0 = RMED/SQRT(C)
2005          RR0 (JA) = R0
2006          RR02(JA) = R0*R0
2007       ENDDO
2008       RETURN
2009       END
2010 C->
2011       SUBROUTINE WOOD_SAXON_INI
2012       COMMON /CWOOD/ RR0(19:56), AA0(19:56), CC0(19:56)
2013       SAVE
2014 
2015       DATA PI /3.1415926/
2016 C...Wood-Saxon parameters from  table 6.2   of Barrett and Jackson
2017       RR0 (19) = 2.59
2018       AA0 (19) = 0.564
2019       RR0 (20) = 2.74
2020       AA0 (20) = 0.569
2021       RR0 (22) = 2.782
2022       AA0 (22) = 0.549
2023       RR0 (24) = 2.99
2024       AA0 (24) = 0.548
2025       RR0 (27) = 2.84
2026       AA0 (27) = 0.569
2027       RR0 (28) = 3.14
2028       AA0 (28) = 0.537
2029       RR0 (29) = 3.77
2030       AA0 (29) = 0.52
2031       RR0 (48) = 3.912
2032       AA0 (48) = 0.5234
2033       RR0 (56) = 3.98
2034       AA0 (56) = 0.569
2035       DO J=19, 56
2036          IF (RR0(J) .LE. 0.)  THEN
2037             RR0(J) = 1.05*FLOAT(J)**0.333333
2038             AA0(J) = 0.545
2039          ENDIF
2040          CC0(J)=3./(4.*PI*RR0(J)**3)/(1.+((AA0(J)*PI)/RR0(J))**2)
2041       ENDDO
2042       RETURN
2043       END
2044 
2045       FUNCTION DENSA (Z)
2046 C....Woods Saxon nuclear density (normalised to 1)
2047 C.   for a nucleus of mass number A.
2048 C.   INPUT z = z coordinate (fm)
2049 C.         JA = integer mass number
2050 C.         B (in common /CC01/)  impact parameter  (fm)
2051 C.  OUTPUT (fm**-3)
2052 C--------------------------------------------------------
2053       SAVE
2054 
2055       COMMON /CC01/  B
2056       COMMON /CCDA/ JA
2057       COMMON /CWOOD/ RR0(19:56), AA0(19:56), CC0(19:56)
2058       R = SQRT (Z*Z + B*B)
2059       DENSA = CC0(JA)/(1.+EXP((R-RR0(JA))/AA0(JA)))
2060       RETURN
2061       END
2062 
2063 C==========================================================================
2064 C. Cross sections
2065 C==========================================================================
2066 
2067       SUBROUTINE SIGMA_AIR (IB,SIG0,SIGEL,NINT,
2068      +                            SIGMA,DSIGMA,SIGQE,DSIGQE)
2069 C...Compute with a montecarlo method the "production"
2070 C.  and "quasi-elastic" cross section for
2071 C.  a nucleus-air  interaction
2072 C.
2073 C.  INPUT : IB            = mass of projectile nucleus
2074 C.          SIG0 (mbarn)  = inelastic pp cross section
2075 C.          NINT            = number  of interactions to generate
2076 C.  OUTPUT : SIGMA (mbarn) = "production" cross section
2077 C.           DSIGMA   "    = error
2078 C.           SIGQE    "    = "quasi-elastic" cross section
2079 C.           DSIGQE   "    = error
2080 C.           additional output is in the common block  /CPROBAB/
2081 C..........................................................................
2082       SAVE
2083 
2084       PARAMETER (IAMAX=56)
2085       PARAMETER (IAMAX2=3136)          ! IAMAX*IAMAX
2086       COMMON  /CPROBAB/ PROBA(IAMAX), DPROBA(IAMAX),
2087      +   PROBB(IAMAX), DPROBB(IAMAX), PROBI(IAMAX2), DPROBI(IAMAX2),
2088      +   P1AEL(0:IAMAX),DP1AEL(0:IAMAX),P1BEL(0:IAMAX), DP1BEL(0:IAMAX),
2089      +   P2AEL(0:IAMAX),DP2AEL(0:IAMAX),P2BEL(0:IAMAX), DP2BEL(0:IAMAX)
2090       COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
2091      +         ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
2092      +         ,JJAEL(IAMAX), JJBEL(IAMAX)
2093       DIMENSION  MMA(0:IAMAX), MMB(0:IAMAX), MMI(0:IAMAX2)
2094       DIMENSION  M1AEL(0:IAMAX), M1BEL(0:IAMAX)
2095       DIMENSION  M2AEL(0:IAMAX), M2BEL(0:IAMAX)
2096       DATA WOX /0.346/
2097       DATA PI /3.1415926/
2098       R2 = 0.1 * SIG0/PI
2099       BMAX = 15.                             ! fm
2100       SIGMA0 = PI*BMAX*BMAX*10.              ! mbarn
2101       IA = 16
2102       DO J=1,IA
2103          MMA(J) = 0
2104          M1AEL(J) = 0
2105          M2AEL(J) = 0
2106       ENDDO
2107       DO J=1,IB
2108          MMB(J) = 0
2109          M1BEL(J) = 0
2110          M2BEL(J) = 0
2111       ENDDO
2112       DO J=1,IA*IB
2113          MMI(J) = 0
2114       ENDDO
2115       NN = 0
2116       M = 0
2117       DO KK=1,NINT
2118          IA = 14 + 2*INT((1.+WOX)*S_RNDM(kk))
2119          CALL INT_NUC (IA, IB, SIG0, SIGEL)
2120          NN = NN + NTRY
2121          MMI(NI) = MMI(NI) + 1
2122          MMA(NA) = MMA(NA)+1
2123          MMB(NB) = MMB(NB)+1
2124          IF (NI .GT. 0)  THEN
2125             M = M+1
2126             M1AEL(NAEL) = M1AEL(NAEL)+1
2127             M1BEL(NBEL) = M1BEL(NBEL)+1
2128          ELSE
2129             M2AEL(NAEL) = M2AEL(NAEL)+1
2130             M2BEL(NBEL) = M2BEL(NBEL)+1
2131          ENDIF
2132       ENDDO
2133       MQE = NINT - M
2134       SIGMA  = SIGMA0 * FLOAT(M)/FLOAT(NN)
2135       DSIGMA = SIGMA0 * SQRT(FLOAT(M))/FLOAT(NN)
2136       SIGQE  = SIGMA0 * FLOAT(MQE)/FLOAT(NN)
2137       DSIGQE = SIGMA0 * SQRT(FLOAT(MQE))/FLOAT(NN)
2138       DO J=1,IA
2139          PROBA(J) = FLOAT(MMA(J))/FLOAT(M)
2140          DPROBA(J) = SQRT(FLOAT(MMA(J)))/FLOAT(M)
2141       ENDDO
2142       DO J=1,IB
2143          PROBB(J) = FLOAT(MMB(J))/FLOAT(M)
2144          DPROBB(J) = SQRT(FLOAT(MMB(J)))/FLOAT(M)
2145       ENDDO
2146       DO J=1,IA*IB
2147          PROBI(J) = FLOAT(MMI(J))/FLOAT(M)
2148          DPROBI(J) = SQRT(FLOAT(MMI(J)))/FLOAT(M)
2149       ENDDO
2150       DO J=0,IA
2151          P1AEL(J) = FLOAT(M1AEL(J))/FLOAT(M)
2152          DP1AEL(J) = SQRT(FLOAT(M1AEL(J)))/FLOAT(M)
2153          P2AEL(J) = FLOAT(M2AEL(J))/FLOAT(MQE)
2154          DP2AEL(J) = SQRT(FLOAT(M2AEL(J)))/FLOAT(MQE)
2155       ENDDO
2156       DO J=0,IB
2157          P1BEL(J) = FLOAT(M1BEL(J))/FLOAT(M)
2158          DP1BEL(J) = SQRT(FLOAT(M1BEL(J)))/FLOAT(M)
2159          P2BEL(J) = FLOAT(M2BEL(J))/FLOAT(MQE)
2160          DP2BEL(J) = SQRT(FLOAT(M2BEL(J)))/FLOAT(MQE)
2161       ENDDO
2162       RETURN
2163       END
2164 C->
2165       SUBROUTINE SIGMA_MC (IA,IB,SIG0,SIGEL,NINT,
2166      +                            SIGMA,DSIGMA,SIGQE,DSIGQE)
2167 C...Compute with a montecarlo method the "production"
2168 C.  and "quasi-elastic" cross section for
2169 C.  a nucleus-nucleus interaction
2170 C.
2171 C.  INPUT : IA            = mass of target nucleus
2172 C.          IB            = mass of projectile nucleus
2173 C.          SIG0 (mbarn)  = inelastic pp cross section
2174 C.          NINT            = number  of interactions to generate
2175 C.  OUTPUT : SIGMA (mbarn) = "production" cross section
2176 C.           DSIGMA   "    = error
2177 C.           SIGQE    "    = "quasi-elastic" cross section
2178 C.           DSIGQE   "    = error
2179 C.           additional output is in the common block  /CPROBAB/
2180 C.           Prob(n_A), Prob(n_B), Prob(n_int)
2181 C..........................................................................
2182       SAVE
2183 
2184       PARAMETER (IAMAX=56)
2185       PARAMETER (IAMAX2=3136)          ! IAMAX*IAMAX
2186       COMMON  /CPROBAB/ PROBA(IAMAX), DPROBA(IAMAX),
2187      +   PROBB(IAMAX), DPROBB(IAMAX), PROBI(IAMAX2), DPROBI(IAMAX2),
2188      +   P1AEL(0:IAMAX),DP1AEL(0:IAMAX),P1BEL(0:IAMAX), DP1BEL(0:IAMAX),
2189      +   P2AEL(0:IAMAX),DP2AEL(0:IAMAX),P2BEL(0:IAMAX), DP2BEL(0:IAMAX)
2190       COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
2191      +         ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
2192      +         ,JJAEL(IAMAX), JJBEL(IAMAX)
2193       DIMENSION  MMA(0:IAMAX), MMB(0:IAMAX), MMI(0:IAMAX2)
2194       DIMENSION  M1AEL(0:IAMAX), M1BEL(0:IAMAX)
2195       DIMENSION  M2AEL(0:IAMAX), M2BEL(0:IAMAX)
2196       DATA PI /3.1415926/
2197       R2 = 0.1 * SIG0/PI
2198       BMAX = 15.                             ! fm
2199       SIGMA0 = PI*BMAX*BMAX*10.              ! mbarn
2200       DO J=1,IA
2201          MMA(J) = 0
2202          M1AEL(J) = 0
2203          M2AEL(J) = 0
2204       ENDDO
2205       DO J=1,IB
2206          MMB(J) = 0
2207          M1BEL(J) = 0
2208          M2BEL(J) = 0
2209       ENDDO
2210       DO J=1,IA*IB
2211          MMI(J) = 0
2212       ENDDO
2213       NN = 0
2214       M = 0
2215       DO KK=1,NINT
2216          CALL INT_NUC (IA, IB, SIG0, SIGEL)
2217          NN = NN + NTRY
2218          MMI(NI) = MMI(NI) + 1
2219          MMA(NA) = MMA(NA)+1
2220          MMB(NB) = MMB(NB)+1
2221          IF (NI .GT. 0)  THEN
2222             M = M+1
2223             M1AEL(NAEL) = M1AEL(NAEL)+1
2224             M1BEL(NBEL) = M1BEL(NBEL)+1
2225          ELSE
2226             M2AEL(NAEL) = M2AEL(NAEL)+1
2227             M2BEL(NBEL) = M2BEL(NBEL)+1
2228          ENDIF
2229       ENDDO
2230       MQE = NINT - M
2231       SIGMA  = SIGMA0 * FLOAT(M)/FLOAT(NN)
2232       DSIGMA = SIGMA0 * SQRT(FLOAT(M))/FLOAT(NN)
2233       SIGQE  = SIGMA0 * FLOAT(MQE)/FLOAT(NN)
2234       DSIGQE = SIGMA0 * SQRT(FLOAT(MQE))/FLOAT(NN)
2235       DO J=1,IA
2236          PROBA(J) = FLOAT(MMA(J))/FLOAT(M)
2237          DPROBA(J) = SQRT(FLOAT(MMA(J)))/FLOAT(M)
2238       ENDDO
2239       DO J=1,IB
2240          PROBB(J) = FLOAT(MMB(J))/FLOAT(M)
2241          DPROBB(J) = SQRT(FLOAT(MMB(J)))/FLOAT(M)
2242       ENDDO
2243       DO J=1,IA*IB
2244          PROBI(J) = FLOAT(MMI(J))/FLOAT(M)
2245          DPROBI(J) = SQRT(FLOAT(MMI(J)))/FLOAT(M)
2246       ENDDO
2247       DO J=0,IA
2248          P1AEL(J) = FLOAT(M1AEL(J))/FLOAT(M)
2249          DP1AEL(J) = SQRT(FLOAT(M1AEL(J)))/FLOAT(M)
2250          P2AEL(J) = FLOAT(M2AEL(J))/FLOAT(MQE)
2251          DP2AEL(J) = SQRT(FLOAT(M2AEL(J)))/FLOAT(MQE)
2252       ENDDO
2253       DO J=0,IB
2254          P1BEL(J) = FLOAT(M1BEL(J))/FLOAT(M)
2255          DP1BEL(J) = SQRT(FLOAT(M1BEL(J)))/FLOAT(M)
2256          P2BEL(J) = FLOAT(M2BEL(J))/FLOAT(MQE)
2257          DP2BEL(J) = SQRT(FLOAT(M2BEL(J)))/FLOAT(MQE)
2258       ENDDO
2259       RETURN
2260       END
2261 
2262 C=============================================================
2263 C.  Cross sections
2264 C=============================================================
2265 
2266 
2267       SUBROUTINE SIG_H_AIR (SSIG, SLOPE, ALPHA,  SIGT, SIGEL, SIGQE)
2268 C...Subroutine to compute hadron-air cross sections
2269 C.  according to:
2270 C.  R.J. Glauber and G.Matthiae  Nucl.Phys. B21, 135, (1970)
2271 C.
2272 C.  Air is a linear combination of Nitrogen and oxygen
2273 C.
2274 C.  INPUT :  SSIG  (mbarn) total pp cross section
2275 C.           SLOPE (GeV**-2)  elastic scattering slope for pp
2276 C.           ALPHA    real/imaginary part of the forward pp elastic
2277 C.                                               scattering amplitude
2278 C.  OUTPUT : SIGT  = Total cross section
2279 C.           SIGEL = Elastic cross section
2280 C.           SIGQEL  = Elastic + Quasi elastic cross section
2281 C......................................................................
2282       SAVE
2283 
2284       DATA  FOX /0.257/
2285       CALL GLAUBER(14,SSIG,SLOPE,ALPHA,SIG1,SIGEL1,SIGQE1)
2286       CALL GLAUBER(16,SSIG,SLOPE,ALPHA,SIG2,SIGEL2,SIGQE2)
2287       SIGT  = (1.-FOX)*SIG1   + FOX*SIG2
2288       SIGEL = (1.-FOX)*SIGEL1 + FOX*SIGEL2
2289       SIGQE = (1.-FOX)*SIGQE1 + FOX*SIGQE2
2290       RETURN
2291       END
2292 
2293       SUBROUTINE GLAUBER(JA,SSIG,SLOPE,ALPHA,SIGT,SIGEL,SIGQEL)
2294 C...Subroutine to compute hadron-Nucleus cross sections
2295 C.  according to:
2296 C.  R.J. Glauber and G.Matthiae  Nucl.Phys. B21, 135, (1970)
2297 C.
2298 C.  This formulas assume that the target nucleus  density is
2299 C.  modeled by a shell-model form.  A reasonable range of models
2300 C   is  4 < JA < 18
2301 C.
2302 C.  INPUT :  A = mass number of the nucleus
2303 C.           SSIG  (mbarn) total pp cross section
2304 C.           SLOPE (GeV**-2)  elastic scattering slope for pp
2305 C.           ALPHA    real/imaginary part of the forward pp elastic
2306 C.                                               scattering amplitude
2307 C.  OUTPUT : SIGT  = Total cross section
2308 C.           SIGEL = Elastic cross section
2309 C.           SIGQEL  = Elastic + Quasi elastic cross section
2310 C.
2311 C. Internally  everything is computed in GeV (length = GeV**-1)
2312 C......................................................................
2313       SAVE
2314 
2315       COMMON /CA0SH/ R0, R02
2316       COMPLEX  ZZ, ZS, ZP, ZC
2317       DIMENSION RR(18)
2318       DATA CMBARN /0.389385/
2319       DATA PI /3.1415926/
2320       DATA BMAX /50./            ! GeV**-1
2321       DATA NB /100/
2322 C...data on Sqrt[<r**2>] (fm). (A=5,8 are not correct).
2323 C   from Barett and Jackson
2324       DATA RR /0.81,2.095,1.88,1.674, 2.56,2.56,2.41,2.5,2.519,2.45
2325      +          ,2.37, 2.460, 2.440, 2.54, 2.58, 2.718, 2.662,2.789 /
2326       A = FLOAT(JA)
2327 C...Parameter of shell model density
2328       R0 = RR(JA)/0.197/SQRT(5./2. - 4./A)         ! GeV**-1
2329       R02 = R0*R0
2330       SIG = SSIG/CMBARN                           ! GeV**-2
2331       DB = BMAX/FLOAT(NB)
2332       SUM = 0.
2333       SUM1 = 0.
2334       SUM2 = 0.
2335       DO JB=1,NB
2336          B = DB*(FLOAT(JB)-0.5)
2337          GS = GLAUBGS (B,SLOPE, SIG)
2338          GP = GLAUBGP (B,SLOPE, SIG)
2339          XS = (1.- GS)
2340          YS =  GS*ALPHA
2341          ZS = CMPLX(XS,YS)
2342          XP = (1.- GP)
2343          YP =  GP*ALPHA
2344          ZP = CMPLX(XP,YP)
2345          ZZ = ZS**4. * ZP**(A-4.)
2346          X = REAL (ZZ)
2347          Y = AIMAG(ZZ)
2348          ZC = CMPLX(X,-Y)
2349          SUM = SUM + (1.-X)*B
2350          SUM1 = SUM1 + ((1.-X)**2 + Y**2)*B
2351          OMS = OMEGAS(B,SIG,SLOPE,ALPHA)
2352          OMP = OMEGAP(B,SIG,SLOPE,ALPHA)
2353          OM = (1.- 2.*GS + OMS)**4. * (1. -2.*GP + OMP)**(A-4.)
2354          SUM2 = SUM2 + (1.-2.*X + OM)*B
2355       ENDDO
2356       SIGT =   SUM  * DB * 4.*PI * CMBARN
2357       SIGEL =  SUM1 * DB * 2.*PI * CMBARN
2358       SIGQEL = SUM2 * DB * 2.*PI * CMBARN
2359       RETURN
2360       END
2361 C->
2362       FUNCTION GLAUBGS (B,SLOPE, SIG)
2363       SAVE
2364 
2365       COMMON /CA0SH/ A0, A02
2366       DATA PI /3.1415926/
2367       GAMMA2 = A02/4. + 0.5*SLOPE
2368       ARG = B**2/(4.*GAMMA2)
2369       GLAUBGS = SIG/(8.*PI*GAMMA2) * EXP(-ARG)
2370       RETURN
2371       END
2372 C->
2373       FUNCTION GLAUBGP (B,SLOPE, SIG)
2374       SAVE
2375 
2376       COMMON /CA0SH/ A0, A02
2377       DATA PI /3.1415926/
2378       GAMMA2 = A02/4. + 0.5*SLOPE
2379       ARG = B**2/(4.*GAMMA2)
2380       C1 = 1.- A02/(6.*GAMMA2)*(1.-ARG)
2381       GLAUBGP = SIG/(8.*PI*GAMMA2) *  C1 * EXP(-ARG)
2382       RETURN
2383       END
2384 C->
2385       FUNCTION OMEGAS (B, SIG, SLOPE, RHO)
2386       SAVE
2387 
2388       COMMON /CA0SH/ A0, A02
2389       DATA PI /3.1415926/
2390       ETA2 = 0.25*(A02 + SLOPE)
2391       F02 = SIG*SIG*(1.+RHO*RHO)/(16.*PI**2)
2392       ARG = -B*B/(4.*ETA2)
2393       OMEGAS = F02/(4.*ETA2*SLOPE) *EXP(ARG)
2394       RETURN
2395       END
2396 C->
2397       FUNCTION OMEGAP (B, SIG, SLOPE, RHO)
2398       SAVE
2399 
2400       COMMON /CA0SH/ A0, A02
2401       DATA PI /3.1415926/
2402       ETA2 = 0.25*(A02 + SLOPE)
2403       F02 = SIG*SIG*(1.+RHO*RHO)/(16.*PI**2)
2404       ARG = -B*B/(4.*ETA2)
2405       OMEGAP=F02/(4.*ETA2*SLOPE)*(1.-A02/(6.*ETA2)*(1.+ARG))*EXP(ARG)
2406       RETURN
2407       END
2408 
2409 C------------------------------------------------------------------------
2410 C.  Fit of Block and Cahn to pp and pbar-p cross sections
2411 C------------------------------------------------------------------------
2412 
2413       SUBROUTINE BLOCK(SQS,SIG1,SIG2,SLOP1,SLOP2,
2414      +                 RHO1,RHO2,SIGEL1,SIGEL2)
2415 C...p-p and pbar-p cross sections
2416 C.  Parametrization of  Block and Cahn
2417 C
2418 C.  INPUT  : SQS   (GeV)  = c.m. energy
2419 C.
2420 C.  OUPUT : SIG1 (mbarn)    = pp  total  cross section
2421 C.          SLOP1 (GeV**2)  = slope of elastic scattering
2422 C.          RHO1            = Real/Imaginary part of the amplitude
2423 C.                            for forward elastic  scattering (pp)
2424 C.          SIGEL1 (mbarn)  = pp  elastic scattering  cross section
2425 C.          [1 -> 2   : pp -> pbar p]
2426 C-----------------------------------------------------------------------
2427       SAVE
2428 
2429       DATA PI /3.1415926/
2430       DATA CMBARN /0.389385/
2431       S = SQS*SQS
2432       CALL FPLUS  (S, FR, FI)
2433       CALL FMINUS (S, GR, GI)
2434       SIG1 = FI-GI
2435       SIG2 = FI+GI
2436       RHO1 = (FR-GR)/(FI-GI)
2437       RHO2 = (FR+GR)/(FI+GI)
2438       CALL SSLOPE (S, BP, BM)
2439       SLOP1 = BP - GI/FI*(BM-BP)
2440       SLOP2 = BP + GI/FI*(BM-BP)
2441       SIGEL1 = SIG1**2*(1.+RHO1**2)/(16.*PI*SLOP1)/CMBARN
2442       SIGEL2 = SIG2**2*(1.+RHO2**2)/(16.*PI*SLOP2)/CMBARN
2443       RETURN
2444       END
2445 
2446       SUBROUTINE FPLUS (S, FR, FI)
2447       SAVE
2448 
2449       COMMON /BLOCKC/ AA, BETA, S0, CC, AMU, DD, ALPHA, A0
2450       COMPLEX Z1, Z2, Z3
2451       DATA PI /3.1415926/
2452       F1 = LOG(S/S0)
2453       Z1 = CMPLX(F1,-PI/2.)
2454       Z1 = Z1*Z1
2455       Z2 = 1. + A0*Z1
2456       Z3 = Z1/Z2
2457       F2 = CC*S**(AMU-1.)
2458       F3 = 0.5*PI*(1.-AMU)
2459       FI = AA + F2*COS(F3) + BETA*REAL(Z3)
2460       FR = -BETA*AIMAG(Z3)+F2*SIN(F3)
2461       RETURN
2462       END
2463 
2464       SUBROUTINE FMINUS (S, FR, FI)
2465       SAVE
2466 
2467       COMMON /BLOCKC/ AA, BETA, S0, CC, AMU, DD, ALPHA, A0
2468       DATA PI /3.1415926/
2469       F1 = S**(ALPHA-1.)
2470       F2 = 0.5*PI*(1.-ALPHA)
2471       FR = -DD*F1*COS(F2)
2472       FI = -DD*F1*SIN(F2)
2473       RETURN
2474       END
2475 
2476       SUBROUTINE SSLOPE (S, BP, BM)
2477       SAVE
2478 
2479       COMMON /BLOCKD/ CP, DP, EP, CM, DM
2480       AL = LOG(S)
2481       BP = CP + DP*AL + EP*AL*AL
2482       BM = CM + DM*AL
2483       RETURN
2484       END
2485 
2486       SUBROUTINE BLOCK_INI
2487 C...Parameters of fit IFIT=1 of Block and Cahn
2488       SAVE
2489 
2490       COMMON /BLOCKC/ AA, BETA, S0, CC, AMU, DD, ALPHA, A0
2491       COMMON /BLOCKD/ CP, DP, EP, CM, DM
2492       AA = 41.74
2493       BETA = 0.66
2494       S0 = 338.5
2495       CC = 0.
2496       AMU = 0.
2497       DD = -39.37
2498       ALPHA = 0.48
2499       A0 = 0.
2500       CP = 10.90
2501       DP = -0.08
2502       EP = 0.043
2503       CM = 23.27
2504       DM = 0.93
2505       RETURN
2506       END
2507 
2508 C=============================================================
2509 C.  Nucleus-nucleus cross sections
2510 C=============================================================
2511 
2512       SUBROUTINE SIGNUC_INI (IA,E0)
2513 C...This subroutine receives in INPUT E0 (TeV)
2514 C.  energy per nucleon and computes the cross sections
2515 C.  and interactions lengths for  all nuclei
2516 C.  with A  between 2 and IA
2517 C.  The output is contained in common block /CLENNN/
2518 C.
2519 C.  Attention: the tabulated cross sections are obtained with
2520 C.  new p-p cross sections as used in SIBYLL 2x,
2521 C.  in addition field dimensions changed (RE 04/2000)
2522 C.
2523 C........................................................
2524       COMMON /CLENNN/ SSIGNUC(60), ALNUC(60)
2525       DIMENSION SIGMA(6,56), SIGQE(6,56)
2526       DIMENSION AA(6)
2527       DATA NE /6/, AMIN /1./, DA /1./
2528       DATA AA /1.,2.,3.,4.,5.,6./
2529       DATA AVOG /6.0221367E-04/
2530       DATA ATARGET /14.514/               ! effective masss of air
2531 C...Data on `inelastic-production' nucleus-air cross section
2532       DATA (SIGMA(J, 2),J=1,6) /  392., 434., 507., 617., 739., 866./
2533       DATA (SIGMA(J, 3),J=1,6) /  459., 501., 579., 691., 819., 943./
2534       DATA (SIGMA(J, 4),J=1,6) /  493., 532., 608., 721., 853., 975./
2535       DATA (SIGMA(J, 5),J=1,6) /  589., 632., 716., 839., 986.,1122./
2536       DATA (SIGMA(J, 6),J=1,6) /  695., 746., 838., 974.,1125.,1275./
2537       DATA (SIGMA(J, 7),J=1,6) /  710., 759., 857., 986.,1141.,1292./
2538       DATA (SIGMA(J, 8),J=1,6) /  751., 808., 902.,1039.,1192.,1344./
2539       DATA (SIGMA(J, 9),J=1,6) /  794., 847., 944.,1082.,1241.,1395./
2540       DATA (SIGMA(J,10),J=1,6) /  808., 857., 958.,1095.,1248.,1411./
2541       DATA (SIGMA(J,11),J=1,6) /  808., 858., 958.,1094.,1255.,1409./
2542       DATA (SIGMA(J,12),J=1,6) /  856., 909.,1001.,1144.,1306.,1462./
2543       DATA (SIGMA(J,13),J=1,6) /  868., 916.,1016.,1166.,1324.,1486./
2544       DATA (SIGMA(J,14),J=1,6) /  913., 961.,1063.,1214.,1386.,1543./
2545       DATA (SIGMA(J,15),J=1,6) /  937., 995.,1101.,1240.,1411.,1573./
2546       DATA (SIGMA(J,16),J=1,6) /  998.,1037.,1157.,1304.,1474.,1644./
2547       DATA (SIGMA(J,17),J=1,6) /  991.,1044.,1154.,1310.,1472.,1647./
2548       DATA (SIGMA(J,18),J=1,6) / 1041.,1101.,1208.,1359.,1544.,1718./
2549       DATA (SIGMA(J,19),J=1,6) / 1094.,1151.,1255.,1431.,1614.,1785./
2550       DATA (SIGMA(J,20),J=1,6) / 1141.,1193.,1314.,1485.,1667.,1842./
2551       DATA (SIGMA(J,21),J=1,6) / 1165.,1222.,1341.,1498.,1694.,1880./
2552       DATA (SIGMA(J,22),J=1,6) / 1159.,1219.,1333.,1506.,1687.,1863./
2553       DATA (SIGMA(J,23),J=1,6) / 1211.,1265.,1381.,1553.,1738.,1918./
2554       DATA (SIGMA(J,24),J=1,6) / 1221.,1285.,1402.,1571.,1757.,1947./
2555       DATA (SIGMA(J,25),J=1,6) / 1245.,1314.,1429.,1595.,1792.,1993./
2556       DATA (SIGMA(J,26),J=1,6) / 1270.,1334.,1447.,1624.,1811.,2005./
2557       DATA (SIGMA(J,27),J=1,6) / 1251.,1314.,1431.,1603.,1790.,1984./
2558       DATA (SIGMA(J,28),J=1,6) / 1289.,1351.,1465.,1654.,1839.,2029./
2559       DATA (SIGMA(J,29),J=1,6) / 1412.,1482.,1617.,1796.,2003.,2204./
2560       DATA (SIGMA(J,30),J=1,6) / 1351.,1409.,1526.,1711.,1895.,2105./
2561       DATA (SIGMA(J,31),J=1,6) / 1361.,1422.,1548.,1734.,1932.,2132./
2562       DATA (SIGMA(J,32),J=1,6) / 1381.,1435.,1574.,1751.,1952.,2140./
2563       DATA (SIGMA(J,33),J=1,6) / 1400.,1466.,1593.,1767.,1971.,2173./
2564       DATA (SIGMA(J,34),J=1,6) / 1414.,1479.,1606.,1796.,1985.,2193./
2565       DATA (SIGMA(J,35),J=1,6) / 1425.,1503.,1621.,1800.,2008.,2203./
2566       DATA (SIGMA(J,36),J=1,6) / 1444.,1516.,1643.,1822.,2031.,2225./
2567       DATA (SIGMA(J,37),J=1,6) / 1460.,1529.,1655.,1844.,2042.,2246./
2568       DATA (SIGMA(J,38),J=1,6) / 1484.,1547.,1676.,1861.,2071.,2269./
2569       DATA (SIGMA(J,39),J=1,6) / 1494.,1557.,1700.,1886.,2085.,2291./
2570       DATA (SIGMA(J,40),J=1,6) / 1510.,1585.,1710.,1897.,2116.,2292./
2571       DATA (SIGMA(J,41),J=1,6) / 1532.,1594.,1724.,1922.,2116.,2317./
2572       DATA (SIGMA(J,42),J=1,6) / 1540.,1617.,1735.,1930.,2144.,2334./
2573       DATA (SIGMA(J,43),J=1,6) / 1559.,1627.,1764.,1944.,2158.,2363./
2574       DATA (SIGMA(J,44),J=1,6) / 1570.,1641.,1774.,1949.,2172.,2367./
2575       DATA (SIGMA(J,45),J=1,6) / 1587.,1651.,1784.,1983.,2187.,2396./
2576       DATA (SIGMA(J,46),J=1,6) / 1609.,1675.,1810.,1997.,2206.,2412./
2577       DATA (SIGMA(J,47),J=1,6) / 1611.,1688.,1814.,2011.,2236.,2436./
2578       DATA (SIGMA(J,48),J=1,6) / 1630.,1700.,1833.,2024.,2240.,2457./
2579       DATA (SIGMA(J,49),J=1,6) / 1645.,1717.,1843.,2047.,2251.,2474./
2580       DATA (SIGMA(J,50),J=1,6) / 1655.,1729.,1858.,2051.,2272.,2488./
2581       DATA (SIGMA(J,51),J=1,6) / 1668.,1742.,1875.,2076.,2288.,2504./
2582       DATA (SIGMA(J,52),J=1,6) / 1684.,1757.,1891.,2091.,2308.,2515./
2583       DATA (SIGMA(J,53),J=1,6) / 1697.,1773.,1908.,2103.,2328.,2544./
2584       DATA (SIGMA(J,54),J=1,6) / 1713.,1785.,1917.,2125.,2335.,2533./
2585       DATA (SIGMA(J,55),J=1,6) / 1726.,1803.,1939.,2131.,2355.,2570./
2586       DATA (SIGMA(J,56),J=1,6) / 1757.,1831.,1970.,2168.,2388.,2620./
2587 C...Data on `quasi-elastic' nucleus-air cross section
2588       DATA (SIGQE(J, 2),J=1,6) /   41.,  39.,  54.,  99., 159., 229./
2589       DATA (SIGQE(J, 3),J=1,6) /   43.,  40.,  55., 101., 167., 238./
2590       DATA (SIGQE(J, 4),J=1,6) /   40.,  40.,  55., 102., 168., 243./
2591       DATA (SIGQE(J, 5),J=1,6) /   45.,  44.,  61., 113., 179., 262./
2592       DATA (SIGQE(J, 6),J=1,6) /   52.,  49.,  67., 121., 195., 278./
2593       DATA (SIGQE(J, 7),J=1,6) /   50.,  48.,  69., 121., 194., 276./
2594       DATA (SIGQE(J, 8),J=1,6) /   52.,  50.,  70., 125., 199., 282./
2595       DATA (SIGQE(J, 9),J=1,6) /   54.,  50.,  69., 126., 206., 289./
2596       DATA (SIGQE(J,10),J=1,6) /   53.,  50.,  70., 127., 204., 288./
2597       DATA (SIGQE(J,11),J=1,6) /   51.,  50.,  69., 125., 203., 285./
2598       DATA (SIGQE(J,12),J=1,6) /   52.,  51.,  70., 127., 210., 293./
2599       DATA (SIGQE(J,13),J=1,6) /   53.,  50.,  71., 129., 210., 293./
2600       DATA (SIGQE(J,14),J=1,6) /   52.,  53.,  72., 134., 212., 301./
2601       DATA (SIGQE(J,15),J=1,6) /   54.,  52.,  75., 135., 219., 306./
2602       DATA (SIGQE(J,16),J=1,6) /   57.,  56.,  76., 136., 222., 313./
2603       DATA (SIGQE(J,17),J=1,6) /   55.,  54.,  76., 137., 221., 311./
2604       DATA (SIGQE(J,18),J=1,6) /   59.,  55.,  79., 142., 223., 320./
2605       DATA (SIGQE(J,19),J=1,6) /   63.,  58.,  81., 143., 231., 323./
2606       DATA (SIGQE(J,20),J=1,6) /   63.,  60.,  82., 149., 237., 332./
2607       DATA (SIGQE(J,21),J=1,6) /   64.,  61.,  84., 150., 236., 333./
2608       DATA (SIGQE(J,22),J=1,6) /   62.,  61.,  84., 145., 237., 329./
2609       DATA (SIGQE(J,23),J=1,6) /   65.,  62.,  83., 149., 240., 341./
2610       DATA (SIGQE(J,24),J=1,6) /   64.,  61.,  84., 152., 244., 338./
2611       DATA (SIGQE(J,25),J=1,6) /   64.,  61.,  86., 153., 249., 343./
2612       DATA (SIGQE(J,26),J=1,6) /   66.,  64.,  88., 156., 251., 344./
2613       DATA (SIGQE(J,27),J=1,6) /   65.,  61.,  84., 155., 244., 345./
2614       DATA (SIGQE(J,28),J=1,6) /   65.,  63.,  87., 153., 249., 347./
2615       DATA (SIGQE(J,29),J=1,6) /   72.,  64.,  91., 161., 260., 362./
2616       DATA (SIGQE(J,30),J=1,6) /   67.,  63.,  90., 157., 252., 355./
2617       DATA (SIGQE(J,31),J=1,6) /   67.,  66.,  90., 161., 256., 353./
2618       DATA (SIGQE(J,32),J=1,6) /   68.,  65.,  89., 161., 256., 355./
2619       DATA (SIGQE(J,33),J=1,6) /   69.,  66.,  89., 161., 257., 360./
2620       DATA (SIGQE(J,34),J=1,6) /   69.,  63.,  89., 160., 262., 359./
2621       DATA (SIGQE(J,35),J=1,6) /   68.,  65.,  90., 166., 259., 360./
2622       DATA (SIGQE(J,36),J=1,6) /   68.,  66.,  92., 167., 257., 362./
2623       DATA (SIGQE(J,37),J=1,6) /   71.,  67.,  92., 166., 261., 365./
2624       DATA (SIGQE(J,38),J=1,6) /   70.,  67.,  92., 167., 265., 367./
2625       DATA (SIGQE(J,39),J=1,6) /   69.,  68.,  96., 163., 261., 370./
2626       DATA (SIGQE(J,40),J=1,6) /   71.,  67.,  92., 167., 266., 371./
2627       DATA (SIGQE(J,41),J=1,6) /   72.,  68.,  93., 167., 267., 373./
2628       DATA (SIGQE(J,42),J=1,6) /   71.,  69.,  96., 167., 271., 370./
2629       DATA (SIGQE(J,43),J=1,6) /   72.,  70.,  95., 168., 265., 369./
2630       DATA (SIGQE(J,44),J=1,6) /   70.,  68.,  98., 170., 270., 372./
2631       DATA (SIGQE(J,45),J=1,6) /   72.,  69.,  97., 167., 272., 379./
2632       DATA (SIGQE(J,46),J=1,6) /   75.,  69.,  98., 173., 272., 370./
2633       DATA (SIGQE(J,47),J=1,6) /   73.,  70.,  93., 170., 275., 377./
2634       DATA (SIGQE(J,48),J=1,6) /   73.,  71.,  96., 170., 275., 379./
2635       DATA (SIGQE(J,49),J=1,6) /   75.,  72.,  96., 171., 275., 379./
2636       DATA (SIGQE(J,50),J=1,6) /   74.,  73.,  97., 174., 275., 380./
2637       DATA (SIGQE(J,51),J=1,6) /   73.,  68.,  97., 173., 281., 388./
2638       DATA (SIGQE(J,52),J=1,6) /   75.,  71.,  99., 175., 275., 390./
2639       DATA (SIGQE(J,53),J=1,6) /   75.,  71.,  98., 180., 279., 383./
2640       DATA (SIGQE(J,54),J=1,6) /   74.,  71.,  99., 177., 280., 386./
2641       DATA (SIGQE(J,55),J=1,6) /   75.,  71., 100., 179., 278., 394./
2642       DATA (SIGQE(J,56),J=1,6) /   76.,  70., 101., 175., 284., 390./
2643 
2644       ASQS = 0.5*LOG10(1.876E+03*E0)
2645       JE = MIN(INT((ASQS-AMIN)/DA)+1,NE-2)
2646       DO JA=2,IA
2647          ABEAM = FLOAT(JA)
2648          S1 = QUAD_INT(ASQS, AA(JE),AA(JE+1),AA(JE+2),
2649      +                   SIGMA(JE,JA),SIGMA(JE+1,JA),SIGMA(JE+2,JA))
2650          S2 = QUAD_INT(ASQS, AA(JE),AA(JE+1),AA(JE+2),
2651      +                   SIGQE(JE,JA),SIGQE(JE+1,JA),SIGQE(JE+2,JA))
2652          SSIGNUC(JA) = S1 + S2
2653          ALNUC(JA) = ATARGET/(AVOG*SSIGNUC(JA))
2654       ENDDO
2655       ALNUC(1) = FPNI(E0, 13)
2656       SSIGNUC(1) = ATARGET/(AVOG*ALNUC(1))
2657 
2658       RETURN
2659       END
2660 
2661 
2662 C=======================================================================
2663 C.  General utilities
2664 C.======================================================================
2665 
2666       FUNCTION QUAD_INT (R,X0,X1,X2,V0,V1,V2)
2667 C...Quadratic interpolation
2668       SAVE
2669 
2670       R0=R-X0
2671       R1=R-X1
2672       R2=R-X2
2673       S0=X0-X1
2674       S1=X0-X2
2675       S2=X1-X2
2676       QUAD_INT = V0*R1*R2/(S0*S1)-V1*R0*R2/(S0*S2)+V2*R0*R1/(S1*S2)
2677       RETURN
2678       END
2679 
2680       FUNCTION GAUSS (FUN, A,B)
2681 C...Returns the  8 points Gauss-Legendre integral
2682 C.  of function FUN from A to B
2683 C...........................................................
2684       SAVE
2685 
2686       DIMENSION X(8), W(8)
2687       DATA X / .0950125098, .2816035507, .4580167776, .6178762444
2688      1          ,.7554044083, .8656312023, .9445750230, .9894009349/
2689       DATA W / .1894506104, .1826034150, .1691565193, .1495959888
2690      1          ,.1246289712, .0951585116, .0622535239, .0271524594/
2691       XM = 0.5*(B+A)
2692       XR = 0.5*(B-A)
2693       SS = 0.
2694       DO J=1,8
2695         DX = XR*X(J)
2696         SS = SS + W(J) * (FUN(XM+DX) + FUN(XM-DX))
2697       ENDDO
2698       GAUSS = XR*SS
2699       RETURN
2700       END
2701 
2702 ctp060203      FUNCTION GASDEV(Idum)
2703       FUNCTION GASDEV(Idum)
2704 C...Gaussian deviation
2705       SAVE
2706 
2707 ctp060123      SAVE GSET
2708       DATA ISET/0/
2709       IF (ISET.EQ.0) THEN
2710 1       V1=2.*S_RNDM(ISET)-1.
2711         V2=2.*S_RNDM(Idum)-1.
2712         R=V1**2+V2**2
2713         IF(R.GE.1.)GO TO 1
2714         FAC=SQRT(-2.*LOG(R)/R)
2715         GSET=V1*FAC
2716         GASDEV=V2*FAC
2717         ISET=1
2718       ELSE
2719         GASDEV=GSET
2720         ISET=0
2721       ENDIF
2722       RETURN
2723       END
2724 
2725       subroutine invert_array (yy, xmin, dx, n, xnew, ymin, dy)
2726 C..    This subroutine receives one   array
2727 C      of n y values in input yy(1:n)
2728 C      that correspond to  equispaced values of x_j = xmin + dx*(j-1)
2729 C
2730 C      and "reverse" the array returning an array of  x values
2731 C      xnew (1:n) that  corresponds to equispaced values of y
2732 C      The relation is assumed monotonous but can be
2733 C      increasing or decreasing
2734 C..............................................................
2735       SAVE
2736 
2737       dimension  yy(n), xnew (n)
2738       ymin = yy(1)
2739       ymax = yy(n)
2740       dy = (ymax - ymin)/float(n-1)
2741       xnew (1) = xmin
2742       xnew (n) = xmin + dx*float(n-1)
2743       k0 = 1
2744       do j=2,n-1
2745          y = ymin + float(j-1)*dy
2746          do k=k0,n
2747             if((yy(k) .gt. y) .eqv. (yy(n) .gt. yy(1))) goto 100
2748          enddo
2749 100      y2 = yy(k)
2750          y1 = yy(k-1)
2751          k0 = k-1
2752          x1 = xmin + dx*float(k-2)
2753          x2 = x1+dx
2754          xnew (j)  = x1 + dx* (y-y1)/(y2-y1)
2755       enddo
2756       return
2757       end
2758 C->
2759       SUBROUTINE SINCO(S,C)
2760       SAVE
2761 
2762       DATA PI /3.1415926/
2763       F = 2.*PI*S_RNDM(0)
2764       C = COS (F)
2765       S = SIN (F)
2766       RETURN
2767       END
2768 
2769 
2770 
2771 C=============================================================
2772 C.  Cross sections for cascade calculations (FPNI)
2773 C=============================================================
2774 
2775 
2776 
2777       SUBROUTINE SIGMA_PP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
2778 C-----------------------------------------------------------------------
2779 C...p-p cross sections
2780 C.
2781 C.  this routine serves the purpose to calculate cascades with different
2782 C.  cross sections
2783 C.
2784 C. INPUT: E0 = Laboratory Energy  (TeV)
2785 C.
2786 C. OUTPUT: SIGT = total cross section
2787 C.         SIGEL = elastic cross section
2788 C.         SIGINEL = inelastic cross section
2789 C.         SLOPE = slope of elastic scattering (GeV**-2)
2790 C.         RHO = Imaginary/Real part of forward elastic amplitude
2791 C.
2792 C.  (old cross section tables end at 10^6 GeV)
2793 C-----------------------------------------------------------------------
2794       SAVE
2795 
2796       DIMENSION SSIG0(51)
2797       DIMENSION SIGDIF(3)
2798       DATA ICSPA /0/
2799       DATA PI /3.1415926/
2800       DATA CMBARN /0.389385/
2801 
2802 C...p-p inelastic cross sections (mbarn)
2803       DATA (SSIG0(J),J=1,51) /
2804      +      32.05,    32.06,    32.08,    32.13,    32.22,    32.36,
2805      +      32.56,    32.85,    33.24,    33.75,    34.37,    35.14,
2806      +      36.05,    37.12,    38.37,    39.78,    41.36,    43.13,
2807      +      45.07,    47.18,    49.47,    51.91,    54.54,    57.28,
2808      +      60.15,    63.15,    66.28,    69.48,    72.80,    76.22,
2809      +      79.71,    83.27,    86.87,    90.55,    94.26,    98.05,
2810      +     101.89,   105.75,   109.71,   113.65,   117.60,   121.55,
2811      +     125.53,   129.56,   133.60,   137.70,   141.77,   145.84,
2812      +     149.92,   154.02,   158.15/
2813 
2814 
2815       SQS = SQRT(2000.*0.938*E0)
2816 
2817 *  old standard NUCLIB/SIBYLL model
2818 
2819       IF(ICSPA.EQ.-1) THEN
2820 
2821         AL = LOG10(SQS)
2822         if(AL.le.1.) then
2823           SIGINEL = SSIG0(1)
2824         else
2825           J1 = (AL - 1.)*10. + 1
2826           J1 = min(J1,50)
2827           T = (AL-1.)*10. - FLOAT(J1-1)
2828           SIGINEL = SSIG0(J1)*(1.-T) + SSIG0(J1+1)*T
2829         endif
2830         CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,SIGEL1,SIGEL2)
2831         R = SIGEL1/SIGT1
2832         RHO = RHO1
2833         SIGT  = SIGINEL/(1.-R)
2834         SIGEL = SIGINEL*R/(1.-R)
2835         SLOPE = SIGT**2/(SIGEL * 16.*PI) * (1.+RHO1**2) /CMBARN
2836 
2837 *  cross section as calculated in SIBYLL
2838 
2839       ELSE IF(ICSPA.EQ.0) THEN
2840 
2841         CALL SIB_SIGMA_HP(1,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
2842 
2843 *  Donnachie-Landshoff  (sig-tot)
2844 
2845       ELSE IF(ICSPA.EQ.1) THEN
2846 
2847         CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,
2848      +             SIGEL1,SIGEL2)
2849         R = SIGEL1/SIGT1
2850         RHO = RHO1
2851 
2852         DELDL = 0.0808
2853         EPSDL = -0.4525
2854         S = SQS*SQS
2855         SIGT = 21.7*S**DELDL+56.08*S**EPSDL
2856         SIGEL = R*SIGT
2857         SIGINEL = SIGT-SIGEL
2858         SLOPE = SIGT**2/(SIGEL * 16.*PI) * (1.+RHO**2) /CMBARN
2859 
2860 *  Donnachie-Landshoff (sig-tot and sig-el)
2861 
2862       ELSE IF(ICSPA.EQ.2) THEN
2863 
2864         DELDL = 0.0808
2865         EPSDL = -0.4525
2866         S = SQS*SQS
2867         SIGT = 21.7*S**DELDL+56.08*S**EPSDL
2868         IMODEL = 1
2869         IF(IMODEL.EQ.1) THEN
2870           ALPHAP = 0.25D0
2871           SLOPE = 8.5D0+2.D0*ALPHAP*LOG(S)
2872         ELSE IF(IMODEL.EQ.2) THEN
2873           ALPHAP = 0.3D0
2874           SLOPE = 8.D0+2.D0*ALPHAP*LOG(S)
2875         ENDIF
2876         SIGEL = SIGT**2/(16.D0*PI*SLOPE*CMBARN)
2877         SIGINEL = SIGT-SIGEL
2878         RHO = 0.
2879 
2880 *  geometrical scaling with Donnachie-Landshoff sig-tot
2881 
2882       ELSE IF(ICSPA.EQ.3) THEN
2883 
2884         R = 0.17D0
2885 
2886         DELDL = 0.0808
2887         EPSDL = -0.4525
2888         S = SQS*SQS
2889         SIGT = 21.7*S**DELDL+56.08*S**EPSDL
2890 
2891         SIGEL = R*SIGT
2892         SIGINEL = SIGT-SIGEL
2893         SLOPE = SIGT**2/(16*PI*SIGEL)/CMBARN
2894         RHO = 0.
2895 
2896       ENDIF
2897 
2898       RETURN
2899       END
2900 
2901 
2902       SUBROUTINE SIGMA_PIP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
2903 C-----------------------------------------------------------------------
2904 C...pi-p cross sections
2905 C.
2906 C.  this routine serves the purpose to calculate cascades with different
2907 C.  cross sections
2908 C.
2909 C. INPUT: E0 = Laboratory Energy  (TeV)
2910 C.
2911 C. OUTPUT: SIGT = total cross section
2912 C.         SIGEL = elastic cross section
2913 C.         SIGINEL = inelastic cross section
2914 C.         SLOPE = slope of elastic scattering (GeV**-2)
2915 C.         RHO = Imaginary/Real part of forward elastic amplitude
2916 C.
2917 C.  (old cross section tables end at 10^6 GeV)
2918 C-----------------------------------------------------------------------
2919       SAVE
2920 
2921       DIMENSION SSIG0(51)
2922       DIMENSION SIGDIF(3)
2923       DATA ICSPA /0/
2924       DATA PI /3.1415926/
2925       DATA CMBARN /0.389385/
2926 
2927 C...pi-p inelastic cross sections (mbarn)
2928       DATA (SSIG0(J),J=1,51) /
2929      +      20.76,    20.78,    20.81,    20.88,    20.98,    21.13,
2930      +      21.33,    21.61,    21.96,    22.39,    22.92,    23.56,
2931      +      24.31,    25.18,    26.18,    27.32,    28.60,    30.04,
2932      +      31.64,    33.40,    35.34,    37.43,    39.72,    42.16,
2933      +      44.77,    47.56,    50.53,    53.66,    56.99,    60.50,
2934      +      64.17,    68.03,    72.05,    76.27,    80.67,    85.27,
2935      +      90.08,    95.04,   100.27,   105.65,   111.21,   116.94,
2936      +     122.87,   129.03,   135.37,   141.93,   148.62,   155.49,
2937      +     162.48,   169.60,   176.94/
2938 
2939       SQS = SQRT(2000.*0.938*E0)
2940 
2941 *  old standard NUCLIB/SIBYLL model
2942 
2943       IF(ICSPA.EQ.-1) THEN
2944 
2945         AL = LOG10(SQS)
2946         if(AL.le.1.) then
2947           SIGINEL = SSIG0(1)
2948         else
2949           J1 = (AL - 1.)*10. + 1
2950           J1 = min(J1,50)
2951           T = (AL-1.)*10. - FLOAT(J1-1)
2952           SIGINEL = SSIG0(J1)*(1.-T) + SSIG0(J1+1)*T
2953         endif
2954         CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,SIGEL1,SIGEL2)
2955         R = SIGEL1/SIGT1
2956         RHO = RHO1
2957         SIGT  = SIGINEL/(1.-R)
2958         SIGEL = SIGINEL*R/(1.-R)
2959         SLOPE = SIGT**2/(SIGEL * 16.*PI) * (1.+RHO1**2) /CMBARN
2960 
2961 *  cross section as calculated in SIBYLL
2962 
2963       ELSE IF(ICSPA.EQ.0) THEN
2964 
2965         CALL SIB_SIGMA_HP(2,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
2966 
2967 *  Donnachie-Landshoff  (sig-tot)
2968 
2969       ELSE IF(ICSPA.EQ.1) THEN
2970 
2971         CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,
2972      +             SIGEL1,SIGEL2)
2973         R = SIGEL1/SIGT1
2974         RHO = RHO1
2975 
2976         DELDL = 0.0808
2977         EPSDL = -0.4525
2978         S = SQS*SQS
2979         SIGT = 13.63*S**DELDL+(36.02+27.56)/2.*S**EPSDL
2980         SIGEL = R*SIGT
2981         SIGINEL = SIGT-SIGEL
2982         SLOPE = SIGT**2/(SIGEL * 16.*PI) * (1.+RHO**2) /CMBARN
2983 
2984 *  Donnachie-Landshoff (sig-tot and sig-el)
2985 
2986       ELSE IF(ICSPA.EQ.2) THEN
2987 
2988         DELDL = 0.0808
2989         EPSDL = -0.4525
2990         S = SQS*SQS
2991         SIGT = 13.63*S**DELDL+(36.02+27.56)/2.*S**EPSDL
2992         IMODEL = 1
2993         IF(IMODEL.EQ.1) THEN
2994           ALPHAP = 0.25D0
2995           SLOPE = 8.5D0+2.D0*ALPHAP*LOG(S)
2996         ELSE IF(IMODEL.EQ.2) THEN
2997           ALPHAP = 0.3D0
2998           SLOPE = 8.D0+2.D0*ALPHAP*LOG(S)
2999         ENDIF
3000         SIGEL = SIGT**2/(16.D0*PI*SLOPE*CMBARN)
3001         SIGINEL = SIGT-SIGEL
3002         RHO = 0.
3003 
3004 *  geometrical scaling with Donnachie-Landshoff sig-tot
3005 
3006       ELSE IF(ICSPA.EQ.3) THEN
3007 
3008         R = 0.17D0
3009 
3010         DELDL = 0.0808
3011         EPSDL = -0.4525
3012         S = SQS*SQS
3013         SIGT = 13.63*S**DELDL+(36.02+27.56)/2.*S**EPSDL
3014 
3015         SIGEL = R*SIGT
3016         SIGINEL = SIGT-SIGEL
3017         SLOPE = SIGT**2/(16*PI*SIGEL)/CMBARN
3018         RHO = 0.
3019 
3020       ENDIF
3021 
3022       RETURN
3023       END
3024 
3025 
3026       SUBROUTINE SIGMA_KP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
3027 C-----------------------------------------------------------------------
3028 C...K-p cross sections
3029 C.
3030 C.  this routine serves the purpose to calculate cascades with different
3031 C.  cross sections
3032 C.
3033 C.  if old cross sections are selected then sigma_pi = sigma_K
3034 C.
3035 C. INPUT: E0 = Laboratory Energy  (TeV)
3036 C.
3037 C. OUTPUT: SIGT = total cross section
3038 C.         SIGEL = elastic cross section
3039 C.         SIGINEL = inelastic cross section
3040 C.         SLOPE = slope of elastic scattering (GeV**-2)
3041 C.         RHO = Imaginary/Real part of forward elastic amplitude
3042 C.
3043 C.  (old cross section tables end at 10^6 GeV)
3044 C-----------------------------------------------------------------------
3045       SAVE
3046 
3047       DIMENSION SSIG0(51)
3048       DIMENSION SIGDIF(3)
3049       DATA ICSPA /0/
3050       DATA PI /3.1415926/
3051       DATA CMBARN /0.389385/
3052 
3053 C...pi-p inelastic cross sections (mbarn)
3054       DATA (SSIG0(J),J=1,51) /
3055      +      20.76,    20.78,    20.81,    20.88,    20.98,    21.13,
3056      +      21.33,    21.61,    21.96,    22.39,    22.92,    23.56,
3057      +      24.31,    25.18,    26.18,    27.32,    28.60,    30.04,
3058      +      31.64,    33.40,    35.34,    37.43,    39.72,    42.16,
3059      +      44.77,    47.56,    50.53,    53.66,    56.99,    60.50,
3060      +      64.17,    68.03,    72.05,    76.27,    80.67,    85.27,
3061      +      90.08,    95.04,   100.27,   105.65,   111.21,   116.94,
3062      +     122.87,   129.03,   135.37,   141.93,   148.62,   155.49,
3063      +     162.48,   169.60,   176.94/
3064 
3065       SQS = SQRT(2000.*0.938*E0)
3066 
3067 *  old standard NUCLIB/SIBYLL model
3068 
3069       IF(ICSPA.EQ.-1) THEN
3070 
3071         AL = LOG10(SQS)
3072         if(AL.le.1.) then
3073           SIGINEL = SSIG0(1)
3074         else
3075           J1 = (AL - 1.)*10. + 1
3076           J1 = min(J1,50)
3077           T = (AL-1.)*10. - FLOAT(J1-1)
3078           SIGINEL = SSIG0(J1)*(1.-T) + SSIG0(J1+1)*T
3079         endif
3080         CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,SIGEL1,SIGEL2)
3081         R = SIGEL1/SIGT1
3082         RHO = RHO1
3083         SIGT  = SIGINEL/(1.-R)
3084         SIGEL = SIGINEL*R/(1.-R)
3085         SLOPE = SIGT**2/(SIGEL * 16.*PI) * (1.+RHO1**2) /CMBARN
3086 
3087 *  cross section as calculated in SIBYLL
3088 
3089       ELSE IF(ICSPA.EQ.0) THEN
3090 
3091         CALL SIB_SIGMA_HP(3,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
3092 
3093 *  Donnachie-Landshoff  (sig-tot)
3094 
3095       ELSE IF(ICSPA.EQ.1) THEN
3096 
3097         CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,
3098      +             SIGEL1,SIGEL2)
3099         R = SIGEL1/SIGT1
3100         RHO = RHO1
3101 
3102         DELDL = 0.0808
3103         EPSDL = -0.4525
3104         S = SQS*SQS
3105         SIGT = 11.82*S**DELDL+(26.36+ 8.15)/2.*S**EPSDL
3106         SIGEL = R*SIGT
3107         SIGINEL = SIGT-SIGEL
3108         SLOPE = SIGT**2/(SIGEL * 16.*PI) * (1.+RHO**2) /CMBARN
3109 
3110 *  Donnachie-Landshoff (sig-tot and sig-el)
3111 
3112       ELSE IF(ICSPA.EQ.2) THEN
3113 
3114         DELDL = 0.0808
3115         EPSDL = -0.4525
3116         S = SQS*SQS
3117         SIGT = 11.82*S**DELDL+(26.36+ 8.15)/2.*S**EPSDL
3118         IMODEL = 1
3119         IF(IMODEL.EQ.1) THEN
3120           ALPHAP = 0.25D0
3121           SLOPE = 8.5D0+2.D0*ALPHAP*LOG(S)
3122         ELSE IF(IMODEL.EQ.2) THEN
3123           ALPHAP = 0.3D0
3124           SLOPE = 8.D0+2.D0*ALPHAP*LOG(S)
3125         ENDIF
3126         SIGEL = SIGT**2/(16.D0*PI*SLOPE*CMBARN)
3127         SIGINEL = SIGT-SIGEL
3128         RHO = 0.
3129 
3130 *  geometrical scaling with Donnachie-Landshoff sig-tot
3131 
3132       ELSE IF(ICSPA.EQ.3) THEN
3133 
3134         R = 0.17D0
3135 
3136         DELDL = 0.0808
3137         EPSDL = -0.4525
3138         S = SQS*SQS
3139         SIGT = 11.82*S**DELDL+(26.36+ 8.15)/2.*S**EPSDL
3140 
3141         SIGEL = R*SIGT
3142         SIGINEL = SIGT-SIGEL
3143         SLOPE = SIGT**2/(16*PI*SIGEL)/CMBARN
3144         RHO = 0.
3145 
3146       ENDIF
3147 
3148       RETURN
3149       END
3150 
3151 
3152 
3153       SUBROUTINE SIGMA_INI
3154 C-----------------------------------------------------------------------
3155 C.  Initialize the cross section and interaction lengths on air
3156 C-----------------------------------------------------------------------
3157       SAVE
3158 
3159       COMMON /CSAIR/ NSQS, ASQSMIN, ASQSMAX, DASQS,
3160      &               SSIG0(61,3),SSIGA(61,3),ALINT(61,3)
3161 
3162       DATA AVOG /6.0221367E-04/
3163       ATARGET = 14.514
3164 
3165       CALL BLOCK_INI
3166 
3167 C...Loop on c.m. energy
3168       NSQS = 61
3169       SQSMIN = 10.
3170       SQSMAX = 1.E+07
3171       ASQSMIN = LOG10(SQSMIN)
3172       ASQSMAX = LOG10(SQSMAX)
3173       DASQS = (ASQSMAX-ASQSMIN)/FLOAT(NSQS-1)
3174       DO J=1,NSQS
3175          ASQS = ASQSMIN + DASQS*FLOAT(J-1)
3176          SQS = 10.**ASQS
3177          E0 = SQS*SQS/(2.*0.938) * 1.E-03       ! TeV
3178 C...p-air
3179          CALL SIGMA_PP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
3180          CALL SIG_H_AIR (SIGT, SLOPE, RHO, SSIGT, SSIGEL, SSIGQE)
3181          SSIGA(J,1) = SSIGT-SSIGQE
3182          SSIG0(J,1) = SIGINEL
3183          ALINT(J,1) = 1./(AVOG*SSIGA(J,1)/ATARGET)
3184 C...pi-air
3185          CALL SIGMA_PIP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
3186          CALL  SIG_H_AIR (SIGT, SLOPE, RHO, SSIGT, SSIGEL, SSIGQE)
3187          SSIGA(J,2) = SSIGT-SSIGQE
3188          SSIG0(J,2) = SIGINEL
3189          ALINT(J,2) = 1./(AVOG*SSIGA(J,2)/ATARGET)
3190 C...K-air
3191          CALL SIGMA_KP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
3192          CALL  SIG_H_AIR (SIGT, SLOPE, RHO, SSIGT, SSIGEL, SSIGQE)
3193          SSIGA(J,3) = SSIGT-SSIGQE
3194          SSIG0(J,3) = SIGINEL
3195          ALINT(J,3) = 1./(AVOG*SSIGA(J,3)/ATARGET)
3196       ENDDO
3197 
3198 *     WRITE(6,'(1X,A)')
3199 *    &  'SIGMA_INI: NUCLIB interaction lengths (p-air, pi-air, K-air)'
3200 *     DO J=1,NSQS
3201 *        SQS = 10.**(ASQSMIN + DASQS*FLOAT(J-1))
3202 *        WRITE(6,'(1X,1P,4E12.3)') SQS,ALINT(J,1),ALINT(J,2),ALINT(J,3)
3203 *     ENDDO
3204 
3205       RETURN
3206       END
3207 
3208 
3209       FUNCTION FPNI (E,Linp)
3210 C-----------------------------------------------------------------------
3211 C...This function  returns the interaction length
3212 C.  of an hadronic particle travelling in air
3213 C.
3214 C.  INPUT:   E (TeV)   particle energy
3215 C.           Linp      particle code
3216 C.  OUTPUT:  FPNI      (g cm-2)
3217 C-----------------------------------------------------------------------
3218       SAVE
3219 
3220       COMMON /CSAIR/ NSQS, ASQSMIN, ASQSMAX, DASQS,
3221      &               SSIG0(61,3),SSIGA(61,3),ALINT(61,3)
3222 
3223       DIMENSION KK(6:14)
3224       DATA KK /3*2, 4*3, 2*1/
3225 
3226       SQS = SQRT(2000.*E*0.937)                        ! GeV
3227       AL = LOG10 (SQS)
3228       L = abs(Linp)
3229       IF (AL .LE. ASQSMIN)  THEN
3230          FPNI = ALINT(1,KK(L))
3231       ELSE
3232          T = (AL-ASQSMIN)/DASQS
3233          J = INT(T)
3234          J = MIN(J,NSQS-2)
3235          T = T-FLOAT(J)
3236          FPNI = ((1.-T)*ALINT(J+1,KK(L)) + T*ALINT(J+2,KK(L)))
3237       ENDIF
3238 
3239       RETURN
3240       END
3241 
3242 
3243 
3244       SUBROUTINE INT_LEN_INI
3245 C-----------------------------------------------------------------------
3246 C...Initialize the interaction lengths from NUCLIB
3247 C-----------------------------------------------------------------------
3248       SAVE
3249 
3250       CALL NUC_GEOM_INI                 ! nucleus profiles
3251       CALL SIGMA_INI                    ! initialize cross sections
3252 
3253       RETURN
3254       END
3255