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=======================================================================