*                                                                      *
*=== fluscw ===========================================================*
*                                                                      *
      DOUBLE PRECISION FUNCTION FLUSCW ( IJ    , PLA   , TXX   , TYY   ,
     &                                   TZZ   , WEE   , XX    , YY    ,
     &                                   ZZ    , NREG  , IOLREG, LLO   ,
     &                                   NSURF )

      INCLUDE 'dblprc.inc'
      INCLUDE 'dimpar.inc'
      INCLUDE 'iounit.inc'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 2003-2019:  CERN & INFN                            *
*     All Rights Reserved.                                             *
*                                                                      *
*     New version of Fluscw for FLUKA9x-FLUKA20xy:                     *
*                                                                      *
*     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     *
*     !!! This is a completely dummy routine for Fluka9x/200x. !!!     *
*     !!! The  name has been kept the same as for older  Fluka !!!     *
*     !!! versions for back-compatibility, even though  Fluscw !!!     *
*     !!! is applied only to estimators which didn't exist be- !!!     *
*     !!! fore Fluka89.                                        !!!     *
*     !!! User  developed versions  can be used for  weighting !!!     *
*     !!! flux-like quantities at runtime                      !!!     *
*     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     *
*                                                                      *
*     Last modified by D. Bozzato (CERN) and R. Froeschl (CERN)        *
*     on may 2023 for the FLUKA.CERN Advanced Course                   *
*                                                                      *
*     Input variables:                                                 *
*                                                                      *
*           Ij = (generalized) particle code (Paprop numbering)        *
*          Pla = particle laboratory momentum (GeV/c) (if > 0),        *
*                or kinetic energy (GeV) (if <0 )                      *
*    Txx,yy,zz = particle direction cosines                            *
*          Wee = particle weight                                       *
*     Xx,Yy,Zz = position                                              *
*         Nreg = (new) region number                                   *
*       Iolreg = (old) region number                                   *
*          Llo = particle generation                                   *
*        Nsurf = transport flag (ignore!)                              *
*                                                                      *
*     Output variables:                                                *
*                                                                      *
*       Fluscw = factor the scored amount will be multiplied by        *
*       Lsczer = logical flag, if true no amount will be scored        *
*                regardless of Fluscw                                  *
*                                                                      *
*     Useful variables (common SCOHLP):                                *
*                                                                      *
*     Flux like binnings/estimators (Fluscw):                          *
*          ISCRNG = 1 --> Boundary crossing estimator                  *
*          ISCRNG = 2 --> Track  length     binning                    *
*          ISCRNG = 3 --> Track  length     estimator                  *
*          ISCRNG = 4 --> Collision density estimator                  *
*          ISCRNG = 5 --> Yield             estimator                  *
*          JSCRNG = # of the binning/estimator                         *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE 'usrbin.inc'
      INCLUDE 'usrtrc.inc'
      INCLUDE 'scohlp.inc'
      INCLUDE 'trackr.inc'
      INCLUDE 'flkmat.inc'
      INCLUDE 'paprop.inc'

      INTEGER ILOW
      CHARACTER*10 SCORNAM

      INTEGER NENERGY
      PARAMETER (NENERGY = 46)
      DOUBLE PRECISION ENERGY(NENERGY), CROSSSEC(NENERGY)
      DOUBLE PRECISION EKIN, COEFF, NORM, HALFLIFE, LAMBDA
      DOUBLE PRECISION LOGINTRP
      SAVE ENERGY, CROSSSEC, NORM, HALFLIFE, LAMBDA

      LOGICAL LFIRST
      DATA LFIRST / .TRUE. /
      SAVE LFIRST


      ! Energy points in GeV
      DATA ENERGY /
     & 5.30900D-03, 5.95700D-03, 6.68300D-03, 7.49900D-03, 8.41400D-03,
     & 9.44100D-03, 1.05900D-02, 1.18900D-02, 1.33400D-02, 1.49700D-02,
     & 1.67900D-02, 1.88300D-02, 2.11300D-02, 2.37200D-02, 2.66100D-02,
     & 2.98500D-02, 3.34900D-02, 3.75800D-02, 4.21700D-02, 4.73200D-02,
     & 5.30900D-02, 5.95700D-02, 6.68300D-02, 7.49900D-02, 8.41400D-02,
     & 9.44100D-02, 1.05900D-01, 1.18900D-01, 1.33400D-01, 1.49700D-01,
     & 1.67900D-01, 1.88300D-01, 2.11300D-01, 2.37200D-01, 2.66100D-01,
     & 2.98500D-01, 3.34900D-01, 3.75800D-01, 4.21700D-01, 4.73200D-01,
     & 5.30900D-01, 5.95700D-01, 6.68300D-01, 7.49900D-01, 8.41400D-01,
     & 9.44100D-01 /


      ! Al-27(n,x)Na-24 production cross section in b
      DATA CROSSSEC /
     & 1.17700D-04, 1.55400D-03, 9.69100D-03, 2.74600D-02, 4.93700D-02,
     & 7.55200D-02, 9.81300D-02, 1.16300D-01, 1.25100D-01, 1.10200D-01,
     & 7.94000D-02, 5.20400D-02, 2.89900D-02, 1.54100D-02, 9.11600D-03,
     & 6.25200D-03, 4.25500D-03, 3.21600D-03, 3.80300D-03, 7.21300D-03,
     & 1.01500D-02, 1.40600D-02, 1.85200D-02, 2.07400D-02, 1.99900D-02,
     & 1.80100D-02, 1.52000D-02, 1.27700D-02, 1.10400D-02, 1.00000D-02,
     & 7.70200D-03, 6.27500D-03, 5.61300D-03, 5.15900D-03, 4.80500D-03,
     & 4.50800D-03, 4.24800D-03, 4.01400D-03, 3.79900D-03, 3.59900D-03,
     & 3.41100D-03, 3.23100D-03, 3.05900D-03, 2.89400D-03, 2.73300D-03,
     & 2.57700D-03 /
     
     
      ! Na-24 half-life (h)
      DATA HALFLIFE / 1.4997D+01 /
     


      !-----------------------------------------------------------------
      ! FIRST CALL
      
      IF ( LFIRST ) THEN
      
         WRITE( LUNOUT, * ) "#"
         WRITE( LUNOUT, * ) "# [fluscw] Requested fluence weighting ",
     &                      " on selected scorings:"
         WRITE( LUNOUT, * ) "# - USRBIN with SDUM(1:6) = 'Na24AS' -> ",
     &                      "specific activity (Bq/g/primary) of Na-24 ",
     &                      "[Al-27(n,x)Na-24] at saturation "
         WRITE( LUNOUT, * ) "# - USRBIN with SDUM(1:6) = 'Na24AC' -> ",
     &                      "specific activity (Bq/g/primary) of Na-24 ",
     &                      "[Al-27(n,x)Na-24] following 10 y ",
     &                      "irradiation and 12 h cooling time."
         WRITE( LUNOUT, * ) "#"

         ! Normalization factor for Al-27(n,x)Na-24 radionuclide
         ! production cross section (b) to Na-24 mass-specific
         ! activity (Bq/g) at saturation
         NORM = 1.0D-24 * AVOGAD / 2.6981539D+01
         
         ! Na-24 decay constant (1/s) to later compute normalizations for
         ! Na-24 mass-specific activity taking into account build-up
         ! and decay
         LAMBDA = LOG( TWOTWO ) / ( HALFLIFE * 3.6D+03 )

         LFIRST = .FALSE.

      ENDIF
      
      !-----------------------------------------------------------------



      !-----------------------------------------------------------------
      ! GENERAL PART

      ! Default is weighting is 1.00D+00
      FLUSCW = ONEONE
      LSCZER = .FALSE.

      ! Get current scoring SDUM for USRBIN
      IF ( ISCRNG .EQ. 2 ) THEN
         SCORNAM = TITUSB(JSCRNG)
      ! Ignore all other scoring options, no weight applied
      ELSE
         FLUSCW = ONEONE
         RETURN
      ENDIF

      ! Get the kinetic energy
      IF ( PLA .LT. ZERZER ) THEN
         EKIN = ABS(PLA)
      ELSEIF ( PLA .GT. ZERZER ) THEN
         EKIN = SQRT( PLA**2 + AM(IJ)**2 ) - AM(IJ)
         EKIN = EKIN
      ! Redundant check
      ELSE
         FLUSCW = ZERZER
         RETURN
      ENDIF
      
      !-----------------------------------------------------------------      
      


      !-----------------------------------------------------------------
      ! APPLY THE CORRECT FLUENCE WEIGHTING
      
      ! Check scoring name
      IF ( SCORNAM(1:5) .EQ. "Na24A" ) THEN
      
         !--------------------------------------------------------------         
         ! Select neutrons only
         
         IF ( IJ .EQ. 8 ) THEN
         
            !-----------------------------------------------------------      
            ! Get the cross section for that energy 
            
            ! Energy below minimum tabulated value
            IF ( EKIN .LT. ENERGY(1) ) THEN
               COEFF = ZERZER
            ! Energy above the maximum tabulated value
            ELSEIF ( EKIN .GE. ENERGY(NENERGY) ) THEN
               COEFF = CROSSSEC(NENERGY)
            ! All other cases: locate the closest energy
            ! with the helper subroutine SEARCHINDEX
            ELSE
               CALL LCT116( ENERGY, NENERGY, EKIN, ILOW )
               ! Particle's energy matches a tabulated one
               IF ( EKIN .EQ. ENERGY(ILOW) ) THEN
                  COEFF = CROSSSEC(ILOW)
               ! Perform a log-log interpolation between values otherwise
               ELSE
                  COEFF = LOGINTRP( ENERGY(ILOW),   ENERGY(ILOW+1),
     &                              CROSSSEC(ILOW), CROSSSEC(ILOW+1),
     &                              EKIN )
               ENDIF
            ENDIF

            !----------------------------------------------------------- 
               
            !-----------------------------------------------------------             
            ! Apply correct normalization
            
            ! Na-24 mass-specific activity (Bq/g) at saturation
            IF ( SCORNAM(1:6) .EQ. "Na24AS" ) THEN
               FLUSCW = COEFF * NORM              
            ! Na-24 mass-specific activity (Bq/g) following 10 y
            ! irradiation and 12 h cooling time
            ELSEIF ( SCORNAM(1:6) .EQ. "Na24AC" ) THEN
               FLUSCW = COEFF * NORM * 
     &                  (ONEONE - EXP(-LAMBDA * 3.15D+08)) *
     &                  EXP(-LAMBDA * (1.2D+01 * 3.6D+03))
            ! No normalization
            ELSE
               FLUSCW = COEFF
            ENDIF
            
            !----------------------------------------------------------- 

         ! If the particle is not a neutron (but scoring name matches), 
         ! apply zero weighting
         ELSE 
            FLUSCW = ZERZER
            RETURN
         ENDIF
         !--------------------------------------------------------------  
             
      ! Scoring name does not match, no weight applied
      ELSE
         FLUSCW = ONEONE
         RETURN
      ENDIF
      
      !-----------------------------------------------------------------
      
      RETURN
