IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 078fb175 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan: Committed by Felix Riehn
Browse files

Resolve "Make CORSIKA 8 fast again"

parent d5a40b5c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Processes/Sibyll/sibyll2.3c.f
+ 195
− 100
View file @ 078fb175
...@@ -425,7 +425,7 @@ C. important information for the generation of events ...@@ -425,7 +425,7 @@ C. important information for the generation of events
C. C.
C PARAMETER (NS_max = 20, NH_max = 80) C PARAMETER (NS_max = 20, NH_max = 80)
C COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2), C COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
C & SSIGN(61,3),SSIGNSD(61,3) ALINT(61,3), ASQSMIN, ASQSMAX, DASQS, NSQS C & SSIGN(61,3,3),SSIGNSD(61,3,3) ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
C. C.
C. NSQS = number of energy points (61 is current version) C. NSQS = number of energy points (61 is current version)
C. ASQSMIN = log_10 [sqrt(s) GeV] minimum value C. ASQSMIN = log_10 [sqrt(s) GeV] minimum value
...@@ -436,8 +436,10 @@ C. ...@@ -436,8 +436,10 @@ C.
C. SSIG(J,1) inelastic cross section for pp interaction C. SSIG(J,1) inelastic cross section for pp interaction
C. at energy: sqrt(s)(GeV) = 10**[ASQSMIN+DASQS*(J-1)] C. at energy: sqrt(s)(GeV) = 10**[ASQSMIN+DASQS*(J-1)]
C. SSIG(J,2) inelastic cross section for pi-p interaction C. SSIG(J,2) inelastic cross section for pi-p interaction
C. SSIGN(J,1) inelastic cross section for p-Air interaction C. SSIGN(J,1,1) inelastic cross section for p-Air interaction
C. SSIGN(J,2) inelastic cross section for pi-Air interaction C. SSIGN(J,2,1) inelastic cross section for pi-Air interaction
C. SSIGN(J,1,2) inelastic cross section for p-Nitrogen interaction
C. SSIGN(J,2,2) inelastic cross section for pi-Nitrogen interaction
C. C.
C. PJETC(n_s,n_j,J,1) Cumulative probability distribution C. PJETC(n_s,n_j,J,1) Cumulative probability distribution
C. for the production of n_s soft interactions and C. for the production of n_s soft interactions and
...@@ -10886,12 +10888,12 @@ c COMMON /S_DEBUG/ Ncall, Ndebug, Lun ...@@ -10886,12 +10888,12 @@ c COMMON /S_DEBUG/ Ncall, Ndebug, Lun
INTEGER NS_max, NH_max INTEGER NS_max, NH_max
PARAMETER (NS_max = 20, NH_max = 80) PARAMETER (NS_max = 20, NH_max = 80)
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,ALINT,ASQSMIN,ASQSMAX, DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,SSIGNEL,ALINT,ASQSMIN,
& DASQS & ASQSMAX,DASQS
INTEGER NSQS INTEGER NSQS
COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2), COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
& SSIGN(61,3), SSIGNSD(61,3), ALINT(61,3), & SSIGN(61,3,3), SSIGNSD(61,3,3), SSIGNEL(61,3,3),
& ASQSMIN, ASQSMAX, DASQS, NSQS & ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
DOUBLE PRECISION STR_mass_val, STR_mass_val_hyp, STR_mass_sea DOUBLE PRECISION STR_mass_val, STR_mass_val_hyp, STR_mass_sea
COMMON /S_CUTOFF/ STR_mass_val, STR_mass_val_hyp, STR_mass_sea COMMON /S_CUTOFF/ STR_mass_val, STR_mass_val_hyp, STR_mass_sea
   
...@@ -10993,12 +10995,12 @@ C----------------------------------------------------------------------- ...@@ -10993,12 +10995,12 @@ C-----------------------------------------------------------------------
INTEGER NS_max, NH_max INTEGER NS_max, NH_max
PARAMETER (NS_max = 20, NH_max = 80) PARAMETER (NS_max = 20, NH_max = 80)
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,ALINT,ASQSMIN,ASQSMAX, DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,SSIGNEL,ALINT,ASQSMIN,
& DASQS & ASQSMAX,DASQS
INTEGER NSQS INTEGER NSQS
COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2), COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
& SSIGN(61,3), SSIGNSD(61,3), ALINT(61,3), & SSIGN(61,3,3), SSIGNSD(61,3,3), SSIGNEL(61,3,3),
& ASQSMIN, ASQSMAX, DASQS, NSQS & ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
DOUBLE PRECISION SSIG_TOT,SSIG_SD1,SSIG_SD2,SSIG_DD,SSIG_B, DOUBLE PRECISION SSIG_TOT,SSIG_SD1,SSIG_SD2,SSIG_DD,SSIG_B,
& SSIG_RHO & SSIG_RHO
COMMON /S_CCSIG2/ SSIG_TOT(61,3),SSIG_SD1(61,3),SSIG_SD2(61,3), COMMON /S_CCSIG2/ SSIG_TOT(61,3),SSIG_SD1(61,3),SSIG_SD2(61,3),
...@@ -11187,12 +11189,12 @@ C----------------------------------------------------------------------- ...@@ -11187,12 +11189,12 @@ C-----------------------------------------------------------------------
INTEGER NS_max, NH_max INTEGER NS_max, NH_max
PARAMETER (NS_max = 20, NH_max = 80) PARAMETER (NS_max = 20, NH_max = 80)
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,ALINT,ASQSMIN,ASQSMAX, DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,SSIGNEL,ALINT,ASQSMIN,
& DASQS & ASQSMAX,DASQS
INTEGER NSQS INTEGER NSQS
COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2), COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
& SSIGN(61,3), SSIGNSD(61,3), ALINT(61,3), & SSIGN(61,3,3), SSIGNSD(61,3,3), SSIGNEL(61,3,3),
& ASQSMIN, ASQSMAX, DASQS, NSQS & ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
DIMENSION PJ(2),PS(2),PW(2) DIMENSION PJ(2),PS(2),PW(2)
   
SAVE SAVE
...@@ -11207,15 +11209,15 @@ C----------------------------------------------------------------------- ...@@ -11207,15 +11209,15 @@ C-----------------------------------------------------------------------
10 FORMAT(//,' Table of cross sections, and average number', 10 FORMAT(//,' Table of cross sections, and average number',
& ' of minijets and wounded nucleons ') & ' of minijets and wounded nucleons ')
15 FORMAT(' [sqrt(s) in GeV, cross sections in mbarn]. ') 15 FORMAT(' [sqrt(s) in GeV, cross sections in mbarn]. ')
16 FORMAT(' sqrt(s) sig(pp) sig(pA) <n_s> <n_j> <n_w>', 16 FORMAT(' sqrt(s) sig(pp) sig(pAir) <n_s> <n_j> <n_w>',
& ' sig(pip) sig(piA) <n_s> <n_j> <n_w>') & ' sig(pip) sig(piAir) <n_s> <n_j> <n_w>')
18 FORMAT(1X,77('-') ) 18 FORMAT(1X,77('-') )
DO J=1,61,1 DO J=1,61,1
SQS = 10.D0**(ASQSMIN + DASQS*DBLE(J-1)) SQS = 10.D0**(ASQSMIN + DASQS*DBLE(J-1))
   
DO K=1,2 DO K=1,2
   
PW(K) = ATARGET*SSIG(J,K)/SSIGN(J,K) PW(K) = ATARGET*SSIG(J,K)/SSIGN(J,K,1)
   
PJ(K) = 0.D0 PJ(K) = 0.D0
PS(K) = 0.D0 PS(K) = 0.D0
...@@ -11235,8 +11237,8 @@ C----------------------------------------------------------------------- ...@@ -11235,8 +11237,8 @@ C-----------------------------------------------------------------------
   
ENDDO ENDDO
   
WRITE(LUN,20) SQS,SSIG(J,1),SSIGN(J,1),PS(1),PJ(1),PW(1) WRITE(LUN,20) SQS,SSIG(J,1),SSIGN(J,1,1),PS(1),PJ(1),PW(1)
& ,SSIG(J,2),SSIGN(J,2),PS(2),PJ(2),PW(2) & ,SSIG(J,2),SSIGN(J,2,1),PS(2),PJ(2),PW(2)
   
ENDDO ENDDO
   
