gnuastro-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[gnuastro-commits] master 5a6cf1d1: CosmicCalculator: allow curved FLRW


From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master 5a6cf1d1: CosmicCalculator: allow curved FLRW models for radial distances
Date: Fri, 2 Feb 2024 12:21:37 -0500 (EST)

branch: master
commit 5a6cf1d158c112613cd56a35e4e1b3aee45f36b0
Author: Boud Roukema <boud@cosmo.torun.pl>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    CosmicCalculator: allow curved FLRW models for radial distances
    
    Until now, the two files 'bin/cosmiccal/ui.c' and 'lib/cosmology.c' would
    prevent the calculation of cosmological distances for non-flat FLRW
    models. If the aim of astcosmiccal is to be limited strictly to LambdaCDM,
    then there should be much stricter limits on all the FLRW local
    cosmological parameters, not just their sum.
    
    With this commit, the error message in 'bin/cosmiccal/ui.c' is modified to
    instead report the current value of the curvature density parameter, and
    the message in 'lib/cosmology.c' is modified to remind the user that the
    model is non-flat. In both cases, the error message does not abort the
    program. To let the user avoid the library warning message, a new 'quiet'
    argument has been added to all the affected 'gal_cosmology_*' functions.
    
    Note that this change is *only* valid for the radial comoving distance,
    i.e. the '--properdistance' at the current epoch. The angular diameter
    distance is wrong for the non-flat cases, because it is not coded in
    astcosmiccalc. Currently, cosmdist [1] can be used for more generic FLRW
    distance/age/redshift calculations.
    
    [1] https://codeberg.org/boud/cosmdist
---
 NEWS                      | 20 +++++++++++-
 bin/cosmiccal/cosmiccal.c | 65 ++++++++++++++++++++++++---------------
 bin/cosmiccal/ui.c        | 13 ++++----
 lib/cosmology.c           | 78 ++++++++++++++++++++++++++++++-----------------
 lib/gnuastro/cosmology.h  | 23 +++++++++-----
 5 files changed, 131 insertions(+), 68 deletions(-)

diff --git a/NEWS b/NEWS
index 1be7789d..f230dcf1 100644
--- a/NEWS
+++ b/NEWS
@@ -213,10 +213,17 @@ See the end of the file for license conditions.
     on Gnuastro's Matrix chat channel (#gnuastro:openastronomy.org).
 
 *** CosmicCalculator
+  - In a non-FLRW model (when the sum of the densities is not 1.0),
+    CosmicCalculator no longer crashes with an error since calculations
+    like the radial comoving distance ('--properdistance') are still valid
+    in a curved cosmology. Instead a warning is printed informing the user
+    that with Gnuastro's current implementation angular diameter based
+    calculations will be wrong. This was implemented by Boud Roukema.
   --angulardiamdist is the new name for the old '--angulardimdist'
     option. This was necessary to avoid confusion of 'dim' in the old name
     with "dimension"; 'diam' makes it clear that this is a "diameter", not
-    'dimension'.
+    'dimension'. This was suggested by Boud Roukema.
+
 *** Crop
   - The 'ICF*' keywords are now written in the 0-th HDU of Crop's outputs
     (the extension with no data, only metadata). Therefore, the WCS is now
@@ -233,6 +240,14 @@ See the end of the file for license conditions.
     directory that was being used by a later call).
 
 *** Library
+  - gal_cosmology_age: a 'quiet' argument is added to disable warnings.
+  - gal_cosmology_proper_distance: similar 'gal_cosmology_age'.
+  - gal_cosmology_comoving_volume: similar 'gal_cosmology_age'.
+  - gal_cosmology_critical_density: similar 'gal_cosmology_age'.
+  - gal_cosmology_angular_distance: similar 'gal_cosmology_age'.
+  - gal_cosmology_luminosity_distance: similar 'gal_cosmology_age'.
+  - gal_cosmology_distance_modulus: similar 'gal_cosmology_age'.
+  - gal_cosmology_to_absolute_mag: similar 'gal_cosmology_age'.
   - gal_fits_key_write: to generalize, the "title" argument has been
     removed (because it is a keyword). It can also create the file if it
     doesn't exist and can optionally free the list of keywords (until now,
