diff --git a/Processes/Sibyll/sibyll2.3c.f b/Processes/Sibyll/sibyll2.3c.f
index e757320c8c8cf72a652bd14b4eb880d92b93f187..9788a8bfaa30738fbaa9f5ff89f58f22f4a689dc 100644
--- a/Processes/Sibyll/sibyll2.3c.f
+++ b/Processes/Sibyll/sibyll2.3c.f
@@ -17840,3560 +17840,6 @@ c     Knuth shuffle..
       KF2=KPL(2)*10+KPL(3)     
       END
             
-C.=========================================================================
-C.  Library of programs for the generation of nucleus-nucleus interactions
-C.  and the study of nucleus-induced cosmic ray showers
-C.
-C.  September 2001  changes in FPNI, and SIGMA_INI,
-C.                  new SIGMA_PP, SIGMA_PPI, SIGMA_KP  (R. Engel)
-C.
-C.  may  1996       small bug  corrected by Dieter Heck in NUC_CONF 
-C.
-C.  march 1996      small modification to the superposition code
-C.
-C.  February 1996   change to FPNI to give an interaction length
-C.                   also  at very low energy  
-C.
-C.  Version 1.01  september 1995 
-C.       (small corrections P.L.)
-C.       the random number generator is called as S_RNDM(0)
-C.  ------------------------------------------------------
-C.  Version 1.00  April 1992
-C.
-C.  Authors:
-C.
-C.     J. Engel
-C.     T.K Gaisser
-C.     P.Lipari
-C.     T. Stanev
-C. 
-C.  This set of routines  when used in  the simulation of cosmic ray
-C.  showers have only three  "contact points" with the "external world"
-C.
-C.    (i) SUBROUTINE NUC_NUC_INI
-C.        (no  calling arguments)         
-C.         to be called once during general initialization
-C.    (ii) SUBROUTINE HEAVY (IA, E0)
-C.         where IA (integer) is the mass number of the projectile
-C.         nucleus  and E0 (TeV) is the energy per nucleon
-C.         The output (positions of first interaction for the IA
-C.         nucleons of the projectile) is  contained in  the common block:
-C.           COMMON /C1STNC/ XX0(60),XX(60),YY(60),AX(60),AY(60)
-C.         In detail:
-C.             XX0(j)   (g cm-2) =  position of interaction
-C.             XX(j) (mm)    x-distance from shower axis
-C.             YY(j) (mm)    y-distance from shower axis
-C.             AX(j) (radiants)  Theta_x with respect to original direction
-C.             AY(j) (radiants)  Theta_y with respect to original direction
-C.      
-C.    (iii)  FUNCTION FPNI (E,L)
-C.           Interaction length in air.
-C.           E (TeV) is the energy of the particle, L is the particle
-C.           code (NOTE: "Sibyll" codes are used : L =1-18) 
-C.           WANRNING : The nucleus-nucleus cross section
-C.           tabulated in the program are "matched" to the p-Air
-C.           cross section calculated  with this FUNCTION, in other words 
-C.           they are both calculated with the same input pp cross section
-C==========================================================================
-
-      SUBROUTINE NUC_NUC_INI
-
-C-----------------------------------------------------------------------
-C...Initialization for the generation of nucleus-nucleus interactions
-C.  INPUT : E0 (TeV) Energy per nucleon of the beam nucleus
-C........................................................................
-      SAVE
-
-      CALL NUC_GEOM_INI                       ! nucleus profiles
-      CALL SIGMA_INI                          ! initialize pp cross sections
-
-      RETURN
-      END
-C=======================================================================
-
-      SUBROUTINE FRAGM1 (IA,NW, NF, IAF)
-
-C-----------------------------------------------------------------------
-C...Nuclear Fragmentation 
-C.  total dissolution of nucleus
-C.......................................................................
-      SAVE
-
-      DIMENSION IAF(60)
-      NF = IA-NW
-      DO J=1,NF
-         IAF(J) = 1
-      ENDDO
-      RETURN
-      END
-C->
-C=======================================================================
-
-      SUBROUTINE FRAGM2 (IA,NW, NF, IAF)
-
-C-----------------------------------------------------------------------
-C...Nuclear Fragmentation 
-C.  Spectator in one single fragment 
-C.......................................................................
-      SAVE
-
-      DIMENSION IAF(60)
-      IF (IA-NW .GT. 0)  THEN
-         NF = 1
-         IAF(1) = IA-NW
-      ELSE
-         NF = 0
-      ENDIF
-      RETURN
-      END
-
-C-----------------------------------------------------------------------
-C...Code of fragmentation  of spectator nucleons
-C.  based on Jon Engel  abrasion-ablation algorithms
-C=======================================================================
-
-      BLOCK DATA FRAG_DATA
-
-C-----------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-
-C...Data for the fragmentation of  nucleus  projectiles
-      COMMON /FRAGMOD/A(10,10,20),AE(10,10,20),ERES(10,10),NFLAGG(10,10)
-      SAVE
-      DATA (NFLAGG(I, 1),I=1,10)  / 
-     +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
-      DATA (NFLAGG(I, 2),I=1,10)  / 
-     +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
-      DATA (NFLAGG(I, 3),I=1,10)  / 
-     +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
-      DATA (NFLAGG(I, 4),I=1,10)  / 
-     +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
-      DATA (NFLAGG(I, 5),I=1,10)  / 
-     +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
-      DATA (NFLAGG(I, 6),I=1,10)  / 
-     +    0,  0,  0,  0,  0,  0,  0,  1,  1,  1 /
-      DATA (NFLAGG(I, 7),I=1,10)  / 
-     +    1,  1,  1,  1,  1,  1,  1,  1,  1,  1 /
-      DATA (NFLAGG(I, 8),I=1,10)  / 
-     +    1,  1,  1,  1,  1,  1,  1,  1,  1,  1 /
-      DATA (NFLAGG(I, 9),I=1,10)  / 
-     +    1,  1,  1,  1,  1,  1,  1,  1,  1,  1 /
-      DATA (NFLAGG(I,10),I=1,10)  / 
-     +    1,  1,  1,  1,  1,  1,  1,  1,  1,  1 /
-      DATA (A(I, 1, 1),I=1,10)  / 
-     +  .438D-01,.172D0  ,.283D0  ,.511D0  ,.715D0  ,.920D0  ,1.19D0  ,
-     +  1.37D0  ,1.65D0  ,2.14D0   /
-      DATA (A(I, 1, 2),I=1,10)  / 
-     +  .147D-01,.249D-01,.439D-01,.592D-01,.776D-01,.886D-01,.108D0  ,
-     +  .117D0  ,.126D0  ,.128D0   /
-      DATA (A(I, 1, 3),I=1,10)  / 
-     +  .216D-02,.627D-02,.834D-02,.108D-01,.144D-01,.152D-01,.196D-01,
-     +  .200D-01,.210D-01,.224D-01 /
-      DATA (A(I, 1, 4),I=1,10)  / 
-     +  .593D-01,.653D-01,.116D0  ,.145D0  ,.184D0  ,.204D0  ,.234D0  ,
-     +  .257D0  ,.271D0  ,.248D0   /
-      DATA (A(I, 1, 5),I=1,10)  / 
-     +  .000D+00,.918D-02,.362D-02,.805D-02,.436D-02,.728D-02,.466D-02,
-     +  .707D-02,.932D-02,.130D-01 /
-      DATA (A(I, 1, 6),I=1,10)  / 
-     +  .000D+00,.180D-02,.247D-02,.208D-02,.224D-02,.214D-02,.226D-02,
-     +  .233D-02,.230D-02,.194D-02 /
-      DATA (A(I, 1, 7),I=1,10)  / 
-     +  .000D+00,.106D-02,.703D-03,.687D-03,.739D-03,.674D-03,.819D-03,
-     +  .768D-03,.756D-03,.720D-03 /
-      DATA (A(I, 1, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,.188D-02,.130D-02,.138D-02,.117D-02,.124D-02,
-     +  .119D-02,.111D-02,.829D-03 /
-      DATA (A(I, 1, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,.302D-03,.258D-03,.249D-03,.208D-03,.248D-03,
-     +  .222D-03,.210D-03,.187D-03 /
-      DATA (A(I, 1,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.235D-03,.222D-03,.172D-03,.181D-03,
-     +  .166D-03,.152D-03,.124D-03 /
-      DATA (A(I, 1,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.238D-03,.179D-03,.145D-03,.156D-03,
-     +  .138D-03,.129D-03,.111D-03 /
-      DATA (A(I, 1,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.368D-03,.400D-03,.255D-03,.262D-03,
-     +  .221D-03,.182D-03,.112D-03 /
-      DATA (A(I, 1,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.753D-04,.712D-04,.527D-04,
-     +  .537D-04,.538D-04,.487D-04 /
-      DATA (A(I, 1,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.103D-03,.589D-04,.578D-04,
-     +  .468D-04,.385D-04,.269D-04 /
-      DATA (A(I, 1,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.444D-04,.372D-04,
-     +  .318D-04,.284D-04,.218D-04 /
-      DATA (A(I, 1,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.487D-04,.473D-04,
-     +  .338D-04,.243D-04,.122D-04 /
-      DATA (A(I, 1,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.121D-04,.117D-04,
-     +  .932D-05,.792D-05,.583D-05 /
-      DATA (A(I, 1,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.147D-04,
-     +  .101D-04,.756D-05,.496D-05 /
-      DATA (A(I, 1,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.755D-05,
-     +  .612D-05,.505D-05,.341D-05 /
-      DATA (A(I, 1,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  .630D-05,.444D-05,.282D-05 /
-      DATA (A(I, 2, 1),I=1,10)  / 
-     +  .269D0  ,.510D0  ,.738D0  ,1.12D0  ,1.46D0  ,1.83D0  ,2.22D0  ,
-     +  2.57D0  ,3.00D0  ,3.67D0   /
-      DATA (A(I, 2, 2),I=1,10)  / 
-     +  .121D0  ,.133D0  ,.190D0  ,.234D0  ,.293D0  ,.332D0  ,.395D0  ,
-     +  .431D0  ,.468D0  ,.502D0   /
-      DATA (A(I, 2, 3),I=1,10)  / 
-     +  .227D-01,.374D-01,.474D-01,.578D-01,.722D-01,.794D-01,.960D-01,
-     +  .102D0  ,.110D0  ,.120D0   /
-      DATA (A(I, 2, 4),I=1,10)  / 
-     +  .287D0  ,.196D0  ,.270D0  ,.314D0  ,.373D0  ,.408D0  ,.462D0  ,
-     +  .498D0  ,.529D0  ,.523D0   /
-      DATA (A(I, 2, 5),I=1,10)  / 
-     +  .000D+00,.433D-01,.218D-01,.384D-01,.263D-01,.385D-01,.298D-01,
-     +  .405D-01,.504D-01,.671D-01 /
-      DATA (A(I, 2, 6),I=1,10)  / 
-     +  .000D+00,.151D-01,.177D-01,.159D-01,.173D-01,.173D-01,.187D-01,
-     +  .196D-01,.201D-01,.191D-01 /
-      DATA (A(I, 2, 7),I=1,10)  / 
-     +  .000D+00,.457D-02,.607D-02,.610D-02,.677D-02,.670D-02,.784D-02,
-     +  .787D-02,.806D-02,.803D-02 /
-      DATA (A(I, 2, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,.702D-02,.536D-02,.558D-02,.510D-02,.554D-02,
-     +  .546D-02,.538D-02,.489D-02 /
-      DATA (A(I, 2, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,.190D-02,.199D-02,.205D-02,.191D-02,.221D-02,
-     +  .214D-02,.213D-02,.204D-02 /
-      DATA (A(I, 2,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.226D-02,.219D-02,.195D-02,.208D-02,
-     +  .204D-02,.203D-02,.194D-02 /
-      DATA (A(I, 2,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.213D-02,.195D-02,.175D-02,.191D-02,
-     +  .183D-02,.179D-02,.166D-02 /
-      DATA (A(I, 2,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.588D-03,.186D-02,.137D-02,.141D-02,
-     +  .128D-02,.117D-02,.947D-03 /
-      DATA (A(I, 2,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.554D-03,.562D-03,.454D-03,
-     +  .485D-03,.505D-03,.509D-03 /
-      DATA (A(I, 2,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.490D-03,.533D-03,.531D-03,
-     +  .476D-03,.437D-03,.369D-03 /
-      DATA (A(I, 2,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.427D-03,.382D-03,
-     +  .358D-03,.340D-03,.294D-03 /
-      DATA (A(I, 2,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.239D-03,.298D-03,
-     +  .238D-03,.196D-03,.134D-03 /
-      DATA (A(I, 2,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.299D-04,.893D-04,
-     +  .796D-04,.744D-04,.683D-04 /
-      DATA (A(I, 2,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.127D-03,
-     +  .107D-03,.916D-04,.720D-04 /
-      DATA (A(I, 2,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.397D-04,
-     +  .630D-04,.565D-04,.461D-04 /
-      DATA (A(I, 2,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  .511D-04,.459D-04,.402D-04 /
-      DATA (A(I, 3, 1),I=1,10)  / 
-     +  .708D0  ,1.02D0  ,1.41D0  ,1.91D0  ,2.42D0  ,3.00D0  ,3.53D0  ,
-     +  4.09D0  ,4.71D0  ,5.57D0   /
-      DATA (A(I, 3, 2),I=1,10)  / 
-     +  .397D0  ,.410D0  ,.539D0  ,.648D0  ,.795D0  ,.910D0  ,1.06D0  ,
-     +  1.17D0  ,1.29D0  ,1.42D0   /
-      DATA (A(I, 3, 3),I=1,10)  / 
-     +  .845D-01,.122D0  ,.157D0  ,.190D0  ,.232D0  ,.262D0  ,.307D0  ,
-     +  .335D0  ,.366D0  ,.402D0   /
-      DATA (A(I, 3, 4),I=1,10)  / 
-     +  .210D0  ,.379D0  ,.450D0  ,.490D0  ,.574D0  ,.636D0  ,.709D0  ,
-     +  .769D0  ,.820D0  ,.849D0   /
-      DATA (A(I, 3, 5),I=1,10)  / 
-     +  .000D+00,.102D0  ,.675D-01,.104D0  ,.858D-01,.115D0  ,.102D0  ,
-     +  .129D0  ,.154D0  ,.194D0   /
-      DATA (A(I, 3, 6),I=1,10)  / 
-     +  .000D+00,.392D-01,.615D-01,.593D-01,.649D-01,.674D-01,.735D-01,
-     +  .779D-01,.817D-01,.828D-01 /
-      DATA (A(I, 3, 7),I=1,10)  / 
-     +  .000D+00,.539D-02,.222D-01,.238D-01,.269D-01,.280D-01,.320D-01,
-     +  .334D-01,.350D-01,.361D-01 /
-      DATA (A(I, 3, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,.838D-02,.130D-01,.133D-01,.131D-01,.141D-01,
-     +  .144D-01,.149D-01,.152D-01 /
-      DATA (A(I, 3, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,.228D-02,.647D-02,.688D-02,.687D-02,.772D-02,
-     +  .786D-02,.811D-02,.824D-02 /
-      DATA (A(I, 3,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.664D-02,.828D-02,.802D-02,.845D-02,
-     +  .869D-02,.902D-02,.930D-02 /
-      DATA (A(I, 3,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.338D-02,.735D-02,.710D-02,.767D-02,
-     +  .767D-02,.776D-02,.756D-02 /
-      DATA (A(I, 3,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.280D-03,.262D-02,.349D-02,.342D-02,
-     +  .322D-02,.312D-02,.291D-02 /
-      DATA (A(I, 3,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.618D-03,.161D-02,.138D-02,
-     +  .148D-02,.155D-02,.166D-02 /
-      DATA (A(I, 3,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.313D-03,.128D-02,.161D-02,
-     +  .150D-02,.144D-02,.134D-02 /
-      DATA (A(I, 3,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.645D-03,.118D-02,
-     +  .115D-02,.111D-02,.103D-02 /
-      DATA (A(I, 3,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.117D-03,.497D-03,
-     +  .581D-03,.501D-03,.401D-03 /
-      DATA (A(I, 3,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.115D-04,.997D-04,
-     +  .202D-03,.203D-03,.206D-03 /
-      DATA (A(I, 3,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.877D-04,
-     +  .242D-03,.263D-03,.226D-03 /
-      DATA (A(I, 3,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.158D-04,
-     +  .881D-04,.152D-03,.136D-03 /
-      DATA (A(I, 3,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  .358D-04,.997D-04,.117D-03 /
-      DATA (A(I, 4, 1),I=1,10)  / 
-     +  .945D0  ,1.29D0  ,1.40D0  ,1.98D0  ,2.73D0  ,3.17D0  ,3.77D0  ,
-     +  4.29D0  ,4.78D0  ,5.54D0   /
-      DATA (A(I, 4, 2),I=1,10)  / 
-     +  .581D0  ,.599D0  ,.645D0  ,.839D0  ,1.10D0  ,1.25D0  ,1.47D0  ,
-     +  1.64D0  ,1.78D0  ,1.99D0   /
-      DATA (A(I, 4, 3),I=1,10)  / 
-     +  .127D0  ,.182D0  ,.202D0  ,.264D0  ,.344D0  ,.387D0  ,.455D0  ,
-     +  .504D0  ,.549D0  ,.611D0   /
-      DATA (A(I, 4, 4),I=1,10)  / 
-     +  .183D0  ,.464D0  ,.351D0  ,.444D0  ,.642D0  ,.659D0  ,.772D0  ,
-     +  .830D0  ,.882D0  ,.930D0   /
-      DATA (A(I, 4, 5),I=1,10)  / 
-     +  .000D+00,.122D0  ,.803D-01,.136D0  ,.134D0  ,.173D0  ,.164D0  ,
-     +  .203D0  ,.239D0  ,.300D0   /
-      DATA (A(I, 4, 6),I=1,10)  / 
-     +  .000D+00,.393D-01,.766D-01,.872D-01,.108D0  ,.111D0  ,.123D0  ,
-     +  .132D0  ,.139D0  ,.145D0   /
-      DATA (A(I, 4, 7),I=1,10)  / 
-     +  .000D+00,.416D-02,.289D-01,.360D-01,.454D-01,.477D-01,.549D-01,
-     +  .583D-01,.618D-01,.654D-01 /
-      DATA (A(I, 4, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,.761D-02,.157D-01,.214D-01,.205D-01,.233D-01,
-     +  .241D-01,.255D-01,.271D-01 /
-      DATA (A(I, 4, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,.238D-02,.803D-02,.123D-01,.123D-01,.140D-01,
-     +  .145D-01,.153D-01,.160D-01 /
-      DATA (A(I, 4,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.695D-02,.150D-01,.154D-01,.166D-01,
-     +  .172D-01,.181D-01,.192D-01 /
-      DATA (A(I, 4,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.355D-02,.104D-01,.143D-01,.156D-01,
-     +  .158D-01,.164D-01,.165D-01 /
-      DATA (A(I, 4,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.112D-03,.276D-02,.568D-02,.736D-02,
-     +  .684D-02,.691D-02,.661D-02 /
-      DATA (A(I, 4,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.740D-03,.222D-02,.339D-02,
-     +  .352D-02,.382D-02,.409D-02 /
-      DATA (A(I, 4,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.369D-03,.160D-02,.322D-02,
-     +  .375D-02,.375D-02,.355D-02 /
-      DATA (A(I, 4,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.750D-03,.190D-02,
-     +  .298D-02,.319D-02,.299D-02 /
-      DATA (A(I, 4,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.260D-03,.673D-03,
-     +  .117D-02,.156D-02,.126D-02 /
-      DATA (A(I, 4,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.283D-05,.131D-03,
-     +  .363D-03,.618D-03,.690D-03 /
-      DATA (A(I, 4,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.205D-03,
-     +  .378D-03,.709D-03,.844D-03 /
-      DATA (A(I, 4,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.654D-05,
-     +  .150D-03,.341D-03,.527D-03 /
-      DATA (A(I, 4,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  .957D-04,.197D-03,.406D-03 /
-      DATA (A(I, 5, 1),I=1,10)  / 
-     +  1.16D0  ,1.70D0  ,2.19D0  ,2.79D0  ,3.33D0  ,3.90D0  ,4.49D0  ,
-     +  5.07D0  ,5.66D0  ,6.38D0   /
-      DATA (A(I, 5, 2),I=1,10)  / 
-     +  .779D0  ,.899D0  ,1.09D0  ,1.28D0  ,1.51D0  ,1.71D0  ,1.96D0  ,
-     +  2.18D0  ,2.39D0  ,2.62D0   /
-      DATA (A(I, 5, 3),I=1,10)  / 
-     +  .167D0  ,.263D0  ,.334D0  ,.408D0  ,.482D0  ,.548D0  ,.632D0  ,
-     +  .700D0  ,.767D0  ,.840D0   /
-      DATA (A(I, 5, 4),I=1,10)  / 
-     +  .203D0  ,.565D0  ,.845D0  ,.867D0  ,.906D0  ,.961D0  ,1.08D0  ,
-     +  1.13D0  ,1.21D0  ,1.25D0   /
-      DATA (A(I, 5, 5),I=1,10)  / 
-     +  .000D+00,.129D0  ,.152D0  ,.237D0  ,.208D0  ,.268D0  ,.258D0  ,
-     +  .312D0  ,.368D0  ,.450D0   /
-      DATA (A(I, 5, 6),I=1,10)  / 
-     +  .000D+00,.460D-01,.126D0  ,.174D0  ,.182D0  ,.188D0  ,.208D0  ,
-     +  .219D0  ,.233D0  ,.239D0   /
-      DATA (A(I, 5, 7),I=1,10)  / 
-     +  .000D+00,.289D-02,.380D-01,.611D-01,.788D-01,.845D-01,.974D-01,
-     +  .103D0  ,.111D0  ,.117D0   /
-      DATA (A(I, 5, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,.137D-01,.223D-01,.374D-01,.436D-01,.488D-01,
-     +  .488D-01,.524D-01,.547D-01 /
-      DATA (A(I, 5, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,.162D-02,.114D-01,.198D-01,.263D-01,.315D-01,
-     +  .323D-01,.348D-01,.364D-01 /
-      DATA (A(I, 5,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.149D-01,.240D-01,.320D-01,.428D-01,
-     +  .436D-01,.469D-01,.493D-01 /
-      DATA (A(I, 5,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.562D-02,.194D-01,.290D-01,.408D-01,
-     +  .460D-01,.492D-01,.500D-01 /
-      DATA (A(I, 5,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.476D-04,.106D-01,.134D-01,.191D-01,
-     +  .227D-01,.264D-01,.253D-01 /
-      DATA (A(I, 5,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.281D-02,.679D-02,.879D-02,
-     +  .123D-01,.165D-01,.190D-01 /
-      DATA (A(I, 5,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.542D-04,.847D-02,.125D-01,
-     +  .144D-01,.173D-01,.192D-01 /
-      DATA (A(I, 5,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.652D-02,.982D-02,
-     +  .129D-01,.159D-01,.192D-01 /
-      DATA (A(I, 5,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.109D-03,.688D-02,
-     +  .751D-02,.845D-02,.905D-02 /
-      DATA (A(I, 5,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.823D-06,.237D-02,
-     +  .318D-02,.446D-02,.569D-02 /
-      DATA (A(I, 5,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.604D-03,
-     +  .610D-02,.673D-02,.827D-02 /
-      DATA (A(I, 5,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.716D-06,
-     +  .412D-02,.519D-02,.617D-02 /
-      DATA (A(I, 5,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  .710D-03,.543D-02,.674D-02 /
-      DATA (A(I, 6, 1),I=1,10)  / 
-     +  1.36D0  ,2.08D0  ,2.67D0  ,3.30D0  ,3.94D0  ,4.62D0  ,5.18D0  ,
-     +  3.60D0  ,3.64D0  ,3.95D0   /
-      DATA (A(I, 6, 2),I=1,10)  / 
-     +  1.07D0  ,1.33D0  ,1.58D0  ,1.82D0  ,2.10D0  ,2.44D0  ,2.74D0  ,
-     +  1.78D0  ,1.73D0  ,1.80D0   /
-      DATA (A(I, 6, 3),I=1,10)  / 
-     +  .158D0  ,.276D0  ,.402D0  ,.506D0  ,.609D0  ,.700D0  ,.802D0  ,
-     +  .638D0  ,.629D0  ,.658D0   /
-      DATA (A(I, 6, 4),I=1,10)  / 
-     +  .308D0  ,.739D0  ,1.02D0  ,1.12D0  ,1.26D0  ,1.35D0  ,1.57D0  ,
-     +  1.94D0  ,1.71D0  ,1.55D0   /
-      DATA (A(I, 6, 5),I=1,10)  / 
-     +  .000D+00,.217D0  ,.183D0  ,.324D0  ,.276D0  ,.395D0  ,.393D0  ,
-     +  .558D0  ,.602D0  ,.681D0   /
-      DATA (A(I, 6, 6),I=1,10)  / 
-     +  .000D+00,.658D-01,.251D0  ,.267D0  ,.299D0  ,.326D0  ,.386D0  ,
-     +  .452D0  ,.475D0  ,.409D0   /
-      DATA (A(I, 6, 7),I=1,10)  / 
-     +  .000D+00,.198D-02,.774D-01,.136D0  ,.149D0  ,.164D0  ,.187D0  ,
-     +  .210D0  ,.238D0  ,.256D0   /
-      DATA (A(I, 6, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,.290D-01,.122D0  ,.139D0  ,.128D0  ,.129D0  ,
-     +  .137D0  ,.147D0  ,.167D0   /
-      DATA (A(I, 6, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,.699D-03,.617D-01,.750D-01,.801D-01,.905D-01,
-     +  .974D-01,.105D0  ,.122D0   /
-      DATA (A(I, 6,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.310D-01,.112D0  ,.127D0  ,.140D0  ,
-     +  .143D0  ,.155D0  ,.176D0   /
-      DATA (A(I, 6,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.277D-02,.889D-01,.143D0  ,.150D0  ,
-     +  .175D0  ,.184D0  ,.208D0   /
-      DATA (A(I, 6,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.202D-04,.343D-01,.959D-01,.109D0  ,
-     +  .115D0  ,.112D0  ,.116D0   /
-      DATA (A(I, 6,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.186D-02,.435D-01,.512D-01,
-     +  .744D-01,.856D-01,.103D0   /
-      DATA (A(I, 6,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.144D-04,.427D-01,.786D-01,
-     +  .911D-01,.993D-01,.108D0   /
-      DATA (A(I, 6,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.466D-02,.518D-01,
-     +  .848D-01,.109D0  ,.119D0   /
-      DATA (A(I, 6,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.655D-05,.330D-01,
-     +  .586D-01,.617D-01,.594D-01 /
-      DATA (A(I, 6,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.228D-06,.328D-02,
-     +  .190D-01,.301D-01,.454D-01 /
-      DATA (A(I, 6,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.218D-04,
-     +  .272D-01,.501D-01,.707D-01 /
-      DATA (A(I, 6,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.146D-06,
-     +  .441D-02,.378D-01,.556D-01 /
-      DATA (A(I, 6,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  .160D-03,.204D-01,.679D-01 /
-      DATA (A(I, 7, 1),I=1,10)  / 
-     +  .522D0  ,.862D0  ,1.14D0  ,1.40D0  ,1.70D0  ,1.94D0  ,2.26D0  ,
-     +  2.48D0  ,2.72D0  ,3.95D0   /
-      DATA (A(I, 7, 2),I=1,10)  / 
-     +  .314D0  ,.450D0  ,.588D0  ,.692D0  ,.834D0  ,.936D0  ,1.09D0  ,
-     +  1.18D0  ,1.28D0  ,1.80D0   /
-      DATA (A(I, 7, 3),I=1,10)  / 
-     +  .814D-01,.147D0  ,.189D0  ,.226D0  ,.272D0  ,.302D0  ,.351D0  ,
-     +  .378D0  ,.406D0  ,.658D0   /
-      DATA (A(I, 7, 4),I=1,10)  / 
-     +  .252D0  ,.864D0  ,1.01D0  ,.851D0  ,.837D0  ,.774D0  ,.763D0  ,
-     +  .757D0  ,.748D0  ,1.55D0   /
-      DATA (A(I, 7, 5),I=1,10)  / 
-     +  .000D+00,.225D0  ,.180D0  ,.276D0  ,.193D0  ,.240D0  ,.190D0  ,
-     +  .228D0  ,.259D0  ,.681D0   /
-      DATA (A(I, 7, 6),I=1,10)  / 
-     +  .000D+00,.485D-01,.272D0  ,.273D0  ,.253D0  ,.216D0  ,.206D0  ,
-     +  .197D0  ,.191D0  ,.409D0   /
-      DATA (A(I, 7, 7),I=1,10)  / 
-     +  .000D+00,.137D-02,.752D-01,.137D0  ,.152D0  ,.134D0  ,.125D0  ,
-     +  .119D0  ,.116D0  ,.256D0   /
-      DATA (A(I, 7, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,.220D-01,.155D0  ,.175D0  ,.155D0  ,.116D0  ,
-     +  .977D-01,.858D-01,.167D0   /
-      DATA (A(I, 7, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,.326D-03,.695D-01,.881D-01,.106D0  ,.897D-01,
-     +  .782D-01,.706D-01,.122D0   /
-      DATA (A(I, 7,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.261D-01,.124D0  ,.131D0  ,.156D0  ,
-     +  .141D0  ,.121D0  ,.176D0   /
-      DATA (A(I, 7,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.785D-03,.864D-01,.130D0  ,.170D0  ,
-     +  .182D0  ,.172D0  ,.208     /
-      DATA (A(I, 7,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.896D-05,.225D-01,.105D0  ,.126D0  ,
-     +  .126D0  ,.135D0  ,.116D0   /
-      DATA (A(I, 7,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.542D-03,.427D-01,.553D-01,
-     +  .744D-01,.980D-01,.103D0   /
-      DATA (A(I, 7,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.515D-05,.377D-01,.831D-01,
-     +  .985D-01,.104D0  ,.108D0   /
-      DATA (A(I, 7,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.285D-02,.495D-01,
-     +  .871D-01,.106D0  ,.119D0   /
-      DATA (A(I, 7,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.110D-05,.284D-01,
-     +  .588D-01,.657D-01,.594D-01 /
-      DATA (A(I, 7,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.722D-07,.176D-02,
-     +  .170D-01,.305D-01,.454D-01 /
-      DATA (A(I, 7,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.148D-05,
-     +  .213D-01,.492D-01,.707D-01 /
-      DATA (A(I, 7,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.323D-07,
-     +  .722D-02,.359D-01,.556D-01 /
-      DATA (A(I, 7,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  .461D-05,.155D-01,.679D-01 /
-      DATA (A(I, 8, 1),I=1,10)  / 
-     +  .630D0  ,.974D0  ,1.29D0  ,1.58D0  ,1.89D0  ,2.16D0  ,2.49D0  ,
-     +  2.75D0  ,3.02D0  ,3.95D0   /
-      DATA (A(I, 8, 2),I=1,10)  / 
-     +  .328D0  ,.459D0  ,.613D0  ,.735D0  ,.879D0  ,.994D0  ,1.15D0  ,
-     +  1.27D0  ,1.38D0  ,1.80D0   /
-      DATA (A(I, 8, 3),I=1,10)  / 
-     +  .748D-01,.121D0  ,.164D0  ,.197D0  ,.235D0  ,.265D0  ,.310D0  ,
-     +  .339D0  ,.370D0  ,.658D0   /
-      DATA (A(I, 8, 4),I=1,10)  / 
-     +  .194D0  ,.211D0  ,.337D0  ,.344D0  ,.339D0  ,.351D0  ,.390    ,
-     +  .419D0  ,.442D0  ,1.55D0   /
-      DATA (A(I, 8, 5),I=1,10)  / 
-     +  .000D+00,.869D-01,.725D-01,.113D0  ,.810D-01,.106D0  ,.951D-01,
-     +  .120D0  ,.143D0  ,.681D0   /
-      DATA (A(I, 8, 6),I=1,10)  / 
-     +  .000D+00,.288D-01,.102D0  ,.922D-01,.857D-01,.845D-01,.932D-01,
-     +  .983D-01,.102D0  ,.409D0   /
-      DATA (A(I, 8, 7),I=1,10)  / 
-     +  .000D+00,.668D-03,.533D-01,.575D-01,.493D-01,.482D-01,.539D-01,
-     +  .558D-01,.582D-01,.256D0   /
-      DATA (A(I, 8, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,.205D-01,.808D-01,.510D-01,.409D-01,.406D-01,
-     +  .394D-01,.389D-01,.167D0   /
-      DATA (A(I, 8, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,.999D-04,.647D-01,.385D-01,.325D-01,.325D-01,
-     +  .316D-01,.314D-01,.122D0   /
-      DATA (A(I, 8,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.169D-01,.834D-01,.611D-01,.565D-01,
-     +  .533D-01,.519D-01,.176D0   /
-      DATA (A(I, 8,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.107D-03,.769D-01,.922D-01,.805D-01,
-     +  .745D-01,.711D-01,.208D0   /
-      DATA (A(I, 8,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.180D-05,.143D-01,.983D-01,.775D-01,
-     +  .627D-01,.541D-01,.116D0   /
-      DATA (A(I, 8,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.157D-04,.346D-01,.507D-01,
-     +  .479D-01,.455D-01,.103D0   /
-      DATA (A(I, 8,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.752D-06,.248D-01,.721D-01,
-     +  .728D-01,.611D-01,.108D0   /
-      DATA (A(I, 8,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.686D-04,.356D-01,
-     +  .731D-01,.791D-01,.119D0   /
-      DATA (A(I, 8,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.838D-07,.151D-01,
-     +  .470D-01,.567D-01,.594D-01 /
-      DATA (A(I, 8,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.759D-08,.400D-04,
-     +  .193D-01,.313D-01,.454D-01 /
-      DATA (A(I, 8,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.385D-07,
-     +  .921D-02,.353D-01,.707D-01 /
-      DATA (A(I, 8,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.219D-08,
-     +  .348D-03,.226D-01,.556D-01 /
-      DATA (A(I, 8,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  .212D-07,.149D-01,.679D-01 /
-      DATA (A(I, 9, 1),I=1,10)  / 
-     +  .736D0  ,1.13D0  ,1.49D0  ,1.82D0  ,2.20D0  ,2.49D0  ,2.86D0  ,
-     +  3.17D0  ,3.49D0  ,3.95D0   /
-      DATA (A(I, 9, 2),I=1,10)  / 
-     +  .339D0  ,.492D0  ,.658D0  ,.789D0  ,.958D0  ,1.08D0  ,1.25D0  ,
-     +  1.37D0  ,1.50D0  ,1.80D0   /
-      DATA (A(I, 9, 3),I=1,10)  / 
-     +  .680D-01,.110D0  ,.150D0  ,.180D0  ,.222D0  ,.247D0  ,.289    ,
-     +  .318D0  ,.349D0  ,.658D0   /
-      DATA (A(I, 9, 4),I=1,10)  / 
-     +  .110D0  ,.104D0  ,.157D0  ,.156D0  ,.210D0  ,.205D0  ,.246D0  ,
-     +  .274D0  ,.300D0  ,1.55D0   /
-      DATA (A(I, 9, 5),I=1,10)  / 
-     +  .000D+00,.379D-01,.347D-01,.477D-01,.486D-01,.576D-01,.569D-01,
-     +  .732D-01,.893D-01,.681D0   /
-      DATA (A(I, 9, 6),I=1,10)  / 
-     +  .000D+00,.223D-01,.354D-01,.312D-01,.436D-01,.400D-01,.489D-01,
-     +  .548D-01,.600D-01,.409D0   /
-      DATA (A(I, 9, 7),I=1,10)  / 
-     +  .000D+00,.338D-03,.149D-01,.142D-01,.215D-01,.188D-01,.248D-01,
-     +  .278D-01,.307D-01,.256D0   /
-      DATA (A(I, 9, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,.553D-02,.862D-02,.150D-01,.106D-01,.145D-01,
-     +  .165D-01,.181D-01,.167D0   /
-      DATA (A(I, 9, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,.375D-04,.641D-02,.111D-01,.792D-02,.112D-01,
-     +  .127D-01,.140D-01,.122D0   /
-      DATA (A(I, 9,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.112D-01,.200D-01,.127D-01,.176D-01,
-     +  .200D-01,.220D-01,.176D0   /
-      DATA (A(I, 9,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.244D-04,.261D-01,.162D-01,.232D-01,
-     +  .263D-01,.287D-01,.208D0   /
-      DATA (A(I, 9,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.455D-06,.635D-02,.121D-01,.186D-01,
-     +  .201D-01,.207D-01,.116D0   /
-      DATA (A(I, 9,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.146D-05,.922D-02,.116D-01,
-     +  .145D-01,.165D-01,.103D0   /
-      DATA (A(I, 9,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.135D-06,.128D-01,.202D-01,
-     +  .215D-01,.220D-01,.108D0   /
-      DATA (A(I, 9,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.237D-05,.229D-01,
-     +  .259D-01,.271D-01,.119D0   /
-      DATA (A(I, 9,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.100D-07,.534D-02,
-     +  .210D-01,.193D-01,.594D-01 /
-      DATA (A(I, 9,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.915D-09,.847D-06,
-     +  .119D-01,.125D-01,.454D-01 /
-      DATA (A(I, 9,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.298D-08,
-     +  .101D-01,.242D-01,.707D-01 /
-      DATA (A(I, 9,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.196D-09,
-     +  .243D-05,.234D-01,.556D-01 /
-      DATA (A(I, 9,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  .575D-09,.364D-02,.679D-01 /
-      DATA (A(I,10, 1),I=1,10)  / 
-     +  .959D0  ,1.46D0  ,1.92D0  ,2.34D0  ,2.80D0  ,3.24D0  ,3.64D0  ,
-     +  4.05D0  ,4.48D0  ,3.95D0   /
-      DATA (A(I,10, 2),I=1,10)  / 
-     +  .343D0  ,.516D0  ,.692D0  ,.836D0  ,1.01D0  ,1.16D0  ,1.31D0  ,
-     +  1.46D0  ,1.61D0  ,1.80D0   /
-      DATA (A(I,10, 3),I=1,10)  / 
-     +  .512D-01,.837D-01,.115D0  ,.138D0  ,.169D0  ,.195D0  ,.220D0  ,
-     +  .245D0  ,.270D0  ,.658D0   /
-      DATA (A(I,10, 4),I=1,10)  / 
-     +  .274D-01,.361D-01,.510D-01,.562D-01,.703D-01,.828D-01,.877D-01,
-     +  .996D-01,.111D0  ,1.55D0   /
-      DATA (A(I,10, 5),I=1,10)  / 
-     +  .000D+00,.850D-02,.875D-02,.118D-01,.124D-01,.170D-01,.154D-01,
-     +  .194D-01,.237D-01,.681D0   /
-      DATA (A(I,10, 6),I=1,10)  / 
-     +  .000D+00,.345D-02,.519D-02,.533D-02,.691D-02,.842D-02,.844D-02,
-     +  .987D-02,.113D-01,.409D0   /
-      DATA (A(I,10, 7),I=1,10)  / 
-     +  .000D+00,.722D-04,.130D-02,.135D-02,.189D-02,.240D-02,.235D-02,
-     +  .281D-02,.331D-02,.256D0   /
-      DATA (A(I,10, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,.283D-03,.272D-03,.394D-03,.557D-03,.480D-03,
-     +  .616D-03,.775D-03,.167D0   /
-      DATA (A(I,10, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,.457D-05,.122D-03,.192D-03,.275D-03,.225D-03,
-     +  .292D-03,.373D-03,.122D0   /
-      DATA (A(I,10,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.119D-03,.185D-03,.278D-03,.201D-03,
-     +  .274D-03,.364D-03,.176D0   /
-      DATA (A(I,10,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.140D-05,.129D-03,.200D-03,.137D-03,
-     +  .188D-03,.252D-03,.208D0   /
-      DATA (A(I,10,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.207D-07,.307D-04,.518D-04,.278D-04,
-     +  .421D-04,.608D-04,.116D0   /
-      DATA (A(I,10,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.306D-07,.252D-04,.111D-04,
-     +  .188D-04,.295D-04,.103D0   /
-      DATA (A(I,10,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.321D-08,.220D-04,.104D-04,
-     +  .162D-04,.243D-04,.108D0   /
-      DATA (A(I,10,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.770D-08,.632D-05,
-     +  .105D-04,.162D-04,.119D0   /
-      DATA (A(I,10,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.117D-09,.199D-05,
-     +  .321D-05,.492D-05,.594D-01 /
-      DATA (A(I,10,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.888E-11,.323D-09,
-     +  .106D-05,.192D-05,.454D-01 /
-      DATA (A(I,10,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.174E-10,
-     +  .131D-05,.218D-05,.707D-01 /
-      DATA (A(I,10,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.994E-12,
-     +  .233D-09,.104D-05,.556D-01 /
-      DATA (A(I,10,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  .144E-11,.724D-06,.679D-01 /
-      DATA (AE(I, 1, 1),I=1,10)  / 
-     +  7.27D0  ,6.29D0  ,7.76D0  ,6.70D0  ,8.17D0  ,7.34D0  ,8.70D0  ,
-     +  8.02D0  ,7.37D0  ,6.18D0   /
-      DATA (AE(I, 1, 2),I=1,10)  / 
-     +  7.41D0  ,7.52D0  ,8.14D0  ,8.20D0  ,8.96D0  ,9.05D0  ,9.96D0  ,
-     +  10.0D0  ,10.1D0  ,9.86D0   /
-      DATA (AE(I, 1, 3),I=1,10)  / 
-     +  7.72D0  ,7.69D0  ,9.17D0  ,8.99D0  ,10.6D0  ,10.5D0  ,12.1D0  ,
-     +  12.1D0  ,12.0D0  ,11.5D0   /
-      DATA (AE(I, 1, 4),I=1,10)  / 
-     +  7.90D0  ,8.48D0  ,9.50D0  ,9.94D0  ,10.8D0  ,11.4D0  ,12.2D0  ,
-     +  12.8D0  ,13.3D0  ,13.8D0   /
-      DATA (AE(I, 1, 5),I=1,10)  / 
-     +  .000D+00,8.52D0  ,9.59D0  ,10.1D0  ,11.1D0  ,11.8D0  ,12.7D0  ,
-     +  13.3D0  ,13.8D0  ,14.4D0   /
-      DATA (AE(I, 1, 6),I=1,10)  / 
-     +  .000D+00,9.00D0  ,10.7D0  ,11.7D0  ,13.2D0  ,14.2D0  ,15.6D0  ,
-     +  16.5D0  ,17.3D0  ,18.0D0   /
-      DATA (AE(I, 1, 7),I=1,10)  / 
-     +  .000D+00,9.01D0  ,11.1D0  ,11.9D0  ,14.3D0  ,15.0D0  ,17.4D0  ,
-     +  18.0D0  ,18.6D0  ,18.8D0   /
-      DATA (AE(I, 1, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,11.2D0  ,12.4D0  ,14.5D0  ,15.7D0  ,17.6D0  ,
-     +  18.8D0  ,19.9D0  ,20.9D0   /
-      DATA (AE(I, 1, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,11.4D0  ,12.7D0  ,15.5D0  ,16.6D0  ,19.3D0  ,
-     +  20.2D0  ,21.1D0  ,21.7D0   /
-      DATA (AE(I, 1,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,13.2D0  ,15.8D0  ,17.3D0  ,19.9D0  ,
-     +  21.2D0  ,22.4D0  ,23.2D0   /
-      DATA (AE(I, 1,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,13.2D0  ,16.3D0  ,17.8D0  ,20.8D0  ,
-     +  22.1D0  ,23.3D0  ,24.2D0   /
-      DATA (AE(I, 1,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,13.4D0  ,16.2D0  ,18.2D0  ,21.0D0  ,
-     +  22.8D0  ,24.4D0  ,25.9D0   /
-      DATA (AE(I, 1,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,16.5D0  ,18.4D0  ,21.6D0  ,
-     +  23.2D0  ,24.8D0  ,26.2D0   /
-      DATA (AE(I, 1,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,16.7D0  ,19.0D0  ,22.3D0  ,
-     +  24.3D0  ,26.1D0  ,27.4D0   /
-      DATA (AE(I, 1,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,19.1D0  ,22.8D0  ,
-     +  24.7D0  ,26.6D0  ,28.2D0   /
-      DATA (AE(I, 1,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,19.2D0  ,23.0D0  ,
-     +  25.3D0  ,27.5D0  ,29.5D0   /
-      DATA (AE(I, 1,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,19.6D0  ,23.3D0  ,
-     +  25.6D0  ,27.8D0  ,29.6D0   /
-      DATA (AE(I, 1,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,23.6D0  ,
-     +  26.2D0  ,28.5D0  ,30.4D0   /
-      DATA (AE(I, 1,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,23.7D0  ,
-     +  26.3D0  ,28.8D0  ,31.0D0   /
-      DATA (AE(I, 1,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  26.5D0  ,29.2D0  ,31.5D0   /
-      DATA (AE(I, 2, 1),I=1,10)  / 
-     +  8.74D0  ,8.16D0  ,9.25D0  ,8.45D0  ,9.46D0  ,8.90D0  ,9.83D0  ,
-     +  9.38D0  ,8.96D0  ,8.15D0   /
-      DATA (AE(I, 2, 2),I=1,10)  / 
-     +  8.96D0  ,9.30D0  ,9.95D0  ,10.0D0  ,10.8D0  ,10.9D0  ,11.7D0  ,
-     +  11.8D0  ,11.9D0  ,11.8D0   /
-      DATA (AE(I, 2, 3),I=1,10)  / 
-     +  9.44D0  ,9.66D0  ,11.0D0  ,11.0D0  ,12.3D0  ,12.5D0  ,13.7D0  ,
-     +  13.9D0  ,14.0D0  ,13.8D0   /
-      DATA (AE(I, 2, 4),I=1,10)  / 
-     +  8.86D0  ,9.81D0  ,10.8D0  ,11.2D0  ,12.0D0  ,12.6D0  ,13.4D0  ,
-     +  14.0D0  ,14.5D0  ,15.1D0   /
-      DATA (AE(I, 2, 5),I=1,10)  / 
-     +  .000D+00,10.2D0  ,11.4D0  ,12.0D0  ,12.9D0  ,13.6D0  ,14.5D0  ,
-     +  15.1D0  ,15.7D0  ,16.3D0   /
-      DATA (AE(I, 2, 6),I=1,10)  / 
-     +  .000D+00,10.7D0  ,12.5D0  ,13.5D0  ,15.1D0  ,16.0D0  ,17.5D0  ,
-     +  18.3D0  ,19.2D0  ,19.9D0   /
-      DATA (AE(I, 2, 7),I=1,10)  / 
-     +  .000D+00,11.5D0  ,12.9D0  ,13.9D0  ,16.1D0  ,17.0D0  ,19.1D0  ,
-     +  19.8D0  ,20.6D0  ,21.0D0   /
-      DATA (AE(I, 2, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,12.4D0  ,13.8D0  ,15.9D0  ,17.2D0  ,19.1D0  ,
-     +  20.3D0  ,21.4D0  ,22.3D0   /
-      DATA (AE(I, 2, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,13.4D0  ,14.5D0  ,17.1D0  ,18.3D0  ,20.9D0  ,
-     +  21.9D0  ,23.0D0  ,23.7D0   /
-      DATA (AE(I, 2,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,14.9D0  ,17.5D0  ,19.1D0  ,21.6D0  ,
-     +  22.9D0  ,24.1D0  ,25.0D0   /
-      DATA (AE(I, 2,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,15.0D0  ,18.0D0  ,19.6D0  ,22.4D0  ,
-     +  23.8D0  ,25.2D0  ,26.2D0   /
-      DATA (AE(I, 2,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,16.2D0  ,17.3D0  ,19.4D0  ,22.2D0  ,
-     +  24.0D0  ,25.7D0  ,27.2D0   /
-      DATA (AE(I, 2,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,17.8D0  ,19.8D0  ,22.9D0  ,
-     +  24.6D0  ,26.2D0  ,27.7D0   /
-      DATA (AE(I, 2,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,19.1D0  ,20.4D0  ,23.7D0  ,
-     +  25.7D0  ,27.6D0  ,29.1D0   /
-      DATA (AE(I, 2,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,20.5D0  ,24.1D0  ,
-     +  26.1D0  ,28.1D0  ,29.9D0   /
-      DATA (AE(I, 2,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,20.9D0  ,23.9D0  ,
-     +  26.4D0  ,28.7D0  ,30.7D0   /
-      DATA (AE(I, 2,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,22.4D0  ,24.2D0  ,
-     +  26.7D0  ,29.0D0  ,30.9D0   /
-      DATA (AE(I, 2,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,24.8D0  ,
-     +  27.3D0  ,29.7D0  ,31.8D0   /
-      DATA (AE(I, 2,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,26.1D0  ,
-     +  27.3D0  ,29.9D0  ,32.3D0   /
-      DATA (AE(I, 2,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  27.4D0  ,30.1D0  ,32.6D0   /
-      DATA (AE(I, 3, 1),I=1,10)  / 
-     +  11.0D0  ,11.0D0  ,11.7D0  ,11.3D0  ,11.9D0  ,11.4D0  ,12.1D0  ,
-     +  11.7D0  ,11.5D0  ,11.0D0   /
-      DATA (AE(I, 3, 2),I=1,10)  / 
-     +  11.2D0  ,12.0D0  ,12.7D0  ,12.9D0  ,13.6D0  ,13.7D0  ,14.4D0  ,
-     +  14.6D0  ,14.7D0  ,14.6D0   /
-      DATA (AE(I, 3, 3),I=1,10)  / 
-     +  12.1D0  ,12.6D0  ,13.7D0  ,13.9D0  ,15.0D0  ,15.2D0  ,16.3D0  ,
-     +  16.5D0  ,16.7D0  ,16.7D0   /
-      DATA (AE(I, 3, 4),I=1,10)  / 
-     +  12.6D0  ,11.3D0  ,12.4D0  ,13.0D0  ,13.8D0  ,14.2D0  ,15.0D0  ,
-     +  15.6D0  ,16.1D0  ,16.6D0   /
-      DATA (AE(I, 3, 5),I=1,10)  / 
-     +  .000D+00,12.6D0  ,13.7D0  ,14.4D0  ,15.3D0  ,16.0D0  ,16.8D0  ,
-     +  17.5D0  ,18.1D0  ,18.6D0   /
-      DATA (AE(I, 3, 6),I=1,10)  / 
-     +  .000D+00,14.0D0  ,14.6D0  ,15.8D0  ,17.4D0  ,18.4D0  ,19.8D0  ,
-     +  20.6D0  ,21.5D0  ,22.2D0   /
-      DATA (AE(I, 3, 7),I=1,10)  / 
-     +  .000D+00,16.0D0  ,15.2D0  ,16.3D0  ,18.3D0  ,19.3D0  ,21.1D0  ,
-     +  22.0D0  ,22.8D0  ,23.5D0   /
-      DATA (AE(I, 3, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,15.6D0  ,15.1D0  ,17.2D0  ,18.6D0  ,20.6D0  ,
-     +  21.8D0  ,22.9D0  ,23.8D0   /
-      DATA (AE(I, 3, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,17.8D0  ,16.3D0  ,18.8D0  ,20.1D0  ,22.5D0  ,
-     +  23.6D0  ,24.7D0  ,25.6D0   /
-      DATA (AE(I, 3,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,17.5D0  ,19.0D0  ,20.7D0  ,23.1D0  ,
-     +  24.5D0  ,25.8D0  ,26.8D0   /
-      DATA (AE(I, 3,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,19.2D0  ,19.4D0  ,21.1D0  ,23.8D0  ,
-     +  25.4D0  ,26.8D0  ,28.0D0   /
-      DATA (AE(I, 3,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,20.7D0  ,19.6D0  ,19.7D0  ,22.4D0  ,
-     +  24.4D0  ,26.2D0  ,27.9D0   /
-      DATA (AE(I, 3,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,21.6D0  ,20.4D0  ,23.2D0  ,
-     +  25.1D0  ,26.9D0  ,28.5D0   /
-      DATA (AE(I, 3,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,23.5D0  ,22.0D0  ,23.8D0  ,
-     +  26.1D0  ,28.1D0  ,29.9D0   /
-      DATA (AE(I, 3,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,23.7D0  ,24.2D0  ,
-     +  26.3D0  ,28.5D0  ,30.4D0   /
-      DATA (AE(I, 3,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,25.4D0  ,24.8D0  ,
-     +  25.6D0  ,28.1D0  ,30.5D0   /
-      DATA (AE(I, 3,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,26.9D0  ,26.8D0  ,
-     +  26.1D0  ,28.4D0  ,30.8D0   /
-      DATA (AE(I, 3,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,28.8D0  ,
-     +  27.6D0  ,29.0D0  ,31.5D0   /
-      DATA (AE(I, 3,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,30.5D0  ,
-     +  29.2D0  ,28.9D0  ,31.5D0   /
-      DATA (AE(I, 3,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  31.0D0  ,30.0D0  ,31.7D0   /
-      DATA (AE(I, 4, 1),I=1,10)  / 
-     +  13.0D0  ,13.2D0  ,14.8D0  ,14.2D0  ,14.2D0  ,14.1D0  ,14.5D0  ,
-     +  14.4D0  ,14.3D0  ,14.0D0   /
-      DATA (AE(I, 4, 2),I=1,10)  / 
-     +  13.5D0  ,14.5D0  ,16.1D0  ,15.9D0  ,16.0D0  ,16.3D0  ,16.8D0  ,
-     +  17.0D0  ,17.1D0  ,17.2D0   /
-      DATA (AE(I, 4, 3),I=1,10)  / 
-     +  14.9D0  ,15.3D0  ,17.2D0  ,17.1D0  ,17.5D0  ,17.8D0  ,18.6D0  ,
-     +  18.9D0  ,19.1D0  ,19.3D0   /
-      DATA (AE(I, 4, 4),I=1,10)  / 
-     +  15.1D0  ,13.5D0  ,16.4D0  ,16.7D0  ,16.4D0  ,17.3D0  ,17.8D0  ,
-     +  18.5D0  ,19.0D0  ,19.6D0   /
-      DATA (AE(I, 4, 5),I=1,10)  / 
-     +  .000D+00,15.6D0  ,17.5D0  ,17.7D0  ,17.8D0  ,18.6D0  ,19.2D0  ,
-     +  19.9D0  ,20.3D0  ,21.1D0   /
-      DATA (AE(I, 4, 6),I=1,10)  / 
-     +  .000D+00,18.0D0  ,18.4D0  ,19.2D0  ,19.8D0  ,20.9D0  ,22.0D0  ,
-     +  23.1D0  ,23.6D0  ,24.7D0   /
-      DATA (AE(I, 4, 7),I=1,10)  / 
-     +  .000D+00,27.4D0  ,19.1D0  ,19.8D0  ,20.7D0  ,21.8D0  ,23.2D0  ,
-     +  24.4D0  ,24.9D0  ,25.9D0   /
-      DATA (AE(I, 4, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,18.9D0  ,18.9D0  ,19.3D0  ,21.1D0  ,22.5D0  ,
-     +  24.0D0  ,24.7D0  ,26.0D0   /
-      DATA (AE(I, 4, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,21.1D0  ,19.7D0  ,20.7D0  ,22.3D0  ,24.0D0  ,
-     +  25.6D0  ,26.3D0  ,27.7D0   /
-      DATA (AE(I, 4,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,21.0D0  ,21.1D0  ,22.9D0  ,24.6D0  ,
-     +  26.5D0  ,27.3D0  ,29.0D0   /
-      DATA (AE(I, 4,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,21.3D0  ,22.4D0  ,23.1D0  ,25.0D0  ,
-     +  27.1D0  ,27.9D0  ,29.8D0   /
-      DATA (AE(I, 4,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,36.6D0  ,21.5D0  ,22.2D0  ,23.1D0  ,
-     +  25.6D0  ,26.8D0  ,29.1D0   /
-      DATA (AE(I, 4,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,22.9D0  ,23.1D0  ,23.7D0  ,
-     +  26.2D0  ,27.3D0  ,29.6D0   /
-      DATA (AE(I, 4,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,30.5D0  ,23.6D0  ,25.0D0  ,
-     +  26.9D0  ,28.2D0  ,30.7D0   /
-      DATA (AE(I, 4,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,25.4D0  ,26.2D0  ,
-     +  27.2D0  ,28.3D0  ,31.0D0   /
-      DATA (AE(I, 4,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,24.5D0  ,25.9D0  ,
-     +  27.4D0  ,27.6D0  ,30.7D0   /
-      DATA (AE(I, 4,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,43.3D0  ,28.4D0  ,
-     +  27.5D0  ,27.9D0  ,30.9D0   /
-      DATA (AE(I, 4,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,27.2D0  ,
-     +  29.1D0  ,29.0D0  ,31.4D0   /
-      DATA (AE(I, 4,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,51.3D0  ,
-     +  30.6D0  ,29.5D0  ,31.4D0   /
-      DATA (AE(I, 4,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  28.8D0  ,30.6D0  ,32.4D0   /
-      DATA (AE(I, 5, 1),I=1,10)  / 
-     +  15.0D0  ,14.9D0  ,15.5D0  ,15.4D0  ,15.9D0  ,15.8D0  ,16.2D0  ,
-     +  16.2D0  ,16.1D0  ,15.9D0   /
-      DATA (AE(I, 5, 2),I=1,10)  / 
-     +  15.4D0  ,16.1D0  ,17.0D0  ,17.4D0  ,18.0D0  ,18.2D0  ,18.7D0  ,
-     +  18.9D0  ,19.0D0  ,19.1D0   /
-      DATA (AE(I, 5, 3),I=1,10)  / 
-     +  17.1D0  ,17.2D0  ,18.3D0  ,18.7D0  ,19.3D0  ,19.6D0  ,20.3D0  ,
-     +  20.6D0  ,20.8D0  ,20.9D0   /
-      DATA (AE(I, 5, 4),I=1,10)  / 
-     +  14.7D0  ,14.8D0  ,15.0D0  ,16.0D0  ,17.0D0  ,17.7D0  ,18.1D0  ,
-     +  19.0D0  ,19.4D0  ,20.0D0   /
-      DATA (AE(I, 5, 5),I=1,10)  / 
-     +  .000D+00,16.7D0  ,17.6D0  ,18.1D0  ,18.6D0  ,19.2D0  ,19.7D0  ,
-     +  20.4D0  ,20.8D0  ,21.2D0   /
-      DATA (AE(I, 5, 6),I=1,10)  / 
-     +  .000D+00,17.8D0  ,18.2D0  ,19.2D0  ,20.0D0  ,21.0D0  ,21.9D0  ,
-     +  23.0D0  ,23.6D0  ,24.3D0   /
-      DATA (AE(I, 5, 7),I=1,10)  / 
-     +  .000D+00,35.2D0  ,18.9D0  ,20.3D0  ,20.6D0  ,21.5D0  ,22.6D0  ,
-     +  23.7D0  ,24.2D0  ,24.7D0   /
-      DATA (AE(I, 5, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,16.4D0  ,18.9D0  ,18.8D0  ,19.6D0  ,20.7D0  ,
-     +  22.3D0  ,23.1D0  ,23.9D0   /
-      DATA (AE(I, 5, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,33.9D0  ,19.8D0  ,20.3D0  ,20.7D0  ,21.9D0  ,
-     +  23.4D0  ,24.1D0  ,24.8D0   /
-      DATA (AE(I, 5,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,18.0D0  ,20.0D0  ,21.4D0  ,22.0D0  ,
-     +  23.8D0  ,24.6D0  ,25.4D0   /
-      DATA (AE(I, 5,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,26.4D0  ,20.4D0  ,21.2D0  ,22.3D0  ,
-     +  23.8D0  ,24.7D0  ,25.5D0   /
-      DATA (AE(I, 5,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,41.7D0  ,18.2D0  ,19.8D0  ,21.1D0  ,
-     +  22.6D0  ,23.4D0  ,24.6D0   /
-      DATA (AE(I, 5,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,22.5D0  ,20.0D0  ,21.7D0  ,
-     +  22.8D0  ,23.7D0  ,24.7D0   /
-      DATA (AE(I, 5,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,54.1D0  ,19.9D0  ,21.9D0  ,
-     +  23.2D0  ,24.3D0  ,25.3D0   /
-      DATA (AE(I, 5,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,21.2D0  ,22.2D0  ,
-     +  23.6D0  ,24.9D0  ,25.5D0   /
-      DATA (AE(I, 5,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,44.9D0  ,21.9D0  ,
-     +  23.8D0  ,25.2D0  ,25.6D0   /
-      DATA (AE(I, 5,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,47.8D0  ,22.7D0  ,
-     +  23.8D0  ,24.9D0  ,26.3D0   /
-      DATA (AE(I, 5,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,35.5D0  ,
-     +  23.9D0  ,25.9D0  ,26.6D0   /
-      DATA (AE(I, 5,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,64.3D0  ,
-     +  24.1D0  ,25.7D0  ,27.1D0   /
-      DATA (AE(I, 5,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  34.0D0  ,25.7D0  ,27.7D0   /
-      DATA (AE(I, 6, 1),I=1,10)  / 
-     +  16.6D0  ,16.5D0  ,16.8D0  ,16.7D0  ,17.0D0  ,16.5D0  ,16.7D0  ,
-     +  18.3D0  ,18.9D0  ,19.0D0   /
-      DATA (AE(I, 6, 2),I=1,10)  / 
-     +  16.2D0  ,16.6D0  ,17.2D0  ,17.4D0  ,17.9D0  ,17.4D0  ,17.7D0  ,
-     +  20.7D0  ,22.0D0  ,22.6D0   /
-      DATA (AE(I, 6, 3),I=1,10)  / 
-     +  18.9D0  ,18.7D0  ,18.8D0  ,18.6D0  ,18.9D0  ,18.6D0  ,18.9D0  ,
-     +  21.0D0  ,22.3D0  ,22.9D0   /
-      DATA (AE(I, 6, 4),I=1,10)  / 
-     +  18.3D0  ,12.7D0  ,14.2D0  ,15.0D0  ,15.7D0  ,16.1D0  ,16.3D0  ,
-     +  16.5D0  ,17.9D0  ,19.0D0   /
-      DATA (AE(I, 6, 5),I=1,10)  / 
-     +  .000D+00,15.7D0  ,15.1D0  ,15.3D0  ,16.5D0  ,16.4D0  ,16.4D0  ,
-     +  17.0D0  ,18.3D0  ,19.4D0   /
-      DATA (AE(I, 6, 6),I=1,10)  / 
-     +  .000D+00,22.9D0  ,14.9D0  ,15.2D0  ,16.2D0  ,16.9D0  ,17.4D0  ,
-     +  18.2D0  ,19.5D0  ,21.1D0   /
-      DATA (AE(I, 6, 7),I=1,10)  / 
-     +  .000D+00,40.7D0  ,18.4D0  ,15.9D0  ,17.1D0  ,17.7D0  ,18.9D0  ,
-     +  19.5D0  ,20.3D0  ,21.1D0   /
-      DATA (AE(I, 6, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,23.3D0  ,16.2D0  ,16.3D0  ,17.3D0  ,18.7D0  ,
-     +  19.5D0  ,20.3D0  ,21.1D0   /
-      DATA (AE(I, 6, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,49.2D0  ,19.0D0  ,19.1D0  ,19.4D0  ,20.2D0  ,
-     +  20.8D0  ,21.6D0  ,22.0D0   /
-      DATA (AE(I, 6,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,27.2D0  ,21.2D0  ,20.8D0  ,21.4D0  ,
-     +  22.3D0  ,22.8D0  ,23.3D0   /
-      DATA (AE(I, 6,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,45.6D0  ,25.0D0  ,22.8D0  ,23.9D0  ,
-     +  23.6D0  ,24.3D0  ,24.4D0   /
-      DATA (AE(I, 6,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,45.8D0  ,29.7D0  ,25.1D0  ,25.3D0  ,
-     +  25.3D0  ,26.0D0  ,26.3D0   /
-      DATA (AE(I, 6,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,42.7D0  ,29.0D0  ,28.0D0  ,
-     +  27.0D0  ,27.2D0  ,27.6D0   /
-      DATA (AE(I, 6,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,62.0D0  ,32.0D0  ,30.0D0  ,
-     +  29.8D0  ,29.5D0  ,29.6D0   /
-      DATA (AE(I, 6,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,44.5D0  ,34.4D0  ,
-     +  32.7D0  ,31.5D0  ,31.8D0   /
-      DATA (AE(I, 6,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,75.6D0  ,37.1D0  ,
-     +  34.6D0  ,34.4D0  ,34.4D0   /
-      DATA (AE(I, 6,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,51.2D0  ,45.2D0  ,
-     +  39.0D0  ,37.5D0  ,36.4D0   /
-      DATA (AE(I, 6,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,74.9D0  ,
-     +  42.3D0  ,39.9D0  ,38.3D0   /
-      DATA (AE(I, 6,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,69.5D0  ,
-     +  50.7D0  ,42.3D0  ,41.4D0   /
-      DATA (AE(I, 6,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  66.3D0  ,48.0D0  ,43.4D0   /
-      DATA (AE(I, 7, 1),I=1,10)  / 
-     +  27.0D0  ,25.8D0  ,26.3D0  ,26.2D0  ,26.7D0  ,26.7D0  ,27.1D0  ,
-     +  27.1D0  ,27.2D0  ,19.0D0   /
-      DATA (AE(I, 7, 2),I=1,10)  / 
-     +  29.1D0  ,28.9D0  ,29.7D0  ,30.3D0  ,31.0D0  ,31.4D0  ,32.0D0  ,
-     +  32.3D0  ,32.7D0  ,22.6D0   /
-      DATA (AE(I, 7, 3),I=1,10)  / 
-     +  31.6D0  ,29.7D0  ,30.9D0  ,31.4D0  ,32.5D0  ,33.1D0  ,34.0D0  ,
-     +  34.6D0  ,35.1D0  ,22.9D0   /
-      DATA (AE(I, 7, 4),I=1,10)  / 
-     +  27.4D0  ,19.9D0  ,20.8D0  ,22.8D0  ,24.6D0  ,26.4D0  ,28.2D0  ,
-     +  29.6D0  ,30.8D0  ,19.0D0   /
-      DATA (AE(I, 7, 5),I=1,10)  / 
-     +  .000D+00,24.6D0  ,24.1D0  ,25.0D0  ,27.2D0  ,28.7D0  ,30.7D0  ,
-     +  31.8D0  ,32.9D0  ,19.4D0   /
-      DATA (AE(I, 7, 6),I=1,10)  / 
-     +  .000D+00,35.6D0  ,25.2D0  ,25.6D0  ,27.9D0  ,30.4D0  ,32.7D0  ,
-     +  34.6D0  ,36.3D0  ,21.1D0   /
-      DATA (AE(I, 7, 7),I=1,10)  / 
-     +  .000D+00,45.4D0  ,30.9D0  ,28.2D0  ,29.0D0  ,31.2D0  ,34.0D0  ,
-     +  35.8D0  ,37.4D0  ,21.1D0   /
-      DATA (AE(I, 7, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,38.2D0  ,29.6D0  ,29.4D0  ,30.3D0  ,33.2D0  ,
-     +  35.5D0  ,37.6D0  ,21.1D0   /
-      DATA (AE(I, 7, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,59.3D0  ,34.5D0  ,33.7D0  ,32.9D0  ,35.4D0  ,
-     +  37.6D0  ,39.6D0  ,22.0D0   /
-      DATA (AE(I, 7,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,44.5D0  ,37.8D0  ,37.5D0  ,37.2D0  ,
-     +  39.0D0  ,41.4D0  ,23.3D0   /
-      DATA (AE(I, 7,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,67.0D0  ,43.6D0  ,42.0D0  ,40.8D0  ,
-     +  41.4D0  ,43.0D0  ,24.4D0   /
-      DATA (AE(I, 7,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,49.9D0  ,50.9D0  ,44.6D0  ,43.9D0  ,
-     +  44.2D0  ,44.2D0  ,26.3D0   /
-      DATA (AE(I, 7,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,67.2D0  ,50.5D0  ,48.7D0  ,
-     +  48.1D0  ,47.2D0  ,27.6D0   /
-      DATA (AE(I, 7,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,68.1D0  ,55.2D0  ,52.3D0  ,
-     +  51.5D0  ,51.6D0  ,29.6D0   /
-      DATA (AE(I, 7,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,68.7D0  ,58.6D0  ,
-     +  56.5D0  ,55.7D0  ,31.8D0   /
-      DATA (AE(I, 7,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,89.3D0  ,62.9D0  ,
-     +  60.0D0  ,59.1D0  ,34.4D0   /
-      DATA (AE(I, 7,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,56.0D0  ,72.9D0  ,
-     +  66.3D0  ,64.2D0  ,36.4D0   /
-      DATA (AE(I, 7,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,105.D0  ,
-     +  71.3D0  ,68.3D0  ,38.3D0   /
-      DATA (AE(I, 7,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,73.4D0  ,
-     +  76.8D0  ,72.4D0  ,41.4D0   /
-      DATA (AE(I, 7,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  107.D0  ,79.9D0  ,43.4D0   /
-      DATA (AE(I, 8, 1),I=1,10)  / 
-     +  35.5D0  ,35.3D0  ,35.7D0  ,35.7D0  ,36.3D0  ,36.3D0  ,36.7D0  ,
-     +  36.7D0  ,36.7D0  ,19.0D0   /
-      DATA (AE(I, 8, 2),I=1,10)  / 
-     +  40.6D0  ,41.4D0  ,41.9D0  ,42.3D0  ,43.2D0  ,43.5D0  ,44.0D0  ,
-     +  44.3D0  ,44.5D0  ,22.6D0   /
-      DATA (AE(I, 8, 3),I=1,10)  / 
-     +  45.4D0  ,45.7D0  ,46.4D0  ,47.0D0  ,48.1D0  ,48.7D0  ,49.4D0  ,
-     +  49.8D0  ,50.2D0  ,22.9D0   /
-      DATA (AE(I, 8, 4),I=1,10)  / 
-     +  43.9D0  ,44.3D0  ,43.4D0  ,45.1D0  ,47.3D0  ,48.7D0  ,49.6D0  ,
-     +  50.5D0  ,51.3D0  ,19.0D0   /
-      DATA (AE(I, 8, 5),I=1,10)  / 
-     +  .000D+00,49.3D0  ,49.6D0  ,50.5D0  ,53.2D0  ,54.2D0  ,55.4D0  ,
-     +  56.1D0  ,56.8D0  ,19.4D0   /
-      DATA (AE(I, 8, 6),I=1,10)  / 
-     +  .000D+00,59.1D0  ,53.0D0  ,55.4D0  ,58.0D0  ,60.0D0  ,61.2D0  ,
-     +  62.5D0  ,63.6D0  ,21.1D0   /
-      DATA (AE(I, 8, 7),I=1,10)  / 
-     +  .000D+00,54.5D0  ,57.1D0  ,59.2D0  ,62.3D0  ,64.4D0  ,66.0D0  ,
-     +  67.3D0  ,68.5D0  ,21.1D0   /
-      DATA (AE(I, 8, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,65.9D0  ,62.1D0  ,65.1D0  ,67.6D0  ,69.4D0  ,
-     +  71.1D0  ,72.6D0  ,21.1D0   /
-      DATA (AE(I, 8, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,72.2D0  ,67.1D0  ,70.5D0  ,73.1D0  ,75.1D0  ,
-     +  76.8D0  ,78.4D0  ,22.0D0   /
-      DATA (AE(I, 8,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,80.1D0  ,75.0D0  ,78.0D0  ,80.0D0  ,
-     +  82.1D0  ,83.9D0  ,23.3D0   /
-      DATA (AE(I, 8,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,94.5D0  ,82.2D0  ,82.8D0  ,85.1D0  ,
-     +  87.3D0  ,89.2D0  ,24.4D0   /
-      DATA (AE(I, 8,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,56.8D0  ,92.5D0  ,87.2D0  ,89.4D0  ,
-     +  91.9D0  ,94.1D0  ,26.3D0   /
-      DATA (AE(I, 8,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,116.D0  ,96.2D0  ,94.4D0  ,
-     +  97.0D0  ,99.2D0  ,27.6D0   /
-      DATA (AE(I, 8,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,78.1D0  ,104.D0  ,102.D0  ,
-     +  102.D0  ,105.D0  ,29.6D0   /
-      DATA (AE(I, 8,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,128.D0  ,111.D0  ,
-     +  109.D0  ,110.D0  ,31.8D0   /
-      DATA (AE(I, 8,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,104.D0  ,118.D0  ,
-     +  117.D0  ,115.D0  ,34.4D0   /
-      DATA (AE(I, 8,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,64.4D0  ,138.D0  ,
-     +  124.D0  ,122.D0  ,36.4D0   /
-      DATA (AE(I, 8,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,133.D0  ,
-     +  133.D0  ,132.D0  ,38.3D0   /
-      DATA (AE(I, 8,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,83.6D0  ,
-     +  146.D0  ,139.D0  ,41.4D0   /
-      DATA (AE(I, 8,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  166.D0  ,147.D0  ,43.4D0   /
-      DATA (AE(I, 9, 1),I=1,10)  / 
-     +  43.3D0  ,43.2D0  ,43.6D0  ,43.8D0  ,44.1D0  ,44.3D0  ,44.7D0  ,
-     +  44.8D0  ,44.8D0  ,19.0D0   /
-      DATA (AE(I, 9, 2),I=1,10)  / 
-     +  50.9D0  ,51.4D0  ,52.0D0  ,52.6D0  ,53.1D0  ,53.6D0  ,54.2D0  ,
-     +  54.5D0  ,54.7D0  ,22.6D0   /
-      DATA (AE(I, 9, 3),I=1,10)  / 
-     +  58.0D0  ,58.4D0  ,59.3D0  ,60.1D0  ,60.7D0  ,61.5D0  ,62.3D0  ,
-     +  62.7D0  ,63.1D0  ,22.9D0   /
-      DATA (AE(I, 9, 4),I=1,10)  / 
-     +  62.0D0  ,63.9D0  ,63.7D0  ,65.7D0  ,65.5D0  ,67.5D0  ,68.2D0  ,
-     +  68.9D0  ,69.7D0  ,19.0D0   /
-      DATA (AE(I, 9, 5),I=1,10)  / 
-     +  .000D+00,72.2D0  ,72.5D0  ,74.2D0  ,74.2D0  ,76.1D0  ,77.0D0  ,
-     +  77.8D0  ,78.6D0  ,19.4D0   /
-      DATA (AE(I, 9, 6),I=1,10)  / 
-     +  .000D+00,80.4D0  ,80.5D0  ,83.1D0  ,83.0D0  ,85.5D0  ,86.8D0  ,
-     +  88.1D0  ,89.2D0  ,21.1D0   /
-      DATA (AE(I, 9, 7),I=1,10)  / 
-     +  .000D+00,63.4D0  ,88.5D0  ,91.3D0  ,91.1D0  ,94.0D0  ,95.8D0  ,
-     +  97.3D0  ,98.6D0  ,21.1D0   /
-      DATA (AE(I, 9, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,98.8D0  ,98.6D0  ,97.8D0  ,102.D0  ,104.D0  ,
-     +  106.D0  ,108.D0  ,21.1D0   /
-      DATA (AE(I, 9, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,84.1D0  ,107.D0  ,107.D0  ,111.D0  ,113.D0  ,
-     +  116.D0  ,117.D0  ,22.0D0   /
-      DATA (AE(I, 9,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,116.D0  ,115.D0  ,119.D0  ,122.D0  ,
-     +  125.D0  ,127.D0  ,23.3D0   /
-      DATA (AE(I, 9,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,111.D0  ,123.D0  ,127.D0  ,131.D0  ,
-     +  134.D0  ,137.D0  ,24.4D0   /
-      DATA (AE(I, 9,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,65.6D0  ,136.D0  ,135.D0  ,140.D0  ,
-     +  143.D0  ,146.D0  ,26.3D0   /
-      DATA (AE(I, 9,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,146.D0  ,144.D0  ,149.D0  ,
-     +  152.D0  ,155.D0  ,27.6D0   /
-      DATA (AE(I, 9,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,88.7D0  ,152.D0  ,158.D0  ,
-     +  162.D0  ,165.D0  ,29.6D0   /
-      DATA (AE(I, 9,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,181.D0  ,167.D0  ,
-     +  171.D0  ,174.D0  ,31.8D0   /
-      DATA (AE(I, 9,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,117.D0  ,174.D0  ,
-     +  180.D0  ,183.D0  ,34.4D0   /
-      DATA (AE(I, 9,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,72.0D0  ,201.D0  ,
-     +  189.D0  ,192.D0  ,36.4D0   /
-      DATA (AE(I, 9,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,151.D0  ,
-     +  198.D0  ,201.D0  ,38.3D0   /
-      DATA (AE(I, 9,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,95.2D0  ,
-     +  220.D0  ,210.D0  ,41.4D0   /
-      DATA (AE(I, 9,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  192.D0  ,217.D0  ,43.4D0   /
-      DATA (AE(I,10, 1),I=1,10)  / 
-     +  62.1D0  ,62.1D0  ,62.6D0  ,62.9D0  ,63.3D0  ,63.3D0  ,64.0D0  ,
-     +  64.0D0  ,64.0D0  ,19.0D0   /
-      DATA (AE(I,10, 2),I=1,10)  / 
-     +  75.1D0  ,75.4D0  ,76.3D0  ,76.8D0  ,77.6D0  ,77.9D0  ,78.8D0  ,
-     +  79.0D0  ,79.3D0  ,22.6D0   /
-      DATA (AE(I,10, 3),I=1,10)  / 
-     +  87.5D0  ,88.3D0  ,89.4D0  ,90.2D0  ,91.3D0  ,91.9D0  ,93.0D0  ,
-     +  93.5D0  ,93.9D0  ,22.9D0   /
-      DATA (AE(I,10, 4),I=1,10)  / 
-     +  104.D0  ,104.D0  ,105.D0  ,106.D0  ,107.D0  ,108.D0  ,109.D0  ,
-     +  110.D0  ,110.D0  ,19.0D0   /
-      DATA (AE(I,10, 5),I=1,10)  / 
-     +  .000D+00,122.D0  ,122.D0  ,123.D0  ,124.D0  ,125.D0  ,126.D0  ,
-     +  127.D0  ,128.D0  ,19.4D0   /
-      DATA (AE(I,10, 6),I=1,10)  / 
-     +  .000D+00,138.D0  ,139.D0  ,140.D0  ,142.D0  ,143.D0  ,144.D0  ,
-     +  146.D0  ,147.D0  ,21.1D0   /
-      DATA (AE(I,10, 7),I=1,10)  / 
-     +  .000D+00,85.3D0  ,158.D0  ,159.D0  ,161.D0  ,162.D0  ,164.D0  ,
-     +  166.D0  ,167.D0  ,21.1D0   /
-      DATA (AE(I,10, 8),I=1,10)  / 
-     +  .000D+00,.000D+00,176.D0  ,177.D0  ,179.D0  ,181.D0  ,183.D0  ,
-     +  184.D0  ,186.D0  ,21.1D0   /
-      DATA (AE(I,10, 9),I=1,10)  / 
-     +  .000D+00,.000D+00,114.D0  ,199.D0  ,201.D0  ,202.D0  ,205.D0  ,
-     +  206.D0  ,207.D0  ,22.0D0   /
-      DATA (AE(I,10,10),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,218.D0  ,219.D0  ,220.D0  ,224.D0  ,
-     +  225.D0  ,226.D0  ,23.3D0   /
-      DATA (AE(I,10,11),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,150.D0  ,238.D0  ,238.D0  ,243.D0  ,
-     +  244.D0  ,245.D0  ,24.4D0   /
-      DATA (AE(I,10,12),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,85.8D0  ,255.D0  ,255.D0  ,261.D0  ,
-     +  262.D0  ,263.D0  ,26.3D0   /
-      DATA (AE(I,10,13),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,195.D0  ,272.D0  ,279.D0  ,
-     +  279.D0  ,280.D0  ,27.6D0   /
-      DATA (AE(I,10,14),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,115.D0  ,290.D0  ,296.D0  ,
-     +  297.D0  ,298.D0  ,29.6D0   /
-      DATA (AE(I,10,15),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,263.D0  ,313.D0  ,
-     +  314.D0  ,315.D0  ,31.8D0   /
-      DATA (AE(I,10,16),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,150.D0  ,330.D0  ,
-     +  331.D0  ,332.D0  ,34.4D0   /
-      DATA (AE(I,10,17),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,90.0D0  ,319.D0  ,
-     +  349.D0  ,349.D0  ,36.4D0   /
-      DATA (AE(I,10,18),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,196.D0  ,
-     +  366.D0  ,367.D0  ,38.3D0   /
-      DATA (AE(I,10,19),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,122.D0  ,
-     +  387.D0  ,384.D0  ,41.4D0   /
-      DATA (AE(I,10,20),I=1,10)  / 
-     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
-     +  247.D0  ,401.D0  ,43.4D0     /
-      DATA (ERES(I, 1),I=1,10)  / 10*0.D0/
-      DATA (ERES(I, 2),I=1,10)  / 10*0.D0/
-      DATA (ERES(I, 3),I=1,10)  / 10*0.D0/
-      DATA (ERES(I, 4),I=1,10)  / 10*0.D0/
-      DATA (ERES(I, 5),I=1,10)  / 10*0.D0/
-      DATA (ERES(I, 6),I=1,10)  / 
-     +    0.000D0, 0.000D0, 0.000D0, 0.000D0, 0.000D0, 0.000D0, 0.000D0,
-     +    2.780D0, 2.880D0, 2.890D0 /
-      DATA (ERES(I, 7),I=1,10)  / 
-     +    1.500D0, 2.460D0, 2.510D0, 2.610D0, 2.700D0, 2.920D0, 3.070D0,
-     +    3.200D0, 3.330D0, 2.890D0 /
-      DATA (ERES(I, 8),I=1,10)  / 
-     +    4.470D0, 4.350D0, 4.390D0, 4.550D0, 4.660D0, 4.890D0, 4.980D0,
-     +    5.100D0, 5.220D0, 2.890D0 /
-      DATA (ERES(I, 9),I=1,10)  / 
-     +    7.480D0, 7.380D0, 7.370D0, 7.480D0, 7.510D0, 7.630D0, 7.660D0,
-     +    7.750D0, 7.820D0, 2.890D0 /
-      DATA (ERES(I,10),I=1,10)  / 
-     +   15.270D0,15.190D0,15.200D0,15.370D0,15.380D0,15.430D0,15.540D0,
-     +   15.590D0,15.630D0, 2.890D0 /
-      END
-C->
-C=======================================================================
-
-      SUBROUTINE FRAGM (IAT,IAP, NW,B, NF, IAF)
-
-C-----------------------------------------------------------------------
-C...Nuclear Fragmentation, Abrasion-ablation model, 
-C...Based on Jon Engel's routines ABRABL 
-C...This most recent version adds for all prefragment
-C...masses > 10 the model calculation for the fragment
-C...mass distribution and the energy carried by the fragment
-C...of W. Friedmann
-C...The average values are used to implement the model
-C...in the montecarlo fashion / TSS, Dec '91
-C.
-C.  INPUT: IAP = mass of incident nucleus
-C.         IAT = mass of target   nucleus
-C.         NW = number of wounded nucleons in the beam nucleus
-C.         B  = impact parameter in the interaction
-C.     
-C.  OUTPUT : NF = number of fragments  of the spectator nucleus
-C.           IAF(1:NF) = mass number of each fragment
-C.           PF(3,60) in common block /FRAGMENTS/ contains
-C.           the three momentum components (MeV/c) of each
-C.           fragment in the projectile frame
-C..............................................................
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-
-      COMMON /FRAGMENTS/ PPP(3,60)
-      COMMON /FRAGMOD/A(10,10,20),AE(10,10,20),ERES(10,10),NFLAGG(10,10)
-      DIMENSION IAF(60)
-      DIMENSION AA(10), EAA(10) 
-      SAVE
-      EXTERNAL GASDEV
-      DATA AA/10.D0,15.D0,20.D0,25.D0,30.D0,35.D0,40.D0,45.D0,50.D0,
-     $        56.D0/
-      DATA EAA/1.D0,2.D0,4.D0,6.D0,8.D0,10.D0,12.D0,16.D0,20.D0,30.D0/
-
-      AP=IAP
-      AT=IAT
-      NPF = IAP - NW
-      IF (NPF .EQ. 0) THEN
-         NF = 0
-         RETURN
-      ENDIF
-
-      EB = ESTAR(AP,AT, B)
-      EBP = ESTARP (NPF, NW)
-C CONTRIBUTION TO E* FROM ENERGY DEPOSITED BY SECONDARIES
-      EB = EB + EBP
-C TOTAL E* IS THE SUM OF THE TWO COMPONENTS
-
-C.....Prefragment transverse momentum (MeV/nucleon)...
-            FK = FERMK(AP)
-C FERMI MOMENTUM OF THE PROJECTILE NUCLEUS
-            IF (NW .LT. IAP) THEN
-            SIG = FK*DSQRT(NW*NPF/(AP-1.D0))/3.162D0
-C GAUSSIAN SIGMA IN ALL THREE DIRECTION
-            ELSE
-            SIG = FK/3.162D0
-C THIS IS NOT CORRECT, TOO LARGE !!!!!!!!!!!!!!
-            ENDIF
-             PPFX = SIG*GASDEV(0)/NPF
-             PPFY = SIG*GASDEV(1)/NPF
-C THREE MOMENTUM COMPONENTS PER NUCLEON FOR THE PREFRAGMENT
-
-C.............Crude model for small prefragment mass .......
-            IF (NPF .LT. 10) THEN
-                 CALL EVAP(NPF, EB, EPS, NNUC, NALP)
-C   EPS IS THE KINETIC ENERGY CARRIED BY THE EVAPORATED NUCLEONS
-               ETOT = 938.D0 + EPS
-                 PP = SQRT((ETOT*ETOT - 8.79844D5)/3.D0)
-C   AVERAGE MOMENTUM OF EVAPORATED NUCLEONS IN EACH DIRECTION
-                 NUC = NPF - NNUC - 4*NALP
-                 NF = 0
-                 IF (NUC .GT. 0) THEN
-                    NF = NF + 1
-                    IAF(NF) = NUC
-                    PPP(1,NF) = NUC*PPFX
-                    PPP(2,NF) = NUC*PPFY
-                 ENDIF
-                 IF (NALP .NE. 0) THEN
-                 DO I=1,NALP
-                   NF = NF + 1
-                    IAF(NF) = 4
-                   CALL SINCO(S1,C1)
-                   CALL SINCO(S2,C2)
-                   PXE = 4.D0*PP*S1*S2
-                   PYE = 4.D0*PP*S1*C2
-                   PPP(1,NF) = 4.D0*PPFX + PXE
-                   PPP(2,NF) = 4.D0*PPFY + PYE
-                   PPP(1,1) = PPP(1,1) - PXE
-                   PPP(2,1) = PPP(2,1) - PYE
-                 ENDDO
-                 ENDIF
-                 IF (NNUC .NE. 0) THEN
-                 DO I=1,NNUC
-                    NF = NF + 1
-                    IAF(NF) = 1
-                    CALL SINCO(S1,C1)
-                    CALL SINCO(S2,C2)
-                    PXE = PP*S1*S2
-                    PYE = PP*S1*C2
-                    PPP(1,NF) = 4.D0*PPFX + PXE
-                    PPP(2,NF) = 4.D0*PPFY + PYE
-                    PPP(1,1) = PPP(1,1) - PXE
-                    PPP(2,1) = PPP(2,1) - PYE
-                 ENDDO
-                 ENDIF
-                 RETURN
-            ENDIF
-
-C.........More refined model calculation .............
-      JA = NPF/5 -1
-      IF (JA .LT. 10) THEN
-      IF ((NPF - AA(JA)) .GT. (AA(JA+1)-NPF)) JA = JA + 1
-      ENDIF
-      ARAT = DBLE(NPF)/AA(JA)
-      DO J=1,10
-      IF (EB .LT. EAA(J)) GO TO 29
-      ENDDO
-      JE = 10
-      GO TO 39
-   29      JE = J
-   39      IF (JE .GT. 1 .AND. JE .NE. 10) THEN
-      IF ((EB - EAA(J-1)) .LT. (EAA(J)-EB)) JE = J - 1
-      ENDIF
-      ERAT = EB/EAA(JE)
-        IF (EB .LT. 1.D0) THEN
-        ERAT = EB
-        ENDIF
-C INTERPOLATE BETWEEN EB=0. (NOTHING HAPPENS) AND EB = 1. MeV
-         IF (JA .EQ. 10 .AND. JE .GT. 6) THEN
-            WRITE(*,*)' JA=',JA,',   JE=',JE
-         ENDIF
-   43    ESUM = 0.D0
-      NSUM = 0
-      JF = 0
-      DO J=20,1,-1
-        FR =  A(JA, JE, J)*ARAT*ERAT
-        N1 = INT(1.D0 + FR)
-        FR1 = FR/DBLE(N1)
-        DO K=1, N1
-          IF (S_RNDM(0) .LT. FR1) THEN
-            JF = JF + 1
-            IAF(JF) = J
-            NSUM = NSUM + J
-            EKIN = ERAT*AE(JA,JE, J)
-            IF (EKIN .GT. 0.D0) THEN
-              ESUM = ESUM + EKIN
-              ETOT = 938.D0*IAF(JF) + EKIN
-              PP = DSQRT(2.D0*(ETOT*ETOT - IAF(JF)**2*8.79844D5)/3.D0)
-              CALL SINCO(S1,C1)
-              CALL SINCO(S2,C2)
-              PPP(1,JF) = PP*S1*S2 + IAF(JF)*PPFX
-              PPP(2,JF) = PP*S1*C2 + IAF(JF)*PPFY
-            ENDIF
-            IF (NSUM .GT. NPF) THEN
-C           WRITE(*,*)' WARNING, NSUM=', NSUM,',  NPF=',NPF
-C           WRITE(*,*)'  ARAT =', ARAT
-              GO TO 43
-            ELSE
-              IF (NSUM .EQ. NPF) THEN
-                GO TO 44
-              ENDIF
-            ENDIF
-          ENDIF
-        ENDDO
-      ENDDO
-      IF (NFLAGG(JA,JE) .EQ. 0) THEN
-C 'THE RESIDUE' IS A NUCLEAR FRAGMENT
-        JF = JF + 1
-        IAF(JF) = NPF - NSUM
-        F1 = NPF*EB - ESUM
-        IF (F1 .LT. 0.D0) F1 = 0.D0
-C GIVE THE REST OF EB TO THE FRAGMENT
-        EKIN = F1
-        IF (EKIN .GT. 0.D0) THEN
-          ETOT = 938.D0*IAF(JF) + EKIN
-          PP = DSQRT(2.D0*(ETOT*ETOT - IAF(JF)**2*8.79844D5)/3.D0)
-          CALL SINCO(S1,C1)
-          CALL SINCO(S2,C2)
-          PPP(1,JF) = PP*S1*S2 + IAF(JF)*PPFX
-          PPP(2,JF) = PP*S1*C2 + IAF(JF)*PPFY
-        ENDIF
-      ELSE
-C 'THE RESIDUE' CONSISTS OF SPECTATOR NUCLEONS
-        N1 = NPF - NSUM
-        DO K=1,N1
-          JF = JF + 1
-          IAF(JF) = 1
-          EKIN = ERAT*ERES(JA,JE)
-          IF (EKIN .GT. 0.D0) THEN
-            ETOT = 938.D0*IAF(JF) + EKIN
-            PP = DSQRT(2.D0*(ETOT*ETOT - IAF(JF)**2*8.79844D5)/3.D0)
-            CALL SINCO(S1,C1)
-            CALL SINCO(S2,C2)
-            PPP(1,JF) = PP*S1*S2 + PPFX
-            PPP(2,JF) = PP*S1*C2 + PPFY
-          ENDIF
-        ENDDO
-      ENDIF
-  44  NF = JF
-      RETURN
-      END
-C->
-C=======================================================================
-
-      FUNCTION ESTARP (NPF, NW)
-
-C-----------------------------------------------------------------------
-C CONTRIBUTION TO E* FROM ENERGY DEPOSITED BY SECONDARIES
-C VERY NAIVE VERSION INCORPORATING HUEFFNER'S IDEAS
-C-----------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      SAVE
-
-      APF = NPF
-      F1 = 15.3D0/APF**0.666666666D0
-C AVERAGE KINETIC ENERGY/NUCLEON IN PREFRAGMENT (MeV)
-C PER PATHLENGTH EQUAL TO THE PREFRAGMENT RADIUS
-      ESTARP = 0.D0
-      DO I=1,NW
-        IF (S_RNDM(0) .GT. 0.5D0) THEN
-          F2 = F1*RDIS(0)
-          ESTARP = ESTARP + F2
-        ENDIF
-      ENDDO
-C SAMPLE RANDOMLY PER WOUNDED NUCLEON, x NW
-      RETURN
-      END
-C=======================================================================
-      
-      FUNCTION RDIS(Idum)
-
-C-----------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-
-      dimension probr(20)
-      SAVE
-      data probr/
-     *      0.10000D0, 0.15748D0, 0.21778D0, 0.28605D0, 0.36060D0,
-     *      0.43815D0, 0.51892D0, 0.60631D0, 0.70002D0, 0.79325D0,
-     *      0.88863D0, 0.98686D0, 1.10129D0, 1.21202D0, 1.32932D0,
-     *      1.44890D0, 1.57048D0, 1.70139D0, 1.83417D0, 2.00000D0/
-
-      rdis = idum
-      nr = INT(20.D0*S_RNDM(0) + 1.D0)
-      if (nr .eq. 1) then
-        f1 = 0.D0
-      else
-        f1 = probr(nr-1)
-      endif
-      dr = probr(nr) - f1
-      rdis = f1 + dr*S_RNDM(1)
-      return
-      end
-
-C=======================================================================
-
-      FUNCTION ESTAR(ap,at,b)
-
-C-----------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-
-      SAVE
-
-c      real*4 ap,at,b,estar
-      sigma=4.5D0  !total n-n cross section in fm**2
-      rt=.82d0*at**.33333333D0 !target radius
-      rp=.82d0*ap**.33333333D0 !projectile radius
-      alpha=rt**2/rp**2
-      beta=b**2/rt**2
-      f=at*sigma/(PI*rt**2)
-      alf = log(f)
-      alalf = log(alpha)
-      gfac=0.d0
-      gfac1=0.d0
-      s1=0.D0
-      s2=0.D0
-      s3=0.D0
-      ii=1
-      do n=0,10 ! This limit may not need to be so high.
-         if(n.ge.2) then
-            gfac1=gfac
-            gfac=gfac+log(float(n)) 
-         endif
-         g0=n*alf -n*beta*alpha/(n+alpha)+alalf
-         g1=g0-log(alpha+n)-gfac
-         g2=(n+2)*log(f)-(n+2)*beta*alpha/(n+2+alpha) 
-     >      +log(n+2+alpha+beta*alpha**2)-3.d0*log(n+2.d0+alpha)-gfac
-         g3=g0-2.d0*log(n+alpha)-gfac1
-         ii=-ii
-         s1=s1+ii*exp(g1)
-         s2=s2+ii*exp(g2)
-         if(n.ge.1) s3=s3+ii*exp(g3)
-      enddo
-
-      pb=s1
-      e1b=197.D0**2/(2.D0*938.d0*rp**2*pb) *s2
-c      a=b*(s3/pb-1)
-c      a=-b*s3/pb
-c      e2b=-.5* 938. * (41./(ap**.333))**2 * a**2 /(197.**2)
-c      estar=e1b+e2b
-      estar = e1b
-      return
-      end
-C=======================================================================
-
-      SUBROUTINE EVAP(npf,eb,eps,nnuc,nalp)
-
-C-----------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      SAVE
-
-      eps=7.5D0+sqrt(8.D0*eb)
-      n=min(npf*int(eb/eps),npf)
-      nalp=n/5
-      nnuc=n-4*nalp
-      return
-      end
-C->
-C=======================================================================
-
-      FUNCTION FERMK(A)
-
-C-----------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      DIMENSION AA(6), FK(6)
-      SAVE
-      DATA AA/4.D0, 6.D0, 12.D0, 24.D0, 40.D0, 57.D0/
-      DATA FK/130.D0,169.D0,221.D0,235.D0,251.D0,260.D0/
-
-      DO I=2,4
-      IF (A .LT. AA(I)) GO TO 25
-      ENDDO
-      I = 5
-   25      F11 = AA(I-1)
-      F12 = AA(I)
-      F13 = AA(I+1)
-      F21 = FK(I-1)
-      F22 = FK(I)
-      F23 = FK(I+1)
-      FERMK = QUAD_INT(A,F11,F12,F13, F21,F22,F23)
-      RETURN
-      END
-
-C*=======================================================================
-C. Multiple interaction structure
-C========================================================================
-
-      SUBROUTINE INT_NUC (IA, IB, SIG0, SIGEL) 
-
-C-----------------------------------------------------------------------
-C...Compute with a montecarlo code  the  "multiple interaction structure"
-C.  of a nucleus-nucleus interaction
-C.
-C.  INPUT : IA            = mass of target nucleus
-C.          IB            = mass of projectile nucleus
-C.          SIG0 (mbarn)  = inelastic pp cross section
-C.          SIGEL(mbarn)  = elastic pp cross section
-C.
-C.  OUTPUT : in common block /CNUCMS/
-C.           B = impact parameter (fm)
-C.           BMAX = maximum impact parameter for generation
-C.           NTRY = number of "trials" before one interaction
-C.           NA = number of wounded nucleons in A
-C.           NB =    "        "        "     in B
-C.           NI = number of nucleon-nucleon inelastic interactions 
-C.           NAEL = number of elastically scattered nucleons in  A 
-C.           NBEL =    "         "           "          "    in  B
-C.           JJA(J)  [J=1:IA]   = number of inelastic interactions 
-C.                                of J-th nucleon of nucleus A
-C.           JJB(J)  [J=1:IB]   = number of inelastic interactions 
-C.                                of J-th nucleon of nucleus B
-C.           JJAEL(J)  [J=1:IA]   = number of elastic interactions 
-C.                                of J-th nucleon of nucleus A
-C.           JJBEL(J)  [J=1:IB]   = number of elastic interactions 
-C.                                of J-th nucleon of nucleus B
-C.           JJINT(J,K)  [J=1:NB, K=1:NA]  (0 = no interaction) 
-C.                                         (1 = interaction )
-C.                                         between nucleon J of A and K of B
-C-----------------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-
-      PARAMETER (IAMAX=56)
-      COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
-     +         ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
-     +         ,JJAEL(IAMAX), JJBEL(IAMAX)
-      DIMENSION XA(IAMAX), YA(IAMAX), XB(IAMAX), YB(IAMAX)
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-      SAVE
-
-      SIGT = SIG0 + SIGEL
-      R2  = 0.1D0 * SIG0/PI
-      R2T = 0.1D0 * SIGT/PI
-      BMAX = 15.D0                             ! fm
-      NTRY = 0
-      CALL NUC_CONF (IA, XA, YA)
-      CALL NUC_CONF (IB, XB, YB)
-      NI = 0
-      NIEL = 0
-      DO JA=1,IA
-         JJA(JA) = 0
-         JJAEL(JA) = 0
-      ENDDO
-      DO JB=1,IB
-         JJB(JB) = 0
-         JJBEL(JB) = 0
-         DO JA=1,IA
-            JJINT(JB,JA) = 0
-         ENDDO
-      ENDDO
-1000  B = BMAX*SQRT(S_RNDM(0))
-      PHI = TWOPI*S_RNDM(1)
-      BX = B*COS(PHI)
-      BY = B*SIN(PHI)
-      NTRY = NTRY+1
-      DO JA=1,IA
-         DO JB=1,IB
-            S = (XA(JA)-XB(JB)-BX)**2 + (YA(JA)-YB(JB)-BY)**2
-            IF (S .LT. R2)  THEN
-               NI = NI + 1
-               JJA(JA) = JJA(JA)+1
-               JJB(JB) = JJB(JB)+1
-               JJINT(JB,JA) = 1
-            ELSE IF (S .LT. R2T)  THEN
-               NIEL = NIEL + 1
-               JJAEL(JA) = JJAEL(JA)+1
-               JJBEL(JB) = JJBEL(JB)+1
-            ENDIF
-         ENDDO
-      ENDDO
-      IF (NI + NIEL .EQ. 0)  GOTO 1000
-      NA = 0
-      NB = 0
-      NAEL = 0
-      NBEL = 0
-      DO JA=1,IA
-         IF (JJA(JA) .GT. 0)  THEN
-            NA = NA + 1
-         ELSE
-            IF (JJAEL(JA) .GT. 0)  NAEL = NAEL+1
-         ENDIF
-      ENDDO
-      DO JB=1,IB
-         IF (JJB(JB) .GT. 0)  THEN
-            NB = NB + 1
-         ELSE
-            IF (JJBEL(JB) .GT. 0)  NBEL = NBEL+1
-         ENDIF
-      ENDDO
-      RETURN
-      END
-C=======================================================================
-
-       SUBROUTINE NUC_CONF (IA, XX, YY)
-
-C-----------------------------------------------------------------------
-C...This routine generates the configuration  of a nucleus 
-C.  need an initialization call to NUC_GEOM_INI
-C.
-C.  INPUT  : IA = mass number of the nucleus
-C.  OUTPUT : XX(1:IA), YY(1:IA) (fm) = position in impact parameter
-C.                                     space of the IA nucleons
-C...................................................................
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-
-      PARAMETER (IAMAX=56)
-      DIMENSION XX(IAMAX), YY(IAMAX)
-      PARAMETER (NB=401)
-      COMMON /CPROFA/ ZMIN, DZ, BBZ(NB,IAMAX)
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-      SAVE
-
-      DO J=1,IA
-         Z = S_RNDM(J)
-         JZ = INT((Z-ZMIN)/DZ)+1
-         JZ = MIN(JZ,400)
-         T = (Z-ZMIN)/DZ - DBLE(JZ-1)
-         B = BBZ(JZ,IA)*(1.D0-T) + BBZ(JZ+1,IA)*T
-         PHI = TWOPI*S_RNDM(J+1)
-         XX(J) = B*COS(PHI)
-         YY(J) = B*SIN(PHI)
-      ENDDO
-      RETURN
-      END
-C=======================================================================
-
-      SUBROUTINE NUC_GEOM_INI
-
-C-----------------------------------------------------------------------
-C...Initialize all nucleus profiles
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-
-      PARAMETER (NB=401)
-      PARAMETER (IAMAX=56)
-      COMMON /CPROF/ DB, BMAX, BB(NB), TB(NB), A
-      COMMON /CPROFA/ ZMIN, DZ, BBZ(NB,IAMAX)
-      DIMENSION FFB(NB), GGB(NB)
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-      SAVE
-
-      CALL SHELL_INI
-      CALL WOOD_SAXON_INI
-      DO IA= 2,IAMAX
-           JA = IA
-         CALL NUC_PROFIL(JA)
-         DO K=1,NB
-           FFB(K) = BB(K)*TB(K) * TWOPI
-         ENDDO            
-         GGB(1) = 0.D0
-         GGB(NB) = 1.D0
-         DO K=2,NB-1
-           GGB(K) = GGB(K-1) + FFB(K-1)*DB
-         ENDDO            
-         CALL INVERT_ARRAY(GGB,0.D0,DB,NB, BBZ(1,IA), ZMIN, DZ)
-      ENDDO
-      RETURN
-      END
-C=======================================================================
-
-      SUBROUTINE NUC_PROFIL (JA)
-
-C-----------------------------------------------------------------------
-C...Compute the profile function T(b)
-C.  normalised as INT[d2b T(b) = 1]
-C.  INPUT : JA = integer mass number of nucleus
-C...............................................
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      PARAMETER (NB=401)
-      EXTERNAL DENSA
-      DOUBLE PRECISION DENSA
-      COMMON /CC01/  B
-      COMMON /CCDA/ JJA
-      COMMON /CPROF/ DB, BMAX, BB(NB), TB(NB), A
-      SAVE
-
-      BMAX = 7.5D0
-      DB = BMAX/DBLE(NB-1)
-      JJA = JA
-      A = JA
-      DO JB=1,NB
-        B = DB*DBLE(JB-1)
-        BB(JB) = B
-        IF (JA .LE. 18)  THEN
-            TB(JB) = PROFNUC (B, JA)
-         ELSE
-            TB(JB) = 2.D0*GAUSS (DENSA,0.D0,BMAX)
-         ENDIF
-      ENDDO
-      RETURN
-      END
-C=======================================================================
-
-      SUBROUTINE NUC1_PROFIL (AA)
-
-C-----------------------------------------------------------------------
-C...Compute the profile function T(b)
-C.  normalised as INT[d2b T(b) = 1]
-C.  INPUT : AA = mass number of nucleus
-C...............................................
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-
-      PARAMETER (NB=401)
-      EXTERNAL DENSA
-      DOUBLE PRECISION DENSA
-      COMMON /CC01/  B
-      COMMON /CPROF/ DB, BMAX, BB(NB), TB(NB), A
-      SAVE
-
-      A = AA
-      IA1 = INT(AA)
-      IA2 = IA1 + 1
-      U = AA - DBLE(IA1)
-      BMAX = 7.5D0
-      DB = BMAX/DBLE(NB-1)
-      DO JB=1,NB
-         B = DB*DBLE(JB-1)
-         BB(JB) = B
-         IF (A .LE. 18.D0)  THEN
-             T1 = PROFNUC (B, IA1)
-             T2 = PROFNUC (B, IA2)
-          ELSE
-             JJA = IA1
-             T1 = 2.D0*GAUSS (DENSA,0.D0,BMAX)
-             JJA = IA2
-             T2 = 2.D0*GAUSS (DENSA,0.D0,BMAX)
-          ENDIF
-          TB(JB) = (1.D0-U)*T1  + U*T2
-      ENDDO
-      RETURN
-      END
-
-C*======================================================================
-C.   Code about nuclear densities
-C=======================================================================
-
-      FUNCTION DENS_NUC (R, JA)
-
-C-----------------------------------------------------------------------
-C....Nuclear density (normalised to 1)
-C.   for a nucleus of mass number JA
-C.   INPUT R = radial coordinate  (fm)
-C.         JA = integer mass number
-C.  OUTPUT (fm**-3)
-C--------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      COMMON /CWOOD/ RR0(19:56), AA0(19:56), CC0(19:56)
-      SAVE
-
-      IF (JA .GT. 18)  THEN
-         DENS_NUC = WOOD_SAXON(R,JA)
-      ELSE IF (JA .NE. 4)  THEN
-         DENS_NUC = HELIUM(R)
-      ELSE
-         DENS_NUC = SHELL(R,JA)
-      ENDIF
-      RETURN
-      END
-C=======================================================================
-
-      FUNCTION WOOD_SAXON (R, JA) 
-
-C-----------------------------------------------------------------------
-C....Wood-Saxon nuclear density (normalised to 1)
-C.   for a nucleus of mass number A.
-C.   INPUT R =  (fm)
-C.         JA = mass number
-C.   OUTPUT (fm**-3)
-C------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      COMMON /CWOOD/ RR0(19:56), AA0(19:56), CC0(19:56)
-      SAVE
-
-      WOOD_SAXON = CC0(JA)/(1.D0+EXP((R-RR0(JA))/AA0(JA)))
-      RETURN
-      END      
-C=======================================================================
-
-      FUNCTION HELIUM (R)
-
-C-----------------------------------------------------------------------
-C... Helium density from Barrett and Jackson
-C.   INPUT R = r coordinate (fm)
-C.   OUTPUT (fm**-3)
-C........................................................
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      SAVE
-      DATA R0 /0.964D0/, CA /0.322D0/   ! fm
-      DATA W /0.517D0/, CC /5.993224D-02/
-
-      HELIUM = CC*(1.D0+W*(R/R0)**2)/(1.D0 + EXP((R-R0)/CA))
-      RETURN
-      END
-C=======================================================================
-
-      FUNCTION SHELL (R,JA)
-
-C-----------------------------------------------------------------------
-C...Density in the shell model
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      COMMON /CSHELL/ RR0(18), RR02(18)
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-      SAVE
-
-      R0 = RR0(JA)
-      C1 = MIN(1.D0,4.D0/DBLE(JA))
-      CS = 1.D0/(R0**3 * PI**1.5D0)
-      CP = 2.D0*CS/3.D0
-      FS = EXP(-(R/R0)**2)
-      FP = (R/R0)**2 * FS
-      SHELL = C1*CS*FS + (1.D0-C1)*CP*FP
-      RETURN
-      END
-C=======================================================================
-
-      FUNCTION PROFNUC (B, JA)
-
-C-----------------------------------------------------------------------
-C...This function return
-C.  the profile T(b) for a nucleus of mass number A
-C.  INPUT B = impact parameter (GeV**-1)
-C.        JA = integer mass number
-C.  OUTPUT  (fm**-2)
-C.
-C.  The  density of the nucleus is the `shell model density'
-C.  the parameter r0 must beinitialized in the common block
-C.............................................................
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      COMMON /CSHELL/ RR0(18), RR02(18)
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-      SAVE
-
-      B2 = B*B
-      ARG = B2/RR02(JA)
-      TS = EXP(-ARG)
-      TP = TS*(2.D0*B2+RR02(JA))/(3.D0*RR02(JA))
-      CS = MIN(1.D0,4.D0/DBLE(JA))
-      PROFNUC = (CS*TS + (1.D0-CS)*TP)/(PI*RR02(JA))
-      RETURN
-      END
-C=======================================================================
-
-      SUBROUTINE SHELL_INI
-
-C-----------------------------------------------------------------------
-C...Initialize the parameter  of the shell model
-C.  for the nuclei with    6 < A < 18
-C..............................................
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-
-      COMMON /CSHELL/ RR0(18), RR02(18)
-      DIMENSION RR(18)
-      SAVE
-C...Data on Sqrt[<r**2>]  in fermi
-      DATA RR /0.81D0, 2.095D0,  1.88D0, 1.674D0,  -1.D0,
-     +         2.56D0, 2.41D0,    -1.D0, 2.519D0, 2.45D0,
-     +         2.37D0, 2.460D0, 2.440D0,  2.54D0, 2.58D0, 
-     +         2.718D0,2.662D0, 2.789D0/
-
-      DO JA=1,18
-         A = DBLE(JA)
-         RMED = RR(JA)
-         IF (RMED .LE. 0.D0)   RMED = 0.5D0*(RR(JA-1) + RR(JA+1))
-         C = MAX(1.5D0,(5.D0/2.D0 - 4.D0/A) )
-         R0 = RMED/SQRT(C)
-         RR0 (JA) = R0
-         RR02(JA) = R0*R0
-      ENDDO
-      RETURN
-      END
-C->
-C=======================================================================
-
-      SUBROUTINE WOOD_SAXON_INI
-
-C-----------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      COMMON /CWOOD/ RR0(19:56), AA0(19:56), CC0(19:56)
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-      SAVE
-
-C...Wood-Saxon parameters from  table 6.2   of Barrett and Jackson
-      RR0 (19) = 2.59D0
-      AA0 (19) = 0.564D0
-      RR0 (20) = 2.74D0
-      AA0 (20) = 0.569D0
-      RR0 (22) = 2.782D0
-      AA0 (22) = 0.549D0
-      RR0 (24) = 2.99D0
-      AA0 (24) = 0.548D0
-      RR0 (27) = 2.84D0
-      AA0 (27) = 0.569D0
-      RR0 (28) = 3.14D0
-      AA0 (28) = 0.537D0
-      RR0 (29) = 3.77D0
-      AA0 (29) = 0.52D0
-      RR0 (48) = 3.912D0
-      AA0 (48) = 0.5234D0
-      RR0 (56) = 3.98D0
-      AA0 (56) = 0.569D0
-      DO J=19, 56
-         IF (RR0(J) .LE. 0.D0)  THEN
-            RR0(J) = 1.05D0*DBLE(J)**0.333333333333D0
-            AA0(J) = 0.545D0
-         ENDIF
-         CC0(J)=3.D0/(4.D0*PI*RR0(J)**3)/(1.D0+((AA0(J)*PI)/RR0(J))**2)
-      ENDDO
-      RETURN
-      END
-C=======================================================================
-
-      FUNCTION DENSA (Z)
-
-C-----------------------------------------------------------------------
-C....Woods Saxon nuclear density (normalised to 1)
-C.   for a nucleus of mass number A.
-C.   INPUT z = z coordinate (fm)
-C.         JA = integer mass number
-C.         B (in common /CC01/)  impact parameter  (fm)
-C.  OUTPUT (fm**-3)
-C--------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      COMMON /CC01/  B
-      COMMON /CCDA/ JA
-      COMMON /CWOOD/ RR0(19:56), AA0(19:56), CC0(19:56)
-      SAVE
-
-      R = SQRT (Z*Z + B*B)
-      DENSA = CC0(JA)/(1.D0+EXP((R-RR0(JA))/AA0(JA)))
-      RETURN
-      END
-
-C*=====================================================================
-C. Cross sections
-C======================================================================
-
-      SUBROUTINE SIGMA_AIR (IB,SIG0,SIGEL,KINT,
-     +                            SIGMA,DSIGMA,SIGQE,DSIGQE)
-
-C-----------------------------------------------------------------------
-C...Compute with a montecarlo method the "production"
-C.  and "quasi-elastic" cross section for  
-C.  a nucleus-air  interaction 
-C.
-C.  INPUT : IB            = mass of projectile nucleus
-C.          SIG0 (mbarn)  = inelastic pp cross section
-C.          KINT            = number  of interactions to generate
-C.  OUTPUT : SIGMA (mbarn) = "production" cross section
-C.           DSIGMA   "    = error
-C.           SIGQE    "    = "quasi-elastic" cross section
-C.           DSIGQE   "    = error
-C.           additional output is in the common block  /CPROBAB/
-C..........................................................................
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-
-      PARAMETER (IAMAX=56)
-      PARAMETER (IAMAX2=3136)          ! IAMAX*IAMAX
-      COMMON  /CPROBAB/ PROBA(IAMAX), DPROBA(IAMAX), 
-     +   PROBB(IAMAX), DPROBB(IAMAX), PROBI(IAMAX2), DPROBI(IAMAX2),
-     +   P1AEL(0:IAMAX),DP1AEL(0:IAMAX),P1BEL(0:IAMAX), DP1BEL(0:IAMAX),
-     +   P2AEL(0:IAMAX),DP2AEL(0:IAMAX),P2BEL(0:IAMAX), DP2BEL(0:IAMAX)
-      COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
-     +         ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
-     +         ,JJAEL(IAMAX), JJBEL(IAMAX)
-      DIMENSION  MMA(0:IAMAX), MMB(0:IAMAX), MMI(0:IAMAX2)
-      DIMENSION  M1AEL(0:IAMAX), M1BEL(0:IAMAX)
-      DIMENSION  M2AEL(0:IAMAX), M2BEL(0:IAMAX)
-      DOUBLE PRECISION FOX
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-      SAVE
-      DATA FOX /0.21522D0/  !atomic percentage of 'non-nitrogen' in air
-
-      R2 = 0.1D0 * SIG0/PI
-      BMAX = 15.D0                           ! fm
-      SIGMA0 = PI*BMAX*BMAX*10.              ! mbarn
-      IA = 16
-      DO J=1,IA
-         MMA(J) = 0
-         M1AEL(J) = 0
-         M2AEL(J) = 0
-      ENDDO
-      DO J=1,IB
-         MMB(J) = 0
-         M1BEL(J) = 0
-         M2BEL(J) = 0
-      ENDDO
-      DO J=1,IA*IB
-         MMI(J) = 0
-      ENDDO
-      NN = 0
-      M = 0
-      DO KK=1,KINT
-c  select target IA from air composition
-         R = S_RNDM(KK)
-         IA = 14
-         IF (R .LT. FOX)  IA = 16
-
-         CALL INT_NUC (IA, IB, SIG0, SIGEL) 
-         NN = NN + NTRY
-         MMI(NI) = MMI(NI) + 1
-         MMA(NA) = MMA(NA)+1
-         MMB(NB) = MMB(NB)+1
-         IF (NI .GT. 0)  THEN
-            M = M+1
-            M1AEL(NAEL) = M1AEL(NAEL)+1
-            M1BEL(NBEL) = M1BEL(NBEL)+1
-         ELSE
-            M2AEL(NAEL) = M2AEL(NAEL)+1
-            M2BEL(NBEL) = M2BEL(NBEL)+1
-         ENDIF
-      ENDDO
-      MQE = KINT - M
-      SIGMA  = SIGMA0 * DBLE(M)/DBLE(NN)
-      DSIGMA = SIGMA0 * SQRT(DBLE(M))/DBLE(NN)
-      SIGQE  = SIGMA0 * DBLE(MQE)/DBLE(NN)
-      DSIGQE = SIGMA0 * SQRT(DBLE(MQE))/DBLE(NN)
-      DO J=1,IA
-         PROBA(J) = DBLE(MMA(J))/DBLE(M)
-         DPROBA(J) = SQRT(DBLE(MMA(J)))/DBLE(M)
-      ENDDO
-      DO J=1,IB
-         PROBB(J) = DBLE(MMB(J))/DBLE(M)
-         DPROBB(J) = SQRT(DBLE(MMB(J)))/DBLE(M)
-      ENDDO
-      DO J=1,IA*IB
-         PROBI(J) = DBLE(MMI(J))/DBLE(M)
-         DPROBI(J) = SQRT(DBLE(MMI(J)))/DBLE(M)
-      ENDDO
-      DO J=0,IA
-         P1AEL(J) = DBLE(M1AEL(J))/DBLE(M)
-         DP1AEL(J) = SQRT(DBLE(M1AEL(J)))/DBLE(M)
-         P2AEL(J) = DBLE(M2AEL(J))/DBLE(MQE)
-         DP2AEL(J) = SQRT(DBLE(M2AEL(J)))/DBLE(MQE)
-      ENDDO
-      DO J=0,IB
-         P1BEL(J) = DBLE(M1BEL(J))/DBLE(M)
-         DP1BEL(J) = SQRT(DBLE(M1BEL(J)))/DBLE(M)
-         P2BEL(J) = DBLE(M2BEL(J))/DBLE(MQE)
-         DP2BEL(J) = SQRT(DBLE(M2BEL(J)))/DBLE(MQE)
-      ENDDO
-      RETURN
-      END
-C->
-C=======================================================================
-
-      SUBROUTINE SIGMA_MC (IA,IB,SIG0,SIGEL,KINT,
-     +                            SIGMA,DSIGMA,SIGQE,DSIGQE)
-
-C-----------------------------------------------------------------------
-C...Compute with a montecarlo method the "production"
-C.  and "quasi-elastic" cross section for  
-C.  a nucleus-nucleus interaction
-C.
-C.  INPUT : IA            = mass of target nucleus
-C.          IB            = mass of projectile nucleus
-C.          SIG0 (mbarn)  = inelastic pp cross section
-C.          KINT            = number  of interactions to generate
-C.  OUTPUT : SIGMA (mbarn) = "production" cross section
-C.           DSIGMA   "    = error
-C.           SIGQE    "    = "quasi-elastic" cross section
-C.           DSIGQE   "    = error
-C.           additional output is in the common block  /CPROBAB/
-C.           Prob(n_A), Prob(n_B), Prob(n_int)
-C..........................................................................
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-
-      PARAMETER (IAMAX=56)
-      PARAMETER (IAMAX2=3136)          ! IAMAX*IAMAX
-      COMMON  /CPROBAB/ PROBA(IAMAX), DPROBA(IAMAX), 
-     +   PROBB(IAMAX), DPROBB(IAMAX), PROBI(IAMAX2), DPROBI(IAMAX2),
-     +   P1AEL(0:IAMAX),DP1AEL(0:IAMAX),P1BEL(0:IAMAX), DP1BEL(0:IAMAX),
-     +   P2AEL(0:IAMAX),DP2AEL(0:IAMAX),P2BEL(0:IAMAX), DP2BEL(0:IAMAX)
-      COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
-     +         ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
-     +         ,JJAEL(IAMAX), JJBEL(IAMAX)
-      DIMENSION  MMA(0:IAMAX), MMB(0:IAMAX), MMI(0:IAMAX2)
-      DIMENSION  M1AEL(0:IAMAX), M1BEL(0:IAMAX)
-      DIMENSION  M2AEL(0:IAMAX), M2BEL(0:IAMAX)
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-      SAVE
-
-      R2 = 0.1D0 * SIG0/PI
-      BMAX = 15.D0                           ! fm
-      SIGMA0 = PI*BMAX*BMAX*10.D0              ! mbarn
-      DO J=1,IA
-         MMA(J) = 0
-         M1AEL(J) = 0
-         M2AEL(J) = 0
-      ENDDO
-      DO J=1,IB
-         MMB(J) = 0
-         M1BEL(J) = 0
-         M2BEL(J) = 0
-      ENDDO
-      DO J=1,IA*IB
-         MMI(J) = 0
-      ENDDO
-      NN = 0
-      M = 0
-      DO KK=1,KINT
-         CALL INT_NUC (IA, IB, SIG0, SIGEL) 
-         NN = NN + NTRY
-         MMI(NI) = MMI(NI) + 1
-         MMA(NA) = MMA(NA)+1
-         MMB(NB) = MMB(NB)+1
-         IF (NI .GT. 0)  THEN
-            M = M+1
-            M1AEL(NAEL) = M1AEL(NAEL)+1
-            M1BEL(NBEL) = M1BEL(NBEL)+1
-         ELSE
-            M2AEL(NAEL) = M2AEL(NAEL)+1
-            M2BEL(NBEL) = M2BEL(NBEL)+1
-         ENDIF
-      ENDDO
-      MQE = KINT - M
-      SIGMA  = SIGMA0 * DBLE(M)/DBLE(NN)
-      DSIGMA = SIGMA0 * SQRT(DBLE(M))/DBLE(NN)
-      SIGQE  = SIGMA0 * DBLE(MQE)/DBLE(NN)
-      DSIGQE = SIGMA0 * SQRT(DBLE(MQE))/DBLE(NN)
-      DO J=1,IA
-         PROBA(J) = DBLE(MMA(J))/DBLE(M)
-         DPROBA(J) = SQRT(DBLE(MMA(J)))/DBLE(M)
-      ENDDO
-      DO J=1,IB
-         PROBB(J) = DBLE(MMB(J))/DBLE(M)
-         DPROBB(J) = SQRT(DBLE(MMB(J)))/DBLE(M)
-      ENDDO
-      DO J=1,IA*IB
-         PROBI(J) = DBLE(MMI(J))/DBLE(M)
-         DPROBI(J) = SQRT(DBLE(MMI(J)))/DBLE(M)
-      ENDDO
-      DO J=0,IA
-         P1AEL(J) = DBLE(M1AEL(J))/DBLE(M)
-         DP1AEL(J) = SQRT(DBLE(M1AEL(J)))/DBLE(M)
-         P2AEL(J) = DBLE(M2AEL(J))/DBLE(MQE)
-         DP2AEL(J) = SQRT(DBLE(M2AEL(J)))/DBLE(MQE)
-      ENDDO
-      DO J=0,IB
-         P1BEL(J) = DBLE(M1BEL(J))/DBLE(M)
-         DP1BEL(J) = SQRT(DBLE(M1BEL(J)))/DBLE(M)
-         P2BEL(J) = DBLE(M2BEL(J))/DBLE(MQE)
-         DP2BEL(J) = SQRT(DBLE(M2BEL(J)))/DBLE(MQE)
-      ENDDO
-      RETURN
-      END
-
-C*=============================================================
-C.  Cross sections
-C*=============================================================
-
-C Glauber h-air cross section calculation moved to inelScreen src file..
-
-C-----------------------------------------------------------------------
-C.  Fit of Block and Cahn to pp and pbar-p cross sections
-C-----------------------------------------------------------------------
-C=======================================================================
-
-      SUBROUTINE BLOCK(SQS,SIG1,SIG2,SLOP1,SLOP2,
-     +                 RHO1,RHO2,SIGEL1,SIGEL2)
-
-C-----------------------------------------------------------------------
-C...p-p and pbar-p cross sections
-C.  Parametrization of  Block and Cahn
-C
-C.  INPUT  : SQS   (GeV)  = c.m. energy
-C.  
-C.  OUPUT : SIG1 (mbarn)    = pp  total  cross section 
-C.          SLOP1 (GeV**2)  = slope of elastic scattering
-C.          RHO1            = Real/Imaginary part of the amplitude
-C.                            for forward elastic  scattering (pp)
-C.          SIGEL1 (mbarn)  = pp  elastic scattering  cross section
-C.          [1 -> 2   : pp -> pbar p]
-C-----------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-      SAVE
-
-      S = SQS*SQS
-      CALL FPLUS  (S, FR, FI)
-      CALL FMINUS (S, GR, GI)
-      SIG1 = FI-GI
-      SIG2 = FI+GI
-      RHO1 = (FR-GR)/(FI-GI)
-      RHO2 = (FR+GR)/(FI+GI)
-      CALL SSLOPE (S, BP, BM)
-      SLOP1 = BP - GI/FI*(BM-BP)
-      SLOP2 = BP + GI/FI*(BM-BP)
-      SIGEL1 = SIG1**2*(1.D0+RHO1**2)/(16.D0*PI*SLOP1)/CMBARN
-      SIGEL2 = SIG2**2*(1.D0+RHO2**2)/(16.D0*PI*SLOP2)/CMBARN
-      RETURN
-      END
-C=======================================================================
-
-      SUBROUTINE FPLUS (S, FR, FI)
-
-C-----------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      COMMON /BLOCKC/ AA, BETA, S0, CC, AMU, DD, ALPHA, A0
-      COMPLEX*16 Z1, Z2, Z3
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-      SAVE
-
-      F1 = LOG(S/S0)
-      Z1 = DCMPLX(F1,-PI/2.D0)
-      Z1 = Z1*Z1
-      Z2 = 1.D0 + A0*Z1
-      Z3 = Z1/Z2
-      F2 = CC*S**(AMU-1.D0)
-      F3 = 0.5D0*PI*(1.-AMU)
-      FI = AA + F2*COS(F3) + BETA*DREAL(Z3)
-      FR = -BETA*DIMAG(Z3)+F2*SIN(F3)
-      RETURN
-      END
-C=======================================================================
-
-      SUBROUTINE FMINUS (S, FR, FI)
-
-C-----------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      COMMON /BLOCKC/ AA, BETA, S0, CC, AMU, DD, ALPHA, A0
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-      SAVE
-
-      F1 = S**(ALPHA-1.D0)
-      F2 = 0.5D0*PI*(1.D0-ALPHA)
-      FR = -DD*F1*COS(F2)
-      FI = -DD*F1*SIN(F2)
-      RETURN
-      END
-C=======================================================================
-
-      SUBROUTINE SSLOPE (S, BP, BM)
-
-C-----------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      COMMON /BLOCKD/ CP, DP, EP, CM, DM
-      SAVE
-
-      AL = LOG(S)
-      BP = CP + DP*AL + EP*AL*AL
-      BM = CM + DM*AL
-      RETURN
-      END
-C=======================================================================
-
-      SUBROUTINE BLOCK_INI
-
-C-----------------------------------------------------------------------
-C...Parameters of fit IFIT=1 of Block and Cahn
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      COMMON /BLOCKC/ AA, BETA, S0, CC, AMU, DD, ALPHA, A0
-      COMMON /BLOCKD/ CP, DP, EP, CM, DM
-      SAVE
-
-      AA = 41.74D0
-      BETA = 0.66D0
-      S0 = 338.5D0
-      CC = 0.D0
-      AMU = 0.D0
-      DD = -39.37D0
-      ALPHA = 0.48D0
-      A0 = 0.D0
-      CP = 10.90D0
-      DP = -0.08D0
-      EP = 0.043D0
-      CM = 23.27D0
-      DM = 0.93D0
-      RETURN
-      END
-
-C*=============================================================
-C.  Nucleus-nucleus cross sections
-C=======================================================================
-
-      SUBROUTINE SIGNUC_INI (IA,E0)
-
-C-----------------------------------------------------------------------
-C...This subroutine receives in INPUT E0 (TeV)
-C.  energy per nucleon and computes the cross sections
-C.  and interactions lengths for  all nuclei
-C.  with A  between 2 and IA
-C.  The output is contained in common block /CLENNN/
-C.
-C.  Attention: the tabulated cross sections are obtained with
-C.  new p-p cross sections as used in SIBYLL 2x,
-C.  in addition field dimensions changed (RE 04/2000)
-C.
-C........................................................
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      COMMON /CLENNN/ SSIGNUC(60), ALNUC(60)
-      DIMENSION SIGMA(6,56), SIGQE(6,56)
-      DIMENSION AA(6)
-      SAVE
-      DATA NE /6/, AMIN /1.D0/, DA /1.D0/
-      DATA AA /1.D0,2.D0,3.D0,4.D0,5.D0,6.D0/
-      DATA AVOG /6.0221367D-04/
-      DATA ATARGET /14.514D0/            ! effective masss of air
-C...Data on `inelastic-production' nucleus-air cross section
-      DATA (SIGMA(J, 2),J=1,6) /
-     &3.842D+02,4.287D+02,4.940D+02,5.887D+02,6.922D+02,7.767D+02/
-      DATA (SIGMA(J, 3),J=1,6) /
-     &4.601D+02,5.149D+02,5.595D+02,6.663D+02,7.641D+02,8.446D+02/
-      DATA (SIGMA(J, 4),J=1,6) /
-     &4.881D+02,5.373D+02,6.005D+02,6.895D+02,7.716D+02,8.967D+02/
-      DATA (SIGMA(J, 5),J=1,6) /
-     &5.874D+02,6.176D+02,7.181D+02,7.993D+02,9.089D+02,1.031D+03/
-      DATA (SIGMA(J, 6),J=1,6) /
-     &7.054D+02,7.399D+02,8.388D+02,9.463D+02,1.080D+03,1.197D+03/
-      DATA (SIGMA(J, 7),J=1,6) /
-     &7.192D+02,7.611D+02,8.449D+02,9.539D+02,1.061D+03,1.176D+03/
-      DATA (SIGMA(J, 8),J=1,6) /
-     &7.550D+02,7.975D+02,9.153D+02,9.944D+02,1.126D+03,1.236D+03/
-      DATA (SIGMA(J, 9),J=1,6) /
-     &7.929D+02,8.392D+02,9.265D+02,1.059D+03,1.167D+03,1.262D+03/
-      DATA (SIGMA(J, 10),J=1,6) /
-     &8.157D+02,8.644D+02,9.512D+02,1.058D+03,1.182D+03,1.298D+03/
-      DATA (SIGMA(J, 11),J=1,6) /
-     &8.039D+02,8.587D+02,9.534D+02,1.055D+03,1.182D+03,1.298D+03/
-      DATA (SIGMA(J, 12),J=1,6) /
-     &8.515D+02,8.957D+02,9.869D+02,1.122D+03,1.253D+03,1.366D+03/
-      DATA (SIGMA(J, 13),J=1,6) /
-     &8.769D+02,9.100D+02,1.018D+03,1.119D+03,1.252D+03,1.341D+03/
-      DATA (SIGMA(J, 14),J=1,6) /
-     &9.058D+02,9.532D+02,1.057D+03,1.171D+03,1.302D+03,1.391D+03/
-      DATA (SIGMA(J, 15),J=1,6) /
-     &9.555D+02,9.799D+02,1.098D+03,1.201D+03,1.342D+03,1.444D+03/
-      DATA (SIGMA(J, 16),J=1,6) /
-     &1.009D+03,1.058D+03,1.149D+03,1.290D+03,1.414D+03,1.520D+03/
-      DATA (SIGMA(J, 17),J=1,6) /
-     &9.907D+02,1.045D+03,1.166D+03,1.290D+03,1.384D+03,1.516D+03/
-      DATA (SIGMA(J, 18),J=1,6) /
-     &1.036D+03,1.121D+03,1.198D+03,1.328D+03,1.470D+03,1.592D+03/
-      DATA (SIGMA(J, 19),J=1,6) /
-     &1.083D+03,1.162D+03,1.250D+03,1.371D+03,1.516D+03,1.661D+03/
-      DATA (SIGMA(J, 20),J=1,6) /
-     &1.146D+03,1.215D+03,1.295D+03,1.443D+03,1.544D+03,1.744D+03/
-      DATA (SIGMA(J, 21),J=1,6) /
-     &1.158D+03,1.234D+03,1.292D+03,1.467D+03,1.618D+03,1.750D+03/
-      DATA (SIGMA(J, 22),J=1,6) /
-     &1.153D+03,1.205D+03,1.329D+03,1.451D+03,1.596D+03,1.734D+03/
-      DATA (SIGMA(J, 23),J=1,6) /
-     &1.210D+03,1.274D+03,1.356D+03,1.493D+03,1.655D+03,1.803D+03/
-      DATA (SIGMA(J, 24),J=1,6) /
-     &1.212D+03,1.273D+03,1.398D+03,1.489D+03,1.641D+03,1.800D+03/
-      DATA (SIGMA(J, 25),J=1,6) /
-     &1.236D+03,1.315D+03,1.423D+03,1.561D+03,1.669D+03,1.855D+03/
-      DATA (SIGMA(J, 26),J=1,6) /
-     &1.279D+03,1.345D+03,1.431D+03,1.595D+03,1.734D+03,1.889D+03/
-      DATA (SIGMA(J, 27),J=1,6) /
-     &1.228D+03,1.304D+03,1.438D+03,1.546D+03,1.714D+03,1.836D+03/
-      DATA (SIGMA(J, 28),J=1,6) /
-     &1.289D+03,1.370D+03,1.451D+03,1.597D+03,1.754D+03,1.913D+03/
-      DATA (SIGMA(J, 29),J=1,6) /
-     &1.411D+03,1.469D+03,1.613D+03,1.777D+03,1.910D+03,2.075D+03/
-      DATA (SIGMA(J, 30),J=1,6) /
-     &1.347D+03,1.401D+03,1.498D+03,1.642D+03,1.816D+03,1.975D+03/
-      DATA (SIGMA(J, 31),J=1,6) /
-     &1.359D+03,1.448D+03,1.551D+03,1.694D+03,1.858D+03,2.007D+03/
-      DATA (SIGMA(J, 32),J=1,6) /
-     &1.358D+03,1.460D+03,1.559D+03,1.698D+03,1.842D+03,1.974D+03/
-      DATA (SIGMA(J, 33),J=1,6) /
-     &1.418D+03,1.448D+03,1.578D+03,1.727D+03,1.872D+03,2.047D+03/
-      DATA (SIGMA(J, 34),J=1,6) /
-     &1.433D+03,1.466D+03,1.605D+03,1.738D+03,1.892D+03,2.019D+03/
-      DATA (SIGMA(J, 35),J=1,6) /
-     &1.430D+03,1.511D+03,1.602D+03,1.752D+03,1.935D+03,2.060D+03/
-      DATA (SIGMA(J, 36),J=1,6) /
-     &1.462D+03,1.499D+03,1.653D+03,1.805D+03,1.920D+03,2.057D+03/
-      DATA (SIGMA(J, 37),J=1,6) /
-     &1.470D+03,1.520D+03,1.656D+03,1.818D+03,1.946D+03,2.131D+03/
-      DATA (SIGMA(J, 38),J=1,6) /
-     &1.470D+03,1.542D+03,1.691D+03,1.800D+03,1.968D+03,2.133D+03/
-      DATA (SIGMA(J, 39),J=1,6) /
-     &1.495D+03,1.588D+03,1.676D+03,1.834D+03,1.969D+03,2.163D+03/
-      DATA (SIGMA(J, 40),J=1,6) /
-     &1.525D+03,1.551D+03,1.722D+03,1.833D+03,2.020D+03,2.192D+03/
-      DATA (SIGMA(J, 41),J=1,6) /
-     &1.526D+03,1.615D+03,1.709D+03,1.899D+03,2.040D+03,2.181D+03/
-      DATA (SIGMA(J, 42),J=1,6) /
-     &1.510D+03,1.567D+03,1.716D+03,1.892D+03,2.056D+03,2.197D+03/
-      DATA (SIGMA(J, 43),J=1,6) /
-     &1.557D+03,1.658D+03,1.776D+03,1.898D+03,2.092D+03,2.200D+03/
-      DATA (SIGMA(J, 44),J=1,6) /
-     &1.556D+03,1.645D+03,1.752D+03,1.920D+03,2.091D+03,2.243D+03/
-      DATA (SIGMA(J, 45),J=1,6) /
-     &1.583D+03,1.663D+03,1.798D+03,1.940D+03,2.051D+03,2.263D+03/
-      DATA (SIGMA(J, 46),J=1,6) /
-     &1.599D+03,1.642D+03,1.799D+03,1.941D+03,2.107D+03,2.268D+03/
-      DATA (SIGMA(J, 47),J=1,6) /
-     &1.611D+03,1.692D+03,1.811D+03,1.956D+03,2.107D+03,2.264D+03/
-      DATA (SIGMA(J, 48),J=1,6) /
-     &1.625D+03,1.706D+03,1.819D+03,1.986D+03,2.139D+03,2.354D+03/
-      DATA (SIGMA(J, 49),J=1,6) /
-     &1.666D+03,1.737D+03,1.854D+03,1.971D+03,2.160D+03,2.318D+03/
-      DATA (SIGMA(J, 50),J=1,6) /
-     &1.648D+03,1.747D+03,1.856D+03,2.023D+03,2.181D+03,2.352D+03/
-      DATA (SIGMA(J, 51),J=1,6) /
-     &1.653D+03,1.763D+03,1.868D+03,2.015D+03,2.203D+03,2.386D+03/
-      DATA (SIGMA(J, 52),J=1,6) /
-     &1.690D+03,1.720D+03,1.902D+03,2.027D+03,2.189D+03,2.357D+03/
-      DATA (SIGMA(J, 53),J=1,6) /
-     &1.690D+03,1.750D+03,1.921D+03,2.059D+03,2.208D+03,2.417D+03/
-      DATA (SIGMA(J, 54),J=1,6) /
-     &1.705D+03,1.781D+03,1.911D+03,2.073D+03,2.242D+03,2.411D+03/
-      DATA (SIGMA(J, 55),J=1,6) /
-     &1.714D+03,1.806D+03,1.896D+03,2.100D+03,2.253D+03,2.411D+03/
-      DATA (SIGMA(J, 56),J=1,6) /
-     &1.774D+03,1.813D+03,1.954D+03,2.098D+03,2.280D+03,2.482D+03/
- 
-      DATA (SIGQE(J, 2),J=1,6) /
-     &4.141D+01,3.708D+01,5.428D+01,8.696D+01,1.403D+02,1.885D+02/
-      DATA (SIGQE(J, 3),J=1,6) /
-     &4.357D+01,3.894D+01,5.177D+01,9.675D+01,1.447D+02,2.029D+02/
-      DATA (SIGQE(J, 4),J=1,6) /
-     &4.123D+01,3.933D+01,6.070D+01,9.482D+01,1.474D+02,2.023D+02/
-      DATA (SIGQE(J, 5),J=1,6) /
-     &4.681D+01,4.287D+01,6.381D+01,1.050D+02,1.519D+02,2.198D+02/
-      DATA (SIGQE(J, 6),J=1,6) /
-     &5.407D+01,5.195D+01,6.723D+01,1.108D+02,1.750D+02,2.368D+02/
-      DATA (SIGQE(J, 7),J=1,6) /
-     &4.975D+01,4.936D+01,6.880D+01,1.162D+02,1.689D+02,2.329D+02/
-      DATA (SIGQE(J, 8),J=1,6) /
-     &5.361D+01,5.027D+01,6.858D+01,1.177D+02,1.759D+02,2.412D+02/
-      DATA (SIGQE(J, 9),J=1,6) /
-     &4.980D+01,5.063D+01,7.210D+01,1.196D+02,1.806D+02,2.299D+02/
-      DATA (SIGQE(J, 10),J=1,6) /
-     &5.170D+01,5.070D+01,7.105D+01,1.182D+02,1.679D+02,2.411D+02/
-      DATA (SIGQE(J, 11),J=1,6) /
-     &4.950D+01,4.950D+01,7.286D+01,1.137D+02,1.769D+02,2.477D+02/
-      DATA (SIGQE(J, 12),J=1,6) /
-     &5.262D+01,5.133D+01,7.110D+01,1.204D+02,1.789D+02,2.501D+02/
-      DATA (SIGQE(J, 13),J=1,6) /
-     &5.320D+01,5.378D+01,6.847D+01,1.200D+02,1.805D+02,2.442D+02/
-      DATA (SIGQE(J, 14),J=1,6) /
-     &5.638D+01,5.271D+01,6.985D+01,1.209D+02,1.867D+02,2.610D+02/
-      DATA (SIGQE(J, 15),J=1,6) /
-     &5.294D+01,5.353D+01,7.435D+01,1.211D+02,1.899D+02,2.612D+02/
-      DATA (SIGQE(J, 16),J=1,6) /
-     &5.668D+01,5.254D+01,7.557D+01,1.269D+02,1.917D+02,2.707D+02/
-      DATA (SIGQE(J, 17),J=1,6) /
-     &5.456D+01,5.721D+01,7.481D+01,1.208D+02,1.859D+02,2.658D+02/
-      DATA (SIGQE(J, 18),J=1,6) /
-     &5.901D+01,5.382D+01,7.591D+01,1.246D+02,1.872D+02,2.874D+02/
-      DATA (SIGQE(J, 19),J=1,6) /
-     &6.328D+01,6.116D+01,8.451D+01,1.318D+02,2.088D+02,2.749D+02/
-      DATA (SIGQE(J, 20),J=1,6) /
-     &5.779D+01,5.924D+01,8.382D+01,1.370D+02,2.062D+02,2.837D+02/
-      DATA (SIGQE(J, 21),J=1,6) /
-     &7.155D+01,5.732D+01,8.231D+01,1.363D+02,2.047D+02,2.820D+02/
-      DATA (SIGQE(J, 22),J=1,6) /
-     &6.699D+01,5.651D+01,8.511D+01,1.477D+02,2.031D+02,2.921D+02/
-      DATA (SIGQE(J, 23),J=1,6) /
-     &6.179D+01,6.269D+01,9.395D+01,1.437D+02,2.195D+02,2.964D+02/
-      DATA (SIGQE(J, 24),J=1,6) /
-     &6.784D+01,6.028D+01,8.622D+01,1.279D+02,2.214D+02,2.867D+02/
-      DATA (SIGQE(J, 25),J=1,6) /
-     &6.589D+01,5.795D+01,8.890D+01,1.385D+02,2.055D+02,2.988D+02/
-      DATA (SIGQE(J, 26),J=1,6) /
-     &6.364D+01,6.325D+01,8.942D+01,1.421D+02,2.128D+02,3.083D+02/
-      DATA (SIGQE(J, 27),J=1,6) /
-     &6.449D+01,6.664D+01,8.986D+01,1.453D+02,2.140D+02,2.932D+02/
-      DATA (SIGQE(J, 28),J=1,6) /
-     &7.284D+01,6.139D+01,8.867D+01,1.425D+02,2.179D+02,2.978D+02/
-      DATA (SIGQE(J, 29),J=1,6) /
-     &7.221D+01,7.085D+01,9.079D+01,1.482D+02,2.277D+02,2.913D+02/
-      DATA (SIGQE(J, 30),J=1,6) /
-     &6.928D+01,6.294D+01,8.935D+01,1.463D+02,2.265D+02,2.834D+02/
-      DATA (SIGQE(J, 31),J=1,6) /
-     &6.611D+01,6.586D+01,9.133D+01,1.461D+02,2.201D+02,2.959D+02/
-      DATA (SIGQE(J, 32),J=1,6) /
-     &6.401D+01,6.177D+01,8.971D+01,1.480D+02,2.155D+02,3.152D+02/
-      DATA (SIGQE(J, 33),J=1,6) /
-     &7.057D+01,6.918D+01,8.410D+01,1.465D+02,2.288D+02,3.088D+02/
-      DATA (SIGQE(J, 34),J=1,6) /
-     &6.453D+01,7.020D+01,9.272D+01,1.517D+02,2.189D+02,2.999D+02/
-      DATA (SIGQE(J, 35),J=1,6) /
-     &6.741D+01,6.295D+01,9.323D+01,1.536D+02,2.190D+02,2.930D+02/
-      DATA (SIGQE(J, 36),J=1,6) /
-     &6.807D+01,7.046D+01,1.025D+02,1.565D+02,2.315D+02,3.090D+02/
-      DATA (SIGQE(J, 37),J=1,6) /
-     &8.082D+01,6.565D+01,9.160D+01,1.572D+02,2.229D+02,3.125D+02/
-      DATA (SIGQE(J, 38),J=1,6) /
-     &6.494D+01,6.964D+01,9.089D+01,1.653D+02,2.336D+02,3.120D+02/
-      DATA (SIGQE(J, 39),J=1,6) /
-     &6.833D+01,6.860D+01,8.933D+01,1.601D+02,2.261D+02,3.167D+02/
-      DATA (SIGQE(J, 40),J=1,6) /
-     &7.021D+01,6.866D+01,8.437D+01,1.588D+02,2.249D+02,2.941D+02/
-      DATA (SIGQE(J, 41),J=1,6) /
-     &7.122D+01,6.205D+01,9.545D+01,1.582D+02,2.335D+02,3.395D+02/
-      DATA (SIGQE(J, 42),J=1,6) /
-     &7.265D+01,6.936D+01,9.486D+01,1.505D+02,2.379D+02,3.248D+02/
-      DATA (SIGQE(J, 43),J=1,6) /
-     &7.048D+01,7.539D+01,9.192D+01,1.566D+02,2.532D+02,3.182D+02/
-      DATA (SIGQE(J, 44),J=1,6) /
-     &6.650D+01,7.139D+01,9.862D+01,1.602D+02,2.289D+02,3.077D+02/
-      DATA (SIGQE(J, 45),J=1,6) /
-     &7.511D+01,6.893D+01,9.245D+01,1.641D+02,2.519D+02,3.381D+02/
-      DATA (SIGQE(J, 46),J=1,6) /
-     &6.437D+01,6.894D+01,8.697D+01,1.544D+02,2.391D+02,3.213D+02/
-      DATA (SIGQE(J, 47),J=1,6) /
-     &7.980D+01,6.958D+01,1.022D+02,1.609D+02,2.408D+02,3.246D+02/
-      DATA (SIGQE(J, 48),J=1,6) /
-     &7.265D+01,7.313D+01,8.989D+01,1.578D+02,2.387D+02,3.235D+02/
-      DATA (SIGQE(J, 49),J=1,6) /
-     &6.959D+01,6.337D+01,9.084D+01,1.656D+02,2.331D+02,3.226D+02/
-      DATA (SIGQE(J, 50),J=1,6) /
-     &7.371D+01,6.807D+01,9.726D+01,1.535D+02,2.445D+02,3.189D+02/
-      DATA (SIGQE(J, 51),J=1,6) /
-     &7.882D+01,6.680D+01,9.377D+01,1.629D+02,2.448D+02,3.297D+02/
-      DATA (SIGQE(J, 52),J=1,6) /
-     &7.223D+01,6.794D+01,9.925D+01,1.738D+02,2.446D+02,3.162D+02/
-      DATA (SIGQE(J, 53),J=1,6) /
-     &7.703D+01,6.971D+01,9.601D+01,1.595D+02,2.484D+02,3.265D+02/
-      DATA (SIGQE(J, 54),J=1,6) /
-     &7.549D+01,7.459D+01,8.984D+01,1.645D+02,2.348D+02,3.201D+02/
-      DATA (SIGQE(J, 55),J=1,6) /
-     &7.891D+01,6.840D+01,1.017D+02,1.698D+02,2.501D+02,3.429D+02/
-      DATA (SIGQE(J, 56),J=1,6) /
-     &7.545D+01,6.673D+01,1.057D+02,1.684D+02,2.424D+02,3.181D+02/
-
-      ASQS = 0.5D0*LOG10(1.876D+03*E0)
-      JE = MIN(INT((ASQS-AMIN)/DA)+1,NE-2)
-      DO JA=2,IA
-         ABEAM = DBLE(JA)
-         S1 = QUAD_INT(ASQS, AA(JE),AA(JE+1),AA(JE+2),
-     +                   SIGMA(JE,JA),SIGMA(JE+1,JA),SIGMA(JE+2,JA))
-         S2 = QUAD_INT(ASQS, AA(JE),AA(JE+1),AA(JE+2),
-     +                   SIGQE(JE,JA),SIGQE(JE+1,JA),SIGQE(JE+2,JA))
-         SSIGNUC(JA) = S1 + S2
-         ALNUC(JA) = ATARGET/(AVOG*SSIGNUC(JA))
-      ENDDO
-      ALNUC(1) = FPNI(E0, 13)
-      SSIGNUC(1) = ATARGET/(AVOG*ALNUC(1))
-
-      RETURN
-      END
-
-
-C*=======================================================================
-C.  General utilities
-C=======================================================================
-
-      FUNCTION QUAD_INT (R,X0,X1,X2,V0,V1,V2)
-
-C-----------------------------------------------------------------------
-C...Quadratic interpolation
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      SAVE
-
-      R0=R-X0
-      R1=R-X1
-      R2=R-X2
-      S0=X0-X1
-      S1=X0-X2
-      S2=X1-X2
-      QUAD_INT = V0*R1*R2/(S0*S1)-V1*R0*R2/(S0*S2)+V2*R0*R1/(S1*S2)
-      RETURN
-      END
-C=======================================================================
-
-      FUNCTION GAUSS (FUN, A,B)
-
-C-----------------------------------------------------------------------
-C...Returns the  8 points Gauss-Legendre integral
-C.  of function FUN from A to B
-C...........................................................
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      DIMENSION X(8), W(8)
-      SAVE
-      DATA X/.0950125098D0, .2816035507D0, .4580167776D0, .6178762444D0,
-     1       .7554044083D0, .8656312023D0, .9445750230D0, .9894009349D0/
-      DATA W/.1894506104D0, .1826034150D0, .1691565193D0, .1495959888D0,
-     1       .1246289712D0, .0951585116D0, .0622535239D0, .0271524594D0/
-
-      XM = 0.5D0*(B+A)
-      XR = 0.5D0*(B-A)
-      SS = 0.D0
-      DO J=1,8
-        DX = XR*X(J)
-        SS = SS + W(J) * (FUN(XM+DX) + FUN(XM-DX))
-      ENDDO
-      GAUSS = XR*SS
-      RETURN
-      END
-C=======================================================================
-
-      SUBROUTINE INVERT_ARRAY (yy, xmin, dx, n, xnew, ymin, dy)
-
-C-----------------------------------------------------------------------
-C..    This subroutine receives one   array
-C      of n y values in input yy(1:n)
-C      that correspond to  equispaced values of x_j = xmin + dx*(j-1)
-C
-C      and "reverse" the array returning an array of  x values
-C      xnew (1:n) that  corresponds to equispaced values of y
-C      The relation is assumed monotonous but can be 
-C      increasing or decreasing
-C..............................................................
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      dimension  yy(n), xnew (n)
-      SAVE
-
-      ymin = yy(1)
-      ymax = yy(n)
-      dy = (ymax - ymin)/float(n-1)
-      xnew (1) = xmin
-      xnew (n) = xmin + dx*float(n-1)
-      k0 = 1
-      do j=2,n-1
-         y = ymin + float(j-1)*dy 
-         do k=k0,n
-            if((yy(k) .gt. y) .eqv. (yy(n) .gt. yy(1))) goto 100
-         enddo
-100      y2 = yy(k)
-         y1 = yy(k-1)
-         k0 = k-1
-         x1 = xmin + dx*float(k-2)
-         x2 = x1+dx
-         xnew (j)  = x1 + dx* (y-y1)/(y2-y1)
-      enddo
-      return
-      end
-C->
-C=======================================================================
-
-      SUBROUTINE SINCO(S,C)
-
-C-----------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-      SAVE
-
-      F = TWOPI*S_RNDM(0)
-      C = COS (F)
-      S = SIN (F)
-      RETURN
-      END
-
-C***********************************************************************
-C.  Cross sections for cascade calculations (FPNI)
-C=======================================================================
-      
-      SUBROUTINE SIGMA_PP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO) 
-
-C-----------------------------------------------------------------------
-C...p-p cross sections
-C.
-C.  this routine serves the purpose to calculate cascades with different 
-C.  cross sections
-C.
-C. INPUT: E0 = Laboratory Energy  (TeV)
-C. 
-C. OUTPUT: SIGT = total cross section
-C.         SIGEL = elastic cross section
-C.         SIGINEL = inelastic cross section
-C.         SLOPE = slope of elastic scattering (GeV**-2)
-C.         RHO = Imaginary/Real part of forward elastic amplitude
-C.   
-C.  (old cross section tables end at 10^6 GeV)
-C-----------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      DIMENSION SSIG0(51)
-      DIMENSION SIGDIF(3)
-      COMMON /CSPA/ ICSPA2(3)
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-      SAVE
-
-C...p-p inelastic cross sections (mbarn)
-      DATA (SSIG0(J),J=1,51) /
-     +      32.05D0,  32.06D0,  32.08D0,  32.13D0,  32.22D0,  32.36D0,
-     +      32.56D0,  32.85D0,  33.24D0,  33.75D0,  34.37D0,  35.14D0,
-     +      36.05D0,  37.12D0,  38.37D0,  39.78D0,  41.36D0,  43.13D0,
-     +      45.07D0,  47.18D0,  49.47D0,  51.91D0,  54.54D0,  57.28D0,
-     +      60.15D0,  63.15D0,  66.28D0,  69.48D0,  72.80D0,  76.22D0,
-     +      79.71D0,  83.27D0,  86.87D0,  90.55D0,  94.26D0,  98.05D0,
-     +     101.89D0, 105.75D0, 109.71D0, 113.65D0, 117.60D0, 121.55D0,
-     +     125.53D0, 129.56D0, 133.60D0, 137.70D0, 141.77D0, 145.84D0,
-     +     149.92D0, 154.02D0, 158.15D0/
-
-      ICSPA = ICSPA2(1)
-
-      SQS = SQRT(2000.D0*0.938D0*E0)
-
-*  pre-LHC SIBYLL2.1 model
-      
-      IF(ICSPA.EQ.-2) THEN
-
-         CALL SIB_SIGMA_EXT(3,SQS,SIGT,SIGEL,SIGINEL,SLOPE,RHO)      
-
-*  old standard NUCLIB/SIBYLL model
-
-      ELSE IF(ICSPA.EQ.-1) THEN
-
-        AL = LOG10(SQS)
-        if(AL.le.1.D0) then
-          SIGINEL = SSIG0(1)
-        else
-          J1 = INT((AL - 1.D0)*10.D0) + 1
-          J1 = min(J1,50)
-          T = (AL-1.D0)*10.D0 - DBLE(J1-1)
-          SIGINEL = SSIG0(J1)*(1.D0-T) + SSIG0(J1+1)*T
-        endif
-        CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,SIGEL1,SIGEL2)
-        R = SIGEL1/SIGT1
-        RHO = RHO1
-        SIGT  = SIGINEL/(1.D0-R)
-        SIGEL = SIGINEL*R/(1.D0-R)
-        SLOPE = SIGT**2/(SIGEL * 16.D0*PI) * (1.D0+RHO1**2) /CMBARN
-
-*  cross section as calculated in SIBYLL
-
-      ELSE IF(ICSPA.EQ.0) THEN
-
-        CALL SIB_SIGMA_HP(1,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
-
-*  Donnachie-Landshoff  (sig-tot)
-
-      ELSE IF(ICSPA.EQ.1) THEN
-
-        CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,
-     +             SIGEL1,SIGEL2)
-        R = SIGEL1/SIGT1
-        RHO = RHO1
-
-        DELDL = 0.0808D0
-        EPSDL = -0.4525D0
-        S = SQS*SQS
-        SIGT = 21.7D0*S**DELDL+56.08D0*S**EPSDL
-        SIGEL = R*SIGT
-        SIGINEL = SIGT-SIGEL
-        SLOPE = SIGT**2/(SIGEL * 16.D0*PI) * (1.D0+RHO**2) /CMBARN
-
-*  Donnachie-Landshoff (sig-tot and sig-el)
-
-      ELSE IF(ICSPA.EQ.2) THEN
-
-        DELDL = 0.0808D0
-        EPSDL = -0.4525D0
-        S = SQS*SQS
-        SIGT = 21.7D0*S**DELDL+56.08D0*S**EPSDL
-        IMODEL = 1
-        IF(IMODEL.EQ.1) THEN
-          ALPHAP = 0.25D0
-          SLOPE = 8.5D0+2.D0*ALPHAP*LOG(S)
-        ELSE IF(IMODEL.EQ.2) THEN
-          ALPHAP = 0.3D0
-          SLOPE = 8.D0+2.D0*ALPHAP*LOG(S)
-        ENDIF
-        SIGEL = SIGT**2/(16.D0*PI*SLOPE*CMBARN)
-        SIGINEL = SIGT-SIGEL
-        RHO = 0.D0
-
-*  geometrical scaling with Donnachie-Landshoff sig-tot
-
-      ELSE IF(ICSPA.EQ.3) THEN
-
-        R = 0.17D0
-
-        DELDL = 0.0808D0
-        EPSDL = -0.4525D0
-        S = SQS*SQS
-        SIGT = 21.7D0*S**DELDL+56.08D0*S**EPSDL
-
-        SIGEL = R*SIGT
-        SIGINEL = SIGT-SIGEL
-        SLOPE = SIGT**2/(16.D0*PI*SIGEL)/CMBARN
-        RHO = 0.D0
-
-c ICSPA=4 reserved for CONEX_EXTENSION
-c      ELSE IF(ICSPA.EQ.4) THEN
-
-*  cross section from 2014 Review of Particle Physics
-        
-      ELSE IF(ICSPA.EQ.5) THEN
-         
-c     elastic slope not included in fit
-c     taking slope parameterization from sigma_pp Donnie.-Landshoff
-         ALPHAP = 0.25D0
-         SLOPE = 8.5D0+4.D0*ALPHAP*LOG(SQS)
-         
-         CALL SIG_RPP2014(1,1,SQS,SLOPE,SIGT,SIGEL,SIGINEL,RHO)
-
-      ENDIF
-
-      RETURN
-      END
-
-C=======================================================================
-
-      SUBROUTINE SIGMA_PIP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO) 
-
-C-----------------------------------------------------------------------
-C...pi-p cross sections
-C.
-C.  this routine serves the purpose to calculate cascades with different 
-C.  cross sections
-C.
-C. INPUT: E0 = Laboratory Energy  (TeV)
-C. 
-C. OUTPUT: SIGT = total cross section
-C.         SIGEL = elastic cross section
-C.         SIGINEL = inelastic cross section
-C.         SLOPE = slope of elastic scattering (GeV**-2)
-C.         RHO = Imaginary/Real part of forward elastic amplitude
-C.
-C.  (old cross section tables end at 10^6 GeV)
-C-----------------------------------------------------------------------
-Cf2py double precision,intent(in) :: e0
-Cf2py double precision,intent(out) :: sigt, sigel, siginel, slope, rho
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      DIMENSION SSIG0(51)
-      DIMENSION SIGDIF(3)
-      COMMON /CSPA/ ICSPA2(3)
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-      SAVE
-
-C...pi-p inelastic cross sections (mbarn)
-      DATA (SSIG0(J),J=1,51) /
-     +      20.76D0,  20.78D0,  20.81D0,  20.88D0,  20.98D0,  21.13D0,
-     +      21.33D0,  21.61D0,  21.96D0,  22.39D0,  22.92D0,  23.56D0,
-     +      24.31D0,  25.18D0,  26.18D0,  27.32D0,  28.60D0,  30.04D0,
-     +      31.64D0,  33.40D0,  35.34D0,  37.43D0,  39.72D0,  42.16D0,
-     +      44.77D0,  47.56D0,  50.53D0,  53.66D0,  56.99D0,  60.50D0,
-     +      64.17D0,  68.03D0,  72.05D0,  76.27D0,  80.67D0,  85.27D0,
-     +      90.08D0,  95.04D0, 100.27D0, 105.65D0, 111.21D0, 116.94D0,
-     +     122.87D0, 129.03D0, 135.37D0, 141.93D0, 148.62D0, 155.49D0,
-     +     162.48D0, 169.60D0, 176.94D0/
-
-      ICSPA = ICSPA2(2)
-
-      SQS = SQRT(2000.D0*0.938D0*E0)
-      
-*  pre-LHC SIBYLL2.1 model
-      
-      IF(ICSPA.EQ.-2) THEN
-
-         CALL SIB_SIGMA_EXT(2,SQS,SIGT,SIGEL,SIGINEL,SLOPE,RHO)
-      
-*  old standard NUCLIB/SIBYLL model
-
-      ELSE IF(ICSPA.EQ.-1) THEN
-
-        AL = LOG10(SQS)
-        if(AL.le.1.D0) then
-          SIGINEL = SSIG0(1)
-        else
-          J1 = INT((AL - 1.D0)*10.D0) + 1
-          J1 = min(J1,50)
-          T = (AL-1.D0)*10.D0 - DBLE(J1-1)
-          SIGINEL = SSIG0(J1)*(1.D0-T) + SSIG0(J1+1)*T
-        endif
-        CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,SIGEL1,SIGEL2)
-        R = SIGEL1/SIGT1
-        RHO = RHO1
-        SIGT  = SIGINEL/(1.D0-R)
-        SIGEL = SIGINEL*R/(1.D0-R)
-        SLOPE = SIGT**2/(SIGEL * 16.D0*PI) * (1.D0+RHO1**2) /CMBARN
-
-*  cross section as calculated in SIBYLL
-
-      ELSE IF(ICSPA.EQ.0) THEN
-
-        CALL SIB_SIGMA_HP(2,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
-
-*  Donnachie-Landshoff  (sig-tot)
-
-      ELSE IF(ICSPA.EQ.1) THEN
-
-        CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,
-     +             SIGEL1,SIGEL2)
-        R = SIGEL1/SIGT1
-        RHO = RHO1
-
-        DELDL = 0.0808D0
-        EPSDL = -0.4525D0
-        S = SQS*SQS
-        SIGT = 13.63D0*S**DELDL+(36.02D0+27.56D0)/2.D0*S**EPSDL
-        SIGEL = R*SIGT
-        SIGINEL = SIGT-SIGEL
-        SLOPE = SIGT**2/(SIGEL * 16.D0*PI) * (1.D0+RHO**2) /CMBARN
-
-*  Donnachie-Landshoff (sig-tot and sig-el)
-
-      ELSE IF(ICSPA.EQ.2) THEN
-
-        DELDL = 0.0808D0
-        EPSDL = -0.4525D0
-        S = SQS*SQS
-        SIGT = 13.63D0*S**DELDL+(36.02D0+27.56D0)/2.D0*S**EPSDL
-        IMODEL = 1
-        IF(IMODEL.EQ.1) THEN
-          ALPHAP = 0.25D0
-          SLOPE = 8.5D0+2.D0*ALPHAP*LOG(S)
-        ELSE IF(IMODEL.EQ.2) THEN
-          ALPHAP = 0.3D0
-          SLOPE = 8.D0+2.D0*ALPHAP*LOG(S)
-        ENDIF
-        SIGEL = SIGT**2/(16.D0*PI*SLOPE*CMBARN)
-        SIGINEL = SIGT-SIGEL
-        RHO = 0.
-
-*  geometrical scaling with Donnachie-Landshoff sig-tot
-
-      ELSE IF(ICSPA.EQ.3) THEN
-
-        R = 0.17D0
-
-        DELDL = 0.0808D0
-        EPSDL = -0.4525D0
-        S = SQS*SQS
-        SIGT = 13.63D0*S**DELDL+(36.02D0+27.56D0)/2.D0*S**EPSDL
-
-        SIGEL = R*SIGT
-        SIGINEL = SIGT-SIGEL
-        SLOPE = SIGT**2/(16.D0*PI*SIGEL)/CMBARN
-        RHO = 0.D0
-
-c ICSPA=4 reserved for CONEX_EXTENSION
-c      ELSE IF(ICSPA.EQ.4) THEN
-
-*  cross section from 2014 Review of Particle Physics
-        
-      ELSE IF(ICSPA.EQ.5) THEN
-         
-c     elastic slope not included in fit
-c     taking slope parameterization from sigma_pp Donnie.-Landshoff
-         ALPHAP = 0.25D0
-         SLOPE = 8.5D0+4.D0*ALPHAP*LOG(SQS)
-         
-         CALL SIG_RPP2014(2,1,SQS,SLOPE,SIGT,SIGEL,SIGINEL,RHO)
-         
-      ENDIF
-
-
-      RETURN
-      END
-
-C=======================================================================
-
-      SUBROUTINE SIGMA_KP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO) 
-
-C-----------------------------------------------------------------------
-C...K-p cross sections
-C.
-C.  this routine serves the purpose to calculate cascades with different 
-C.  cross sections
-C.
-C.  if old cross sections are selected then sigma_pi = sigma_K
-C.
-C. INPUT: E0 = Laboratory Energy  (TeV)
-C. 
-C. OUTPUT: SIGT = total cross section
-C.         SIGEL = elastic cross section
-C.         SIGINEL = inelastic cross section
-C.         SLOPE = slope of elastic scattering (GeV**-2)
-C.         RHO = Imaginary/Real part of forward elastic amplitude
-C.
-C.  (old cross section tables end at 10^6 GeV)
-C-----------------------------------------------------------------------
-Cf2py double precision,intent(in) :: e0
-Cf2py double precision,intent(out) :: sigt, sigel, siginel, slope, rho
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      DIMENSION SSIG0(51)
-      DIMENSION SIGDIF(3)
-      COMMON /CSPA/ ICSPA2(3)
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-      SAVE
-
-C...pi-p inelastic cross sections (mbarn)
-      DATA (SSIG0(J),J=1,51) /
-     +      20.76D0,  20.78D0,  20.81D0,  20.88D0,  20.98D0,  21.13D0,
-     +      21.33D0,  21.61D0,  21.96D0,  22.39D0,  22.92D0,  23.56D0,
-     +      24.31D0,  25.18D0,  26.18D0,  27.32D0,  28.60D0,  30.04D0,
-     +      31.64D0,  33.40D0,  35.34D0,  37.43D0,  39.72D0,  42.16D0,
-     +      44.77D0,  47.56D0,  50.53D0,  53.66D0,  56.99D0,  60.50D0,
-     +      64.17D0,  68.03D0,  72.05D0,  76.27D0,  80.67D0,  85.27D0,
-     +      90.08D0,  95.04D0, 100.27D0, 105.65D0, 111.21D0, 116.94D0,
-     +     122.87D0, 129.03D0, 135.37D0, 141.93D0, 148.62D0, 155.49D0,
-     +     162.48D0, 169.60D0, 176.94D0/
-
-      ICSPA = ICSPA2(3)
-      
-      SQS = SQRT(2000.D0*0.938D0*E0)      
-
-*  pre-LHC SIBYLL2.1 model
-      
-      IF(ICSPA.EQ.-2) THEN
-
-         CALL SIB_SIGMA_EXT(3,SQS,SIGT,SIGEL,SIGINEL,SLOPE,RHO)
-      
-*  old standard NUCLIB/SIBYLL model
-
-      ELSE IF(ICSPA.EQ.-1) THEN
-
-        AL = LOG10(SQS)
-        if(AL.le.1.D0) then
-          SIGINEL = SSIG0(1)
-        else
-          J1 = INT((AL - 1.D0)*10.D0) + 1
-          J1 = min(J1,50)
-          T = (AL-1.D0)*10.D0 - DBLE(J1-1)
-          SIGINEL = SSIG0(J1)*(1.D0-T) + SSIG0(J1+1)*T
-        endif
-        CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,SIGEL1,SIGEL2)
-        R = SIGEL1/SIGT1
-        RHO = RHO1
-        SIGT  = SIGINEL/(1.D0-R)
-        SIGEL = SIGINEL*R/(1.D0-R)
-        SLOPE = SIGT**2/(SIGEL * 16.D0*PI) * (1.D0+RHO1**2) /CMBARN
-
-*  cross section as calculated in SIBYLL
-
-      ELSE IF(ICSPA.EQ.0) THEN
-
-        CALL SIB_SIGMA_HP(3,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
-
-*  Donnachie-Landshoff  (sig-tot)
-
-      ELSE IF(ICSPA.EQ.1) THEN
-
-        CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,
-     +             SIGEL1,SIGEL2)
-        R = SIGEL1/SIGT1
-        RHO = RHO1
-
-        DELDL = 0.0808D0
-        EPSDL = -0.4525D0
-        S = SQS*SQS
-        SIGT = 11.82D0*S**DELDL+(26.36D0+ 8.15D0)/2.D0*S**EPSDL
-        SIGEL = R*SIGT
-        SIGINEL = SIGT-SIGEL
-        SLOPE = SIGT**2/(SIGEL * 16.D0*PI) * (1.D0+RHO**2) /CMBARN
-
-*  Donnachie-Landshoff (sig-tot and sig-el)
-
-      ELSE IF(ICSPA.EQ.2) THEN
-
-        DELDL = 0.0808D0
-        EPSDL = -0.4525D0
-        S = SQS*SQS
-        SIGT = 11.82D0*S**DELDL+(26.36D0+ 8.15D0)/2.D0*S**EPSDL
-        IMODEL = 1
-        IF(IMODEL.EQ.1) THEN
-          ALPHAP = 0.25D0
-          SLOPE = 8.5D0+2.D0*ALPHAP*LOG(S)
-        ELSE IF(IMODEL.EQ.2) THEN
-          ALPHAP = 0.3D0
-          SLOPE = 8.D0+2.D0*ALPHAP*LOG(S)
-        ENDIF
-        SIGEL = SIGT**2/(16.D0*PI*SLOPE*CMBARN)
-        SIGINEL = SIGT-SIGEL
-        RHO = 0.D0
-
-*  geometrical scaling with Donnachie-Landshoff sig-tot
-
-      ELSE IF(ICSPA.EQ.3) THEN
-
-        R = 0.17D0
-
-        DELDL = 0.0808D0
-        EPSDL = -0.4525D0
-        S = SQS*SQS
-        SIGT = 11.82D0*S**DELDL+(26.36D0+ 8.15D0)/2.D0*S**EPSDL
-
-        SIGEL = R*SIGT
-        SIGINEL = SIGT-SIGEL
-        SLOPE = SIGT**2/(16.D0*PI*SIGEL)/CMBARN
-        RHO = 0.D0
-        
-c ICSPA=4 reserved for CONEX_EXTENSION
-c      ELSE IF(ICSPA.EQ.4) THEN
-
-
-*  cross section from 2014 Review of Particle Physics
-        
-      ELSE IF(ICSPA.EQ.5) THEN
-         
-c     elastic slope not included in fit
-c     taking slope parameterization from sigma_pp Donnie.-Landshoff
-         ALPHAP = 0.25D0
-         SLOPE = 8.5D0+4.D0*ALPHAP*LOG(SQS)
-         
-         CALL SIG_RPP2014(3,1,SQS,SLOPE,SIGT,SIGEL,SIGINEL,RHO)
-
-      ENDIF
-
-      RETURN
-      END
-
-C=======================================================================
-
-      SUBROUTINE SIGMA_INI 
-
-C-----------------------------------------------------------------------
-C.  Initialize the cross section and interaction lengths in air
-C.  cross section model can be chosen, per particle, by setting ICSPA2()
-C.  default is Sibyll cross section (0,0,0)      
-C-----------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      COMMON /CSAIR/ ASQSMIN, ASQSMAX, DASQS,
-     &     SSIG0(61,3),SSIGA(61,3),ALINT(61,3),NSQS
-      
-      COMMON /CSPA/ ICSPA2(3)
-
-      INTEGER NCALL, NDEBUG, LUN
-      COMMON /S_DEBUG/ NCALL, NDEBUG, LUN
-
-C--------------------------------------------------------------------
-C     SIBYLL utility common blocks containing constants       \FR'14
-C--------------------------------------------------------------------
-      DOUBLE PRECISION EPS3,EPS5,EPS8,EPS10
-      COMMON /SIB_EPS/ EPS3,EPS5,EPS8,EPS10
-
-      DOUBLE PRECISION PI,TWOPI,CMBARN
-      COMMON /SIB_CST/ PI,TWOPI,CMBARN
-
-      DOUBLE PRECISION FACN
-      DIMENSION FACN(3:10)
-      COMMON /SIB_FAC/ FACN
-      SAVE
-      DATA ICSPA2 /0,0,0/
-      DATA AVOG /6.0221367D-04/
-      DATA ATARGET /14.514D0/            ! effective masss of air
-
-      IF(NDEBUG.gt.0)
-     &     write(lun,*) ' SIGMA_INI: using cross section model no.',
-     &     (ICSPA2(i),i=1,3)
-
-      CALL BLOCK_INI
-
-C...Loop on c.m. energy 
-      NSQS = 61
-      SQSMIN = 10.D0
-      SQSMAX = 1.d+07
-      ASQSMIN = LOG10(SQSMIN)
-      ASQSMAX = LOG10(SQSMAX)
-      DASQS = (ASQSMAX-ASQSMIN)/DBLE(NSQS-1)
-      DO J=1,NSQS
-         ASQS = ASQSMIN + DASQS*DBLE(J-1)
-         SQS = 10.D0**ASQS
-         E0 = SQS*SQS/(2.D0*0.938D0) * 1.D-03       ! TeV
-C...p-air
-         CALL SIGMA_PP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
-C     using parametrization by Goulianos for diff. cross section
-c     (depends on elastic cross section)
-c     used to determine coupling to intermediate resonances in Glauber calc (ALAM)
-c     assumed to be universal, i.e. same coupling used for proton, pion and kaons
-         CALL SIB_HADCS1(1,SQS,SIGT1,SIGEL1,SIGINEL1,SLOPE1,RHO1)
-         SIGEFF = 0.68D0*(1.D0+36.D0/SQS**2)
-     &        *LOG(0.6D0+0.02D0/1.5D0*SQS**2)
-         SIGEFF = MAX(0.D0,SIGEFF)
-         ALAM = sqrt(SIGEFF/SIGEL1)
-         SSIGSD = 2.D0 * SIGEFF        
-         CALL SIG_H_AIR (SIGT, SLOPE, RHO, ALAM,
-     &        SSIGT, SSIGEL, SSIGQE, SIGSD, SIGQSD )
-         SSIGA(J,1) = SSIGT-SSIGQE ! had-air production cross section
-         SSIG0(J,1) = SIGINEL   ! had-nucleon inel. cross section
-         ALINT(J,1) = 1.D0/(AVOG*SSIGA(J,1)/ATARGET) ! interaction length in air
-C...pi-air
-         CALL SIGMA_PIP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO) 
-         CALL  SIG_H_AIR (SIGT, SLOPE, RHO, ALAM,
-     &        SSIGT, SSIGEL, SSIGQE, SIGSD, SIGQSD )
-         SSIGA(J,2) = SSIGT-SSIGQE
-         SSIG0(J,2) = SIGINEL
-         ALINT(J,2) = 1.D0/(AVOG*SSIGA(J,2)/ATARGET)
-C...K-air
-         CALL SIGMA_KP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO) 
-         CALL  SIG_H_AIR (SIGT, SLOPE, RHO, ALAM,
-     &        SSIGT, SSIGEL, SSIGQE, SIGSD, SIGQSD )
-         SSIGA(J,3) = SSIGT-SSIGQE
-         SSIG0(J,3) = SIGINEL
-         ALINT(J,3) = 1.D0/(AVOG*SSIGA(J,3)/ATARGET)
-      ENDDO
-
-      if (ndebug .gt. 0 ) THEN
-        WRITE(LUN,'(1X,A)') 
-     &  ' SIGMA_INI: NUCLIB interaction lengths [g/cm**2]'
-        WRITE(LUN,'(1X,A)') 
-     &  '     sqs,       p-air,      pi-air,     K-air'
-      DO J=1,NSQS
-         SQS = 10.D0**(ASQSMIN + DASQS*DBLE(J-1))
-         WRITE(LUN,'(1X,1P,4E12.3)') 
-     &        SQS,ALINT(J,1),ALINT(J,2),ALINT(J,3)
-        ENDDO
-      endif
-
-      RETURN
-      END
-
-C=======================================================================
-
-      FUNCTION FPNI (E,Linp)
-
-C-----------------------------------------------------------------------
-C...This function  returns the interaction length 
-C.  of an hadronic particle travelling in air
-C.
-C.  INPUT:   E (TeV)   particle energy
-C.           Linp      particle code
-C.  OUTPUT:  FPNI      (g cm-2)
-C-----------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-            
-      COMMON /CSAIR/ ASQSMIN, ASQSMAX, DASQS,
-     &     SSIG0(61,3),SSIGA(61,3),ALINT(61,3),NSQS
-
-      DIMENSION KK(6:14)
-      SAVE
-      DATA KK /3*2, 4*3, 2*1/
-
-      SQS = SQRT(2000.D0*E*0.937D0)                        ! GeV
-      AL = LOG10 (SQS)
-      L = abs(Linp)
-      IF (AL .LE. ASQSMIN)  THEN
-         FPNI = ALINT(1,KK(L))
-      ELSE
-         T = (AL-ASQSMIN)/DASQS
-         J = INT(T)
-         J = MIN(J,NSQS-2)
-         T = T-DBLE(J)
-         FPNI = ((1.D0-T)*ALINT(J+1,KK(L)) + T*ALINT(J+2,KK(L)))
-      ENDIF
-      RETURN
-      END
-
-C=======================================================================
-      
-      FUNCTION FSIGHAIR (E,Linp)
-
-C-----------------------------------------------------------------------
-C...This function returns the production cross section
-C.  of an hadronic particle with air calculated in NUCLIB (SIGMA_INI)     
-C.
-C.  INPUT:   E (TeV)   particle energy
-C.           Linp      particle code
-C.  OUTPUT:  SIG_PROD  (mb)
-C-----------------------------------------------------------------------
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-            
-      COMMON /CSAIR/ ASQSMIN, ASQSMAX, DASQS,
-     &     SSIG0(61,3),SSIGA(61,3),ALINT(61,3),NSQS
-
-      DIMENSION KK(6:14)
-      SAVE
-      DATA KK /3*2, 4*3, 2*1/
-
-      SQS = SQRT(2000.D0*E*0.937D0)                        ! GeV
-      AL = LOG10 (SQS)
-      L = abs(Linp)
-      IF (AL .LE. ASQSMIN)  THEN
-         FSIGHAIR = SSIGA(1,KK(L))
-      ELSE
-         T = (AL-ASQSMIN)/DASQS
-         J = INT(T)
-         J = MIN(J,NSQS-2)
-         T = T-DBLE(J)
-         FSIGHAIR = ((1.D0-T)*SSIGA(J+1,KK(L)) + T*SSIGA(J+2,KK(L)))
-      ENDIF     
-      RETURN
-      END
-
-C=======================================================================
-
-      SUBROUTINE INT_LEN_INI
-
-C-----------------------------------------------------------------------
-C...Initialize the interaction lengths from NUCLIB
-C-----------------------------------------------------------------------
-      SAVE
-      
-      CALL NUC_GEOM_INI                 ! nucleus profiles
-      CALL SIGMA_INI                    ! initialize cross sections
-
-      RETURN
-      END
-C=======================================================================
 
       SUBROUTINE TRANSFONSHELL(ECM,XM1in,XM2in,XMAX,IMOD,P1,P2,LBAD)