[gnuastro-commits] master 0a1b8bd: Library (cosmology.h) and CosmicCalulator: velocity to z conversion

[Mohammad Akhlaghi]

branch: master
commit 0a1b8bde36441951b96b82685fe01a4ef72b6ab8
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Library (cosmology.h) and CosmicCalulator: velocity to z conversion
    
    Until now the only inputs to CosmicCalculator were the redshift and the
    observed wavelength. But in some scenarios, we have the velocity of the
    object instead and although the conversion isn't too complicated, it is
    annoying!
    
    With this commit, it is now possible to give the input redshift of
    CosmicCalculator through the '--velocity' option (value will be
    interpretted to be in units of km/s) and the used velocity in that run can
    be printed out with the new '--usedvelocity' option. Also, in the general
    information printed, we now also print the velocity.
---
 NEWS                      |  4 ++++
 bin/cosmiccal/args.h      | 27 +++++++++++++++++++++++++++
 bin/cosmiccal/cosmiccal.c |  9 ++++++++-
 bin/cosmiccal/main.h      |  3 ++-
 bin/cosmiccal/ui.c        | 32 ++++++++++++++++++++++----------
 bin/cosmiccal/ui.h        |  6 ++++--
 doc/gnuastro.texi         | 30 ++++++++++++++++++++++++++----
 lib/cosmology.c           | 24 ++++++++++++++++++++++++
 lib/gnuastro/cosmology.h  |  6 ++++++
 9 files changed, 123 insertions(+), 18 deletions(-)

diff --git a/NEWS b/NEWS
index a1b57c6..17893cc 100644
--- a/NEWS
+++ b/NEWS
@@ -26,6 +26,8 @@ See the end of the file for license conditions.
      in Table.
 
   CosmicCalculator:
+   --velocity: Velocity (in km/s) to use instead of input redshift.
+   --usedvelocity: Print the velocity (in km/s) at input redshift.
    --listlinesatz: Print the wavelength of all pre-defined spectral lines
      at given redshift as a simple table with the line names. This is very
      convenient and can be used in conjunction with '--obsline' for example
@@ -72,6 +74,8 @@ See the end of the file for license conditions.
    - Spectral lines library: SiIII, OIII, CIV, NV and rest of Lyman series.
    - GAL_CONFIG_HAVE_WCSLIB_DIS_H: if the host's WCSLIB supports distortions.
    - GAL_CONFIG_HAVE_WCSLIB_MJDREF: if the host's WCSLIB reads MJDREF keyword.
+   - gal_cosmology_velocity_from_z: Calculate velocity from redshift.
+   - gal_cosmology_z_from_velocity: Calculate redshift from velocity.
    - gal_data_array_ptr_calloc: Allocate array of pointers to gal_data_t
    - gal_data_array_ptr_free: Free all the datasets within the array and 
itself.
    - gal_fits_key_list_title_add: Add a title key word to the list.
diff --git a/bin/cosmiccal/args.h b/bin/cosmiccal/args.h
index 947098c..12e9eed 100644
--- a/bin/cosmiccal/args.h
+++ b/bin/cosmiccal/args.h
@@ -45,6 +45,19 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET
     },
     {
+      "velocity",
+      UI_KEY_VELOCITY,
+      "FLT",
+      0,
+      "Velocity of interest in km/s.",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->velocity,
+      GAL_TYPE_FLOAT64,
+      GAL_OPTIONS_RANGE_GE_0,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
       "obsline",
       UI_KEY_OBSLINE,
       "STR,FLT",
@@ -303,6 +316,20 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET,
       ui_add_to_single_value,
     },
+    {
+      "usedvelocity",
+      UI_KEY_USEDVELOCITY,
+      0,
+      0,
+      "Used velocity (in km/s) for this run.",
+      UI_GROUP_BASIC,
+      &p->specific,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_add_to_single_value,
+    },
 
 
 
diff --git a/bin/cosmiccal/cosmiccal.c b/bin/cosmiccal/cosmiccal.c
index f63d3ec..70bd485 100644
--- a/bin/cosmiccal/cosmiccal.c
+++ b/bin/cosmiccal/cosmiccal.c
@@ -72,7 +72,7 @@ cosmiccal_print_input(struct cosmiccalparams *p)
 static void
 cosmiccal_printall(struct cosmiccalparams *p)
 {
-  double ad, ld, vz, pd, absmagconv;
+  double ad, ld, vz, pd, vel, absmagconv;
   double curage, ccritd, distmod, outage, zcritd;
 
   /* The user wants everything, do all the calculations and print
@@ -104,6 +104,8 @@ cosmiccal_printall(struct cosmiccalparams *p)
   zcritd=gal_cosmology_critical_density(p->redshift, p->H0, p->olambda,
                                         p->omatter, p->oradiation);
 
+  vel=gal_cosmology_velocity_from_z(p->redshift);
+
   vz=gal_cosmology_comoving_volume(p->redshift, p->H0, p->olambda, p->omatter,
                                    p->oradiation);
 
@@ -115,6 +117,7 @@ cosmiccal_printall(struct cosmiccalparams *p)
   printf(    " ------------\n");
   printf(FLTFORMAT, "Age of Universe now (Ga*):", curage);
   printf(EXPFORMAT, "Critical density now (g/cm^3):",  ccritd);
+  printf(FLTFORMAT, "Velocity at z (km/s):", vel);
   printf(FLTFORMAT, "Proper distance to z (Mpc):", pd);
   printf(FLTFORMAT, "Angular diameter distance to z (Mpc):", ad);
   printf(FLTFORMAT, "Tangential distance covered by 1 arcsec at z (Kpc):",
@@ -256,6 +259,10 @@ cosmiccal(struct cosmiccalparams *p)
                                                           p->oradiation));
               break;
 
+            case UI_KEY_USEDVELOCITY:
+              printf("%g", gal_cosmology_velocity_from_z(p->redshift));
+              break;
+
             case UI_KEY_LINEATZ:
               printf("%g", gal_list_f64_pop(&p->specific_arg)*(1+p->redshift));
               break;
diff --git a/bin/cosmiccal/main.h b/bin/cosmiccal/main.h
index 8fb700c..bb32c4e 100644
--- a/bin/cosmiccal/main.h
+++ b/bin/cosmiccal/main.h
@@ -47,11 +47,12 @@ struct cosmiccalparams
 
   /* Input: */
   double              redshift; /* Redshift of interest.                */
+  gal_data_t          *obsline; /* Observed wavelength of a line.       */
+  double              velocity; /* Velocity of interest.                */
   double                    H0; /* Current expansion rate (km/sec/Mpc). */
   double               olambda; /* Current cosmological constant dens.  */
   double               omatter; /* Current matter density.              */
   double            oradiation; /* Current radiation density.           */
-  gal_data_t          *obsline; /* Observed wavelength of a line.       */
   uint8_t            listlines; /* List the known spectral lines.       */
   uint8_t         listlinesatz; /* List the known spectral lines.       */
 
diff --git a/bin/cosmiccal/ui.c b/bin/cosmiccal/ui.c
index c5d73c5..5d77c93 100644
--- a/bin/cosmiccal/ui.c
+++ b/bin/cosmiccal/ui.c
@@ -116,6 +116,7 @@ ui_initialize_options(struct cosmiccalparams *p,
 
   /* Program specific initializations. */
   p->redshift            = NAN;
+  p->velocity            = NAN;
 
   /* Modify the common options. */
   for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
@@ -378,6 +379,9 @@ ui_parse_obsline(struct argp_option *option, char *arg,
 static void
 ui_read_check_only_options(struct cosmiccalparams *p)
 {
+  int hasobsline=p->obsline!=NULL;
+  int hasredshift=!isnan(p->redshift);
+  int hasvelocity=!isnan(p->velocity);
   double sum = p->olambda + p->omatter + p->oradiation;
 
   /* Check if the density fractions add up to 1 (within floating point
@@ -394,17 +398,19 @@ ui_read_check_only_options(struct cosmiccalparams *p)
     error(EXIT_FAILURE, 0, "'--listlines' and '--listlinesatz' can't be "
           "called together");
 
-  /* Make sure that atleast one of '--redshift' or '--obsline' are given
-     for the running redshift. Except when '--listlines' is used (that
-     doesn't need any redshift). */
-  if(isnan(p->redshift) && p->obsline==NULL && p->listlines==0)
-    error(EXIT_FAILURE, 0, "no redshift specified! Please use '--redshift' "
-          "or '--obsline' to specify a redshift");
+  /* Make sure that atleast one of '--redshift', '--obsline', or velocity
+     are given for the running redshift. Except when '--listlines' is used
+     (that doesn't need any redshift). */
+  if(isnan(p->redshift) && isnan(p->velocity) && p->obsline==NULL
+     && p->listlines==0)
+    error(EXIT_FAILURE, 0, "no redshift/velocity specified! Please use "
+          "'--redshift', '--velocity' (in km/s), or '--obsline' to specify "
+          "a redshift, run with '--help' for more");
 
   /* Make sure that '--redshift' and '--obsline' aren't called together. */
-  if(!isnan(p->redshift) && p->obsline)
-    error(EXIT_FAILURE, 0, "'--redshift' and '--obsline' cannot be called "
-          "together");
+  if( (hasredshift + hasvelocity + hasobsline) > 1 )
+    error(EXIT_FAILURE, 0, "only one of '--redshift', '--velocity', or "
+          "'--obsline' can be called in each run");
 }
 
 
@@ -475,13 +481,19 @@ ui_preparations(struct cosmiccalparams *p)
   if(p->listlines)
     ui_list_lines(p);
 
-  /* If '--obsline' has been given, set the redshift based on it. */
+  /* If '--obsline' has been given, set the redshift based on it (it can't
+     be called with '--velocity'). */
   if(p->obsline)
     p->redshift = ( (p->obsline->status==GAL_SPECLINES_INVALID)
                     ? gal_speclines_line_redshift(obsline[0], obsline[1])
                     : gal_speclines_line_redshift_code(obsline[0],
                                                        p->obsline->status) );
 
+  /* If '--velocity' has been given, set the redshift based on it (it can't
+     be called with '--obsline'). */
+  if( !isnan(p->velocity) )
+    p->redshift = gal_cosmology_z_from_velocity(p->velocity);
+
   /* Currently GSL will fail for z=0. So if a value of zero is given (bug
      #56299). As a work-around, in such cases, we'll change it to an
      extremely small value. NOTE: This has to be after the 'obsline'
diff --git a/bin/cosmiccal/ui.h b/bin/cosmiccal/ui.h
index 142d79f..6f0c034 100644
--- a/bin/cosmiccal/ui.h
+++ b/bin/cosmiccal/ui.h
@@ -43,13 +43,14 @@ enum program_args_groups
 
 /* Available letters for short options:
 
-   f j k n p t w x y
-   B E J Q R W X Y
+   f j k n p t w x
+   B E J Q R W X
 */
 enum option_keys_enum
 {
   /* With short-option version. */
   UI_KEY_REDSHIFT            = 'z',
+  UI_KEY_VELOCITY            = 'y',
   UI_KEY_OBSLINE             = 'O',
   UI_KEY_H0                  = 'H',
   UI_KEY_OLAMBDA             = 'l',
@@ -69,6 +70,7 @@ enum option_keys_enum
   UI_KEY_LOOKBACKTIME        = 'b',
   UI_KEY_CRITICALDENSITY     = 'c',
   UI_KEY_VOLUME              = 'v',
+  UI_KEY_USEDVELOCITY        = 'Y',
   UI_KEY_LINEATZ             = 'i',
 
   /* Only with long version (start with a value 1000, the rest will be set
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index f7e1b37..de352fc 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -17877,8 +17877,17 @@ The inputs to CosmicCalculator can be specified with 
the following options:
 @item -z FLT
 @itemx --redshift=FLT
 The redshift of interest.
-You can also specify the target redshift from spectral lines and their 
observed wavelengths, for that see @option{--obsline}.
-Note that this option cannot be called with @option{--obsline}.
+There are two other ways that you can specify the target redshift:
+1) Spectral lines and their observed wavelengths, see @option{--obsline}.
+2) Velocity, see @option{--velocity}.
+Hence this option cannot be called with @option{--obsline} or 
@option{--velocity}.
+
+@item -y FLT
+@itemx --velocity=FLT
+Input velocity in km/s.
+The given value will be converted to redshift internally, and used in any 
subsequent calculation.
+This option is thus an alternative to @code{--redshift} or @code{--obsline}, 
it cannot be used with them.
+The conversion will be done with the more general and accurate relativistic 
equation of @mymath{1+z=\sqrt{(c+v)/(c-v)}}, not the simplified 
@mymath{z\approx v/c}.
 
 @item -H FLT
 @itemx --H0=FLT
@@ -17901,7 +17910,7 @@ Radiation density divided by the critical density in 
the current Universe (@myma
 @cindex Rest-frame wavelength
 @cindex Wavelength, rest-frame
 Find the redshift to use in next steps based on the rest-frame and observed 
wavelengths of a line.
-This option can be used instead of @option{--redshift}.
+This option is thus an alternative to @code{--redshift} or @code{--velocity}, 
it cannot be used with them.
 Wavelengths are assumed to be in Angstroms.
 The first argument identifies the line.
 It can be one of the standard names below, or any rest-frame wavelength in 
Angstroms.
@@ -18128,7 +18137,13 @@ In case you have forgot the units, you can use the 
@option{--help} option which
 @itemx --usedredshift
 The redshift that was used in this run.
 In many cases this is the main input parameter to CosmicCalculator, but it is 
useful in others.
-For example in combination with @option{--obsline} (where you give an observed 
and rest-frame wavelength and would like to know the redshift), or if you want 
to run CosmicCalculator in a loop while changing the redshift and you want to 
keep the redshift value.
+For example in combination with @option{--obsline} (where you give an observed 
and rest-frame wavelength and would like to know the redshift) or with 
@option{--velocity} (where you specify the velocity instead of redshift).
+Another example is when you run CosmicCalculator in a loop, while changing the 
redshift and you want to keep the redshift value with the resulting calculation.
+
+@item -Y
+@itemx --usedvelocity
+The velocity (in km/s) that was used in this run.
+The conversion from redshift will be done with the more general and accurate 
relativistic equation of @mymath{1+z=\sqrt{(c+v)/(c-v)}}, not the simplified 
@mymath{z\approx v/c}.
 
 @item -G
 @itemx --agenow
@@ -25783,6 +25798,13 @@ 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}.
 @end deftypefun
 
+@deftypefun double gal_cosmology_velocity_from_z (double @code{z})
+Return the velocity (in km/s) corresponding to the given redshift (@code{z}).
+@end deftypefun
+
+@deftypefun double gal_cosmology_z_from_velocity (double @code{v})
+Return the redshift corresponding to the given velocity (@code{v} in km/s).
+@end deftypefun
 
 @node Library demo programs,  , Gnuastro library, Library
 @section Library demo programs
diff --git a/lib/cosmology.c b/lib/cosmology.c
index f02261d..873b6b3 100644
--- a/lib/cosmology.c
+++ b/lib/cosmology.c
@@ -310,3 +310,27 @@ gal_cosmology_to_absolute_mag(double z, double H0, double 
o_lambda_0,
                                            o_radiation_0);
   return dm-2.5*log10(1.0+z);
 }
+
+
+
+
+
+/* Velocity at given redshift in units of km/s. */
+double
+gal_cosmology_velocity_from_z(double z)
+{
+  double c=GSL_CONST_MKSA_SPEED_OF_LIGHT;
+  return c * ( (1+z)*(1+z) - 1 ) / ( (1+z)*(1+z) + 1 ) / 1000;
+}
+
+
+
+
+
+/* Redshift at given velocity (in units of km/s). */
+double
+gal_cosmology_z_from_velocity(double v)
+{
+  double c=GSL_CONST_MKSA_SPEED_OF_LIGHT/1000;
+  return sqrt( (c+v)/(c-v) ) - 1;
+}
diff --git a/lib/gnuastro/cosmology.h b/lib/gnuastro/cosmology.h
index e49f607..9b48fb8 100644
--- a/lib/gnuastro/cosmology.h
+++ b/lib/gnuastro/cosmology.h
@@ -88,6 +88,12 @@ double
 gal_cosmology_to_absolute_mag(double z, double H0, double o_lambda_0,
                               double o_matter_0, double o_radiation_0);
 
+double
+gal_cosmology_velocity_from_z(double z);
+
+double
+gal_cosmology_z_from_velocity(double v);
+
 __END_C_DECLS    /* From C++ preparations */
 
 #endif           /* __GAL_COSMOLOGY_H__ */