The ParaMonte Documentation Website
Current view: top level - kernel - Batse_mod.f90 (source / functions) Hit Total Coverage
Test: ParaMonte 1.5.1 :: MPI Parallel Kernel - Code Coverage Report Lines: 72 72 100.0 %
Date: 2021-01-08 13:07:16 Functions: 8 8 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
      44             : !> This module contains the procedures for modeling and manipulating the BATSE Catalog.
      45             : !> \author
      46             : !> Amir Shahmoradi, Tuesday April 30, 2019, 12:58 PM, SEIR, UTA
      47             : 
      48             : module Batse_mod
      49             : 
      50             :     use Constants_mod, only: IK, RK, PI, LN10
      51             : 
      52             :     implicit none
      53             : 
      54             :     character(*), parameter :: MODULE_NAME = "@Batse_mod"
      55             : 
      56             :     integer(IK) , parameter :: NLGRB = 1366_IK, NSGRB = 565_IK, NVAR = 4_IK
      57             :     real(RK)    , parameter :: MAX_LOG10PH53_4_LOGPBOLZERO  = 6.318167895318538_RK
      58             :     real(RK)    , parameter :: MIN_LOG10PH53_4_LOGPBOLZERO  = 4.92_RK
      59             :     real(RK)    , parameter :: MAX_LOGPH53_4_LOGPBOLZERO    = MAX_LOG10PH53_4_LOGPBOLZERO * LN10
      60             :     real(RK)    , parameter :: MIN_LOGPH53_4_LOGPBOLZERO    = MIN_LOG10PH53_4_LOGPBOLZERO * LN10
      61             :     real(RK)    , parameter :: DIF_LOGPH53_4_LOGPBOLZERO    = MAX_LOGPH53_4_LOGPBOLZERO - MIN_LOGPH53_4_LOGPBOLZERO
      62             : 
      63             :     ! Parameters for the effective peak flux of SGRBs (P64ms conversion to P1024):
      64             : 
      65             : !   ! The scale of the change in BATSE efficiency for different GRB durations.
      66             : !   real(RK)    , parameter :: THRESH_ERFC_BASE             = +0.146314238936889_RK
      67             : !
      68             : !   ! The scale of the change in BATSE efficiency for different GRB durations.
      69             : !   real(RK)    , parameter :: THRESH_ERFC_AMP              = +0.570285374263156_RK
      70             : !
      71             : !   ! Mean duration in the Error function used to model the connection between the peak fluxes in 64 and 1024 ms.
      72             : !   real(RK)    , parameter :: THRESH_ERFC_AVG              = -0.480811530417719_RK
      73             : !
      74             : !   ! scale of the duration in the Error function used to model the connection between the peak fluxes in 64 and 1024 ms.
      75             : !   real(RK)    , parameter :: THRESH_ERFC_STD              = +1.292443569094922_RK
      76             : 
      77             :     ! The scale of the change in BATSE efficiency for different GRB durations.
      78             :     real(RK)    , parameter :: THRESH_ERFC_BASE             = +0.146314238936889_RK * LN10
      79             : 
      80             :     ! The scale of the change in BATSE efficiency for different GRB durations.
      81             :     real(RK)    , parameter :: THRESH_ERFC_AMP              = +0.282313526464596_RK * LN10  ! Serfc
      82             : 
      83             :     ! Mean duration in the Error function used to model the connection between the peak fluxes in 64 and 1024 ms.
      84             :     real(RK)    , parameter :: THRESH_ERFC_AVG              = -0.483553339256463_RK * LN10  ! meandur
      85             : 
      86             :     ! scale of the duration in the Error function used to model the connection between the peak fluxes in 64 and 1024 ms.
      87             :     real(RK)    , parameter :: THRESH_ERFC_STD              = 1.0514698984694800_RK * LN10  ! scaledur
      88             : 
      89             :     ! inverse scale of the duration in the Error function used to model the connection between the peak fluxes in 64 and 1024 ms.
      90             :     real(RK)    , parameter :: THRESH_ERFC_STD_INV          = 1._RK / THRESH_ERFC_STD       ! inverse scaledur
      91             : 
      92             :     ! The height of the ERFC function.
      93             :     real(RK)    , parameter :: THRESH_ERFC_HEIGHT           = 2_IK * THRESH_ERFC_AMP
      94             : 
      95             :     ! correction that must be added to logPbol64ms to convert it to effective peak flux.
      96             :     ! Effective LogPbol limit above which trigger efficiency is 100%, for any Log(Epk) and Log(dur). It is equivalent to maximum Log(Pbol) at very long durations.
      97             :     real(RK)    , parameter :: THRESH_LOGPBOL64_CORRECTION  = DIF_LOGPH53_4_LOGPBOLZERO - MIN_LOGPH53_4_LOGPBOLZERO + THRESH_ERFC_HEIGHT ! equivalent to lpb_correction
      98             : 
      99             :     ! GRB attributes
     100             :     type :: Event_type
     101             :         real(RK) :: logPbol, logEpk, logSbol, logT90, logPF53
     102             :     end type Event_type
     103             : 
     104             :     type :: GRB_type
     105             :         integer(IK) :: count
     106             :         type(Event_type), allocatable :: Event(:)
     107             :         !type(Event_type) :: Event(NLGRB)
     108             :     end type GRB_type
     109             : 
     110             : #if defined CAF_ENABLED
     111             :     !> The `GRB_type` class containing GRB prompt attributes.
     112             :     type(GRB_type) :: GRB[*]
     113             : #else
     114             :     !> The `GRB_type` class containing GRB prompt attributes.
     115             :     type(GRB_type) :: GRB
     116             : #endif
     117             : 
     118             :     integer(IK), allocatable :: Trigger(:)
     119             :     !integer(IK) :: Trigger(NLGRB)
     120             :     !integer(IK) :: TriggerSGRB(NSGRB)
     121             : 
     122             :     interface getLogEffectivePeakPhotonFluxCorrection
     123             :         module procedure :: getLogEffectivePeakPhotonFluxCorrection_SPR
     124             :         module procedure :: getLogEffectivePeakPhotonFluxCorrection_DPR
     125             :     end interface getLogEffectivePeakPhotonFluxCorrection
     126             : 
     127             :     interface getLogEffectivePeakPhotonFlux
     128             :         module procedure :: getLogEffectivePeakPhotonFlux_SPR
     129             :         module procedure :: getLogEffectivePeakPhotonFlux_DPR
     130             :     end interface getLogEffectivePeakPhotonFlux
     131             : 
     132             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     133             : 
     134             : contains
     135             : 
     136             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     137             : 
     138             :     !> Return all log of data in natural (Neper) base.
     139             :     !>
     140             :     !> \param[in]   inFilePath  :   The path to the input BATSE file.
     141             :     !> \param[in]   outFilePath :   The path to the output BATSE file.
     142             :     !> \param[in]   isLgrb      :   A logical flag indicating what type of input file is being processed.
     143           6 :     subroutine readDataGRB(inFilePath,outFilePath,isLgrb)
     144             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     145             :         !DEC$ ATTRIBUTES DLLEXPORT :: readDataGRB
     146             : #endif
     147             : 
     148             :         use Parallelism_mod, only: Image_type
     149             :         use Constants_mod, only: IK, RK
     150             :         implicit none
     151             :         character(*), intent(in)            :: inFilePath, outFilePath
     152             :         integer(IK)                         :: inFileUnit, outFileUnit, igrb
     153             :         logical     , intent(in)            :: isLgrb
     154           6 :         type(Image_type)                    :: Image
     155             : 
     156           6 :         call Image%query()
     157             : 
     158           6 :         if (isLgrb) then
     159           3 :             GRB%count = NLGRB
     160             :         else
     161           3 :             GRB%count = NSGRB
     162             :         end if
     163             : 
     164           6 :         if (allocated(GRB%Event)) deallocate(GRB%Event); allocate(GRB%Event(GRB%count))
     165           6 :         if (allocated(Trigger)) deallocate(Trigger); allocate(Trigger(GRB%count))
     166             : 
     167             :         open( newunit = inFileUnit & ! LCOV_EXCL_LINE
     168             :             , file = inFilePath & ! LCOV_EXCL_LINE
     169             :             , status = "old" & ! LCOV_EXCL_LINE
     170             : #if defined INTEL_COMPILER_ENABLED && defined OS_IS_WINDOWS
     171             :             , SHARED & ! LCOV_EXCL_LINE
     172             : #endif
     173           6 :             )
     174             : 
     175           6 :         if (Image%isFirst) then
     176             :         ! LCOV_EXCL_START
     177             :             open(newunit=outFileUnit,file=outFilePath,status="replace")
     178             :             write(outFileUnit,"(9a30)"  ) "trigger"             &
     179             :                                         , "logPbol_1eV_20MeV"   &
     180             :                                         , "logSbol_1eV_20MeV"   &
     181             :                                         , "logEpk"              &
     182             :                                         , "logEPR1024"          &
     183             :                                         , "logEFR"              &
     184             :                                         , "logFPR1024"          &
     185             :                                         , "logT90"              &
     186             :                                         , "logEffPF53"
     187             :         end if
     188             :         ! LCOV_EXCL_STOP
     189             : 
     190             :         ! skip the header row in the input file
     191           6 :         read(inFileUnit,*)
     192             : 
     193             :         ! read BATSE GRB data
     194        5799 :         do igrb = 1, GRB%count
     195             : 
     196             :             read(inFileUnit,*   ) Trigger(igrb)             & ! LCOV_EXCL_LINE
     197        5793 :                                 , GRB%Event(igrb)%logPF53   &
     198        5793 :                                 , GRB%Event(igrb)%logEpk    &
     199        5793 :                                 , GRB%Event(igrb)%logSbol   &
     200       11586 :                                 , GRB%Event(igrb)%logT90
     201             : 
     202             :             ! convert all values to logarithm in base Neper
     203             : 
     204        5793 :             GRB%Event(igrb)%logPF53 = LN10 * GRB%Event(igrb)%logPF53
     205        5793 :             GRB%Event(igrb)%logEpk  = LN10 * GRB%Event(igrb)%logEpk
     206        5793 :             GRB%Event(igrb)%logSbol = LN10 * GRB%Event(igrb)%logSbol
     207        5793 :             GRB%Event(igrb)%logT90  = LN10 * GRB%Event(igrb)%logT90
     208             : 
     209             :             ! convert photon count data to energy in units of ergs
     210             : 
     211        5793 :             GRB%Event(igrb)%logPbol = getLogPbol( GRB%Event(igrb)%logEpk, GRB%Event(igrb)%logPF53 )
     212        5793 :             if (isLgrb) then
     213        4098 :                 GRB%Event(igrb)%logSbol = getLogPbol( GRB%Event(igrb)%logEpk, GRB%Event(igrb)%logSbol )
     214             :             else
     215        1695 :                 GRB%Event(igrb)%logPF53 = GRB%Event(igrb)%logPF53 - THRESH_ERFC_AMP * erfc( (GRB%Event(igrb)%logT90-THRESH_ERFC_AVG) * THRESH_ERFC_STD_INV )
     216             :             end if
     217             : 
     218             :             ! write the converted data to output file
     219             : 
     220        5799 :             if (Image%isFirst) then
     221             :             ! LCOV_EXCL_START
     222             :                 write(outFileUnit,"(I30,8E30.6)") Trigger(igrb)                                     &
     223             :                                                 , GRB%Event(igrb)%logPbol                           &
     224             :                                                 , GRB%Event(igrb)%logSbol                           &
     225             :                                                 , GRB%Event(igrb)%logEpk                            &
     226             :                                                 , GRB%Event(igrb)%logEpk-GRB%Event(igrb)%logPbol    &
     227             :                                                 , GRB%Event(igrb)%logEpk-GRB%Event(igrb)%logSbol    &
     228             :                                                 , GRB%Event(igrb)%logSbol-GRB%Event(igrb)%logPbol   &
     229             :                                                 , GRB%Event(igrb)%logT90                            &
     230             :                                                 , GRB%Event(igrb)%logPF53
     231             : 
     232             :             end if
     233             :             ! LCOV_EXCL_STOP
     234             : 
     235             :         end do
     236             : 
     237           6 :         if (Image%isFirst) close(outFileUnit)
     238             : 
     239           6 :         close(inFileUnit)
     240             : 
     241             : !#if defined CAF_ENABLED
     242             : !        sync images(*)
     243             : !
     244             : !    else
     245             : !
     246             : !        sync images(1)
     247             : !        do igrb = 1, GRB%count
     248             : !            GRB%Event(igrb)%logPbol = GRB[1]%Event(igrb)%logPbol
     249             : !            GRB%Event(igrb)%logSbol = GRB[1]%Event(igrb)%logSbol
     250             : !            GRB%Event(igrb)%logPF53 = GRB[1]%Event(igrb)%logPF53
     251             : !            GRB%Event(igrb)%logEpk  = GRB[1]%Event(igrb)%logEpk
     252             : !            GRB%Event(igrb)%logT90  = GRB[1]%Event(igrb)%logT90
     253             : !        end do
     254             : !
     255             : !    end if
     256             : !#endif
     257             : 
     258           6 :     end subroutine readDataGRB
     259             : 
     260             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     261             : 
     262         774 :     pure function getLog10PF53(log10epk,log10pbol) result(log10PF53)
     263             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     264             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLog10PF53
     265             : #endif
     266             :         ! Given Log10(Epk [KeV]) of an LGRB and its bolometric (0.0001-20000 KeV) peak flux [in units of Ergs/s], Log(Pbol),
     267             :         ! this function calculates the corresponding Log10(peak photon flux) in the BATSE detection energy range [50,300] KeV.
     268             :         ! Amir Shahmoradi, Wednesday June 27, 2012, 7:15 PM, IFS, The University of Texas at Austin.
     269             :         ! Amir Shahmoradi, Wednesday June 27, 2012, 9:29 PM, IFS, The University of Texas at Austin.
     270           6 :         use Constants_mod, only: RK
     271             :         implicit none
     272             :         real(RK), intent(in)    :: log10epk,log10pbol
     273             :         real(RK)                :: log10PF53
     274         774 :         if (log10epk<-2.915056638230699_RK) then
     275           9 :             log10PF53 = log10pbol    + 4.9200000000000000_RK
     276         765 :         elseif (log10epk<1.5_RK) then
     277             :             log10PF53 = log10pbol   + 5.7361200000000000_RK + log10epk * &
     278             :                                     ( 0.3093600000000000_RK + log10epk * &
     279             :                                     ( 0.0045600000000000_RK + log10epk * &
     280             :                                     ( 0.0015900000000000_RK + log10epk * &
     281             :                                     ( 0.0001533360000000_RK - log10epk   &
     282         396 :                                     * 0.0003574800000000_RK ))))
     283         369 :         elseif (log10epk<2.5_RK) then
     284             :             log10PF53 = log10pbol   + 1.9112800000000000_RK + log10epk * &
     285             :                                     ( 39.710390000000000_RK - log10epk * &
     286             :                                     ( 96.606280000000000_RK - log10epk * &
     287             :                                     ( 109.24696000000000_RK - log10epk * &
     288             :                                     ( 67.271800000000000_RK - log10epk * &
     289             :                                     ( 23.402390000000000_RK - log10epk * &
     290             :                                     ( 4.3454400000000000_RK - log10epk   &
     291          90 :                                     * 0.3360600000000000_RK ))))))
     292         279 :         elseif (log10epk<4._RK) then
     293             :             log10PF53 = log10pbol   + 2.8020600000000000_RK + log10epk * &
     294             :                                     ( 4.5690700000000000_RK - log10epk * &
     295             :                                     ( 1.9277200000000000_RK - log10epk * &
     296             :                                     ( 0.2938100000000000_RK - log10epk   &
     297         135 :                                     * 0.0148900000000000_RK )))
     298         144 :         elseif (log10epk<5.4093868613659435_RK) then
     299             :             log10PF53 = log10pbol   - 10.465330000000000_RK + log10epk * &
     300             :                                     ( 26.706370000000000_RK - log10epk * &
     301             :                                     ( 14.476310000000000_RK - log10epk * &
     302             :                                     ( 3.5404100000000000_RK - log10epk * &
     303             :                                     ( 0.4095700000000000_RK - log10epk   &
     304         135 :                                     * 0.0183100000000000_RK ))))
     305             :         else ! if (log10epk>=5.4093868613659435_RK) then
     306           9 :             log10PF53 = log10pbol   + 4.9200000000000000_RK
     307             :         end if
     308         774 :     end function getLog10PF53
     309             : 
     310             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     311             : 
     312       10407 :     pure function getLogPF53(logEpk,logPbol) result(logPF53)
     313             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     314             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogPF53
     315             : #endif
     316             :         ! Given Log(Epk [KeV]) of an LGRB and its bolometric (0.0001-20000 KeV) peak flux [in units of Ergs/s], Log(Pbol),
     317             :         ! this function calculates the corresponding Log(peak photon flux) in the BATSE detection energy range [50,300] KeV.
     318             :         ! Amir Shahmoradi, Wednesday June 27, 2012, 7:15 PM, IFS, The University of Texas at Austin.
     319             :         ! Amir Shahmoradi, Wednesday June 27, 2012, 9:29 PM, IFS, The University of Texas at Austin.
     320         774 :         use Constants_mod, only: RK
     321             :         implicit none
     322             :         real(RK), intent(in)    :: logEpk,logPbol
     323             :         real(RK)                :: logPF53
     324       10407 :         if (logEpk<-6.712165960423344_RK) then
     325           6 :             logPF53 = logPbol    + 11.32871865753070600_RK
     326       10401 :         elseif (logEpk<3.453877639491069_RK) then
     327             :             logPF53 = logPbol   + 13.20790440362500600_RK + logEpk * &
     328             :                                 ( 0.309360000000000000_RK + logEpk * &
     329             :                                 ( 0.001980382837478830_RK + logEpk * &
     330             :                                 ( 0.000299892598248466_RK + logEpk * &
     331             :                                 ( 1.25602147173493e-05_RK - logEpk   &
     332         330 :                                 * 1.27171265917873e-05_RK ))))
     333       10071 :         elseif (logEpk<5.756462732485115_RK) then
     334             :             logPF53 = logPbol   + 4.400884836537660000_RK + logEpk * &
     335             :                                 ( 39.71039000000000000_RK - logEpk * &
     336             :                                 ( 41.95557432120050000_RK - logEpk * &
     337             :                                 ( 20.60525451895990000_RK - logEpk * &
     338             :                                 ( 5.510436247342930000_RK - logEpk * &
     339             :                                 ( 0.832525333390336000_RK - logEpk * &
     340             :                                 ( 0.067135977132248900_RK - logEpk   &
     341        7362 :                                 * 0.002254876138523550_RK ))))))
     342        2709 :         elseif (logEpk<9.210340371976184_RK) then
     343             :             logPF53 = logPbol   + 6.451981585674900000_RK + logEpk * &
     344             :                                 ( 4.569070000000000000_RK - logEpk * &
     345             :                                 ( 0.837198158654537000_RK - logEpk * &
     346             :                                 ( 0.055416002698982300_RK - logEpk   &
     347        2613 :                                 * 0.001219684856402480_RK )))
     348          96 :         elseif (logEpk<12.455573549219071_RK) then
     349             :             logPF53 = logPbol   - 24.09731285126340000_RK + logEpk * &
     350             :                                 ( 26.70637000000000000_RK - logEpk * &
     351             :                                 ( 6.286981551320860000_RK - logEpk * &
     352             :                                 ( 0.667762738216888000_RK - logEpk * &
     353             :                                 ( 0.033549115287895400_RK - logEpk   &
     354          90 :                                 * 0.000651366755890191_RK ))))
     355             :         else
     356           6 :             logPF53 = logPbol   + 11.32871865753070600_RK
     357             :         end if
     358             :         !write(*,"(*(g0.13,:,', '))") "logEpk, logPbol, logPF53<0.0: ", logEpk, logPbol, logPF53
     359       10407 :     end function getLogPF53
     360             : 
     361             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     362             : 
     363       10149 :     pure function getLogPbol(logEpk,logPF53) result(logPbol)
     364             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     365             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogPbol
     366             : #endif
     367             :         ! Given Log10(Epk [KeV]) of an LGRB and its bolometric (0.0001-20000 KeV) peak flux [in units of Ergs/s], Log(Pbol),
     368             :         ! this function calculates the corresponding Log10(peak photon flux) in the BATSE detection energy range [50,300] KeV.
     369             :         ! Amir Shahmoradi, Wednesday June 27, 2012, 7:15 PM, IFS, The University of Texas at Austin.
     370             :         ! Amir Shahmoradi, Wednesday June 27, 2012, 9:29 PM, IFS, The University of Texas at Austin.
     371       10407 :         use Constants_mod, only: RK
     372             :         implicit none
     373             :         real(RK), intent(in)    :: logEpk, logPF53
     374             :         real(RK)                :: logPbol
     375       10149 :         logPbol = logPF53 - getLogPF53(logEpk,0._RK)
     376       20298 :     end function getLogPbol
     377             : 
     378             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     379             : 
     380           3 :     pure function getLogEffectivePeakPhotonFlux_SPR(logPeakPhotonFlux64ms,logT90) result(logEffectivePeakPhotonFlux)
     381             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     382             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogEffectivePeakPhotonFlux_SPR
     383             : #endif
     384             :         ! Converts an input natural-log peak photon flux in 64ms timescale to an effective triggering peak photon flux.
     385             :         ! To do so, the observed T90 duration of the event is also necessary as input.
     386             :         ! Reference: Eqn A4 of Shahmoradi and Nemiroff 2015, MNRAS, Short versus long gamma-ray bursts.
     387             :         use, intrinsic :: iso_fortran_env, only: RK => real32
     388             :         implicit none
     389             :         real(RK), intent(in)    :: logPeakPhotonFlux64ms, logT90
     390             :         real(RK)                :: logEffectivePeakPhotonFlux
     391           3 :         logEffectivePeakPhotonFlux  = logPeakPhotonFlux64ms - getLogEffectivePeakPhotonFluxCorrection(logT90)
     392       10149 :     end function getLogEffectivePeakPhotonFlux_SPR
     393             : 
     394             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     395             : 
     396           3 :     pure function getLogEffectivePeakPhotonFlux_DPR(logPeakPhotonFlux64ms,logT90) result(logEffectivePeakPhotonFlux)
     397             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     398             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogEffectivePeakPhotonFlux_DPR
     399             : #endif
     400             :         ! Converts an input natural-log peak photon flux in 64ms timescale to an effective triggering peak photon flux.
     401             :         ! To do so, the observed T90 duration of the event is also necessary as input.
     402             :         ! Reference: Eqn A4 of Shahmoradi and Nemiroff 2015, MNRAS, Short versus long gamma-ray bursts.
     403             :         use, intrinsic :: iso_fortran_env, only: RK => real64
     404             :         implicit none
     405             :         real(RK), intent(in)    :: logPeakPhotonFlux64ms, logT90
     406             :         real(RK)                :: logEffectivePeakPhotonFlux
     407           3 :         logEffectivePeakPhotonFlux  = logPeakPhotonFlux64ms - getLogEffectivePeakPhotonFluxCorrection(logT90)
     408           3 :     end function getLogEffectivePeakPhotonFlux_DPR
     409             : 
     410             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     411             : 
     412           6 :     pure function getLogEffectivePeakPhotonFluxCorrection_SPR(logT90) result(logEffectivePeakPhotonFluxCorrection)
     413             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     414             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogEffectivePeakPhotonFluxCorrection_SPR
     415             : #endif
     416             :         ! Converts an input natural-log peak photon flux in 64ms timescale to an effective triggering peak photon flux.
     417             :         ! To do so, the observed T90 duration of the event is also necessary as input.
     418             :         ! Reference: Eqn A4 of Shahmoradi and Nemiroff 2015, MNRAS, Short versus long gamma-ray bursts.
     419             :         use, intrinsic :: iso_fortran_env, only: RK => real32
     420             :         implicit none
     421             :         real(RK), intent(in)    :: logT90
     422             :         real(RK)                :: logEffectivePeakPhotonFluxCorrection
     423           6 :         logEffectivePeakPhotonFluxCorrection    = real(THRESH_ERFC_AMP,RK) * erfc(real((logT90-THRESH_ERFC_AVG)/THRESH_ERFC_STD,kind=RK))
     424             :                                               ! + THRESH_ERFC_BASE ! adding this term will make the effective peak flux equivalent to PF1024ms
     425           3 :     end function getLogEffectivePeakPhotonFluxCorrection_SPR
     426             : 
     427             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     428             : 
     429           6 :     pure function getLogEffectivePeakPhotonFluxCorrection_DPR(logT90) result(logEffectivePeakPhotonFluxCorrection)
     430             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     431             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogEffectivePeakPhotonFluxCorrection_DPR
     432             : #endif
     433             :         ! Converts an input natural-log peak photon flux in 64ms timescale to an effective triggering peak photon flux.
     434             :         ! To do so, the observed T90 duration of the event is also necessary as input.
     435             :         ! Reference: Eqn A4 of Shahmoradi and Nemiroff 2015, MNRAS, Short versus long gamma-ray bursts.
     436             :         use, intrinsic :: iso_fortran_env, only: RK => real64
     437             :         implicit none
     438             :         real(RK), intent(in)    :: logT90
     439             :         real(RK)                :: logEffectivePeakPhotonFluxCorrection
     440           6 :         logEffectivePeakPhotonFluxCorrection    = THRESH_ERFC_AMP * erfc(real((logT90-THRESH_ERFC_AVG)/THRESH_ERFC_STD,kind=RK))
     441             :                                               ! + THRESH_ERFC_BASE ! adding this term will make the effective peak flux equivalent to PF1024ms
     442           6 :     end function getLogEffectivePeakPhotonFluxCorrection_DPR
     443             : 
     444             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     445             : 
     446             : end module Batse_mod ! LCOV_EXCL_LINE

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