...@@ -11252,7 +11254,8 @@ C======================================================================= ...@@ -11252,7 +11254,8 @@ C=======================================================================
SUBROUTINE SIG_AIR_INI SUBROUTINE SIG_AIR_INI
   
C----------------------------------------------------------------------- C-----------------------------------------------------------------------
C...Initialize the cross section and interaction lengths on air C... Initialize the cross section and interaction lengths on air,
C. nitrogen and oxygen
C. (this version initializes p-air, pi-air, and K-air cross sections) C. (this version initializes p-air, pi-air, and K-air cross sections)
C. C.
C. also calculates the low mass beam diffraction cross section in hAir \FR C. also calculates the low mass beam diffraction cross section in hAir \FR
...@@ -11264,12 +11267,12 @@ C----------------------------------------------------------------------- ...@@ -11264,12 +11267,12 @@ C-----------------------------------------------------------------------
INTEGER NS_max, NH_max INTEGER NS_max, NH_max
PARAMETER (NS_max = 20, NH_max = 80) PARAMETER (NS_max = 20, NH_max = 80)
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,ALINT,ASQSMIN,ASQSMAX, DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,SSIGNEL,ALINT,ASQSMIN,
& DASQS & ASQSMAX,DASQS
INTEGER NSQS INTEGER NSQS
COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2), COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
& SSIGN(61,3), SSIGNSD(61,3), ALINT(61,3), & SSIGN(61,3,3), SSIGNSD(61,3,3), SSIGNEL(61,3,3),
& ASQSMIN, ASQSMAX, DASQS, NSQS & ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
COMMON /GLAUB_SCR/ XI_MAX , ALAM(61) COMMON /GLAUB_SCR/ XI_MAX , ALAM(61)
   
INTEGER NCALL, NDEBUG, LUN INTEGER NCALL, NDEBUG, LUN
...@@ -11280,15 +11283,15 @@ C----------------------------------------------------------------------- ...@@ -11280,15 +11283,15 @@ C-----------------------------------------------------------------------
INTEGER IPAR INTEGER IPAR
COMMON /S_CFLAFR/ PAR(NPAR_max), IPAR(NIPAR_max) COMMON /S_CFLAFR/ PAR(NPAR_max), IPAR(NIPAR_max)
DIMENSION SIGDIF(3) DIMENSION SIGDIF(3)
DIMENSION ITARGC(3)
CHARACTER*3 TARGN
DIMENSION TARGN(3)
SAVE SAVE
DATA AVOG /6.0221367D-04/ DATA AVOG /6.0221367D-04/
DATA ATARGET /14.514D0/ DATA ATARGET /14.514D0/
DATA ITARGC /0,14,16/
DATA TARGN /'air','nit','oxy'/
   
c PRINT *,'inel. screening in hadron - nucleus interactions'
c ALAM = 0.5
c PRINT *,'const. coupling: ', ALAM
IF ( IPAR(12).GT.0 ) THEN IF ( IPAR(12).GT.0 ) THEN
if (ndebug.gt.0) then if (ndebug.gt.0) then
WRITE(LUN,*) ' SIG_AIR_INI:' WRITE(LUN,*) ' SIG_AIR_INI:'
...@@ -11298,61 +11301,81 @@ c PRINT *,'const. coupling: ', ALAM ...@@ -11298,61 +11301,81 @@ c PRINT *,'const. coupling: ', ALAM
if (ndebug.gt.0)WRITE(LUN,*)' low mass Xi_max: ' , XI_MAX if (ndebug.gt.0)WRITE(LUN,*)' low mass Xi_max: ' , XI_MAX
ENDIF ENDIF
   
C...target loop (air, N, O)
DO IK=1,3
IAT = ITARGC(IK)
WRITE(6,*) 'SIG_AIR_INI: initializing target: (i,A)',
& ik, IAT, TARGN(IK) , '..'
C...particle loop (p, pi, K) C...particle loop (p, pi, K)
DO K=1,3 DO K=1,3
if (NDEBUG .gt. 0 ) then if (NDEBUG .gt. 0 ) then
WRITE(LUN,'(/,1X,A,A)') WRITE(6,'(/,1X,A,A)')
& 'Table: J, sqs, SIGtot, SIGprod, SIG_SD,', & 'Table: J, IK, sqs, SIGtot, SIGprod, SIG_SD,',
& ' Lambda ' &' Lambda '
WRITE(LUN,*) WRITE(6,*)
& '-------------------------------------------------', & '-------------------------------------------------',
& '-------------' & '-------------'
endif endif
DO J=1,NSQS DO J=1,NSQS
   
ASQS = ASQSMIN + DASQS*DBLE(J-1) ASQS = ASQSMIN + DASQS*DBLE(J-1)
SQS = 10.D0**ASQS SQS = 10.D0**ASQS
   
IF (K.EQ.1) THEN IF (K.EQ.1) THEN
c Goulianos param. from GAP-2012-056, Mx**2s = 0.02 c Goulianos param. from GAP-2012-056, Mx**2s = 0.02
c against PDG elastic cross section c against PDG elastic cross section
CALL SIB_HADCS1(K,SQS,SIGT1,SIGEL1,SIGINEL1,SLOPE1,RHO1) CALL SIB_HADCS1
SIGEFF = 0.68D0*(1.D0+36.D0/SQS**2)* & (K,SQS,SIGT1,SIGEL1,SIGINEL1,SLOPE1,RHO1)
& dlog(0.6D0+XI_MAX/1.5D0*SQS**2) SIGEFF = 0.68D0*(1.D0+36.D0/SQS**2)*
ALAM(J) = dSQRT(SIGEFF/SIGEL1) & dlog(0.6D0+XI_MAX/1.5D0*SQS**2)
ENDIF ALAM(J) = dSQRT(SIGEFF/SIGEL1)
ENDIF
c CALL SIB_HADCSL(k,SQS, CALL SIB_SIGMA_HP(K,SQS,
c & SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO) & SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
IF(IK.eq.1)THEN
CALL SIB_SIGMA_HP(K,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO) c fixed O-N mixture
CALL SIG_H_AIR CALL SIG_H_AIR
& (SIGT, SLOPE, RHO, ALAM(J), & (SIGT, SLOPE, RHO, ALAM(J),
& SSIGT, SSIGEL, SSIGQE, SSIGSD, SSIGQSD) & SSIGT, SSIGEL, SSIGQE, SSIGSD, SSIGQSD)
ELSE
if (ndebug .gt. 0 ) WRITE(LUN,'(1X,I2,1P,5E12.3)') CALL SIG_H_NUC
& K,SQS,SSIGT,SSIGT-SSIGQE,SSIGQSD,ALAM(J) & (IAT, SIGT, SLOPE, RHO, ALAM(J),
C particle production cross section & SSIGT, SSIGEL, SSIGQE, SSIGSD, SSIGQSD)
SSIGN(J,K) = SSIGT-SSIGQE ENDIF
SSIGNSD(J,K) = SSIGQSD if (ndebug .gt. 0 ) WRITE(6,'(1X,I2,1P,5E12.3)')
ALINT(J,K) = 1.D0/(AVOG*SSIGn(j,K)/ATARGET) & K,SQS,SSIGT,SSIGT-SSIGQE,SSIGQSD,ALAM(J)
ENDDO C particle production cross section
SSIGN(J,K,IK) = SSIGT-SSIGQE
c diffractive cross section
SSIGNSD(J,K,IK) = SSIGQSD
c elastic cross section
SSIGNEL(J,K,IK) = SSIGEL
c interaction length
IF(IK.eq.1)then
ALINT(J,K,IK) = 1.D0/(AVOG*SSIGn(j,K,IK)/ATARGET)
else
ALINT(J,K,IK) = 1.D0/(AVOG*SSIGn(j,K,IK)/IAT)
endif
ENDDO
ENDDO
if (ndebug .gt. 0 ) then
WRITE(6,'(/,1X,A)')
& ' SIG_AIR_INI: NUCLIB interaction lengths [g/cm**2]'
WRITE(6,*) 'target:', TARGN(IK)
WRITE(6,'(1X,A)')
& ' sqs, p-targ, pi-targ, K-targ'
DO J=1,NSQS
ASQS = ASQSMIN + DASQS*DBLE(J-1)
SQS = 10.D0**ASQS
WRITE(6,'(1X,1P,4E12.3)')
& SQS,ALINT(J,1,IK),ALINT(J,2,IK),ALINT(J,3,IK)
ENDDO
endif
ENDDO ENDDO
if (ndebug .gt. 0 ) then
WRITE(LUN,'(/,1X,A)')
& ' SIG_AIR_INI: NUCLIB interaction lengths [g/cm**2]'
WRITE(LUN,'(1X,A)')
& ' sqs, p-air, pi-air, K-air'
DO J=1,NSQS
ASQS = ASQSMIN + DASQS*DBLE(J-1)
SQS = 10.D0**ASQS
WRITE(LUN,'(1X,1P,4E12.3)')
& SQS,ALINT(J,1),ALINT(J,2),ALINT(J,3)
ENDDO
endif
END END
C======================================================================= C=======================================================================
   
SUBROUTINE SAMPLE_TARGET(NW,XCHG,KRMNT,XJET,Irec,IREJ) SUBROUTINE SAMPLE_TARGET(NW,XCHG,KRMNT,XJET,Irec,IREJ)
...@@ -11873,6 +11896,54 @@ C...................................................................... ...@@ -11873,6 +11896,54 @@ C......................................................................
SIGQSD = (1.D0-FOX)*SIGQSD1 + FOX*SIGQSD2 SIGQSD = (1.D0-FOX)*SIGQSD1 + FOX*SIGQSD2
RETURN RETURN
END END
C=======================================================================
SUBROUTINE SIG_H_NUC
+ (IAT,SSIG,SLOPE,ALPHA,ALAM,SIGT,SIGEL,SIGQE,SIGSD,SIGQSD)
C-----------------------------------------------------------------------
C**********************************************************************
C...Subroutine to compute hadron-nucleus cross sections
C. according to:
C. R.J. Glauber and G.Matthiae Nucl.Phys. B21, 135, (1970)
C.
C. INPUT : IAT nucleon number in target nucleus
C. SSIG (mbarn) total pp cross section
C. SLOPE (GeV**-2) elastic scattering slope for pp
C. ALPHA real/imaginary part of the forward pp elastic
C. scattering amplitude
C. OUTPUT : SIGT = Total cross section
C. SIGEL = Elastic cross section
C. SIGQEL = Elastic + Quasi elastic cross section
C. SIGSD = single diff. cross section (beam)
C. SIGQSD = Elastic + Quasi elastic SD cross section (beam)
C.
C......................................................................
IMPLICIT NONE
INTEGER IAT
DOUBLE PRECISION SSIG,SLOPE,ALPHA,ALAM
DOUBLE PRECISION SIG1,SIGEL1,SIGQE1,SIGSD1,SIGQSD1
DOUBLE PRECISION SIGT,SIGEL,SIGQE,SIGSD,SIGQSD
INTEGER NCALL, NDEBUG, LUN
COMMON /S_DEBUG/ NCALL, NDEBUG, LUN
SAVE
IF(IAT.eq.0.or.IAT.gt.18) THEN
WRITE(LUN,'(//,1X,A)')
& ' SIG_H_NUC: number of target nucleons too large!',
& ' (1<=IAT<=18)'
SIGT = -1.D0
STOP
ENDIF
CALL GLAUBER2
+ (IAT,SSIG,SLOPE,ALPHA,ALAM,SIG1,SIGEL1,SIGQE1,SIGSD1,SIGQSD1)
SIGT = SIG1
SIGEL = SIGEL1
SIGQE = SIGQE1
SIGSD = SIGSD1
SIGQSD = SIGQSD1
RETURN
END
   
C======================================================================= C=======================================================================
   
...@@ -13785,12 +13856,12 @@ c external types ...@@ -13785,12 +13856,12 @@ c external types
INTEGER NS_max, NH_max INTEGER NS_max, NH_max
PARAMETER (NS_max = 20, NH_max = 80) PARAMETER (NS_max = 20, NH_max = 80)
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,ALINT,ASQSMIN,ASQSMAX, DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,SSIGNEL,ALINT,ASQSMIN,
& DASQS & ASQSMAX,DASQS
INTEGER NSQS INTEGER NSQS
COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2), COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
& SSIGN(61,3), SSIGNSD(61,3), ALINT(61,3), & SSIGN(61,3,3), SSIGNSD(61,3,3), SSIGNEL(61,3,3),
& ASQSMIN, ASQSMAX, DASQS, NSQS & ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
DOUBLE PRECISION SSIG_TOT,SSIG_SD1,SSIG_SD2,SSIG_DD,SSIG_B, DOUBLE PRECISION SSIG_TOT,SSIG_SD1,SSIG_SD2,SSIG_DD,SSIG_B,
& SSIG_RHO & SSIG_RHO
COMMON /S_CCSIG2/ SSIG_TOT(61,3),SSIG_SD1(61,3),SSIG_SD2(61,3), COMMON /S_CCSIG2/ SSIG_TOT(61,3),SSIG_SD1(61,3),SSIG_SD2(61,3),
...@@ -13892,12 +13963,12 @@ c external types ...@@ -13892,12 +13963,12 @@ c external types
INTEGER NS_max, NH_max INTEGER NS_max, NH_max
PARAMETER (NS_max = 20, NH_max = 80) PARAMETER (NS_max = 20, NH_max = 80)
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,ALINT,ASQSMIN,ASQSMAX, DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,SSIGNEL,ALINT,ASQSMIN,
& DASQS & ASQSMAX,DASQS
INTEGER NSQS INTEGER NSQS
COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2), COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
& SSIGN(61,3), SSIGNSD(61,3), ALINT(61,3), & SSIGN(61,3,3), SSIGNSD(61,3,3), SSIGNEL(61,3,3),
& ASQSMIN, ASQSMAX, DASQS, NSQS & ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
DOUBLE PRECISION SSIG_TOT,SSIG_SD1,SSIG_SD2,SSIG_DD,SSIG_B, DOUBLE PRECISION SSIG_TOT,SSIG_SD1,SSIG_SD2,SSIG_DD,SSIG_B,
& SSIG_RHO & SSIG_RHO
COMMON /S_CCSIG2/ SSIG_TOT(61,3),SSIG_SD1(61,3),SSIG_SD2(61,3), COMMON /S_CCSIG2/ SSIG_TOT(61,3),SSIG_SD1(61,3),SSIG_SD2(61,3),
...@@ -13970,12 +14041,12 @@ Cf2py double precision, intent(out) :: SIGprod,SIGbdif ...@@ -13970,12 +14041,12 @@ Cf2py double precision, intent(out) :: SIGprod,SIGbdif
INTEGER NS_max, NH_max INTEGER NS_max, NH_max
PARAMETER (NS_max = 20, NH_max = 80) PARAMETER (NS_max = 20, NH_max = 80)
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,ALINT,ASQSMIN,ASQSMAX, DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,SSIGNEL,ALINT,ASQSMIN,
& DASQS & ASQSMAX,DASQS
INTEGER NSQS INTEGER NSQS
COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2), COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
& SSIGN(61,3), SSIGNSD(61,3), ALINT(61,3), & SSIGN(61,3,3), SSIGNSD(61,3,3), SSIGNEL(61,3,3),
& ASQSMIN, ASQSMAX, DASQS, NSQS & ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
   
c external c external
DOUBLE PRECISION SQS,SIGPROD,SIGBDIF DOUBLE PRECISION SQS,SIGPROD,SIGBDIF
...@@ -14004,8 +14075,8 @@ c internal ...@@ -14004,8 +14075,8 @@ c internal
J1 = max(J1,1) J1 = max(J1,1)
endif endif
T = (AL-1.D0)*10.D0 - DBLE(J1-1) T = (AL-1.D0)*10.D0 - DBLE(J1-1)
SIGprod = SSIGN(J1,L)*(1.D0-T) + SSIGN(J1+1,L)*T SIGprod = SSIGN(J1,L,1)*(1.D0-T) + SSIGN(J1+1,L,1)*T
SIGbdif = SSIGNSD(J1,L)*(1.D0-T) + SSIGNSD(J1+1,L)*T SIGbdif = SSIGNSD(J1,L,1)*(1.D0-T) + SSIGNSD(J1+1,L,1)*T
   
END END
C======================================================================= C=======================================================================
...@@ -14024,6 +14095,7 @@ C SQS sqrt(s) ...@@ -14024,6 +14095,7 @@ C SQS sqrt(s)
C C
C output: SIGprod particle production cross section (mb) C output: SIGprod particle production cross section (mb)
C SIGbdif q.ela and ela beam diff. cross section C SIGbdif q.ela and ela beam diff. cross section
C SIGela elastic cross section
C----------------------------------------------------------------------- C-----------------------------------------------------------------------
Cf2py integer, intent(in) :: L,IAT Cf2py integer, intent(in) :: L,IAT
Cf2py double precision, intent(in) :: SQS Cf2py double precision, intent(in) :: SQS
...@@ -14035,12 +14107,12 @@ Cf2py double precision, intent(out) :: SIGprod,SIGbdif ...@@ -14035,12 +14107,12 @@ Cf2py double precision, intent(out) :: SIGprod,SIGbdif
INTEGER NCALL, NDEBUG, LUN INTEGER NCALL, NDEBUG, LUN
COMMON /S_DEBUG/ NCALL, NDEBUG, LUN COMMON /S_DEBUG/ NCALL, NDEBUG, LUN
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,ALINT,ASQSMIN,ASQSMAX, DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,SSIGNEL,ALINT,ASQSMIN,
& DASQS & ASQSMAX,DASQS
INTEGER NSQS INTEGER NSQS
COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2), COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
& SSIGN(61,3), SSIGNSD(61,3), ALINT(61,3), & SSIGN(61,3,3), SSIGNSD(61,3,3), SSIGNEL(61,3,3),
& ASQSMIN, ASQSMAX, DASQS, NSQS & ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
   
