[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 3b6c15d 036/113: Merged with recent work in ma
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 3b6c15d 036/113: Merged with recent work in master, no conflicts |
|
Date: |
Fri, 16 Apr 2021 10:33:38 -0400 (EDT) |
branch: master
commit 3b6c15d1574e33e9e04096951835e26cc42ba285
Merge: 06798a5 49d6c034
Author: Mohammad Akhlaghi <akhlaghi@gnu.org>
Commit: Mohammad Akhlaghi <akhlaghi@gnu.org>
Merged with recent work in master, no conflicts
All the recent work in master is now merged with this branch. There were no
conflicts.
---
.mailmap | 1 +
NEWS | 29 ++
bin/arithmetic/arithmetic.c | 428 ++++++++++++++--
bin/arithmetic/arithmetic.h | 15 +
bin/cosmiccal/cosmiccal.c | 240 ++-------
bin/cosmiccal/cosmiccal.h | 5 -
bin/cosmiccal/main.h | 5 -
bin/cosmiccal/ui.c | 11 +-
bin/mkcatalog/args.h | 68 ++-
bin/mkcatalog/columns.c | 181 +++++--
bin/mkcatalog/main.h | 4 +-
bin/mkcatalog/mkcatalog.c | 31 +-
bin/mkcatalog/ui.c | 51 +-
bin/mkcatalog/ui.h | 18 +-
doc/gnuastro.texi | 553 +++++++++++++--------
doc/plotsrc/Makefile | 2 +-
doc/plotsrc/all.tex | 2 +-
.../tex/{sphericalplane.tex => sphereandplane.tex} | 0
lib/Makefile.am | 25 +-
lib/arithmetic.c | 2 +
lib/cosmology.c | 312 ++++++++++++
lib/gnuastro/arithmetic.h | 2 +
lib/gnuastro/cosmology.h | 96 ++++
lib/statistics.c | 2 +-
lib/wcs.c | 11 +-
25 files changed, 1519 insertions(+), 575 deletions(-)
diff --git a/.mailmap b/.mailmap
index 1da35be..bccc8e5 100644
--- a/.mailmap
+++ b/.mailmap
@@ -1 +1,2 @@
+Boud Roukema <boud@cosmo.torun.pl>
<akhlaghi@gnu.org> <makhlaghi@gmail.com>
\ No newline at end of file
diff --git a/NEWS b/NEWS
index e4f4b4b..af7058f 100644
--- a/NEWS
+++ b/NEWS
@@ -8,6 +8,10 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
All programs: a value of `0' to the `--numthreads' option will use the
number of threads available to the system at run time.
+ Arithmetic: The new operators `filter-median' and `filter-mean' can be
+ used to filter (smooth) the input. The size of the filter can be set as
+ the other operands to these operators.
+
BuildProgram: The new `--la' option allows the identification of a
different Libtool `.la' file for Libtool linking information.
@@ -16,6 +20,20 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
keywords in that HDU will be printed (as if only `--printallkeys' was
called).
+ MakeCatalog: physical nature agnostic WCS column names. Previously the
+ first WCS axis was always assumed to be RA and the second DEC. So for
+ example even if you had a spectrum (with X and wavelength as the two WCS
+ dimensions), you would have to ask for `--ra' and `--dec'. The new `--w1'
+ and `--w2' options are now generic and don't assume any particular type
+ only their order in the FITS header. MakeCatalog now also uses the CTYPE
+ and CUNIT keywords to set the names and units of its output columns. The
+ `--ra' and `--dec' options are now just internal aliases for `--w1' or
+ `--w2' which will be determined based on the input's CTYPE keyword. Also
+ the new `--geow1', `--geow2', `--clumpsw1', `--clumpsw2',
+ `--clumpsgeow1', `--clumpsgeow2' options replace the old options
+ `--geora', `--geodec', `--clumpsra', `--clumpsdec', `--clumpsgeora',
+ `--clumpsgeodec'. No alias is currently defined for the latter group.
+
NoiseChisel: with the new `--widekernel' option it is now possible to use
a wider kernel to identify which tiles contain signal. The rest of the
steps (identifying the quantile threshold on the selected tiles and etc)
@@ -36,11 +54,19 @@ 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.
** Removed features
+ MakeCatalog: `--zeropoint' option doesn't have a short option name any
+ more. Previously it was `-z' which was confusing because `-x' and `-y'
+ were used to refer to image coordinate positions.
+
NoiseChisel: The `--dilate' and `--dilatengb' options have been
removed. Growing of true detections is no longer done through dilation
but through the `--detgrowquant' and `--detgrowmaxholesize' options (see
@@ -62,6 +88,9 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
to limit the range of header keywords to read the WCS, similar to how
they are used in `gal_wcs_read'.
+ `gal_wcs_pixel_area_arcsec2' will return NaN (instead of aborting) when
+ input is unreasonable (not two dimensions or not in units of degrees).
+
** Bug fixes
ConvertType crash when changing values (bug #52010).
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index b213e7d..36ed627 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -31,6 +31,9 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/wcs.h>
#include <gnuastro/fits.h>
+#include <gnuastro/threads.h>
+#include <gnuastro/dimension.h>
+#include <gnuastro/statistics.h>
#include <gnuastro/arithmetic.h>
#include <gnuastro-internal/checkset.h>
@@ -60,7 +63,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
CTYPE a=*(CTYPE *)(data->array); if(a>0) return a; }
static size_t
-set_number_of_operands(struct arithmeticparams *p, gal_data_t *data,
+pop_number_of_operands(struct arithmeticparams *p, gal_data_t *data,
char *token_string)
{
/* Check if its a number. */
@@ -121,6 +124,313 @@ set_number_of_operands(struct arithmeticparams *p,
gal_data_t *data,
+/**********************************************************************/
+/**************** Filtering operators *****************/
+/**********************************************************************/
+#define ARITHMETIC_FILTER_DIM 10
+
+struct arithmetic_filter_p
+{
+ int operator; /* The type of filtering. */
+ size_t *fsize; /* Filter size. */
+ size_t *hpfsize; /* Positive Half-filter size. */
+ size_t *hnfsize; /* Negative Half-filter size. */
+ gal_data_t *input; /* Input dataset. */
+ gal_data_t *out; /* Output dataset. */
+
+ int hasblank; /* If the dataset has blank values.*/
+};
+
+
+
+
+
+/* Main filtering work function. */
+static void *
+arithmetic_filter(void *in_prm)
+{
+ struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
+ struct arithmetic_filter_p *afp=(struct arithmetic_filter_p *)tprm->params;
+ gal_data_t *input=afp->input;
+
+ size_t ind,index;
+ int out_type_checked=0;
+ gal_data_t *result=NULL;
+ size_t *hpfsize=afp->hpfsize, *hnfsize=afp->hnfsize;
+ size_t *tsize, *dsize=input->dsize, *fsize=afp->fsize;
+ size_t i, j, coord[ARITHMETIC_FILTER_DIM], ndim=input->ndim;
+ size_t start[ARITHMETIC_FILTER_DIM], end[ARITHMETIC_FILTER_DIM];
+ gal_data_t *tile=gal_data_alloc(NULL, input->type, ndim, afp->fsize, NULL,
+ 0, -1, NULL, NULL, NULL);
+
+ /* Prepare the tile. */
+ free(tile->array);
+ tsize=tile->dsize;
+ tile->block=input;
+
+
+ /* Go over all the pixels that were assigned to this thread. */
+ for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
+ {
+ /* Get the coordinate of the pixel. */
+ ind=tprm->indexs[i];
+ gal_dimension_index_to_coord(ind, ndim, dsize, coord);
+
+ /* See which dimensions need trimming. */
+ tile->size=1;
+ for(j=0;j<ndim;++j)
+ {
+ /* Estimate the coordinate of the filter's starting point. Note
+ that we are dealing with size_t (unsigned int) type here, so
+ there are no negatives. A negative result will produce an
+ extremely large number, so instead of checking for negative,
+ we can just see if the result of a subtraction is less than
+ the width of the input. */
+ if( (coord[j] - hnfsize[j] > dsize[j])
+ || (coord[j] + hpfsize[j] >= dsize[j]) )
+ {
+ start[j] = ( (coord[j] - hnfsize[j] > dsize[j])
+ ? 0 : coord[j] - hnfsize[j] );
+ end[j] = ( (coord[j] + hpfsize[j] >= dsize[j])
+ ? dsize[j]
+ : coord[j] + hpfsize[j] + 1);
+ tsize[j] = end[j] - start[j];
+ }
+ else /* We are NOT on the edge (given requested filter width). */
+ {
+ tsize[j] = fsize[j];
+ start[j] = coord[j] - hnfsize[j];
+ }
+ tile->size *= tsize[j];
+ }
+
+ /* Set the tile's starting pointer. */
+ index=gal_dimension_coord_to_index(ndim, dsize, start);
+ tile->array=gal_data_ptr_increment(input->array, index, input->type);
+
+ /* Do the necessary calculation. */
+ switch(afp->operator)
+ {
+ case ARITHMETIC_OP_FILTER_MEDIAN:
+ result=gal_statistics_median(tile, 0);
+ break;
+
+ case ARITHMETIC_OP_FILTER_MEAN:
+ result=gal_statistics_mean(tile);
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. `afp->operator' code %d is not "
+ "recognized", PACKAGE_BUGREPORT, __func__, afp->operator);
+ }
+
+ /* Make sure the output array type and result's type are the same. We
+ only need to do this once, but we'll suffice to once for each
+ thread for simplicify of the code, it is too negligible to have
+ any real affect. */
+ if( out_type_checked==0)
+ {
+ if(result->type!=afp->out->type )
+ error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s so "
+ "we can address the problem. The tyes of `result' and "
+ "`out' aren't the same, they are respectively: `%s' and "
+ "`%s'", __func__, PACKAGE_BUGREPORT,
+ gal_type_name(result->type, 1),
+ gal_type_name(afp->out->type, 1));
+ out_type_checked=1;
+ }
+
+ /* Copy the result into the output array. */
+ memcpy(gal_data_ptr_increment(afp->out->array, ind, afp->out->type),
+ result->array, gal_type_sizeof(afp->out->type));
+
+ /* Clean up. */
+ gal_data_free(result);
+ }
+
+
+ /* Clean up. */
+ tile->array=NULL;
+ tile->block=NULL;
+ gal_data_free(tile);
+
+
+ /* Wait for all the other threads to finish, then return. */
+ if(tprm->b) pthread_barrier_wait(tprm->b);
+ return NULL;
+}
+
+
+
+
+
+static void
+wrapper_for_filter(struct arithmeticparams *p, char *token, int operator)
+{
+ size_t i=0, ndim, one=1;
+ int type=GAL_TYPE_INVALID;
+ struct arithmetic_filter_p afp={0};
+ size_t fsize[ARITHMETIC_FILTER_DIM];
+ gal_data_t *tmp, *tmp2, *zero, *comp, *fsize_list=NULL;
+ size_t hnfsize[ARITHMETIC_FILTER_DIM], hpfsize[ARITHMETIC_FILTER_DIM];
+
+
+ /* Get the input's number of dimensions. */
+ afp.input=operands_pop(p, token);
+ afp.operator=operator;
+ ndim=afp.input->ndim;
+ afp.hnfsize=hnfsize;
+ afp.hpfsize=hpfsize;
+ afp.fsize=fsize;
+
+
+ /* A small sanity check. */
+ if(ndim>ARITHMETIC_FILTER_DIM)
+ error(EXIT_FAILURE, 0, "%s: currently only datasets with less than "
+ "%d dimensions are acceptable. The input has %zu dimensions",
+ __func__, ARITHMETIC_FILTER_DIM, ndim);
+
+
+ /* A zero value for checking the value of input widths. */
+ zero=gal_data_alloc(NULL, GAL_TYPE_INT32, 1, &one, NULL, 1, -1, NULL,
+ NULL, NULL);
+
+
+ /* Based on the dimensions of the popped operand, pop the necessary
+ number of operands. */
+ for(i=0;i<ndim;++i)
+ gal_list_data_add(&fsize_list, operands_pop(p, token));
+
+
+ /* Make sure the filter size only has single values. */
+ i=0;
+ for(tmp=fsize_list; tmp!=NULL; tmp=tmp->next)
+ {
+ ++i;
+ if(tmp->size!=1)
+ error(EXIT_FAILURE, 0, "the filter length values given to the "
+ "filter operators can only be numbers. Value number %zu has "
+ "%zu elements, so its an array", ndim-i-1, tmp->size);
+ }
+
+
+ /* If the input only has one element, filtering makes no sense, so don't
+ waste time, just add the input onto the stack. */
+ if(afp.input->size==1) afp.out=afp.input;
+ else
+ {
+ /* Allocate an array for the size of the filter and fill it in. The
+ values must be written in the inverse order since the user gives
+ dimensions with the FITS standard. */
+ i=ndim-1;
+ for(tmp=fsize_list; tmp!=NULL; tmp=tmp->next)
+ {
+ /* Make sure the user has given an integer type. */
+ if(tmp->type==GAL_TYPE_FLOAT32 || tmp->type==GAL_TYPE_FLOAT64)
+ error(EXIT_FAILURE, 0, "lengths of filter along dimensions "
+ "must be integer values, not floats. The given length "
+ "along dimension %zu is a float", ndim-i);
+
+ /* Make sure it isn't negative. */
+ comp=gal_arithmetic(GAL_ARITHMETIC_OP_GT, 0, tmp, zero);
+ if( *(uint8_t *)(comp->array) == 0 )
+ error(EXIT_FAILURE, 0, "lengths of filter along dimensions "
+ "must be positive. The given length in dimension %zu"
+ "is either zero or negative", ndim-i);
+
+ /* Convert the input into size_t and put it into the array that
+ keeps the filter size. */
+ tmp2=gal_data_copy_to_new_type(tmp, GAL_TYPE_SIZE_T);
+ fsize[ i ] = *(size_t *)(tmp2->array);
+ gal_data_free(tmp2);
+
+ /* If the width is larger than the input's size, change the width
+ to the input's size. */
+ if( fsize[i] > afp.input->dsize[i] )
+ error(EXIT_FAILURE, 0, "%s: the filter size along dimension %zu "
+ "(%zu) is greater than the input's length in that "
+ "dimension (%zu)", __func__, i, fsize[i],
+ afp.input->dsize[i]);
+
+ /* Go onto the previous dimension. */
+ --i;
+ }
+
+
+ /* Set the half filter sizes. Note that when the size is an odd
+ number, the number of pixels before and after the actual pixel are
+ equal, but for an even number, we will look into one element more
+ when looking before than the ones after. */
+ for(i=0;i<ndim;++i)
+ {
+ if( fsize[i]%2 )
+ hnfsize[i]=hpfsize[i]=fsize[i]/2;
+ else
+ { hnfsize[i]=fsize[i]/2; hpfsize[i]=fsize[i]/2-1; }
+ }
+
+
+ /* See if the input has blank pixels. */
+ afp.hasblank=gal_blank_present(afp.input, 1);
+
+ /* Set the type of the output dataset. */
+ switch(operator)
+ {
+ case ARITHMETIC_OP_FILTER_MEDIAN:
+ type=afp.input->type;
+ break;
+
+ case ARITHMETIC_OP_FILTER_MEAN:
+ type=GAL_TYPE_FLOAT64;
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s to fix "
+ "the problem. The `operator' code %d is not recognized",
+ PACKAGE_BUGREPORT, __func__, operator);
+ }
+
+ /* Allocate the output dataset. Note that filtering doesn't change
+ the units of the dataset. */
+ afp.out=gal_data_alloc(NULL, type, ndim, afp.input->dsize,
+ afp.input->wcs, 0, afp.input->minmapsize,
+ NULL, afp.input->unit, NULL);
+
+ /* Spin off threads for each pixel. */
+ gal_threads_spin_off(arithmetic_filter, &afp, afp.input->size,
+ p->cp.numthreads);
+ }
+
+
+ /* Add the output to the top of the stack. */
+ operands_add(p, NULL, afp.out);
+
+
+ /* Clean up and add the output on top of the stack */
+ gal_data_free(afp.input);
+ gal_list_data_free(fsize_list);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
/***************************************************************/
/************* Reverse Polish algorithm *************/
/***************************************************************/
@@ -274,6 +584,13 @@ reversepolish(struct arithmeticparams *p)
else if (!strcmp(token->v, "float64"))
{ op=GAL_ARITHMETIC_OP_TO_FLOAT64; nop=1; }
+ /* Filters. */
+ else if (!strcmp(token->v, "filter-median"))
+ { op=ARITHMETIC_OP_FILTER_MEDIAN; nop=0; }
+ else if (!strcmp(token->v, "filter-mean"))
+ { op=ARITHMETIC_OP_FILTER_MEAN; nop=0; }
+
+
/* Finished checks with known operators */
else
error(EXIT_FAILURE, 0, "the argument \"%s\" could not be "
@@ -281,51 +598,76 @@ reversepolish(struct arithmeticparams *p)
token->v);
- /* Pop the necessary number of operators. Note that the operators
- are poped from a linked list (which is last-in-first-out). So
- for the operators which need a specific order, the first poped
- operand is actally the last (right most, in in-fix notation)
- input operand.*/
- switch(nop)
+ /* See if the arithmetic library must be called or not. */
+ if(nop)
{
- case 1:
- d1=operands_pop(p, token->v);
- break;
-
- case 2:
- d2=operands_pop(p, token->v);
- d1=operands_pop(p, token->v);
- break;
-
- case 3:
- d3=operands_pop(p, token->v);
- d2=operands_pop(p, token->v);
- d1=operands_pop(p, token->v);
- break;
-
- case -1:
- /* This case is when the number of operands is itself an
- operand. So the first popped operand must be an integer
- number, we will use that to construct a linked list of any
- number of operands within the single `d1' pointer. */
- d1=NULL;
- numop=set_number_of_operands(p, operands_pop(p, token->v),
- token->v);
- for(i=0;i<numop;++i)
- gal_list_data_add(&d1, operands_pop(p, token->v));
- break;
-
- default:
- error(EXIT_FAILURE, 0, "No operators need %d operands", nop);
+ /* Pop the necessary number of operators. Note that the
+ operators are poped from a linked list (which is
+ last-in-first-out). So for the operators which need a
+ specific order, the first poped operand is actally the
+ last (right most, in in-fix notation) input operand.*/
+ switch(nop)
+ {
+ case 1:
+ d1=operands_pop(p, token->v);
+ break;
+
+ case 2:
+ d2=operands_pop(p, token->v);
+ d1=operands_pop(p, token->v);
+ break;
+
+ case 3:
+ d3=operands_pop(p, token->v);
+ d2=operands_pop(p, token->v);
+ d1=operands_pop(p, token->v);
+ break;
+
+ case -1:
+ /* This case is when the number of operands is itself an
+ operand. So the first popped operand must be an
+ integer number, we will use that to construct a linked
+ list of any number of operands within the single `d1'
+ pointer. */
+ d1=NULL;
+ numop=pop_number_of_operands(p, operands_pop(p, token->v),
+ token->v);
+ for(i=0;i<numop;++i)
+ gal_list_data_add(&d1, operands_pop(p, token->v));
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "no operators, `%s' needs %d "
+ "operand(s)", token->v, nop);
+ }
+
+
+ /* Run the arithmetic operation. Note that `gal_arithmetic'
+ is a variable argument function (like printf). So the
+ number of arguments it uses depend on the operator. So
+ when the operator doesn't need three operands, the extra
+ arguments will be ignored. */
+ operands_add(p, NULL, gal_arithmetic(op, flags, d1, d2, d3));
}
-
- /* Run the arithmetic operation. Note that `gal_arithmetic' is a
- variable argument function (like printf). So the number of
- arguments it uses depend on the operator. So when the operator
- doesn't need three operands, the extra arguments will be
- ignored. */
- operands_add(p, NULL, gal_arithmetic(op, flags, d1, d2, d3));
+ /* No need to call the arithmetic library, call the proper
+ wrappers directly. */
+ else
+ {
+ switch(op)
+ {
+ case ARITHMETIC_OP_FILTER_MEAN:
+ case ARITHMETIC_OP_FILTER_MEDIAN:
+ wrapper_for_filter(p, token->v, op);
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! please contact us at "
+ "%s to fix the problem. The code %d is not "
+ "recognized for `op'", __func__, PACKAGE_BUGREPORT,
+ op);
+ }
+ }
}
}
diff --git a/bin/arithmetic/arithmetic.h b/bin/arithmetic/arithmetic.h
index e7f2ea0..833c5bf 100644
--- a/bin/arithmetic/arithmetic.h
+++ b/bin/arithmetic/arithmetic.h
@@ -23,6 +23,21 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#ifndef IMGARITH_H
#define IMGARITH_H
+
+#include <gnuastro/arithmetic.h>
+
+
+/* These are operator codes for functions that aren't in the arithmetic
+ library. */
+enum arithmetic_prog_operators
+{
+ ARITHMETIC_OP_FILTER_MEDIAN=GAL_ARITHMETIC_OP_LAST_CODE,
+ ARITHMETIC_OP_FILTER_MEAN,
+};
+
+
+
+
void
imgarith(struct arithmeticparams *p);
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/bin/mkcatalog/args.h b/bin/mkcatalog/args.h
index 94e841b..30db17b 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -534,7 +534,7 @@ struct argp_option program_options[] =
UI_KEY_RA,
0,
0,
- "Right ascension of flux weighted center.",
+ "Flux weighted center right ascension.",
ARGS_GROUP_COLUMNS,
0,
GAL_TYPE_INVALID,
@@ -548,7 +548,7 @@ struct argp_option program_options[] =
UI_KEY_DEC,
0,
0,
- "Declination of flux weighted center.",
+ "Flux weighted center declination.",
ARGS_GROUP_COLUMNS,
0,
GAL_TYPE_INVALID,
@@ -558,11 +558,11 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
- "geora",
- UI_KEY_GEORA,
+ "w1",
+ UI_KEY_W1,
0,
0,
- "Right ascension of geometric center.",
+ "Flux weighted center in first WCS axis.",
ARGS_GROUP_COLUMNS,
0,
GAL_TYPE_INVALID,
@@ -572,11 +572,11 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
- "geodec",
- UI_KEY_GEODEC,
+ "w2",
+ UI_KEY_W2,
0,
0,
- "Declination of geometric center.",
+ "Flux weighted center in second WCS axis.",
ARGS_GROUP_COLUMNS,
0,
GAL_TYPE_INVALID,
@@ -586,11 +586,11 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
- "clumpsra",
- UI_KEY_CLUMPSRA,
+ "geow1",
+ UI_KEY_GEOW1,
0,
0,
- "Right ascension of f.wht center of all clumps.",
+ "Geometric center in first WCS axis.",
ARGS_GROUP_COLUMNS,
0,
GAL_TYPE_INVALID,
@@ -600,11 +600,11 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
- "clumpsdec",
- UI_KEY_CLUMPSDEC,
+ "geow2",
+ UI_KEY_GEOW2,
0,
0,
- "Declination of f.wht center of all clumps.",
+ "Geometric center in second WCS axis.",
ARGS_GROUP_COLUMNS,
0,
GAL_TYPE_INVALID,
@@ -614,11 +614,11 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
- "clumpsgeora",
- UI_KEY_CLUMPSGEORA,
+ "clumpsw1",
+ UI_KEY_CLUMPSW1,
0,
0,
- "Right ascension of geo. center of all clumps.",
+ "Flux.wht center of all clumps in 1st WCS.",
ARGS_GROUP_COLUMNS,
0,
GAL_TYPE_INVALID,
@@ -628,11 +628,39 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
- "clumpsgeodec",
- UI_KEY_CLUMPSGEODEC,
+ "clumpsw2",
+ UI_KEY_CLUMPSW2,
0,
0,
- "Declination of geometric center of all clumps.",
+ "Flux.wht center of all clumps in 2nd WCS.",
+ ARGS_GROUP_COLUMNS,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
+ "clumpsgeow1",
+ UI_KEY_CLUMPSGEOW1,
+ 0,
+ 0,
+ "Geometric center of all clumps in 1st WCS.",
+ ARGS_GROUP_COLUMNS,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
+ "clumpsgeow2",
+ UI_KEY_CLUMPSGEOW2,
+ 0,
+ 0,
+ "Geometric center of all clumps in 2nd WCS.",
ARGS_GROUP_COLUMNS,
0,
GAL_TYPE_INVALID,
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index 809f673..a3e08f4 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -27,8 +27,11 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <errno.h>
#include <error.h>
#include <stdlib.h>
+#include <string.h>
#include <pthread.h>
+#include <gnuastro-internal/checkset.h>
+
#include "main.h"
#include "mkcatalog.h"
@@ -196,6 +199,74 @@ columns_alloc_clumpsgeoradec(struct mkcatalogparams *p)
/******************************************************************/
/********** Column definition/allocation ***************/
/******************************************************************/
+static void
+columns_wcs_preparation(struct mkcatalogparams *p)
+{
+ size_t i;
+ gal_list_i32_t *colcode;
+ int continue_wcs_check=1;
+
+ /* Make sure a WCS structure is present if we need it. */
+ for(colcode=p->columnids; colcode!=NULL; colcode=colcode->next)
+ {
+ if(continue_wcs_check)
+ {
+ switch(colcode->v)
+ {
+ /* High-level. */
+ case UI_KEY_RA:
+ case UI_KEY_DEC:
+
+ /* Low-level. */
+ case UI_KEY_W1:
+ case UI_KEY_W2:
+ case UI_KEY_GEOW1:
+ case UI_KEY_GEOW2:
+ case UI_KEY_CLUMPSW1:
+ case UI_KEY_CLUMPSW2:
+ case UI_KEY_CLUMPSGEOW1:
+ case UI_KEY_CLUMPSGEOW2:
+ if(p->input->wcs)
+ continue_wcs_check=0;
+ else
+ error(EXIT_FAILURE, 0, "input doesn't have WCS meta-data for "
+ "defining world coordinates (like RA and Dec). Atleast "
+ "one of the requested columns requires this "
+ "information");
+ break;
+ }
+ }
+ else
+ break;
+ }
+
+ /* Convert the high-level WCS columns to low-level ones. */
+ for(colcode=p->columnids; colcode!=NULL; colcode=colcode->next)
+ switch(colcode->v)
+ {
+ case UI_KEY_RA:
+ case UI_KEY_DEC:
+ /* Check all the CTYPES. */
+ for(i=0;i<p->input->ndim;++i)
+ if( !strcmp(p->ctype[i], colcode->v==UI_KEY_RA ? "RA" : "DEC") )
+ {
+ colcode->v = i==0 ? UI_KEY_W1 : UI_KEY_W2;
+ break;
+ }
+
+ /* Make sure it actually existed. */
+ if(i==p->input->ndim)
+ error(EXIT_FAILURE, 0, "%s (hdu: %s): %s not present in any of "
+ "the WCS axis types (CTYPE)", p->inputname, p->cp.hdu,
+ colcode->v==UI_KEY_RA ? "RA" : "DEC");
+ break;
+ }
+}
+
+
+
+
+
/* Set the necessary parameters for each output column and allocate the
space necessary to keep the values. */
void
@@ -207,6 +278,13 @@ columns_define_alloc(struct mkcatalogparams *p)
char *name=NULL, *unit=NULL, *ocomment=NULL, *ccomment=NULL;
uint8_t otype=GAL_TYPE_INVALID, ctype=GAL_TYPE_INVALID, *oiflag, *ciflag;
+ /* If there is any columns that need WCS, the input image needs to have a
+ WCS in its headers. So before anything, we need to check if a WCS is
+ present or not. This can't be done after the initial setting of column
+ properties because the WCS-related columns use information that is
+ based on it (for units and names). */
+ columns_wcs_preparation(p);
+
/* Allocate the array for which intermediate parameters are
necessary. The basic issue is that higher-level calculations require a
smaller domain of raw measurements. So to avoid having to calculate
@@ -427,10 +505,10 @@ columns_define_alloc(struct mkcatalogparams *p)
oiflag[ OCOL_C_GY ] = 1;
break;
- case UI_KEY_RA:
- name = "RA";
- unit = "degrees";
- ocomment = "Flux weighted center right ascension.";
+ case UI_KEY_W1:
+ name = p->ctype[0];
+ unit = p->input->wcs->cunit[0];
+ ocomment = "Flux weighted center in first WCS axis.";
ccomment = ocomment;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_FLOAT64;
@@ -442,10 +520,10 @@ columns_define_alloc(struct mkcatalogparams *p)
columns_alloc_radec(p);
break;
- case UI_KEY_DEC:
- name = "DEC";
- unit = "degrees";
- ocomment = "Flux weighted center declination.";
+ case UI_KEY_W2:
+ name = p->ctype[1];
+ unit = p->input->wcs->cunit[1];
+ ocomment = "Flux weighted center in second WCS axis.";
ccomment = ocomment;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_FLOAT64;
@@ -457,10 +535,10 @@ columns_define_alloc(struct mkcatalogparams *p)
columns_alloc_radec(p);
break;
- case UI_KEY_GEORA:
- name = "GEO_RA";
- unit = "degrees";
- ocomment = "Geometric center right ascension.";
+ case UI_KEY_GEOW1:
+ name = gal_checkset_malloc_cat("GEO_", p->ctype[0]);
+ unit = p->input->wcs->cunit[0];
+ ocomment = "Geometric center in first WCS axis.";
ccomment = ocomment;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_FLOAT64;
@@ -474,10 +552,10 @@ columns_define_alloc(struct mkcatalogparams *p)
columns_alloc_georadec(p);
break;
- case UI_KEY_GEODEC:
- name = "GEO_DEC";
- unit = "degrees";
- ocomment = "Geometric center declination.";
+ case UI_KEY_GEOW2:
+ name = gal_checkset_malloc_cat("GEO_", p->ctype[1]);
+ unit = p->input->wcs->cunit[1];
+ ocomment = "Geometric center in second WCS axis.";
ccomment = ocomment;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_FLOAT64;
@@ -491,10 +569,10 @@ columns_define_alloc(struct mkcatalogparams *p)
columns_alloc_georadec(p);
break;
- case UI_KEY_CLUMPSRA:
- name = "CLUMPS_RA";
- unit = "degrees";
- ocomment = "RA of all clumps flux weighted center.";
+ case UI_KEY_CLUMPSW1:
+ name = gal_checkset_malloc_cat("CLUMPS_", p->ctype[0]);
+ unit = p->input->wcs->cunit[0];
+ ocomment = "Flux.wht center of all clumps in 1st WCS axis.";
ccomment = NULL;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_INVALID;
@@ -506,10 +584,10 @@ columns_define_alloc(struct mkcatalogparams *p)
columns_alloc_clumpsradec(p);
break;
- case UI_KEY_CLUMPSDEC:
- name = "CLUMPS_DEC";
- unit = "degrees";
- ocomment = "Declination of all clumps flux weighted center.";
+ case UI_KEY_CLUMPSW2:
+ name = gal_checkset_malloc_cat("CLUMPS_", p->ctype[1]);
+ unit = p->input->wcs->cunit[1];
+ ocomment = "Flux.wht center of all clumps in 2nd WCS axis.";
ccomment = NULL;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_INVALID;
@@ -521,10 +599,10 @@ columns_define_alloc(struct mkcatalogparams *p)
columns_alloc_clumpsradec(p);
break;
- case UI_KEY_CLUMPSGEORA:
- name = "CLUMPS_RA";
- unit = "degrees";
- ocomment = "RA of all clumps geometric center.";
+ case UI_KEY_CLUMPSGEOW1:
+ name = gal_checkset_malloc_cat("CLUMPS_GEO", p->ctype[0]);
+ unit = p->input->wcs->cunit[0];
+ ocomment = "Geometric center of all clumps in 1st WCS axis.";
ccomment = NULL;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_INVALID;
@@ -536,10 +614,10 @@ columns_define_alloc(struct mkcatalogparams *p)
columns_alloc_clumpsgeoradec(p);
break;
- case UI_KEY_CLUMPSGEODEC:
- name = "CLUMPS_DEC";
- unit = "degrees";
- ocomment = "Declination of all clumps geometric center.";
+ case UI_KEY_CLUMPSGEOW2:
+ name = gal_checkset_malloc_cat("CLUMPS_GEO", p->ctype[1]);
+ unit = p->input->wcs->cunit[1];
+ ocomment = "Geometric center of all clumps in 2nd WCS axis.";
ccomment = NULL;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_INVALID;
@@ -650,8 +728,7 @@ columns_define_alloc(struct mkcatalogparams *p)
disp_fmt = GAL_TABLE_DISPLAY_FMT_GENERAL;
disp_width = 8;
disp_precision = 3;
- p->upperlimit = 1;
- /* Upper limit measurement doesn't need per-pixel calculations. */
+ p->upperlimit = 1; /* Doesn't need per-pixel calculations. */
break;
case UI_KEY_UPPERLIMITMAG:
@@ -665,8 +742,7 @@ columns_define_alloc(struct mkcatalogparams *p)
disp_width = 8;
disp_precision = 3;
p->upperlimit = 1;
- p->hasmag = 1;
- /* Upper limit magnitude doesn't need per-pixel calculations. */
+ p->hasmag = 1; /* Doesn't need per-pixel calculations. */
break;
case UI_KEY_RIVERAVE:
@@ -896,6 +972,7 @@ columns_define_alloc(struct mkcatalogparams *p)
"column code", __func__, PACKAGE_BUGREPORT, colcode->v);
}
+
/* If this is an objects column, add it to the list of columns. We
will be using the `status' element to keep the MakeCatalog code
for the columns. */
@@ -1246,40 +1323,44 @@ columns_fill(struct mkcatalog_passparams *pp)
oi[OCOL_C_NUMALL] );
break;
- case UI_KEY_RA:
- case UI_KEY_DEC:
+ case UI_KEY_W1:
+ case UI_KEY_W2:
p->rd_vo[0][oind] = POS_V_G(oi, OCOL_SUMWHT, OCOL_NUMALL,
OCOL_VX, OCOL_GX);
p->rd_vo[1][oind] = POS_V_G(oi, OCOL_SUMWHT, OCOL_NUMALL,
OCOL_VY, OCOL_GY);
break;
- case UI_KEY_GEORA:
- case UI_KEY_GEODEC:
+ case UI_KEY_GEOW1:
+ case UI_KEY_GEOW2:
p->rd_go[0][oind] = MKC_RATIO( oi[OCOL_GX], oi[OCOL_NUMALL] );
p->rd_go[1][oind] = MKC_RATIO( oi[OCOL_GY], oi[OCOL_NUMALL] );
break;
- case UI_KEY_CLUMPSRA:
- case UI_KEY_CLUMPSDEC:
+ case UI_KEY_CLUMPSW1:
+ case UI_KEY_CLUMPSW2:
p->rd_vcc[0][oind] = POS_V_G(oi, OCOL_C_SUMWHT, OCOL_C_NUMALL,
OCOL_C_VX, OCOL_C_GX);
p->rd_vcc[1][oind] = POS_V_G(oi, OCOL_C_SUMWHT, OCOL_C_NUMALL,
OCOL_C_VY, OCOL_C_GY);
break;
- case UI_KEY_CLUMPSGEORA:
- case UI_KEY_CLUMPSGEODEC:
+ case UI_KEY_CLUMPSGEOW1:
+ case UI_KEY_CLUMPSGEOW2:
p->rd_gcc[0][oind] = MKC_RATIO( oi[OCOL_C_GX], oi[OCOL_C_NUMALL] );
p->rd_gcc[1][oind] = MKC_RATIO( oi[OCOL_C_GY], oi[OCOL_C_NUMALL] );
break;
case UI_KEY_BRIGHTNESS:
- ((float *)colarr)[oind] = oi[ OCOL_NUM ]>0.0f ? oi[ OCOL_SUM ] : NAN;
+ ((float *)colarr)[oind] = ( oi[ OCOL_NUM ]>0.0f
+ ? oi[ OCOL_SUM ]
+ : NAN );
break;
case UI_KEY_CLUMPSBRIGHTNESS:
- ((float *)colarr)[oind] = oi[ OCOL_C_NUM ]>0.0f ?oi[ OCOL_C_SUM
]:NAN;
+ ((float *)colarr)[oind] = ( oi[ OCOL_C_NUM ]>0.0f
+ ? oi[ OCOL_C_SUM ]
+ : NAN );
break;
case UI_KEY_MAGNITUDE:
@@ -1407,16 +1488,16 @@ columns_fill(struct mkcatalog_passparams *pp)
ci[CCOL_NUMALL] );
break;
- case UI_KEY_RA:
- case UI_KEY_DEC:
+ case UI_KEY_W1:
+ case UI_KEY_W2:
p->rd_vc[0][cind] = POS_V_G(ci, CCOL_SUMWHT, CCOL_NUMALL,
CCOL_VX, CCOL_GX);
p->rd_vc[1][cind] = POS_V_G(ci, CCOL_SUMWHT, CCOL_NUMALL,
CCOL_VY, CCOL_GY);
break;
- case UI_KEY_GEORA:
- case UI_KEY_GEODEC:
+ case UI_KEY_GEOW1:
+ case UI_KEY_GEOW2:
p->rd_gc[0][cind] = MKC_RATIO( ci[CCOL_GX], ci[CCOL_NUMALL] );
p->rd_gc[1][cind] = MKC_RATIO( ci[CCOL_GY], ci[CCOL_NUMALL] );
break;
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index 0dd53e9..39db6f6 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -136,7 +136,7 @@ struct mkcatalogparams
{
/* From command-line */
struct gal_options_common_params cp; /* Common parameters. */
- gal_list_i32_t *columnids; /* The desired column codes. */
+ gal_list_i32_t *columnids; /* The desired column codes. */
char *inputname; /* Input filename. */
char *objectsfile; /* File name of objects file. */
char *objectshdu; /* HDU of objects image. */
@@ -195,6 +195,8 @@ struct mkcatalogparams
double **rd_vcc; /* All clumps RA-Dec flx. wht. X, Y. */
double **rd_gcc; /* All clumps RA-Dec geometric X, Y. */
+ char **ctype; /* Type of WCS axis. */
+
uint8_t hasblank; /* Dataset has blank values. */
uint8_t hasmag; /* Catalog has magnitude columns. */
uint8_t upperlimit; /* Calculate upper limit magnitude. */
diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index 5a90129..e41278d 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -604,14 +604,14 @@ mkcatalog_wcs_conversion(struct mkcatalogparams *p)
probably columns that don't need any correction. */
switch(column->status)
{
- case UI_KEY_RA: c=p->rd_vo[0]; break;
- case UI_KEY_DEC: c=p->rd_vo[1]; break;
- case UI_KEY_GEORA: c=p->rd_go[0]; break;
- case UI_KEY_GEODEC: c=p->rd_go[1]; break;
- case UI_KEY_CLUMPSRA: c=p->rd_vcc[0]; break;
- case UI_KEY_CLUMPSDEC: c=p->rd_vcc[1]; break;
- case UI_KEY_CLUMPSGEORA: c=p->rd_gcc[0]; break;
- case UI_KEY_CLUMPSGEODEC: c=p->rd_gcc[1]; break;
+ case UI_KEY_W1: c=p->rd_vo[0]; break;
+ case UI_KEY_W2: c=p->rd_vo[1]; break;
+ case UI_KEY_GEOW1: c=p->rd_go[0]; break;
+ case UI_KEY_GEOW2: c=p->rd_go[1]; break;
+ case UI_KEY_CLUMPSW1: c=p->rd_vcc[0]; break;
+ case UI_KEY_CLUMPSW2: c=p->rd_vcc[1]; break;
+ case UI_KEY_CLUMPSGEOW1: c=p->rd_gcc[0]; break;
+ case UI_KEY_CLUMPSGEOW2: c=p->rd_gcc[1]; break;
}
/* Copy the elements. */
@@ -630,10 +630,10 @@ mkcatalog_wcs_conversion(struct mkcatalogparams *p)
probably columns that don't need any correction. */
switch(column->status)
{
- case UI_KEY_RA: c=p->rd_vc[0]; break;
- case UI_KEY_DEC: c=p->rd_vc[1]; break;
- case UI_KEY_GEORA: c=p->rd_gc[0]; break;
- case UI_KEY_GEODEC: c=p->rd_gc[1]; break;
+ case UI_KEY_W1: c=p->rd_vc[0]; break;
+ case UI_KEY_W2: c=p->rd_vc[1]; break;
+ case UI_KEY_GEOW1: c=p->rd_gc[0]; break;
+ case UI_KEY_GEOW2: c=p->rd_gc[1]; break;
}
/* Copy the elements. */
@@ -718,8 +718,11 @@ mkcatalog_outputs_same_start(struct mkcatalogparams *p,
int o0c1,
if(p->input->wcs)
{
pixarea=gal_wcs_pixel_area_arcsec2(p->input->wcs);
- asprintf(&str, "Pixel area (arcsec^2): %g", pixarea);
- gal_list_str_add(&comments, str, 0);
+ if( isnan(pixarea)==0 )
+ {
+ asprintf(&str, "Pixel area (arcsec^2): %g", pixarea);
+ gal_list_str_add(&comments, str, 0);
+ }
}
if(p->hasmag)
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index 4e89812..9b86483 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -29,7 +29,6 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <string.h>
#include <inttypes.h>
-#include <gnuastro/wcs.h>
#include <gnuastro/fits.h>
#include <gnuastro/blank.h>
#include <gnuastro/threads.h>
@@ -325,6 +324,42 @@ ui_check_options_and_arguments(struct mkcatalogparams *p)
/**************************************************************/
/*************** Preparations *******************/
/**************************************************************/
+/* If a WCS structure is present, then read its basic information to use in
+ the table meta-data. */
+static void
+ui_wcs_info(struct mkcatalogparams *p)
+{
+ char *c;
+ size_t i;
+
+ /* Read the basic WCS information. */
+ if(p->input->wcs)
+ {
+ /* Allocate space for the array of strings. */
+ errno=0;
+ p->ctype=malloc(p->input->ndim * sizeof *p->ctype);
+ if(p->ctype==NULL)
+ error(EXIT_FAILURE, 0, "%s: %zu bytes for `p->ctype'", __func__,
+ p->input->ndim * sizeof *p->ctype);
+
+ /* Fill in the values. */
+ for(i=0;i<p->input->ndim;++i)
+ {
+ /* CTYPE might contain `-' characters, we just want the first
+ non-dash characters. The loop will either stop either at the end
+ or where there is a dash. So we can just replace it with an
+ end-of-string character. */
+ gal_checkset_allocate_copy(p->input->wcs->ctype[i], &p->ctype[i]);
+ for(c=p->ctype[i]; *c!='\0' && *c!='-'; ++c) c=c;
+ *c='\0';
+ }
+ }
+}
+
+
+
+
+
static void
ui_preparations_read_inputs(struct mkcatalogparams *p)
{
@@ -343,9 +378,13 @@ ui_preparations_read_inputs(struct mkcatalogparams *p)
GAL_TYPE_FLOAT32, p->cp.minmapsize, 0,0);
+ /* Read basic WCS information for final table meta-data. */
+ ui_wcs_info(p);
+
+
/* Currently MakeCatalog is only implemented for 2D images. */
if(p->input->ndim!=2)
- error(EXIT_FAILURE, 0, "%s (%s) has %zu dimensions, MakeCatalog "
+ error(EXIT_FAILURE, 0, "%s (hdu %s) has %zu dimensions, MakeCatalog "
"currently only supports 2D inputs", p->inputname, p->cp.hdu,
p->input->ndim);
@@ -994,6 +1033,14 @@ ui_free_report(struct mkcatalogparams *p, struct timeval
*t1)
if(p->rd_gcc) free(p->rd_gcc);
}
+ /* Free the types of the WCS coordinates (for catalog meta-data). */
+ if(p->ctype)
+ {
+ for(d=0;d<p->input->ndim;++d)
+ free(p->ctype[d]);
+ free(p->ctype);
+ }
+
/* If a random number generator was allocated, free it. */
if(p->rng) gsl_rng_free(p->rng);
diff --git a/bin/mkcatalog/ui.h b/bin/mkcatalog/ui.h
index c5ab69a..5173056 100644
--- a/bin/mkcatalog/ui.h
+++ b/bin/mkcatalog/ui.h
@@ -29,7 +29,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Available letters for short options:
- f g k l u v w x y
+ f g k l u v w x y z
H J L W X Y
*/
enum option_keys_enum
@@ -39,7 +39,6 @@ enum option_keys_enum
UI_KEY_CLUMPSFILE = 'C',
UI_KEY_SKYFILE = 's',
UI_KEY_STDFILE = 't',
- UI_KEY_ZEROPOINT = 'z',
UI_KEY_SKYSUBTRACTED = 'E',
UI_KEY_THRESHOLD = 'R',
UI_KEY_ENVSEED = 'e',
@@ -68,6 +67,7 @@ enum option_keys_enum
UI_KEY_CLUMPSHDU,
UI_KEY_SKYHDU,
UI_KEY_STDHDU,
+ UI_KEY_ZEROPOINT,
UI_KEY_SFMAGNSIGMA,
UI_KEY_SFMAGAREA,
UI_KEY_UPMASKFILE,
@@ -86,12 +86,14 @@ enum option_keys_enum
UI_KEY_CLUMPSY,
UI_KEY_CLUMPSGEOX,
UI_KEY_CLUMPSGEOY,
- UI_KEY_GEORA,
- UI_KEY_GEODEC,
- UI_KEY_CLUMPSRA,
- UI_KEY_CLUMPSDEC,
- UI_KEY_CLUMPSGEORA,
- UI_KEY_CLUMPSGEODEC,
+ UI_KEY_W1,
+ UI_KEY_W2,
+ UI_KEY_GEOW1,
+ UI_KEY_GEOW2,
+ UI_KEY_CLUMPSW1,
+ UI_KEY_CLUMPSW2,
+ UI_KEY_CLUMPSGEOW1,
+ UI_KEY_CLUMPSGEOW2,
UI_KEY_CLUMPSBRIGHTNESS,
UI_KEY_NORIVERBRIGHTNESS,
UI_KEY_CLUMPSMAGNITUDE,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index af0f025..2868fec 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})
@@ -2088,9 +2089,9 @@ Warp started on Mon Apr 6 16:51:59 953
Using 8 CPU threads.
Input: cat_convolved.fits (hdu: 1)
matrix:
- 0.2000 0.0000 0.4000
- 0.0000 0.2000 0.4000
- 0.0000 0.0000 1.0000
+ 0.2000 0.0000 0.4000
+ 0.0000 0.2000 0.4000
+ 0.0000 0.0000 1.0000
$ ls
0_cat.fits cat_convolved_scaled.fits cat.txt
@@ -8691,6 +8692,49 @@ standard deviation of the respective pixels in all
operands in the stack.
Similar to @command{min}, but the pixels of the output will contain
the median of the respective pixels in all operands in the stack.
+@item filter-mean
+Apply mean filtering (or @url{https://en.wikipedia.org/wiki/Moving_average,
+moving average}) on the input dataset. During mean filtering, each pixel
+(data element) is replaced by the mean value of all its surrounding pixels
+(excluding blank values). The number of surrounding pixels in each
+dimension (to calculate the mean) is determined through the earlier
+operands that have been pushed onto the stack prior to the input
+dataset. The number of necessary operands is determined by the
+dimensionality of the input dataset (first popped operand). The order of
+the dimensions on the command-line is the order in FITS format. Here is one
+example:
+
+@example
+$ astarithmetic 5 4 image.fits filter-mean
+@end example
+
+@noindent
+In this example, each pixel is replaced by the mean of a 5 by 4 box around
+it. The box is 5 pixels along the first FITS dimension (horizontal when
+viewed in ds9) and 4 pixels along the second FITS dimension (vertical).
+
+Each pixel will be placed in the center of the box that the mean is
+calculated on. If the given width along a dimension is even, then the
+center is assumed to be between the pixels (not in the center of a
+pixel). When the pixel is close to the center, the pixels of the box that
+fall outside the image are ignored. Therefore, on the edge, less points
+will be used in calculating the mean.
+
+The final effect of mean filtering is to smooth the input image, it is
+essentially a convolution with a kernel that has identical values for all
+its pixels (is flat), see @ref{Convolution process}.
+
+@item filter-median
+Apply @url{https://en.wikipedia.org/wiki/Median_filter, median filtering}
+on the input dataset. This is very similar to @command{filter-mean}, except
+that instead of the mean value of the box pixels, the median value is used
+to replace a pixel value. For more on how to use this operator, please see
+@command{filter-mean}.
+
+The median is less susceptible to outliers compared to the mean. As a
+result, after median filtering, the pixel values will be more discontinuous
+than mean filtering.
+
@item lt
Less than: If the second popped (or left operand in infix notation, see
@ref{Reverse polish notation}) value is smaller than the first popped
@@ -14213,7 +14257,7 @@ only one of them. The latter cases are explicitly
marked with [Objects] or
@item --i
@itemx --ids
-This is a unique option it can add multiple columns to the final
+This is a unique option which can add multiple columns to the final
catalog(s). Calling this option will put the object IDs (@option{--objid})
in the objects catalog and host-object-ID (@option{--hostobjid}) and
ID-in-host-object (@option{--idinhostobj}) into the clumps catalog. Hence
@@ -14289,35 +14333,61 @@ the second FITS axis. See @option{--geox}.
@item -r
@itemx --ra
-Flux weighted right ascension of all objects or clumps, see @option{--x}.
+Flux weighted right ascension of all objects or clumps, see
+@option{--x}. This is just an alias for one of the lower-level
+@option{--w1} or @option{--w2} options. Using the FITS WCS keywords
+(@code{CTYPE}), MakeCatalog will determine which axis corresponds to the
+right ascension. If no @code{CTYPE} keywords start with @code{RA}, an error
+will be printed when requesting this column and MakeCatalog will abort.
@item -d
@itemx --dec
-Flux weighted declination of all objects or clumps, see @option{--x}.
+Flux weighted declination of all objects or clumps, see @option{--x}. This
+is just an alias for one of the lower-level @option{--w1} or @option{--w2}
+options. Using the FITS WCS keywords (@code{CTYPE}), MakeCatalog will
+determine which axis corresponds to the declination. If no @code{CTYPE}
+keywords start with @code{DEC}, an error will be printed when requesting
+this column and MakeCatalog will abort.
+
+@item --w1
+Flux weighted first WCS axis of all objects or clumps, see
+@option{--x}. The first WCS axis is commonly used as right ascension in
+images.
-@item --geora
-Geometric center right ascension of all objects or clumps, see
-@option{--geox}.
+@item --w2
+Flux weighted second WCS axis of all objects or clumps, see
+@option{--x}. The second WCS axis is commonly used as declination in
+images.
-@item --geodec
-Geometric center declination of all objects or clumps, see
-@option{--geox}.
+@item --geow1
+Geometric center in first WCS axis of all objects or clumps, see
+@option{--geox}. The first WCS axis is commonly used as right ascension in
+images.
-@item --clumpsra
-[Objects] Flux weighted right ascension of all clumps in this object,
-see @option{--x}.
+@item --geow2
+Geometric center in second WCS axis of all objects or clumps, see
+@option{--geox}. The second WCS axis is commonly used as declination in
+images.
-@item --clumpsdec
+@item --clumpsw1
+[Objects] Flux weighted center in first WCS axis of all clumps in this
+object, see @option{--x}. The first WCS axis is commonly used as right
+ascension in images.
+
+@item --clumpsw2
[Objects] Flux weighted declination of all clumps in this object, see
-@option{--x}.
+@option{--x}. The second WCS axis is commonly used as declination in
+images.
-@item --clumpsgeora
-[Objects] Geometric center right ascension of all clumps in this
-object, see @option{--geox}.
+@item --clumpsgeow1
+[Objects] Geometric center right ascension of all clumps in this object,
+see @option{--geox}. The first WCS axis is commonly used as right ascension
+in images.
-@item --clumpsgeodec
-[Objects] Geometric center declination of all clumps in this object,
-see @option{--geox}.
+@item --clumpsgeow2
+[Objects] Geometric center declination of all clumps in this object, see
+@option{--geox}. The second WCS axis is commonly used as declination in
+images.
@item -b
@itemx --brightness
@@ -16182,8 +16252,8 @@ One line examples:
## Add noise with a standard deviation of 100 to image:
$ astmknoise --sigma=100 image.fits
-## Add noise to input image assuming a background magnitude (with zeropoint
-## magnitude of 0) and a certain instrumental noise:
+## Add noise to input image assuming a background magnitude (with
+## zeropoint magnitude of 0) and a certain instrumental noise:
$ astmknoise --background=-10 -z0 --instrumental=20 mockimage.fits
@end example
@@ -16305,36 +16375,44 @@ interested readers can study those books.
@node Distance on a 2D curved space, Extending distance concepts to 3D,
CosmicCalculator, CosmicCalculator
@subsection Distance on a 2D curved space
-The observations to date (for example the Plank 2013 results), have
-not measured the presence of a significant curvature in the
-universe. However to be generic (and allow its measurement if it does
-in fact exist), it is very important to create a framework that allows
-curvature. As 3D beings, it is impossible for us to mentally create
-(visualize) a picture of the curvature of a 3D volume in a 4D
-space. Hence, here we will assume a 2D surface and discuss distances
-on that 2D surface when it is flat, or when the 2D surface is curved
-(in a 3D space). Once the concepts have been created/visualized here,
-in @ref{Extending distance concepts to 3D}, we will extend them to the
-real 3D universe we live in and hope to study.
-
-To be more understandable (actively discuss from an observer's point
-of view) let's assume we have an imaginary 2D friend living on the 2D
-space (which @emph{might} be curved in 3D). So here we will be working
-with it in its efforts to analyze distances on its 2D universe. The
-start of the analysis might seem too mundane, but since it is
-impossible to imagine a 3D curved space, it is important to review all
-the very basic concepts thoroughly for an easy transition to a
-universe we cannot visualize any more (a curved 3D space in 4D).
+The observations to date (for example the Planck 2015 results), have not
+measured@footnote{The observations are interpeted under the assumption of
+uniform curvature. For a relativistic alternative to dark energy (and maybe
+also some part of dark matter), non-uniform curvature may be even be more
+critical, but that is beyond the scope of this brief explanation.} the
+presence of significant curvature in the universe. However to be generic
+(and allow its measurement if it does in fact exist), it is very important
+to create a framework that allows non-zero uniform curvature. However, this
+section is not intended to be a fully thorough and mathematically complete
+derivation of these concepts. There are many references available for such
+reviews that go deep into the abstract mathematical proofs. The emphasis
+here is on visualization of the concepts for a beginner.
+
+As 3D beings, it is difficult for us to mentally create (visualize) a
+picture of the curvature of a 3D volume. Hence, here we will assume a 2D
+surface/space and discuss distances on that 2D surface when it is flat and
+when it is curved. Once the concepts have been created/visualized here, we
+will extend them, in @ref{Extending distance concepts to 3D}, to a real 3D
+spatial @emph{slice} of the Universe we live in and hope to study.
+
+To be more understandable (actively discuss from an observer's point of
+view) let's assume there's an imaginary 2D creature living on the 2D space
+(which @emph{might} be curved in 3D). Here, we will be working with this
+creature in its efforts to analyze distances in its 2D universe. The start
+of the analysis might seem too mundane, but since it is difficult to
+imagine a 3D curved space, it is important to review all the very basic
+concepts thoroughly for an easy transition to a universe that is more
+difficult to visualize (a curved 3D space embedded in 4D).
To start, let's assume a static (not expanding or shrinking), flat 2D
-surface similar to @ref{flatplane} and that our 2D friend is observing its
-universe from point @mymath{A}. One of the most basic ways to parametrize
-this space is through the Cartesian coordinates (@mymath{x},
+surface similar to @ref{flatplane} and that the 2D creature is observing
+its universe from point @mymath{A}. One of the most basic ways to
+parametrize this space is through the Cartesian coordinates (@mymath{x},
@mymath{y}). In @ref{flatplane}, the basic axes of these two coordinates
are plotted. An infinitesimal change in the direction of each axis is
written as @mymath{dx} and @mymath{dy}. For each point, the infinitesimal
changes are parallel with the respective axes and are not shown for
-clarity. Another very useful way of parameterizing this space is through
+clarity. Another very useful way of parametrizing this space is through
polar coordinates. For each point, we define a radius (@mymath{r}) and
angle (@mymath{\phi}) from a fixed (but arbitrary) reference axis. In
@ref{flatplane} the infinitesimal changes for each polar coordinate are
@@ -16348,69 +16426,72 @@ the same radius.
plane.}
@end float
-Assuming a certain position, which can be parameterized as @mymath{(x,y)},
-or @mymath{(r,\phi)}, a general infinitesimal change change in its position
-will place it in the coordinates @mymath{(x+dx,y+dy)} and
-@mymath{(r+dr,\phi+d\phi)}. The distance (on the flat 2D surface) that is
-covered by this infinitesimal change in the static universe (@mymath{ds_s},
-the subscript signifies the static nature of this universe) can be written
-as:
+Assuming an object is placed at a certain position, which can be
+parameterized as @mymath{(x,y)}, or @mymath{(r,\phi)}, a general
+infinitesimal change in its position will place it in the coordinates
+@mymath{(x+dx,y+dy)} and @mymath{(r+dr,\phi+d\phi)}. The distance (on the
+flat 2D surface) that is covered by this infinitesimal change in the static
+universe (@mymath{ds_s}, the subscript signifies the static nature of this
+universe) can be written as:
@dispmath{ds_s=dx^2+dy^2=dr^2+r^2d\phi^2}
-The main question is this: how can our 2D friend incorporate the (possible)
-curvature in its universe when it is calculating distances? The universe it
-lives in might equally be a locally flat but globally curved surface like
-@ref{sphericalplane}. The answer to this question but for a 3D being (us)
-is the whole purpose to this discussion. So here we want to give our 2D
-friend (and later, ourselves) the tools to measure distances if the space
+The main question is this: how can the 2D creature incorporate the
+(possible) curvature in its universe when it's calculating distances? The
+universe that it lives in might equally be a curved surface like
+@ref{sphereandplane}. The answer to this question but for a 3D being (us)
+is the whole purpose to this discussion. Here, we want to give the 2D
+creature (and later, ourselves) the tools to measure distances if the space
(that hosts the objects) is curved.
-@ref{sphericalplane} assumes a spherical shell with radius @mymath{R}
-as the curved 2D plane for simplicity. The spherical shell is tangent
-to the 2D plane and only touches it at @mymath{A}. The result will be
-generalized afterwards. The first step in measuring the distance in a
-curved space is to imagine a third dimension along the @mymath{z} axis
-as shown in @ref{sphericalplane}. For simplicity, the @mymath{z} axis
-is assumed to pass through the center of the spherical shell. Our
-imaginary 2D friend cannot visualize the third dimension or a curved
-2D surface within it, so the remainder of this discussion is purely
-abstract for it (similar to us being unable to visualize a 3D curved
-space in 4D). But since we are 3D creatures, we have the advantage of
-visualizing the following steps. Fortunately our 2D friend knows our
-mathematics, so it can follow along with us.
-
-With the third axis added, a generic infinitesimal change over
-@emph{the full} 3D space corresponds to the distance:
-@dispmath{ds_s^2=dx^2+dy^2+dz^2=dr^2+r^2d\phi^2+dz^2.}It is very
-important to recognize that this change of distance is for @emph{any}
-point in the 3D space, not just those changes that occur on the 2D
-spherical shell of @ref{sphericalplane}. Recall that our 2D friend can
-only do measurements in the 2D spherical shell, not the full 3D
-space. So we have to constrain this general change to any change on
-the 2D spherical shell. To do that, let's look at the arbitrary point
-@mymath{P} on the 2D spherical shell. Its image (@mymath{P'}) on the
-flat plain is also displayed. From the dark triangle, we see that
-
-@float Figure,sphericalplane
-@center@image{gnuastro-figures/sphericalplane, 10cm, , }
-
-@caption{2D spherical plane (centered on @mymath{O}) and flat plane
-(gray) tangent to it at point @mymath{A}.}
+@ref{sphereandplane} assumes a spherical shell with radius @mymath{R} as
+the curved 2D plane for simplicity. The 2D plane is tangent to the
+spherical shell and only touches it at @mymath{A}. This idea will be
+generalized later. The first step in measuring the distance in a curved
+space is to imagine a third dimension along the @mymath{z} axis as shown in
+@ref{sphereandplane}. For simplicity, the @mymath{z} axis is assumed to
+pass through the center of the spherical shell. Our imaginary 2D creature
+cannot visualize the third dimension or a curved 2D surface within it, so
+the remainder of this discussion is purely abstract for it (similar to us
+having difficulty in visualizing a 3D curved space in 4D). But since we are
+3D creatures, we have the advantage of visualizing the following
+steps. Fortunately the 2D creature is already familiar with our
+mathematical constructs, so it can follow our reasoning.
+
+With the third axis added, a generic infinitesimal change over @emph{the
+full} 3D space corresponds to the distance:
+
+@dispmath{ds_s^2=dx^2+dy^2+dz^2=dr^2+r^2d\phi^2+dz^2.}
+
+@float Figure,sphereandplane
+@center@image{gnuastro-figures/sphereandplane, 10cm, , }
+
+@caption{2D spherical shell (centered on @mymath{O}) and flat plane (light
+gray) tangent to it at point @mymath{A}.}
@end float
+It is very important to recognize that this change of distance is for
+@emph{any} point in the 3D space, not just those changes that occur on the
+2D spherical shell of @ref{sphereandplane}. Recall that our 2D friend can
+only do measurements on the 2D surfaces, not the full 3D space. So we have
+to constrain this general change to any change on the 2D spherical
+shell. To do that, let's look at the arbitrary point @mymath{P} on the 2D
+spherical shell. Its image (@mymath{P'}) on the flat plain is also
+displayed. From the dark gray triangle, we see that
+
@dispmath{\sin\theta={r\over R},\quad\cos\theta={R-z\over R}.}These
-relations allow our 2D friend to find the value of @mymath{z} (an
-abstract dimension for it) as a function of r (distance on a flat 2D
-plane, which it can visualize) and thus eliminate @mymath{z}. From
-@mymath{\sin^2\theta+\cos^2\theta=1}, we get @mymath{z^2-2Rz+r^2=0}
-and solving for @mymath{z}, we find:
-@dispmath{z=R\left(1\pm\sqrt{1-{r^2\over R^2}}\right).}The
-@mymath{\pm} can be understood from @ref{sphericalplane}: For each
-@mymath{r}, there are two points on the sphere, one in the upper
-hemisphere and one in the lower hemisphere. An infinitesimal change in
-@mymath{r}, will create the following infinitesimal change in
-@mymath{z}:
+relations allow the 2D creature to find the value of @mymath{z} (an
+abstract dimension for it) as a function of r (distance on a flat 2D plane,
+which it can visualize) and thus eliminate @mymath{z}. From
+@mymath{\sin^2\theta+\cos^2\theta=1}, we get @mymath{z^2-2Rz+r^2=0} and
+solving for @mymath{z}, we find:
+
+@dispmath{z=R\left(1\pm\sqrt{1-{r^2\over R^2}}\right).}
+
+The @mymath{\pm} can be understood from @ref{sphereandplane}: For each
+@mymath{r}, there are two points on the sphere, one in the upper hemisphere
+and one in the lower hemisphere. An infinitesimal change in @mymath{r},
+will create the following infinitesimal change in @mymath{z}:
@dispmath{dz={\mp r\over R}\left(1\over
\sqrt{1-{r^2/R^2}}\right)dr.}Using the positive signed equation
@@ -16418,45 +16499,48 @@ instead of @mymath{dz} in the @mymath{ds_s^2}
equation above, we get:
@dispmath{ds_s^2={dr^2\over 1-r^2/R^2}+r^2d\phi^2.}
-The derivation above was done for a spherical shell of radius
-@mymath{R} as a curved 2D surface. To generalize it to any surface, we
-can define @mymath{K=1/R^2} as the curvature parameter. Then the
-general infinitesimal change in a static universe can be written as:
-@dispmath{ds_s^2={dr^2\over 1-Kr^2}+r^2d\phi^2.}Therefore, we see that
-a positive @mymath{K} represents a real @mymath{R} which signifies a
-closed 2D spherical shell like @ref{sphericalplane}. When
-@mymath{K=0}, we have a flat plane (@ref{flatplane}) and a negative
-@mymath{K} will correspond to an imaginary @mymath{R}. The latter two
-cases are open universes (where @mymath{r} can extend to infinity).
-However, when @mymath{K>0}, we have a closed universe, where
-@mymath{r} cannot become larger than @mymath{R} as in
-@ref{sphericalplane}.
+The derivation above was done for a spherical shell of radius @mymath{R} as
+a curved 2D surface. To generalize it to any surface, we can define
+@mymath{K=1/R^2} as the curvature parameter. Then the general infinitesimal
+change in a static universe can be written as:
+
+@dispmath{ds_s^2={dr^2\over 1-Kr^2}+r^2d\phi^2.}
+
+Therefore, when @mymath{K>0} (and curvature is the same everywhere), we
+have a finite universe, where @mymath{r} cannot become larger than
+@mymath{R} as in @ref{sphereandplane}. When @mymath{K=0}, we have a flat
+plane (@ref{flatplane}) and a negative @mymath{K} will correspond to an
+imaginary @mymath{R}. The latter two cases may be infinite in area (which
+is not a simple concept, but mathematically can be modelled with @mymath{r}
+extending infinitely), or finite-area (like a cylinder is flat everywhere
+with @mymath{ds_s^2={dx^2 + dy^2}}, but finite in one direction in size).
@cindex Proper distance
-A very important issue that can be discussed now (while we are still
-in 2D and can actually visualize things) is that
-@mymath{\overrightarrow{r}} is tangent to the curved space at the
-observer's position. In other words, it is on the gray flat surface of
-@ref{sphericalplane}, even when the universe if curved:
-@mymath{\overrightarrow{r}=P'-A}. Therefore for the point @mymath{P}
-on a curved space, the raw coordinate @mymath{r} is the distance to
-@mymath{P'}, not @mymath{P}. The distance to the point @mymath{P} (at
-a specific coordinate @mymath{r} on the flat plane) on the curved
-surface (thick line in @ref{sphericalplane}) is called the
-@emph{proper distance} and is displayed with @mymath{l}. For the
-specific example of @ref{sphericalplane}, the proper distance can be
-calculated with: @mymath{l=R\theta} (@mymath{\theta} is in
-radians). using the @mymath{\sin\theta} relation found above, we can
-find @mymath{l} as a function of @mymath{r}:
+A very important issue that can be discussed now (while we are still in 2D
+and can actually visualize things) is that @mymath{\overrightarrow{r}} is
+tangent to the curved space at the observer's position. In other words, it
+is on the gray flat surface of @ref{sphereandplane}, even when the universe
+if curved: @mymath{\overrightarrow{r}=P'-A}. Therefore for the point
+@mymath{P} on a curved space, the raw coordinate @mymath{r} is the distance
+to @mymath{P'}, not @mymath{P}. The distance to the point @mymath{P} (at a
+specific coordinate @mymath{r} on the flat plane) over the curved surface
+(thick line in @ref{sphereandplane}) is called the @emph{proper distance}
+and is displayed with @mymath{l}. For the specific example of
+@ref{sphereandplane}, the proper distance can be calculated with:
+@mymath{l=R\theta} (@mymath{\theta} is in radians). using the
+@mymath{\sin\theta} relation found above, we can find @mymath{l} as a
+function of @mymath{r}:
@dispmath{\theta=\sin^{-1}\left({r\over R}\right)\quad\rightarrow\quad
-l(r)=R\sin^{-1}\left({r\over R}\right)}@mymath{R} is just an arbitrary
-constant and can be directly found from @mymath{K}, so for cleaner
-equations, it is common practice to set @mymath{R=1}, which gives:
-@mymath{l(r)=\sin^{-1}r}. Also note that if @mymath{R=1}, then
-@mymath{l=\theta}. Generally, depending on the the curvature, in a
-@emph{static} universe the proper distance can be written as a
-function of the coordinate @mymath{r} as (from now on we are assuming
+l(r)=R\sin^{-1}\left({r\over R}\right)}
+
+
+@mymath{R} is just an arbitrary constant and can be directly found from
+@mymath{K}, so for cleaner equations, it is common practice to set
+@mymath{R=1}, which gives: @mymath{l(r)=\sin^{-1}r}. Also note that when
+@mymath{R=1}, then @mymath{l=\theta}. Generally, depending on the the
+curvature, in a @emph{static} universe the proper distance can be written
+as a function of the coordinate @mymath{r} as (from now on we are assuming
@mymath{R=1}):
@dispmath{l(r)=\sin^{-1}(r)\quad(K>0),\quad\quad
@@ -16467,71 +16551,91 @@ more simpler and abstract form of
@dispmath{ds_s^2=dl^2+r^2d\phi^2.}
@cindex Comoving distance
-Until now, we had assumed a static universe (not changing with
-time). But our observations so far appear to indicate that the
-universe is expanding (isn't static). Since there is no reason to
-expect the observed expansion is unique to our particular position of
-the universe, we expect the universe to be expanding at all points
-with the same rate at the same time. Therefore, to add a time
-dependence to our distance measurements, we can simply add a
-multiplicative scaling factor, which is a function of time:
+Until now, we had assumed a static universe (not changing with time). But
+our observations so far appear to indicate that the universe is expanding
+(it isn't static). Since there is no reason to expect the observed
+expansion is unique to our particular position of the universe, we expect
+the universe to be expanding at all points with the same rate at the same
+time. Therefore, to add a time dependence to our distance measurements, we
+can include a multiplicative scaling factor, which is a function of time:
@mymath{a(t)}. The functional form of @mymath{a(t)} comes from the
-cosmology and the physics we assume for it: general relativity.
-
-With this scaling factor, the proper distance will also depend on
-time. As the universe expands (moves), the distance will also move to
-larger values. We thus define a distance measure, or coordinate, that
-is independent of time and thus doesn't `move' which we call the
-@emph{comoving distance} and display with @mymath{\chi} such that:
-@mymath{l(r,t)=\chi(r)a(t)}. We thus shift the @mymath{r} dependence
-of the proper distance we derived above for a static universe to the
-comoving distance:
+cosmology, the physics we assume for it: general relativity, and the choice
+of whether the universe is uniform (`homogeneous') in density and curvature
+or inhomogeneous. In this section, the functional form of @mymath{a(t)} is
+irrelevant, so we can aviod these issues.
+
+With this scaling factor, the proper distance will also depend on time. As
+the universe expands, the distance between two given points will shift to
+larger values. We thus define a distance measure, or coordinate, that is
+independent of time and thus doesn’t `move'. We call it the @emph{comoving
+distance} and display with @mymath{\chi} such that:
+@mymath{l(r,t)=\chi(r)a(t)}. We have therefore, shifted the @mymath{r}
+dependence of the proper distance we derived above for a static universe to
+the comoving distance:
@dispmath{\chi(r)=\sin^{-1}(r)\quad(K>0),\quad\quad
\chi(r)=r\quad(K=0),\quad\quad \chi(r)=\sinh^{-1}(r)\quad(K<0).}
-Therefore @mymath{\chi(r)} is the proper distance of an object at a
-specific reference time: @mymath{t=t_r} (the @mymath{r} subscript
-signifies ``reference'') when @mymath{a(t_r)=1}. At any arbitrary
-moment (@mymath{t\neq{t_r}}) before or after @mymath{t_r}, the proper
-distance to the object can simply be scaled with
-@mymath{a(t)}. Measuring the change of distance in a time-dependent
-(expanding) universe will also involve the speed of the object
-changing positions. Hence, let's assume that we are only thinking
-about the change in distance caused by something (light) moving at the
-speed of light. This speed is postulated as the only constant and
-frame-of-reference-independent speed in the universe, making our
-calculations easier, light is also the major source of information we
-receive from the universe, so this is a reasonable assumption for most
-extra-galactic studies. We can thus parametrize the change in distance
-as
-
-@dispmath{ds^2=c^2dt^2-a^2(t)ds_s^2 =
-c^2dt^2-a^2(t)(d\chi^2+r^2d\phi^2).}
+Therefore, @mymath{\chi(r)} is the proper distance to an object at a
+specific reference time: @mymath{t=t_r} (the @mymath{r} subscript signifies
+``reference'') when @mymath{a(t_r)=1}. At any arbitrary moment
+(@mymath{t\neq{t_r}}) before or after @mymath{t_r}, the proper distance to
+the object can be scaled with @mymath{a(t)}.
+
+Measuring the change of distance in a time-dependent (expanding) universe
+only makes sense if we can add up space and time@footnote{In other words,
+making our spacetime consistent with Minkowski spacetime geometry. In this
+geometry, different observers at a given point (event) in spacetime split
+up spacetime into `space' and `time' in different ways, just like people at
+the same spatial position can make different choices of splitting up a map
+into `left--right' and `up--down'. This model is well supported by
+twentieth and twenty-first century observations.}. But we can only add bits
+of space and time together if we measure them in the same units: with a
+conversion constant (similar to how 1000 is used to convert a kilometer
+into meters). Experimentally, we find strong support for the hypothesis
+that this conversion constant is the speed of light (or gravitational
+waves@footnote{The speed of gravitational waves was recently found to be
+very similar to that of light in vaccum, see
+@url{https://arxiv.org/abs/1710.05834, arXiv:1710.05834}.}) in a
+vacuum. This speed is postulated to be constant@footnote{In @emph{natural
+units}, speed is measured in units of the speed of light in vaccum.} and is
+almost always written as @mymath{c}. We can thus parametrize the change in
+distance on an expanding 2D surface as
+
+@dispmath{ds^2=c^2dt^2-a^2(t)ds_s^2 = c^2dt^2-a^2(t)(d\chi^2+r^2d\phi^2).}
@node Extending distance concepts to 3D, Invoking astcosmiccal, Distance on a
2D curved space, CosmicCalculator
@subsection Extending distance concepts to 3D
-The concepts of @ref{Distance on a 2D curved space} are here extended
-to a 3D space that @emph{might} be curved in a 4D space. We can start
-with the generic infinitesimal distance in a static 3D universe, but
-this time not in spherical coordinates instead of polar coordinates.
-@mymath{\theta} is shown in @ref{sphericalplane}, but here we are 3D
-beings, positioned on @mymath{O} (the center of the sphere) and the
-point @mymath{O} is tangent to a 4D-sphere. In our 3D space, a generic
-infinitesimal displacement will have the distance:
-@dispmath{ds_s^2=dx^2+dy^2+dz^2=dr^2+r^2(d\theta^2+\sin^2{\theta}d\phi^2).}Like
-our 2D friend before, we now have to assume an abstract dimension
-which we cannot visualize. Let's call the fourth dimension @mymath{w},
-then the general change in coordinates in the @emph{full} four
+The concepts of @ref{Distance on a 2D curved space} are here extended to a
+3D space that @emph{might} be curved. We can start with the generic
+infinitesimal distance in a static 3D universe, but this time in spherical
+coordinates instead of polar coordinates. @mymath{\theta} is shown in
+@ref{sphereandplane}, but here we are 3D beings, positioned on @mymath{O}
+(the center of the sphere) and the point @mymath{O} is tangent to a
+4D-sphere. In our 3D space, a generic infinitesimal displacement will
+correspond to the following distance in spherical coordinates:
+
+@dispmath{ds_s^2=dx^2+dy^2+dz^2=dr^2+r^2(d\theta^2+\sin^2{\theta}d\phi^2).}
+
+Like the 2D creature before, we now have to assume an abstract dimension
+which we cannot visualize easily. Let's call the fourth dimension
+@mymath{w}, then the general change in coordinates in the @emph{full} four
dimensional space will be:
-@dispmath{ds_s^2=dr^2+r^2(d\theta^2+\sin^2{\theta}d\phi^2)+dw^2.}But we
-can only work on a 3D curved space, so following exactly the same
+
+@dispmath{ds_s^2=dr^2+r^2(d\theta^2+\sin^2{\theta}d\phi^2)+dw^2.}
+
+@noindent
+But we can only work on a 3D curved space, so following exactly the same
steps and conventions as our 2D friend, we arrive at:
-@dispmath{ds_s^2={dr^2\over
-1-Kr^2}+r^2(d\theta^2+\sin^2{\theta}d\phi^2).}In a non-static universe
-(with a scale factor a(t), the distance can be written as:
+
+@dispmath{ds_s^2={dr^2\over 1-Kr^2}+r^2(d\theta^2+\sin^2{\theta}d\phi^2).}
+
+@noindent
+In a non-static universe (with a scale factor a(t)), the distance can be
+written as:
+
@dispmath{ds^2=c^2dt^2-a^2(t)[d\chi^2+r^2(d\theta^2+\sin^2{\theta}d\phi^2)].}
@@ -17515,6 +17619,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
@@ -17742,7 +17847,7 @@ number of threads on the system and spin-off threads.
You can see a
demonstration of using these functions in @ref{Library demo -
multi-threaded operation}.
-@deftp {C @code{struct})} gal_threads_params
+@deftp {C @code{struct}} gal_threads_params
Structure keeping the parameters of each thread. When each thread is
created, a pointer to this structure is passed to it. The @code{params}
element can be the pointer to a structure defined by the user which
@@ -20342,7 +20447,9 @@ for each dimension.
@end deftypefun
@deftypefun double gal_wcs_pixel_area_arcsec2 (struct wcsprm @code{*wcs})
-Return the pixel area of @code{wcs} in arcsecond squared.
+Return the pixel area of @code{wcs} in arcsecond squared. If the input WCS
+structure is not two dimensional and the units (@code{CUNIT} keywords) are
+not @code{deg} (for degrees), then this function will return a NaN.
@end deftypefun
@deftypefun void gal_wcs_world_to_img (struct wcsprm @code{*wcs}, double
@code{*ra}, double @code{*dec}, double @code{**x}, double @code{**y}, size_t
@code{size})
@@ -22372,7 +22479,7 @@ checking that is the most CPU intensive part will only
be done once.
@end deftypefun
-@node Git wrappers, , Interpolation, Gnuastro library
+@node Git wrappers, Cosmology library, Interpolation, Gnuastro library
@subsection Git wrappers (@file{git.h})
@cindex Git
@@ -22409,6 +22516,62 @@ controlled directory, then the output will be the
@code{NULL} pointer.
@end deftypefun
+@node Cosmology library, , Git wrappers, Gnuastro library
+@subsection 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},
+@code{o_radiation_0} 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
+@ref{CosmicCalculator}.
+
+@deftypefun 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.
+@end deftypefun
+
+@deftypefun 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.
+@end deftypefun
+
+@deftypefun 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.
+@end deftypefun
+
+@deftypefun 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
+@mymath{g/cm^3}.
+@end deftypefun
+
+@deftypefun 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.
+@end deftypefun
+
+@deftypefun 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.
+@end deftypefun
+
+@deftypefun 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).
+@end deftypefun
+
+@deftypefun 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}.
+@end deftypefun
+
+
@node Library demo programs, , Gnuastro library, Library
@section Library demo programs
diff --git a/doc/plotsrc/Makefile b/doc/plotsrc/Makefile
index 3759b9c..f377300 100644
--- a/doc/plotsrc/Makefile
+++ b/doc/plotsrc/Makefile
@@ -42,7 +42,7 @@ all.pdf: all.tex ./tex/*.tex ./conversions.sh
cp tikz/all-figure0.eps ../gnuastro-figures/iandtime.eps
cp tikz/all-figure1.eps ../gnuastro-figures/samplingfreq.eps
cp tikz/all-figure2.eps ../gnuastro-figures/flatplane.eps
- cp tikz/all-figure3.eps ../gnuastro-figures/sphericalplane.eps
+ cp tikz/all-figure3.eps ../gnuastro-figures/sphereandplane.eps
# Make all the conversions:
./conversions.sh ../gnuastro-figures/
diff --git a/doc/plotsrc/all.tex b/doc/plotsrc/all.tex
index 8d6351d..1700910 100644
--- a/doc/plotsrc/all.tex
+++ b/doc/plotsrc/all.tex
@@ -148,6 +148,6 @@ appropriate directory.
\input{tex/flatplane.tex}
-\input{tex/sphericalplane.tex}
+\input{tex/sphereandplane.tex}
\end{document}
diff --git a/doc/plotsrc/tex/sphericalplane.tex
b/doc/plotsrc/tex/sphereandplane.tex
similarity index 100%
rename from doc/plotsrc/tex/sphericalplane.tex
rename to doc/plotsrc/tex/sphereandplane.tex
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/arithmetic.c b/lib/arithmetic.c
index 6e20539..db61e9c 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -1484,6 +1484,8 @@ gal_arithmetic(int operator, unsigned char flags, ...)
out=arithmetic_abs(flags, d1);
break;
+
+ /* Multi-operand operators */
case GAL_ARITHMETIC_OP_MIN:
case GAL_ARITHMETIC_OP_MAX:
case GAL_ARITHMETIC_OP_NUM:
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 <akhlaghi@gnu.org>
+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/arithmetic.h b/lib/gnuastro/arithmetic.h
index 5d50518..ba20c6a 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -132,6 +132,8 @@ enum gal_arithmetic_operators
GAL_ARITHMETIC_OP_TO_INT64, /* Convert to int64_t. */
GAL_ARITHMETIC_OP_TO_FLOAT32, /* Convert to float32. */
GAL_ARITHMETIC_OP_TO_FLOAT64, /* Convert to float64. */
+
+ GAL_ARITHMETIC_OP_LAST_CODE, /* Last code of the library operands. */
};
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 <akhlaghi@gnu.org>
+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__ */
diff --git a/lib/statistics.c b/lib/statistics.c
index 6d5e013..a85ce88 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -1805,7 +1805,7 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
meanstd=gal_statistics_mean_std(nbs);
/* Put the three final values in usable (with a type) pointers. */
- med = median_d->array;
+ med = median_d->array;
mean = meanstd->array;
std = &((double *)(meanstd->array))[1];
diff --git a/lib/wcs.c b/lib/wcs.c
index 80167a1..0649e57 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -581,9 +581,13 @@ gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs)
/* A small sanity check. Later, when higher dimensions are necessary, we
can find which ones correlate to RA and Dec and use them to find the
pixel area in arcsec^2. */
- if(wcs->naxis!=2)
- error(EXIT_FAILURE, 0, "%s: currently only 2D datasets supported. "
- "The input WCS has %d dimensions", __func__, wcs->naxis);
+ if(wcs->naxis!=2) return NAN;
+
+ /* Check if the units of the axis are degrees or not. Currently all FITS
+ images I have worked with use `deg' for degrees. If other alternatives
+ exist, we can add corrections later. */
+ if( strcmp("deg", wcs->cunit[0]) || strcmp("deg", wcs->cunit[1]) )
+ return NAN;
/* Get the pixel scales along each axis in degrees, then multiply. */
pixscale=gal_wcs_pixel_scale(wcs);
@@ -612,6 +616,7 @@ gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs)
+
/**************************************************************/
/********** Array conversion ************/
/**************************************************************/
- [gnuastro-commits] master 3e76e1f 010/113: MakeProfiles --kernel builds 3D kernels also, (continued)
- [gnuastro-commits] master 3e76e1f 010/113: MakeProfiles --kernel builds 3D kernels also, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 81f3f65 023/113: More --coordcol options acceptable in Crop, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master e47f8db 024/113: Merged recent work in master, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 8fa5ff1 026/113: Minor edit in book (part added in last commit), Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 40f0a56 013/113: Minor corrections to MakeProfiles continued, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 8372486 020/113: NoiseChisel's detection complete in 3D, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 6cc3d25 027/113: No 3D projections in function to inspect NoiseChisel outputs, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master e990860 029/113: Merged recent updates in master, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 6bdc5d6 030/113: Oversegmentation connectivity one less for 3D, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master be3bfd8 033/113: NoiseChisel configuration file in 3D updated, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 3b6c15d 036/113: Merged with recent work in master, no conflicts,
Mohammad Akhlaghi <=
- [gnuastro-commits] master fb3660f 037/113: MakeCatalog works in 3D, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 73fdf0c 006/113: MakeProfiles builds 3D ellipsoids, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 03d73fe 011/113: Some minor corrections in code and comments, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 87ab805 014/113: Identifiers for integer constants in BZERO comparisons, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 1483201 009/113: Noised cube created in make check, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master bdeaba9 012/113: Corrected name of 3D catalog for tarball inclusion, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 6b53c5e 015/113: Bug fixes in master: commits 1ff1c25 and 540af65, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master adeb1c6 016/113: NoiseChisel's convolution step in 3D completed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master ace4160 017/113: New --mcolisbrightness option for MakeProfiles, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master f06acc5 018/113: Merged recent work from master, Mohammad Akhlaghi, 2021/04/16