The ParaMonte Documentation Website
Current view: top level - kernel - StarFormation_mod.f90 (source / functions) Hit Total Coverage
Test: ParaMonte 1.5.1 :: MPI Parallel Kernel - Code Coverage Report Lines: 179 190 94.2 %
Date: 2021-01-08 13:07:16 Functions: 23 23 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       2             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       3             : !!!!
       4             : !!!!   MIT License
       5             : !!!!
       6             : !!!!   ParaMonte: plain powerful parallel Monte Carlo library.
       7             : !!!!
       8             : !!!!   Copyright (C) 2012-present, The Computational Data Science Lab
       9             : !!!!
      10             : !!!!   This file is part of the ParaMonte library.
      11             : !!!!
      12             : !!!!   Permission is hereby granted, free of charge, to any person obtaining a
      13             : !!!!   copy of this software and associated documentation files (the "Software"),
      14             : !!!!   to deal in the Software without restriction, including without limitation
      15             : !!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
      16             : !!!!   and/or sell copies of the Software, and to permit persons to whom the
      17             : !!!!   Software is furnished to do so, subject to the following conditions:
      18             : !!!!
      19             : !!!!   The above copyright notice and this permission notice shall be
      20             : !!!!   included in all copies or substantial portions of the Software.
      21             : !!!!
      22             : !!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
      23             : !!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
      24             : !!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
      25             : !!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
      26             : !!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
      27             : !!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
      28             : !!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
      29             : !!!!
      30             : !!!!   ACKNOWLEDGMENT
      31             : !!!!
      32             : !!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
      33             : !!!!   As per the ParaMonte library license agreement terms, if you use any parts of
      34             : !!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
      35             : !!!!   work (education/research/industry/development/...) by citing the ParaMonte
      36             : !!!!   library as described on this page:
      37             : !!!!
      38             : !!!!       https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
      39             : !!!!
      40             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      41             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      42             : 
      43             : !> \brief This module contains procedures for computing the Cosmic Star Formation Rate in the Universe.
      44             : !> \author Amir Shahmoradi
      45             : 
      46             : module StarFormation_mod
      47             : 
      48             :     use Constants_mod, only: IK, RK, PI, NEGINF_RK
      49             : 
      50             :     implicit none
      51             : 
      52             :     character(*), parameter :: MODULE_NAME = "@StarFormation_mod"
      53             : 
      54             :     real(RK)    , parameter :: ZMIN = 0.1e0_RK
      55             :     real(RK)    , parameter :: ZMAX = 1.0e2_RK
      56             : 
      57             :     abstract interface
      58             :     function getLogRate_proc(zplus1,logzplus1,twiceLogLumDisMpc) result(logRate)
      59             :         import :: RK, IK
      60             :         real(RK), intent(in)    :: zplus1, logzplus1, twiceLogLumDisMpc
      61             :         real(RK)                :: logRate
      62             :     end function getLogRate_proc
      63             :     end interface
      64             : 
      65             :     abstract interface
      66             :     function getRateDensity_proc(zplus1) result(rateDensity)
      67             :         import :: RK, IK
      68             :         real(RK), intent(in)    :: zplus1
      69             :         real(RK)                :: rateDensity
      70             :     end function getRateDensity_proc
      71             :     end interface
      72             : 
      73             :     abstract interface
      74             :     function getMergerDelayTimePDF_proc(mergerDelayTime) result(mergerDelayTimeProb)
      75             :         import :: RK, IK
      76             :         real(RK), intent(in)    :: mergerDelayTime
      77             :         real(RK)                :: mergerDelayTimeProb
      78             :     end function getMergerDelayTimePDF_proc
      79             :     end interface
      80             : 
      81             : #if defined OS_IS_WSL
      82             :         procedure(getMergerDelayTimePDF_proc), pointer  :: getMergerDelayTimePDF_WSL        !< This madness bypasses the Microsoft Subsystem for Linux Internal Function call GFortran Segmentation Fault error.
      83             :         procedure(getRateDensity_proc), pointer         :: getStarFormationRateDensity_WSL  !< This madness bypasses the Microsoft Subsystem for Linux Internal Function call GFortran Segmentation Fault error.
      84             :         real(RK)                                        :: maxRelativeErrorDefault_WSL      !< This madness bypasses the Microsoft Subsystem for Linux Internal Function call GFortran Segmentation Fault error.
      85             :         integer(IK)                                     :: nRefinementDefault_WSL           !< This madness bypasses the Microsoft Subsystem for Linux Internal Function call GFortran Segmentation Fault error.
      86             :         real(RK)                                        :: lookBackTimeRef_WSL              !< This madness bypasses the Microsoft Subsystem for Linux Internal Function call GFortran Segmentation Fault error.
      87             : #endif
      88             : 
      89             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      90             : 
      91             : contains
      92             : 
      93             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      94             : 
      95             :     !> \brief
      96             :     !> Return GRBFR density based on the formation rate estimates of Petrosian et al (2015).
      97           3 :     pure function getLogRateDensityP15(logzplus1) result(logDensitySFR)
      98             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
      99             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateDensityP15
     100             : #endif
     101             :         use Constants_mod, only: RK
     102             :         implicit none
     103             :         real(RK), intent(in)    :: logzplus1
     104             :         real(RK), parameter     :: logz1plus1 = log(1._RK+4.50_RK)
     105             :         real(RK), parameter     :: exponentHighZ = -7.8_RK
     106             :         real(RK), parameter     :: logNormFac2 = -exponentHighZ * logz1plus1
     107             :         real(RK)                :: logDensitySFR
     108           3 :         if (logzplus1<0._RK) then
     109           0 :             logDensitySFR = NEGINF_RK
     110           3 :         elseif (logzplus1<logz1plus1) then
     111           0 :             logDensitySFR = 0._RK
     112             :         else
     113           3 :             logDensitySFR = logzplus1*exponentHighZ + logNormFac2
     114             :         end if
     115           3 :     end function getLogRateDensityP15
     116             : 
     117             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     118             : 
     119             :     !> \brief
     120             :     !> Return GRBFR density based on the formation rate estimates of Hopkins and Beacom (2007).
     121           3 :     pure function getLogRateDensityH06(logzplus1) result(logDensitySFR)
     122             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     123             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateDensityH06
     124             : #endif
     125           3 :         use Constants_mod, only: RK
     126             :         implicit none
     127             :         real(RK), intent(in)    :: logzplus1
     128             :         real(RK), parameter     :: logz0plus1 = log(1._RK+0.97_RK)
     129             :         real(RK), parameter     :: logz1plus1 = log(1._RK+4.50_RK)
     130             :         real(RK), parameter     :: g0 = +3.4_RK
     131             :         real(RK), parameter     :: g1 = -0.3_RK
     132             :         real(RK), parameter     :: g2 = -7.8_RK
     133             :         real(RK), parameter     :: logNormFac1 = logz0plus1*(g0-g1)
     134             :         real(RK), parameter     :: logNormFac2 = logz1plus1*(g1-g2) + logNormFac1
     135             :         real(RK)                :: logDensitySFR
     136           3 :         if (logzplus1<0._RK) then
     137           0 :             logDensitySFR = NEGINF_RK
     138           3 :         elseif (logzplus1<logz0plus1) then
     139           0 :             logDensitySFR = logzplus1*g0
     140           3 :         elseif (logzplus1<logz1plus1) then
     141           0 :             logDensitySFR = logzplus1*g1 + logNormFac1
     142             :         else
     143           3 :             logDensitySFR = logzplus1*g2 + logNormFac2
     144             :         end if
     145           3 :     end function getLogRateDensityH06
     146             : 
     147             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     148             : 
     149             :     !> \brief
     150             :     !> Return GRBFR density based on the formation rate estimates of Li (2008).
     151           3 :     pure function getLogRateDensityL08(logzplus1) result(logDensitySFR)
     152             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     153             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateDensityL08
     154             : #endif
     155           3 :         use Constants_mod, only: RK
     156             :         implicit none
     157             :         real(RK), intent(in)    :: logzplus1
     158             :         real(RK), parameter     :: logz0plus1 = log(1._RK+0.993_RK)
     159             :         real(RK), parameter     :: logz1plus1 = log(1._RK+3.800_RK)
     160             :         real(RK), parameter     :: g0 = +3.3000_RK
     161             :         real(RK), parameter     :: g1 = +0.0549_RK
     162             :         real(RK), parameter     :: g2 = -4.4600_RK
     163             :         real(RK), parameter     :: logNormFac1 = logz0plus1*(g0-g1)
     164             :         real(RK), parameter     :: logNormFac2 = logz1plus1*(g1-g2) + logNormFac1
     165             :         real(RK)                :: logDensitySFR
     166           3 :         if (logzplus1<0._RK) then
     167           0 :             logDensitySFR = NEGINF_RK
     168           3 :         elseif (logzplus1<logz0plus1) then
     169           0 :             logDensitySFR = logzplus1*g0
     170           3 :         elseif (logzplus1<logz1plus1) then
     171           0 :             logDensitySFR = logzplus1*g1 + logNormFac1
     172             :         else
     173           3 :             logDensitySFR = logzplus1*g2 + logNormFac2
     174             :         end if
     175           3 :     end function getLogRateDensityL08
     176             : 
     177             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     178             : 
     179             :     !> \brief
     180             :     !> Return GRBFR density based on the formation rate estimates of Butler et al. (2010).
     181             :     !>
     182             :     !> \remark
     183             :     !> The formation estimate of B10 takes into acount the metalicity-correction to the SFR.
     184        3891 :     pure function getLogRateDensityB10(logzplus1) result(logDensitySFR)
     185             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     186             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateDensityB10
     187             : #endif
     188           3 :         use Constants_mod, only: IK, RK
     189             :         implicit none
     190             :         real(RK), intent(in)    :: logzplus1
     191             :         real(RK), parameter     :: logz0plus1 = log(1._RK+0.97_RK)
     192             :         real(RK), parameter     :: logz1plus1 = log(1._RK+4.00_RK)
     193             :         real(RK), parameter     :: g0 = +3.14_RK
     194             :         real(RK), parameter     :: g1 = +1.36_RK
     195             :         real(RK), parameter     :: g2 = -2.92_RK
     196             :         real(RK), parameter     :: logNormFac1 = logz0plus1*(g0-g1)
     197             :         real(RK), parameter     :: logNormFac2 = logz1plus1*(g1-g2) + logNormFac1
     198             :         real(RK)                :: logDensitySFR
     199        3891 :         if (logzplus1<0._RK) then
     200           0 :             logDensitySFR = NEGINF_RK
     201        3891 :         elseif (logzplus1<logz0plus1) then
     202           0 :             logDensitySFR = logzplus1*g0
     203        3891 :         elseif (logzplus1<logz1plus1) then
     204        1194 :             logDensitySFR = logzplus1*g1 + logNormFac1
     205             :         else
     206        2697 :             logDensitySFR = logzplus1*g2 + logNormFac2
     207             :         end if
     208        3891 :     end function getLogRateDensityB10
     209             : 
     210             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     211             : 
     212             :     !> \brief
     213             :     !> Return the Comoving Star Formation Rate Density according to Eqn 15 of Madau 2014: Cosmic Star-Formation History.\n
     214             :     !> `densitySFR(z) = 0.015 * (1+z)^2.7 / ( 1 + [(1+z)/2.9]^5.6 )`
     215           3 :     pure function getLogRateDensityM14(zplus1,logzplus1) result(logDensitySFR)
     216             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     217             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateDensityM14
     218             : #endif
     219        3891 :         use Constants_mod, only: RK
     220             :         implicit none
     221             :         real(RK), intent(in)    :: zplus1, logzplus1
     222             :         real(RK)                :: logDensitySFR
     223             :         real(RK), parameter     :: logAmplitude = log(0.015_RK)
     224             :         real(RK), parameter     :: lowerExp = 2.7_RK
     225             :         real(RK), parameter     :: upperExp = 5.6_RK
     226             :         real(RK), parameter     :: zplus1Break = 2.9_RK
     227             :         real(RK), parameter     :: zplus1Coeff = 1._RK / (zplus1Break**upperExp)
     228           3 :         logDensitySFR = logAmplitude + lowerExp*logzplus1 - log( 1._RK + zplus1Coeff * zplus1**upperExp )
     229           6 :     end function getLogRateDensityM14
     230             : 
     231             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     232             : 
     233             :     !> \brief
     234             :     !> Return the Comoving Star Formation Rate Density according to Eqn 1 of Madau 2017: Cosmic Star-Formation History
     235             :     !> `densitySFR(z) = 0.01 * (1+z)^2.6 / ( 1 + [(1+z)/3.2]^6.2 )`
     236           3 :     pure function getLogRateDensityM17(zplus1,logzplus1) result(logDensitySFR)
     237             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     238             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateDensityM17
     239             : #endif
     240           3 :         use Constants_mod, only: RK
     241             :         implicit none
     242             :         real(RK), intent(in)    :: zplus1, logzplus1
     243             :         real(RK)                :: logDensitySFR
     244             :         real(RK), parameter     :: logAmplitude = log(0.01_RK)
     245             :         real(RK), parameter     :: lowerExp = 2.6_RK
     246             :         real(RK), parameter     :: upperExp = 6.2_RK
     247             :         real(RK), parameter     :: zplus1Break = 3.2_RK
     248             :         real(RK), parameter     :: zplus1Coeff = 1._RK / (zplus1Break**upperExp)
     249           3 :         logDensitySFR = logAmplitude + lowerExp*logzplus1 - log( 1._RK + zplus1Coeff * zplus1**upperExp )
     250           6 :     end function getLogRateDensityM17
     251             : 
     252             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     253             : 
     254             :     !> \brief
     255             :     !> Return the Mordau Comoving Star Formation Rate Density with updated parameters from Fermi 2018.\n
     256             :     ! `densitySFR(z) = 0.013 * (1+z)^2.99 / ( 1 + [(1+z)/2.63]^6.19 )`
     257           3 :     pure function getLogRateDensityF18(zplus1,logzplus1) result(logDensitySFR)
     258             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     259             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateDensityF18
     260             : #endif
     261           3 :         use Constants_mod, only: RK
     262             :         implicit none
     263             :         real(RK), intent(in)    :: zplus1, logzplus1
     264             :         real(RK)                :: logDensitySFR
     265             :         real(RK), parameter     :: logAmplitude = log(0.013_RK)
     266             :         real(RK), parameter     :: lowerExp = 2.99_RK
     267             :         real(RK), parameter     :: upperExp = 6.19_RK
     268             :         real(RK), parameter     :: zplus1Break = 2.63_RK
     269             :         real(RK), parameter     :: zplus1Coeff = 1._RK / (zplus1Break**upperExp)
     270           3 :         logDensitySFR = logAmplitude + lowerExp*logzplus1 - log( 1._RK + zplus1Coeff * zplus1**upperExp )
     271           6 :     end function getLogRateDensityF18
     272             : 
     273             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     274             : 
     275             :     !> \brief
     276             :     !> Return the cosmic formation rate according to the work of Petrosian et al. (2015).
     277           3 :     pure function getLogRateP15(zplus1,logzplus1,twiceLogLumDisMpc) result(logRate)
     278             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     279             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateP15
     280             : #endif
     281           3 :         use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
     282             :         use Constants_mod, only: RK, PI
     283             :         implicit none
     284             :         real(RK), intent(in)    :: zplus1, logzplus1, twiceLogLumDisMpc
     285             :         real(RK), parameter     :: LOG_COEF = log(4._RK*PI*LS2HC)
     286             :         real(RK)                :: logRate
     287             :         logRate = LOG_COEF + twiceLogLumDisMpc - ( 3._RK*logzplus1 + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) ) &
     288           3 :                 + getLogRateDensityP15(logzplus1)
     289           6 :     end function getLogRateP15
     290             : 
     291             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     292             : 
     293             :     !> \brief
     294             :     !> Return the cosmic star formation rate according to the work of Hopkins and Beacom (2007).
     295           3 :     pure function getLogRateH06(zplus1,logzplus1,twiceLogLumDisMpc) result(logRate)
     296             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     297             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateH06
     298             : #endif
     299           3 :         use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
     300             :         use Constants_mod, only: RK, PI
     301             :         implicit none
     302             :         real(RK), intent(in)    :: zplus1, logzplus1, twiceLogLumDisMpc
     303             :         real(RK), parameter     :: LOG_COEF = log(4._RK*PI*LS2HC)
     304             :         real(RK)                :: logRate
     305             :         logRate = LOG_COEF + twiceLogLumDisMpc - ( 3._RK*logzplus1 + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) ) &
     306           3 :                 + getLogRateDensityH06(logzplus1)
     307           6 :     end function getLogRateH06
     308             : 
     309             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     310             : 
     311             :     !> \brief
     312             :     !> Return the cosmic star formation rate according to the work of Li (2008).
     313           3 :     pure function getLogRateL08(zplus1,logzplus1,twiceLogLumDisMpc) result(logRate)
     314             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     315             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateL08
     316             : #endif
     317           3 :         use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
     318             :         use Constants_mod, only: RK, PI
     319             :         implicit none
     320             :         real(RK), intent(in)    :: zplus1, logzplus1, twiceLogLumDisMpc
     321             :         real(RK), parameter     :: LOG_COEF = log(4._RK*PI*LS2HC)
     322             :         real(RK)                :: logRate
     323             :         logRate = LOG_COEF + twiceLogLumDisMpc - ( 3._RK*logzplus1 + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) ) &
     324           3 :                 + getLogRateDensityL08(logzplus1)
     325           6 :     end function getLogRateL08
     326             : 
     327             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     328             : 
     329             :     !> \brief
     330             :     !> Return the cosmic GRB formation rate according to the work of Butler et al. (2010).
     331           3 :     pure function getLogRateB10(zplus1,logzplus1,twiceLogLumDisMpc) result(logRate)
     332             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     333             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateB10
     334             : #endif
     335           3 :         use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
     336             :         use Constants_mod, only: RK, PI
     337             :         implicit none
     338             :         real(RK), intent(in)    :: zplus1, logzplus1, twiceLogLumDisMpc
     339             :         real(RK), parameter     :: LOG_COEF = log(4._RK*PI*LS2HC)
     340             :         real(RK)                :: logRate
     341             :         logRate = LOG_COEF + twiceLogLumDisMpc - ( 3._RK*logzplus1 + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) ) &
     342           3 :                 + getLogRateDensityB10(logzplus1)
     343           6 :     end function getLogRateB10
     344             : 
     345             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     346             : 
     347             :     !> \brief
     348             :     !> Return the cosmic star formation rate according to the work of Madau et al. (2014).
     349           3 :     pure function getLogRateM14(zplus1,logzplus1,twiceLogLumDisMpc) result(logRate)
     350             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     351             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateM14
     352             : #endif
     353           3 :         use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
     354             :         use Constants_mod, only: RK, PI
     355             :         implicit none
     356             :         real(RK), intent(in)    :: zplus1, logzplus1, twiceLogLumDisMpc
     357             :         real(RK), parameter     :: LOG_COEF = log(4._RK*PI*LS2HC)
     358             :         real(RK)                :: logRate
     359             :         logRate = LOG_COEF + twiceLogLumDisMpc - ( 3._RK*logzplus1 + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) ) &
     360           3 :                 + getLogRateDensityM14(zplus1,logzplus1)
     361           6 :     end function getLogRateM14
     362             : 
     363             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     364             : 
     365             :     !> \brief
     366             :     !> Return the cosmic star formation rate according to the work of Madau et al. (2017).
     367           3 :     pure function getLogRateM17(zplus1,logzplus1,twiceLogLumDisMpc) result(logRate)
     368             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     369             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateM17
     370             : #endif
     371           3 :         use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
     372             :         use Constants_mod, only: RK, PI
     373             :         implicit none
     374             :         real(RK), intent(in)    :: zplus1, logzplus1, twiceLogLumDisMpc
     375             :         real(RK), parameter     :: LOG_COEF = log(4._RK*PI*LS2HC)
     376             :         real(RK)                :: logRate
     377             :         logRate = LOG_COEF + twiceLogLumDisMpc - ( 3._RK*logzplus1 + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) ) &
     378           3 :                 + getLogRateDensityM17(zplus1,logzplus1)
     379           6 :     end function getLogRateM17
     380             : 
     381             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     382             : 
     383             :     !> \brief
     384             :     !> Return the cosmic star formation rate according to the work of the Fermi collaboration (2018).
     385           3 :     pure function getLogRateF18(zplus1,logzplus1,twiceLogLumDisMpc) result(logRate)
     386             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     387             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateF18
     388             : #endif
     389           3 :         use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
     390             :         use Constants_mod, only: RK, PI
     391             :         implicit none
     392             :         real(RK), intent(in)    :: zplus1, logzplus1, twiceLogLumDisMpc
     393             :         real(RK), parameter     :: LOG_COEF = log(4._RK*PI*LS2HC)
     394             :         real(RK)                :: logRate
     395             :         logRate = LOG_COEF + twiceLogLumDisMpc - ( 3._RK*logzplus1 + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) ) &
     396           3 :                 + getLogRateDensityF18(zplus1,logzplus1)
     397           6 :     end function getLogRateF18
     398             : 
     399             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     400             : 
     401             :     ! Amir Shahmoradi, Wednesday 5:43 PM, December 25, 2013, IFS, UT Austin
     402             :     ! compute the binary merger rate as a function of z, according to Shahmoradi and Nemiroff (2015)
     403             :     ! equivalent to delayed_rate_Belz_Li(z) in S15
     404             :     ! returns 0, if z>6.501_RK or z<0.09_RK
     405          12 :     pure function getBinaryMergerRateS15(z) result(binaryMergerRateS15)
     406             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     407             :         !DEC$ ATTRIBUTES DLLEXPORT :: getBinaryMergerRateS15
     408             : #endif
     409           3 :         use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
     410             :         use Constants_mod, only: RK, PI
     411             :         implicit none
     412             :         real(RK)    , intent(in)    :: z
     413             :         real(RK)                    :: binaryMergerRateS15
     414          12 :         if (z>2.5_RK .and. z<=6.501_RK) then
     415             :             binaryMergerRateS15 = &
     416             :                                 - 2.09118024744342000_RK &
     417             :                                 + 5.15382361299299000_RK * z &
     418             :                                 - 5.46442271664195000_RK * z**2 &
     419             :                                 + 3.29445310883082000_RK * z**3 &
     420             :                                 - 1.24547016168265000_RK * z**4 &
     421             :                                 + 0.30628893690508400_RK * z**5 &
     422             :                                 - 0.04904403249641820_RK * z**6 &
     423             :                                 + 0.00493757380504717_RK * z**7 &
     424             :                                 - 2.84061971928750e-4_RK * z**8 &
     425           3 :                                 + 7.12674138757750e-6_RK * z**9
     426           9 :         elseif (z>1.0_RK .and. z<=2.5_RK) then
     427             :             binaryMergerRateS15 = &
     428             :                                 - 0.86022576265904100_RK &
     429             :                                 + 4.22669545558817000_RK * z &
     430             :                                 - 8.86086728534670000_RK * z**2 &
     431             :                                 + 10.4863792284648000_RK * z**3 &
     432             :                                 - 7.64722909221129000_RK * z**4 &
     433             :                                 + 3.51616699500767000_RK * z**5 &
     434             :                                 - 0.99555474471022000_RK * z**6 &
     435             :                                 + 0.15876893754371900_RK * z**7 &
     436           3 :                                 - 0.01092541997736420_RK * z**8
     437           6 :         elseif (z<=1._RK .and. z>=0.09_RK) then
     438             :             binaryMergerRateS15 = &
     439             :                                 + 1.92595299989370e-4_RK &
     440             :                                 - 0.00345273599582578_RK * z &
     441             :                                 + 0.03157500615320920_RK * z**2 &
     442             :                                 - 0.04470545521198460_RK * z**3 &
     443             :                                 + 0.06812481521281660_RK * z**4 &
     444           0 :                                 - 0.03846033416253570_RK * z**5
     445             :         else
     446           6 :             binaryMergerRateS15 = 0._RK
     447             :         end if
     448          12 :     end function getBinaryMergerRateS15
     449             : 
     450             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     451             : 
     452             :     ! computes the natural log of the binary merger rate as a function of logzplus1.
     453             :     ! However, note that the computed rate is dN/dz, even though the input is logzplus1.
     454             :     ! the merger delay time distribution is the same as that of Shahmoradi and Nemiroff (2015), but with the logMean:log(0.1_RK) and sigma: 0.9612813_RK..
     455             :     ! returns 0, if z>19.929999999999882_RK or z<0.03_RK
     456          21 :     pure function getLogBinaryMergerRateLognormH06(logzplus1) result(logBinaryMergerRate)
     457             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     458             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogBinaryMergerRateLognormH06
     459             : #endif
     460          12 :         use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
     461             :         use Constants_mod, only: RK, PI
     462             :         implicit none
     463             :         real(RK)    , intent(in)    :: logzplus1
     464             :         real(RK)                    :: logBinaryMergerRate
     465          21 :         if (logzplus1>0.02955880224154443_RK .and. logzplus1<=0.1441003439737565_RK) then
     466             :             logBinaryMergerRate = &
     467             :                                 - 14.26464149493092000_RK &
     468             :                                 + 84.73477757043948000_RK * logzplus1 &
     469             :                                 - 488.5893985602366500_RK * logzplus1**2 &
     470           3 :                                 + 1154.414655194473900_RK * logzplus1**3
     471          18 :         elseif (logzplus1>0.1441003439737565_RK .and. logzplus1<=0.6575200029167926_RK) then
     472             :             logBinaryMergerRate = &
     473             :                                 - 11.19700066921606300_RK &
     474             :                                 + 20.46712963401572300_RK * logzplus1 &
     475             :                                 - 24.31794334813894300_RK * logzplus1**2 &
     476           3 :                                 + 12.21213317590724400_RK * logzplus1**3
     477          15 :         elseif (logzplus1>0.6575200029167926_RK .and. logzplus1<=1.5591966959973538_RK) then
     478             :             logBinaryMergerRate = &
     479             :                                 - 9.094912666461765000_RK &
     480             :                                 + 15.23119806754538900_RK * logzplus1 &
     481             :                                 - 18.77526325204311800_RK * logzplus1**2 &
     482             :                                 + 9.941360355936961000_RK * logzplus1**3 &
     483           3 :                                 - 2.077370913197473000_RK * logzplus1**4
     484          12 :         elseif (logzplus1>1.5591966959973538_RK .and. logzplus1<=1.7056567701746455_RK) then
     485             :             logBinaryMergerRate = &
     486             :                                 - 2392.907733171019000_RK &
     487             :                                 + 6210.872744126407000_RK * logzplus1 &
     488             :                                 - 6054.866136454215000_RK * logzplus1**2 &
     489             :                                 + 2622.628785434413700_RK * logzplus1**3 &
     490           3 :                                 - 426.0273477222719000_RK * logzplus1**4
     491           9 :         elseif (logzplus1>1.7056567701746455_RK .and. logzplus1<=3.0411835364579027_RK) then
     492             :             logBinaryMergerRate = &
     493             :                                 + 9.538876239886940000_RK &
     494             :                                 - 8.753418172517534000_RK * logzplus1 &
     495             :                                 - 0.159980818030374640_RK * logzplus1**2 &
     496           3 :                                 - 0.088551503657680930_RK * logzplus1**3
     497           6 :         elseif (logzplus1<=0.02955880224154443_RK .or. logzplus1>3.0411835364579027_RK) then
     498           6 :             logBinaryMergerRate = 0._RK
     499             :         end if
     500          21 :     end function getLogBinaryMergerRateLognormH06
     501             : 
     502             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     503             : 
     504             :     ! computes the natural log of the binary merger rate as a function of logzplus1.
     505             :     ! However, note that the computed rate is dN/dz, even though the input is logzplus1.
     506             :     ! the merger delay time distribution is the same as that of Shahmoradi and Nemiroff (2015), but with the logMean:log(0.1_RK) and sigma: 0.9612813_RK..
     507             :     ! returns 0, if z>19.929999999999882_RK or z<0.03_RK
     508          21 :     pure function getLogBinaryMergerRateLognormL08(logzplus1) result(logBinaryMergerRate)
     509             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     510             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogBinaryMergerRateLognormL08
     511             : #endif
     512          21 :         use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
     513             :         use Constants_mod, only: RK, PI
     514             :         implicit none
     515             :         real(RK)    , intent(in)    :: logzplus1
     516             :         real(RK)                    :: logBinaryMergerRate
     517          21 :         if (logzplus1>0.02955880224154443_RK .and. logzplus1<=0.20701416938432557_RK) then
     518             :             logBinaryMergerRate = &
     519             :                                 - 14.53696144309043900_RK &
     520             :                                 + 94.70274747509626000_RK * logzplus1 &
     521             :                                 - 687.3663996060040000_RK * logzplus1**2 &
     522             :                                 + 2695.421036673770700_RK * logzplus1**3 &
     523           3 :                                 - 4077.601561165490000_RK * logzplus1**4
     524          18 :         elseif (logzplus1>0.20701416938432557_RK .and. logzplus1<=0.8241754429663476_RK) then
     525             :             logBinaryMergerRate = &
     526             :                                 - 13.5104005566057670_RK &
     527             :                                 + 49.6443928683743600_RK * logzplus1 &
     528             :                                 - 164.286063097338630_RK * logzplus1**2 &
     529             :                                 + 315.721394966368100_RK * logzplus1**3 &
     530             :                                 - 300.345052726248640_RK * logzplus1**4 &
     531           3 :                                 + 108.470535327547080_RK * logzplus1**5
     532          15 :         elseif (logzplus1>0.8241754429663476_RK .and. logzplus1<=1.4243124283074096_RK) then
     533             :             logBinaryMergerRate = &
     534             :                                 - 8.77634469738400500_RK &
     535             :                                 + 13.1999684738558810_RK * logzplus1 &
     536             :                                 - 15.8698236818922140_RK * logzplus1**2 &
     537             :                                 + 8.48676936452957000_RK * logzplus1**3 &
     538           3 :                                 - 1.83190451512279620_RK * logzplus1**4
     539          12 :         elseif (logzplus1>1.4243124283074096_RK .and. logzplus1<=1.6154199841116488_RK) then
     540             :             logBinaryMergerRate = &
     541             :                                 + 4158.29353781047900_RK &
     542             :                                 - 10954.1105856433040_RK * logzplus1 &
     543             :                                 + 10789.3451136201870_RK * logzplus1**2 &
     544             :                                 - 4713.80244702217800_RK * logzplus1**3 &
     545           3 :                                 + 770.488645040204600_RK * logzplus1**4
     546           9 :         elseif (logzplus1>1.6154199841116488_RK .and. logzplus1<=3.0411835364579027_RK) then
     547             :             logBinaryMergerRate = &
     548             :                                 + 0.37742655174185624_RK &
     549             :                                 + 0.30883738015163340_RK * logzplus1 &
     550             :                                 - 4.04937550957291800_RK * logzplus1**2 &
     551             :                                 + 1.11680537027038170_RK * logzplus1**3 &
     552           3 :                                 - 0.13770838345089523_RK * logzplus1**4
     553           6 :         elseif (logzplus1<=0.02955880224154443_RK .or. logzplus1>3.0411835364579027_RK) then
     554           6 :             logBinaryMergerRate = 0._RK
     555             :         end if
     556          21 :     end function getLogBinaryMergerRateLognormL08
     557             : 
     558             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     559             : 
     560             :     ! computes the natural log of the binary merger rate as a function of logzplus1.
     561             :     ! However, note that the computed rate is dN/dz, even though the input is logzplus1.
     562             :     ! the merger delay time distribution is the same as that of Shahmoradi and Nemiroff (2015), but with the logMean:log(0.1_RK) and sigma: 0.9612813_RK.
     563             :     ! returns 0, if z>19.929999999999882_RK or z<0.03_RK
     564          21 :     pure function getLogBinaryMergerRateLognormB10(logzplus1) result(logBinaryMergerRate)
     565             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     566             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogBinaryMergerRateLognormB10
     567             : #endif
     568          21 :         use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
     569             :         use Constants_mod, only: RK, PI
     570             :         implicit none
     571             :         real(RK)    , intent(in)    :: logzplus1
     572             :         real(RK)                    :: logBinaryMergerRate
     573          21 :         if (logzplus1>0.02955880224154443_RK .and. logzplus1<=0.20701416938432557_RK) then
     574             :             logBinaryMergerRate = &
     575             :                                 - 15.27802857671202000_RK &
     576             :                                 + 94.54179164991284000_RK * logzplus1 &
     577             :                                 - 687.3676159275769000_RK * logzplus1**2 &
     578             :                                 + 2695.420977251770600_RK * logzplus1**3 &
     579           3 :                                 - 4077.601406650646000_RK * logzplus1**4
     580          18 :         elseif (logzplus1>0.20701416938432557_RK .and. logzplus1<=0.8241754429663476_RK) then
     581             :             logBinaryMergerRate = &
     582             :                                 - 13.5066182170954650_RK &
     583             :                                 + 40.1985219822299200_RK * logzplus1 &
     584             :                                 - 121.506350703598660_RK * logzplus1**2 &
     585             :                                 + 224.621285123736100_RK * logzplus1**3 &
     586             :                                 - 210.878836655472500_RK * logzplus1**4 &
     587           3 :                                 + 76.3335749498628400_RK * logzplus1**5
     588          15 :         elseif (logzplus1>0.8241754429663476_RK .and. logzplus1<=1.4243124283074096_RK) then
     589             :             logBinaryMergerRate = &
     590             :                                 - 10.0515447816113700_RK &
     591             :                                 + 12.6659826494097970_RK * logzplus1 &
     592             :                                 - 13.2268991886238200_RK * logzplus1**2 &
     593             :                                 + 6.84523627043807100_RK * logzplus1**3 &
     594           3 :                                 - 1.44645280124922220_RK * logzplus1**4
     595          12 :         elseif (logzplus1>1.4243124283074096_RK .and. logzplus1<=1.6104374127671848_RK) then
     596             :             logBinaryMergerRate = &
     597             :                                 - 1187.90539057029950_RK &
     598             :                                 + 3240.19327021926350_RK * logzplus1 &
     599             :                                 - 3330.70645904271000_RK * logzplus1**2 &
     600             :                                 + 1522.87499612399850_RK * logzplus1**3 &
     601           3 :                                 - 261.341408956542300_RK * logzplus1**4
     602           9 :         elseif (logzplus1>1.6104374127671848_RK .and. logzplus1<=3.0411835364579027_RK) then
     603             :             logBinaryMergerRate = &
     604             :                                 - 1.43934839576471260_RK &
     605             :                                 + 1.72951867017028120_RK * logzplus1 &
     606             :                                 - 4.06729555225025000_RK * logzplus1**2 &
     607             :                                 + 1.18253386764330200_RK * logzplus1**3 &
     608           3 :                                 - 0.15201156018584210_RK * logzplus1**4
     609           6 :         elseif (logzplus1<=0.02955880224154443_RK .or. logzplus1>3.0411835364579027_RK) then
     610           6 :             logBinaryMergerRate = 0._RK
     611             :         end if
     612          21 :     end function getLogBinaryMergerRateLognormB10
     613             : 
     614             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     615             : 
     616             :     ! computes the natural log of the binary merger rate as a function of logzplus1.
     617             :     ! However, note that the computed rate is dN/dz, even though the input is logzplus1.
     618             :     ! the merger delay time distribution is the same as that of Shahmoradi and Nemiroff (2015), but with the logMean:log(0.1_RK) and sigma: 0.9612813_RK..
     619             :     ! returns 0, if z>19.929999999999882_RK or z<0.03_RK
     620          21 :     pure function getLogBinaryMergerRateLognormM14(logzplus1) result(logBinaryMergerRate)
     621             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     622             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogBinaryMergerRateLognormM14
     623             : #endif
     624          21 :         use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
     625             :         use Constants_mod, only: RK, PI
     626             :         implicit none
     627             :         real(RK)    , intent(in)    :: logzplus1
     628             :         real(RK)                    :: logBinaryMergerRate
     629          21 :         if (logzplus1>0.02955880224154443_RK .and. logzplus1<=0.16551443847757297_RK) then
     630             :             logBinaryMergerRate = &
     631             :                                 - 13.91129314580349600_RK &
     632             :                                 + 78.88963489621422000_RK * logzplus1 &
     633             :                                 - 420.9801740859396700_RK * logzplus1**2 &
     634           3 :                                 + 902.4783078800951000_RK * logzplus1**3
     635          18 :         elseif (logzplus1>0.16551443847757297_RK .and. logzplus1<=0.9282193027394269_RK) then
     636             :             logBinaryMergerRate = &
     637             :                                 - 11.00951046136480500_RK &
     638             :                                 + 21.38817515748999000_RK * logzplus1 &
     639             :                                 - 33.29451048508970000_RK * logzplus1**2 &
     640             :                                 + 29.32158835244860400_RK * logzplus1**3 &
     641           3 :                                 - 10.97377449040440000_RK * logzplus1**4
     642          15 :         elseif (logzplus1>0.9282193027394269_RK .and. logzplus1<=1.3937663759585892_RK) then
     643             :             logBinaryMergerRate = &
     644             :                                 - 8.254476015464371000_RK &
     645             :                                 + 3.620963332444886000_RK * logzplus1 &
     646             :                                 + 6.734585433384001000_RK * logzplus1**2 &
     647             :                                 - 9.151412394211048000_RK * logzplus1**3 &
     648           3 :                                 + 2.516171777428496000_RK * logzplus1**4
     649          12 :         elseif (logzplus1>1.3937663759585892_RK .and. logzplus1<=3.0411835364579027_RK) then
     650             :             logBinaryMergerRate = &
     651             :                                 - 6.539697727782377000_RK &
     652             :                                 + 8.522331572601950000_RK * logzplus1 &
     653             :                                 - 8.242990979412244000_RK * logzplus1**2 &
     654             :                                 + 2.316632169715435300_RK * logzplus1**3 &
     655           6 :                                 - 0.266462340853027450_RK * logzplus1**4
     656           6 :         elseif (logzplus1<=0.02955880224154443_RK .or. logzplus1>3.0411835364579027_RK) then
     657           6 :             logBinaryMergerRate = 0._RK
     658             :         end if
     659          21 :     end function getLogBinaryMergerRateLognormM14
     660             : 
     661             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     662             : 
     663             :     ! computes the natural log of the binary merger rate as a function of logzplus1.
     664             :     ! However, note that the computed rate is dN/dz, even though the input is logzplus1.
     665             :     ! the merger delay time distribution is the same as that of Shahmoradi and Nemiroff (2015), but with the logMean:log(0.1_RK) and sigma: 0.9612813_RK..
     666             :     ! returns 0, if z>19.929999999999882_RK or z<0.03_RK
     667          21 :     pure function getLogBinaryMergerRateLognormM17(logzplus1) result(logBinaryMergerRate)
     668             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     669             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogBinaryMergerRateLognormM17
     670             : #endif
     671          21 :         use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
     672             :         use Constants_mod, only: RK, PI
     673             :         implicit none
     674             :         real(RK)    , intent(in)    :: logzplus1
     675             :         real(RK)                    :: logBinaryMergerRate
     676          21 :         if (logzplus1>0.02955880224154443_RK .and. logzplus1<=0.16551443847757297_RK) then
     677             :             logBinaryMergerRate = &
     678             :                                 - 14.01939141013502300_RK &
     679             :                                 + 78.80010843737509000_RK * logzplus1 &
     680             :                                 - 420.9593253775164000_RK * logzplus1**2 &
     681           3 :                                 + 902.5668042795056000_RK * logzplus1**3
     682          18 :         elseif (logzplus1>0.16551443847757297_RK .and. logzplus1<=0.9282193027394269_RK) then
     683             :             logBinaryMergerRate = &
     684             :                                 - 11.12915953695671500_RK &
     685             :                                 + 21.43217730985805500_RK * logzplus1 &
     686             :                                 - 33.75904206577289000_RK * logzplus1**2 &
     687             :                                 + 30.03916282499634200_RK * logzplus1**3 &
     688           3 :                                 - 11.12086545981264500_RK * logzplus1**4
     689          15 :         elseif (logzplus1>0.9282193027394269_RK .and. logzplus1<=1.3937663759585892_RK) then
     690             :             logBinaryMergerRate = &
     691             :                                 - 1.802362223155230800_RK &
     692             :                                 - 20.58526172567768200_RK * logzplus1 &
     693             :                                 + 38.93828966743146000_RK * logzplus1**2 &
     694             :                                 - 27.19863916580484500_RK * logzplus1**3 &
     695           3 :                                 + 6.138928143113263000_RK * logzplus1**4
     696          12 :         elseif (logzplus1>1.3937663759585892_RK .and. logzplus1<=3.0411835364579027_RK) then
     697             :             logBinaryMergerRate = &
     698             :                                 - 7.711815956299844000_RK &
     699             :                                 + 11.68891993486079700_RK * logzplus1 &
     700             :                                 - 10.62908897824095300_RK * logzplus1**2 &
     701             :                                 + 2.945685425778300700_RK * logzplus1**3 &
     702           6 :                                 - 0.327069839977957850_RK * logzplus1**4
     703           6 :         elseif (logzplus1<=0.02955880224154443_RK .or. logzplus1>3.0411835364579027_RK) then
     704           6 :             logBinaryMergerRate = 0._RK
     705             :         end if
     706          21 :     end function getLogBinaryMergerRateLognormM17
     707             : 
     708             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     709             : 
     710             :     ! computes the natural log of the binary merger rate as a function of logzplus1.
     711             :     ! However, note that the computed rate is dN/dz, even though the input is logzplus1.
     712             :     ! the merger delay time distribution is the same as that of Shahmoradi and Nemiroff (2015), but with the logMean:log(0.1_RK) and sigma: 0.9612813_RK..
     713             :     ! returns 0, if z>19.929999999999882_RK or z<0.03_RK
     714          21 :     pure function getLogBinaryMergerRateLognormF18(logzplus1) result(logBinaryMergerRate)
     715             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     716             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogBinaryMergerRateLognormF18
     717             : #endif
     718          21 :         use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
     719             :         use Constants_mod, only: RK, PI
     720             :         implicit none
     721             :         real(RK)    , intent(in)    :: logzplus1
     722             :         real(RK)                    :: logBinaryMergerRate
     723          21 :         if (logzplus1>0.02955880224154443_RK .and. logzplus1<=0.16551443847757297_RK) then
     724             :             logBinaryMergerRate = &
     725             :                                 - 13.80128475140318000_RK &
     726             :                                 + 79.17963739241087000_RK * logzplus1 &
     727             :                                 - 420.9808813943490700_RK * logzplus1**2 &
     728           3 :                                 + 902.4149755380632000_RK * logzplus1**3
     729          18 :         elseif (logzplus1>0.16551443847757297_RK .and. logzplus1<=0.9282193027394269_RK) then
     730             :             logBinaryMergerRate = &
     731             :                                 - 10.89131941401880800_RK &
     732             :                                 + 21.59483273763076400_RK * logzplus1 &
     733             :                                 - 33.07662054750123000_RK * logzplus1**2 &
     734             :                                 + 29.23608723915102600_RK * logzplus1**3 &
     735           3 :                                 - 11.33984422193848700_RK * logzplus1**4
     736          15 :         elseif (logzplus1>0.9282193027394269_RK .and. logzplus1<=1.3937663759585892_RK) then
     737             :             logBinaryMergerRate = &
     738             :                                 - 14.02518830695731000_RK &
     739             :                                 + 24.92009885816913300_RK * logzplus1 &
     740             :                                 - 20.04762612951797000_RK * logzplus1**2 &
     741             :                                 + 4.885281389900379500_RK * logzplus1**3 &
     742           3 :                                 - 0.168402813838908260_RK * logzplus1**4
     743          12 :         elseif (logzplus1>1.3937663759585892_RK .and. logzplus1<=3.0411835364579027_RK) then
     744             :             logBinaryMergerRate = &
     745             :                                 - 4.348081430972656000_RK &
     746             :                                 + 4.815143234949144000_RK * logzplus1 &
     747             :                                 - 6.143880845780776000_RK * logzplus1**2 &
     748             :                                 + 1.738835623950871300_RK * logzplus1**3 &
     749           6 :                                 - 0.206972882929076480_RK * logzplus1**4
     750           6 :         elseif (logzplus1<=0.02955880224154443_RK .or. logzplus1>3.0411835364579027_RK) then
     751           6 :             logBinaryMergerRate = 0._RK
     752             :         end if
     753          21 :     end function getLogBinaryMergerRateLognormF18
     754             : 
     755             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     756             : 
     757          24 :     function getBinaryMergerRateDensity ( zplus1 &
     758             :                                         , zplus1Max &
     759             :                                         , nRefinement &
     760             :                                         , maxRelativeError &
     761             :                                         , getMergerDelayTimePDF &
     762             :                                         , getStarFormationRateDensity &
     763          24 :                                         ) result(binaryMergerRateDensity)
     764             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     765             :         !DEC$ ATTRIBUTES DLLEXPORT :: getBinaryMergerRateDensity
     766             : #endif
     767             :         use, intrinsic :: iso_fortran_env, only: output_unit
     768          21 :         use Cosmology_mod, only: getLookBackTime
     769             :         use Constants_mod, only: RK, HUGE_RK
     770             :         use Integration_mod, only: doQuadRombOpen, ErrorMessage!, midinf
     771             :         use Integration_mod, only: midexp
     772             :        !use Integration_mod, only: midinf
     773             :         implicit none
     774             :         real(RK)    , intent(in)                :: zplus1
     775             :         real(RK)    , intent(in), optional      :: zplus1Max, maxRelativeError
     776             :         integer(IK) , intent(in), optional      :: nRefinement
     777             :         procedure(getRateDensity_proc)          :: getStarFormationRateDensity
     778             :         procedure(getMergerDelayTimePDF_proc)   :: getMergerDelayTimePDF
     779             :         integer(IK)                             :: neval, ierr, nRefinementDefault
     780          24 :         real(RK)                                :: zplus1MaxDefault, maxRelativeErrorDefault, relerr
     781          24 :         real(RK)                                :: binaryMergerRateDensity, lookBackTimeRef
     782             : 
     783          21 :         nRefinementDefault = 7_IK; if (present(nRefinement)) nRefinementDefault = nRefinement
     784          24 :         zplus1MaxDefault = HUGE_RK; if (present(zplus1Max)) zplus1MaxDefault = zplus1Max
     785          24 :         maxRelativeErrorDefault = 1.e-6_RK; if (present(maxRelativeError)) maxRelativeErrorDefault = maxRelativeError
     786             : 
     787             :         lookBackTimeRef = getLookBackTime   ( zplus1 = zplus1 &
     788             :                                             , maxRelativeError = maxRelativeErrorDefault &
     789             :                                             , nRefinement = nRefinementDefault &
     790          24 :                                             )
     791             : 
     792             : #if defined OS_IS_WSL
     793             :         getMergerDelayTimePDF_WSL => getMergerDelayTimePDF
     794             :         getStarFormationRateDensity_WSL => getStarFormationRateDensity
     795             :         maxRelativeErrorDefault_WSL = maxRelativeErrorDefault
     796             :         nRefinementDefault_WSL = nRefinementDefault
     797             :         lookBackTimeRef_WSL = lookBackTimeRef
     798             : #endif
     799             : 
     800             :         call doQuadRombOpen ( getFunc           = getBinaryMergerRateDensityIntegrand   &
     801             :                             , integrate         = midexp                                &
     802             :                            !, integrate         = midinf                                &
     803             :                             , lowerLim          = zplus1                                &
     804             :                             , upperLim          = zplus1MaxDefault                      &
     805             :                             , maxRelativeError  = maxRelativeErrorDefault               &
     806             :                             , nRefinement       = nRefinementDefault                    &
     807             :                             , integral          = binaryMergerRateDensity               &
     808             :                             , relativeError     = relerr                                &
     809             :                             , numFuncEval       = neval                                 &
     810             :                             , ierr              = ierr                                  &
     811          24 :                             )
     812          24 :         if (ierr/=0_IK) then
     813             :             ! LCOV_EXCL_START
     814             :             write(output_unit,"(A)") ErrorMessage(ierr)
     815             :             error stop
     816             :             ! LCOV_EXCL_STOP
     817             :         end if
     818             : 
     819             : #if defined OS_IS_WSL
     820             :         nullify(getStarFormationRateDensity_WSL)
     821             :         nullify(getMergerDelayTimePDF_WSL)
     822             : #else
     823             :     contains
     824             : 
     825        3888 :         function getBinaryMergerRateDensityIntegrand(zplus1) result(binaryMergerRateIntegrand)
     826             : 
     827          24 :             use Cosmology_mod, only: getUniverseAgeDerivative
     828             :             implicit none
     829             :             real(RK)    , intent(in)    :: zplus1
     830             :             real(RK)                    :: binaryMergerRateIntegrand !,lognormpdf
     831        3888 :             real(RK)                    :: mergerDelayTime
     832             : 
     833             :             ! note that zp<z always, so that delay>0.
     834             :             mergerDelayTime = getLookBackTime   ( zplus1 = zplus1 &
     835             :                                                 , maxRelativeError = maxRelativeErrorDefault &
     836             :                                                 , nRefinement = nRefinementDefault &
     837        3888 :                                                 )
     838        3888 :             mergerDelayTime = mergerDelayTime - lookBackTimeRef
     839        3888 :             if (mergerDelayTime<=0._RK) then
     840             :                 ! LCOV_EXCL_START
     841             :                 write(output_unit,"(A)") "The mergerDelayTime is non-positive in getBinaryMergerRateDensityIntegrand(): (zplus1, mergerDelayTime) = ", zplus1, mergerDelayTime
     842             :                 error stop
     843             :                 ! LCOV_EXCL_STOP
     844             :             end if
     845             : 
     846             :             binaryMergerRateIntegrand   = getMergerDelayTimePDF(mergerDelayTime) &
     847             :                                         * getStarFormationRateDensity(zplus1) &
     848        3888 :                                         * getUniverseAgeDerivative(zplus1)
     849             : 
     850        3888 :         end function getBinaryMergerRateDensityIntegrand
     851             : #endif
     852             :     end function getBinaryMergerRateDensity
     853             : 
     854             : #if defined OS_IS_WSL
     855             :         ! This madness bypasses the Microsoft Subsystem for Linux Internal Function call GFortran Segmentation Fault error.
     856             :         function getBinaryMergerRateDensityIntegrand(zplus1) result(binaryMergerRateIntegrand)
     857             :             use, intrinsic :: iso_fortran_env, only: output_unit
     858             :             use Cosmology_mod, only: getUniverseAgeDerivative
     859             :             use Cosmology_mod, only: getLookBackTime
     860             :             implicit none
     861             :             real(RK)    , intent(in)    :: zplus1
     862             :             real(RK)                    :: binaryMergerRateIntegrand !,lognormpdf
     863             :             real(RK)                    :: mergerDelayTime
     864             :             ! note that zp<z always, so that delay>0.
     865             :             mergerDelayTime = getLookBackTime   ( zplus1 = zplus1 &
     866             :                                                 , maxRelativeError = maxRelativeErrorDefault_WSL &
     867             :                                                 , nRefinement = nRefinementDefault_WSL &
     868             :                                                 )
     869             :             mergerDelayTime = mergerDelayTime - lookBackTimeRef_WSL
     870             :             if (mergerDelayTime<=0._RK) then
     871             :             ! LCOV_EXCL_START
     872             :                 write(output_unit,"(A)") "The mergerDelayTime is non-positive in getBinaryMergerRateDensityIntegrand(): (zplus1, mergerDelayTime) = ", zplus1, mergerDelayTime
     873             :                 error stop
     874             :             end if
     875             :             ! LCOV_EXCL_START
     876             : 
     877             :             binaryMergerRateIntegrand   = getMergerDelayTimePDF_WSL(mergerDelayTime) &
     878             :                                         * getStarFormationRateDensity_WSL(zplus1) &
     879             :                                         * getUniverseAgeDerivative(zplus1)
     880             : 
     881             :         end function getBinaryMergerRateDensityIntegrand
     882             : #endif
     883             : 
     884             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     885             : 
     886             :     function getBinaryMergerRate( zplus1 &
     887             :                                 , zplus1Max &
     888             :                                 , nRefinement &
     889             :                                 , maxRelativeError &
     890             :                                 , getMergerDelayTimePDF &
     891             :                                 , getStarFormationRateDensity &
     892             :                                 ) result(binaryMergerRate)
     893             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     894             :         !DEC$ ATTRIBUTES DLLEXPORT :: getBinaryMergerRate
     895             : #endif
     896             :         use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE, getLogLumDisWicMpc
     897             :         use Constants_mod, only: RK, PI
     898             :         implicit none
     899             :         real(RK)    , intent(in)                :: zplus1
     900             :         real(RK)    , intent(in), optional      :: zplus1Max, maxRelativeError
     901             :         integer(IK) , intent(in), optional      :: nRefinement
     902             :         procedure(getRateDensity_proc)          :: getStarFormationRateDensity
     903             :         procedure(getMergerDelayTimePDF_proc)   :: getMergerDelayTimePDF
     904             :         real(RK)                                :: binaryMergerRate
     905             :         real(RK), parameter                     :: LOG_COEF = log(4._RK*PI*LS2HC)
     906             : 
     907             :         binaryMergerRate    = exp( LOG_COEF + 2_IK*getLogLumDisWicMpc(zplus1) - ( 3._RK*log(zplus1) + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) ) ) &
     908             :                             * getBinaryMergerRateDensity( zplus1 &
     909             :                                                         , zplus1Max &
     910             :                                                         , nRefinement &
     911             :                                                         , maxRelativeError &
     912             :                                                         , getMergerDelayTimePDF &
     913             :                                                         , getStarFormationRateDensity &
     914             :                                                         )
     915             : 
     916             :     end function getBinaryMergerRate
     917             : 
     918             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     919             : 
     920             : end module StarFormation_mod

ParaMonte: Plain Powerful Parallel Monte Carlo Library 
The Computational Data Science Lab
© Copyright 2012 - 2021