@@ -285,6 +300,9 @@ See the end of the file for license conditions.
     '--frac-max'; reported by Helena Domínguez Sánchez.
   - bug #65194: --quiet does not disable the warning for un-given
     minmapsize; reported by Boud Roukema.
+  - bug #65195: CosmicCalculator prevents non-flat FLRW models and angular
+    diameter distance only correct in flat models. Found and fixed by Boud
+    Roukema.
 
 
 
diff --git a/bin/cosmiccal/cosmiccal.c b/bin/cosmiccal/cosmiccal.c
index 90d11ac9..de6edc5c 100644
--- a/bin/cosmiccal/cosmiccal.c
+++ b/bin/cosmiccal/cosmiccal.c
@@ -72,38 +72,41 @@ cosmiccal_printall(struct cosmiccalparams *p)
   double curage, ccritd, distmod, outage, zcritd;
 
   /* The user wants everything, do all the calculations and print
-     everything with full descriptions. */
+     everything with full descriptions.
+
+     Note that 'quiet' is activated here because the same sanity check is
+     already done in CosmicCalculator's 'ui.c'. */
   curage=gal_cosmology_age(0.0f, p->H0, p->olambda, p->omatter,
-                           p->oradiation);
+                           p->oradiation, 1);
 
   ccritd=gal_cosmology_critical_density(0.0f, p->H0, p->olambda, p->omatter,
-                                        p->oradiation);
+                                        p->oradiation, 1);
 
   pd=gal_cosmology_proper_distance(p->redshift, p->H0, p->olambda, p->omatter,
-                                   p->oradiation);
+                                   p->oradiation, 1);
 
   ad=gal_cosmology_angular_distance(p->redshift, p->H0, p->olambda,
-                                    p->omatter, p->oradiation);
+                                    p->omatter, p->oradiation, 1);
 
   ld=gal_cosmology_luminosity_distance(p->redshift, p->H0, p->olambda,
-                                       p->omatter, p->oradiation);
+                                       p->omatter, p->oradiation, 1);
 
   distmod=gal_cosmology_distance_modulus(p->redshift, p->H0, p->olambda,
-                                         p->omatter, p->oradiation);
+                                         p->omatter, p->oradiation, 1);
 
   absmagconv=gal_cosmology_to_absolute_mag(p->redshift, p->H0, p->olambda,
-                                           p->omatter, p->oradiation);
+                                           p->omatter, p->oradiation, 1);
 
   outage=gal_cosmology_age(p->redshift, p->H0, p->olambda, p->omatter,
-                           p->oradiation);
+                           p->oradiation, 1);
 
   zcritd=gal_cosmology_critical_density(p->redshift, p->H0, p->olambda,
-                                        p->omatter, p->oradiation);
+                                        p->omatter, p->oradiation, 1);
 
   vel=gal_cosmology_velocity_from_z(p->redshift);
 
   vz=gal_cosmology_comoving_volume(p->redshift, p->H0, p->olambda, p->omatter,
-                                   p->oradiation);
+                                   p->oradiation, 1);
 
   /* Print out results: */
   cosmiccal_print_input(p);
@@ -163,7 +166,10 @@ cosmiccal(struct cosmiccalparams *p)
     }
 
   /* In case the user just wants one number, only print that and
-     return. */
+     return.
+
+     Note that 'quiet' is activated here because the same sanity check is
+     already done in CosmicCalculator's 'ui.c'. */
   if(p->specific)
     {
       for(tmp=p->specific;tmp!=NULL;tmp=tmp->next)
@@ -177,14 +183,15 @@ cosmiccal(struct cosmiccalparams *p)
 
             case UI_KEY_AGENOW:
               printf("%f", gal_cosmology_age(0.0f, p->H0, p->olambda,
-                                             p->omatter, p->oradiation));
+                                             p->omatter, p->oradiation, 1));
               break;
 
             case UI_KEY_CRITICALDENSITYNOW:
               printf("%e", gal_cosmology_critical_density(0.0f, p->H0,
                                                           p->olambda,
                                                           p->omatter,
-                                                          p->oradiation));
+                                                          p->oradiation,
+                                                          1));
               break;
 
             case UI_KEY_PROPERDISTANCE:
@@ -192,7 +199,8 @@ cosmiccal(struct cosmiccalparams *p)
                                                          p->H0,
                                                          p->olambda,
                                                          p->omatter,
-                                                         p->oradiation));
+                                                         p->oradiation,
+                                                         1));
               break;
 
             case UI_KEY_ANGULARDIAMDIST:
@@ -200,7 +208,8 @@ cosmiccal(struct cosmiccalparams *p)
                                                           p->H0,
                                                           p->olambda,
                                                           p->omatter,
-                                                          p->oradiation));
+                                                          p->oradiation,
+                                                          1));
               break;
 
             case UI_KEY_ARCSECTANDIST:
@@ -208,7 +217,8 @@ cosmiccal(struct cosmiccalparams *p)
                                                             p->H0,
                                                             p->olambda,
                                                             p->omatter,
-                                                            p->oradiation)
+                                                            p->oradiation,
+                                                            1)
                              * 1000 * M_PI / 3600 / 180 ) );
               break;
 
@@ -217,7 +227,8 @@ cosmiccal(struct cosmiccalparams *p)
                                                              p->H0,
                                                              p->olambda,
                                                              p->omatter,
-                                                             p->oradiation));
+                                                             p->oradiation,
+                                                             1));
               break;
 
             case UI_KEY_DISTANCEMODULUS:
@@ -225,26 +236,28 @@ cosmiccal(struct cosmiccalparams *p)
                                                           p->H0,
                                                           p->olambda,
                                                           p->omatter,
-                                                          p->oradiation));
+                                                          p->oradiation,
+                                                          1));
               break;
 
             case UI_KEY_ABSMAGCONV:
               printf("%f", gal_cosmology_to_absolute_mag(p->redshift, p->H0,
                                                          p->olambda,
                                                          p->omatter,
-                                                         p->oradiation));
+                                                         p->oradiation,
+                                                         1));
               break;
 
             case UI_KEY_AGE:
               printf("%f", gal_cosmology_age(p->redshift, p->H0, p->olambda,
-                                             p->omatter, p->oradiation));
+                                             p->omatter, p->oradiation, 1));
               break;
 
             case UI_KEY_LOOKBACKTIME:
               curage=gal_cosmology_age(0.0f, p->H0, p->olambda, p->omatter,
-                                       p->oradiation);
+                                       p->oradiation, 1);
               zage=gal_cosmology_age(p->redshift, p->H0, p->olambda,
-                                     p->omatter, p->oradiation);
+                                     p->omatter, p->oradiation, 1);
               printf("%f", curage-zage);
               break;
 
@@ -253,14 +266,16 @@ cosmiccal(struct cosmiccalparams *p)
                                                           p->H0,
                                                           p->olambda,
                                                           p->omatter,
-                                                          p->oradiation));
+                                                          p->oradiation,
+                                                          1));
               break;
 
             case UI_KEY_VOLUME:
               printf("%f", gal_cosmology_comoving_volume(p->redshift, p->H0,
                                                          p->olambda,
                                                          p->omatter,
-                                                         p->oradiation));
+                                                         p->oradiation,
+                                                         1));
               break;
 
             case UI_KEY_USEDVELOCITY:
diff --git a/bin/cosmiccal/ui.c b/bin/cosmiccal/ui.c
index 3ca6c35f..45a1ed41 100644
--- a/bin/cosmiccal/ui.c
+++ b/bin/cosmiccal/ui.c
@@ -407,12 +407,13 @@ ui_check_only_options(struct cosmiccalparams *p)
 
   /* Check if the density fractions add up to 1 (within floating point
      error). */
-  if( sum > (1+1e-8) || sum < (1-1e-8) )
-    error(EXIT_FAILURE, 0, "sum of fractional densities is not 1, but "
-          "%.8f. The cosmological constant ('olambda'), matter "
-          "('omatter') and radiation ('oradiation') densities are given "
-          "as %.8f, %.8f, %.8f", sum, p->olambda, p->omatter,
-          p->oradiation);
+  if( (sum > (1+1e-8) || sum < (1-1e-8)) && p->cp.quiet==0)
+    error(EXIT_SUCCESS, 0, "WARNING: non-flat FLRW model: the curvature "
+          "density parameter is %.8f; therefore angular diameter based "
+          "distances like '--angulardiamdist' or '--arcsectandist' will "
+          "be wrong in Gnuastro's current implementation; see "
+          "https://savannah.gnu.org/bugs/?65195. This warning message "
+          "can be disabled with '--quiet'", 1.0-sum);
 
   /* Make sure that '--listlines' and '--listlinesatz' aren't called
      together. */
diff --git a/lib/cosmology.c b/lib/cosmology.c
index b5ad954d..5d21a872 100644
--- a/lib/cosmology.c
+++ b/lib/cosmology.c
@@ -86,7 +86,7 @@ struct cosmology_integrand_t
     sum should not exceed 1. */
 static void
 cosmology_density_check(double o_lambda_0, double o_matter_0,
-                        double o_radiation_0)
+                        double o_radiation_0, int quiet)
 {
   double sum = o_lambda_0 + o_matter_0 + o_radiation_0;
 
@@ -111,13 +111,12 @@ cosmology_density_check(double o_lambda_0, double 
o_matter_0,
 
   /* Check if the density fractions add up to 1 (within floating point
       error). */
-  if( sum > (1+1e-8) || sum < (1-1e-8) )
-    error(EXIT_FAILURE, 0, "sum of fractional densities is not 1, "
-          "but %.8f. The cosmological constant ('olambda'), matter "
-          "('omatter') and radiation ('oradiation') densities are given "
-          "as %.8f, %.8f, %.8f", sum, o_lambda_0, o_matter_0,
-          o_radiation_0);
-
+  if( (sum > (1+1e-8) || sum < (1-1e-8)) && quiet==0 )
+    error(EXIT_SUCCESS, 0, "WARNING: non-flat FLRW model: the curvature "
+          "density parameter is %.8f; therefore angular diameter based "
+          "distances like will be wrong in Gnuastro's current "
+          "implementation; see https://savannah.gnu.org/bugs/?65195. "
+          "This warning message can be disabled with '--quiet'", 1.0-sum);
 }
 
 
@@ -220,7 +219,7 @@ cosmology_integrand_comoving_volume(double z, void *params)
    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)
+                  double o_radiation_0, int quiet)
 {
   gsl_function F;
   double result, error;
@@ -232,7 +231,7 @@ gal_cosmology_age(double z, double H0, double o_lambda_0, 
double o_matter_0,
 
   /* Basic sanity check (no problem with the usage of these variables in
      the definitions above; they are just  */
-  cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0);
+  cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0, quiet);
 
   /* Set the GSL function parameters. */
   F.params=&p;
@@ -250,7 +249,8 @@ gal_cosmology_age(double z, double H0, double o_lambda_0, 
double o_matter_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)
+                              double o_matter_0, double o_radiation_0,
+                              int quiet)
 {
   int status;
   size_t neval;
@@ -265,7 +265,7 @@ gal_cosmology_proper_distance(double z, double H0, double 
o_lambda_0,
 
   /* Basic sanity check (no problem with the usage of these variables in
      the definitions above; they are just  */
-  cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0);
+  cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0, quiet);
 
   /* Set the GSL function parameters */
   F.params=&p;
@@ -306,7 +306,8 @@ gal_cosmology_proper_distance(double z, double H0, double 
o_lambda_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)
+                              double o_matter_0, double o_radiation_0,
+                              int quiet)
 {
   int status;
   size_t neval;
@@ -323,7 +324,7 @@ gal_cosmology_comoving_volume(double z, double H0, double 
o_lambda_0,
 
   /* Basic sanity check (no problem with the usage of these variables in
      the definitions above; they are just  */
-  cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0);
+  cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0, quiet);
 
   /* Set the GSL function parameters */
   F.params=&p;
@@ -364,7 +365,8 @@ gal_cosmology_comoving_volume(double z, double H0, double 
o_lambda_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)
+                               double o_matter_0, double o_radiation_0,
+                               int quiet)
 {
   double H;
   double H0s=H0/1000/GSL_CONST_MKSA_PARSEC;     /* H0 in units of seconds. */
@@ -374,7 +376,7 @@ gal_cosmology_critical_density(double z, double H0, double 
o_lambda_0,
 
   /* Basic sanity check (no problem with the usage of these variables in
      the definitions above; they are just  */
-  cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0);
+  cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0, quiet);
 
   /* Set the place holder, then return the result. */
   H = H0s * cosmology_integrand_Ez(z, &p);
@@ -388,11 +390,16 @@ gal_cosmology_critical_density(double z, double H0, 
double o_lambda_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)
+                               double o_matter_0, double o_radiation_0,
+                               int quiet)
 {
-  cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0);
+  /* Sanity checks. */
+  cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0, quiet);
+
+  /* Note that here 'quiet' is activated because we already do the
+     density check once. */
   return gal_cosmology_proper_distance(z, H0, o_lambda_0, o_matter_0,
-                                       o_radiation_0) / (1+z);
+                                       o_radiation_0, 1) / (1+z);
 }
 
 
@@ -402,11 +409,16 @@ gal_cosmology_angular_distance(double z, double H0, 
double o_lambda_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)
+                                  double o_matter_0, double o_radiation_0,
+                                  int quiet)
 {
-  cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0);
+  /* Sanity checks. */
+  cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0, quiet);
+
+  /* Note that here 'quiet' is activated because we already do the density
+     check once. */
   return gal_cosmology_proper_distance(z, H0, o_lambda_0, o_matter_0,
-                                       o_radiation_0) * (1+z);
+                                       o_radiation_0, 1) * (1+z);
 }
 
 
@@ -416,11 +428,16 @@ gal_cosmology_luminosity_distance(double z, double H0, 
double o_lambda_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)
+                               double o_matter_0, double o_radiation_0,
+                               int quiet)
 {
-  cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0);
+  /* Sanity checks. */
+  cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0, quiet);
+
+  /* Note that here 'quiet' is activated because we already do the density
+     check once. */
   double ld=gal_cosmology_luminosity_distance(z, H0, o_lambda_0, o_matter_0,
-                                              o_radiation_0);
+                                              o_radiation_0, 1);
   return 5*(log10(ld*1000000)-1);
 }
 
@@ -431,11 +448,16 @@ gal_cosmology_distance_modulus(double z, double H0, 
double o_lambda_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)
+                              double o_matter_0, double o_radiation_0,
+                              int quiet)
 {
-  cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0);
+  /* Sanity checks. */
+  cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0, quiet);
+
+  /* Note that here 'quiet' is activated because we already do the density
+     check once. */
   double dm=gal_cosmology_distance_modulus(z, H0, o_lambda_0, o_matter_0,
-                                           o_radiation_0);
+                                           o_radiation_0, 1);
   return dm-2.5*log10(1.0+z);
 }
 
diff --git a/lib/gnuastro/cosmology.h b/lib/gnuastro/cosmology.h
index 7e338142..9f00168c 100644
--- a/lib/gnuastro/cosmology.h
+++ b/lib/gnuastro/cosmology.h
@@ -51,42 +51,49 @@ __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);
+                  double o_radiation_0, int quiet);
 
 /* 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);
+                              double o_matter_0, double o_radiation_0,
+                              int quiet);
 
 /* 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);
+                              double o_matter_0, double o_radiation_0,
+                              int quiet);
 
 /* 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 o_matter_0, double o_radiation_0,
+                               int quiet);
 
 /* 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);
+                               double o_matter_0, double o_radiation_0,
+                               int quiet);
 
 /* 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);
+                                  double o_matter_0, double o_radiation_0,
+                                  int quiet);
 
 /* 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 o_matter_0, double o_radiation_0,
+                               int quiet);
 
 /* 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 o_matter_0, double o_radiation_0,
+                              int quiet);
 
 double
 gal_cosmology_velocity_from_z(double z);



reply via email to

[Prev in Thread] Current Thread [Next in Thread]