Line data    Source code

       1             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       2             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       3             : !!!!                                                                                                                            !!!!
       4             : !!!!    ParaMonte: Parallel Monte Carlo and Machine Learning Library.                                                           !!!!
       5             : !!!!                                                                                                                            !!!!
       6             : !!!!    Copyright (C) 2012-present, The Computational Data Science Lab                                                          !!!!
       7             : !!!!                                                                                                                            !!!!
       8             : !!!!    This file is part of the ParaMonte library.                                                                             !!!!
       9             : !!!!                                                                                                                            !!!!
      10             : !!!!    LICENSE                                                                                                                 !!!!
      11             : !!!!                                                                                                                            !!!!
      12             : !!!!       https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md                                                          !!!!
      13             : !!!!                                                                                                                            !!!!
      14             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      15             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      16             : 
      17             : !>  \brief
      18             : !>  This file contains implementations of procedures [pm_cosmicRate](@ref pm_cosmicRate).
      19             : !>
      20             : !>  \finmain
      21             : !>
      22             : !>  \author
      23             : !>  \AmirShahmoradi, September 1, 2017, 12:00 AM, Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin
      24             : 
      25             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      26             : 
      27             :         ! Define the array indexing rules.
      28             : #if     D1_ENABLED
      29             :         integer(IK) :: i
      30             : #define GET_ELEMENT(Array) Array(i)
      31             : #elif   D0_ENABLED
      32             : #define GET_ELEMENT(Array) Array
      33             : #define ALL
      34             : #else
      35             : #error  "Unrecognized interface."
      36             : #endif
      37             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      38             : #if     getLogRateDensity_ENABLED && (H06_ENABLED || L08_ENABLED || B10_ENABLED || P15_ENABLED)
      39             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      40             : 
      41             : #if     H06_ENABLED
      42             :         real(RKC), parameter :: G0 = +3.4_RKC
      43             :         real(RKC), parameter :: G1 = -0.3_RKC
      44             :         real(RKC), parameter :: G2 = -7.8_RKC
      45             :         real(RKC), parameter :: LOGZ0PLUS1 = log(1._RKC + 0.97_RKC)
      46             :         real(RKC), parameter :: LOGZ1PLUS1 = log(1._RKC + 4.50_RKC)
      47             :         real(RKC), parameter :: LOG_NORM_FAC_1 = LOGZ0PLUS1 * (G0 - G1)
      48             :         real(RKC), parameter :: LOG_NORM_FAC_2 = LOGZ1PLUS1 * (G1 - G2) + LOG_NORM_FAC_1
      49             : #elif   L08_ENABLED
      50             :         real(RKC), parameter :: G0 = +3.3000_RKC
      51             :         real(RKC), parameter :: G1 = +0.0549_RKC
      52             :         real(RKC), parameter :: G2 = -4.4600_RKC
      53             :         real(RKC), parameter :: LOGZ0PLUS1 = log(1._RKC + 0.993_RKC)
      54             :         real(RKC), parameter :: LOGZ1PLUS1 = log(1._RKC + 3.800_RKC)
      55             :         real(RKC), parameter :: LOG_NORM_FAC_1 = LOGZ0PLUS1 * (G0 - G1)
      56             :         real(RKC), parameter :: LOG_NORM_FAC_2 = LOGZ1PLUS1 * (G1 - G2) + LOG_NORM_FAC_1
      57             : #elif   B10_ENABLED
      58             :         real(RKC), parameter :: G0 = +3.14_RKC
      59             :         real(RKC), parameter :: G1 = +1.36_RKC
      60             :         real(RKC), parameter :: G2 = -2.92_RKC
      61             :         real(RKC), parameter :: LOGZ0PLUS1 = log(1._RKC + 0.97_RKC)
      62             :         real(RKC), parameter :: LOGZ1PLUS1 = log(1._RKC + 4.00_RKC)
      63             :         real(RKC), parameter :: LOG_NORM_FAC_1 = LOGZ0PLUS1 * (G0 - G1)
      64             :         real(RKC), parameter :: LOG_NORM_FAC_2 = LOGZ1PLUS1 * (G1 - G2) + LOG_NORM_FAC_1
      65             : #elif   P15_ENABLED
      66             :         real(RKC), parameter :: EXPONENT_HIGH_Z = -7.8_RKC
      67             :         real(RKC), parameter :: LOGZ1PLUS1 = log(1._RKC + 4.5_RKC)
      68             :         real(RKC), parameter :: LOG_NORM_FAC_2 = -EXPONENT_HIGH_Z * LOGZ1PLUS1
      69             : #else
      70             : #error  "Unrecognized interface."
      71             : #endif
      72       45455 :         CHECK_ASSERTION(__LINE__, all([logzplus1 >= 0._RKC]), SK_": The condition `logzplus1 >= 0.` must hold. logzplus1 = "//getStr(logzplus1)) ! fpp
      73             : #if     D1_ENABLED
      74             :         do concurrent(i = 1_IK : size(logzplus1, kind = IK))
      75             : #endif
      76             : #if         H06_ENABLED || L08_ENABLED || B10_ENABLED
      77       34983 :             if (GET_ELEMENT(logzplus1) < LOGZ0PLUS1) then
      78        3475 :                 GET_ELEMENT(logRateDensity) = GET_ELEMENT(logzplus1) * G0
      79       31508 :             elseif (GET_ELEMENT(logzplus1) < LOGZ1PLUS1) then
      80       11651 :                 GET_ELEMENT(logRateDensity) = GET_ELEMENT(logzplus1) * G1 + LOG_NORM_FAC_1
      81             :             else
      82       19857 :                 GET_ELEMENT(logRateDensity) = GET_ELEMENT(logzplus1) * G2 + LOG_NORM_FAC_2
      83             :             end if
      84             : #elif       P15_ENABLED
      85       10472 :             if (GET_ELEMENT(logzplus1) < LOGZ1PLUS1) then
      86         276 :                 GET_ELEMENT(logRateDensity) = 0._RKC
      87             :             else
      88        4720 :                 GET_ELEMENT(logRateDensity) = GET_ELEMENT(logzplus1) * EXPONENT_HIGH_Z + LOG_NORM_FAC_2
      89             :             end if
      90             : #endif
      91             : #if     D1_ENABLED
      92             :         end do
      93             : #endif
      94             : 
      95             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      96             : #elif   getLogRateDensity_ENABLED && (M14_ENABLED || M17_ENABLED || F18_ENABLED)
      97             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      98             : 
      99             : #if     M14_ENABLED
     100             :         real(RKC), parameter :: LOG_AMPLITUDE = log(0.015_RKC)
     101             :         real(RKC), parameter :: EXPLONENT_LOWER = 2.7_RKC
     102             :         real(RKC), parameter :: EXPLONENT_UPPER = 5.6_RKC
     103             :         real(RKC), parameter :: ZPLUS1_BREAK = 2.9_RKC
     104             :         real(RKC), parameter :: ZPLUS1_COEFF = 1._RKC / (ZPLUS1_BREAK**EXPLONENT_UPPER)
     105             : #elif   M17_ENABLED
     106             :         real(RKC), parameter :: LOG_AMPLITUDE = log(0.01_RKC)
     107             :         real(RKC), parameter :: EXPLONENT_LOWER = 2.6_RKC
     108             :         real(RKC), parameter :: EXPLONENT_UPPER = 6.2_RKC
     109             :         real(RKC), parameter :: ZPLUS1_BREAK = 3.2_RKC
     110             :         real(RKC), parameter :: ZPLUS1_COEFF = 1._RKC / (ZPLUS1_BREAK**EXPLONENT_UPPER)
     111             : #elif   F18_ENABLED
     112             :         real(RKC), parameter :: LOG_AMPLITUDE = log(0.013_RKC)
     113             :         real(RKC), parameter :: EXPLONENT_LOWER = 2.99_RKC
     114             :         real(RKC), parameter :: EXPLONENT_UPPER = 6.19_RKC
     115             :         real(RKC), parameter :: ZPLUS1_BREAK = 2.63_RKC
     116             :         real(RKC), parameter :: ZPLUS1_COEFF = 1._RKC / (ZPLUS1_BREAK**EXPLONENT_UPPER)
     117             : #endif
     118       36305 :         CHECK_ASSERTION(__LINE__, zplus1 >= 1._RKC, SK_": The condition `zplus1 >= 1.` must hold. logzplus1 = "//getStr(zplus1)) ! fpp
     119       36305 :         CHECK_ASSERTION(__LINE__, logzplus1 >= 0._RKC, SK_": The condition `logzplus1 >= 0.` must hold. logzplus1 = "//getStr(logzplus1)) ! fpp
     120       36305 :         logRateDensity = LOG_AMPLITUDE + EXPLONENT_LOWER * logzplus1 - log(1._RKC + ZPLUS1_COEFF * zplus1**EXPLONENT_UPPER)
     121             : #else
     122             :         !%%%%%%%%%%%%%%%%%%%%%%%%
     123             : #error  "Unrecognized interface."
     124             :         !%%%%%%%%%%%%%%%%%%%%%%%%
     125             : #endif
     126             : #undef GET_ELEMENT