[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master c40e069: Cosmology functions moved into the li
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master c40e069: Cosmology functions moved into the library |
|
Date: |
Tue, 3 Oct 2017 10:43:01 -0400 (EDT) |
branch: master
commit c40e0694b638b86dcc57d025ed70d7e92588174f
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Cosmology functions moved into the library
Until now, Gnuastro's cosmology functions were only available with the
CosmicCalculator program. But these functions can be useful in many other
contexts (other future Gnuastro programs or for users of the library).
With this commit, these functions are now defined in the library and
available for easy usage.
This was done with the help of Boud Roukema.
This completes task #14665.
---
NEWS | 4 +
bin/cosmiccal/cosmiccal.c | 240 ++++++-----------------------------
bin/cosmiccal/cosmiccal.h | 5 -
bin/cosmiccal/main.h | 5 -
bin/cosmiccal/ui.c | 11 +-
doc/gnuastro.texi | 60 ++++++++-
lib/Makefile.am | 25 ++--
lib/cosmology.c | 312 ++++++++++++++++++++++++++++++++++++++++++++++
lib/gnuastro/cosmology.h | 96 ++++++++++++++
9 files changed, 521 insertions(+), 237 deletions(-)
diff --git a/NEWS b/NEWS
index 67992d1..6ea4787 100644
--- a/NEWS
+++ b/NEWS
@@ -36,6 +36,10 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
now possible to detect signal out to much lower surface brightness limits
and the detections don't look boxy any more.
+ Cosmology library: A new set of cosmology functions are now included in
+ the library (declared in `gnuastro/cosmology.h'). These functions are
+ also used in the CosmicCalculator program.
+
`gal_fits_key_img_blank': returns the value that must be used in the
BLANK keyword for the given type as defined by the FITS standard.
diff --git a/bin/cosmiccal/cosmiccal.c b/bin/cosmiccal/cosmiccal.c
index d040154..e1755fa 100644
--- a/bin/cosmiccal/cosmiccal.c
+++ b/bin/cosmiccal/cosmiccal.c
@@ -26,202 +26,13 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdio.h>
#include <stdlib.h>
-#include <gsl/gsl_const_mksa.h>
-#include <gsl/gsl_integration.h>
+#include <gnuastro/cosmology.h>
#include "main.h"
#include "cosmiccal.h"
-/**************************************************************/
-/************ Integrand functions *************/
-/**************************************************************/
-/* In these functions, z as a separate argument, this is not necessarily
- the same z as the redshift in cosmiccalparams. */
-
-double
-Ez(double z, void *params)
-{
- struct cosmiccalparams *p=(struct cosmiccalparams *)params;
- return sqrt( p->olambda
- + p->ocurv * (1+z) * (1+z)
- + p->omatter * (1+z) * (1+z) * (1+z)
- + p->oradiation * (1+z) * (1+z) * (1+z) * (1+z));
-}
-
-
-
-
-
-double
-age(double z, void *params)
-{
- return 1 / ( (1+z)*Ez(z,params) );
-}
-
-
-
-
-
-double
-propdist(double z, void *params)
-{
- return 1 / ( Ez(z,params) );
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/**************************************************************/
-/************ Integrators *************/
-/**************************************************************/
-/* Estimate the age of the universe, note that z might be different
- from the desired redshift. */
-double
-ageofuniverse(struct cosmiccalparams *p, double z)
-{
- gsl_function F;
- double result, error;
- gsl_integration_workspace *w=gsl_integration_workspace_alloc(GSLILIMIT);
-
- /* Set the GSL function parameters */
- F.params=p;
- F.function=&age;
-
- gsl_integration_qagiu(&F, z, GSLIEPSABS, GSLIEPSREL, GSLILIMIT, w,
- &result, &error);
-
- return result / p->H0s / (365*GSL_CONST_MKSA_DAY) / 1e9;
-}
-
-
-
-
-
-/* Return the proper distance to a source at z in units of mega
- parsecs */
-double
-properdistance(struct cosmiccalparams *p, double z)
-{
- size_t neval;
- gsl_function F;
- double result, error;
-
- /* Set the GSL function parameters */
- F.params=p;
- F.function=&propdist;
-
- gsl_integration_qng(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL,
- &result, &error, &neval);
-
- return result * p->c / p->H0s / (1e6 * GSL_CONST_MKSA_PARSEC);
-}
-
-
-
-
-
-double
-covolume(double z, void *params)
-{
- size_t neval;
- gsl_function F;
- double result, error;
-
- /* Set the GSL function parameters */
- F.params=params;
- F.function=&propdist;
-
- gsl_integration_qng(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL,
- &result, &error, &neval);
-
- return result * result / ( Ez(z,params) );
-}
-
-
-
-
-
-double
-comovingvolume(struct cosmiccalparams *p, double z)
-{
- size_t neval;
- gsl_function F;
- double result, error;
- double cH = p->c / p->H0s / (1e6 * GSL_CONST_MKSA_PARSEC);
-
- /* Set the GSL function parameters */
- F.params=p;
- F.function=&covolume;
-
- gsl_integration_qng(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL,
- &result, &error, &neval);
-
- return result * 4 * M_PI * cH*cH*cH;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/**************************************************************/
-/************ Intermediary functions *************/
-/**************************************************************/
-/* Critical density at redshift z in units of gram/cm^3*/
-double
-criticaldensity(struct cosmiccalparams *p, double z)
-{
- double H = p->H0s*Ez(z,p);
- return 3*H*H/(8*M_PI*GSL_CONST_MKSA_GRAVITATIONAL_CONSTANT)/1000;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
@@ -243,11 +54,14 @@ cosmiccal(struct cosmiccalparams *p)
/* In case the user just wants one number, only print that and
return. */
if(p->onlyvolume){
- printf("%f\n", comovingvolume(p,p->redshift));
+ printf("%f\n", gal_cosmology_comoving_volume(p->redshift, p->H0,
+ p->olambda, p->omatter,
+ p->oradiation));
return;
}
if(p->onlyabsmagconv){
- pd=properdistance(p, p->redshift);
+ pd=gal_cosmology_proper_distance(p->redshift, p->H0, p->olambda,
+ p->omatter, p->oradiation);
ld=pd*(1+p->redshift);
distmod=5*(log10(ld*1000000)-1);
absmagconv=distmod-2.5*log10(1+p->redshift);
@@ -257,17 +71,35 @@ cosmiccal(struct cosmiccalparams *p)
/* The user wants everything, do all the calculations and print
everything with full descriptions. */
- curage=ageofuniverse(p, 0.0f);
- ccritd=criticaldensity(p, 0.0f);
- vz=comovingvolume(p,p->redshift);
- pd=properdistance(p, p->redshift);
- outage=ageofuniverse(p, p->redshift);
- zcritd=criticaldensity(p, p->redshift);
-
- ad=pd/(1+p->redshift);
- ld=pd*(1+p->redshift);
- distmod=5*(log10(ld*1000000)-1);
- absmagconv=distmod-2.5*log10(1+p->redshift);
+ curage=gal_cosmology_age(0.0f, p->H0, p->olambda, p->omatter,
+ p->oradiation);
+
+ ccritd=gal_cosmology_critical_density(0.0f, p->H0, p->olambda, p->omatter,
+ p->oradiation);
+
+ vz=gal_cosmology_comoving_volume(p->redshift, p->H0, p->olambda, p->omatter,
+ p->oradiation);
+
+ pd=gal_cosmology_proper_distance(p->redshift, p->H0, p->olambda, p->omatter,
+ p->oradiation);
+
+ outage=gal_cosmology_age(p->redshift, p->H0, p->olambda, p->omatter,
+ p->oradiation);
+
+ zcritd=gal_cosmology_critical_density(p->redshift, p->H0, p->olambda,
+ p->omatter, p->oradiation);
+
+ ad=gal_cosmology_angular_distance(p->redshift, p->H0, p->olambda, p->omatter,
+ p->oradiation);
+
+ ld=gal_cosmology_luminosity_distance(p->redshift, p->H0, p->olambda,
+ p->omatter, p->oradiation);
+
+ distmod=gal_cosmology_distance_modulus(p->redshift, p->H0, p->olambda,
+ p->omatter, p->oradiation);
+
+ absmagconv=gal_cosmology_to_absolute_mag(p->redshift, p->H0, p->olambda,
+ p->omatter, p->oradiation);
/* Print out results: */
printf("%s\n", PROGRAM_STRING);
@@ -280,7 +112,7 @@ cosmiccal(struct cosmiccalparams *p)
printf(FLTFORMAT, "Matter fractional density, now:", p->omatter);
printf(EXPFORMAT, "Radiation fractional density, now:", p->oradiation);
printf(EXPFORMAT, "Curvatue fractional density (from the above):",
- p->ocurv);
+ 1 - ( p->olambda + p->omatter + p->oradiation ));
printf("\n\n Universe now\n");
diff --git a/bin/cosmiccal/cosmiccal.h b/bin/cosmiccal/cosmiccal.h
index 99e2800..946a3a7 100644
--- a/bin/cosmiccal/cosmiccal.h
+++ b/bin/cosmiccal/cosmiccal.h
@@ -23,11 +23,6 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#ifndef COSMICCAL_H
#define COSMICCAL_H
-/* For the GSL integrations */
-#define GSLILIMIT 1000
-#define GSLIEPSABS 0
-#define GSLIEPSREL 1e-7
-
/* Format of final outputs */
# define FLTFORMAT " - %-55s%f\n"
# define EXPFORMAT " - %-55s%e\n"
diff --git a/bin/cosmiccal/main.h b/bin/cosmiccal/main.h
index 248b348..b8a04c7 100644
--- a/bin/cosmiccal/main.h
+++ b/bin/cosmiccal/main.h
@@ -58,11 +58,6 @@ struct cosmiccalparams
uint8_t onlyabsmagconv; /* Only print conversion to abs. mag. */
/* Internal: */
- double K; /* Curvature constant. */
- double c; /* Speed of light. */
- double H0s; /* Current expansion rate (1/sec). */
- double ocurv; /* Curvature density today. */
-
time_t rawtime; /* Starting time of the program. */
};
diff --git a/bin/cosmiccal/ui.c b/bin/cosmiccal/ui.c
index a7d8d98..0e697a6 100644
--- a/bin/cosmiccal/ui.c
+++ b/bin/cosmiccal/ui.c
@@ -107,9 +107,6 @@ ui_initialize_options(struct cosmiccalparams *p,
cp->program_authors = PROGRAM_AUTHORS;
cp->coptions = gal_commonopts_options;
- /* Speed of light: */
- p->c = GSL_CONST_MKSA_SPEED_OF_LIGHT;
-
/* Modify the common options. */
for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
{
@@ -195,7 +192,7 @@ parse_opt(int key, char *arg, struct argp_state *state)
static void
ui_read_check_only_options(struct cosmiccalparams *p)
{
- double sum=p->olambda + p->omatter + p->oradiation;
+ double sum = p->olambda + p->omatter + p->oradiation;
/* Check if the density fractions add up to 1 (within floating point
error). */
@@ -204,12 +201,6 @@ ui_read_check_only_options(struct cosmiccalparams *p)
"The cosmological constant (`olambda'), matter (`omatter') "
"and radiation (`oradiation') densities are given as %.8f, %.8f, "
"%.8f", sum, p->olambda, p->omatter, p->oradiation);
-
- /* The curvature fractional density: */
- p->ocurv=1-sum;
-
- /* Convert H0 from km/sec/Mpc to 1/sec: */
- p->H0s=p->H0/1000/GSL_CONST_MKSA_PARSEC;
}
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index babd660..305dcfc 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -555,6 +555,7 @@ Gnuastro library
* Convolution functions:: Library functions to do convolution.
* Interpolation:: Interpolate (over blank values possibly).
* Git wrappers:: Wrappers for functions in libgit2.
+* Cosmology library:: Cosmological calculations.
Multithreaded programming (@file{threads.h})
@@ -17252,6 +17253,7 @@ documentation will correspond to your installed version.
* Convolution functions:: Library functions to do convolution.
* Interpolation:: Interpolate (over blank values possibly).
* Git wrappers:: Wrappers for functions in libgit2.
+* Cosmology library:: Cosmological calculations.
@end menu
@node Configuration information, Multithreaded programming, Gnuastro library,
Gnuastro library
@@ -22084,7 +22086,7 @@ checking that is the most CPU intensive part will only
be done once.
@end deftypefun
address@hidden Git wrappers, , Interpolation, Gnuastro library
address@hidden Git wrappers, Cosmology library, Interpolation, Gnuastro library
@subsection Git wrappers (@file{git.h})
@cindex Git
@@ -22121,6 +22123,62 @@ controlled directory, then the output will be the
@code{NULL} pointer.
@end deftypefun
address@hidden Cosmology library, , Git wrappers, Gnuastro library
address@hidden Cosmology library (@file{cosmology.h})
+
+This library does the main cosmological calculations that are commonly
+necessary in extra-galactic astronomical studies. The main variable in this
+context is the redshift (@mymath{z}). The cosmological input parameters in
+the functions below are @code{H0}, @code{o_lambda_0}, @code{o_matter_0},
address@hidden which respectively represent the current (at redshift
+0) expansion rate (Hubble constant in units of km/sec/Mpc), cosmological
+constant (@mymath{\Lambda}), matter and radiation densities.
+
+All these functions are declared in @file{gnuastro/cosmology.h}. For a more
+extended introduction/discussion of the cosmological parameters, please see
address@hidden
+
address@hidden double gal_cosmology_age (double @code{z}, double @code{H0},
double @code{o_lambda_0}, double @code{o_matter_0}, double @code{o_radiation_0})
+Returns the age of the universe at redshift @code{z} in units of Giga
+years.
address@hidden deftypefun
+
address@hidden double gal_cosmology_proper_distance (double @code{z}, double
@code{H0}, double @code{o_lambda_0}, double @code{o_matter_0}, double
@code{o_radiation_0})
+Returns the proper distance to an object at redshift @code{z} in units of
+Mega parsecs.
address@hidden deftypefun
+
address@hidden double gal_cosmology_comoving_volume (double @code{z}, double
@code{H0}, double @code{o_lambda_0}, double @code{o_matter_0}, double
@code{o_radiation_0})
+Returns the comoving volume over 4pi stradian to @code{z} in units of Mega
+parsecs cube.
address@hidden deftypefun
+
address@hidden double gal_cosmology_critical_density (double @code{z}, double
@code{H0}, double @code{o_lambda_0}, double @code{o_matter_0}, double
@code{o_radiation_0})
+Returns the critical density at redshift @code{z} in units of
address@hidden/cm^3}.
address@hidden deftypefun
+
address@hidden double gal_cosmology_angular_distance (double @code{z}, double
@code{H0}, double @code{o_lambda_0}, double @code{o_matter_0}, double
@code{o_radiation_0})
+Return the angular diameter distance to an object at redshift @code{z} in
+units of Mega parsecs.
address@hidden deftypefun
+
address@hidden double gal_cosmology_luminosity_distance (double @code{z},
double @code{H0}, double @code{o_lambda_0}, double @code{o_matter_0}, double
@code{o_radiation_0})
+Return the luminosity diameter distance to an object at redshift @code{z}
+in units of Mega parsecs.
address@hidden deftypefun
+
address@hidden double gal_cosmology_distance_modulus (double @code{z}, double
@code{H0}, double @code{o_lambda_0}, double @code{o_matter_0}, double
@code{o_radiation_0})
+Return the distance modulus at redshift @code{z} (with no units).
address@hidden deftypefun
+
address@hidden double gal_cosmology_to_absolute_mag (double @code{z}, double
@code{H0}, double @code{o_lambda_0}, double @code{o_matter_0}, double
@code{o_radiation_0})
+Return the conversion from apparent to absolute magnitdue for an object at
+redshift @code{z}. This value has to be added to the apparent magnitude to
+give the absolute magnitude of an object at redshift @code{z}.
address@hidden deftypefun
+
+
@node Library demo programs, , Gnuastro library, Library
@section Library demo programs
diff --git a/lib/Makefile.am b/lib/Makefile.am
index db56754..9e2bf3c 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -50,11 +50,11 @@ libgnuastro_la_LIBADD =
$(top_builddir)/bootstrapped/lib/libgnu.la
# Specify the library .c files
-libgnuastro_la_SOURCES = arithmetic.c arithmetic-binary.c \
- arithmetic-onlyint.c binary.c blank.c box.c checkset.c convolve.c data.c \
- fits.c git.c interpolate.c list.c options.c permutation.c polygon.c \
- qsort.c dimension.c statistics.c table.c tableintern.c threads.c tile.c \
- timing.c txt.c type.c wcs.c
+libgnuastro_la_SOURCES = arithmetic.c arithmetic-binary.c \
+ arithmetic-onlyint.c binary.c blank.c box.c checkset.c convolve.c \
+ cosmology.c data.c fits.c git.c interpolate.c list.c options.c \
+ permutation.c polygon.c qsort.c dimension.c statistics.c table.c \
+ tableintern.c threads.c tile.c timing.c txt.c type.c wcs.c
@@ -66,13 +66,14 @@ libgnuastro_la_SOURCES = arithmetic.c arithmetic-binary.c
\
# in the $(headersdir) directory. Some of the header files don't need to be
# installed.
headersdir=$(top_srcdir)/lib/gnuastro
-pkginclude_HEADERS = gnuastro/config.h $(headersdir)/arithmetic.h \
- $(headersdir)/binary.h $(headersdir)/blank.h $(headersdir)/box.h \
- $(headersdir)/convolve.h $(headersdir)/data.h $(headersdir)/dimension.h \
- $(headersdir)/fits.h $(headersdir)/git.h $(headersdir)/interpolate.h \
- $(headersdir)/list.h $(headersdir)/permutation.h $(headersdir)/polygon.h \
- $(headersdir)/qsort.h $(headersdir)/statistics.h $(headersdir)/table.h \
- $(headersdir)/threads.h $(headersdir)/tile.h $(headersdir)/txt.h \
+pkginclude_HEADERS = gnuastro/config.h $(headersdir)/arithmetic.h \
+ $(headersdir)/binary.h $(headersdir)/blank.h $(headersdir)/box.h \
+ $(headersdir)/convolve.h $(headersdir)/cosmology.h $(headersdir)/data.h \
+ $(headersdir)/dimension.h $(headersdir)/fits.h $(headersdir)/git.h \
+ $(headersdir)/interpolate.h $(headersdir)/list.h \
+ $(headersdir)/permutation.h $(headersdir)/polygon.h \
+ $(headersdir)/qsort.h $(headersdir)/statistics.h $(headersdir)/table.h \
+ $(headersdir)/threads.h $(headersdir)/tile.h $(headersdir)/txt.h \
$(headersdir)/type.h $(headersdir)/wcs.h
diff --git a/lib/cosmology.c b/lib/cosmology.c
new file mode 100644
index 0000000..15532cf
--- /dev/null
+++ b/lib/cosmology.c
@@ -0,0 +1,312 @@
+/*********************************************************************
+Cosmological calculations.
+This is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+ Mohammad Akhlaghi <address@hidden>
+Contributing author(s):
+Copyright (C) 2015, Free Software Foundation, Inc.
+
+Gnuastro is free software: you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation, either version 3 of the License, or (at your
+option) any later version.
+
+Gnuastro is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
+**********************************************************************/
+#include <config.h>
+
+#include <time.h>
+#include <errno.h>
+#include <error.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <gsl/gsl_const_mksa.h>
+#include <gsl/gsl_integration.h>
+
+
+
+
+/**************************************************************/
+/************ Definitions *************/
+/**************************************************************/
+/* These are basic definitions that commonly go into the header files. But
+ because this is a library and the user imports the header file, it is
+ easier to just have them here in the main C file to avoid filling up the
+ user's name-space with junk. */
+struct cosmology_integrand_t
+{
+ double o_lambda_0;
+ double o_curv_0;
+ double o_matter_0;
+ double o_radiation_0;
+};
+
+
+
+
+
+/* For the GSL integrations */
+#define GSLILIMIT 1000
+#define GSLIEPSABS 0
+#define GSLIEPSREL 1e-7
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/**************************************************************/
+/************ Integrand functions *************/
+/**************************************************************/
+/* These are integrands, they won't be giving the final value. */
+static double
+cosmology_integrand_Ez(double z, void *params)
+{
+ struct cosmology_integrand_t *p=(struct cosmology_integrand_t *)params;
+ return sqrt( p->o_lambda_0
+ + p->o_curv_0 * (1+z) * (1+z)
+ + p->o_matter_0 * (1+z) * (1+z) * (1+z)
+ + p->o_radiation_0 * (1+z) * (1+z) * (1+z) * (1+z));
+}
+
+
+
+
+
+static double
+cosmology_integrand_age(double z, void *params)
+{
+ return 1 / ( (1.0 + z) * cosmology_integrand_Ez(z,params) );
+}
+
+
+
+
+
+static double
+cosmology_integrand_proper_dist(double z, void *params)
+{
+ return 1 / ( cosmology_integrand_Ez(z,params) );
+}
+
+
+
+
+
+static double
+cosmology_integrand_comoving_volume(double z, void *params)
+{
+ size_t neval;
+ gsl_function F;
+ double result, error;
+
+ /* Set the GSL function parameters */
+ F.params=params;
+ F.function=&cosmology_integrand_proper_dist;
+
+ gsl_integration_qng(&F, 0.0, z, GSLIEPSABS, GSLIEPSREL,
+ &result, &error, &neval);
+
+ return result * result / ( cosmology_integrand_Ez(z,params) );
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/**************************************************************/
+/************ Final functions *************/
+/**************************************************************/
+/* Age of the universe (in Gyrs). H0 is in units of (km/sec/Mpc) and the
+ fractional densities must add up to 1. */
+double
+gal_cosmology_age(double z, double H0, double o_lambda_0, double o_matter_0,
+ double o_radiation_0)
+{
+ gsl_function F;
+ double result, error;
+ double o_curv_0 = 1.0 - ( o_lambda_0 + o_matter_0 + o_radiation_0 );
+ double H0s=H0/1000/GSL_CONST_MKSA_PARSEC; /* H0 in units of seconds. */
+ gsl_integration_workspace *w=gsl_integration_workspace_alloc(GSLILIMIT);
+ struct cosmology_integrand_t p={o_lambda_0, o_curv_0, o_matter_0,
+ o_radiation_0};
+
+ /* Set the GSL function parameters. */
+ F.params=&p;
+ F.function=&cosmology_integrand_age;
+ gsl_integration_qagiu(&F, z, GSLIEPSABS, GSLIEPSREL, GSLILIMIT, w,
+ &result, &error);
+
+ return result / H0s / (365*GSL_CONST_MKSA_DAY) / 1e9;
+}
+
+
+
+
+
+/* Proper distance to z (Mpc). */
+double
+gal_cosmology_proper_distance(double z, double H0, double o_lambda_0,
+ double o_matter_0, double o_radiation_0)
+{
+ size_t neval;
+ gsl_function F;
+ double result, error, c=GSL_CONST_MKSA_SPEED_OF_LIGHT;
+ double o_curv_0 = 1.0 - ( o_lambda_0 + o_matter_0 + o_radiation_0 );
+ double H0s=H0/1000/GSL_CONST_MKSA_PARSEC; /* H0 in units of seconds. */
+ struct cosmology_integrand_t p={o_lambda_0, o_curv_0, o_matter_0,
+ o_radiation_0};
+
+ /* Set the GSL function parameters */
+ F.params=&p;
+ F.function=&cosmology_integrand_proper_dist;
+
+ /* Do the integration. */
+ gsl_integration_qng(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL, &result,
+ &error, &neval);
+
+ /* Return the result. */
+ return result * c / H0s / (1e6 * GSL_CONST_MKSA_PARSEC);
+}
+
+
+
+
+
+/* Comoving volume over 4pi stradian to z (Mpc^3). */
+double
+gal_cosmology_comoving_volume(double z, double H0, double o_lambda_0,
+ double o_matter_0, double o_radiation_0)
+{
+ size_t neval;
+ gsl_function F;
+ double result, error;
+ double c=GSL_CONST_MKSA_SPEED_OF_LIGHT;
+ double H0s=H0/1000/GSL_CONST_MKSA_PARSEC; /* H0 in units of seconds. */
+ double cH = c / H0s / (1e6 * GSL_CONST_MKSA_PARSEC);
+ double o_curv_0 = 1.0 - ( o_lambda_0 + o_matter_0 + o_radiation_0 );
+ struct cosmology_integrand_t p={o_lambda_0, o_curv_0, o_matter_0,
+ o_radiation_0};
+
+ /* Set the GSL function parameters */
+ F.params=&p;
+ F.function=&cosmology_integrand_comoving_volume;
+
+ /* Do the integration. */
+ gsl_integration_qng(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL,
+ &result, &error, &neval);
+
+ /* Return the result. */
+ return result * 4 * M_PI * cH*cH*cH;
+}
+
+
+
+
+
+/* Critical density at redshift z in units of g/cm^3. */
+double
+gal_cosmology_critical_density(double z, double H0, double o_lambda_0,
+ double o_matter_0, double o_radiation_0)
+{
+ double H;
+ double H0s=H0/1000/GSL_CONST_MKSA_PARSEC; /* H0 in units of seconds. */
+ double o_curv_0 = 1.0 - ( o_lambda_0 + o_matter_0 + o_radiation_0 );
+ struct cosmology_integrand_t p={o_lambda_0, o_curv_0, o_matter_0,
+ o_radiation_0};
+
+ /* Set the place holder, then return the result. */
+ H = H0s * cosmology_integrand_Ez(z, &p);
+ return 3*H*H/(8*M_PI*GSL_CONST_MKSA_GRAVITATIONAL_CONSTANT)/1000;
+}
+
+
+
+
+
+/* Angular diameter distance to z (Mpc). */
+double
+gal_cosmology_angular_distance(double z, double H0, double o_lambda_0,
+ double o_matter_0, double o_radiation_0)
+{
+ return gal_cosmology_proper_distance(z, H0, o_lambda_0, o_matter_0,
+ o_radiation_0) / (1+z);
+}
+
+
+
+
+
+/* Luminosity distance to z (Mpc). */
+double
+gal_cosmology_luminosity_distance(double z, double H0, double o_lambda_0,
+ double o_matter_0, double o_radiation_0)
+{
+ return gal_cosmology_proper_distance(z, H0, o_lambda_0, o_matter_0,
+ o_radiation_0) * (1+z);
+}
+
+
+
+
+
+/* Distance modulus at z (no units). */
+double
+gal_cosmology_distance_modulus(double z, double H0, double o_lambda_0,
+ double o_matter_0, double o_radiation_0)
+{
+ double ld=gal_cosmology_luminosity_distance(z, H0, o_lambda_0, o_matter_0,
+ o_radiation_0);
+ return 5*(log10(ld*1000000)-1);
+}
+
+
+
+
+
+/* Convert apparent to absolute magnitude. */
+double
+gal_cosmology_to_absolute_mag(double z, double H0, double o_lambda_0,
+ double o_matter_0, double o_radiation_0)
+{
+ double dm=gal_cosmology_distance_modulus(z, H0, o_lambda_0, o_matter_0,
+ o_radiation_0);
+ return dm-2.5*log10(1.0+z);
+}
diff --git a/lib/gnuastro/cosmology.h b/lib/gnuastro/cosmology.h
new file mode 100644
index 0000000..44d5d52
--- /dev/null
+++ b/lib/gnuastro/cosmology.h
@@ -0,0 +1,96 @@
+/*********************************************************************
+Cosmological calculations.
+This is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+ Mohammad Akhlaghi <address@hidden>
+Contributing author(s):
+Copyright (C) 2015, Free Software Foundation, Inc.
+
+Gnuastro is free software: you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation, either version 3 of the License, or (at your
+option) any later version.
+
+Gnuastro is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
+**********************************************************************/
+#ifndef __GAL_COSMOLOGY_H__
+#define __GAL_COSMOLOGY_H__
+
+/* Include other headers if necessary here. Note that other header files
+ must be included before the C++ preparations below */
+
+
+
+/* C++ Preparations */
+#undef __BEGIN_C_DECLS
+#undef __END_C_DECLS
+#ifdef __cplusplus
+# define __BEGIN_C_DECLS extern "C" {
+# define __END_C_DECLS }
+#else
+# define __BEGIN_C_DECLS /* empty */
+# define __END_C_DECLS /* empty */
+#endif
+/* End of C++ preparations */
+
+
+
+/* Actual header contants (the above were for the Pre-processor). */
+__BEGIN_C_DECLS /* From C++ preparations */
+
+
+
+
+/* Age of the universe (in Gyrs). */
+double
+gal_cosmology_age(double z, double H0, double o_lambda_0, double o_matter_0,
+ double o_radiation_0);
+
+/* Proper distance to z (Mpc). */
+double
+gal_cosmology_proper_distance(double z, double H0, double o_lambda_0,
+ double o_matter_0, double o_radiation_0);
+
+/* Comoving volume over 4pi stradian to z (Mpc^3). */
+double
+gal_cosmology_comoving_volume(double z, double H0, double o_lambda_0,
+ double o_matter_0, double o_radiation_0);
+
+/* Critical density at redshift z in units of g/cm^3. */
+double
+gal_cosmology_critical_density(double z, double H0, double o_lambda_0,
+ double o_matter_0, double o_radiation_0);
+
+/* Angular diameter distance to z (Mpc). */
+double
+gal_cosmology_angular_distance(double z, double H0, double o_lambda_0,
+ double o_matter_0, double o_radiation_0);
+
+/* Luminosity distance to z (Mpc). */
+double
+gal_cosmology_luminosity_distance(double z, double H0, double o_lambda_0,
+ double o_matter_0, double o_radiation_0);
+
+/* Distance modulus at z (no units). */
+double
+gal_cosmology_distance_modulus(double z, double H0, double o_lambda_0,
+ double o_matter_0, double o_radiation_0);
+
+/* Convert apparent to absolute magnitude. */
+double
+gal_cosmology_to_absolute_mag(double z, double H0, double o_lambda_0,
+ double o_matter_0, double o_radiation_0);
+
+
+
+
+__END_C_DECLS /* From C++ preparations */
+
+#endif /* __GAL_COSMOLOGY_H__ */
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master c40e069: Cosmology functions moved into the library,
Mohammad Akhlaghi <=