*=== End of function Fluscw ===========================================*
      END


*                                                                      *
*=== searchindex ======================================================*
*                                                                      *
      ! Helper subroutine
      ! Given an ordered array XX of size N and a value X
      ! (with X >= XX(1) and  X < XX(N)) locate the index J
      ! such that XX(J) <= X < XX(J+1)


c       SUBROUTINE SEARCHINDEX( XX, N, X, J )
c
c       INTEGER J, N
c       DOUBLE PRECISION X, XX(N)
c       INTEGER JL, JM, JU
c
c       JL = 0
c       JU = N + 1
c   10  IF (JU - JL .GT. 1) THEN
c         JM = (JU + JL)/2
c         IF ((XX(N) .GE. XX(1)) .EQV. (X .GE. XX(JM))) THEN
c          JL = JM
c         ELSE
c          JU = JM
c         ENDIF
c         GOTO 10
c       ENDIF
c       IF(X .EQ. XX(1))THEN
c         J = 1
c       ELSE IF(X .EQ. XX(N))THEN
c         J = N - 1
c       ELSE
c         J = JL
c       ENDIF
c
c       RETURN
c       END
c


*                                                                      *
*=== loglog ===========================================================*
*                                                                      *
      ! Helper function for log10-log10 interpolation
      DOUBLE PRECISION FUNCTION LOGINTRP( XMIN, XMAX, YMIN, YMAX, XVAL )

      DOUBLE PRECISION XMIN, XMAX, YMIN, YMAX, XVAL
      DOUBLE PRECISION LXMIN, LXMAX, LYMIN, LYMAX, LXVAL, LINTERP

      LXMIN = LOG10( XMIN )
      LXMAX = LOG10( XMAX )
      LYMIN = LOG10( YMIN )
      LYMAX = LOG10( YMAX )
      LXVAL = LOG10( XVAL )
      LINTERP = LYMIN + (LXVAL - LXMIN)*(LYMAX - LYMIN)/(LXMAX - LXMIN)
      LOGINTRP = 10.0D+00 ** LINTERP

      RETURN
      END

