1
!! Copyright 2009, Andrew Benson <abenson@caltech.edu>
3
!! This file is part of Galacticus.
5
!! Galacticus is free software: you can redistribute it and/or modify
6
!! it under the terms of the GNU General Public License as published by
7
!! the Free Software Foundation, either version 3 of the License, or
8
!! (at your option) any later version.
10
!! Galacticus is distributed in the hope that it will be useful,
11
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
12
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13
!! GNU General Public License for more details.
15
!! You should have received a copy of the GNU General Public License
16
!! along with Galacticus. If not, see <http://www.gnu.org/licenses/>.
21
module CDM_Primordial_Power_Spectrum
22
use ISO_Varying_String
25
public :: Primordial_Power_Spectrum_CDM, CDM_Primordial_Power_Spectrum_State_Retrieve
27
! Flag to indicate if this module has been initialized.
28
logical :: powerSpectrumInitialized=.false., tablesInitialized=.false.
30
! Variables to hold the tabulated power spectrum data.
31
integer :: powerSpectrumNumberPoints=-1
32
double precision, allocatable, dimension(:) :: powerSpectrumLogWavenumber,powerSpectrumLogP
33
type(fgsl_interp) :: interpolationObject
34
type(fgsl_interp_accel) :: interpolationAccelerator
35
logical :: resetInterpolation=.true.
37
! Name of power spectrum method used.
38
type(varying_string) :: powerSpectrumMethod
40
! Pointer to the subroutine that tabulates the transfer function and template interface for that subroutine.
41
procedure(Power_Spectrum_Tabulate_Template), pointer :: Power_Spectrum_Tabulate => null()
43
subroutine Power_Spectrum_Tabulate_Template(logWavenumber,powerSpectrumNumberPoints,powerSpectrumLogWavenumber &
45
double precision, intent(in) :: logWavenumber
46
double precision, allocatable, dimension(:), intent(inout) :: powerSpectrumLogWavenumber,powerSpectrumLogP
47
integer, intent(out) :: powerSpectrumNumberPoints
48
end subroutine Power_Spectrum_Tabulate_Template
53
double precision function Primordial_Power_Spectrum_CDM(wavenumber)
54
!% Return the CDM primordial power spectrum for $k=${\tt wavenumber} [Mpc$^{-1}$].
55
use Numerical_Interpolation
57
double precision, intent(in) :: wavenumber
58
double precision :: logWavenumber
60
! Get logarithm of wavenumber.
61
logWavenumber=dlog(wavenumber)
63
!$omp critical(Power_Spectrum_Initialization)
64
! Initialize if necessary.
65
if (.not.(powerSpectrumInitialized.and.tablesInitialized)) then
66
call Power_Spectrum_Initialize(logWavenumber)
67
call Interpolate_Done(interpolationObject,interpolationAccelerator,resetInterpolation)
68
resetInterpolation=.true.
71
! If wavenumber is out of range, attempt to remake the table.
72
if (logWavenumber<powerSpectrumLogWavenumber(1) .or. logWavenumber>powerSpectrumLogWavenumber(powerSpectrumNumberPoints) )&
74
call Power_Spectrum_Tabulate(logWavenumber,powerSpectrumNumberPoints,powerSpectrumLogWavenumber,powerSpectrumLogP)
75
call Interpolate_Done(interpolationObject,interpolationAccelerator,resetInterpolation)
76
resetInterpolation=.true.
78
!$omp end critical(Power_Spectrum_Initialization)
80
! Interpolate in the tabulated function and return a value.
81
Primordial_Power_Spectrum_CDM=dexp(Interpolate(powerSpectrumNumberPoints,powerSpectrumLogWavenumber,powerSpectrumLogP &
82
&,interpolationObject,interpolationAccelerator,logWavenumber,reset=resetInterpolation))
85
end function Primordial_Power_Spectrum_CDM
87
subroutine Power_Spectrum_Initialize(logWavenumber)
88
!% Initializes the transfer function module.
91
!# <include directive="powerSpectrumMethod" type="moduleUse">
92
include 'structure_formation.CDM.power_spectrum.primordial.modules.inc'
95
double precision, intent(in) :: logWavenumber
97
if (.not.powerSpectrumInitialized) then
98
! Get the primordial power spectrum method parameter.
100
!@ <name>powerSpectrumMethod</name>
101
!@ <defaultValue>power law</defaultValue>
102
!@ <attachedTo>module</attachedTo>
104
!@ The name of the method to be used for computing the primordial power spectrum.
107
call Get_Input_Parameter('powerSpectrumMethod',powerSpectrumMethod,defaultValue='power law')
108
! Include file that makes calls to all available method initialization routines.
109
!# <include directive="powerSpectrumMethod" type="code" action="subroutine">
110
!# <subroutineArgs>powerSpectrumMethod,Power_Spectrum_Tabulate</subroutineArgs>
111
include 'structure_formation.CDM.power_spectrum.inc'
113
if (.not.associated(Power_Spectrum_Tabulate)) call Galacticus_Error_Report('Power_Spectrum_Initialize','method '&
114
&//char(powerSpectrumMethod)//' is unrecognized')
116
! Call routine to initialize the transfer function.
117
call Power_Spectrum_Tabulate(logWavenumber,powerSpectrumNumberPoints,powerSpectrumLogWavenumber,powerSpectrumLogP)
118
tablesInitialized=.true.
119
! Flag that the module is now initialized.
120
powerSpectrumInitialized=.true.
122
end subroutine Power_Spectrum_Initialize
124
!# <galacticusStateRetrieveTask>
125
!# <unitName>CDM_Primordial_Power_Spectrum_State_Retrieve</unitName>
126
!# </galacticusStateRetrieveTask>
127
subroutine CDM_Primordial_Power_Spectrum_State_Retrieve(stateFile,fgslStateFile)
128
!% Reset the tabulation if state is to be retrieved. This will force tables to be rebuilt.
129
use Memory_Management
131
integer, intent(in) :: stateFile
132
type(fgsl_file), intent(in) :: fgslStateFile
134
powerSpectrumNumberPoints=0
135
if (allocated(powerSpectrumLogWavenumber)) call Dealloc_Array(powerSpectrumLogWavenumber)
136
if (allocated(powerSpectrumLogP )) call Dealloc_Array(powerSpectrumLogP )
137
tablesInitialized=.false.
139
end subroutine CDM_Primordial_Power_Spectrum_State_Retrieve
141
end module CDM_Primordial_Power_Spectrum
added
removed
source/structure_formation.CDM.power_spectrum.primordial.F90
