From e538fcce75f589dce53dabc947eaa4bfb3ab78ba Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Tue, 11 Feb 2020 15:48:12 +0100
Subject: [PATCH] updated sibyll2.3c.f to 2.3c03
---
Processes/Sibyll/sibyll2.3c.f | 3634 ++++++++++++++++++++++++++++++++-
1 file changed, 3598 insertions(+), 36 deletions(-)
diff --git a/Processes/Sibyll/sibyll2.3c.f b/Processes/Sibyll/sibyll2.3c.f
index 1a5a80f02..e757320c8 100644
--- a/Processes/Sibyll/sibyll2.3c.f
+++ b/Processes/Sibyll/sibyll2.3c.f
@@ -7,7 +7,7 @@ C SSSSSS IIIIIII BBBBB YY LLLLLLL LLLLLLL
C=======================================================================
C Code for SIBYLL: hadronic interaction Monte Carlo event generator
C=======================================================================
-C Version 2.3c02 (Jun-01-2017, modified Dec-11-2017)
+C Version 2.3c03 (Jun-01-2017, modified Aug-22-2019)
C
C with CHARM production
C
@@ -32,6 +32,8 @@ C friehn@lip.pt
C stanev@bartol.udel.edu
C
C last changes relative to Sibyll 2.3c:
+C * added cross section tables for hadron-nitrogen and hadron-oxygen
+C (changed S_CCSIG common)
C * no remnant in high mass diff. events (pi0-had scattering)
C * repaired had-nuc. cross section routine for kaon beams
C routine remains inactive in ordinary calls.
@@ -469,7 +471,7 @@ C-----------------------------------------------------------------------
* /,' ','| F. RIEHN et al., Proc. 35th Int. Cosmic Ray Conf.|',
* /,' ','| Bexco, Busan, Korea, cont. 301 (2017) |',
* /,' ','| |',
- * /,' ','| last modifications: F. Riehn (12/11/2017) |',
+ * /,' ','| last modifications: F. Riehn (08/22/2019) |',
* /,' ','====================================================',
* /)
@@ -480,7 +482,7 @@ C-----------------------------------------------------------------------
CALL BLOCK_INI
CALL NUC_GEOM_INI
CALL SIG_AIR_INI
-c CALL DEC_INI
+ CALL DEC_INI
c... charm frag. normalisation
CALL ZNORMAL
@@ -4453,8 +4455,11 @@ C--------------------------------------------------------------------
DOUBLE PRECISION FACN
DIMENSION FACN(3:10)
COMMON /SIB_FAC/ FACN
+
+ DOUBLE PRECISION EPS100
SAVE
-
+ DATA EPS100 /1.D-100/
+
DO J=0,NH_max
DO I=0,NS_max
P_int(I,J) = 0.D0
@@ -4535,16 +4540,20 @@ C--------------------------------------------------------------------
chi2_hard_12 = (1.D0-al1+ga1)*(1.D0-al2-ga2)*chi2_hard
chi2_hard_21 = (1.D0-al1-ga1)*(1.D0-al2+ga2)*chi2_hard
+ ef_11 = max(-0.5D0*(chi2_soft_11+chi2_hard_11),log(EPS100))
+ ef_22 = max(-0.5D0*(chi2_soft_22+chi2_hard_22),log(EPS100))
+ ef_12 = max(-0.5D0*(chi2_soft_12+chi2_hard_12),log(EPS100))
+ ef_21 = max(-0.5D0*(chi2_soft_21+chi2_hard_21),log(EPS100))
+
+ ef_11 = dexp(ef_11)
+ ef_22 = dexp(ef_22)
+ ef_12 = dexp(ef_12)
+ ef_21 = dexp(ef_21)
- ef_11 = dexp(-0.5D0*(chi2_soft_11+chi2_hard_11))
- ef_22 = dexp(-0.5D0*(chi2_soft_22+chi2_hard_22))
- ef_12 = dexp(-0.5D0*(chi2_soft_12+chi2_hard_12))
- ef_21 = dexp(-0.5D0*(chi2_soft_21+chi2_hard_21))
-
- esf_11 = ef_11**2
- esf_22 = ef_22**2
- esf_12 = ef_12**2
- esf_21 = ef_21**2
+ esf_11 = max(ef_11,EPS100)**2
+ esf_22 = max(ef_22,EPS100)**2
+ esf_12 = max(ef_12,EPS100)**2
+ esf_21 = max(ef_21,EPS100)**2
F_ine = B*(1.D0 - fe_11*esf_11 - fe_12*esf_12
& - fe_21*esf_21 - fe_22*esf_22)
@@ -4578,24 +4587,24 @@ C--------------------------------------------------------------------
soft_rec_22 = 1.D0/chi2_soft_22
soft_rec_12 = 1.D0/chi2_soft_12
soft_rec_21 = 1.D0/chi2_soft_21
- chi2_hard_11 = max(chi2_hard_11,EPS10)
- chi2_hard_22 = max(chi2_hard_22,EPS10)
- chi2_hard_12 = max(chi2_hard_12,EPS10)
- chi2_hard_21 = max(chi2_hard_21,EPS10)
+ chi2_hard_11 = max(chi2_hard_11,EPS100)
+ chi2_hard_22 = max(chi2_hard_22,EPS100)
+ chi2_hard_12 = max(chi2_hard_12,EPS100)
+ chi2_hard_21 = max(chi2_hard_21,EPS100)
DO I=0,I0MAX
- soft_rec_11 = soft_rec_11*chi2_soft_11
- soft_rec_22 = soft_rec_22*chi2_soft_22
- soft_rec_12 = soft_rec_12*chi2_soft_12
- soft_rec_21 = soft_rec_21*chi2_soft_21
- hard_rec_11 = 1.D0/chi2_hard_11
- hard_rec_22 = 1.D0/chi2_hard_22
- hard_rec_12 = 1.D0/chi2_hard_12
- hard_rec_21 = 1.D0/chi2_hard_21
+ soft_rec_11 = max(soft_rec_11*chi2_soft_11,EPS100)
+ soft_rec_22 = max(soft_rec_22*chi2_soft_22,EPS100)
+ soft_rec_12 = max(soft_rec_12*chi2_soft_12,EPS100)
+ soft_rec_21 = max(soft_rec_21*chi2_soft_21,EPS100)
+ hard_rec_11 = max(1.D0/chi2_hard_11,EPS100)
+ hard_rec_22 = max(1.D0/chi2_hard_22,EPS100)
+ hard_rec_12 = max(1.D0/chi2_hard_12,EPS100)
+ hard_rec_21 = max(1.D0/chi2_hard_21,EPS100)
DO J=0,J0MAX
- hard_rec_11 = hard_rec_11*chi2_hard_11
- hard_rec_22 = hard_rec_22*chi2_hard_22
- hard_rec_12 = hard_rec_12*chi2_hard_12
- hard_rec_21 = hard_rec_21*chi2_hard_21
+ hard_rec_11 = max(hard_rec_11*chi2_hard_11,EPS100)
+ hard_rec_22 = max(hard_rec_22*chi2_hard_22,EPS100)
+ hard_rec_12 = max(hard_rec_12*chi2_hard_12,EPS100)
+ hard_rec_21 = max(hard_rec_21*chi2_hard_21,EPS100)
P_int(I,J) = P_int(I,J)
& + fe_11*soft_rec_11*hard_rec_11*fac_11
& + fe_22*soft_rec_22*hard_rec_22*fac_22
@@ -5535,7 +5544,7 @@ c read hadron-nucleon cross section from table
IF(IPAR(12).eq.3)THEN
c distinguish between nuclear cross sections..
IF(IAFLG.eq.0)THEN
-c if target is nucleus calc. hadron-nucleus cross section (slow)
+c if target is nucleus calc. hadron-nucleus cross section
CALL SIB_SIGMA_HNUC(L,IA,SQS,SIGprod,SIGbdif,SIGela)
ELSE
c if target is air read hadron-air cross section from table
@@ -11209,8 +11218,8 @@ C-----------------------------------------------------------------------
10 FORMAT(//,' Table of cross sections, and average number',
& ' of minijets and wounded nucleons ')
15 FORMAT(' [sqrt(s) in GeV, cross sections in mbarn]. ')
-16 FORMAT(' sqrt(s) sig(pp) sig(pAir) <n_s> <n_j> <n_w>',
- & ' sig(pip) sig(piAir) <n_s> <n_j> <n_w>')
+16 FORMAT(' sqrt(s) sig(pp) sig(pA) <n_s> <n_j> <n_w>',
+ & ' sig(pip) sig(piA) <n_s> <n_j> <n_w>')
18 FORMAT(1X,77('-') )
DO J=1,61,1
SQS = 10.D0**(ASQSMIN + DASQS*DBLE(J-1))
@@ -11251,7 +11260,7 @@ C-----------------------------------------------------------------------
C=======================================================================
- SUBROUTINE SIG_AIR_INI
+ SUBROUTINE SIG_AIR_INI
C-----------------------------------------------------------------------
C... Initialize the cross section and interaction lengths on air,
@@ -11375,7 +11384,6 @@ c interaction length
endif
ENDDO
END
-
C=======================================================================
SUBROUTINE SAMPLE_TARGET(NW,XCHG,KRMNT,XJET,Irec,IREJ)
@@ -11870,6 +11878,7 @@ C. INPUT : SSIG (mbarn) total pp cross section
C. SLOPE (GeV**-2) elastic scattering slope for pp
C. ALPHA real/imaginary part of the forward pp elastic
C. scattering amplitude
+C. ALAM coupling to inel. intermediat states
C. OUTPUT : SIGT = Total cross section
C. SIGEL = Elastic cross section
C. SIGQEL = Elastic + Quasi elastic cross section
@@ -11896,6 +11905,7 @@ C......................................................................
SIGQSD = (1.D0-FOX)*SIGQSD1 + FOX*SIGQSD2
RETURN
END
+
C=======================================================================
SUBROUTINE SIG_H_NUC
@@ -11946,7 +11956,7 @@ C......................................................................
END
C=======================================================================
-
+
SUBROUTINE GLAUBER2
+ (JA,SSIG,SLOPE,ALPHA,ALAM,SIGT,SIGEL,SIGQEL,SIGSD,SIGQSD)
@@ -14080,7 +14090,6 @@ c internal
END
C=======================================================================
-
SUBROUTINE SIB_SIGMA_HNUC (L,IAT,SQS,SIGprod,SIGbdif,SIGela)
C-----------------------------------------------------------------------
@@ -17831,8 +17840,3561 @@ 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)
C-----------------------------------------------------------------------
--
GitLab