C-------------------------------------------------------------------- C--------------------------------------------------------------------
C SIBYLL utility common blocks containing constants \FR'14 C SIBYLL utility common blocks containing constants \FR'14
...@@ -14065,17 +14137,40 @@ c external ...@@ -14065,17 +14137,40 @@ c external
INTEGER L,IAT INTEGER L,IAT
   
c internal c internal
DOUBLE PRECISION ALAM DOUBLE PRECISION ALAM,T,AL
INTEGER IPARM,ICSMOD INTEGER IPARM,ICSMOD,IK,J1
SAVE SAVE
IF(NSQS.LE.0) THEN IF(NSQS.LE.0) THEN
WRITE(LUN,'(//,1X,A)') WRITE(LUN,'(//,1X,A)')
& ' SIB_SIGMA_HNUC: interpolation table not initialized.' & ' SIB_SIGMA_HNUC: interpolation table not initialized.'
STOP STOP
ENDIF ENDIF
IF(IAT.eq.0.or.IAT.eq.14.or.IAT.eq.16)THEN
IF(IAT.ge.0.and.IAT.lt.19)THEN c read cross section from table
IF(IAT.eq.0) THEN
IK=1
ELSEIF(IAT.eq.14)THEN
IK=2
ELSE
IK=3
ENDIF
AL = LOG10(SQS)
J1 = INT((AL - 1.D0)*10.D0 + 1)
if((j1.lt.1).or.(j1.gt.NSQS)) then
if (ndebug .gt. 0)
& write (LUN,'(1x,a,i3,1p,e12.3)')
& ' SIB_SIGMA_HNUC: energy out of range ',L,sqs
endif
if((j1.lt.1).or.(j1.ge.NSQS)) then
J1 = min(J1,NSQS-1)
J1 = max(J1,1)
endif
T = (AL-1.D0)*10.D0 - DBLE(J1-1)
SIGprod = SSIGN(J1,L,IK)*(1.D0-T) + SSIGN(J1+1,L,IK)*T
SIGbdif = SSIGNSD(J1,L,IK)*(1.D0-T) + SSIGNSD(J1+1,L,IK)*T
SIGela = SSIGNEL(J1,L,IK)*(1.D0-T) + SSIGNEL(J1+1,L,IK)*T
ELSEIF(IAT.lt.19)THEN
c calculate cross section
IF(ndebug.gt.0)THEN IF(ndebug.gt.0)THEN
WRITE(LUN,'(1X,A,2I4,F8.2)') WRITE(LUN,'(1X,A,2I4,F8.2)')
& 'SIB_SIGMA_HNUC: L,IAT,SQS:',L,IAT,SQS & 'SIB_SIGMA_HNUC: L,IAT,SQS:',L,IAT,SQS
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum