diff --git a/doc/Coding.tex b/doc/Coding.tex index 2d42cbf03d..15241bfab6 100644 --- a/doc/Coding.tex +++ b/doc/Coding.tex @@ -570,7 +570,7 @@ \subsection{Finding Roots of Equations}\index{root finding}\index{numerical algo \end{lstlisting} The above example begins by importing the {\normalfont \ttfamily Root\_Finder} module and then creating a {\normalfont \ttfamily rootFinder} object called {\normalfont \ttfamily finder}. This is made OpenMP {\normalfont \ttfamily threadprivate} so that it may be used simultaneously by all threads. The first step is to initialize {\normalfont \ttfamily finder}---the {\normalfont \ttfamily isInitialized} method tells us if this has already happened. The most important step is to specify the function that will evaluate $f(x)$. This is done via the {\normalfont \ttfamily rootFunction} method---once done, the {\normalfont \ttfamily rootFinder} object is marked as initialized (and the {\normalfont \ttfamily isInitialized} method will return {\normalfont \ttfamily true}). All other initialization steps are optional. In this example, we use the {\normalfont \ttfamily type} method to specify that the {\normalfont \ttfamily Brent} algorithm should be used for root finding. Any valid \href{http://www.gnu.org/software/gsl/manual/html_node/Root-Bracketing-Algorithms.html}{GSL-supported root finding algorithm} can be used. We then use the {\normalfont \ttfamily tolerance} method to specify both the absolute and relative tolerances in the $x$ variable that must be attained to declare the root to be found. Both arguments are optional---default values of $10^{-10}$ will be used if either tolerance is not specified. -The final step of initialization is to call the {\normalfont \ttfamily rangeExpand} method. This specifies how the initial guessed value or range for $x$ should be expanded to bracket the root. If you plan to always specify an initial range, and know that it will always bracket the root, you do not need to specify how the range should be expanded. In this case we've specified that range expansion is multiplicative---that is, the lower and upper values of $x$ defining the range will be multiplied by fixed factors until the root is bracketed---via the {\normalfont \ttfamily rangeExpandType=rangeExpandMultiplicative} option. Alternatively, additive expansion is possible using {\normalfont \ttfamily rangeExpandType=rangeExpandAdditive}. The factors by which to multiply the lower and upper bounds of the range (or the factor to add in the case of additive expansion) are specified by the {\normalfont \ttfamily rangeExpandDownward} and {\normalfont \ttfamily rangeExpandUpward} options. Iit is possible to specify absolute lower/upper limits to the range via the {\normalfont \ttfamily rangeDownwardLimit} and {\normalfont \ttfamily rangeUpwardLimit} options. The range will not be expanded beyond these limits---if the root cannot be bracketed without exceeding these limits an error condition will occur. Finally, it is possible to indicate the expected sign of $f(x)$ at the lower and/or upper limits via the {\normalfont \ttfamily rangeExpandDownwardSignExpect} and {\normalfont \ttfamily rangeExpandUpwardSignExpect} options. Valid settings are {\normalfont \ttfamily rangeExpandSignExpectNegative}, {\normalfont \ttfamily rangeExpandSignExpectPositive}, and {\normalfont \ttfamily rangeExpandSignExpectNone} (the default---implying that there is no expectation for the sign). If the sign of $f(x)$ is specified, then range expansion will stop once the expected sign is found. This can often improve efficiency, by allowing the range expander to expand the range in only one direction, resulting in a narrower range in which to search for the root. +The final step of initialization is to call the {\normalfont \ttfamily rangeExpand} method. This specifies how the initial guessed value or range for $x$ should be expanded to bracket the root. If you plan to always specify an initial range, and know that it will always bracket the root, you do not need to specify how the range should be expanded. In this case we've specified that range expansion is multiplicative---that is, the lower and upper values of $x$ defining the range will be multiplied by fixed factors until the root is bracketed---via the {\normalfont \ttfamily rangeExpandType=rangeExpandMultiplicative} option. Alternatively, additive expansion is possible using {\normalfont \ttfamily rangeExpandType=rangeExpandAdditive}. The factors by which to multiply the lower and upper bounds of the range (or the factor to add in the case of additive expansion) are specified by the {\normalfont \ttfamily rangeExpandDownward} and {\normalfont \ttfamily rangeExpandUpward} options. It is possible to specify absolute lower/upper limits to the range via the {\normalfont \ttfamily rangeDownwardLimit} and {\normalfont \ttfamily rangeUpwardLimit} options. The range will not be expanded beyond these limits---if the root cannot be bracketed without exceeding these limits an error condition will occur. Finally, it is possible to indicate the expected sign of $f(x)$ at the lower and/or upper limits via the {\normalfont \ttfamily rangeExpandDownwardSignExpect} and {\normalfont \ttfamily rangeExpandUpwardSignExpect} options. Valid settings are {\normalfont \ttfamily rangeExpandSignExpectNegative}, {\normalfont \ttfamily rangeExpandSignExpectPositive}, and {\normalfont \ttfamily rangeExpandSignExpectNone} (the default---implying that there is no expectation for the sign). If the sign of $f(x)$ is specified, then range expansion will stop once the expected sign is found. This can often improve efficiency, by allowing the range expander to expand the range in only one direction, resulting in a narrower range in which to search for the root. Finally, we use the {\normalfont \ttfamily find} method to return the value of the root. The first argument to {\normalfont \ttfamily find} is the name of the function that evaluates $f(x)$. Additionally, we must supply either {\normalfont \ttfamily rootGuess} (a scalar value guess to use as the initial value for both the lower and upper values of the range---note that range expansion must be allowed in this case), or {\normalfont \ttfamily rootRange} (a two-element array to use as the initial lower and upper values of the range bracketing the root). @@ -580,7 +580,7 @@ \section{Computation Dependencies and Data Files}\label{sec:codeUniqueLabels}\in In many situations, some module in \glc\ might want to perform a calculation and then store the results to a file so that they can be reused later. A good example is the \refPhysics{transferFunctionCAMB} transfer function model, which computes a transfer function using {\normalfont \scshape CAMB} and stores this function in a file so that it can be re-read next time, avoiding the need to recompute the transfer function. A problem arises in such cases as the calculation may depend on the values of parameters (in our example, the transfer function will depend on cosmological parameters for example). We would like to record which parameter values this calculation refers to, perhaps encoding these into the file name, so that we can reuse these data in a future run only if the parameter values are unchanged. Given the modular nature of \glc\ it is impossible to know in advance which parameters will be relevant (e.g. does the cosmological parameter implementation have a parameter that describes a time varying equation of state for dark energy?). -To addres this problem, \glc\ provides a mechanism to generate a unique descriptor for a given object. This descriptor encodes the parameter used to construct the object, and recursively includes the parameters used to construct any other object which is composited. A long-form (human readable) descriptor is returned by the {\normalfont \ttfamily descriptor} method associated with all {\normalfont \ttfamily functionClass} objects. Additionally, the {\normalfont \ttfamily hashedDescriptor} method will return an MD5 hash of the descriptor, which will be unique (up to collisions) and can be used to idetify the object both internally and, for example, when used as a suffix to file names. If the optional {\normalfont \ttfamily includeSourceDigest} argument is set to true in the {\normalfont \ttfamily hashedDescriptor} method then the hashed descriptor will include a hash of the source code of the object (and all composited objects) such that the descriptor will change should the source code be changed. +To address this problem, \glc\ provides a mechanism to generate a unique descriptor for a given object. This descriptor encodes the parameter used to construct the object, and recursively includes the parameters used to construct any other object which is composited. A long-form (human readable) descriptor is returned by the {\normalfont \ttfamily descriptor} method associated with all {\normalfont \ttfamily functionClass} objects. Additionally, the {\normalfont \ttfamily hashedDescriptor} method will return an MD5 hash of the descriptor, which will be unique (up to collisions) and can be used to identify the object both internally and, for example, when used as a suffix to file names. If the optional {\normalfont \ttfamily includeSourceDigest} argument is set to true in the {\normalfont \ttfamily hashedDescriptor} method then the hashed descriptor will include a hash of the source code of the object (and all composited objects) such that the descriptor will change should the source code be changed. \section{Optimization}\label{sec:Optimization}\index{optimization} @@ -626,6 +626,21 @@ \section{Enumerations} \input{autoEnumerationDefinitions} +\section{Defined Constants} + +\glc\ defines numerous constants, including mathematical constants (e.g. $\pi$), physical constants (e.g. the speed of light), unit conversions (e.g. Angstroms to meters), and prefixes (e.g. ``kilo'', ``mega'', etc.). These should be used whenever a constant is needed in the code---it is bad practice to use the numerical value of a constant directly in the code\footnote{Both because it is prone to mistakes (the more times a numerical value is used directly, the more chances there are for typos and other errors), and because using a named constant makes it much easier to understand \emph{what} the code is doing and \emph{why}.}. + +All defined constants are described below, along with references to their source, and the name of the module in which the constant is defined. To import a constant into a function, you would add a {\normalfont \ttfamily use} statement. For example, for the constant $\pi$, you would use: +\begin{verbatim} + use :: Numerical_Constants_Math, only : Pi +\end{verbatim} +after which you can use this constant, e.g.: +\begin{verbatim} + volumeSphere=4.0d0*Pi/3.0d0*radiusSphere**3 +\end{verbatim} + +\input{constants} + \section{Classes} The return type of each method, and the interfaces (i.e. the types and names of its arguments) are specified for each method of each object. A ``{\normalfont \ttfamily void}'' return type indicates a subroutine. @@ -696,5 +711,5 @@ \subsection{Methods available to all {\normalfont \ttfamily functionClass}es}\la \end{itemize} \end{description} - + \input{dataMethods} diff --git a/doc/Galacticus.bib b/doc/Galacticus.bib index a9002192f9..23be184c2f 100644 --- a/doc/Galacticus.bib +++ b/doc/Galacticus.bib @@ -3665,6 +3665,19 @@ @article{inoue_updated_2014 adsnote = {Provided by the SAO/NASA Astrophysics Data System} } +@book{jackson_classical_1999, + address = {New York}, + edition = {3. ed}, + title = {Classical electrodynamics}, + isbn = {978-0-471-30932-1}, + language = {eng}, + publisher = {Wiley}, + author = {Jackson, John David}, + year = {1999}, + note = {OCLC: 925677836}, + keywords = {Electrodynamics, Électrodynamique, Elektrodynamik}, +} + @article{jassal_wmap_2005, title = {{WMAP} constraints on low redshift evolution of dark energy}, volume = {356}, @@ -5755,6 +5768,25 @@ @article{obreschkow_simulation_2009 pages = {1467--1484} } +@article{oke_secondary_1983, + title = {Secondary standard stars for absolute spectrophotometry.}, + volume = {266}, + issn = {0004-637X}, + url = {https://ui.adsabs.harvard.edu/abs/1983ApJ...266..713O}, + doi = {10.1086/160817}, + abstract = {Based on an adopted absolute spectral energy distribution for the primary standard star Alpha Lyrae, absolute fluxes are given for the four very metal-deficient F type subdwarfs HD 19445, HD 84937, BD + 26.2606 deg, and BD + 17.4703 deg. Somewhat inferior data are also given for HD 140283. The data are given for 40-A bands and cover the wavelength range from 3080 A to 12,000 A. The four stars, all near magnitude 9 and distributed around the sky, are intended as secondary standards for absolute spectrophotometry.}, + urldate = {2025-01-29}, + journal = {The Astrophysical Journal}, + author = {Oke, J. B. and Gunn, J. E.}, + month = mar, + year = {1983}, + note = {Publisher: IOP +ADS Bibcode: 1983ApJ...266..713O}, + keywords = {Astronomy, Calibrating, Line Spectra, Reference Stars, Spectral Energy Distribution, Stellar Spectrophotometry, Subdwarf Stars}, + pages = {713--717}, + file = {Full Text PDF:/home/abensonca/.mozilla/firefox/f54gqgdx.default/zotero/storage/FBTQED33/Oke and Gunn - 1983 - Secondary standard stars for absolute spectrophoto.pdf:application/pdf}, +} + @article{ondaro-mallea_non-universality_2022, title = {Non-universality of the mass function: dependence on the growth rate and power spectrum shape}, volume = {509}, diff --git a/perl/Galacticus/Build/SourceTree.pm b/perl/Galacticus/Build/SourceTree.pm index 32e3eec1cd..94ac8446e3 100644 --- a/perl/Galacticus/Build/SourceTree.pm +++ b/perl/Galacticus/Build/SourceTree.pm @@ -41,7 +41,7 @@ use Galacticus::Build::SourceTree::Process::DebugHDF5; use Galacticus::Build::SourceTree::Process::DebugMPI; use Galacticus::Build::SourceTree::Process::ProfileOpenMP; use Galacticus::Build::SourceTree::Process::ThreadSafeIO; -use Galacticus::Build::SourceTree::Process::GSLConstants; +use Galacticus::Build::SourceTree::Process::Constants; use Galacticus::Build::SourceTree::Process::HDF5FCInterop; use Galacticus::Build::SourceTree::Process::Constructors; use Galacticus::Build::SourceTree::Process::ConditionalCall; diff --git a/perl/Galacticus/Build/SourceTree/Process/Constants.pm b/perl/Galacticus/Build/SourceTree/Process/Constants.pm new file mode 100644 index 0000000000..676e820f45 --- /dev/null +++ b/perl/Galacticus/Build/SourceTree/Process/Constants.pm @@ -0,0 +1,145 @@ +# Contains a Perl module which implements physical/numerical constants. + +package Galacticus::Build::SourceTree::Process::Constants; +use strict; +use warnings; +use utf8; +use Cwd; +use lib $ENV{'GALACTICUS_EXEC_PATH'}."/perl"; +use File::Temp qw/ tempfile /; +use XML::Simple; + +# Insert hooks for our functions. +$Galacticus::Build::SourceTree::Hooks::processHooks{'constant'} = \&Process_Constant; + +sub Process_Constant { + # Get the tree. + my $tree = shift(); + # Initialize structure to store lists of constants. + my $constants; + my @constantsOrphaned; + my $fileName; + # Walk the tree. + my $node = $tree; + my $depth = 0; + while ( $node ) { + if ( $node->{'type'} eq "file" ) { + $fileName = $node->{'name'}; + } elsif ( $node->{'type'} eq "module" ) { + my $moduleName = $node->{'name'}; + if ( @constantsOrphaned ) { + foreach my $constant ( @constantsOrphaned ) { + $constant->{'module'} = $moduleName; + } + push(@{$constants->{'constant'}},@constantsOrphaned); + undef(@constantsOrphaned); + } + } elsif ( $node->{'type'} eq "constant" && ! $node->{'directive'}->{'processed'} ) { + my $code; + my $type = exists($node->{'directive'}->{'type'}) ? $node->{'directive'}->{'type'} : "double precision"; + if ( exists($node->{'directive'}->{'gslSymbol'}) ) { + # Constant value is to be extracted from GSL. + # Generate, compile, and execute a minimal C code which simply outputs the relevant GSL constant. + if ( $type eq "enum" ) { + # An enum type. + my @members = split(/\s*,\s*/,$node->{'directive'}->{'members'}); + (undef, my $tmpSourceFileName) = tempfile('tempXXXXX', SUFFIX => '.c', OPEN => 0, DIR => $ENV{'BUILDPATH'}); + (undef, my $tmpExecFileName ) = tempfile('tempXXXXX', SUFFIX => '.x', OPEN => 0, DIR => $ENV{'BUILDPATH'}); + open(my $tmpSource,">".$tmpSourceFileName); + print $tmpSource "\#include \n"; + print $tmpSource "\#include \n"; + print $tmpSource "\#include {'directive'}->{'gslHeader'}.".h>\n"; + print $tmpSource "int main () {\n"; + print $tmpSource " printf(\"".join(" ",("%i") x scalar(@members))."\", ".$node->{'directive'}->{'members'}.");\n"; + print $tmpSource "}\n"; + close($tmpSource); + system($ENV{'CCOMPILER'}." -o ".$tmpExecFileName." ".$tmpSourceFileName." ".$ENV{'CFLAGS'}); + die("Galacticus::Build::SourceTree::Process::GSLConstants: failed to compile") + unless ( $? == 0 ); + open(my $gslPipe,$tmpExecFileName."|"); + my $gslEnum = <$gslPipe>; + close($gslPipe); + unlink($tmpSourceFileName,$tmpExecFileName); + my @values = split(" ",$gslEnum); + my @definitions; + for(my $i=0;$i '.c', OPEN => 0, DIR => $ENV{'BUILDPATH'}); + (undef, my $tmpExecFileName ) = tempfile('tempXXXXX', SUFFIX => '.x', OPEN => 0, DIR => $ENV{'BUILDPATH'}); + open(my $tmpSource,">".$tmpSourceFileName); + print $tmpSource "\#include \n"; + print $tmpSource "\#include \n"; + print $tmpSource "\#include {'directive'}->{'gslHeader'}.".h>\n"; + if ( $type eq "double precision" ) { + print $tmpSource "\#ifdef LDBL_DECIMAL_DIG\n"; + print $tmpSource " \#define OP_LDBL_Digs (LDBL_DECIMAL_DIG)\n"; + print $tmpSource "\#else \n"; + print $tmpSource " \#ifdef DECIMAL_DIG\n"; + print $tmpSource " \#define OP_LDBL_Digs (DECIMAL_DIG)\n"; + print $tmpSource " \#else \n"; + print $tmpSource " \#define OP_LDBL_Digs (LDBL_DIG + 3)\n"; + print $tmpSource " \#endif\n"; + print $tmpSource "\#endif\n"; + print $tmpSource "int main () {\n"; + print $tmpSource " printf(\"%.*e\", OP_LDBL_Digs - 1, ".$node->{'directive'}->{'gslSymbol'}.");\n"; + print $tmpSource "}\n"; + } elsif ( $type eq "integer" ) { + print $tmpSource "int main () {\n"; + print $tmpSource " printf(\"%i\", ".$node->{'directive'}->{'gslSymbol'}.");\n"; + print $tmpSource "}\n"; + } + close($tmpSource); + system($ENV{'CCOMPILER'}." -o ".$tmpExecFileName." ".$tmpSourceFileName." ".$ENV{'CFLAGS'}); + die("Galacticus::Build::SourceTree::Process::GSLConstants: failed to compile") + unless ( $? == 0 ); + open(my $gslPipe,$tmpExecFileName."|"); + my $gslConstant = <$gslPipe>; + close($gslPipe); + unlink($tmpSourceFileName,$tmpExecFileName); + $gslConstant =~ s/e/d/; + $node->{'directive'}->{'value'} = $gslConstant; + $code = $type.", parameter, public :: ".$node->{'directive'}->{'variable'}."=".$gslConstant."\n"; + } + } else { + # Constant value is defined internally. + $code = $type.", parameter, public :: ".$node->{'directive'}->{'variable'}."=".$node->{'directive'}->{'value'}."\n"; + } + # Insert a new code node to define the constant. + my $newNode = + { + type => "code", + content => $code, + firstChild => undef() + }; + &Galacticus::Build::SourceTree::InsertAfterNode($node,[$newNode]); + # Accumulate to a list of constants from this module. + push(@constantsOrphaned,$node->{'directive'}); + # Mark the directive as processed. + $node->{'directive'}->{'processed'} = 1; + } + # Walk to the next node in the tree. + $node = &Galacticus::Build::SourceTree::Walk_Tree($node,\$depth); + } + # Add file name attribute. + foreach my $constant ( @{$constants->{'constant'}} ) { + $constant->{'fileName'} = $fileName; + } + # If building documentation, output accumulated constants to file. + if ( defined($constants) && exists($ENV{'GALACTICUS_BUILD_DOCS'}) && $ENV{'GALACTICUS_BUILD_DOCS'} eq "yes" ) { + my $xml = new XML::Simple(); + (my $constantsName = $fileName) =~ s/\.F90$/.constants.xml/; + open(my $constantsFile,">",$ENV{'BUILDPATH'}."/".$constantsName); + print $constantsFile $xml->XMLout($constants,RootName => "constants"); + close($constantsFile); + } +} + +1; diff --git a/perl/Galacticus/Build/SourceTree/Process/GSLConstants.pm b/perl/Galacticus/Build/SourceTree/Process/GSLConstants.pm deleted file mode 100644 index 5b2a9fa58b..0000000000 --- a/perl/Galacticus/Build/SourceTree/Process/GSLConstants.pm +++ /dev/null @@ -1,107 +0,0 @@ -# Contains a Perl module which implements extraction of constants from GSL. - -package Galacticus::Build::SourceTree::Process::GSLConstants; -use strict; -use warnings; -use utf8; -use Cwd; -use lib $ENV{'GALACTICUS_EXEC_PATH'}."/perl"; -use File::Temp qw/ tempfile /; - -# Insert hooks for our functions. -$Galacticus::Build::SourceTree::Hooks::processHooks{'gslConstant'} = \&Process_GSLConstant; - -sub Process_GSLConstant { - # Get the tree. - my $tree = shift(); - # Walk the tree. - my $node = $tree; - my $depth = 0; - while ( $node ) { - if ( $node->{'type'} eq "gslConstant" && ! $node->{'directive'}->{'processed'} ) { - # Generate, compile, and execute a minimal C code which simply outputs the relevant GSL constant. - my $type = exists($node->{'directive'}->{'type'}) ? $node->{'directive'}->{'type'} : "double precision"; - my $code; - if ( $type eq "enum" ) { - # An enum type. - my @members = split(/\s*,\s*/,$node->{'directive'}->{'members'}); - (undef, my $tmpSourceFileName) = tempfile('tempXXXXX', SUFFIX => '.c', OPEN => 0, DIR => $ENV{'BUILDPATH'}); - (undef, my $tmpExecFileName ) = tempfile('tempXXXXX', SUFFIX => '.x', OPEN => 0, DIR => $ENV{'BUILDPATH'}); - open(my $tmpSource,">".$tmpSourceFileName); - print $tmpSource "\#include \n"; - print $tmpSource "\#include \n"; - print $tmpSource "\#include {'directive'}->{'gslHeader'}.".h>\n"; - print $tmpSource "int main () {\n"; - print $tmpSource " printf(\"".join(" ",("%i") x scalar(@members))."\", ".$node->{'directive'}->{'members'}.");\n"; - print $tmpSource "}\n"; - close($tmpSource); - system($ENV{'CCOMPILER'}." -o ".$tmpExecFileName." ".$tmpSourceFileName." ".$ENV{'CFLAGS'}); - die("Galacticus::Build::SourceTree::Process::GSLConstants: failed to compile") - unless ( $? == 0 ); - open(my $gslPipe,$tmpExecFileName."|"); - my $gslEnum = <$gslPipe>; - close($gslPipe); - unlink($tmpSourceFileName,$tmpExecFileName); - my @values = split(" ",$gslEnum); - my @definitions; - for(my $i=0;$i '.c', OPEN => 0, DIR => $ENV{'BUILDPATH'}); - (undef, my $tmpExecFileName ) = tempfile('tempXXXXX', SUFFIX => '.x', OPEN => 0, DIR => $ENV{'BUILDPATH'}); - open(my $tmpSource,">".$tmpSourceFileName); - print $tmpSource "\#include \n"; - print $tmpSource "\#include \n"; - print $tmpSource "\#include {'directive'}->{'gslHeader'}.".h>\n"; - if ( $type eq "double precision" ) { - print $tmpSource "\#ifdef LDBL_DECIMAL_DIG\n"; - print $tmpSource " \#define OP_LDBL_Digs (LDBL_DECIMAL_DIG)\n"; - print $tmpSource "\#else \n"; - print $tmpSource " \#ifdef DECIMAL_DIG\n"; - print $tmpSource " \#define OP_LDBL_Digs (DECIMAL_DIG)\n"; - print $tmpSource " \#else \n"; - print $tmpSource " \#define OP_LDBL_Digs (LDBL_DIG + 3)\n"; - print $tmpSource " \#endif\n"; - print $tmpSource "\#endif\n"; - print $tmpSource "int main () {\n"; - print $tmpSource " printf(\"%.*e\", OP_LDBL_Digs - 1, ".$node->{'directive'}->{'gslSymbol'}.");\n"; - print $tmpSource "}\n"; - } elsif ( $type eq "integer" ) { - print $tmpSource "int main () {\n"; - print $tmpSource " printf(\"%i\", ".$node->{'directive'}->{'gslSymbol'}.");\n"; - print $tmpSource "}\n"; - } - close($tmpSource); - system($ENV{'CCOMPILER'}." -o ".$tmpExecFileName." ".$tmpSourceFileName." ".$ENV{'CFLAGS'}); - die("Galacticus::Build::SourceTree::Process::GSLConstants: failed to compile") - unless ( $? == 0 ); - open(my $gslPipe,$tmpExecFileName."|"); - my $gslConstant = <$gslPipe>; - close($gslPipe); - unlink($tmpSourceFileName,$tmpExecFileName); - $gslConstant =~ s/e/d/; - $code = $type.", parameter, public :: ".$node->{'directive'}->{'variable'}."=".$gslConstant."\n"; - } - my $newNode = - { - type => "code", - content => $code, - firstChild => undef() - }; - &Galacticus::Build::SourceTree::InsertAfterNode($node,[$newNode]); - # Mark the directive as processed. - $node->{'directive'}->{'processed'} = 1; - } - # Walk to the next node in the tree. - $node = &Galacticus::Build::SourceTree::Walk_Tree($node,\$depth); - } -} - -1; diff --git a/schema/constant.xsd b/schema/constant.xsd new file mode 100644 index 0000000000..4e3a9d559a --- /dev/null +++ b/schema/constant.xsd @@ -0,0 +1,22 @@ + + + + + + + + + + + + + + + + + + + + + + diff --git a/scripts/doc/Build_Documentation.sh b/scripts/doc/Build_Documentation.sh index e768ce1d57..79c9edaceb 100755 --- a/scripts/doc/Build_Documentation.sh +++ b/scripts/doc/Build_Documentation.sh @@ -7,7 +7,7 @@ PPN=1 FORCE=yes SUFFIX= -DIR=./work/build +DIR=./work/build/ CLEAN=no # Get options. @@ -62,6 +62,13 @@ if [ $? -ne 0 ]; then exit 1 fi +# Extract constants data. +./scripts/doc/constants.pl $DIR doc/constants.tex +if [ $? -ne 0 ]; then + echo Failed to extract constants data + exit 1 +fi + # Move to the documentation folder. cd doc diff --git a/scripts/doc/constants.pl b/scripts/doc/constants.pl new file mode 100755 index 0000000000..c2d0655221 --- /dev/null +++ b/scripts/doc/constants.pl @@ -0,0 +1,105 @@ +#!/usr/bin/env perl +use strict; +use warnings; +use lib $ENV{'GALACTICUS_EXEC_PATH'}."/perl"; +use XML::Simple; +use File::Find; +use List::ExtraUtils; +use Scalar::Util qw(looks_like_number); + +# Format constant definitions for incorporation into the documentation. +# Andrew Benson (29-January-2025) + +# Get arguments. +die("Usage: constants.pl ") + unless ( scalar(@ARGV) == 2 ); +my $buildPath = $ARGV[0]; +my $outputFile = $ARGV[1]; + +# Initialize a list to contain constants. +my @constants; + +# Process all files. +my @directories = ( $buildPath ); +find(\&process,@directories); + +# Define known groups. +my %knownGroups = + ( + "astrophysical" => "Astrophysical constants", + "atomic" => "Atomic physics constants", + "physical" => "Physical constants", + "math" => "Mathematical constants", + "units" => "Units", + "prefixes" => "SI Prefixes", + "GSL" => "Gnu Scientific Library constants", + "misc" => "Miscellaneous constants" + ); + +# Find all groups. +my $groups; +foreach my $constant ( @constants ) { + next + unless ( exists($constant->{'value'}) ); + if ( exists($constant->{'group'}) ) { + my @groupNames = split(":",$constant->{'group'}); + foreach my $groupName ( @groupNames ) { + die("Unknown group name '".$groupName."'") + unless ( exists($knownGroups{$groupName}) ); + push(@{$groups->{$groupName}->{'constant'}},$constant); + } + } else { + push(@{$groups->{'misc'}->{'constant'}},$constant); + } +} + +# Iterate over groups and constants, formatting into LaTeX. +open(my $output,">",$outputFile); +foreach my $groupName ( sort(keys(%{$groups})) ) { + print $output "\\subsection{".$knownGroups{$groupName}."}\n"; + my @constants = sort {$a->{'variable'} cmp $b->{'variable'}} @{$groups->{$groupName}->{'constant'}}; + foreach my $constant ( @constants ) { + next + unless ( exists($constant->{'value'}) ); + my $reference; + if ( exists($constant->{'referenceURL'}) ) { + $reference = "\\href{".$constant->{'referenceURL'}."}{".$constant->{'reference'}."}"; + } else { + $reference = $constant->{'reference'}; + } + $reference =~ s/&/\\&/; + (my $variable = $constant->{'variable'}) =~ s/_/\\_/g; + my $external = exists($constant->{'externalDescription'}) ? " (See \\href{".$constant->{'externalDescription'}."}{here}.)" : ""; + my $symbol = exists($constant->{'symbol'}) ? "\$".$constant->{'symbol'}."\$: " : ""; + my $units = exists($constant->{'units' }) && $constant->{'units'} ne "dimensionless" ? " [".$constant->{'units'}."]" : ""; + (my $value = $constant->{'value'}) =~ s/d([\+\-0-9]+)$/e$1/; + $value =~ s/_[_a-zA-Z0-9]+$//; + (my $module = $constant->{'module'}) =~ s/_/\\_/g; + (my $fileName = $constant->{'fileName'}) =~ s/\./_/g; + my $moduleURL = "https://github.com/galacticusorg/galacticus/releases/download/masterRelease/Galacticus_Source.pdf\\#source.".$fileName.":".lc($constant->{'module'}); + print $output "\\noindent {\\normalfont \\ttfamily ".$variable." = ".$value."}".$units."\\\\\n"; + print $output "\\indent ".$symbol.$constant->{'description'}.$external."\\\\\n"; + print $output "\\indent Module: \\href{".$moduleURL."}{\\normalfont \\ttfamily ".$module."}\\\\\n"; + print $output "\\indent Reference: ".$reference."\n\n\\medskip\n"; + } +} +close($output); + +exit; + +sub process { + # Process constant definition files, accumulating all constants. + my $fileName = $_; + my $fullName = $File::Find::name; + # Ignore files not containing constant definitions. + return + unless ( $fileName =~ m/\.constants.xml$/ ); + # Parse the file. + my $xml = new XML::Simple; + my $constantsFromFile = $xml->XMLin($fileName); + # Append to our list. + push(@constants,&List::ExtraUtils::as_array($constantsFromFile->{'constant'})); +} + + + diff --git a/source/atomic.radiation.gaunt_factors.vanHoof2014.F90 b/source/atomic.radiation.gaunt_factors.vanHoof2014.F90 index 164bcb5d0d..3ff3a2b075 100644 --- a/source/atomic.radiation.gaunt_factors.vanHoof2014.F90 +++ b/source/atomic.radiation.gaunt_factors.vanHoof2014.F90 @@ -136,7 +136,7 @@ double precision function vanHoof2014Total(self,atomicNumber,electronNumber,temp !!} use :: Error , only : Error_Report use :: Numerical_Constants_Physical, only : boltzmannsConstant - use :: Numerical_Constants_Units , only : rydbergs + use :: Numerical_Constants_Units , only : rydberg implicit none class (gauntFactorVanHoof2014), intent(inout) :: self integer , intent(in ) :: atomicNumber, electronNumber @@ -161,7 +161,7 @@ double precision function vanHoof2014Total(self,atomicNumber,electronNumber,temp & -electronNumber & & +1 & & )**2 & - & *rydbergs & + & *rydberg & & /boltzmannsConstant & & /temperature g =log10(gammaSquared) diff --git a/source/interface.GSL.F90 b/source/interface.GSL.F90 index 478fac7a2c..458f5691b1 100644 --- a/source/interface.GSL.F90 +++ b/source/interface.GSL.F90 @@ -145,27 +145,27 @@ Type for GSL special function results. ! Error codes. !![ - - - - - - - - - - - - - - + + + + + + + + + + + + + + !!] ! Precision modes. !![ - - - + + + !!] contains diff --git a/source/intergalactic_medium.filtering_mass.Gnedin2000.F90 b/source/intergalactic_medium.filtering_mass.Gnedin2000.F90 index aae84fee89..6efe0b04bc 100644 --- a/source/intergalactic_medium.filtering_mass.Gnedin2000.F90 +++ b/source/intergalactic_medium.filtering_mass.Gnedin2000.F90 @@ -541,7 +541,8 @@ Return the coefficients of the fitting function for the filtering mass at early class (intergalacticMediumFilteringMassGnedin2000), intent(inout) :: self double precision , dimension(4) :: coefficients double precision , intent(in ) :: time - double precision :: omegaMatter , expansionFactor, & + logical , save :: warnedOmega =.false., warnedRedshift =.false. + double precision :: omegaMatter , expansionFactor , & & redshift ! Extract matter density and redshift. @@ -549,8 +550,26 @@ Return the coefficients of the fitting function for the filtering mass at early expansionFactor=self%cosmologyFunctions_ %expansionFactor (time ) redshift =self%cosmologyFunctions_ %redshiftFromExpansionFactor(expansionFactor) ! Validate input. - if (omegaMatter < 0.25d0 .or. omegaMatter > 0.40d0) call Warn('gnedin2000CoefficientsEarlyEpoch: matter density outside validated range of fitting function; 0.25 ≤ Ωₘ ≤ 0.40') - if (redshift < 7.00d0 .or. redshift > 150.00d0) call Warn('gnedin2000CoefficientsEarlyEpoch: redshift outside validated range of fitting function; 7 ≤ z ≤ 150' ) + if (omegaMatter < 0.25d0 .or. omegaMatter > 0.40d0) then + if (.not.warnedOmega) then + !$omp critical (intergalacticMediumFilteringMassGnedin2000WarnedOmega) + if (.not.warnedOmega) then + call Warn('gnedin2000CoefficientsEarlyEpoch: matter density outside validated range of fitting function; 0.25 ≤ Ωₘ ≤ 0.40') + warnedOmega=.true. + end if + !$omp end critical (intergalacticMediumFilteringMassGnedin2000WarnedOmega) + end if + end if + if (redshift < 7.00d0 .or. redshift > 150.00d0) then + if (.not.warnedRedshift) then + !$omp critical (intergalacticMediumFilteringMassGnedin2000WarnedRedshift) + if (.not.warnedRedshift) then + call Warn('gnedin2000CoefficientsEarlyEpoch: redshift outside validated range of fitting function; 7 ≤ z ≤ 150' ) + warnedRedshift=.true. + end if + !$omp end critical (intergalacticMediumFilteringMassGnedin2000WarnedRedshift) + end if + end if ! Evaluate fitting function. coefficients(1)=-0.38d0*(omegaMatter**2)+ 0.41d0*omegaMatter- 0.16d0 coefficients(2)=+3.30d0*(omegaMatter**2)- 3.38d0*omegaMatter+ 1.15d0 diff --git a/source/math.linear_algebra.F90 b/source/math.linear_algebra.F90 index 261d2cc24b..3539870e0a 100644 --- a/source/math.linear_algebra.F90 +++ b/source/math.linear_algebra.F90 @@ -466,9 +466,9 @@ end function gsl_blas_dgemm ! Extract enums needed for BLAS functions. !![ - - - + + + !!] contains diff --git a/source/numerical.constants.astronomical.F90 b/source/numerical.constants.astronomical.F90 index 82000a77d5..ba77aee6ef 100644 --- a/source/numerical.constants.astronomical.F90 +++ b/source/numerical.constants.astronomical.F90 @@ -33,67 +33,70 @@ module Numerical_Constants_Astronomical implicit none public - ! Solar mass (in kg). + ! Physical properties of the Sun. !![ - + + + !!] - ! Solar radius (in m; Allen's Astrophysical Quantities, page 340). - double precision, parameter :: radiusSolar =6.95508d8 - - ! Solar luminosity (in W; Allen's Astrophysical Quantities, page 340). - double precision, parameter :: luminositySolar =3.845d26 - - ! Solar composition (Allen's Astrophysical Quantities, page 28). - double precision, parameter :: hydrogenByMassSolar =0.7070d0 - double precision, parameter :: heliumByMassSolar =0.2740d0 - double precision, parameter :: metallicitySolar =0.0188d0 - double precision, parameter :: heliumToHydrogenAbundanceSolar = (heliumByMassSolar /atomicMassHelium ) & - & /(hydrogenByMassSolar /atomicMassHydrogen) + ! Solar composition. + !![ + + + + + !!] - ! Primordial composition (Cyburt, Fields, & Olive; 2008; JCAP; 11; 12; https://ui.adsabs.harvard.edu/abs/2008JCAP...11..012C). - double precision, parameter :: hydrogenByMassPrimordial = 0.7514d+00 - double precision, parameter :: heliumByMassPrimordial = 0.2486d+00 ! Theoretical expectation based on WMAP results. - double precision, parameter :: lithiumToHydrogenAbundancePrimordial= 5.2400d-10 - double precision, parameter :: heliumToHydrogenAbundancePrimordial = (heliumByMassPrimordial /atomicMassHelium ) & - & /(hydrogenByMassPrimordial/atomicMassHydrogen) - double precision, parameter :: metallicityPrimordial = lithiumToHydrogenAbundancePrimordial*atomicMassLithium7/atomicMassHydrogen - double precision, parameter :: meanAtomicMassPrimordial = 1.0d0/(2.0d0*hydrogenByMassPrimordial/atomicMassHydrogen+3.0d0*heliumByMassPrimordial/atomicMassHelium) + ! Primordial composition. + !![ + + + + + + + !!] - ! Parsec and related quantities (in m). + ! Astrophysical scales. !![ - + + + + + !!] - double precision, parameter :: kiloParsec =kilo*parsec - double precision, parameter :: megaParsec =mega*parsec ! Newton's gravitational constant (in Galacticus' M_Solar, Mpc, km/s unit system). - double precision, parameter :: gravitationalConstantGalacticus=gravitationalConstant*massSolar/(kilo**2)/megaParsec - - ! Years and related quantities (in s). - double precision, parameter :: year =31558149.8d0 ! Sidereal year. - double precision, parameter :: gigaYear =giga*year + !![ + + !!] ! Conversion from Mpc/(km/s) to Gyr. - double precision, parameter :: Mpc_per_km_per_s_To_Gyr =megaParsec/kilo/gigaYear + !![ + + !!] ! Conversion from magnitudes to optical depth. - double precision, parameter :: magnitudesPerOpticalDepth = 2.5d0/log(10.0d0) + !![ + + !!] - ! AB magnitude system: - ! The AB magnitude system is defined such that: - ! m = -2.5log₁₀(F_nu/[ergs/s/cm²/Hz])-48.57 - ! Computing the flux at 10pc gives us the zero point for the absolute magnitude scale (units of W/Hz). - double precision, parameter :: offsetAB =48.57d0 - double precision, parameter :: luminosityZeroPointAB =(10.0d0**(-offsetAB/2.5d0))*4.0d0*Pi*((10.0d0*parsec*hecto)**2)*ergs + ! AB magnitude system. + !![ + + + !!] ! Angular conversions. - double precision, parameter :: arcminutesToDegrees = 1.0d0 / 60.0d0 - double precision, parameter :: arcsecondsToDegrees = 1.0d0 /3600.0d0 - double precision, parameter :: hoursToDegrees =360.0d0 / 24.0d0 - double precision, parameter :: minutesToDegrees =360.0d0 / 24.0d0/ 60.0d0 - double precision, parameter :: secondsToDegrees =360.0d0 / 24.0d0/3660.0d0 - double precision, parameter :: hoursToRadians = 2.0d0*Pi/ 24.0d0 - double precision, parameter :: degreesToRadians = 2.0d0*Pi/ 360.0d0 + !![ + + + + + + + + !!] end module Numerical_Constants_Astronomical diff --git a/source/numerical.constants.atomic.F90 b/source/numerical.constants.atomic.F90 index 75272871cb..51799e9d77 100644 --- a/source/numerical.constants.atomic.F90 +++ b/source/numerical.constants.atomic.F90 @@ -26,32 +26,23 @@ module Numerical_Constants_Atomic Contains various useful atomic constants. !!} use :: Numerical_Constants_Physical, only : electronMass , plancksConstant, speedLight - use :: Numerical_Constants_Units , only : angstromsPerMeter, rydbergs + use :: Numerical_Constants_Units , only : angstromsPerMeter, rydberg implicit none public - ! Atomic mass unit (in kg). + ! Atomic masses. !![ - + + + + + + !!] - ! Atomic masses. - double precision, parameter :: atomicMassHydrogen =1.007825d0 - double precision, parameter :: atomicMassHelium =4.002602d0 - double precision, parameter :: atomicMassLithium7 =7.016004d0 - - ! Mass of hydrogen and helium atom (in kg). - double precision, parameter :: massHydrogenAtom =atomicMassHydrogen*atomicMassUnit - double precision, parameter :: massHeliumAtom =atomicMassHelium *atomicMassUnit - ! Hydrogen Lyman series limit wavelength including correction for finite mass of the atom. - double precision, parameter :: lymanSeriesLimitWavelengthHydrogen=+plancksConstant & - & *speedLight & - & /rydbergs & - & *angstromsPerMeter & - & /( & - & +1.0d0 & - & -electronMass & - & /massHydrogenAtom & - & ) + !![ + + !!] + end module Numerical_Constants_Atomic diff --git a/source/numerical.constants.math.F90 b/source/numerical.constants.math.F90 index 0e40c82d83..6a9f704f84 100644 --- a/source/numerical.constants.math.F90 +++ b/source/numerical.constants.math.F90 @@ -25,41 +25,18 @@ module Numerical_Constants_Math !!{ Contains various useful mathematical constants. !!} - use :: Kind_Numbers , only : kind_quad + use :: Kind_Numbers, only : kind_quad implicit none - ! e. !![ - + + + + + + + + !!] - ! Pi. - !![ - - !!] - real(kind=kind_quad), public, parameter :: PiQuadPrecision=3.141592653589793238462643383279502884197_kind_quad - - ! ! Natural logarithm of 10. - !![ - - !!] - - ! Natural logarithm of 2. - !![ - - !!] - - ! Euler's constant. - !![ - - !!] - - ! Riemann zeta-function values. - !! ζ(3) - https://oeis.org/A002117 - double precision, public, parameter :: riemannZeta3=1.20205690315959428539973816151144999076d0 - - ! Catalan's constant - !! G - https://oeis.org/A006752 - double precision, public, parameter :: catalan =0.91596559417721901505460351493238411077d0 - end module Numerical_Constants_Math diff --git a/source/numerical.constants.physical.F90 b/source/numerical.constants.physical.F90 index 4980cc0adf..2cbef95b69 100644 --- a/source/numerical.constants.physical.F90 +++ b/source/numerical.constants.physical.F90 @@ -29,60 +29,19 @@ module Numerical_Constants_Physical implicit none public - ! Speed of light (m/s). !![ - + + + + + + + + + + + + !!] - ! Newton's gravitational constant (in SI units). - !![ - - !!] - - ! Stefan-Boltzmann constant (in units of J/s/m²/K⁴). - !![ - - !!] - - ! Radiation constant (in units of J/m³/K⁴). - double precision, parameter :: radiationConstant =4.0d0*stefanBoltzmannConstant/speedLight - - ! Boltzmann's constant (in units of J/K). - !![ - - !!] - - ! Thomson cross section (in units of m²). - !![ - - !!] - - ! Electron mass (in units of kg). - !![ - - !!] - - ! Planck's constant (in units of J s). - !![ - - !!] - - ! Electron Charge (in units of C). - !![ - - !!] - - ! Permitivity of free space (in SI units). - !![ - - !!] - - ! Fine structure constant (unitless) - !![ - - !!] - - ! Classical electron radius (m) - double precision, parameter :: electronRadius = 1.0d0 / (4.0d0 * Pi * eps0) * electronCharge**2 / (electronMass * speedLight**2) - end module Numerical_Constants_Physical diff --git a/source/numerical.constants.prefixes.F90 b/source/numerical.constants.prefixes.F90 index 70a11d1855..de8c488439 100644 --- a/source/numerical.constants.prefixes.F90 +++ b/source/numerical.constants.prefixes.F90 @@ -27,32 +27,35 @@ module Numerical_Constants_Prefixes !!} implicit none public - - double precision , parameter :: quetta =1.0d+30 - double precision , parameter :: ronna =1.0d+27 - double precision , parameter :: yotta =1.0d+24 - double precision , parameter :: zetta =1.0d+21 - double precision , parameter :: exa =1.0d+18 - double precision , parameter :: peta =1.0d+15 - double precision , parameter :: tera =1.0d+12 - double precision , parameter :: giga =1.0d+09 - double precision , parameter :: mega =1.0d+06 - double precision , parameter :: kilo =1.0d+03 - double precision , parameter :: hecto =1.0d+02 - double precision , parameter :: deca =1.0d+01 - double precision , parameter :: deci =1.0d-01 - double precision , parameter :: centi =1.0d-02 - double precision , parameter :: milli =1.0d-03 - double precision , parameter :: micro =1.0d-06 - double precision , parameter :: nano =1.0d-09 - double precision , parameter :: pico =1.0d-12 - double precision , parameter :: femto =1.0d-15 - double precision , parameter :: atto =1.0d-18 - double precision , parameter :: zepto =1.0d-21 - double precision , parameter :: yocto =1.0d-24 - double precision , parameter :: ronto =1.0d-27 - double precision , parameter :: quecto =1.0d-30 - character (len=2), parameter, dimension(-10:10) :: siPrefix=['q ','r ','y ','z ','a ','f ','p ','n ','μ','m ',' ','k ','M ','G ','T ','P ','E ','Z ','Y ','R ','Q '] + + !![ + + + + + + + + + + + + + + + + + + + + + + + + + !!] + + character(len=2), parameter, dimension(-10:10) :: siPrefix=['q ','r ','y ','z ','a ','f ','p ','n ','μ','m ',' ','k ','M ','G ','T ','P ','E ','Z ','Y ','R ','Q '] contains diff --git a/source/numerical.constants.units.F90 b/source/numerical.constants.units.F90 index c09649cb06..2c5b9692fa 100644 --- a/source/numerical.constants.units.F90 +++ b/source/numerical.constants.units.F90 @@ -29,32 +29,15 @@ module Numerical_Constants_Units implicit none public - ! Ergs in Joules. - double precision, parameter :: ergs =1.0d-7 - - ! Rydberg in Joules. !![ - + + + + + + + + !!] - ! Angstroms in microns. - double precision, parameter :: angstromsPerMicron=1.0d4 - - ! Angstroms in meters. - double precision, parameter :: angstromsPerMeter =1.0d10 - - ! Electron volt (in units of Joules). - !![ - - !!] - - ! Barn (cross section unit, in units of m²). - double precision, parameter :: barn =1.0d-28 - - ! Degree (in units of radians). - double precision, parameter :: degree =Pi/180.0d0 - - ! Arcsecond (in units of radians). - double precision, parameter :: arcsecond =Pi/180.0d0/60.0d0/60.0d0 - end module Numerical_Constants_Units diff --git a/source/numerical.integration.F90 b/source/numerical.integration.F90 index e78e952023..1af5e2cc78 100644 --- a/source/numerical.integration.F90 +++ b/source/numerical.integration.F90 @@ -36,12 +36,12 @@ module Numerical_Integration ! Integrator types. !![ - - - - - - + + + + + + !!] !![ diff --git a/testSuite/test-star-formation-histories.py b/testSuite/test-star-formation-histories.py index b8694e8db3..b624316a77 100755 --- a/testSuite/test-star-formation-histories.py +++ b/testSuite/test-star-formation-histories.py @@ -46,7 +46,7 @@ # Find Solar metallicity from the Galacticus source file. metallicitySolar = None for line in open(os.environ['GALACTICUS_EXEC_PATH']+'/source/numerical.constants.astronomical.F90'): - match = re.match(r'.*metallicitySolar\s*=\s*([0-9\.\+\-]+)d([0-9\+\-]+)',line) + match = re.match(r'.*variable="metallicitySolar"\s+value="([0-9\.\+\-]+)d([0-9\+\-]+)"',line) if match: metallicitySolar = float(match.group(1))*10**int(match.group(2)) if metallicitySolar is None: