1This is gnuastro.info, produced by makeinfo version 6.8 from
2gnuastro.texi.
3
4This book documents version 0.16 of the GNU Astronomy Utilities
5(Gnuastro).  Gnuastro provides various programs and libraries for
6astronomical data manipulation and analysis.
7
8   Copyright © 2015-2021, Free Software Foundation, Inc.
9
10     Permission is granted to copy, distribute and/or modify this
11     document under the terms of the GNU Free Documentation License,
12     Version 1.3 or any later version published by the Free Software
13     Foundation; with no Invariant Sections, no Front-Cover Texts, and
14     no Back-Cover Texts.  A copy of the license is included in the
15     section entitled “GNU Free Documentation License”.
16INFO-DIR-SECTION Astronomy
17START-INFO-DIR-ENTRY
18* Gnuastro: (gnuastro).       GNU Astronomy Utilities.
19* libgnuastro: (gnuastro)Gnuastro library. Full Gnuastro library doc.
20
21* help-gnuastro: (gnuastro)help-gnuastro mailing list. Getting help.
22
23* bug-gnuastro: (gnuastro)Report a bug. How to report bugs
24
25* Arithmetic: (gnuastro)Arithmetic. Arithmetic operations on pixels.
26* astarithmetic: (gnuastro)Invoking astarithmetic. Options to Arithmetic.
27
28* BuildProgram: (gnuastro)BuildProgram. Compile and run programs using Gnuastro’s library.
29* astbuildprog: (gnuastro)Invoking astbuildprog. Options to BuildProgram.
30
31* ConvertType: (gnuastro)ConvertType. Convert different file types.
32* astconvertt: (gnuastro)Invoking astconvertt. Options to ConvertType.
33
34* Convolve: (gnuastro)Convolve. Convolve an input file with kernel.
35* astconvolve: (gnuastro)Invoking astconvolve. Options to Convolve.
36
37* CosmicCalculator: (gnuastro)CosmicCalculator. For cosmological params.
38* astcosmiccal: (gnuastro)Invoking astcosmiccal. Options to CosmicCalculator.
39
40* Crop: (gnuastro)Crop. Crop region(s) from image(s).
41* astcrop: (gnuastro)Invoking astcrop. Options to Crop.
42
43* Fits: (gnuastro)Fits. View and manipulate FITS extensions and keywords.
44* astfits: (gnuastro)Invoking astfits. Options to Fits.
45
46* MakeCatalog: (gnuastro)MakeCatalog. Make a catalog from labeled image.
47* astmkcatalog: (gnuastro)Invoking astmkcatalog. Options to MakeCatalog.
48
49* MakeNoise: (gnuastro)MakeNoise. Make (add) noise to an image.
50* astmknoise: (gnuastro)Invoking astmknoise. Options to MakeNoise.
51
52* MakeProfiles: (gnuastro)MakeProfiles. Make mock profiles.
53* astmkprof: (gnuastro)Invoking astmkprof. Options to MakeProfiles.
54
55* Match: (gnuastro)Match. Match two separate catalogs.
56* astmatch: (gnuastro)Invoking astmatch. Options to Match.
57
58* NoiseChisel: (gnuastro)NoiseChisel. Detect signal in noise.
59* astnoisechisel: (gnuastro)Invoking astnoisechisel. Options to NoiseChisel.
60
61* Segment: (gnuastro)Segment. Segment detections based on signal structure.
62* astsegment: (gnuastro)Invoking astsegment. Options to Segment.
63
64* Query: (gnuastro)Query. Access remote databases for downloading data.
65* astquery: (gnuastro)Invoking astquery. Options to Query.
66
67* Statistics: (gnuastro)Statistics. Get image Statistics.
68* aststatistics: (gnuastro)Invoking aststatistics. Options to Statistics.
69
70* Table: (gnuastro)Table. Read and write FITS binary or ASCII tables.
71* asttable: (gnuastro)Invoking asttable. Options to Table.
72
73* Warp: (gnuastro)Warp. Warp a dataset to a new grid.
74* astwarp: (gnuastro)Invoking astwarp. Options to Warp.
75
76* astscript: (gnuastro)Installed scripts. Gnuastro’s installed scripts.
77* astscript-sort-by-night: (gnuastro)Invoking astscript-sort-by-night. Options to this script
78* astscript-radial-profile: (gnuastro)Invoking astscript-radial-profile. Options to this script
79* astscript-ds9-region: (gnuastro)Invoking astscript-ds9-region. Options to this script
80
81END-INFO-DIR-ENTRY
82
83
84File: gnuastro.info,  Node: Surface brightness limit of image,  Next: Upper limit magnitude of image,  Prev: Upper limit magnitude of each detection,  Up: Quantifying measurement limits
85
867.4.3.5 Surface brightness limit of image
87.........................................
88
89As we make more observations on one region of the sky and add/combine
90the observations into one dataset, both the signal and the noise
91increase.  However, the signal increases much faster than the noise:
92Assuming you add $N$ datasets with equal exposure times, the signal will
93increases as a multiple of $N$, while noise increases as $\sqrt{N}$.
94Therefore the signal-to-noise ratio increases by a factor of $\sqrt{N}$.
95Visually, fainter (per pixel) parts of the objects/signal in the image
96will become more visible/detectable.  The noise-level is known as the
97dataset’s surface brightness limit.
98
99   You can think of the noise as muddy water that is completely covering
100a flat ground(1).  The signal (coming from astronomical objects in real
101data) will be summits/hills that start from the flat sky level (under
102the muddy water) and their summits can sometimes reach above the muddy
103water.  Let’s assume that in your first observation the muddy water has
104just been stirred and except a few small peaks, you can’t see anything
105through the mud.  As you wait and make more observations/exposures, the
106mud settles down and the _depth_ of the transparent water increases.  As
107a result, more and more summits become visible and the lower parts of
108the hills (parts with lower surface brightness) can be seen more
109clearly.  In this analogy(2), height (from the ground) is the _surface
110brightness_ and the height of the muddy water at the moment you combine
111your data, is your _surface brightness limit_ for that moment.
112
113   The outputs of NoiseChisel include the Sky standard deviation
114($\sigma$) on every group of pixels (a tile) that were calculated from
115the undetected pixels in each tile, see *note Tessellation:: and *note
116NoiseChisel output::.  Let’s take $\sigma_m$ as the median $\sigma$ over
117the successful meshes in the image (prior to interpolation or
118smoothing).  It is recorded in the ‘MEDSTD’ keyword of the ‘SKY_STD’
119extension of NoiseChisel’s output.
120
121   On different instruments, pixels cover different spatial angles over
122the sky.  For example, the width of each pixel on the ACS camera on the
123Hubble Space Telescope (HST) is roughly 0.05 seconds of arc, while the
124pixels of SDSS are each 0.396 seconds of arc (almost eight times
125wider(3)).  Nevertheless, irrespective of its sky coverage, a pixel is
126our unit of data collection.
127
128   To start with, we define the low-level Surface brightness limit or
129_depth_, in units of magnitude/pixel with the equation below (assuming
130the image has zero point magnitude $z$ and we want the $n$th multiple of
131$\sigma_m$).
132
133$$SB_{n\sigma,\rm pixel}=-2.5\times\log_{10}{(n\sigma_m)}+z \quad\quad [mag/pixel]$$
134
135   As an example, the XDF survey covers part of the sky that the HST has
136observed the most (for 85 orbits) and is consequently very small
137($\sim4$ minutes of arc, squared).  On the other hand, the CANDELS
138survey, is one of the widest multi-color surveys done by the HST
139covering several fields (about 720 arcmin$^2$) but its deepest fields
140have only 9 orbits observation.  The $1\sigma$ depth of the XDF and
141CANDELS-deep surveys in the near infrared WFC3/F160W filter are
142respectively 34.40 and 32.45 magnitudes/pixel.  In a single orbit image,
143this same field has a $1\sigma$ depth of 31.32 magnitudes/pixel.  Recall
144that a larger magnitude corresponds to less brightness, see *note
145Brightness flux magnitude::.
146
147   The low-level magnitude/pixel measurement above is only useful when
148all the datasets you want to use, or compare, have the same pixel size.
149However, you will often find yourself using, or comparing, datasets from
150various instruments with different pixel scales (projected pixel width,
151in arc-seconds).  If we know the pixel scale, we can obtain a more
152easily comparable surface brightness limit in units of:
153magnitude/arcsec$^2$.  But another complication is that astronomical
154objects are usually larger than 1 arcsec$^2$.  As a result, it is common
155to measure the surface brightness limit over a larger (but fixed,
156depending on context) area.
157
158   Let’s assume that every pixel is $p$ arcsec$^2$ and we want the
159surface brightness limit for an object covering A arcsec$^2$ (so $A/p$
160is the number of pixels that cover an area of $A$ arcsec$^2$).  On the
161other hand, noise is added in RMS(4), hence the noise level in $A$
162arcsec$^2$ is $n\sigma_m\sqrt{A/p}$.  But we want the result in units of
163arcsec$^2$, so we should divide this by $A$ arcsec$^2$:
164$n\sigma_m\sqrt{A/p}/A=n\sigma_m\sqrt{A/(pA^2)}=n\sigma_m/\sqrt{pA}$.
165Plugging this into the magnitude equation, we get the $n\sigma$ surface
166brightness limit, over an area of A arcsec$^2$, in units of
167magnitudes/arcsec$^2$:
168
169$$SB_{{n\sigma,\rm A arcsec}^2}=-2.5\times\log_{10}{\left(n\sigma_m\over \sqrt{pA}\right)+z} \quad\quad [mag/arcsec^2]$$
170
171   MakeCatalog will calculate the input dataset’s $SB_{n\sigma,\rm
172pixel}$ and $SB_{{n\sigma,\rm A arcsec}^2}$ and will write them as the
173‘SBLMAGPIX’ and ‘SBLMAG’ keywords the output catalog(s), see *note
174MakeCatalog output::.  You can set your desired $n$-th multiple of
175$\sigma$ and the $A$ arcsec$^2$ area using the following two options
176respectively: ‘--sfmagnsigma’ and ‘--sfmagarea’ (see *note MakeCatalog
177output::).  Just note that $SB_{{n\sigma,\rm A arcsec}^2}$ is only
178calculated if the input has World Coordinate System (WCS). Without WCS,
179the pixel scale cannot be derived.
180
181   As you saw in its derivation, the calculation above extrapolates the
182noise in one pixel over all the input’s pixels!  It therefore implicitly
183assumes that the noise is the same in all of the pixels.  But this only
184happens in individual exposures: reduced data will have correlated noise
185because they are a stack of many individual exposures that have been
186warped (thus mixing the pixel values).  A more accurate measure which
187will provide a realistic value for every labeled region is known as the
188_upper-limit magnitude_, which is discussed below.
189
190   ---------- Footnotes ----------
191
192   (1) The ground is the sky value in this analogy, see *note Sky
193value::.  Note that this analogy only holds for a flat sky value across
194the surface of the image or ground.
195
196   (2) Note that this muddy water analogy is not perfect, because while
197the water-level remains the same all over a peak, in data analysis, the
198Poisson noise increases with the level of data.
199
200   (3) Ground-based instruments like the SDSS suffer from strong
201smoothing due to the atmosphere.  Therefore, increasing the pixel
202resolution (or decreasing the width of a pixel) won’t increase the
203received information).
204
205   (4) If you add three datasets with noise $\sigma_1$, $\sigma_2$ and
206$\sigma_3$, the resulting noise level is
207$\sigma_t=\sqrt{\sigma_1^2+\sigma_2^2+\sigma_3^2}$, so when
208$\sigma_1=\sigma_2=\sigma_3\equiv\sigma$, then
209$\sigma_t=\sigma\sqrt{3}$.  In this case, the area $A$ is covered by
210$A/p$ pixels, so the noise level is $\sigma_t=\sigma\sqrt{A/p}$.
211
212
213File: gnuastro.info,  Node: Upper limit magnitude of image,  Prev: Surface brightness limit of image,  Up: Quantifying measurement limits
214
2157.4.3.6 Upper limit magnitude of image
216......................................
217
218As mentioned above, the upper-limit magnitude will depend on the shape
219of each object’s footprint.  Therefore we can measure the dataset’s
220upper-limit magnitude using standard shapes.  Traditionally a circular
221aperture of a fixed size (in arcseconds) has been used.  For a full
222example of implementing this, see the respective section in the tutorial
223(*note Image surface brightness limit::).
224
225
226File: gnuastro.info,  Node: Measuring elliptical parameters,  Next: Adding new columns to MakeCatalog,  Prev: Quantifying measurement limits,  Up: MakeCatalog
227
2287.4.4 Measuring elliptical parameters
229-------------------------------------
230
231The shape or morphology of a target is one of the most commonly desired
232parameters of a target.  Here, we will review the derivation of the most
233basic/simple morphological parameters: the elliptical parameters for a
234set of labeled pixels.  The elliptical parameters are: the (semi-)major
235axis, the (semi-)minor axis and the position angle along with the
236central position of the profile.  The derivations below follow the
237SExtractor manual derivations with some added explanations for easier
238reading.
239
240   Let’s begin with one dimension for simplicity: Assume we have a set
241of $N$ values $B_i$ (for example showing the spatial distribution of a
242target’s brightness), each at position $x_i$.  The simplest parameter we
243can define is the geometric center of the object ($x_g$) (ignoring the
244brightness values): $x_g=(\sum_ix_i)/N$.  _Moments_ are defined to
245incorporate both the value (brightness) and position of the data.  The
246first moment can be written as:
247
248            $$\overline{x}={\sum_iB_ix_i \over \sum_iB_i}$$
249
250This is essentially the weighted (by $B_i$) mean position.  The
251geometric center ($x_g$, defined above) is a special case of this with
252all $B_i=1$.  The second moment is essentially the variance of the
253distribution:
254
255      $$\overline{x^2}\equiv{\sum_iB_i(x_i-\overline{x})^2 \over
256            \sum_iB_i} = {\sum_iB_ix_i^2 \over \sum_iB_i} -
257      2\overline{x}{\sum_iB_ix_i\over\sum_iB_i} + \overline{x}^2
258         ={\sum_iB_ix_i^2 \over \sum_iB_i} - \overline{x}^2$$
259
260The last step was done from the definition of $\overline{x}$.  Hence,
261the square root of $\overline{x^2}$ is the spatial standard deviation
262(along the one-dimension) of this particular brightness distribution
263($B_i$).  Crudely (or qualitatively), you can think of its square root
264as the distance (from $\overline{x}$) which contains a specific amount
265of the flux (depending on the $B_i$ distribution).  Similar to the first
266moment, the geometric second moment can be found by setting all $B_i=1$.
267So while the first moment quantified the position of the brightness
268distribution, the second moment quantifies how that brightness is
269dispersed about the first moment.  In other words, it quantifies how
270“sharp” the object’s image is.
271
272   Before continuing to two dimensions and the derivation of the
273elliptical parameters, let’s pause for an important implementation
274technicality.  You can ignore this paragraph and the next two if you
275don’t want to implement these concepts.  The basic definition (first
276definition of $\overline{x^2}$ above) can be used without any major
277problem.  However, using this fraction requires two runs over the data:
278one run to find $\overline{x}$ and another run to find $\overline{x^2}$
279from $\overline{x}$, this can be slow.  The advantage of the last
280fraction above, is that we can estimate both the first and second
281moments in one run (since the $-\overline{x}^2$ term can easily be added
282later).
283
284   The logarithmic nature of floating point number digitization creates
285a complication however: suppose the object is located between pixels
28610000 and 10020.  Hence the target’s pixels are only distributed over 20
287pixels (with a standard deviation $<20$), while the mean has a value of
288$\sim10000$.  The $\sum_iB_i^2x_i^2$ will go to very very large values
289while the individual pixel differences will be orders of magnitude
290smaller.  This will lower the accuracy of our calculation due to the
291limited accuracy of floating point operations.  The variance only
292depends on the distance of each point from the mean, so we can shift all
293position by a constant/arbitrary $K$ which is much closer to the mean:
294$\overline{x-K}=\overline{x}-K$.  Hence we can calculate the second
295order moment using:
296
297        $$\overline{x^2}={\sum_iB_i(x_i-K)^2 \over \sum_iB_i} -
298                         (\overline{x}-K)^2 $$
299
300The closer $K$ is to $\overline{x}$, the better (the sums of squares
301will involve smaller numbers), as long as $K$ is within the object
302limits (in the example above: $10000\leq{K}\leq10020$), the floating
303point error induced in our calculation will be negligible.  For the most
304simplest implementation, MakeCatalog takes $K$ to be the smallest
305position of the object in each dimension.  Since $K$ is arbitrary and an
306implementation/technical detail, we will ignore it for the remainder of
307this discussion.
308
309   In two dimensions, the mean and variances can be written as:
310
311             $$\overline{x}={\sum_iB_ix_i\over B_i}, \quad
312           \overline{x^2}={\sum_iB_ix_i^2 \over \sum_iB_i} -
313                           \overline{x}^2$$
314             $$\overline{y}={\sum_iB_iy_i\over B_i}, \quad
315           \overline{y^2}={\sum_iB_iy_i^2 \over \sum_iB_i} -
316                           \overline{y}^2$$
317            $$\quad\quad\quad\quad\quad\quad\quad\quad\quad
318           \overline{xy}={\sum_iB_ix_iy_i \over \sum_iB_i} -
319                   \overline{x}\times\overline{y}$$
320
321   If an elliptical profile’s major axis exactly lies along the $x$
322axis, then $\overline{x^2}$ will be directly proportional with the
323profile’s major axis, $\overline{y^2}$ with its minor axis and
324$\overline{xy}=0$.  However, in reality we are not that lucky and
325(assuming galaxies can be parameterized as an ellipse) the major axis of
326galaxies can be in any direction on the image (in fact this is one of
327the core principles behind weak-lensing by shear estimation).  So the
328purpose of the remainder of this section is to define a strategy to
329measure the position angle and axis ratio of some randomly positioned
330ellipses in an image, using the raw second moments that we have
331calculated above in our image coordinates.
332
333   Let’s assume we have rotated the galaxy by $\theta$, the new second
334order moments are:
335
336        $$\overline{x_\theta^2} = \overline{x^2}\cos^2\theta +
337                     \overline{y^2}\sin^2\theta -
338                 2\overline{xy}\cos\theta\sin\theta $$
339        $$\overline{y_\theta^2} = \overline{x^2}\sin^2\theta +
340                     \overline{y^2}\cos^2\theta +
341                 2\overline{xy}\cos\theta\sin\theta$$
342     $$\overline{xy_\theta} = \overline{x^2}\cos\theta\sin\theta -
343                 \overline{y^2}\cos\theta\sin\theta +
344              \overline{xy}(\cos^2\theta-\sin^2\theta)$$
345
346The best $\theta$ ($\theta_0$, where major axis lies along the
347$x_\theta$ axis) can be found by:
348
349$$\left.{\partial \overline{x_\theta^2} \over \partial \theta}\right|_{\theta_0}=0$$
350   Taking the derivative, we get:
351     $$2\cos\theta_0\sin\theta_0(\overline{y^2}-\overline{x^2}) +
352        2(\cos^2\theta_0-\sin^2\theta_0)\overline{xy}=0$$ When
353   $\overline{x^2}\neq\overline{y^2}$, we can write:
354                           $$\tan2\theta_0 =
355        2{\overline{xy} \over \overline{x^2}-\overline{y^2}}.$$
356
357MakeCatalog uses the standard C math library’s ‘atan2’ function to
358estimate $\theta_0$, which we define as the position angle of the
359ellipse.  To recall, this is the angle of the major axis of the ellipse
360with the $x$ axis.  By definition, when the elliptical profile is
361rotated by $\theta_0$, then $\overline{xy_{\theta_0}}=0$,
362$\overline{x_{\theta_0}^2}$ will be the extent of the maximum variance
363and $\overline{y_{\theta_0}^2}$ the extent of the minimum variance
364(which are perpendicular for an ellipse).  Replacing $\theta_0$ in the
365equations above for $\overline{x_\theta}$ and $\overline{y_\theta}$, we
366can get the semi-major ($A$) and semi-minor ($B$) lengths:
367
368        $$A^2\equiv\overline{x_{\theta_0}^2}= {\overline{x^2} +
369\overline{y^2} \over 2} + \sqrt{\left({\overline{x^2}-\overline{y^2} \over 2}\right)^2 + \overline{xy}^2}$$
370
371        $$B^2\equiv\overline{y_{\theta_0}^2}= {\overline{x^2} +
372\overline{y^2} \over 2} - \sqrt{\left({\overline{x^2}-\overline{y^2} \over 2}\right)^2 + \overline{xy}^2}$$
373
374   As a summary, it is important to remember that the units of $A$ and
375$B$ are in pixels (the standard deviation of a positional distribution)
376and that they represent the spatial light distribution of the object in
377both image dimensions (rotated by $\theta_0$).  When the object cannot
378be represented as an ellipse, this interpretation breaks down:
379$\overline{xy_{\theta_0}}\neq0$ and $\overline{y_{\theta_0}^2}$ will not
380be the direction of minimum variance.
381
382
383File: gnuastro.info,  Node: Adding new columns to MakeCatalog,  Next: Invoking astmkcatalog,  Prev: Measuring elliptical parameters,  Up: MakeCatalog
384
3857.4.5 Adding new columns to MakeCatalog
386---------------------------------------
387
388MakeCatalog is designed to allow easy addition of different measurements
389over a labeled image (see Akhlaghi [2016]
390(https://arxiv.org/abs/1611.06387v1)).  A check-list style description
391of necessary steps to do that is described in this section.  The common
392development characteristics of MakeCatalog and other Gnuastro programs
393is explained in *note Developing::.  We strongly encourage you to have a
394look at that chapter to greatly simplify your navigation in the code.
395After adding and testing your column, you are most welcome (and
396encouraged) to share it with us so we can add to the next release of
397Gnuastro for everyone else to also benefit from your efforts.
398
399   MakeCatalog will first pass over each label’s pixels two times and do
400necessary raw/internal calculations.  Once the passes are done, it will
401use the raw information for filling the final catalog’s columns.  In the
402first pass it will gather mainly object information and in the second
403run, it will mainly focus on the clumps, or any other measurement that
404needs an output from the first pass.  These two passes are designed to
405be raw summations: no extra processing.  This will allow parallel
406processing and simplicity/clarity.  So if your new calculation, needs
407new raw information from the pixels, then you will need to also modify
408the respective ‘mkcatalog_first_pass’ and ‘mkcatalog_second_pass’
409functions (both in ‘bin/mkcatalog/mkcatalog.c’) and define new raw table
410columns in ‘main.h’ (hopefully the comments in the code are clear
411enough).
412
413   In all these different places, the final columns are sorted in the
414same order (same order as *note Invoking astmkcatalog::).  This allows a
415particular column/option to be easily found in all steps.  Therefore in
416adding your new option, be sure to keep it in the same relative place in
417the list in all the separate places (it doesn’t necessarily have to be
418in the end), and near conceptually similar options.
419
420main.h421     The ‘objectcols’ and ‘clumpcols’ enumerated variables (‘enum’)
422     define the raw/internal calculation columns.  If your new column
423     requires new raw calculations, add a row to the respective list.
424     If your calculation requires any other settings parameters, you
425     should add a variable to the ‘mkcatalogparams’ structure.
426
427ui.c428     If the new column needs raw calculations (an entry was added in
429     ‘objectcols’ and ‘clumpcols’), specify which inputs it needs in
430     ‘ui_necessary_inputs’, similar to the other options.  Afterwards,
431     if your column includes any particular settings (you needed to add
432     a variable to the ‘mkcatalogparams’ structure in ‘main.h’), you
433     should do the sanity checks and preparations for it here.
434
435ui.h436     The ‘option_keys_enum’ associates a unique value for each option to
437     MakeCatalog.  The options that have a short option version, the
438     single character short comment is used for the value.  Those that
439     don’t have a short option version, get a large integer
440     automatically.  You should add a variable here to identify your
441     desired column.
442
443args.h444     This file specifies all the parameters for the GNU C library, Argp
445     structure that is in charge of reading the user’s options.  To
446     define your new column, just copy an existing set of parameters and
447     change the first, second and 5th values (the only ones that differ
448     between all the columns), you should use the macro you defined in
449ui.h’ here.
450
451columns.c452     This file contains the main definition and high-level calculation
453     of your new column through the ‘columns_define_alloc’ and
454     ‘columns_fill’ functions.  In the first, you specify the basic
455     information about the column: its name, units, comments, type (see
456     *note Numeric data types::) and how it should be printed if the
457     output is a text file.  You should also specify the raw/internal
458     columns that are necessary for this column here as the many
459     existing examples show.  Through the types for objects and rows,
460     you can specify if this column is only for clumps, objects or both.
461
462     The second main function (‘columns_fill’) writes the final value
463     into the appropriate column for each object and clump.  As you can
464     see in the many existing examples, you can define your processing
465     on the raw/internal calculations here and save them in the output.
466
467mkcatalog.c468     This file contains the low-level parsing functions.  To be
469     optimized, the parsing is done in parallel through the
470     ‘mkcatalog_single_object’ function.  This function initializes the
471     necessary arrays and calls the lower-level ‘parse_objects’ and
472     ‘parse_clumps’ for actually going over the pixels.  They are all
473     heavily commented, so you should be able to follow where to add
474     your necessary low-level calculations.
475
476doc/gnuastro.texi477     Update this manual and add a description for the new column.
478
479
480File: gnuastro.info,  Node: Invoking astmkcatalog,  Prev: Adding new columns to MakeCatalog,  Up: MakeCatalog
481
4827.4.6 Invoking MakeCatalog
483--------------------------
484
485MakeCatalog will do measurements and produce a catalog from a labeled
486dataset and optional values dataset(s).  The executable name is
487‘astmkcatalog’ with the following general template
488
489     $ astmkcatalog [OPTION ...] InputImage.fits
490
491One line examples:
492
493     ## Create catalog with RA, Dec, Magnitude and Magnitude error,
494     ## from Segment's output:
495     $ astmkcatalog --ra --dec --magnitude --magnitudeerr seg-out.fits
496
497     ## Same catalog as above (using short options):
498     $ asmkcatalog -rdmG seg-out.fits
499
500     ## Write the catalog to a text table:
501     $ astmkcatalog -mpQ seg-out.fits --output=cat.txt
502
503     ## Output columns specified in `columns.conf':
504     $ astmkcatalog --config=columns.conf seg-out.fits
505
506     ## Use object and clump labels from a K-band image, but pixel values
507     ## from an i-band image.
508     $ astmkcatalog K_segmented.fits --hdu=DETECTIONS --clumpscat     \
509                    --clumpsfile=K_segmented.fits --clumpshdu=CLUMPS  \
510                    --valuesfile=i_band.fits
511
512If MakeCatalog is to do processing (not printing help or option values),
513an input labeled image should be provided.  The options described in
514this section are those that are particular to MakeProfiles.  For
515operations that MakeProfiles shares with other programs (mainly
516involving input/output or general processing steps), see *note Common
517options::.  Also see *note Common program behavior:: for some general
518characteristics of all Gnuastro programs including MakeCatalog.
519
520   The various measurements/columns of MakeCatalog are requested as
521options, either on the command-line or in configuration files, see *note
522Configuration files::.  The full list of available columns is available
523in *note MakeCatalog measurements::.  Depending on the requested
524columns, MakeCatalog needs more than one input dataset, for more
525details, please see *note MakeCatalog inputs and basic settings::.  The
526upper-limit measurements in particular need several configuration
527options which are thoroughly discussed in *note Upper-limit settings::.
528Finally, in *note MakeCatalog output:: the output file(s) created by
529MakeCatalog are discussed.
530
531* Menu:
532
533* MakeCatalog inputs and basic settings::  Input files and basic settings.
534* Upper-limit settings::        Settings for upper-limit measurements.
535* MakeCatalog measurements::    Available measurements in MakeCatalog.
536* MakeCatalog output::          File names of MakeCatalog’s output table.
537
538
539File: gnuastro.info,  Node: MakeCatalog inputs and basic settings,  Next: Upper-limit settings,  Prev: Invoking astmkcatalog,  Up: Invoking astmkcatalog
540
5417.4.6.1 MakeCatalog inputs and basic settings
542.............................................
543
544MakeCatalog works by using a localized/labeled dataset (see *note
545MakeCatalog::).  This dataset maps/labels pixels to a specific target
546(row number in the final catalog) and is thus the only necessary input
547dataset to produce a minimal catalog in any situation.  Because it only
548has labels/counters, it must have an integer type (see *note Numeric
549data types::), see below if your labels are in a floating point
550container.  When the requested measurements only need this dataset (for
551example ‘--geox’, ‘--geoy’, or ‘--geoarea’), MakeCatalog won’t read any
552more datasets.
553
554   Low-level measurements that only use the labeled image are rarely
555sufficient for any high-level science case.  Therefore necessary input
556datasets depend on the requested columns in each run.  For example,
557let’s assume you want the brightness/magnitude and signal-to-noise ratio
558of your labeled regions.  For these columns, you will also need to
559provide an extra dataset containing values for every pixel of the
560labeled input (to measure brightness) and another for the Sky standard
561deviation (to measure error).  All such auxiliary input files have to
562have the same size (number of pixels in each dimension) as the input
563labeled image.  Their numeric data type is irrelevant (they will be
564converted to 32-bit floating point internally).  For the full list of
565available measurements, see *note MakeCatalog measurements::.
566
567   The “values” dataset is used for measurements like
568brightness/magnitude, or flux-weighted positions.  If it is a real
569image, by default it is assumed to be already Sky-subtracted prior to
570running MakeCatalog.  If it isn’t, you use the ‘--subtractsky’ option
571to, so MakeCatalog reads and subtracts the Sky dataset before any
572processing.  To obtain the Sky value, you can use the ‘--sky’ option of
573*note Statistics::, but the best recommended method is *note
574NoiseChisel::, see *note Sky value::.
575
576   MakeCatalog can also do measurements on sub-structures of detections.
577In other words, it can produce two catalogs.  Following the nomenclature
578of Segment (see *note Segment::), the main labeled input dataset is
579known as “object” labels and the (optional) sub-structure input dataset
580is known as “clumps”.  If MakeCatalog is run with the ‘--clumpscat’
581option, it will also need a labeled image containing clumps, similar to
582what Segment produces (see *note Segment output::).  Since clumps are
583defined within detected regions (they exist over signal, not noise),
584MakeCatalog uses their boundaries to subtract the level of signal under
585them.
586
587   There are separate options to explicitly request a file name and
588HDU/extension for each of the required input datasets as fully described
589below (with the ‘--*file’ format).  When each dataset is in a separate
590file, these options are necessary.  However, one great advantage of the
591FITS file format (that is heavily used in astronomy) is that it allows
592the storage of multiple datasets in one file.  So in most situations
593(for example if you are using the outputs of *note NoiseChisel:: or
594*note Segment::), all the necessary input datasets can be in one file.
595
596   When none of the ‘--*file’ options are given, MakeCatalog will assume
597the necessary input datasets are in the file given as its argument
598(without any option).  When the Sky or Sky standard deviation datasets
599are necessary and the only ‘--*file’ option called is ‘--valuesfile’,
600MakeCatalog will search for these datasets (with the default/given HDUs)
601in the file given to ‘--valuesfile’ (before looking into the main
602argument file).
603
604   When the clumps image (necessary with the ‘--clumpscat’ option) is
605used, MakeCatalog looks into the (possibly existing) ‘NUMLABS’ keyword
606for the total number of clumps in the image (irrespective of how many
607objects there are).  If its not present, it will count them and possibly
608re-label the clumps so the clump labels always start with 1 and finish
609with the total number of clumps in each object.  The re-labeled clumps
610image will be stored with the ‘-clumps-relab.fits’ suffix.  This can
611slightly slow-down the run.
612
613   Note that ‘NUMLABS’ is automatically written by Segment in its
614outputs, so if you are feeding Segment’s clump labels, you can benefit
615from the improved speed.  Otherwise, if you are creating the clumps
616label dataset manually, it may be good to include the ‘NUMLABS’ keyword
617in its header and also be sure that there is no gap in the clump labels.
618For example if an object has three clumps, they are labeled as 1, 2, 3.
619If they are labeled as 1, 3, 4, or any other combination of three
620positive integers that aren’t an increment of the previous, you might
621get unknown behavior.
622
623   It may happen that your labeled objects image was created with a
624program that only outputs floating point files.  However, you know it
625only has integer valued pixels that are stored in a floating point
626container.  In such cases, you can use Gnuastro’s Arithmetic program
627(see *note Arithmetic::) to change the numerical data type of the image
628(‘float.fits’) to an integer type image (‘int.fits’) with a command like
629below:
630
631     $ astarithmetic float.fits int32 --output=int.fits
632
633   To summarize: if the input file to MakeCatalog is the default/full
634output of Segment (see *note Segment output::) you don’t have to worry
635about any of the ‘--*file’ options below.  You can just give Segment’s
636output file to MakeCatalog as described in *note Invoking
637astmkcatalog::.  To feed NoiseChisel’s output into MakeCatalog, just
638change the labeled dataset’s header (with ‘--hdu=DETECTIONS’).  The full
639list of input dataset options and general setting options are described
640below.
641
642‘-l FITS’
643‘--clumpsfile=FITS’
644     The FITS file containing the labeled clumps dataset when
645     ‘--clumpscat’ is called (see *note MakeCatalog output::).  When
646     ‘--clumpscat’ is called, but this option isn’t, MakeCatalog will
647     look into the main input file (given as an argument) for the
648     required extension/HDU (value to ‘--clumpshdu’).
649
650‘--clumpshdu=STR’
651     The HDU/extension of the clump labels dataset.  Only pixels with
652     values above zero will be considered.  The clump labels dataset has
653     to be an integer data type (see *note Numeric data types::) and
654     only pixels with a value larger than zero will be used.  See *note
655     Segment output:: for a description of the expected format.
656
657‘-v FITS’
658‘--valuesfile=FITS’
659     The file name of the (sky-subtracted) values dataset.  When any of
660     the columns need values to associate with the input labels (for
661     example to measure the brightness/magnitude of a galaxy),
662     MakeCatalog will look into a “values” for the respective pixel
663     values.  In most common processing, this is the actual astronomical
664     image that the labels were defined, or detected, over.  The
665     HDU/extension of this dataset in the given file can be specified
666     with ‘--valueshdu’.  If this option is not called, MakeCatalog will
667     look for the given extension in the main input file.
668
669‘--valueshdu=STR/INT670     The name or number (counting from zero) of the extension containing
671     the “values” dataset, see the descriptions above and those in
672     ‘--valuesfile’ for more.
673
674‘-s FITS/FLT675‘--insky=FITS/FLT676     Sky value as a single number, or the file name containing a dataset
677     (different values per pixel or tile).  The Sky dataset is only
678     necessary when ‘--subtractsky’ is called or when a column directly
679     related to the Sky value is requested (currently ‘--sky’).  This
680     dataset may be a tessellation, with one element per tile (see
681     ‘--oneelempertile’ of NoiseChisel’s *note Processing options::).
682
683     When the Sky dataset is necessary but this option is not called,
684     MakeCatalog will assume it is an HDU/extension (specified by
685     ‘--skyhdu’) in one of the already given files.  First it will look
686     for it in the ‘--valuesfile’ (if it is given) and then the main
687     input file (given as an argument).
688
689     By default the values dataset is assumed to be already Sky
690     subtracted, so this dataset is not necessary for many of the
691     columns.
692
693‘--skyhdu=STR’
694     HDU/extension of the Sky dataset, see ‘--skyfile’.
695
696‘--subtractsky’
697     Subtract the sky value or dataset from the values file prior to any
698     processing.
699
700‘-t STR/FLT701‘--instd=STR/FLT702     Sky standard deviation value as a single number, or the file name
703     containing a dataset (different values per pixel or tile).  With
704     the ‘--variance’ option you can tell MakeCatalog to interpret this
705     value/dataset as a variance image, not standard deviation.
706
707     *Important note:* This must only be the SKY standard deviation or
708     variance (not including the signal’s contribution to the error).
709     In other words, the final standard deviation of a pixel depends on
710     how much signal there is in it.  MakeCatalog will find the amount
711     of signal within each pixel (while subtracting the Sky, if
712     ‘--subtractsky’ is called) and account for the extra error due to
713     it’s value (signal).  Therefore if the input standard deviation (or
714     variance) image also contains the contribution of signal to the
715     error, then the final error measurements will be over-estimated.
716
717‘--stdhdu=STR’
718     The HDU of the Sky value standard deviation image.
719
720‘--variance’
721     The dataset given to ‘--instd’ (and ‘--stdhdu’ has the Sky variance
722     of every pixel, not the Sky standard deviation.
723
724‘--forcereadstd’
725     Read the input STD image even if it is not required by any of the
726     requested columns.  This is because some of the output catalog’s
727     metadata may need it, for example to calculate the dataset’s
728     surface brightness limit (see *note Quantifying measurement
729     limits::, configured with ‘--sfmagarea’ and ‘--sfmagnsigma’ in
730     *note MakeCatalog output::).
731
732     Furthermore, if the input STD image doesn’t have the ‘MEDSTD’
733     keyword (that is meant to contain the representative standard
734     deviation of the full image), with this option, the median will be
735     calculated and used for the surface brightness limit.
736
737‘-z FLT’
738‘--zeropoint=FLT’
739     The zero point magnitude for the input image, see *note Brightness
740     flux magnitude::.
741
742‘--sigmaclip FLT,FLT’
743     The sigma-clipping parameters when any of the sigma-clipping
744     related columns are requested (for example ‘--sigclip-median’ or
745     ‘--sigclip-number’).
746
747     This option takes two values: the first is the multiple of
748     $\sigma$, and the second is the termination criteria.  If the
749     latter is larger than 1, it is read as an integer number and will
750     be the number of times to clip.  If it is smaller than 1, it is
751     interpreted as the tolerance level to stop clipping.  See *note
752     Sigma clipping:: for a complete explanation.
753
754‘--fracmax=FLT[,FLT]’
755     The fractions (one or two) of maximum value in objects or clumps to
756     be used in the related columns, for example ‘--fracmaxarea1’,
757     ‘--fracmaxsum1’ or ‘--fracmaxradius1’, see *note MakeCatalog
758     measurements::.  For the maximum value, see the description of
759     ‘--maximum’ column below.  The value(s) of this option must be
760     larger than 0 and smaller than 1 (they are a fraction).  When only
761     ‘--fracmaxarea1’ or ‘--fracmaxsum1’ is requested, one value must be
762     given to this option, but if ‘--fracmaxarea2’ or ‘--fracmaxsum2’
763     are also requested, two values must be given to this option.  The
764     values can be written as simple floating point numbers, or as
765     fractions, for example ‘0.25,0.75’ and ‘0.25,3/4’ are the same.
766
767‘--spatialresolution=FLT’
768     The error in measuring spatial properties (for example the area) in
769     units of pixels.  You can think of this as the FWHM of the
770     dataset’s PSF and is used in measurements like the error in surface
771     brightness (‘--sberror’, see *note MakeCatalog measurements::).
772     Ideally, images are taken in the optimal Nyquist sampling *note
773     Sampling theorem::, so the default value for this option is 2.  But
774     in practice real images my be over-sampled (usually ground-based
775     images, where you will need to increase the default value) or
776     undersampled (some space-based images, where you will need to
777     decrease the default value).
778
779
780File: gnuastro.info,  Node: Upper-limit settings,  Next: MakeCatalog measurements,  Prev: MakeCatalog inputs and basic settings,  Up: Invoking astmkcatalog
781
7827.4.6.2 Upper-limit settings
783............................
784
785The upper-limit magnitude was discussed in *note Quantifying measurement
786limits::.  Unlike other measured values/columns in MakeCatalog, the
787upper limit magnitude needs several extra parameters which are discussed
788here.  All the options specific to the upper-limit measurements start
789with ‘up’ for “upper-limit”.  The only exception is ‘--envseed’ that is
790also present in other programs and is general for any job requiring
791random number generation in Gnuastro (see *note Generating random
792numbers::).
793
794   One very important consideration in Gnuastro is reproducibility.
795Therefore, the values to all of these parameters along with others (like
796the random number generator type and seed) are also reported in the
797comments of the final catalog when the upper limit magnitude column is
798desired.  The random seed that is used to define the random positions
799for each object or clump is unique and set based on the (optionally)
800given seed, the total number of objects and clumps and also the labels
801of the clumps and objects.  So with identical inputs, an identical
802upper-limit magnitude will be found.  However, even if the seed is
803identical, when the ordering of the object/clump labels differs between
804different runs, the result of upper-limit measurements will not be
805identical.
806
807   MakeCatalog will randomly place the object/clump footprint over the
808dataset.  When the randomly placed footprint doesn’t fall on any object
809or masked region (see ‘--upmaskfile’) it will be used in the final
810distribution.  Otherwise that particular random position will be ignored
811and another random position will be generated.  Finally, when the
812distribution has the desired number of successfully measured random
813samples (‘--upnum’) the distribution’s properties will be measured and
814placed in the catalog.
815
816   When the profile is very large or the image is significantly covered
817by detections, it might not be possible to find the desired number of
818samplings in a reasonable time.  MakeProfiles will continue searching
819until it is unable to find a successful position (since the last
820successful measurement(1)), for a large multiple of ‘--upnum’
821(currently(2) this is 10).  If ‘--upnum’ successful samples cannot be
822found until this limit is reached, MakeCatalog will set the upper-limit
823magnitude for that object to NaN (blank).
824
825   MakeCatalog will also print a warning if the range of positions
826available for the labeled region is smaller than double the size of the
827region.  In such cases, the limited range of random positions can
828artificially decrease the standard deviation of the final distribution.
829If your dataset can allow it (it is large enough), it is recommended to
830use a larger range if you see such warnings.
831
832‘--upmaskfile=FITS’
833     File name of mask image to use for upper-limit calculation.  In
834     some cases (especially when doing matched photometry), the object
835     labels specified in the main input and mask image might not be
836     adequate.  In other words they do not necessarily have to cover
837     _all_ detected objects: the user might have selected only a few of
838     the objects in their labeled image.  This option can be used to
839     ignore regions in the image in these situations when estimating the
840     upper-limit magnitude.  All the non-zero pixels of the image
841     specified by this option (in the ‘--upmaskhdu’ extension) will be
842     ignored in the upper-limit magnitude measurements.
843
844     For example, when you are using labels from another image, you can
845     give NoiseChisel’s objects image output for this image as the value
846     to this option.  In this way, you can be sure that regions with
847     data do not harm your distribution.  See *note Quantifying
848     measurement limits:: for more on the upper limit magnitude.
849
850‘--upmaskhdu=STR’
851     The extension in the file specified by ‘--upmask’.
852
853‘--upnum=INT’
854     The number of random samples to take for all the objects.  A larger
855     value to this option will give a more accurate result
856     (asymptotically), but it will also slow down the process.  When a
857     randomly positioned sample overlaps with a detected/masked pixel it
858     is not counted and another random position is found until the
859     object completely lies over an undetected region.  So you can be
860     sure that for each object, this many samples over undetected
861     objects are made.  See the upper limit magnitude discussion in
862     *note Quantifying measurement limits:: for more.
863
864‘--uprange=INT,INT’
865     The range/width of the region (in pixels) to do random sampling
866     along each dimension of the input image around each object’s
867     position.  This is not a mandatory option and if not given (or
868     given a value of zero in a dimension), the full possible range of
869     the dataset along that dimension will be used.  This is useful when
870     the noise properties of the dataset vary gradually.  In such cases,
871     using the full range of the input dataset is going to bias the
872     result.  However, note that decreasing the range of available
873     positions too much will also artificially decrease the standard
874     deviation of the final distribution (and thus bias the upper-limit
875     measurement).
876
877‘--envseed’
878     Read the random number generator type and seed value from the
879     environment (see *note Generating random numbers::).  Random
880     numbers are used in calculating the random positions of different
881     samples of each object.
882
883‘--upsigmaclip=FLT,FLT’
884     The raw distribution of random values will not be used to find the
885     upper-limit magnitude, it will first be $\sigma$-clipped (see *note
886     Sigma clipping::) to avoid outliers in the distribution (mainly the
887     faint undetected wings of bright/large objects in the image).  This
888     option takes two values: the first is the multiple of $\sigma$, and
889     the second is the termination criteria.  If the latter is larger
890     than 1, it is read as an integer number and will be the number of
891     times to clip.  If it is smaller than 1, it is interpreted as the
892     tolerance level to stop clipping.  See *note Sigma clipping:: for a
893     complete explanation.
894
895‘--upnsigma=FLT’
896     The multiple of the final ($\sigma$-clipped) standard deviation (or
897     $\sigma$) used to measure the upper-limit brightness or magnitude.
898
899‘--checkuplim=INT[,INT]’
900     Print a table of positions and measured values for all the full
901     random distribution used for one particular object or clump.  If
902     only one integer is given to this option, it is interpreted to be
903     an object’s label.  If two values are given, the first is the
904     object label and the second is the ID of requested clump within it.
905
906     The output is a table with three columns (its type is determined
907     with the ‘--tableformat’ option, see *note Input output options::).
908     The first two columns are the position of the first pixel in each
909     random sampling of this particular object/clump.  The the third
910     column is the measured flux over that region.  If the region
911     overlapped with a detection or masked pixel, then its measured
912     value will be a NaN (not-a-number).  The total number of rows is
913     thus unknown, but you can be sure that the number of rows with
914     non-NaN measurements is the number given to the ‘--upnum’ option.
915
916   ---------- Footnotes ----------
917
918   (1) The counting of failed positions restarts on every successful
919measurement.
920
921   (2) In Gnuastro’s source, this constant number is defined as the
922‘MKCATALOG_UPPERLIMIT_MAXFAILS_MULTIP’ macro in ‘bin/mkcatalog/main.h’,
923see *note Downloading the source::.
924
925
926File: gnuastro.info,  Node: MakeCatalog measurements,  Next: MakeCatalog output,  Prev: Upper-limit settings,  Up: Invoking astmkcatalog
927
9287.4.6.3 MakeCatalog measurements
929................................
930
931The final group of options particular to MakeCatalog are those that
932specify which measurements/columns should be written into the final
933output table.  The current measurements in MakeCatalog are those which
934only produce one final value for each label (for example its total
935brightness: a single number).  All the different label’s measurements
936can be written as one column in a final table/catalog that contains
937other columns for other similar single-number measurements.
938
939   In this case, all the different label’s measurements can be written
940as one column in a final table/catalog that contains other columns for
941other similar single-number measurements.  The majority of this section
942is devoted to MakeCatalog’s single-valued measurements.  However,
943MakeCatalog can also do measurements that produce more than one value
944for each label.  Currently the only such measurement is generation of
945spectra from 3D cubes with the ‘--spectrum’ option and it is discussed
946in the end of this section.
947
948   Command-line options are used to identify which measurements you want
949in the final catalog(s) and in what order.  If any of the options below
950is called on the command line or in any of the configuration files, it
951will be included as a column in the output catalog.  The order of the
952columns is in the same order as the options were seen by MakeCatalog
953(see *note Configuration file precedence::).  Some of the columns apply
954to both “objects” and “clumps” and some are particular to only one of
955them (for the definition of “objects” and “clumps”, see *note
956Segment::).  Columns/options that are unique to one catalog (only
957objects, or only clumps), are explicitly marked with [Objects] or
958[Clumps] to specify the catalog they will be placed in.
959
960‘--i’
961‘--ids’
962     This is a unique option which can add multiple columns to the final
963     catalog(s).  Calling this option will put the object IDs
964     (‘--objid’) in the objects catalog and host-object-ID
965     (‘--hostobjid’) and ID-in-host-object (‘--idinhostobj’) into the
966     clumps catalog.  Hence if only object catalogs are required, it has
967     the same effect as ‘--objid’.
968
969‘--objid’
970     [Objects] ID of this object.
971
972‘-j’
973‘--hostobjid’
974     [Clumps] The ID of the object which hosts this clump.
975
976‘--idinhostobj’
977     [Clumps] The ID of this clump in its host object.
978
979‘-x’
980‘--x’
981     The flux weighted center of all objects and clumps along the first
982     FITS axis (horizontal when viewed in SAO DS9), see $\overline{x}$
983     in *note Measuring elliptical parameters::.  The weight has to have
984     a positive value (pixel value larger than the Sky value) to be
985     meaningful!  Specially when doing matched photometry, this might
986     not happen: no pixel value might be above the Sky value.  For such
987     detections, the geometric center will be reported in this column
988     (see ‘--geox’).  You can use ‘--weightarea’ to see which was used.
989
990‘-y’
991‘--y’
992     The flux weighted center of all objects and clumps along the second
993     FITS axis (vertical when viewed in SAO DS9).  See ‘--x’.
994
995‘-z’
996‘--z’
997     The flux weighted center of all objects and clumps along the third
998     FITS axis.  See ‘--x’.
999
1000‘--geox’
1001     The geometric center of all objects and clumps along the first FITS
1002     axis axis.  The geometric center is the average pixel positions
1003     irrespective of their pixel values.
1004
1005‘--geoy’
1006     The geometric center of all objects and clumps along the second
1007     FITS axis axis, see ‘--geox’.
1008
1009‘--geoz’
1010     The geometric center of all objects and clumps along the third FITS
1011     axis axis, see ‘--geox’.
1012
1013‘--minvx’
1014     Position of pixel with minimum value in objects and clumps, along
1015     the first FITS axis.
1016
1017‘--maxvx’
1018     Position of pixel with minimum value in objects and clumps, along
1019     the first FITS axis.
1020
1021‘--minvy’
1022     Position of pixel with minimum value in objects and clumps, along
1023     the first FITS axis.
1024
1025‘--maxvy’
1026     Position of pixel with minimum value in objects and clumps, along
1027     the first FITS axis.
1028
1029‘--minvz’
1030     Position of pixel with minimum value in objects and clumps, along
1031     the first FITS axis.
1032
1033‘--maxvz’
1034     Position of pixel with minimum value in objects and clumps, along
1035     the first FITS axis.
1036
1037‘--minx’
1038     The minimum position of all objects and clumps along the first FITS
1039     axis.
1040
1041‘--maxx’
1042     The maximum position of all objects and clumps along the first FITS
1043     axis.
1044
1045‘--miny’
1046     The minimum position of all objects and clumps along the second
1047     FITS axis.
1048
1049‘--maxy’
1050     The maximum position of all objects and clumps along the second
1051     FITS axis.
1052
1053‘--minz’
1054     The minimum position of all objects and clumps along the third FITS
1055     axis.
1056
1057‘--maxz’
1058     The maximum position of all objects and clumps along the third FITS
1059     axis.
1060
1061‘--clumpsx’
1062     [Objects] The flux weighted center of all the clumps in this object
1063     along the first FITS axis.  See ‘--x’.
1064
1065‘--clumpsy’
1066     [Objects] The flux weighted center of all the clumps in this object
1067     along the second FITS axis.  See ‘--x’.
1068
1069‘--clumpsz’
1070     [Objects] The flux weighted center of all the clumps in this object
1071     along the third FITS axis.  See ‘--x’.
1072
1073‘--clumpsgeox’
1074     [Objects] The geometric center of all the clumps in this object
1075     along the first FITS axis.  See ‘--geox’.
1076
1077‘--clumpsgeoy’
1078     [Objects] The geometric center of all the clumps in this object
1079     along the second FITS axis.  See ‘--geox’.
1080
1081‘--clumpsgeoz’
1082     [Objects] The geometric center of all the clumps in this object
1083     along the third FITS axis.  See ‘--geoz’.
1084
1085‘-r’
1086‘--ra’
1087     Flux weighted right ascension of all objects or clumps, see ‘--x’.
1088     This is just an alias for one of the lower-level ‘--w1’ or ‘--w2’
1089     options.  Using the FITS WCS keywords (‘CTYPE’), MakeCatalog will
1090     determine which axis corresponds to the right ascension.  If no
1091     ‘CTYPE’ keywords start with ‘RA’, an error will be printed when
1092     requesting this column and MakeCatalog will abort.
1093
1094‘-d’
1095‘--dec’
1096     Flux weighted declination of all objects or clumps, see ‘--x’.
1097     This is just an alias for one of the lower-level ‘--w1’ or ‘--w2’
1098     options.  Using the FITS WCS keywords (‘CTYPE’), MakeCatalog will
1099     determine which axis corresponds to the declination.  If no ‘CTYPE’
1100     keywords start with ‘DEC’, an error will be printed when requesting
1101     this column and MakeCatalog will abort.
1102
1103‘--w1’
1104     Flux weighted first WCS axis of all objects or clumps, see ‘--x’.
1105     The first WCS axis is commonly used as right ascension in images.
1106
1107‘--w2’
1108     Flux weighted second WCS axis of all objects or clumps, see ‘--x’.
1109     The second WCS axis is commonly used as declination in images.
1110
1111‘--w3’
1112     Flux weighted third WCS axis of all objects or clumps, see ‘--x’.
1113     The third WCS axis is commonly used as wavelength in integral field
1114     unit data cubes.
1115
1116‘--geow1’
1117     Geometric center in first WCS axis of all objects or clumps, see
1118     ‘--geox’.  The first WCS axis is commonly used as right ascension
1119     in images.
1120
1121‘--geow2’
1122     Geometric center in second WCS axis of all objects or clumps, see
1123     ‘--geox’.  The second WCS axis is commonly used as declination in
1124     images.
1125
1126‘--geow3’
1127     Geometric center in third WCS axis of all objects or clumps, see
1128     ‘--geox’.  The third WCS axis is commonly used as wavelength in
1129     integral field unit data cubes.
1130
1131‘--clumpsw1’
1132     [Objects] Flux weighted center in first WCS axis of all clumps in
1133     this object, see ‘--x’.  The first WCS axis is commonly used as
1134     right ascension in images.
1135
1136‘--clumpsw2’
1137     [Objects] Flux weighted declination of all clumps in this object,
1138     see ‘--x’.  The second WCS axis is commonly used as declination in
1139     images.
1140
1141‘--clumpsw3’
1142     [Objects] Flux weighted center in third WCS axis of all clumps in
1143     this object, see ‘--x’.  The third WCS axis is commonly used as
1144     wavelength in integral field unit data cubes.
1145
1146‘--clumpsgeow1’
1147     [Objects] Geometric center right ascension of all clumps in this
1148     object, see ‘--geox’.  The first WCS axis is commonly used as right
1149     ascension in images.
1150
1151‘--clumpsgeow2’
1152     [Objects] Geometric center declination of all clumps in this
1153     object, see ‘--geox’.  The second WCS axis is commonly used as
1154     declination in images.
1155
1156‘--clumpsgeow3’
1157     [Objects] Geometric center in third WCS axis of all clumps in this
1158     object, see ‘--geox’.  The third WCS axis is commonly used as
1159     wavelength in integral field unit data cubes.
1160
1161‘-b’
1162‘--brightness’
1163     The brightness (sum of all pixel values), see *note Brightness flux
1164     magnitude::.  For clumps, the ambient brightness (flux of river
1165     pixels around the clump multiplied by the area of the clump) is
1166     removed, see ‘--riverave’.  So the sum of all the clumps brightness
1167     in the clump catalog will be smaller than the total clump
1168     brightness in the ‘--clumpbrightness’ column of the objects
1169     catalog.
1170
1171     If no usable pixels are present over the clump or object (for
1172     example they are all blank), the returned value will be NaN (note
1173     that zero is meaningful).
1174
1175‘--brightnesserr’
1176     The ($1\sigma$) error in measuring the brightness of a label
1177     (objects or clumps).
1178
1179     The returned value will be NaN when the label covers only NaN
1180     pixels in the values image, or a pixel is NaN in the ‘--instd’
1181     image, but non-NaN in the values image.  The latter situation
1182     usually happens when there is a bug in the previous steps of your
1183     analysis, and is important because those pixels with a NaN in the
1184     ‘--instd’ image may contribute significantly to the final error.
1185     If you want to ignore those pixels in the error measurement, set
1186     them to zero (which is a meaningful number in such scenarios).
1187
1188‘--clumpbrightness’
1189     [Objects] The total brightness of the clumps within an object.
1190     This is simply the sum of the pixels associated with clumps in the
1191     object.  If no usable pixels are present over the clump or object
1192     (for example they are all blank), the stored value will be NaN
1193     (note that zero is meaningful).
1194
1195‘--brightnessnoriver’
1196     [Clumps] The Sky (not river) subtracted clump brightness.  By
1197     definition, for the clumps, the average brightness of the rivers
1198     surrounding it are subtracted from it for a first order accounting
1199     for contamination by neighbors.  In cases where you will be
1200     calculating the flux brightness difference later (one example
1201     below) the contamination will be (mostly) removed at that stage,
1202     which is why this column was added.
1203
1204     If no usable pixels are present over the clump or object (for
1205     example they are all blank), the stored value will be NaN (note
1206     that zero is meaningful).
1207
1208‘--mean’
1209     The mean sky subtracted value of pixels within the object or clump.
1210     For clumps, the average river flux is subtracted from the sky
1211     subtracted mean.
1212
1213‘--median’
1214     The median sky subtracted value of pixels within the object or
1215     clump.  For clumps, the average river flux is subtracted from the
1216     sky subtracted median.
1217
1218‘--maximum’
1219     The maximum value of pixels within the object or clump.  When the
1220     label (object or clump) is larger than three pixels, the maximum is
1221     actually derived by the mean of the brightest three pixels, not the
1222     largest pixel value of the same label.  This is because noise
1223     fluctuations can be very strong in the extreme values of the
1224     objects/clumps due to Poisson noise (which gets stronger as the
1225     mean gets higher).  Simply using the maximum pixel value will
1226     create a strong scatter in results that depend on the maximum (for
1227     example the ‘--fwhm’ option also uses this value internally).
1228
1229‘--sigclip-number’
1230     The number of elements/pixels in the dataset after sigma-clipping
1231     the object or clump.  The sigma-clipping parameters can be set with
1232     the ‘--sigmaclip’ option described in *note MakeCatalog inputs and
1233     basic settings::.  For more on Sigma-clipping, see *note Sigma
1234     clipping::.
1235
1236‘--sigclip-median’
1237     The sigma-clipped median value of the object of clump’s pixel
1238     distribution.  For more on sigma-clipping and how to define it, see
1239     ‘--sigclip-number’.
1240
1241‘--sigclip-mean’
1242     The sigma-clipped mean value of the object of clump’s pixel
1243     distribution.  For more on sigma-clipping and how to define it, see
1244     ‘--sigclip-number’.
1245
1246‘--sigclip-std’
1247     The sigma-clipped standard deviation of the object of clump’s pixel
1248     distribution.  For more on sigma-clipping and how to define it, see
1249     ‘--sigclip-number’.
1250
1251‘-m’
1252‘--magnitude’
1253     The magnitude of clumps or objects, see ‘--brightness’.
1254
1255‘-e’
1256‘--magnitudeerr’
1257     The magnitude error of clumps or objects.  The magnitude error is
1258     calculated from the signal-to-noise ratio (see ‘--sn’ and *note
1259     Quantifying measurement limits::).  Note that until now this error
1260     assumes uncorrelated pixel values and also does not include the
1261     error in estimating the aperture (or error in generating the
1262     labeled image).
1263
1264     For now these factors have to be found by other means.  Task 14124
1265     (https://savannah.gnu.org/task/index.php?14124) has been defined
1266     for work on adding these sources of error too.
1267
1268     The returned value will be NaN when the label covers only NaN
1269     pixels in the values image, or a pixel is NaN in the ‘--instd’
1270     image, but non-NaN in the values image.  The latter situation
1271     usually happens when there is a bug in the previous steps of your
1272     analysis, and is important because those pixels with a NaN in the
1273     ‘--instd’ image may contribute significantly to the final error.
1274     If you want to ignore those pixels in the error measurement, set
1275     them to zero (which is a meaningful number in such scenarios).
1276
1277‘--clumpsmagnitude’
1278     [Objects] The magnitude of all clumps in this object, see
1279     ‘--clumpbrightness’.
1280
1281‘--upperlimit’
1282     The upper limit value (in units of the input image) for this object
1283     or clump.  This is the sigma-clipped standard deviation of the
1284     random distribution, multiplied by the value of ‘--upnsigma’).  See
1285     *note Quantifying measurement limits:: and *note Upper-limit
1286     settings:: for a complete explanation.  This is very important for
1287     the fainter and smaller objects in the image where the measured
1288     magnitudes are not reliable.
1289
1290‘--upperlimitmag’
1291     The upper limit magnitude for this object or clump.  See *note
1292     Quantifying measurement limits:: and *note Upper-limit settings::
1293     for a complete explanation.  This is very important for the fainter
1294     and smaller objects in the image where the measured magnitudes are
1295     not reliable.
1296
1297‘--upperlimitsb’
1298     The upper-limit surface brightness (in units of mag/arcsec$^2$) of
1299     this labeled region (object or clump).  This is just a simple
1300     wrapper over lower-level columns: setting B and A as the value in
1301     the columns ‘--upperlimit’ and ‘--areaarcsec2’, we fill this column
1302     by simply use the surface brightness equation of *note Brightness
1303     flux magnitude::.
1304
1305‘--upperlimitonesigma’
1306     The $1\sigma$ upper limit value (in units of the input image) for
1307     this object or clump.  See *note Quantifying measurement limits::
1308     and *note Upper-limit settings:: for a complete explanation.  When
1309     ‘--upnsigma=1’, this column’s values will be the same as
1310     ‘--upperlimit’.
1311
1312‘--upperlimitsigma’
1313     The position of the total brightness measured within the
1314     distribution of randomly placed upperlimit measurements in units of
1315     the distribution’s $\sigma$ or standard deviation.  See *note
1316     Quantifying measurement limits:: and *note Upper-limit settings::
1317     for a complete explanation.
1318
1319‘--upperlimitquantile’
1320     The position of the total brightness measured within the
1321     distribution of randomly placed upperlimit measurements as a
1322     quantile (value between 0 or 1).  See *note Quantifying measurement
1323     limits:: and *note Upper-limit settings:: for a complete
1324     explanation.  If the object is brighter than the brightest randomly
1325     placed profile, a value of ‘inf’ is returned.  If it is less than
1326     the minimum, a value of ‘-inf’ is reported.
1327
1328‘--upperlimitskew’
1329     This column contains the non-parametric skew of the
1330     $\sigma$-clipped random distribution that was used to estimate the
1331     upper-limit magnitude.  Taking $\mu$ as the mean, $\nu$ as the
1332     median and $\sigma$ as the standard deviation, the traditional
1333     definition of skewness is defined as: $(\mu-\nu)/\sigma$.
1334
1335     This can be a good measure to see how much you can trust the random
1336     measurements, or in other words, how accurately the regions with
1337     signal have been masked/detected.  If the skewness is strong (and
1338     to the positive), then you can tell that you have a lot of
1339     undetected signal in the dataset, and therefore that the
1340     upper-limit measurement (and other measurements) are not reliable.
1341
1342‘--riverave’
1343     [Clumps] The average brightness of the river pixels around this
1344     clump.  River pixels were defined in Akhlaghi and Ichikawa 2015.
1345     In short they are the pixels immediately outside of the clumps.
1346     This value is used internally to find the brightness (or magnitude)
1347     and signal to noise ratio of the clumps.  It can generally also be
1348     used as a scale to gauge the base (ambient) flux surrounding the
1349     clump.  In case there was no river pixels, then this column will
1350     have the value of the Sky under the clump.  So note that this value
1351     is _not_ sky subtracted.
1352
1353‘--rivernum’
1354     [Clumps] The number of river pixels around this clump, see
1355     ‘--riverave’.
1356
1357‘-n’
1358‘--sn’
1359     The Signal to noise ratio (S/N) of all clumps or objects.  See
1360     Akhlaghi and Ichikawa (2015) for the exact equations used.
1361
1362     The returned value will be NaN when the label covers only NaN
1363     pixels in the values image, or a pixel is NaN in the ‘--instd’
1364     image, but non-NaN in the values image.  The latter situation
1365     usually happens when there is a bug in the previous steps of your
1366     analysis, and is important because those pixels with a NaN in the
1367     ‘--instd’ image may contribute significantly to the final error.
1368     If you want to ignore those pixels in the error measurement, set
1369     them to zero (which is a meaningful number).
1370
1371‘--sky’
1372     The sky flux (per pixel) value under this object or clump.  This is
1373     actually the mean value of all the pixels in the sky image that lie
1374     on the same position as the object or clump.
1375
1376‘--std’
1377     The sky value standard deviation (per pixel) for this clump or
1378     object.  This is the square root of the mean variance under the
1379     object, or the root mean square.
1380
1381‘-C’
1382‘--numclumps’
1383     [Objects] The number of clumps in this object.
1384
1385‘-a’
1386‘--area’
1387     The raw area (number of pixels/voxels) in any clump or object
1388     independent of what pixel it lies over (if it is NaN/blank or
1389     unused for example).
1390
1391‘--areaarcsec2’
1392     The used (non-blank in values image) area of the labeled region in
1393     units of arc-seconds squared.  This column is just the value of the
1394     ‘--area’ column, multiplied by the area of each pixel in the input
1395     image (in units of arcsec^2).  Similar to the ‘--ra’ or ‘--dec’
1396     columns, for this option to work, the objects extension used has to
1397     have a WCS structure.
1398
1399‘--areaminv’
1400     The number of pixels that are equal to the minimum value of the
1401     labeled region (clump or object).
1402
1403‘--areamaxv’
1404     The number of pixels that are equal to the maximum value of the
1405     labeled region (clump or object).
1406
1407‘--surfacebrightness’
1408     The surface brightness (in units of mag/arcsec$^2$) of the labeled
1409     region (objects or clumps).  For more on the definition of the
1410     surface brightness, see *note Brightness flux magnitude::.
1411
1412‘--sberror’
1413     Error in measuring the surface brightness (the
1414     ‘--surfacebrightness’ column).  This column will use the value
1415     given to ‘--spatialresolution’ in the processing (in pixels).  For
1416     more on ‘--spatialresolution’, see see *note MakeCatalog inputs and
1417     basic settings:: and for the equation used to derive the surface
1418     brightness error, see *note Surface brightness error of each
1419     detection::.
1420
1421‘--areaxy’
1422     Similar to ‘--area’, when the clump or object is projected onto the
1423     first two dimensions.  This is only available for 3-dimensional
1424     datasets.  When working with Integral Field Unit (IFU) datasets,
1425     this projection onto the first two dimensions would be a
1426     narrow-band image.
1427
1428‘--fwhm’
1429     The full width at half maximum (in units of pixels, along the
1430     semi-major axis) of the labeled region (object or clump).  The
1431     maximum value is estimated from the mean of the top-three pixels
1432     with the highest values, see the description under ‘--maximum’.
1433     The number of pixels that have half the value of that maximum are
1434     then found (value in the ‘--halfmaxarea’ column) and a radius is
1435     estimated from the area.  See the description under
1436     ‘--halfsumradius’ for more on converting area to radius along major
1437     axis.
1438
1439     Because of its non-parametric nature, this column is most reliable
1440     on clumps and should only be used in objects with great caution.
1441     This is because objects can have more than one clump (peak with
1442     true signal) and multiple peaks are not treated separately in
1443     objects, so the result of this column will be biased.
1444
1445     Also, because of its non-parametric nature, this FWHM it does not
1446     account for the PSF, and it will be strongly affected by noise if
1447     the object is faint.  So when half the maximum value (which can be
1448     requested using the ‘--maximum’ column) is too close to the local
1449     noise level (which can be requested using the ‘--std’ column), the
1450     value returned in this column is meaningless (its just noise peaks
1451     which are randomly distributed over the area).  You can therefore
1452     use the ‘--maximum’ and ‘--std’ columns to remove unreliable FWHMs.
1453     For example if a labeled region’s maximum is less than 2 times the
1454     sky standard deviation, the value will certainly be unreliable
1455     (half of that is $1\sigma$!).  For a more reliable value, this
1456     fraction should be around 4 (so half the maximum is 2$\sigma$).
1457
1458‘--halfmaxarea’
1459     The number of pixels with values larger than half the maximum flux
1460     within the labeled region.  This option is used to estimate
1461     ‘--fwhm’, so please read the notes there for the caveats and
1462     necessary precautions.
1463
1464‘--halfmaxradius’
1465     The radius of region containing half the maximum flux within the
1466     labeled region.  This is just half the value reported by ‘--fwhm’.
1467
1468‘--halfmaxsum’
1469     The sum of the pixel values containing half the maximum flux within
1470     the labeled region (or those that are counted in ‘--halfmaxarea’).
1471     This option uses the pixels within ‘--fwhm’, so please read the
1472     notes there for the caveats and necessary precautions.
1473
1474‘--halfmaxsb’
1475     The surface brightness (in units of mag/arcsec$^2$) within the
1476     region that contains half the maximum value of the labeled region.
1477     For more on the definition of the surface brightness, see *note
1478     Brightness flux magnitude::.
1479
1480‘--halfsumarea’
1481     The number of pixels that contain half the object or clump’s total
1482     sum of pixels (half the value in the ‘--brightness’ column).  To
1483     count this area, all the non-blank values associated with the given
1484     label (object or clump) will be sorted and summed in order
1485     (starting from the maximum), until the sum becomes larger than half
1486     the total sum of the label’s pixels.
1487
1488     This option is thus good for clumps (which are defined to have a
1489     single peak in their morphology), but for objects you should be
1490     careful: if the object includes multiple peaks/clumps at roughly
1491     the same level, then the area reported by this option will be
1492     distributed over all the peaks.
1493
1494‘--halfsumsb’
1495     Surface brightness (in units of mag/arcsec$^2$) within the area
1496     that contains half the total sum of the label’s pixels (object or
1497     clump).  For more on the definition of the surface brightness, see
1498     *note Brightness flux magnitude::.
1499
1500     This column just plugs in the values of half the value of the
1501     ‘--brightness’ column and the ‘--halfsumarea’ column, into the
1502     surface brightness equation.  Therefore please see the description
1503     in ‘--halfsumarea’ to understand the systematics of this column and
1504     potential biases.
1505
1506‘--halfsumradius’
1507     Radius (in units of pixels) derived from the area that contains
1508     half the total sum of the label’s pixels (value reported by
1509     ‘--halfsumarea’).  If the area is $A_h$ and the axis ratio is $q$,
1510     then the value returned in this column is $\sqrt{A_h/({\pi}q)}$.
1511     This option is a good measure of the concentration of the
1512     _observed_ (after PSF convolution and noisy) object or clump, But
1513     as described below it underestimates the effective radius.  Also,
1514     it should be used in caution with objects that may have multiple
1515     clumps.  It is most reliable with clumps or objects that have one
1516     or zero clumps, see the note under ‘--halfsumarea’.
1517
1518     Recall that in general, for an ellipse with semi-major axis $a$,
1519     semi-minor axis $b$, and axis ratio $q=b/a$ the area ($A$) is
1520     $A={\pi}ab={\pi}qa^2$.  For a circle (where $q=1$), this simplifies
1521     to the familiar $A={\pi}a^2$.
1522
1523     This option should not be confused with the _effective radius_ for
1524     Sérsic profiles, commonly written as $r_e$.  For more on the Sérsic
1525     profile and $r_e$, please see *note Galaxies::.  Therefore, when
1526     $r_e$ is meaningful for the target (the target is elliptically
1527     symmetric and can be parameterized as a Sérsic profile), $r_e$
1528     should be derived from fitting the profile with a Sérsic function
1529     which has been convolved with the PSF. But from the equation above,
1530     you see that this radius is derived from the raw image’s labeled
1531     values (after convolution, with no parametric profile), so this
1532     column’s value will generally be (much) smaller than $r_e$,
1533     depending on the PSF, depth of the dataset, the morphology, or if a
1534     fraction of the profile falls on the edge of the image.
1535
1536     In other words, this option can only be interpreted as an effective
1537     radius if there is no noise and no PSF and the profile within the
1538     image extends to infinity (or a very large multiple of the
1539     effective radius) and it not near the edge of the image.
1540
1541‘--fracmaxarea1’
1542‘--fracmaxarea2’
1543     Number of pixels brighter than the given fraction(s) of the maximum
1544     pixel value.  For the maximum value, see the description of
1545     ‘--maximum’ column.  The fraction(s) are given through the
1546     ‘--fracmax’ option (that can take two values) and is described in
1547     *note MakeCatalog inputs and basic settings::.  Recall that in
1548     ‘--halfmaxarea’, the fraction is fixed to 0.5.  Hence, added with
1549     these two columns, you can sample three parts of the profile area.
1550
1551‘--fracmaxsum1’
1552‘--fracmaxsum2’
1553     Sum of pixels brighter than the given fraction(s) of the maximum
1554     pixel value.  For the maximum value, see the description of
1555     ‘--maximum’ column below.  The fraction(s) are given through the
1556     ‘--fracmax’ option (that can take two values) and is described in
1557     *note MakeCatalog inputs and basic settings::.  Recall that in
1558     ‘--halfmaxsum’, the fraction is fixed to 0.5.  Hence, added with
1559     these two columns, you can sample three parts of the profile’s sum
1560     of pixels.
1561
1562‘--fracmaxradius1’
1563‘--fracmaxradius2’
1564     Radius (in units of pixels) derived from the area that contains the
1565     given fractions of the maximum valued pixel(s) of the label’s
1566     pixels (value reported by ‘--fracmaxarea1’ or ‘--fracmaxarea2’).
1567     For the maximum value, see the description of ‘--maximum’ column
1568     below.  The fractions are given through the ‘--fracmax’ option
1569     (that can take two values) and is described in *note MakeCatalog
1570     inputs and basic settings::.  Recall that in ‘--fwhm’, the fraction
1571     is fixed to 0.5.  Hence, added with these two columns, you can
1572     sample three parts of the profile’s radius.
1573
1574‘--clumpsarea’
1575     [Objects] The total area of all the clumps in this object.
1576
1577‘--weightarea’
1578     The area (number of pixels) used in the flux weighted position
1579     calculations.
1580
1581‘--geoarea’
1582     The area of all the pixels labeled with an object or clump.  Note
1583     that unlike ‘--area’, pixel values are completely ignored in this
1584     column.  For example, if a pixel value is blank, it won’t be
1585     counted in ‘--area’, but will be counted here.
1586
1587‘--geoareaxy’
1588     Similar to ‘--geoarea’, when the clump or object is projected onto
1589     the first two dimensions.  This is only available for 3-dimensional
1590     datasets.  When working with Integral Field Unit (IFU) datasets,
1591     this projection onto the first two dimensions would be a
1592     narrow-band image.
1593
1594‘-A’
1595‘--semimajor’
1596     The pixel-value weighted root mean square (RMS) along the
1597     semi-major axis of the profile (assuming it is an ellipse) in units
1598     of pixels.  See *note Measuring elliptical parameters::.
1599
1600‘-B’
1601‘--semiminor’
1602     The pixel-value weighted root mean square (RMS) along the
1603     semi-minor axis of the profile (assuming it is an ellipse) in units
1604     of pixels.  See *note Measuring elliptical parameters::.
1605
1606‘--axisratio’
1607     The pixel-value weighted axis ratio (semi-minor/semi-major) of the
1608     object or clump.
1609
1610‘-p’
1611‘--positionangle’
1612     The pixel-value weighted angle of the semi-major axis with the
1613     first FITS axis in degrees.  See *note Measuring elliptical
1614     parameters::.
1615
1616‘--geosemimajor’
1617     The geometric (ignoring pixel values) root mean square (RMS) along
1618     the semi-major axis of the profile, assuming it is an ellipse, in
1619     units of pixels.
1620
1621‘--geosemiminor’
1622     The geometric (ignoring pixel values) root mean square (RMS) along
1623     the semi-minor axis of the profile, assuming it is an ellipse, in
1624     units of pixels.
1625
1626‘--geoaxisratio’
1627     The geometric (ignoring pixel values) axis ratio of the profile,
1628     assuming it is an ellipse.
1629
1630‘--geopositionangle’
1631     The geometric (ignoring pixel values) angle of the semi-major axis
1632     with the first FITS axis in degrees.
1633
1634   Above, all of MakeCatalog’s single-valued measurements were listed.
1635As mentioned in the start of this section, MakeCatalog can also do
1636multi-valued measurements per label.  Currently the only such
1637measurement is the creation of spectra from 3D data cubes as discussed
1638below:
1639
1640‘--spectrum’
1641     Generate a spectrum (measurement along the first two FITS
1642     dimensions) for each label when the input dataset is a 3D data
1643     cube.  With this option, a seprate table/spectrum will be generated
1644     for every label.  If the output is a FITS file, each label’s
1645     spectrum will be written into an extension of that file with a
1646     standard name of ‘SPECTRUM_NN’ (the label will be replaced with
1647     ‘NN’).  If the output is a plain text file, each label’s spectrum
1648     will be written into a separate file with the suffix ‘spec-NN.txt’.
1649     See *note MakeCatalog output:: for more on specifying MakeCatalog’s
1650     output file.
1651
1652     The spectra will contain one row for every slice (third FITS
1653     dimension) of the cube.  Since the physical nature of the third
1654     dimension is different, two types of spectra (along with their
1655     errors) are measured: 1) Sum of values in each slice that only have
1656     the requested label.  2) Sum of values on the 2D projection of the
1657     whole label (the area of this projection can be requested with the
1658     ‘--areaxy’ column above).
1659
1660     Labels can overlap when they are projected onto the first two FITS
1661     dimensions (the spatial domain).  To help separate them,
1662     MakeCatalog does a third measurement on each slice: the area, sum
1663     of values and error of all pixels that belong to other labels but
1664     overlap with the 2D projection.  This can be used to see how
1665     reliable the emission line measurement is (on the projected
1666     spectra) and also if multiple lines (labeled regions) belong to the
1667     same physical object.
1668
1669‘--inbetweenints’
1670     Output will contain one row for all integers between 1 and the
1671     largest label in the input (irrespective of their existance in the
1672     input image).  By default, MakeCatalog’s output will only contain
1673     rows with integers that actually corresponded to at least one pixel
1674     in the input dataset.
1675
1676     For example if the input’s only labeled pixel values are 11 and 13,
1677     MakeCatalog’s default output will only have two rows.  If you use
1678     this option, it will have 13 rows and all the columns corresponding
1679     to integer identifiers that didn’t correspond to any pixel will be
1680     0 or NaN (depending on context).
1681
1682
1683File: gnuastro.info,  Node: MakeCatalog output,  Prev: MakeCatalog measurements,  Up: Invoking astmkcatalog
1684
16857.4.6.4 MakeCatalog output
1686..........................
1687
1688After it has completed all the requested measurements (see *note
1689MakeCatalog measurements::), MakeCatalog will store its measurements in
1690table(s).  If an output filename is given (see ‘--output’ in *note Input
1691output options::), the format of the table will be deduced from the
1692name.  When it isn’t given, the input name will be appended with a
1693‘_cat’ suffix (see *note Automatic output::) and its format will be
1694determined from the ‘--tableformat’ option, which is also discussed in
1695*note Input output options::.  ‘--tableformat’ is also necessary when
1696the requested output name is a FITS table (recall that FITS can accept
1697ASCII and binary tables, see *note Table::).
1698
1699   By default (when ‘--spectrum’ or ‘--clumpscat’ aren’t called) only a
1700single catalog/table will be created for the labeled “objects”.
1701
1702   • if ‘--clumpscat’ is called, a secondary catalog/table will also be
1703     created for “clumps” (one of the outputs of the Segment program,
1704     for more on “objects” and “clumps”, see *note Segment::).  In
1705     short, if you only have one labeled image, you don’t have to worry
1706     about clumps and just ignore this.
1707   • When ‘--spectrum’ is called, it is not mandatory to specify any
1708     single-valued measurement columns.  In this case, the output will
1709     only be the spectra of each labeled region within a 3D datacube.
1710     For more, see the description of ‘--spectrum’ in *note MakeCatalog
1711     measurements::.
1712
1713   When possible, MakeCatalog will also measure the full input’s noise
1714level (also known as surface brightness limit, see *note Quantifying
1715measurement limits::).  Since these measurements are related to the
1716noise and not any particular labeled object, they are stored as keywords
1717in the output table.  Furthermore, they are only possible when a
1718standard deviation image has been loaded (done automatically for any
1719column measurement that involves noise, for example ‘--sn’ or ‘--std’).
1720But if you just want the surface brightness limit and no noise-related
1721column, you can use ‘--forcereadstd’.  All these keywords start with
1722‘SBL’ (for “surface brightness limit”) and are described below:
1723
1724‘SBLSTD’
1725     Per-pixel standard deviation.  If a ‘MEDSTD’ keyword exists in the
1726     standard deviation dataset, then that value is directly used.
1727
1728‘SBLNSIG’
1729     Sigma multiple for surface brightness limit (value you gave to
1730     ‘--sfmagnsigma’), used for ‘SBLMAGPX’ and ‘SBLMAG’.
1731
1732‘SBLMAGPX’
1733     Per-pixel surface brightness limit (in units of magnitudes/pixel).
1734
1735‘SBLAREA’
1736     Area (in units of arcsec$^2$) used in ‘SBLMAG’ (value you gave to
1737     ‘--sfmagarea’).
1738
1739‘SBLMAG’
1740     Surface brightness limit of data calculated over ‘SBLAREA’ (in
1741     units of mag/arcsec$^2$).
1742
1743   When any of the upper-limit measurements are requested, the input
1744parameters for the upper-limit measurement are stored in the keywords
1745starting with ‘UP’: ‘UPNSIGMA’, ‘UPNUMBER’, ‘UPRNGNAM’, ‘UPRNGSEE’,
1746‘UPSCMLTP’, ‘UPSCTOL’.  These are primarily input arguments, so they
1747correspond to the options with a similar name.
1748
1749   The full list of MakeCatalog’s options relating to the output file
1750format and keywords are listed below.  See *note MakeCatalog
1751measurements:: for specifying which columns you want in the final
1752catalog.
1753
1754‘-C’
1755‘--clumpscat’
1756     Do measurements on clumps and produce a second catalog (only
1757     devoted to clumps).  When this option is given, MakeCatalog will
1758     also look for a secondary labeled dataset (identifying
1759     substructure) and produce a catalog from that.  For more on the
1760     definition on “clumps”, see *note Segment::.
1761
1762     When the output is a FITS file, the objects and clumps
1763     catalogs/tables will be stored as multiple extensions of one FITS
1764     file.  You can use *note Table:: to inspect the column meta-data
1765     and contents in this case.  However, in plain text format (see
1766     *note Gnuastro text table format::), it is only possible to keep
1767     one table per file.  Therefore, if the output is a text file, two
1768     output files will be created, ending in ‘_o.txt’ (for objects) and
1769_c.txt’ (for clumps).
1770
1771‘--noclumpsort’
1772     Don’t sort the clumps catalog based on object ID (only relevant
1773     with ‘--clumpscat’).  This option will benefit the performance(1)
1774     of MakeCatalog when it is run on multiple threads _and_ the
1775     position of the rows in the clumps catalog is irrelevant (for
1776     example you just want the number-counts).
1777
1778     MakeCatalog does all its measurements on each _object_
1779     independently and in parallel.  As a result, while it is writing
1780     the measurements on each object’s clumps, it doesn’t know how many
1781     clumps there were in previous objects.  Each thread will just fetch
1782     the first available row and write the information of clumps (in
1783     order) starting from that row.  After all the measurements are
1784     done, by default (when this option isn’t called), MakeCatalog will
1785     reorder/permute the clumps catalog to have both the object and
1786     clump ID in an ascending order.
1787
1788     If you would like to order the catalog later (when its a plain text
1789     file), you can run the following command to sort the rows by object
1790     ID (and clump ID within each object), assuming they are
1791     respectively the first and second columns:
1792
1793          $ awk '!/^#/' out_c.txt | sort -g -k1,1 -k2,2
1794
1795‘--sfmagnsigma=FLT’
1796     Value to multiply with the median standard deviation (from a
1797     ‘MEDSTD’ keyword in the Sky standard deviation image) for
1798     estimating the surface brightness limit.  Note that the surface
1799     brightness limit is only reported when a standard deviation image
1800     is read, in other words a column using it is requested (for example
1801     ‘--sn’) or ‘--forcereadstd’ is called.
1802
1803     This value is a per-pixel value, not per object/clump and is not
1804     found over an area or aperture, like the common $5\sigma$ values
1805     that are commonly reported as a measure of depth or the upper-limit
1806     measurements (see *note Quantifying measurement limits::).
1807
1808‘--sfmagarea=FLT’
1809     Area (in arc-seconds squared) to convert the per-pixel estimation
1810     of ‘--sfmagnsigma’ in the comments section of the output tables.
1811     Note that the surface brightness limit is only reported when a
1812     standard deviation image is read, in other words a column using it
1813     is requested (for example ‘--sn’) or ‘--forcereadstd’ is called.
1814
1815     Note that this is just a unit conversion using the World Coordinate
1816     System (WCS) information in the input’s header.  It does not
1817     actually do any measurements on this area.  For random measurements
1818     on any area, please use the upper-limit columns of MakeCatalog (see
1819     the discussion on upper-limit measurements in *note Quantifying
1820     measurement limits::).
1821
1822   ---------- Footnotes ----------
1823
1824   (1) The performance boost due to ‘--noclumpsort’ can only be felt
1825when there are a huge number of objects.  Therefore, by default the
1826output is sorted to avoid miss-understandings or bugs in the user’s
1827scripts when the user forgets to sort the outputs.
1828
1829
1830File: gnuastro.info,  Node: Match,  Prev: MakeCatalog,  Up: Data analysis
1831
18327.5 Match
1833=========
1834
1835Data can come come from different telescopes, filters, software and even
1836different configurations for a single software.  As a result, one of the
1837primary things to do after generating catalogs from each of these
1838sources (for example with *note MakeCatalog::), is to find which sources
1839in one catalog correspond to which in the other(s).  In other words, to
1840‘match’ the two catalogs with each other.
1841
1842   Gnuastro’s Match program is in charge of such operations.  The
1843nearest objects in the two catalogs, within the given aperture, will be
1844found and given as output.  The aperture can be a circle or an ellipse
1845with any orientation.
1846
1847* Menu:
1848
1849* Invoking astmatch::           Inputs, outputs and options of Match
1850
1851
1852File: gnuastro.info,  Node: Invoking astmatch,  Prev: Match,  Up: Match
1853
18547.5.1 Invoking Match
1855--------------------
1856
1857When given two catalogs, Match finds the rows that are nearest to each
1858other within an input aperture.  The executable name is ‘astmatch’ with
1859the following general template
1860
1861     $ astmatch [OPTION ...] input-1 input-2
1862
1863One line examples:
1864
1865     ## 1D wavelength match (within 5 angstroms) of the two inputs.
1866     ## The wavelengths are in the 5th and 10th columns respectively.
1867     $ astmatch --aperture=5e-10 --ccol1=5 --ccol2=10 in1.fits in2.txt
1868
1869     ## Match the two catalogs with a circular aperture of width 2.
1870     ## (Units same as given positional columns).
1871     ## (By default two columns are given for `--ccol1' and `--ccol2',
1872     ##  The number of values to these determines the dimensions).
1873     $ astmatch --aperture=2 input1.txt input2.fits
1874
1875     ## Similar to before, but the output is created by merging various
1876     ## columns from the two inputs: columns 1, RA, DEC from the first
1877     ## input, followed by all columns starting with `MAG' and the `BRG'
1878     ## column from second input and finally the 10th from first input.
1879     $ astmatch --aperture=2 input1.txt input2.fits                   \
1880                --outcols=a1,aRA,aDEC,b/^MAG/,bBRG,a10
1881
1882     ## Assuming both inputs have the same column metadata (same name
1883     ## and numeric type), the output will contain all the rows of the
1884     ## first input, appended with the non-matching rows of the second
1885     ## input (good when you need to merge multiple catalogs that
1886     ## may have matching items, which you don't want to repeat).
1887     $ astmatch input1.fits input2.fits --ccol1=RA,DEC --ccol2=RA,DEC \
1888                --aperture=1/3600 --notmatched --outcols=_all
1889
1890     ## Match the two catalogs within an elliptical aperture of 1 and 2
1891     ## arc-seconds along RA and Dec respectively.
1892     $ astmatch --aperture=1/3600,2/3600 in1.fits in2.txt
1893
1894     ## Match the RA and DEC columns of the first input with the RA_D
1895     ## and DEC_D columns of the second within a 0.5 arcseconds aperture.
1896     $ astmatch --ccol1=RA,DEC --ccol2=RA_D,DEC_D --aperture=0.5/3600  \
1897                in1.fits in2.fits
1898
1899     ## Match in 3D (RA, Dec and Wavelength).
1900     $ astmatch --ccol1=2,3,4 --ccol2=2,3,4 -a0.5/3600,0.5/3600,5e-10 \
1901                in1.fits in2.txt
1902
1903   Match will find the rows that are nearest to each other in two
1904catalogs (given some coordinate columns).  Therefore two catalogs are
1905necessary for input.  However, they don’t necessarily have to be files:
19061) the first catalog can also come from the standard input (for example
1907a pipe, see *note Standard input::); 2) when only one point is needed,
1908you can use the ‘--coord’ option to avoid creating a file for the second
1909catalog.  When the inputs are files, they can be plain text tables or
1910FITS tables, for more see *note Tables::.
1911
1912   Match follows the same basic behavior of all Gnuastro programs as
1913fully described in *note Common program behavior::.  If the first input
1914is a FITS file, the common ‘--hdu’ option (see *note Input output
1915options::) should be used to identify the extension.  When the second
1916input is FITS, the extension must be specified with ‘--hdu2’.
1917
1918   When ‘--quiet’ is not called, Match will print the number of matches
1919found in standard output (on the command-line).  When matches are found,
1920by default, the output file(s) will be the re-arranged input tables such
1921that the rows match each other: both output tables will have the same
1922number of rows which are matched with each other.  If ‘--outcols’ is
1923called, the output is a single table with rows chosen from either of the
1924two inputs in any order.  If the ‘--logasoutput’ option is called, the
1925output will be a single table with the contents of the log file, see
1926below.  If no matches are found, the columns of the output table(s) will
1927have zero rows (with proper meta-data).
1928
1929   If no output file name is given with the ‘--output’ option, then
1930automatic output *note Automatic output:: will be used to determine the
1931output name(s).  Depending on ‘--tableformat’ (see *note Input output
1932options::), the output will then be a (possibly multi-extension) FITS
1933file or (possibly two) plain text file(s).  When the output is a FITS
1934file, the default re-arranged inputs will be two extensions of the
1935output FITS file.  With ‘--outcols’ and ‘--logasoutput’, the FITS output
1936will be a single table (in one extension).
1937
1938   When the ‘--log’ option is called (see *note Operating mode
1939options::), and there was a match, Match will also create a file named
1940astmatch.fits’ (or ‘astmatch.txt’, depending on ‘--tableformat’, see
1941*note Input output options::) in the directory it is run in.  This log
1942table will have three columns.  The first and second columns show the
1943matching row/record number (counting from 1) of the first and second
1944input catalogs respectively.  The third column is the distance between
1945the two matched positions.  The units of the distance are the same as
1946the given coordinates (given the possible ellipticity, see description
1947of ‘--aperture’ below).  When ‘--logasoutput’ is called, no log file
1948(with a fixed name) will be created.  In this case, the output file
1949(possibly given by the ‘--output’ option) will have the contents of this
1950log file.
1951
1952*‘--log’ isn’t thread-safe*: As described above, when ‘--logasoutput’ is
1953not called, the Log file has a fixed name for all calls to Match.
1954Therefore if a separate log is requested in two simultaneous calls to
1955Match in the same directory, Match will try to write to the same file.
1956This will cause problems like unreasonable log file, undefined behavior,
1957or a crash.
1958
1959‘-H STR’
1960‘--hdu2=STR’
1961     The extension/HDU of the second input if it is a FITS file.  When
1962     it isn’t a FITS file, this option’s value is ignored.  For the
1963     first input, the common option ‘--hdu’ must be used.
1964
1965‘--outcols=STR[,STR,[...]]’
1966     Columns (from both inputs) to write into a single matched table
1967     output.  The value to ‘--outcols’ must be a comma-separated list of
1968     column identifiers (number or name, see *note Selecting table
1969     columns::).  The expected format depends on ‘--notmatched’ and
1970     explained below.  By default (when ‘--nomatched’ is not called),
1971     the number of rows in the output will be equal to the number of
1972     matches.  However, when ‘--notmatched’ is called, all the rows
1973     (from the requested columns) of the first input are placed in the
1974     output, and the not-matched rows of the second input are inserted
1975     afterwards (useful when you want to merge unique entries of
1976     multiple catalogs into one).
1977
1978     Default (only matching rows)
1979          The first character of each string specifies the input
1980          catalog: ‘a’ for the first and ‘b’ for the second.  The rest
1981          of the characters of the string will be directly used to
1982          identify the proper column(s) in the respective table.  See
1983          *note Selecting table columns:: for how columns can be
1984          specified in Gnuastro.
1985
1986          For example the output of ‘--outcols=a1,bRA,bDEC’ will have
1987          three columns: the first column of the first input, along with
1988          the ‘RA’ and ‘DEC’ columns of the second input.
1989
1990          If the string after ‘a’ or ‘b’ is ‘_all’, then all the columns
1991          of the respective input file will be written in the output.
1992          For example the command below will print all the input columns
1993          from the first catalog along with the 5th column from the
1994          second:
1995
1996               $ astmatch a.fits b.fits --outcols=a_all,b5
1997
1998          ‘_all’ can be used multiple times, possibly on both inputs.
1999          Tip: if an input’s column is called ‘_all’ (an unlikely name!)
2000          and you don’t want all the columns from that table the output,
2001          use its column number to avoid confusion.
2002
2003          Another example is given in the one-line examples above.
2004          Compared to the default case (where two tables with all their
2005          columns) are saved separately, using this option is much
2006          faster: it will only read and re-arrange the necessary columns
2007          and it will write a single output table.  Combined with
2008          regular expressions in large tables, this can be a very
2009          powerful and convenient way to merge various tables into one.
2010
2011          When ‘--coord’ is given, no second catalog will be read.  The
2012          second catalog will be created internally based on the values
2013          given to ‘--coord’.  So column names aren’t defined and you
2014          can only request integer column numbers that are less than the
2015          number of coordinates given to ‘--coord’.  For example if you
2016          want to find the row matching RA of 1.2345 and Dec of 6.7890,
2017          then you should use ‘--coord=1.2345,6.7890’.  But when using
2018          ‘--outcols’, you can’t give ‘bRA’, or ‘b25’.
2019
2020     With ‘--notmatched’
2021          Only the column names/numbers should be given (for example
2022          ‘--outcols=RA,DEC,MAGNITUDE’).  It is assumed that both input
2023          tables have the requested column(s) and that the numerical
2024          data types of each column in each input (with same name) is
2025          the same as the corresponding column in the other.  Therefore
2026          if one input has a ‘MAGNITUDE’ column with a 32-bit floating
2027          point type, but the ‘MAGNITUDE’ column of the other is 64-bit
2028          floating point, Match will crash with an error.  The metadata
2029          of the columns will come from the first input.
2030
2031          As an example, let’s assume ‘input1.txt’ and ‘input2.fits2032          each have a different number of columns and rows.  However,
2033          they both have the ‘RA’ (64-bit floating point), ‘DEC’ (64-bit
2034          floating point) and ‘MAGNITUDE’ (32-bit floating point)
2035          columns.  If ‘input1.txt’ has 100 rows and ‘input2.fits’ has
2036          300 rows (such that 50 of them match within 1 arcsec of the
2037          first), then the output of the command above will have
2038          $100+(300-50)=350$ rows and only three columns.  Other columns
2039          in each catalog, which may be different, are ignored.
2040
2041               $ astmatch input1.txt  --ccol1=RA,DEC \
2042                          input2.fits --ccol2=RA,DEC \
2043                          --aperture=1/3600 \
2044                          --notmatched --outcols=RA,DEC,MAGNITUDE
2045
2046‘-l’
2047‘--logasoutput’
2048     The output file will have the contents of the log file: indexes in
2049     the two catalogs that match with each other along with their
2050     distance, see description of the log file above.
2051
2052     When this option is called, a separate log file will not be created
2053     and the output will not contain any of the input columns (either as
2054     two tables containing the re-arranged columns of each input, or a
2055     single table mixing columns), only their indices in the log format.
2056
2057‘--notmatched’
2058     Write the non-matching rows into the outputs, not the matched ones.
2059     By default, this will produce two output tables, that will not
2060     necessarily have the same number of rows.  However, when called
2061     with ‘--outcols’, its possible to import non-matching rows of the
2062     second into the first.  See the description of ‘--outcols’ for
2063     more.
2064
2065‘-c INT/STR[,INT/STR]’
2066‘--ccol1=INT/STR[,INT/STR]’
2067     The coordinate columns of the first input.  The number of
2068     dimensions for the match is determined by the number of
2069     comma-separated values given to this option.  The values can be the
2070     column number (counting from 1), exact column name or a regular
2071     expression.  For more, see *note Selecting table columns::.  See
2072     the one-line examples above for some usages of this option.
2073
2074‘-C INT/STR[,INT/STR]’
2075‘--ccol2=INT/STR[,INT/STR]’
2076     The coordinate columns of the second input.  See the example in
2077     ‘--ccol1’ for more.
2078
2079‘-d FLT[,FLT]’
2080‘--coord=FLT[,FLT]’
2081     Manually specify the coordinates to match against the given
2082     catalog.  With this option, Match will not look for a second input
2083     file/table and will directly use the coordinates given to this
2084     option.
2085
2086     When this option is called, the output changes in the following
2087     ways: 1) when ‘--outcols’ is specified, for the second input, it
2088     can only accept integer numbers that are less than the number of
2089     values given to this option, see description of that option for
2090     more.  2) By default (when ‘--outcols’ isn’t used), only the
2091     matching row of the first table will be output (a single file), not
2092     two separate files (one for each table).
2093
2094     This option is good when you have a (large) catalog and only want
2095     to match a single coordinate to it (for example to find the nearest
2096     catalog entry to your desired point).  With this option, you can
2097     write the coordinates on the command-line and thus avoid the need
2098     to make a single-row file.
2099
2100‘-a FLT[,FLT[,FLT]]’
2101‘--aperture=FLT[,FLT[,FLT]]’
2102     Parameters of the aperture for matching.  The values given to this
2103     option can be fractions, for example when the position columns are
2104     in units of degrees, ‘1/3600’ can be used to ask for one
2105     arc-second.  The interpretation of the values depends on the
2106     requested dimensions (determined from ‘--ccol1’ and ‘--ccol2’) and
2107     how many values are given to this option.
2108
2109     When multiple objects are found within the aperture, the match is
2110     defined as the nearest one.  In a multi-dimensional dataset, when
2111     the aperture is a general ellipse or ellipsoid (and not a circle or
2112     sphere), the distance is calculated in the elliptical space along
2113     the major axis.  For the defintion of this distance, see $r_{el}$
2114     in *note Defining an ellipse and ellipsoid::.
2115
2116     1D match
2117          The aperture/interval can only take one value: half of the
2118          interval around each point (maximum distance from each point).
2119
2120     2D match
2121          In a 2D match, the aperture can be a circle, an ellipse
2122          aligned in the axes or an ellipse with a rotated major axis.
2123          To simply the usage, you can determine the shape based on the
2124          number of free parameters for each.
2125
2126          1 number
2127               For example ‘--aperture=2’.  The aperture will be a
2128               circle of the given radius.  The value will be in the
2129               same units as the columns in ‘--ccol1’ and ‘--ccol2’).
2130
2131          2 numbers
2132               For example ‘--aperture=3,4e-10’.  The aperture will be
2133               an ellipse (if the two numbers are different) with the
2134               respective value along each dimension.  The numbers are
2135               in units of the first and second axis.  In the example
2136               above, the semi-axis value along the first axis will be 3
2137               (in units of the first coordinate) and along the second
2138               axis will be $4\times10^{-10}$ (in units of the second
2139               coordinate).  Such values can happen if you are comparing
2140               catalogs of a spectra for example.  If more than one
2141               object exists in the aperture, the nearest will be found
2142               along the major axis as described in *note Defining an
2143               ellipse and ellipsoid::.
2144
2145          3 numbers
2146               For example ‘--aperture=2,0.6,30’.  The aperture will be
2147               an ellipse (if the second value is not 1).  The first
2148               number is the semi-major axis, the second is the axis
2149               ratio and the third is the position angle (in degrees).
2150               If multiple matches are found within the ellipse, the
2151               distance (to find the nearest) is calculated along the
2152               major axis in the elliptical space, see *note Defining an
2153               ellipse and ellipsoid::.
2154
2155     3D match
2156          The aperture (matching volume) can be a sphere, an ellipsoid
2157          aligned on the three axises or a genenral ellipsoid rotated in
2158          any direction.  To simplifythe usage, the shape can be
2159          determined based on the number of values given to this option.
2160
2161          1 number
2162               For example ‘--aperture=3’.  The matching volume will be
2163               a sphere of the given radius.  The value is in the same
2164               units as the input coordinates.
2165
2166          3 numbers
2167               For example ‘--aperture=4,5,6e-10’.  The aperture will be
2168               a general ellipsoid with the respective extent along each
2169               dimension.  The numbers must be in the same units as each
2170               axis.  This is very similar to the two number case of 2D
2171               inputs.  See there for more.
2172
2173          6 numbers
2174               For example ‘--aperture=4,0.5,0.6,10,20,30’.  The numbers
2175               represent the full general ellipsoid definition (in any
2176               orientation).  For the definition of a general ellipsoid,
2177               see *note Defining an ellipse and ellipsoid::.  The first
2178               number is the semi-major axis.  The second and third are
2179               the two axis ratios.  The last three are the three Euler
2180               angles in units of degrees in the ZXZ order as fully
2181               described in *note Defining an ellipse and ellipsoid::.
2182
2183
2184File: gnuastro.info,  Node: Modeling and fittings,  Next: High-level calculations,  Prev: Data analysis,  Up: Top
2185
21868 Modeling and fitting
2187**********************
2188
2189In order to fully understand observations after initial analysis on the
2190image, it is very important to compare them with the existing models to
2191be able to further understand both the models and the data.  The tools
2192in this chapter create model galaxies and will provide 2D fittings to be
2193able to understand the detections.
2194
2195* Menu:
2196
2197* MakeProfiles::                Making mock galaxies and stars.
2198* MakeNoise::                   Make (add) noise to an image.
2199
2200
2201File: gnuastro.info,  Node: MakeProfiles,  Next: MakeNoise,  Prev: Modeling and fittings,  Up: Modeling and fittings
2202
22038.1 MakeProfiles
2204================
2205
2206MakeProfiles will create mock astronomical profiles from a catalog,
2207either individually or together in one output image.  In data analysis,
2208making a mock image can act like a calibration tool, through which you
2209can test how successfully your detection technique is able to detect a
2210known set of objects.  There are commonly two aspects to detecting: the
2211detection of the fainter parts of bright objects (which in the case of
2212galaxies fade into the noise very slowly) or the complete detection of
2213an over-all faint object.  Making mock galaxies is the most accurate
2214(and idealistic) way these two aspects of a detection algorithm can be
2215tested.  You also need mock profiles in fitting known functional
2216profiles with observations.
2217
2218   MakeProfiles was initially built for extra galactic studies, so
2219currently the only astronomical objects it can produce are stars and
2220galaxies.  We welcome the simulation of any other astronomical object.
2221The general outline of the steps that MakeProfiles takes are the
2222following:
2223
2224  1. Build the full profile out to its truncation radius in a possibly
2225     over-sampled array.
2226
2227  2. Multiply all the elements by a fixed constant so its total
2228     magnitude equals the desired total magnitude.
2229
2230  3. If ‘--individual’ is called, save the array for each profile to a
2231     FITS file.
2232
2233  4. If ‘--nomerged’ is not called, add the overlapping pixels of all
2234     the created profiles to the output image and abort.
2235
2236   Using input values, MakeProfiles adds the World Coordinate System
2237(WCS) headers of the FITS standard to all its outputs (except PSF
2238images!).  For a simple test on a set of mock galaxies in one image,
2239there is no need for the third step or the WCS information.
2240
2241   However in complicated simulations like weak lensing simulations,
2242where each galaxy undergoes various types of individual transformations
2243based on their position, those transformations can be applied to the
2244different individual images with other programs.  After all the
2245transformations are applied, using the WCS information in each
2246individual profile image, they can be merged into one output image for
2247convolution and adding noise.
2248
2249* Menu:
2250
2251* Modeling basics::             Astronomical modeling basics.
2252* If convolving afterwards::    Considerations for convolving later.
2253* Profile magnitude::           Definition of total profile magnitude.
2254* Invoking astmkprof::          Inputs and Options for MakeProfiles.
2255
2256
2257File: gnuastro.info,  Node: Modeling basics,  Next: If convolving afterwards,  Prev: MakeProfiles,  Up: MakeProfiles
2258
22598.1.1 Modeling basics
2260---------------------
2261
2262In the subsections below, first a review of some very basic information
2263and concepts behind modeling a real astronomical image is given.  You
2264can skip this subsection if you are already sufficiently familiar with
2265these concepts.
2266
2267* Menu:
2268
2269* Defining an ellipse and ellipsoid::  Definition of these important shapes.
2270* PSF::                         Radial profiles for the PSF.
2271* Stars::                       Making mock star profiles.
2272* Galaxies::                    Radial profiles for galaxies.
2273* Sampling from a function::    Sample a function on a pixelated canvas.
2274* Oversampling::                Oversampling the model.
2275
2276
2277File: gnuastro.info,  Node: Defining an ellipse and ellipsoid,  Next: PSF,  Prev: Modeling basics,  Up: Modeling basics
2278
22798.1.1.1 Defining an ellipse and ellipsoid
2280.........................................
2281
2282The PSF, see *note PSF::, and galaxy radial profiles are generally
2283defined on an ellipse.  Therefore, in this section we’ll start defining
2284an ellipse on a pixelated 2D surface.  Labeling the major axis of an
2285ellipse $a$, and its minor axis with $b$, the _axis ratio_ is defined
2286as: $q\equiv b/a$.  The major axis of an ellipse can be aligned in any
2287direction, therefore the angle of the major axis with respect to the
2288horizontal axis of the image is defined to be the _position angle_ of
2289the ellipse and in this book, we show it with $\theta$.
2290
2291   Our aim is to put a radial profile of any functional form $f(r)$ over
2292an ellipse.  Hence we need to associate a radius/distance to every point
2293in space.  Let’s define the radial distance $r_{el}$ as the distance on
2294the major axis to the center of an ellipse which is located at $i_c$ and
2295$j_c$ (in other words $r_{el}\equiv{a}$).  We want to find $r_{el}$ of a
2296point located at $(i,j)$ (in the image coordinate system) from the
2297center of the ellipse with axis ratio $q$ and position angle $\theta$.
2298First the coordinate system is rotated(1) by $\theta$ to get the new
2299rotated coordinates of that point $(i_r,j_r)$:
2300
2301           $$i_r(i,j)=+(i_c-i)\cos\theta+(j_c-j)\sin\theta$$
2302           $$j_r(i,j)=-(i_c-i)\sin\theta+(j_c-j)\cos\theta$$
2303
2304Recall that an ellipse is defined by $(i_r/a)^2+(j_r/b)^2=1$ and that we
2305defined $r_{el}\equiv{a}$.  Hence, multiplying all elements of the
2306ellipse definition with $r_{el}^2$ we get the elliptical distance at
2307this point point located: $r_{el}=\sqrt{i_r^2+(j_r/q)^2}$.  To place the
2308radial profiles explained below over an ellipse, $f(r_{el})$ is
2309calculated based on the functional radial profile desired.
2310
2311   An ellipse in 3D, or an ellipsoid
2312(https://en.wikipedia.org/wiki/Ellipsoid), can be defined following
2313similar principles as before.  Labeling the major (largest) axis length
2314as $a$, the second and third (in a right-handed coordinate system) axis
2315lengths can be labeled as $b$ and $c$.  Hence we have two axis ratios:
2316$q_1\equiv{b/a}$ and $q_2\equiv{c/a}$.  The orientation of the ellipsoid
2317can be defined from the orientation of its major axis.  There are many
2318ways to define 3D orientation and order matters.  So to be clear, here
2319we use the ZXZ (or $Z_1X_2Z_3$) proper Euler angles
2320(https://en.wikipedia.org/wiki/Euler_angles) to define the 3D
2321orientation.  In short, when a point is rotated in this order, we first
2322rotate it around the Z axis (third axis) by $\alpha$, then about the
2323(rotated) X axis by $\beta$ and finally about the (rotated) Z axis by
2324$\gamma$.
2325
2326   Following the discussion in *note Merging multiple warpings::, we can
2327define the full rotation with the following matrix multiplication.
2328However, here we are rotating the coordinates, not the point.
2329Therefore, both the rotation angles and rotation order are reversed.  We
2330are also not using homogeneous coordinates (see *note Warping basics::)
2331since we aren’t concerned with translation in this context:
2332
2333              $$\left[\matrix{i_r\cr j_r\cr k_r}\right] =
2334\left[\matrix{cos\gamma&sin\gamma&0\cr -sin\gamma&cos\gamma&0\cr 0&0&1}\right]
2335\left[\matrix{1&0&0\cr 0&cos\beta&sin\beta\cr 0&-sin\beta&cos\beta }\right]
2336\left[\matrix{cos\alpha&sin\alpha&0\cr -sin\alpha&cos\alpha&0\cr 0&0&1}\right]
2337           \left[\matrix{i_c-i\cr j_c-j\cr k_c-k}\right] $$
2338
2339Recall that an ellipsoid can be characterized with
2340$(i_r/a)^2+(j_r/b)^2+(k_r/c)^2=1$, so similar to before
2341($r_{el}\equiv{a}$), we can find the ellipsoidal radius at pixel
2342$(i,j,k)$ as: $r_{el}=\sqrt{i_r^2+(j_r/q_1)^2+(k_r/q_2)^2}$.
2343
2344   MakeProfiles builds the profile starting from the nearest element
2345(pixel in an image) in the dataset to the profile center.  The profile
2346value is calculated for that central pixel using monte carlo
2347integration, see *note Sampling from a function::.  The next pixel is
2348the next nearest neighbor to the central pixel as defined by $r_{el}$.
2349This process goes on until the profile is fully built upto the
2350truncation radius.  This is done fairly efficiently using a breadth
2351first parsing strategy(2) which is implemented through an ordered linked
2352list.
2353
2354   Using this approach, we build the profile by expanding the
2355circumference.  Not one more extra pixel has to be checked (the
2356calculation of $r_{el}$ from above is not cheap in CPU terms).  Another
2357consequence of this strategy is that extending MakeProfiles to three
2358dimensions becomes very simple: only the neighbors of each pixel have to
2359be changed.  Everything else after that (when the pixel index and its
2360radial profile have entered the linked list) is the same, no matter the
2361number of dimensions we are dealing with.
2362
2363   ---------- Footnotes ----------
2364
2365   (1) Do not confuse the signs of $sin$ with the rotation matrix
2366defined in *note Warping basics::.  In that equation, the point is
2367rotated, here the coordinates are rotated and the point is fixed.
2368
2369   (2) <http://en.wikipedia.org/wiki/Breadth-first_search>
2370
2371
2372File: gnuastro.info,  Node: PSF,  Next: Stars,  Prev: Defining an ellipse and ellipsoid,  Up: Modeling basics
2373
23748.1.1.2 Point spread function
2375.............................
2376
2377Assume we have a ‘point’ source, or a source that is far smaller than
2378the maximum resolution (a pixel).  When we take an image of it, it will
2379‘spread’ over an area.  To quantify that spread, we can define a
2380‘function’.  This is how the point spread function or the PSF of an
2381image is defined.  This ‘spread’ can have various causes, for example in
2382ground based astronomy, due to the atmosphere.  In practice we can never
2383surpass the ‘spread’ due to the diffraction of the lens aperture.
2384Various other effects can also be quantified through a PSF. For example,
2385the simple fact that we are sampling in a discrete space, namely the
2386pixels, also produces a very small ‘spread’ in the image.
2387
2388   Convolution is the mathematical process by which we can apply a
2389‘spread’ to an image, or in other words blur the image, see *note
2390Convolution process::.  The Brightness of an object should remain
2391unchanged after convolution, see *note Brightness flux magnitude::.
2392Therefore, it is important that the sum of all the pixels of the PSF be
2393unity.  The PSF image also has to have an odd number of pixels on its
2394sides so one pixel can be defined as the center.  In MakeProfiles, the
2395PSF can be set by the two methods explained below.
2396
2397Parametric functions
2398     A known mathematical function is used to make the PSF. In this
2399     case, only the parameters to define the functions are necessary and
2400     MakeProfiles will make a PSF based on the given parameters for each
2401     function.  In both cases, the center of the profile has to be
2402     exactly in the middle of the central pixel of the PSF (which is
2403     automatically done by MakeProfiles).  When talking about the PSF,
2404     usually, the full width at half maximum or FWHM is used as a scale
2405     of the width of the PSF.
2406
2407     ‘Gaussian’
2408          In the older papers, and to a lesser extent even today, some
2409          researchers use the 2D Gaussian function to approximate the
2410          PSF of ground based images.  In its most general form, a
2411          Gaussian function can be written as:
2412
2413      $$f(r)=a \exp \left( -(x-\mu)^2 \over 2\sigma^2 \right)+d$$
2414
2415          Since the center of the profile is pre-defined, $\mu$ and $d$
2416          are constrained.  $a$ can also be found because the function
2417          has to be normalized.  So the only important parameter for
2418          MakeProfiles is the $\sigma$.  In the Gaussian function we
2419          have this relation between the FWHM and $\sigma$:
2420
2421      $$\rm{FWHM}_g=2\sqrt{2\ln{2}}\sigma \approx 2.35482\sigma$$
2422
2423     ‘Moffat’
2424          The Gaussian profile is much sharper than the images taken
2425          from stars on photographic plates or CCDs.  Therefore in 1969,
2426          Moffat proposed this functional form for the image of stars:
2427
2428  $$f(r)=a \left[ 1+\left( r\over \alpha \right)^2 \right]^{-\beta}$$
2429
2430          Again, $a$ is constrained by the normalization, therefore two
2431          parameters define the shape of the Moffat function: $\alpha$
2432          and $\beta$.  The radial parameter is $\alpha$ which is
2433          related to the FWHM by
2434
2435              $$\rm{FWHM}_m=2\alpha\sqrt{2^{1/\beta}-1}$$
2436
2437          Comparing with the PSF predicted from atmospheric turbulence
2438          theory with a Moffat function, Trujillo et al.(1)  claim that
2439          $\beta$ should be 4.765.  They also show how the Moffat PSF
2440          contains the Gaussian PSF as a limiting case when
2441          $\beta\to\infty$.
2442
2443An input FITS image
2444     An input image file can also be specified to be used as a PSF. If
2445     the sum of its pixels are not equal to 1, the pixels will be
2446     multiplied by a fraction so the sum does become 1.
2447
2448   While the Gaussian is only dependent on the FWHM, the Moffat function
2449is also dependent on $\beta$.  Comparing these two functions with a
2450fixed FWHM gives the following results:
2451
2452   • Within the FWHM, the functions don’t have significant differences.
2453   • For a fixed FWHM, as $\beta$ increases, the Moffat function becomes
2454     sharper.
2455   • The Gaussian function is much sharper than the Moffat functions,
2456     even when $\beta$ is large.
2457
2458   ---------- Footnotes ----------
2459
2460   (1) Trujillo, I., J. A. L. Aguerri, J. Cepa, and C. M. Gutierrez
2461(2001).  “The effects of seeing on Sérsic profiles - II. The Moffat
2462PSF”.  In: MNRAS 328, pp.  977—985.
2463
2464
2465File: gnuastro.info,  Node: Stars,  Next: Galaxies,  Prev: PSF,  Up: Modeling basics
2466
24678.1.1.3 Stars
2468.............
2469
2470In MakeProfiles, stars are generally considered to be a point source.
2471This is usually the case for extra galactic studies, were nearby stars
2472are also in the field.  Since a star is only a point source, we assume
2473that it only fills one pixel prior to convolution.  In fact, exactly for
2474this reason, in astronomical images the light profiles of stars are one
2475of the best methods to understand the shape of the PSF and a very large
2476fraction of scientific research is preformed by assuming the shapes of
2477stars to be the PSF of the image.
2478
2479
2480File: gnuastro.info,  Node: Galaxies,  Next: Sampling from a function,  Prev: Stars,  Up: Modeling basics
2481
24828.1.1.4 Galaxies
2483................
2484
2485Today, most practitioners agree that the flux of galaxies can be modeled
2486with one or a few generalized de Vaucouleur’s (or Sérsic) profiles.
2487
2488$$I(r) = I_e \exp \left ( -b_n \left[ \left( r \over r_e \right)^{1/n} -1 \right] \right )$$
2489
2490   Gérard de Vaucouleurs (1918-1995) was first to show in 1948 that this
2491function resembles the galaxy light profiles, with the only difference
2492that he held $n$ fixed to a value of 4.  Twenty years later in 1968, J.
2493L. Sérsic showed that $n$ can have a variety of values and does not
2494necessarily need to be 4.  This profile depends on the effective radius
2495($r_e$) which is defined as the radius which contains half of the
2496profile brightness (see *note Profile magnitude::).  $I_e$ is the flux
2497at the effective radius.  The Sérsic index $n$ is used to define the
2498concentration of the profile within $r_e$ and $b_n$ is a constant
2499dependent on $n$.  MacArthur et al.(1)  show that for $n>0.35$, $b_n$
2500can be accurately approximated using this equation:
2501
2502$$b_n=2n - {1\over 3} + {4\over 405n} + {46\over 25515n^2} + {131\over 1148175n^3}-{2194697\over 30690717750n^4}$$
2503
2504   ---------- Footnotes ----------
2505
2506   (1) MacArthur, L. A., S. Courteau, and J. A. Holtzman (2003).
2507“Structure of Disk-dominated Galaxies.  I. Bulge/Disk Parameters,
2508Simulations, and Secular Evolution”.  In: ApJ 582, pp.  689—722.
2509
2510
2511File: gnuastro.info,  Node: Sampling from a function,  Next: Oversampling,  Prev: Galaxies,  Up: Modeling basics
2512
25138.1.1.5 Sampling from a function
2514................................
2515
2516A pixel is the ultimate level of accuracy to gather data, we can’t get
2517any more accurate in one image, this is known as sampling in signal
2518processing.  However, the mathematical profiles which describe our
2519models have infinite accuracy.  Over a large fraction of the area of
2520astrophysically interesting profiles (for example galaxies or PSFs), the
2521variation of the profile over the area of one pixel is not too
2522significant.  In such cases, the elliptical radius ($r_{el}$ of the
2523center of the pixel can be assigned as the final value of the pixel, see
2524*note Defining an ellipse and ellipsoid::).
2525
2526   As you approach their center, some galaxies become very sharp (their
2527value significantly changes over one pixel’s area).  This sharpness
2528increases with smaller effective radius and larger Sérsic values.  Thus
2529rendering the central value extremely inaccurate.  The first method that
2530comes to mind for solving this problem is integration.  The functional
2531form of the profile can be integrated over the pixel area in a 2D
2532integration process.  However, unfortunately numerical integration
2533techniques also have their limitations and when such sharp profiles are
2534needed they can become extremely inaccurate.
2535
2536   The most accurate method of sampling a continuous profile on a
2537discrete space is by choosing a large number of random points within the
2538boundaries of the pixel and taking their average value (or Monte Carlo
2539integration).  This is also, generally speaking, what happens in
2540practice with the photons on the pixel.  The number of random points can
2541be set with ‘--numrandom’.
2542
2543   Unfortunately, repeating this Monte Carlo process would be extremely
2544time and CPU consuming if it is to be applied to every pixel.  In order
2545to not loose too much accuracy, in MakeProfiles, the profile is built
2546using both methods explained below.  The building of the profile begins
2547from its central pixel and continues (radially) outwards.  Monte Carlo
2548integration is first applied (which yields $F_r$), then the central
2549pixel value ($F_c$) is calculated on the same pixel.  If the fractional
2550difference ($|F_r-F_c|/F_r$) is lower than a given tolerance level
2551(specified with ‘--tolerance’) MakeProfiles will stop using Monte Carlo
2552integration and only use the central pixel value.
2553
2554   The ordering of the pixels in this inside-out construction is based
2555on $r=\sqrt{(i_c-i)^2+(j_c-j)^2}$, not $r_{el}$, see *note Defining an
2556ellipse and ellipsoid::.  When the axis ratios are large (near one) this
2557is fine.  But when they are small and the object is highly elliptical,
2558it might seem more reasonable to follow $r_{el}$ not $r$.  The problem
2559is that the gradient is stronger in pixels with smaller $r$ (and larger
2560$r_{el}$) than those with smaller $r_{el}$.  In other words, the
2561gradient is strongest along the minor axis.  So if the next pixel is
2562chosen based on $r_{el}$, the tolerance level will be reached sooner and
2563lots of pixels with large fractional differences will be missed.
2564
2565   Monte Carlo integration uses a random number of points.  Thus, every
2566time you run it, by default, you will get a different distribution of
2567points to sample within the pixel.  In the case of large profiles, this
2568will result in a slight difference of the pixels which use Monte Carlo
2569integration each time MakeProfiles is run.  To have a deterministic
2570result, you have to fix the random number generator properties which is
2571used to build the random distribution.  This can be done by setting the
2572‘GSL_RNG_TYPE’ and ‘GSL_RNG_SEED’ environment variables and calling
2573MakeProfiles with the ‘--envseed’ option.  To learn more about the
2574process of generating random numbers, see *note Generating random
2575numbers::.
2576
2577   The seed values are fixed for every profile: with ‘--envseed’, all
2578the profiles have the same seed and without it, each will get a
2579different seed using the system clock (which is accurate to within one
2580microsecond).  The same seed will be used to generate a random number
2581for all the sub-pixel positions of all the profiles.  So in the former,
2582the sub-pixel points checked for all the pixels undergoing Monte carlo
2583integration in all profiles will be identical.  In other words, the
2584sub-pixel points in the first (closest to the center) pixel of all the
2585profiles will be identical with each other.  All the second pixels
2586studied for all the profiles will also receive an identical (different
2587from the first pixel) set of sub-pixel points and so on.  As long as the
2588number of random points used is large enough or the profiles are not
2589identical, this should not cause any systematic bias.
2590
2591
2592File: gnuastro.info,  Node: Oversampling,  Prev: Sampling from a function,  Up: Modeling basics
2593
25948.1.1.6 Oversampling
2595....................
2596
2597The steps explained in *note Sampling from a function:: do give an
2598accurate representation of a profile prior to convolution.  However, in
2599an actual observation, the image is first convolved with or blurred by
2600the atmospheric and instrument PSF in a continuous space and then it is
2601sampled on the discrete pixels of the camera.
2602
2603   In order to more accurately simulate this process, the unconvolved
2604image and the PSF are created on a finer pixel grid.  In other words,
2605the output image is a certain odd-integer multiple of the desired size,
2606we can call this ‘oversampling’.  The user can specify this multiple as
2607a command-line option.  The reason this has to be an odd number is that
2608the PSF has to be centered on the center of its image.  An image with an
2609even number of pixels on each side does not have a central pixel.
2610
2611   The image can then be convolved with the PSF (which should also be
2612oversampled on the same scale).  Finally, image can be sub-sampled to
2613get to the initial desired pixel size of the output image.  After this,
2614mock noise can be added as explained in the next section.  This is
2615because unlike the PSF, the noise occurs in each output pixel, not on a
2616continuous space like all the prior steps.
2617
2618
2619File: gnuastro.info,  Node: If convolving afterwards,  Next: Profile magnitude,  Prev: Modeling basics,  Up: MakeProfiles
2620
26218.1.2 If convolving afterwards
2622------------------------------
2623
2624In case you want to convolve the image later with a given point spread
2625function, make sure to use a larger image size.  After convolution, the
2626profiles become larger and a profile that is normally completely outside
2627of the image might fall within it.
2628
2629   On one axis, if you want your final (convolved) image to be $m$
2630pixels and your PSF is $2n+1$ pixels wide, then when calling
2631MakeProfiles, set the axis size to $m+2n$, not $m$.  You also have to
2632shift all the pixel positions of the profile centers on the that axis by
2633$n$ pixels to the positive.
2634
2635   After convolution, you can crop the outer $n$ pixels with the section
2636crop box specification of Crop: ‘--section=n:*-n,n:*-n’ assuming your
2637PSF is a square, see *note Crop section syntax::.  This will also remove
2638all discrete Fourier transform artifacts (blurred sides) from the final
2639image.  To facilitate this shift, MakeProfiles has the options
2640‘--xshift’, ‘--yshift’ and ‘--prepforconv’, see *note Invoking
2641astmkprof::.
2642
2643
2644File: gnuastro.info,  Node: Profile magnitude,  Next: Invoking astmkprof,  Prev: If convolving afterwards,  Up: MakeProfiles
2645
26468.1.3 Profile magnitude
2647-----------------------
2648
2649To find the profile brightness or its magnitude, (see *note Brightness
2650flux magnitude::), it is customary to use the 2D integration of the flux
2651to infinity.  However, in MakeProfiles we do not follow this idealistic
2652approach and apply a more realistic method to find the total brightness
2653or magnitude: the sum of all the pixels belonging to a profile within
2654its predefined truncation radius.  Note that if the truncation radius is
2655not large enough, this can be significantly different from the total
2656integrated light to infinity.
2657
2658   An integration to infinity is not a realistic condition because no
2659galaxy extends indefinitely (important for high Sérsic index profiles),
2660pixelation can also cause a significant difference between the actual
2661total pixel sum value of the profile and that of integration to
2662infinity, especially in small and high Sérsic index profiles.  To be
2663safe, you can specify a large enough truncation radius for such compact
2664high Sérsic index profiles.
2665
2666   If oversampling is used then the brightness is calculated using the
2667over-sampled image, see *note Oversampling:: which is much more
2668accurate.  The profile is first built in an array completely bounding it
2669with a normalization constant of unity (see *note Galaxies::).  Taking
2670$B$ to be the desired brightness and $S$ to be the sum of the pixels in
2671the created profile, every pixel is then multiplied by $B/S$ so the sum
2672is exactly $B$.
2673
2674   If the ‘--individual’ option is called, this same array is written to
2675a FITS file.  If not, only the overlapping pixels of this array and the
2676output image are kept and added to the output array.
2677
2678
2679File: gnuastro.info,  Node: Invoking astmkprof,  Prev: Profile magnitude,  Up: MakeProfiles
2680
26818.1.4 Invoking MakeProfiles
2682---------------------------
2683
2684MakeProfiles will make any number of profiles specified in a catalog
2685either individually or in one image.  The executable name is ‘astmkprof’
2686with the following general template
2687
2688     $ astmkprof [OPTION ...] [Catalog]
2689
2690One line examples:
2691
2692     ## Make an image with profiles in catalog.txt (with default size):
2693     $ astmkprof catalog.txt
2694
2695     ## Make the profiles in catalog.txt over image.fits:
2696     $ astmkprof --background=image.fits catalog.txt
2697
2698     ## Make a Moffat PSF with FWHM 3pix, beta=2.8, truncation=5
2699     $ astmkprof --kernel=moffat,3,2.8,5 --oversample=1
2700
2701     ## Make profiles in catalog, using RA and Dec in the given column:
2702     $ astmkprof --ccol=RA_CENTER --ccol=DEC_CENTER --mode=wcs catalog.txt
2703
2704     ## Make a 1500x1500 merged image (oversampled 500x500) image along
2705     ## with an individual image for all the profiles in catalog:
2706     $ astmkprof --individual --oversample 3 --mergedsize=500,500 cat.txt
2707
2708The parameters of the mock profiles can either be given through a
2709catalog (which stores the parameters of many mock profiles, see *note
2710MakeProfiles catalog::), or the ‘--kernel’ option (see *note
2711MakeProfiles output dataset::).  The catalog can be in the FITS ASCII,
2712FITS binary format, or plain text formats (see *note Tables::).  A plain
2713text catalog can also be provided using the Standard input (see *note
2714Standard input::).  The columns related to each parameter can be
2715determined both by number, or by match/search criteria using the column
2716names, units, or comments, with the options ending in ‘col’, see below.
2717
2718   Without any file given to the ‘--background’ option, MakeProfiles
2719will make a zero-valued image and build the profiles on that (its size
2720and main WCS parameters can also be defined through the options
2721described in *note MakeProfiles output dataset::).  Besides the
2722main/merged image containing all the profiles in the catalog, it is also
2723possible to build individual images for each profile (only enclosing one
2724full profile to its truncation radius) with the ‘--individual’ option.
2725
2726   If an image is given to the ‘--background’ option, the pixels of that
2727image are used as the background value for every pixel hence flux value
2728of each profile pixel will be added to the pixel in that background
2729value.  You can disable this with the ‘--clearcanvas’ option (which will
2730initialize the background to zero-valued pixels and build profiles over
2731that).  With the ‘--background’ option, the values to all options
2732relating to the “canvas” (output size and WCS) will be ignored if
2733specified, for example ‘--oversample’, ‘--mergedsize’, and
2734‘--prepforconv’.
2735
2736   The sections below discuss the options specific to MakeProfiles based
2737on context: the input catalog settings which can have many rows for
2738different profiles are discussed in *note MakeProfiles catalog::, in
2739*note MakeProfiles profile settings::, we discuss how you can set
2740general profile settings (that are the same for all the profiles in the
2741catalog).  Finally *note MakeProfiles output dataset:: and *note
2742MakeProfiles log file:: discuss the outputs of MakeProfiles and how you
2743can configure them.  Besides these, MakeProfiles also supports all the
2744common Gnuastro program options that are discussed in *note Common
2745options::, so please flip through them is well for a more comfortable
2746usage.
2747
2748   When building 3D profiles, there are more degrees of freedom.  Hence,
2749more columns are necessary and all the values related to dimensions (for
2750example size of dataset in each dimension and the WCS properties) must
2751also have 3 values.  To allow having an independent set of default
2752values for creating 3D profiles, MakeProfiles also installs a
2753astmkprof-3d.conf’ configuration file (see *note Configuration
2754files::).  You can use this for default 3D profile values.  For example,
2755if you installed Gnuastro with the prefix ‘/usr/local’ (the default
2756location, see *note Installation directory::), you can benefit from this
2757configuration file by running MakeProfiles like the example below.  As
2758with all configuration files, if you want to customize a given option,
2759call it before the configuration file.
2760
2761     $ astmkprof --config=/usr/local/etc/astmkprof-3d.conf catalog.txt
2762
2763   To further simplify the process, you can define a shell alias in any
2764startup file (for example ‘~/.bashrc’, see *note Installation
2765directory::).  Assuming that you installed Gnuastro in ‘/usr/local’, you
2766can add this line to the startup file (you may put it all in one line,
2767it is broken into two lines here for fitting within page limits).
2768
2769     alias astmkprof-3d="astmkprof --config=/usr/local/etc/astmkprof-3d.conf"
2770
2771Using this alias, you can call MakeProfiles with the name ‘astmkprof-3d’
2772(instead of ‘astmkprof’).  It will automatically load the 3D specific
2773configuration file first, and then parse any other arguments, options or
2774configuration files.  You can change the default values in this 3D
2775configuration file by calling them on the command-line as you do with
2776‘astmkprof’(1).
2777
2778   Please see *note Sufi simulates a detection:: for a very complete
2779tutorial explaining how one could use MakeProfiles in conjunction with
2780other Gnuastro’s programs to make a complete simulated image of a mock
2781galaxy.
2782
2783* Menu:
2784
2785* MakeProfiles catalog::        Required catalog properties.
2786* MakeProfiles profile settings::  Configuration parameters for all profiles.
2787* MakeProfiles output dataset::  The canvas/dataset to build profiles over.
2788* MakeProfiles log file::       A description of the optional log file.
2789
2790   ---------- Footnotes ----------
2791
2792   (1) Recall that for single-invocation options, the last command-line
2793invocation takes precedence over all previous invocations (including
2794those in the 3D configuration file).  See the description of ‘--config’
2795in *note Operating mode options::.
2796
2797
2798File: gnuastro.info,  Node: MakeProfiles catalog,  Next: MakeProfiles profile settings,  Prev: Invoking astmkprof,  Up: Invoking astmkprof
2799
28008.1.4.1 MakeProfiles catalog
2801............................
2802
2803The catalog containing information about each profile can be in the FITS
2804ASCII, FITS binary, or plain text formats (see *note Tables::).  The
2805latter can also be provided using standard input (see *note Standard
2806input::).  Its columns can be ordered in any desired manner.  You can
2807specify which columns belong to which parameters using the set of
2808options discussed below.  For example through the ‘--rcol’ and ‘--tcol’
2809options, you can specify the column that contains the radial parameter
2810for each profile and its truncation respectively.  See *note Selecting
2811table columns:: for a thorough discussion on the values to these
2812options.
2813
2814   The value for the profile center in the catalog (the ‘--ccol’ option)
2815can be a floating point number so the profile center can be on any
2816sub-pixel position.  Note that pixel positions in the FITS standard
2817start from 1 and an integer is the pixel center.  So a 2D image actually
2818starts from the position (0.5, 0.5), which is the bottom-left corner of
2819the first pixel.  When a ‘--background’ image with WCS information is
2820provided or you specify the WCS parameters with the respective options,
2821you may also use RA and Dec to identify the center of each profile (see
2822the ‘--mode’ option below).
2823
2824   In MakeProfiles, profile centers do not have to be in (overlap with)
2825the final image.  Even if only one pixel of the profile within the
2826truncation radius overlaps with the final image size, the profile is
2827built and included in the final image image.  Profiles that are
2828completely out of the image will not be created (unless you explicitly
2829ask for it with the ‘--individual’ option).  You can use the output log
2830file (created with ‘--log’ to see which profiles were within the image,
2831see *note Common options::.
2832
2833   If PSF profiles (Moffat or Gaussian, see *note PSF::) are in the
2834catalog and the profiles are to be built in one image (when
2835‘--individual’ is not used), it is assumed they are the PSF(s) you want
2836to convolve your created image with.  So by default, they will not be
2837built in the output image but as separate files.  The sum of pixels of
2838these separate files will also be set to unity (1) so you are ready to
2839convolve, see *note Convolution process::.  As a summary, the position
2840and magnitude of PSF profile will be ignored.  This behavior can be
2841disabled with the ‘--psfinimg’ option.  If you want to create all the
2842profiles separately (with ‘--individual’) and you want the sum of the
2843PSF profile pixels to be unity, you have to set their magnitudes in the
2844catalog to the zero point magnitude and be sure that the central
2845positions of the profiles don’t have any fractional part (the PSF center
2846has to be in the center of the pixel).
2847
2848   The list of options directly related to the input catalog columns is
2849shown below.
2850
2851‘--ccol=STR/INT2852     Center coordinate column for each dimension.  This option must be
2853     called two times to define the center coordinates in an image.  For
2854     example ‘--ccol=RA’ and ‘--ccol=DEC’ (along with ‘--mode=wcs’) will
2855     inform MakeProfiles to look into the catalog columns named ‘RA’ and
2856     ‘DEC’ for the Right Ascension and Declination of the profile
2857     centers.
2858
2859‘--fcol=INT/STR2860     The functional form of the profile with one of the values below
2861     depending on the desired profile.  The column can contain either
2862     the numeric codes (for example ‘‘1’’) or string characters (for
2863     example ‘‘sersic’’).  The numeric codes are easier to use in
2864     scripts which generate catalogs with hundreds or thousands of
2865     profiles.
2866
2867     The string format can be easier when the catalog is to be
2868     written/checked by hand/eye before running MakeProfiles.  It is
2869     much more readable and provides a level of documentation.  All
2870     Gnuastro’s recognized table formats (see *note Recognized table
2871     formats::) accept string type columns.  To have string columns in a
2872     plain text table/catalog, see *note Gnuastro text table format::.
2873
2874        • Sérsic profile with ‘‘sersic’’ or ‘‘1’’.
2875
2876        • Moffat profile with ‘‘moffat’’ or ‘‘2’’.
2877
2878        • Gaussian profile with ‘‘gaussian’’ or ‘‘3’’.
2879
2880        • Point source with ‘‘point’’ or ‘‘4’’.
2881
2882        • Flat profile with ‘‘flat’’ or ‘‘5’’.
2883
2884        • Circumference profile with ‘‘circum’’ or ‘‘6’’.  A fixed value
2885          will be used for all pixels less than or equal to the
2886          truncation radius ($r_t$) and greater than $r_t-w$ ($w$ is the
2887          value to the ‘--circumwidth’).
2888
2889        • Radial distance profile with ‘‘distance’’ or ‘‘7’’.  At the
2890          lowest level, each pixel only has an elliptical radial
2891          distance given the profile’s shape and orientation (see *note
2892          Defining an ellipse and ellipsoid::).  When this profile is
2893          chosen, the pixel’s elliptical radial distance from the
2894          profile center is written as its value.  For this profile, the
2895          value in the magnitude column (‘--mcol’) will be ignored.
2896
2897          You can use this for checks or as a first approximation to
2898          define your own higher-level radial function.  In the latter
2899          case, just note that the central values are going to be
2900          incorrect (see *note Sampling from a function::).
2901
2902        • Custom profile with ‘‘custom’’ or ‘‘8’’.  The values to use
2903          for each radial interval should be in the table given to
2904          ‘--customtable’.  For more, see *note MakeProfiles profile
2905          settings::.
2906
2907        • Azimuthal angle profile with ‘‘azimuth’’ or ‘‘9’’.  Every
2908          pixel within the truncation radius will be given its azimuthal
2909          angle (in degrees, from 0 to 360) from the major axis.  In
2910          combination with the radial distance profile, you can now
2911          create complex features in polar coordinates, such as tidal
2912          tails or tidal shocks (using the Arithmetic program to mix the
2913          radius and azimuthal angle through a function to create your
2914          desired features).
2915
2916‘--rcol=STR/INT2917     The radius parameter of the profiles.  Effective radius ($r_e$) if
2918     Sérsic, FWHM if Moffat or Gaussian.
2919
2920‘--ncol=STR/INT2921     The Sérsic index ($n$) or Moffat $\beta$.
2922
2923‘--pcol=STR/INT2924     The position angle (in degrees) of the profiles relative to the
2925     first FITS axis (horizontal when viewed in SAO DS9).  When building
2926     a 3D profile, this is the first Euler angle: first rotation of the
2927     ellipsoid major axis from the first FITS axis (rotating about the
2928     third axis).  See *note Defining an ellipse and ellipsoid::.
2929
2930‘--p2col=STR/INT2931     Second Euler angle (in degrees) when building a 3D ellipsoid.  This
2932     is the second rotation of the ellipsoid major axis (following
2933     ‘--pcol’) about the (rotated) X axis.  See *note Defining an
2934     ellipse and ellipsoid::.  This column is ignored when building a 2D
2935     profile.
2936
2937‘--p3col=STR/INT2938     Third Euler angle (in degrees) when building a 3D ellipsoid.  This
2939     is the third rotation of the ellipsoid major axis (following
2940     ‘--pcol’ and ‘--p2col’) about the (rotated) Z axis.  See *note
2941     Defining an ellipse and ellipsoid::.  This column is ignored when
2942     building a 2D profile.
2943
2944‘--qcol=STR/INT2945     The axis ratio of the profiles (minor axis divided by the major
2946     axis in a 2D ellipse).  When building a 3D ellipse, this is the
2947     ratio of the major axis to the semi-axis length of the second
2948     dimension (in a right-handed coordinate system).  See $q1$ in *note
2949     Defining an ellipse and ellipsoid::.
2950
2951‘--q2col=STR/INT2952     The ratio of the ellipsoid major axis to the third semi-axis length
2953     (in a right-handed coordinate system) of a 3D ellipsoid.  See $q1$
2954     in *note Defining an ellipse and ellipsoid::.  This column is
2955     ignored when building a 2D profile.
2956
2957‘--mcol=STR/INT2958     The total pixelated magnitude of the profile within the truncation
2959     radius, see *note Profile magnitude::.
2960
2961‘--tcol=STR/INT2962     The truncation radius of this profile.  By default it is in units
2963     of the radial parameter of the profile (the value in the ‘--rcol’
2964     of the catalog).  If ‘--tunitinp’ is given, this value is
2965     interpreted in units of pixels (prior to oversampling) irrespective
2966     of the profile.
2967
2968
2969File: gnuastro.info,  Node: MakeProfiles profile settings,  Next: MakeProfiles output dataset,  Prev: MakeProfiles catalog,  Up: Invoking astmkprof
2970
29718.1.4.2 MakeProfiles profile settings
2972.....................................
2973
2974The profile parameters that differ between each created profile are
2975specified through the columns in the input catalog and described in
2976*note MakeProfiles catalog::.  Besides those there are general settings
2977for some profiles that don’t differ between one profile and another,
2978they are a property of the general process.  For example how many random
2979points to use in the monte-carlo integration, this value is fixed for
2980all the profiles.  The options described in this section are for
2981configuring such properties.
2982
2983‘--mode=STR’
2984     Interpret the center position columns (‘--ccol’ in *note
2985     MakeProfiles catalog::) in image or WCS coordinates.  This option
2986     thus accepts only two values: ‘img’ and ‘wcs’.  It is mandatory
2987     when a catalog is being used as input.
2988
2989‘-r’
2990‘--numrandom’
2991     The number of random points used in the central regions of the
2992     profile, see *note Sampling from a function::.
2993
2994‘-e’
2995‘--envseed’
2996     Use the value to the ‘GSL_RNG_SEED’ environment variable to
2997     generate the random Monte Carlo sampling distribution, see *note
2998     Sampling from a function:: and *note Generating random numbers::.
2999
3000‘-t FLT’
3001‘--tolerance=FLT’
3002     The tolerance to switch from Monte Carlo integration to the central
3003     pixel value, see *note Sampling from a function::.
3004
3005‘-p’
3006‘--tunitinp’
3007     The truncation column of the catalog is in units of pixels.  By
3008     default, the truncation column is considered to be in units of the
3009     radial parameters of the profile (‘--rcol’).  Read it as
3010     ‘t-unit-in-p’ for ‘truncation unit in pixels’.
3011
3012‘-f’
3013‘--mforflatpix’
3014     When making fixed value profiles (flat and circumference, see
3015     ‘‘--fcol’’), don’t use the value in the column specified by
3016     ‘‘--mcol’’ as the magnitude.  Instead use it as the exact value
3017     that all the pixels of these profiles should have.  This option is
3018     irrelevant for other types of profiles.  This option is very useful
3019     for creating masks, or labeled regions in an image.  Any integer,
3020     or floating point value can used in this column with this option,
3021     including ‘NaN’ (or ‘‘nan’’, or ‘‘NAN’’, case is irrelevant), and
3022     infinities (‘inf’, ‘-inf’, or ‘+inf’).
3023
3024     For example, with this option if you set the value in the magnitude
3025     column (‘--mcol’) to ‘NaN’, you can create an elliptical or
3026     circular mask over an image (which can be given as the argument),
3027     see *note Blank pixels::.  Another useful application of this
3028     option is to create labeled elliptical or circular apertures in an
3029     image.  To do this, set the value in the magnitude column to the
3030     label you want for this profile.  This labeled image can then be
3031     used in combination with NoiseChisel’s output (see *note
3032     NoiseChisel output::) to do aperture photometry with MakeCatalog
3033     (see *note MakeCatalog::).
3034
3035     Alternatively, if you want to mark regions of the image (for
3036     example with an elliptical circumference) and you don’t want to use
3037     NaN values (as explained above) for some technical reason, you can
3038     get the minimum or maximum value in the image (1) using Arithmetic
3039     (see *note Arithmetic::), then use that value in the magnitude
3040     column along with this option for all the profiles.
3041
3042     Please note that when using MakeProfiles on an already existing
3043     image, you have to set ‘‘--oversample=1’’.  Otherwise all the
3044     profiles will be scaled up based on the oversampling scale in your
3045     configuration files (see *note Configuration files::) unless you
3046     have accounted for oversampling in your catalog.
3047
3048‘--mcolisbrightness’
3049     The value given in the “magnitude column” (specified by ‘--mcol’,
3050     see *note MakeProfiles catalog::) must be interpreted as
3051     brightness, not magnitude.  The zero point magnitude (value to the
3052     ‘--zeropoint’ option) is ignored and the given value must have the
3053     same units as the input dataset’s pixels.
3054
3055     Recall that the total profile magnitude or brightness that is
3056     specified with in the ‘--mcol’ column of the input catalog is not
3057     an integration to infinity, but the actual sum of pixels in the
3058     profile (until the desired truncation radius).  See *note Profile
3059     magnitude:: for more on this point.
3060
3061‘--magatpeak’
3062     The magnitude column in the catalog (see *note MakeProfiles
3063     catalog::) will be used to find the brightness only for the peak
3064     profile pixel, not the full profile.  Note that this is the flux of
3065     the profile’s peak pixel in the final output of MakeProfiles.  So
3066     beware of the oversampling, see *note Oversampling::.
3067
3068     This option can be useful if you want to check a mock profile’s
3069     total magnitude at various truncation radii.  Without this option,
3070     no matter what the truncation radius is, the total magnitude will
3071     be the same as that given in the catalog.  But with this option,
3072     the total magnitude will become brighter as you increase the
3073     truncation radius.
3074
3075     In sharper profiles, sometimes the accuracy of measuring the peak
3076     profile flux is more than the overall object brightness.  In such
3077     cases, with this option, the final profile will be built such that
3078     its peak has the given magnitude, not the total profile.
3079
3080     *CAUTION:* If you want to use this option for comparing with
3081     observations, please note that MakeProfiles does not do
3082     convolution.  Unless you have de-convolved your data, your images
3083     are convolved with the instrument and atmospheric PSF, see *note
3084     PSF::.  Particularly in sharper profiles, the flux in the peak
3085     pixel is strongly decreased after convolution.  Also note that in
3086     such cases, besides de-convolution, you will have to set
3087     ‘--oversample=1’ otherwise after resampling your profile with Warp
3088     (see *note Warp::), the peak flux will be different.
3089
3090‘--customtable FITS/TXT3091     The filename of the table to use in the custom profiles (see
3092     description of ‘--fcol’ in *note MakeProfiles catalog::.  This can
3093     be a plain-text table, or FITS table, see *note Tables::, if it is
3094     a FITS table, you can use ‘--customtablehdu’ to specify which HDU
3095     should be used (described below).
3096
3097     A custom profile can have any value you want for a given radial
3098     profile (including NaN/blank values).  Each interval is defined by
3099     its minimum (inclusive) and maximum (exclusive) radius, when a
3100     pixel center falls within a radius interval, the value specified
3101     for that interval will be used.  If a pixel is not in the given
3102     intervals, a value of 0.0 will be used for that pixel.
3103
3104     The table should have 3 columns as shown below.  If the intervals
3105     are contiguous (the maximum value of the previous interval is equal
3106     to the minimum value of an interval) and the intervals all have the
3107     same size (difference between minimum and maximum values) the
3108     creation of these profiles will be fast.  However, if the intervals
3109     are not sorted and contiguous, Makeprofiles will parse the
3110     intervals from the top of the table and use the first interval that
3111     contains the pixel center.
3112
3113     Column 1:
3114          The interval’s minimum radius.
3115     Column 2:
3116          The interval’s maximum radius.
3117     Column 3:
3118          The value to be used for pixels within the given interval
3119          (including NaN/blank).
3120
3121     For example let’s assume you have the radial profile below in a
3122     file called ‘radial.txt’.  The first column is the larger interval
3123     radius (in units of pixels) and the second column is the value in
3124     that interval:
3125
3126          1    100
3127          2    90
3128          3    50
3129          4    10
3130          5    2
3131          6    0.1
3132          7    0.05
3133
3134     You can construct the table to give to ‘--customtable’ with the
3135     command below (using Gnuastro’s *note Column arithmetic::).
3136
3137          asttable radial.fits -c'arith $1 1 -' -c1,2 -ocustom.fits
3138
3139     In case the intervals are different from 1 (for example 0.5),
3140     change the ‘$1 1 -’ to ‘$1 0.5 -’.  On a side-note, Gnuastro has
3141     features to extract the radial profile of an object from the image,
3142     see *note Generate radial profile::.
3143
3144‘--customtablehdu INT/STR3145     The HDU/extension in the FITS file given to ‘--customtable’.
3146
3147‘-X INT,INT’
3148‘--shift=INT,INT’
3149     Shift all the profiles and enlarge the image along each dimension.
3150     To better understand this option, please see $n$ in *note If
3151     convolving afterwards::.  This is useful when you want to convolve
3152     the image afterwards.  If you are using an external PSF, be sure to
3153     oversample it to the same scale used for creating the mock images.
3154     If a background image is specified, any possible value to this
3155     option is ignored.
3156
3157‘-c’
3158‘--prepforconv’
3159     Shift all the profiles and enlarge the image based on half the
3160     width of the first Moffat or Gaussian profile in the catalog,
3161     considering any possible oversampling see *note If convolving
3162     afterwards::.  ‘--prepforconv’ is only checked and possibly
3163     activated if ‘--xshift’ and ‘--yshift’ are both zero (after reading
3164     the command-line and configuration files).  If a background image
3165     is specified, any possible value to this option is ignored.
3166
3167‘-z FLT’
3168‘--zeropoint=FLT’
3169     The zero point magnitude of the input.  For more on the zero point
3170     magnitude, see *note Brightness flux magnitude::.
3171
3172‘-w FLT’
3173‘--circumwidth=FLT’
3174     The width of the circumference if the profile is to be an
3175     elliptical circumference or annulus.  See the explanations for this
3176     type of profile in ‘--fcol’.
3177
3178‘-R’
3179‘--replace’
3180     Do not add the pixels of each profile over the background, or other
3181     profiles.  But replace the values.
3182
3183     By default, when two profiles overlap, the final pixel value is the
3184     sum of all the profiles that overlap on that pixel.  This is the
3185     expected situation when dealing with physical object profiles like
3186     galaxies or stars/PSF. However, when MakeProfiles is used to build
3187     integer labeled images (for example in *note Aperture
3188     photometry::), this is not the expected situation: the sum of two
3189     labels will be a new label.  With this option, the pixels are not
3190     added but the largest (maximum) value over that pixel is used.
3191     Because the maximum operator is independent of the order of values,
3192     the output is also thread-safe.
3193
3194   ---------- Footnotes ----------
3195
3196   (1) The minimum will give a better result, because the maximum can be
3197too high compared to most pixels in the image, making it harder to
3198display.
3199
3200
3201File: gnuastro.info,  Node: MakeProfiles output dataset,  Next: MakeProfiles log file,  Prev: MakeProfiles profile settings,  Up: Invoking astmkprof
3202
32038.1.4.3 MakeProfiles output dataset
3204...................................
3205
3206MakeProfiles takes an input catalog uses basic properties that are
3207defined there to build a dataset, for example a 2D image containing the
3208profiles in the catalog.  In *note MakeProfiles catalog:: and *note
3209MakeProfiles profile settings::, the catalog and profile settings were
3210discussed.  The options of this section, allow you to configure the
3211output dataset (or the canvas that will host the built profiles).
3212
3213‘-k FITS’
3214‘--background=FITS’
3215     A background image FITS file to build the profiles on.  The
3216     extension that contains the image should be specified with the
3217     ‘--backhdu’ option, see below.  When a background image is
3218     specified, it will be used to derive all the information about the
3219     output image.  Hence, the following options will be ignored:
3220     ‘--mergedsize’, ‘--oversample’, ‘--crpix’, ‘--crval’ (generally,
3221     all other WCS related parameters) and the output’s data type (see
3222     ‘--type’ in *note Input output options::).
3223
3224     The image will act like a canvas to build the profiles on: profile
3225     pixel values will be summed with the background image pixel values.
3226     With the ‘--replace’ option you can disable this behavior and
3227     replace the profile pixels with the background pixels.  If you want
3228     to use all the image information above, except for the pixel values
3229     (you want to have a blank canvas to build the profiles on, based on
3230     an input image), you can call ‘--clearcanvas’, to set all the input
3231     image’s pixels to zero before starting to build the profiles over
3232     it (this is done in memory after reading the input, so nothing will
3233     happen to your input file).
3234
3235‘-B STR/INT3236‘--backhdu=STR/INT3237     The header data unit (HDU) of the file given to ‘--background’.
3238
3239‘-C’
3240‘--clearcanvas’
3241     When an input image is specified (with the ‘--background’ option,
3242     set all its pixels to 0.0 immediately after reading it into memory.
3243     Effectively, this will allow you to use all its properties
3244     (described under the ‘--background’ option), without having to
3245     worry about the pixel values.
3246
3247     ‘--clearcanvas’ can come in handy in many situations, for example
3248     if you want to create a labeled image (segmentation map) for
3249     creating a catalog (see *note MakeCatalog::).  In other cases, you
3250     might have modeled the objects in an image and want to create them
3251     on the same frame, but without the original pixel values.
3252
3253‘-E STR/INT,FLT[,FLT,[...]]’
3254‘--kernel=STR/INT,FLT[,FLT,[...]]’
3255     Only build one kernel profile with the parameters given as the
3256     values to this option.  The different values must be separated by a
3257     comma (<,>).  The first value identifies the radial function of the
3258     profile, either through a string or through a number (see
3259     description of ‘--fcol’ in *note MakeProfiles catalog::).  Each
3260     radial profile needs a different total number of parameters: Sérsic
3261     and Moffat functions need 3 parameters: radial, Sérsic index or
3262     Moffat $\beta$, and truncation radius.  The Gaussian function needs
3263     two parameters: radial and truncation radius.  The point function
3264     doesn’t need any parameters and flat and circumference profiles
3265     just need one parameter (truncation radius).
3266
3267     The PSF or kernel is a unique (and highly constrained) type of
3268     profile: the sum of its pixels must be one, its center must be the
3269     center of the central pixel (in an image with an odd number of
3270     pixels on each side), and commonly it is circular, so its axis
3271     ratio and position angle are one and zero respectively.  Kernels
3272     are commonly necessary for various data analysis and data
3273     manipulation steps (for example see *note Convolve::, and *note
3274     NoiseChisel::.  Because of this it is inconvenient to define a
3275     catalog with one row and many zero valued columns (for all the
3276     non-necessary parameters).  Hence, with this option, it is possible
3277     to create a kernel with MakeProfiles without the need to create a
3278     catalog.  Here are some examples:
3279
3280     ‘--kernel=moffat,3,2.8,5’
3281          A Moffat kernel with FWHM of 3 pixels, $\beta=2.8$ which is
3282          truncated at 5 times the FWHM.
3283
3284     ‘--kernel=gaussian,2,3’
3285          A circular Gaussian kernel with FWHM of 2 pixels and truncated
3286          at 3 times the FWHM.
3287
3288     This option may also be used to create a 3D kernel.  To do that,
3289     two small modifications are necessary: add a ‘-3d’ (or ‘-3D’) to
3290     the profile name (for example ‘moffat-3d’) and add a number
3291     (axis-ratio along the third dimension) to the end of the parameters
3292     for all profiles except ‘point’.  The main reason behind providing
3293     an axis ratio in the third dimension is that in 3D astronomical
3294     datasets, commonly the third dimension doesn’t have the same nature
3295     (units/sampling) as the first and second.
3296
3297     For example in IFU (optical) or Radio datacubes, the first and
3298     second dimensions are commonly spatial/angular positions (like RA
3299     and Dec) but the third dimension is wavelength or frequency (in
3300     units of Angstroms for Herz).  Because of this different nature
3301     (which also affects the processing), it may be necessary for the
3302     kernel to have a different extent in that direction.
3303
3304     If the 3rd dimension axis ratio is equal to $1.0$, then the kernel
3305     will be a spheroid.  If its smaller than $1.0$, the kernel will be
3306     button-shaped: extended less in the third dimension.  However, when
3307     it islarger than $1.0$, the kernel will be bullet-shaped: extended
3308     more in the third dimension.  In the latter case, the radial
3309     parameter will correspond to the length along the 3rd dimension.
3310     For example, let’s have a look at the two examples above but in 3D:
3311
3312     ‘--kernel=moffat-3d,3,2.8,5,0.5’
3313          An ellipsoid Moffat kernel with FWHM of 3 pixels, $\beta=2.8$
3314          which is truncated at 5 times the FWHM. The ellipsoid is
3315          circular in the first two dimensions, but in the third
3316          dimension its extent is half the first two.
3317
3318     ‘--kernel=gaussian-3d,2,3,1’
3319          A spherical Gaussian kernel with FWHM of 2 pixels and
3320          truncated at 3 times the FWHM.
3321
3322     Ofcourse, if a specific kernel is needed that doesn’t fit the
3323     constraints imposed by this option, you can always use a catalog to
3324     define any arbitrary kernel.  Just call the ‘--individual’ and
3325     ‘--nomerged’ options to make sure that it is built as a separate
3326     file (individually) and no “merged” image of the input profiles is
3327     created.
3328
3329‘-x INT,INT’
3330‘--mergedsize=INT,INT’
3331     The number of pixels along each axis of the output, in FITS order.
3332     This is before over-sampling.  For example if you call MakeProfiles
3333     with ‘--mergedsize=100,150 --oversample=5’ (assuming no shift due
3334     for later convolution), then the final image size along the first
3335     axis will be 500 by 750 pixels.  Fractions are acceptable as values
3336     for each dimension, however, they must reduce to an integer, so
3337     ‘--mergedsize=150/3,300/3’ is acceptable but
3338     ‘--mergedsize=150/4,300/4’ is not.
3339
3340     When viewing a FITS image in DS9, the first FITS dimension is in
3341     the horizontal direction and the second is vertical.  As an
3342     example, the image created with the example above will have 500
3343     pixels horizontally and 750 pixels vertically.
3344
3345     If a background image is specified, this option is ignored.
3346
3347‘-s INT’
3348‘--oversample=INT’
3349     The scale to over-sample the profiles and final image.  If not an
3350     odd number, will be added by one, see *note Oversampling::.  Note
3351     that this ‘--oversample’ will remain active even if an input image
3352     is specified.  If your input catalog is based on the background
3353     image, be sure to set ‘--oversample=1’.
3354
3355‘--psfinimg’
3356     Build the possibly existing PSF profiles (Moffat or Gaussian) in
3357     the catalog into the final image.  By default they are built
3358     separately so you can convolve your images with them, thus their
3359     magnitude and positions are ignored.  With this option, they will
3360     be built in the final image like every other galaxy profile.  To
3361     have a final PSF in your image, make a point profile where you want
3362     the PSF and after convolution it will be the PSF.
3363
3364‘-i’
3365‘--individual’
3366     If this option is called, each profile is created in a separate
3367     FITS file within the same directory as the output and the row
3368     number of the profile (starting from zero) in the name.  The file
3369     for each row’s profile will be in the same directory as the final
3370     combined image of all the profiles and will have the final image’s
3371     name as a suffix.  So for example if the final combined image is
3372     named ‘./out/fromcatalog.fits’, then the first profile that will be
3373     created with this option will be named ‘./out/0_fromcatalog.fits’.
3374
3375     Since each image only has one full profile out to the truncation
3376     radius the profile is centered and so, only the sub-pixel position
3377     of the profile center is important for the outputs of this option.
3378     The output will have an odd number of pixels.  If there is no
3379     oversampling, the central pixel will contain the profile center.
3380     If the value to ‘--oversample’ is larger than unity, then the
3381     profile center is on any of the central ‘--oversample’’d pixels
3382     depending on the fractional value of the profile center.
3383
3384     If the fractional value is larger than half, it is on the bottom
3385     half of the central region.  This is due to the FITS definition of
3386     a real number position: The center of a pixel has fractional value
3387     $0.00$ so each pixel contains these fractions: .5 – .75 – .00
3388     (pixel center) – .25 – .5.
3389
3390‘-m’
3391‘--nomerged’
3392     Don’t make a merged image.  By default after making the profiles,
3393     they are added to a final image with side lengths specified by
3394     ‘--mergedsize’ if they overlap with it.
3395
3396The options below can be used to define the world coordinate system
3397(WCS) properties of the MakeProfiles outputs.  The option names are
3398deliberately chosen to be the same as the FITS standard WCS keywords.
3399See Section 8 of Pence et al [2010]
3400(https://doi.org/10.1051/0004-6361/201015362) for a short introduction
3401to WCS in the FITS standard(1).
3402
3403   If you look into the headers of a FITS image with WCS for example you
3404will see all these names but in uppercase and with numbers to represent
3405the dimensions, for example ‘CRPIX1’ and ‘PC2_1’.  You can see the FITS
3406headers with Gnuastro’s *note Fits:: program using a command like this:
3407‘$ astfits -p image.fits’.
3408
3409   If the values given to any of these options does not correspond to
3410the number of dimensions in the output dataset, then no WCS information
3411will be added.
3412
3413‘--crpix=FLT,FLT’
3414     The pixel coordinates of the WCS reference point.  Fractions are
3415     acceptable for the values of this option.
3416
3417‘--crval=FLT,FLT’
3418     The WCS coordinates of the Reference point.  Fractions are
3419     acceptable for the values of this option.
3420
3421‘--cdelt=FLT,FLT’
3422     The resolution (size of one data-unit or pixel in WCS units) of the
3423     non-oversampled dataset.  Fractions are acceptable for the values
3424     of this option.
3425
3426‘--pc=FLT,FLT,FLT,FLT’
3427     The PC matrix of the WCS rotation, see the FITS standard (link
3428     above) to better understand the PC matrix.
3429
3430‘--cunit=STR,STR’
3431     The units of each WCS axis, for example ‘deg’.  Note that these
3432     values are part of the FITS standard (link above).  MakeProfiles
3433     won’t complain if you use non-standard values, but later usage of
3434     them might cause trouble.
3435
3436‘--ctype=STR,STR’
3437     The type of each WCS axis, for example ‘RA---TAN’ and ‘DEC--TAN’.
3438     Note that these values are part of the FITS standard (link above).
3439     MakeProfiles won’t complain if you use non-standard values, but
3440     later usage of them might cause trouble.
3441
3442   ---------- Footnotes ----------
3443
3444   (1) The world coordinate standard in FITS is a very beautiful and
3445powerful concept to link/associate datasets with the outside world
3446(other datasets).  The description in the FITS standard (link above)
3447only touches the tip of the ice-burg.  To learn more please see Greisen
3448and Calabretta [2002] (https://doi.org/10.1051/0004-6361:20021326),
3449Calabretta and Greisen [2002]
3450(https://doi.org/10.1051/0004-6361:20021327), Greisen et al.  [2006]
3451(https://doi.org/10.1051/0004-6361:20053818), and Calabretta et al.
3452(http://www.atnf.csiro.au/people/mcalabre/WCS/dcs_20040422.pdf)
3453
3454
3455File: gnuastro.info,  Node: MakeProfiles log file,  Prev: MakeProfiles output dataset,  Up: Invoking astmkprof
3456
34578.1.4.4 MakeProfiles log file
3458.............................
3459
3460Besides the final merged dataset of all the profiles, or the individual
3461datasets (see *note MakeProfiles output dataset::), if the ‘--log’
3462option is called MakeProfiles will also create a log file in the current
3463directory (where you run MockProfiles).  See *note Common options:: for
3464a full description of ‘--log’ and other options that are shared between
3465all Gnuastro programs.  The values for each column are explained in the
3466first few commented lines of the log file (starting with ‘#’ character).
3467Here is a more complete description.
3468
3469   • An ID (row number of profile in input catalog).
3470
3471   • The total magnitude of the profile in the output dataset.  When the
3472     profile does not completely overlap with the output dataset, this
3473     will be different from your input magnitude.
3474
3475   • The number of pixels (in the oversampled image) which used Monte
3476     Carlo integration and not the central pixel value, see *note
3477     Sampling from a function::.
3478
3479   • The fraction of flux in the Monte Carlo integrated pixels.
3480
3481   • If an individual image was created, this column will have a value
3482     of ‘1’, otherwise it will have a value of ‘0’.
3483
3484
3485File: gnuastro.info,  Node: MakeNoise,  Prev: MakeProfiles,  Up: Modeling and fittings
3486
34878.2 MakeNoise
3488=============
3489
3490Real data are always buried in noise, therefore to finalize a simulation
3491of real data (for example to test our observational algorithms) it is
3492essential to add noise to the mock profiles created with MakeProfiles,
3493see *note MakeProfiles::.  Below, the general principles and concepts to
3494help understand how noise is quantified is discussed.  MakeNoise options
3495and argument are then discussed in *note Invoking astmknoise::.
3496
3497* Menu:
3498
3499* Noise basics::                Noise concepts and definitions.
3500* Invoking astmknoise::         Options and arguments to MakeNoise.
3501
3502
3503File: gnuastro.info,  Node: Noise basics,  Next: Invoking astmknoise,  Prev: MakeNoise,  Up: MakeNoise
3504
35058.2.1 Noise basics
3506------------------
3507
3508Deep astronomical images, like those used in extragalactic studies,
3509seriously suffer from noise in the data.  Generally speaking, the
3510sources of noise in an astronomical image are photon counting noise and
3511Instrumental noise which are discussed in *note Photon counting noise::
3512and *note Instrumental noise::.  This review finishes with *note
3513Generating random numbers:: which is a short introduction on how random
3514numbers are generated.  We will see that while software random number
3515generators are not perfect, they allow us to obtain a reproducible
3516series of random numbers through setting the random number generator
3517function and seed value.  Therefore in this section, we’ll also discuss
3518how you can set these two parameters in Gnuastro’s programs (including
3519MakeNoise).
3520
3521* Menu:
3522
3523* Photon counting noise::       Poisson noise
3524* Instrumental noise::          Readout, dark current and other sources.
3525* Final noised pixel value::    How the final noised value is calculated.
3526* Generating random numbers::   How random numbers are generated.
3527
3528
3529File: gnuastro.info,  Node: Photon counting noise,  Next: Instrumental noise,  Prev: Noise basics,  Up: Noise basics
3530
35318.2.1.1 Photon counting noise
3532.............................
3533
3534With the very accurate electronics used in today’s detectors, photon
3535counting noise(1) is the most significant source of uncertainty in most
3536datasets.  To understand this noise (error in counting), we need to take
3537a closer look at how a distribution produced by counting can be modeled
3538as a parametric function.
3539
3540   Counting is an inherently discrete operation, which can only produce
3541positive (including zero) integer outputs.  For example we can’t count
3542$3.2$ or $-2$ of anything.  We only count $0$, $1$, $2$, $3$ and so on.
3543The distribution of values, as a result of counting efforts is formally
3544known as the Poisson distribution
3545(https://en.wikipedia.org/wiki/Poisson_distribution).  It is associated
3546to Siméon Denis Poisson, because he discussed it while working on the
3547number of wrongful convictions in court cases in his 1837 book(2).
3548
3549   Let’s take $\lambda$ to represent the expected mean count of
3550something.  Furthermore, let’s take $k$ to represent the result of one
3551particular counting attempt.  The probability density function of
3552getting $k$ counts (in each attempt, given the expected/mean count of
3553$\lambda$) can be written as:
3554
3555$$f(k)={\lambda^k \over k!} e^{-\lambda},\quad k\in {0, 1, 2, 3, \dots }$$
3556
3557   Because the Poisson distribution is only applicable to positive
3558values (note the factorial operator, which only applies to non-negative
3559integers), naturally it is very skewed when $\lambda$ is near zero.  One
3560qualitative way to understand this behavior is that there simply aren’t
3561enough integers smaller than $\lambda$, than integers that are larger
3562than it.  Therefore to accommodate all possibilities/counts, it has to
3563be strongly skewed when $\lambda$ is small.
3564
3565   As $\lambda$ becomes larger, the distribution becomes more and more
3566symmetric.  A very useful property of the Poisson distribution is that
3567the mean value is also its variance.  When $\lambda$ is very large, say
3568$\lambda>1000$, then the Normal (Gaussian) distribution
3569(https://en.wikipedia.org/wiki/Normal_distribution), is an excellent
3570approximation of the Poisson distribution with mean $\mu=\lambda$ and
3571standard deviation $\sigma=\sqrt{\lambda}$.  In other words, a Poisson
3572distribution (with a sufficiently large $\lambda$) is simply a Gaussian
3573that only has one free parameter ($\mu=\lambda$ and
3574$\sigma=\sqrt{\lambda}$), instead of the two parameters (independent
3575$\mu$ and $\sigma$) that it originally has.
3576
3577   In real situations, the photons/flux from our targets are added to a
3578certain background flux (observationally, the _Sky_ value).  The Sky
3579value is defined to be the average flux of a region in the dataset with
3580no targets.  Its physical origin can be the brightness of the atmosphere
3581(for ground-based instruments), possible stray light within the imaging
3582instrument, the average flux of undetected targets, etc.  The Sky value
3583is thus an ideal definition, because in real datasets, what lies deep in
3584the noise (far lower than the detection limit) is never known(3).  To
3585account for all of these, the sky value is defined to be the average
3586count/value of the undetected regions in the image.  In a mock
3587image/dataset, we have the luxury of setting the background (Sky) value.
3588
3589   In each element of the dataset (pixel in an image), the flux is the
3590sum of contributions from various sources (after convolution by the PSF,
3591see *note PSF::).  Let’s name the convolved sum of possibly overlapping
3592objects, $I_{nn}$.  $nn$ representing ‘no noise’.  For now, let’s assume
3593the background ($B$) is constant and sufficiently high for the Poisson
3594distribution to be approximated by a Gaussian.  Then the flux after
3595adding noise is a random value taken from a Gaussian distribution with
3596the following mean ($\mu$) and standard deviation ($\sigma$):
3597
3598            $$\mu=B+I_{nn}, \quad \sigma=\sqrt{B+I_{nn}}$$
3599
3600   Since this type of noise is inherent in the objects we study, it is
3601usually measured on the same scale as the astronomical objects, namely
3602the magnitude system, see *note Brightness flux magnitude::.  It is then
3603internally converted to the flux scale for further processing.
3604
3605   ---------- Footnotes ----------
3606
3607   (1) In practice, we are actually counting the electrons that are
3608produced by each photon, not the actual photons.
3609
3610   (2) [From Wikipedia] Poisson’s result was also derived in a previous
3611study by Abraham de Moivre in 1711.  Therefore some people suggest it
3612should rightly be called the de Moivre distribution.
3613
3614   (3) In a real image, a relatively large number of very faint objects
3615can been fully buried in the noise and never detected.  These undetected
3616objects will bias the background measurement to slightly larger values.
3617Our best approximation is thus to simply assume they are uniform, and
3618consider their average effect.  See Figure 1 (a.1 and a.2) and Section
36192.2 in Akhlaghi and Ichikawa [2015] (https://arxiv.org/abs/1505.01664).
3620
3621
3622File: gnuastro.info,  Node: Instrumental noise,  Next: Final noised pixel value,  Prev: Photon counting noise,  Up: Noise basics
3623
36248.2.1.2 Instrumental noise
3625..........................
3626
3627While taking images with a camera, a dark current is fed to the pixels,
3628the variation of the value of this dark current over the pixels, also
3629adds to the final image noise.  Another source of noise is the readout
3630noise that is produced by the electronics in the detector.
3631Specifically, the parts that attempt to digitize the voltage produced by
3632the photo-electrons in the analog to digital converter.  With the
3633current generation of instruments, this source of noise is not as
3634significant as the noise due to the background Sky discussed in *note
3635Photon counting noise::.
3636
3637   Let $C$ represent the combined standard deviation of all these
3638instrumental sources of noise.  When only this source of noise is
3639present, the noised pixel value would be a random value chosen from a
3640Gaussian distribution with
3641
3642            $$\mu=I_{nn}, \quad \sigma=\sqrt{C^2+I_{nn}}$$
3643
3644   This type of noise is independent of the signal in the dataset, it is
3645only determined by the instrument.  So the flux scale (and not magnitude
3646scale) is most commonly used for this type of noise.  In practice, this
3647value is usually reported in analog-to-digital units or ADUs, not flux
3648or electron counts.  The gain value of the device can be used to convert
3649between these two, see *note Brightness flux magnitude::.
3650
3651
3652File: gnuastro.info,  Node: Final noised pixel value,  Next: Generating random numbers,  Prev: Instrumental noise,  Up: Noise basics
3653
36548.2.1.3 Final noised pixel value
3655................................
3656
3657Based on the discussions in *note Photon counting noise:: and *note
3658Instrumental noise::, depending on the values you specify for $B$ and
3659$C$ from the above, the final noised value for each pixel is a random
3660value chosen from a Gaussian distribution with
3661
3662          $$\mu=B+I_{nn}, \quad \sigma=\sqrt{C^2+B+I_{nn}}$$
3663
3664
3665File: gnuastro.info,  Node: Generating random numbers,  Prev: Final noised pixel value,  Up: Noise basics
3666
36678.2.1.4 Generating random numbers
3668.................................
3669
3670As discussed above, to generate noise we need to make random samples of
3671a particular distribution.  So it is important to understand some
3672general concepts regarding the generation of random numbers.  For a very
3673complete and nice introduction we strongly advise reading Donald Knuth’s
3674“The art of computer programming”, volume 2, chapter 3(1).  Quoting from
3675the GNU Scientific Library manual, “If you don’t own it, you should stop
3676reading right now, run to the nearest bookstore, and buy it”(2)!
3677
3678   Using only software, we can only produce what is called a
3679psuedo-random sequence of numbers.  A true random number generator is a
3680hardware (let’s assume we have made sure it has no systematic biases),
3681for example throwing dice or flipping coins (which have remained from
3682the ancient times).  More modern hardware methods use atmospheric noise,
3683thermal noise or other types of external electromagnetic or quantum
3684phenomena.  All pseudo-random number generators (software) require a
3685seed to be the basis of the generation.  The advantage of having a seed
3686is that if you specify the same seed for multiple runs, you will get an
3687identical sequence of random numbers which allows you to reproduce the
3688same final noised image.
3689
3690   The programs in GNU Astronomy Utilities (for example MakeNoise or
3691MakeProfiles) use the GNU Scientific Library (GSL) to generate random
3692numbers.  GSL allows the user to set the random number generator through
3693environment variables, see *note Installation directory:: for an
3694introduction to environment variables.  In the chapter titled “Random
3695Number Generation” they have fully explained the various random number
3696generators that are available (there are a lot of them!).  Through the
3697two environment variables ‘GSL_RNG_TYPE’ and ‘GSL_RNG_SEED’ you can
3698specify the generator and its seed respectively.
3699
3700   If you don’t specify a value for ‘GSL_RNG_TYPE’, GSL will use its
3701default random number generator type.  The default type is sufficient
3702for most general applications.  If no value is given for the
3703‘GSL_RNG_SEED’ environment variable and you have asked Gnuastro to read
3704the seed from the environment (through the ‘--envseed’ option), then GSL
3705will use the default value of each generator to give identical outputs.
3706If you don’t explicitly tell Gnuastro programs to read the seed value
3707from the environment variable, then they will use the system time
3708(accurate to within a microsecond) to generate (apparently random)
3709seeds.  In this manner, every time you run the program, you will get a
3710different random number distribution.
3711
3712   There are two ways you can specify values for these environment
3713variables.  You can call them on the same command-line for example:
3714
3715     $ GSL_RNG_TYPE="taus" GSL_RNG_SEED=345 astmknoise input.fits
3716
3717In this manner the values will only be used for this particular
3718execution of MakeNoise.  Alternatively, you can define them for the full
3719period of your terminal session or script length, using the shell’s
3720‘export’ command with the two separate commands below (for a script
3721remove the ‘$’ signs):
3722
3723     $ export GSL_RNG_TYPE="taus"
3724     $ export GSL_RNG_SEED=345
3725
3726The subsequent programs which use GSL’s random number generators will
3727hence forth use these values in this session of the terminal you are
3728running or while executing this script.  In case you want to set fixed
3729values for these parameters every time you use the GSL random number
3730generator, you can add these two lines to your ‘.bashrc’ startup
3731script(3), see *note Installation directory::.
3732
3733   *IMPORTANT NOTE:* If the two environment variables ‘GSL_RNG_TYPE’ and
3734‘GSL_RNG_SEED’ are defined, GSL will report them by default, even if you
3735don’t use the ‘--envseed’ option.  For example, see this call to
3736MakeProfiles:
3737
3738     $ export GSL_RNG_TYPE=taus
3739     $ export GSL_RNG_SEED=345
3740     $ astmkprof -s1 --kernel=gaussian,2,5
3741     GSL_RNG_TYPE=taus
3742     GSL_RNG_SEED=345
3743     MakeProfiles V.VV started on DDD MMM DDD HH:MM:SS YYYY
3744       - Building one gaussian kernel
3745       - Random number generator (RNG) type: taus
3746       - Basic RNG seed: 1618960836
3747       ---- ./kernel.fits created.
3748       -- Output: ./kernel.fits
3749     MakeProfiles finished in 0.068945 seconds
3750
3751The first two output lines (showing the names and values of the GSL
3752environment variables) are printed by GSL before MakeProfiles actually
3753starts generating random numbers.  Gnuastro’s programs will report the
3754actual values they use independently (after the name of the program),
3755you should check them for the final values used, not GSL’s printed
3756values.  In the example above, did you notice how the random number
3757generator seed above is different between GSL and MakeProfiles?
3758However, if ‘--envseed’ was given, both printed seeds would be the same.
3759
3760   ---------- Footnotes ----------
3761
3762   (1) Knuth, Donald.  1998.  The art of computer programming.
3763Addison–Wesley.  ISBN 0-201-89684-2
3764
3765   (2) For students, running to the library might be more affordable!
3766
3767   (3) Don’t forget that if you are going to give your scripts (that use
3768the GSL random number generator) to others you have to make sure you
3769also tell them to set these environment variable separately.  So for
3770scripts, it is best to keep all such variable definitions within the
3771script, even if they are within your ‘.bashrc’.
3772
3773
3774File: gnuastro.info,  Node: Invoking astmknoise,  Prev: Noise basics,  Up: MakeNoise
3775
37768.2.2 Invoking MakeNoise
3777------------------------
3778
3779MakeNoise will add noise to an existing image.  The executable name is
3780‘astmknoise’ with the following general template
3781
3782     $ astmknoise [OPTION ...] InputImage.fits
3783
3784One line examples:
3785
3786     ## Add noise with a standard deviation of 100 to image.
3787     ## (this is independent of the pixel value: not Poission noise)
3788     $ astmknoise --sigma=100 image.fits
3789
3790     ## Add noise to input image assuming a background magnitude (with
3791     ## zero point magnitude of 0) and a certain instrumental noise:
3792     $ astmknoise --background=-10 -z0 --instrumental=20 mockimage.fits
3793
3794If actual processing is to be done, the input image is a mandatory
3795argument.  The full list of options common to all the programs in
3796Gnuastro can be seen in *note Common options::.  The type (see *note
3797Numeric data types::) of the output can be specified with the ‘--type’
3798option, see *note Input output options::.  The header of the output FITS
3799file keeps all the parameters that were influential in making it.  This
3800is done for future reproducibility.
3801
3802‘-b FLT’
3803‘--background=FLT’
3804     The background value (per pixel) that will be added to each pixel
3805     value (internally) to estimate Poisson noise, see *note Photon
3806     counting noise::.  By default the units of this value are assumed
3807     to be in magnitudes, hence a ‘--zeropoint’ is also necessary.  But
3808     if the background is in units of brightness, you need add
3809     ‘--bgisbrightness’, see *note Brightness flux magnitude::
3810
3811     Internally, the value given to this option will be converted to
3812     brightness ($b$, when ‘--bgisbrightness’ is called, the value will
3813     be used directly).  Assuming the pixel value is $p$, the random
3814     value for that pixel will be taken from a Gaussian distribution
3815     with mean of $p+b$ and standard deviation of $\sqrt{p+b}$.  With
3816     this option, the noise will therefore be dependent on the pixel
3817     values: according to the Poission noise model, as the pixel value
3818     becomes larger, its noise will also become larger.  This is thus a
3819     realistic way to model noise, see *note Photon counting noise::.
3820
3821‘-B’
3822‘--bgisbrightness’
3823     The value given to ‘--background’ should be interpreted as
3824     brightness, not as a magnitude.
3825
3826‘-z FLT’
3827‘--zeropoint=FLT’
3828     The zero point magnitude used to convert the value of
3829     ‘--background’ (in units of magnitude) to flux, see *note
3830     Brightness flux magnitude::.
3831
3832‘-i FLT’
3833‘--instrumental=FLT’
3834     The instrumental noise which is in units of flux, see *note
3835     Instrumental noise::.
3836
3837‘-s FLT’
3838‘--sigma=FLT’
3839     The total noise sigma in the same units as the pixel values.  With
3840     this option, the ‘--background’, ‘--zeropoint’ and ‘--instrumental’
3841     will be ignored.  With this option, the noise will be independent
3842     of the pixel values (which is not realistic, see *note Photon
3843     counting noise::).  Hence it is only useful if you are working on
3844     low surface brightness regions where the change in pixel value (and
3845     thus real noise) is insignificant.
3846
3847     Generally, *usage of this option is discouraged* unless you
3848     understand the risks of not simulating real noise.  This is because
3849     with this option, you will not get Poisson noise (the common noise
3850     model for astronomical imaging), where the noise varies based on
3851     pixel value.  Use ‘--background’ for adding Poission noise.
3852
3853‘-e’
3854‘--envseed’
3855     Use the ‘GSL_RNG_SEED’ environment variable for the seed used in
3856     the random number generator, see *note Generating random numbers::.
3857     With this option, the output image noise is always going to be
3858     identical (or reproducible).
3859
3860‘-d’
3861‘--doubletype’
3862     Save the output in the double precision floating point format that
3863     was used internally.  This option will be most useful if the input
3864     images were of integer types.
3865
3866
3867File: gnuastro.info,  Node: High-level calculations,  Next: Installed scripts,  Prev: Modeling and fittings,  Up: Top
3868
38699 High-level calculations
3870*************************
3871
3872After the reduction of raw data (for example with the programs in *note
3873Data manipulation::) you will have reduced images/data ready for
3874processing/analyzing (for example with the programs in *note Data
3875analysis::).  But the processed/analyzed data (or catalogs) are still
3876not enough to derive any scientific result.  Even higher-level analysis
3877is still needed to convert the observed magnitudes, sizes or volumes
3878into physical quantities that we associate with each catalog entry or
3879detected object which is the purpose of the tools in this section.
3880
3881* Menu:
3882
3883* CosmicCalculator::            Calculate cosmological variables
3884
3885
3886File: gnuastro.info,  Node: CosmicCalculator,  Prev: High-level calculations,  Up: High-level calculations
3887
38889.1 CosmicCalculator
3889====================
3890
3891To derive higher-level information regarding our sources in
3892extra-galactic astronomy, cosmological calculations are necessary.  In
3893Gnuastro, CosmicCalculator is in charge of such calculations.  Before
3894discussing how CosmicCalculator is called and operates (in *note
3895Invoking astcosmiccal::), it is important to provide a rough but mostly
3896self sufficient review of the basics and the equations used in the
3897analysis.  In *note Distance on a 2D curved space:: the basic idea of
3898understanding distances in a curved and expanding 2D universe (which we
3899can visualize) are reviewed.  Having solidified the concepts there, in
3900*note Extending distance concepts to 3D::, the formalism is extended to
3901the 3D universe we are trying to study in our research.
3902
3903   The focus here is obtaining a physical insight into these equations
3904(mainly for the use in real observational studies).  There are many
3905books thoroughly deriving and proving all the equations with all
3906possible initial conditions and assumptions for any abstract universe,
3907interested readers can study those books.
3908
3909* Menu:
3910
3911* Distance on a 2D curved space::  Distances in 2D for simplicity
3912* Extending distance concepts to 3D::  Going to 3D (our real universe).
3913* Invoking astcosmiccal::       How to run CosmicCalculator
3914
3915
3916File: gnuastro.info,  Node: Distance on a 2D curved space,  Next: Extending distance concepts to 3D,  Prev: CosmicCalculator,  Up: CosmicCalculator
3917
39189.1.1 Distance on a 2D curved space
3919-----------------------------------
3920
3921The observations to date (for example the Planck 2015 results), have not
3922measured(1) the presence of significant curvature in the universe.
3923However to be generic (and allow its measurement if it does in fact
3924exist), it is very important to create a framework that allows non-zero
3925uniform curvature.  However, this section is not intended to be a fully
3926thorough and mathematically complete derivation of these concepts.
3927There are many references available for such reviews that go deep into
3928the abstract mathematical proofs.  The emphasis here is on visualization
3929of the concepts for a beginner.
3930
3931   As 3D beings, it is difficult for us to mentally create (visualize) a
3932picture of the curvature of a 3D volume.  Hence, here we will assume a
39332D surface/space and discuss distances on that 2D surface when it is
3934flat and when it is curved.  Once the concepts have been
3935created/visualized here, we will extend them, in *note Extending
3936distance concepts to 3D::, to a real 3D spatial _slice_ of the Universe
3937we live in and hope to study.
3938
3939   To be more understandable (actively discuss from an observer’s point
3940of view) let’s assume there’s an imaginary 2D creature living on the 2D
3941space (which _might_ be curved in 3D). Here, we will be working with
3942this creature in its efforts to analyze distances in its 2D universe.
3943The start of the analysis might seem too mundane, but since it is
3944difficult to imagine a 3D curved space, it is important to review all
3945the very basic concepts thoroughly for an easy transition to a universe
3946that is more difficult to visualize (a curved 3D space embedded in 4D).
3947
3948   To start, let’s assume a static (not expanding or shrinking), flat 2D
3949surface similar to *note Figure 9.1: flatplane. and that the 2D creature
3950is observing its universe from point $A$.  One of the most basic ways to
3951parameterize this space is through the Cartesian coordinates ($x$, $y$).
3952In *note Figure 9.1: flatplane, the basic axes of these two coordinates
3953are plotted.  An infinitesimal change in the direction of each axis is
3954written as $dx$ and $dy$.  For each point, the infinitesimal changes are
3955parallel with the respective axes and are not shown for clarity.
3956Another very useful way of parameterizing this space is through polar
3957coordinates.  For each point, we define a radius ($r$) and angle
3958($\phi$) from a fixed (but arbitrary) reference axis.  In *note Figure
39599.1: flatplane. the infinitesimal changes for each polar coordinate are
3960plotted for a random point and a dashed circle is shown for all points
3961with the same radius.
3962
3963�[image src="gnuastro-figures/flatplane.png" text="../gnuastro-figures//flatplane.eps"�]
3964
3965
3966Figure 9.1: Two dimensional Cartesian and polar coordinates on a flat
3967plane.
3968
3969   Assuming an object is placed at a certain position, which can be
3970parameterized as $(x,y)$, or $(r,\phi)$, a general infinitesimal change
3971in its position will place it in the coordinates $(x+dx,y+dy)$ and
3972$(r+dr,\phi+d\phi)$.  The distance (on the flat 2D surface) that is
3973covered by this infinitesimal change in the static universe ($ds_s$, the
3974subscript signifies the static nature of this universe) can be written
3975as:
3976
3977                  $$ds_s=dx^2+dy^2=dr^2+r^2d\phi^2$$
3978
3979   The main question is this: how can the 2D creature incorporate the
3980(possible) curvature in its universe when it’s calculating distances?
3981The universe that it lives in might equally be a curved surface like
3982*note Figure 9.2: sphereandplane.  The answer to this question but for a
39833D being (us) is the whole purpose to this discussion.  Here, we want to
3984give the 2D creature (and later, ourselves) the tools to measure
3985distances if the space (that hosts the objects) is curved.
3986
3987   *note Figure 9.2: sphereandplane. assumes a spherical shell with
3988radius $R$ as the curved 2D plane for simplicity.  The 2D plane is
3989tangent to the spherical shell and only touches it at $A$.  This idea
3990will be generalized later.  The first step in measuring the distance in
3991a curved space is to imagine a third dimension along the $z$ axis as
3992shown in *note Figure 9.2: sphereandplane.  For simplicity, the $z$ axis
3993is assumed to pass through the center of the spherical shell.  Our
3994imaginary 2D creature cannot visualize the third dimension or a curved
39952D surface within it, so the remainder of this discussion is purely
3996abstract for it (similar to us having difficulty in visualizing a 3D
3997curved space in 4D). But since we are 3D creatures, we have the
3998advantage of visualizing the following steps.  Fortunately the 2D
3999creature is already familiar with our mathematical constructs, so it can
4000follow our reasoning.
4001
4002   With the third axis added, a generic infinitesimal change over _the
4003full_ 3D space corresponds to the distance:
4004
4005            $$ds_s^2=dx^2+dy^2+dz^2=dr^2+r^2d\phi^2+dz^2.$$
4006
4007�[image src="gnuastro-figures/sphereandplane.png" text="../gnuastro-figures//sphereandplane.eps"�]
4008
4009
4010Figure 9.2: 2D spherical shell (centered on $O$) and flat plane (light
4011gray) tangent to it at point $A$.
4012
4013   It is very important to recognize that this change of distance is for
4014_any_ point in the 3D space, not just those changes that occur on the 2D
4015spherical shell of *note Figure 9.2: sphereandplane.  Recall that our 2D
4016friend can only do measurements on the 2D surfaces, not the full 3D
4017space.  So we have to constrain this general change to any change on the
40182D spherical shell.  To do that, let’s look at the arbitrary point $P$
4019on the 2D spherical shell.  Its image ($P'$) on the flat plain is also
4020displayed.  From the dark gray triangle, we see that
4021
4022$$\sin\theta={r\over R},\quad\cos\theta={R-z\over R}.$$These relations allow the 2D creature to find the value of $z$ (an abstract dimension for it) as a function of r (distance on a flat 2D plane, which it can visualize) and thus eliminate $z$.
4023   From $\sin^2\theta+\cos^2\theta=1$, we get $z^2-2Rz+r^2=0$ and
4024solving for $z$, we find:
4025
4026           $$z=R\left(1\pm\sqrt{1-{r^2\over R^2}}\right).$$
4027
4028   The $\pm$ can be understood from *note Figure 9.2: sphereandplane.:
4029For each $r$, there are two points on the sphere, one in the upper
4030hemisphere and one in the lower hemisphere.  An infinitesimal change in
4031$r$, will create the following infinitesimal change in $z$:
4032
4033                    $$dz={\mp r\over R}\left(1\over
4034\sqrt{1-{r^2/R^2}}\right)dr.$$Using the positive signed equation instead of $dz$ in the $ds_s^2$ equation above, we get:
4035
4036             $$ds_s^2={dr^2\over 1-r^2/R^2}+r^2d\phi^2.$$
4037
4038   The derivation above was done for a spherical shell of radius $R$ as
4039a curved 2D surface.  To generalize it to any surface, we can define
4040$K=1/R^2$ as the curvature parameter.  Then the general infinitesimal
4041change in a static universe can be written as:
4042
4043               $$ds_s^2={dr^2\over 1-Kr^2}+r^2d\phi^2.$$
4044
4045   Therefore, when $K>0$ (and curvature is the same everywhere), we have
4046a finite universe, where $r$ cannot become larger than $R$ as in *note
4047Figure 9.2: sphereandplane.  When $K=0$, we have a flat plane (*note
4048Figure 9.1: flatplane.) and a negative $K$ will correspond to an
4049imaginary $R$.  The latter two cases may be infinite in area (which is
4050not a simple concept, but mathematically can be modeled with $r$
4051extending infinitely), or finite-area (like a cylinder is flat
4052everywhere with $ds_s^2={dx^2 + dy^2}$, but finite in one direction in
4053size).
4054
4055   A very important issue that can be discussed now (while we are still
4056in 2D and can actually visualize things) is that $\overrightarrow{r}$ is
4057tangent to the curved space at the observer’s position.  In other words,
4058it is on the gray flat surface of *note Figure 9.2: sphereandplane, even
4059when the universe if curved: $\overrightarrow{r}=P'-A$.  Therefore for
4060the point $P$ on a curved space, the raw coordinate $r$ is the distance
4061to $P'$, not $P$.  The distance to the point $P$ (at a specific
4062coordinate $r$ on the flat plane) over the curved surface (thick line in
4063*note Figure 9.2: sphereandplane.) is called the _proper distance_ and
4064is displayed with $l$.  For the specific example of *note Figure 9.2:
4065sphereandplane, the proper distance can be calculated with: $l=R\theta$
4066($\theta$ is in radians).  Using the $\sin\theta$ relation found above,
4067we can find $l$ as a function of $r$:
4068
4069    $$\theta=\sin^{-1}\left({r\over R}\right)\quad\rightarrow\quad
4070               l(r)=R\sin^{-1}\left({r\over R}\right)$$
4071
4072   $R$ is just an arbitrary constant and can be directly found from $K$,
4073so for cleaner equations, it is common practice to set $R=1$, which
4074gives: $l(r)=\sin^{-1}r$.  Also note that when $R=1$, then $l=\theta$.
4075Generally, depending on the curvature, in a _static_ universe the proper
4076distance can be written as a function of the coordinate $r$ as (from now
4077on we are assuming $R=1$):
4078
4079               $$l(r)=\sin^{-1}(r)\quad(K>0),\quad\quad
4080    l(r)=r\quad(K=0),\quad\quad l(r)=\sinh^{-1}(r)\quad(K<0).$$With
4081   $l$, the infinitesimal change of distance can be written in a more
4082simpler and abstract form of
4083
4084                      $$ds_s^2=dl^2+r^2d\phi^2.$$
4085
4086   Until now, we had assumed a static universe (not changing with time).
4087But our observations so far appear to indicate that the universe is
4088expanding (it isn’t static).  Since there is no reason to expect the
4089observed expansion is unique to our particular position of the universe,
4090we expect the universe to be expanding at all points with the same rate
4091at the same time.  Therefore, to add a time dependence to our distance
4092measurements, we can include a multiplicative scaling factor, which is a
4093function of time: $a(t)$.  The functional form of $a(t)$ comes from the
4094cosmology, the physics we assume for it: general relativity, and the
4095choice of whether the universe is uniform (‘homogeneous’) in density and
4096curvature or inhomogeneous.  In this section, the functional form of
4097$a(t)$ is irrelevant, so we can avoid these issues.
4098
4099   With this scaling factor, the proper distance will also depend on
4100time.  As the universe expands, the distance between two given points
4101will shift to larger values.  We thus define a distance measure, or
4102coordinate, that is independent of time and thus doesn’t ‘move’.  We
4103call it the _comoving distance_ and display with $\chi$ such that:
4104$l(r,t)=\chi(r)a(t)$.  We have therefore, shifted the $r$ dependence of
4105the proper distance we derived above for a static universe to the
4106comoving distance:
4107
4108              $$\chi(r)=\sin^{-1}(r)\quad(K>0),\quad\quad
4109   \chi(r)=r\quad(K=0),\quad\quad \chi(r)=\sinh^{-1}(r)\quad(K<0).$$
4110
4111   Therefore, $\chi(r)$ is the proper distance to an object at a
4112specific reference time: $t=t_r$ (the $r$ subscript signifies
4113“reference”) when $a(t_r)=1$.  At any arbitrary moment ($t\neq{t_r}$)
4114before or after $t_r$, the proper distance to the object can be scaled
4115with $a(t)$.
4116
4117   Measuring the change of distance in a time-dependent (expanding)
4118universe only makes sense if we can add up space and time(2).  But we
4119can only add bits of space and time together if we measure them in the
4120same units: with a conversion constant (similar to how 1000 is used to
4121convert a kilometer into meters).  Experimentally, we find strong
4122support for the hypothesis that this conversion constant is the speed of
4123light (or gravitational waves(3)) in a vacuum.  This speed is postulated
4124to be constant(4) and is almost always written as $c$.  We can thus
4125parameterize the change in distance on an expanding 2D surface as
4126
4127  $$ds^2=c^2dt^2-a^2(t)ds_s^2 = c^2dt^2-a^2(t)(d\chi^2+r^2d\phi^2).$$
4128
4129   ---------- Footnotes ----------
4130
4131   (1) The observations are interpreted under the assumption of uniform
4132curvature.  For a relativistic alternative to dark energy (and maybe
4133also some part of dark matter), non-uniform curvature may be even be
4134more critical, but that is beyond the scope of this brief explanation.
4135
4136   (2) In other words, making our space-time consistent with Minkowski
4137space-time geometry.  In this geometry, different observers at a given
4138point (event) in space-time split up space-time into ‘space’ and ‘time’
4139in different ways, just like people at the same spatial position can
4140make different choices of splitting up a map into ‘left–right’ and
4141‘up–down’.  This model is well supported by twentieth and twenty-first
4142century observations.
4143
4144   (3) The speed of gravitational waves was recently found to be very
4145similar to that of light in vacuum, see arXiv:1710.05834
4146(https://arxiv.org/abs/1710.05834).
4147
4148   (4) In _natural units_, speed is measured in units of the speed of
4149light in vacuum.
4150
4151
4152File: gnuastro.info,  Node: Extending distance concepts to 3D,  Next: Invoking astcosmiccal,  Prev: Distance on a 2D curved space,  Up: CosmicCalculator
4153
41549.1.2 Extending distance concepts to 3D
4155---------------------------------------
4156
4157The concepts of *note Distance on a 2D curved space:: are here extended
4158to a 3D space that _might_ be curved.  We can start with the generic
4159infinitesimal distance in a static 3D universe, but this time in
4160spherical coordinates instead of polar coordinates.  $\theta$ is shown
4161in *note Figure 9.2: sphereandplane, but here we are 3D beings,
4162positioned on $O$ (the center of the sphere) and the point $O$ is
4163tangent to a 4D-sphere.  In our 3D space, a generic infinitesimal
4164displacement will correspond to the following distance in spherical
4165coordinates:
4166
4167 $$ds_s^2=dx^2+dy^2+dz^2=dr^2+r^2(d\theta^2+\sin^2{\theta}d\phi^2).$$
4168
4169   Like the 2D creature before, we now have to assume an abstract
4170dimension which we cannot visualize easily.  Let’s call the fourth
4171dimension $w$, then the general change in coordinates in the _full_ four
4172dimensional space will be:
4173
4174      $$ds_s^2=dr^2+r^2(d\theta^2+\sin^2{\theta}d\phi^2)+dw^2.$$
4175
4176But we can only work on a 3D curved space, so following exactly the same
4177steps and conventions as our 2D friend, we arrive at:
4178
4179  $$ds_s^2={dr^2\over 1-Kr^2}+r^2(d\theta^2+\sin^2{\theta}d\phi^2).$$
4180
4181In a non-static universe (with a scale factor a(t)), the distance can be
4182written as:
4183
4184$$ds^2=c^2dt^2-a^2(t)[d\chi^2+r^2(d\theta^2+\sin^2{\theta}d\phi^2)].$$
4185
4186
4187File: gnuastro.info,  Node: Invoking astcosmiccal,  Prev: Extending distance concepts to 3D,  Up: CosmicCalculator
4188
41899.1.3 Invoking CosmicCalculator
4190-------------------------------
4191
4192CosmicCalculator will calculate cosmological variables based on the
4193input parameters.  The executable name is ‘astcosmiccal’ with the
4194following general template
4195
4196     $ astcosmiccal [OPTION...] ...
4197
4198One line examples:
4199
4200     ## Print basic cosmological properties at redshift 2.5:
4201     $ astcosmiccal -z2.5
4202
4203     ## Only print Comoving volume over 4pi stradian to z (Mpc^3):
4204     $ astcosmiccal --redshift=0.8 --volume
4205
4206     ## Print redshift and age of universe when Lyman-alpha line is
4207     ## at 6000 angstrom (another way to specify redshift).
4208     $ astcosmiccal --obsline=lyalpha,6000 --age
4209
4210     ## Print luminosity distance, angular diameter distance and age
4211     ## of universe in one row at redshift 0.4
4212     $ astcosmiccal -z0.4 -LAg
4213
4214     ## Assume Lambda and matter density of 0.7 and 0.3 and print
4215     ## basic cosmological parameters for redshift 2.1:
4216     $ astcosmiccal -l0.7 -m0.3 -z2.1
4217
4218     ## Print wavelength of all pre-defined spectral lines when
4219     ## Lyman-alpha is observed at 4000 Angstroms.
4220     $ astcosmiccal --obsline=lyalpha,4000 --listlinesatz
4221
4222   The input parameters (for example current matter density, etc) can be
4223given as command-line options or in the configuration files, see *note
4224Configuration files::.  For a definition of the different parameters,
4225please see the sections prior to this.  If no redshift is given,
4226CosmicCalculator will just print its input parameters and abort.  For a
4227full list of the input options, please see *note CosmicCalculator input
4228options::.
4229
4230   Without any particular output requested (and only a given redshift),
4231CosmicCalculator will print all basic cosmological calculations (one per
4232line) with some explanations before each.  This can be good when you
4233want a general feeling of the conditions at a specific redshift.
4234Alternatively, if any specific calculation(s) are requested (its
4235possible to call more than one), only the requested value(s) will be
4236calculated and printed with one character space between them.  In this
4237case, no description or units will be printed.  See *note
4238CosmicCalculator basic cosmology calculations:: for the full list of
4239these options along with some explanations how when/how they can be
4240useful.
4241
4242   Another common operation in observational cosmology is dealing with
4243spectral lines at different redshifts.  CosmicCalculator also has
4244features to help in such situations, please see *note CosmicCalculator
4245spectral line calculations::.
4246
4247* Menu:
4248
4249* CosmicCalculator input options::  Options to specify input conditions.
4250* CosmicCalculator basic cosmology calculations::  Like distance modulus, distances and etc.
4251* CosmicCalculator spectral line calculations::  How they get affected by redshift.
4252
4253
4254File: gnuastro.info,  Node: CosmicCalculator input options,  Next: CosmicCalculator basic cosmology calculations,  Prev: Invoking astcosmiccal,  Up: Invoking astcosmiccal
4255
42569.1.3.1 CosmicCalculator input options
4257......................................
4258
4259The inputs to CosmicCalculator can be specified with the following
4260options:
4261
4262‘-z FLT’
4263‘--redshift=FLT’
4264     The redshift of interest.  There are two other ways that you can
4265     specify the target redshift: 1) Spectral lines and their observed
4266     wavelengths, see ‘--obsline’.  2) Velocity, see ‘--velocity’.
4267     Hence this option cannot be called with ‘--obsline’ or
4268     ‘--velocity’.
4269
4270‘-y FLT’
4271‘--velocity=FLT’
4272     Input velocity in km/s.  The given value will be converted to
4273     redshift internally, and used in any subsequent calculation.  This
4274     option is thus an alternative to ‘--redshift’ or ‘--obsline’, it
4275     cannot be used with them.  The conversion will be done with the
4276     more general and accurate relativistic equation of
4277     $1+z=\sqrt{(c+v)/(c-v)}$, not the simplified $z\approx v/c$.
4278
4279‘-H FLT’
4280‘--H0=FLT’
4281     Current expansion rate (in km sec$^{-1}$ Mpc$^{-1}$).
4282
4283‘-l FLT’
4284‘--olambda=FLT’
4285     Cosmological constant density divided by the critical density in
4286     the current Universe ($\Omega_{\Lambda,0}$).
4287
4288‘-m FLT’
4289‘--omatter=FLT’
4290     Matter (including massive neutrinos) density divided by the
4291     critical density in the current Universe ($\Omega_{m,0}$).
4292
4293‘-r FLT’
4294‘--oradiation=FLT’
4295     Radiation density divided by the critical density in the current
4296     Universe ($\Omega_{r,0}$).
4297
4298‘-O STR/FLT,FLT’
4299‘--obsline=STR/FLT,FLT’
4300     Find the redshift to use in next steps based on the rest-frame and
4301     observed wavelengths of a line.  This option is thus an alternative
4302     to ‘--redshift’ or ‘--velocity’, it cannot be used with them.
4303     Wavelengths are assumed to be in Angstroms.  The first argument
4304     identifies the line.  It can be one of the standard names below, or
4305     any rest-frame wavelength in Angstroms.  The second argument is the
4306     observed wavelength of that line.  For example
4307     ‘--obsline=lyalpha,6000’ is the same as ‘--obsline=1215.64,6000’.
4308
4309     The pre-defined names are listed below, sorted from red (longer
4310     wavelength) to blue (shorter wavelength).  You can get this list on
4311     the command-line with the ‘--listlines’.
4312
4313     ‘siired’
4314          [6731Å] SII doublet’s redder line.
4315
4316     ‘sii’
4317          [6724Å] SII doublet’s mean center at .
4318
4319     ‘siiblue’
4320          [6717Å] SII doublet’s bluer line.
4321
4322     ‘niired’
4323          [6584Å] NII doublet’s redder line.
4324
4325     ‘nii’
4326          [6566Å] NII doublet’s mean center.
4327
4328     ‘halpha’
4329          [6562.8Å] H-$\alpha$ line.
4330
4331     ‘niiblue’
4332          [6548Å] NII doublet’s bluer line.
4333
4334     ‘oiiired-vis’
4335          [5007Å] OIII doublet’s redder line in the visible.
4336
4337     ‘oiii-vis’
4338          [4983Å] OIII doublet’s mean center in the visible.
4339
4340     ‘oiiiblue-vis’
4341          [4959Å] OIII doublet’s bluer line in the visible.
4342
4343     ‘hbeta’
4344          [4861.36Å] H-$\beta$ line.
4345
4346     ‘heii-vis’
4347          [4686Å] HeII doublet’s redder line in the visible.
4348
4349     ‘hgamma’
4350          [4340.46Å] H-$\gamma$ line.
4351
4352     ‘hdelta’
4353          [4101.74Å] H-$\delta$ line.
4354
4355     ‘hepsilon’
4356          [3970.07Å] H-$\epsilon$ line.
4357
4358     ‘neiii’
4359          [3869Å] NEIII line.
4360
4361     ‘oiired’
4362          [3729Å] OII doublet’s redder line.
4363
4364     ‘oii’
4365          [3727.5Å] OII doublet’s mean center.
4366
4367     ‘oiiblue’
4368          [3726Å] OII doublet’s bluer line.
4369
4370     ‘blimit’
4371          [3646Å] Balmer limit.
4372
4373     ‘mgiired’
4374          [2803Å] MgII doublet’s redder line.
4375
4376     ‘mgii’
4377          [2799.5Å] MgII doublet’s mean center.
4378
4379     ‘mgiiblue’
4380          [2796Å] MgII doublet’s bluer line.
4381
4382     ‘ciiired’
4383          [1909Å] CIII doublet’s redder line.
4384
4385     ‘ciii’
4386          [1908Å] CIII doublet’s mean center.
4387
4388     ‘ciiiblue’
4389          [1907Å] CIII doublet’s bluer line.
4390
4391     ‘si_iiired’
4392          [1892Å] SiIII doublet’s redder line.
4393
4394     ‘si_iii’
4395          [1887.5Å] SiIII doublet’s mean center.
4396
4397     ‘si_iiiblue’
4398          [1883Å] SiIII doublet’s bluer line.
4399
4400     ‘oiiired-uv’
4401          [1666Å] OIII doublet’s redder line in the ultra-violet.
4402
4403     ‘oiii-uv’
4404          [1663.5Å] OIII doublet’s mean center in the ultra-violet.
4405
4406     ‘oiiiblue-uv’
4407          [1661Å] OIII doublet’s bluer line in the ultra-violet.
4408
4409     ‘heii-uv’
4410          [1640Å] HeII doublet’s bluer line in the ultra-violet.
4411
4412     ‘civred’
4413          [1551Å] CIV doublet’s redder line.
4414
4415     ‘civ’
4416          [1549Å] CIV doublet’s mean center.
4417
4418     ‘civblue’
4419          [1548Å] CIV doublet’s bluer line.
4420
4421     ‘nv’
4422          [1240Å] NV (four times ionized Sodium).
4423
4424     ‘lyalpha’
4425          [1215.67Å] Lyman-$\alpha$ line.
4426
4427     ‘lybeta’
4428          [1025.7Å] Lyman-$\beta$ line.
4429
4430     ‘lygamma’
4431          [972.54Å] Lyman-$\gamma$ line.
4432
4433     ‘lydelta’
4434          [949.74Å] Lyman-$\delta$ line.
4435
4436     ‘lyepsilon’
4437          [937.80Å] Lyman-$\epsilon$ line.
4438
4439     ‘lylimit’
4440          [912Å] Lyman limit.
4441
4442
4443File: gnuastro.info,  Node: CosmicCalculator basic cosmology calculations,  Next: CosmicCalculator spectral line calculations,  Prev: CosmicCalculator input options,  Up: Invoking astcosmiccal
4444
44459.1.3.2 CosmicCalculator basic cosmology calculations
4446.....................................................
4447
4448By default, when no specific calculations are requested,
4449CosmicCalculator will print a complete set of all its calculators (one
4450line for each calculation, see *note Invoking astcosmiccal::).  The full
4451list of calculations can be useful when you don’t want any specific
4452value, but just a general view.  In other contexts (for example in a
4453batch script or during a discussion), you know exactly what you want and
4454don’t want to be distracted by all the extra information.
4455
4456   You can use any number of the options described below in any order.
4457When any of these options are requested, CosmicCalculator’s output will
4458just be a single line with a single space between the (possibly)
4459multiple values.  In the example below, only the tangential distance
4460along one arc-second (in kpc), absolute magnitude conversion, and age of
4461the universe at redshift 2 are printed (recall that you can merge short
4462options together, see *note Options::).
4463
4464     $ astcosmiccal -z2 -sag
4465     8.585046 44.819248 3.289979
4466
4467   Here is one example of using this feature in scripts: by adding the
4468following two lines in a script to keep/use the comoving volume with
4469varying redshifts:
4470
4471     z=3.12
4472     vol=$(astcosmiccal --redshift=$z --volume)
4473
4474In a script, this operation might be necessary for a large number of
4475objects (several of galaxies in a catalog for example).  So the fact
4476that all the other default calculations are ignored will also help you
4477get to your result faster.
4478
4479   If you are indeed dealing with many (for example thousands) of
4480redshifts, using CosmicCalculator is not the best/fastest solution.
4481Because it has to go through all the configuration files and
4482preparations for each invocation.  To get the best efficiency (least
4483overhead), we recommend using Gnuastro’s cosmology library (see *note
4484Cosmology library::).  CosmicCalculator also calls the library functions
4485defined there for its calculations, so you get the same result with no
4486overhead.  Gnuastro also has libraries for easily reading tables into a
4487C program, see *note Table input output::.  Afterwards, you can easily
4488build and run your C program for the particular processing with *note
4489BuildProgram::.
4490
4491   If you just want to inspect the value of a variable visually, the
4492description (which comes with units) might be more useful.  In such
4493cases, the following command might be better.  The other calculations
4494will also be done, but they are so fast that you will not notice on
4495modern computers (the time it takes your eye to focus on the result is
4496usually longer than the processing: a fraction of a second).
4497
4498     $ astcosmiccal --redshift=0.832 | grep volume
4499
4500   The full list of CosmicCalculator’s specific calculations is present
4501below in two groups: basic cosmology calculations and those related to
4502spectral lines.  In case you have forgot the units, you can use the
4503‘--help’ option which has the units along with a short description.
4504
4505‘-e’
4506‘--usedredshift’
4507     The redshift that was used in this run.  In many cases this is the
4508     main input parameter to CosmicCalculator, but it is useful in
4509     others.  For example in combination with ‘--obsline’ (where you
4510     give an observed and rest-frame wavelength and would like to know
4511     the redshift) or with ‘--velocity’ (where you specify the velocity
4512     instead of redshift).  Another example is when you run
4513     CosmicCalculator in a loop, while changing the redshift and you
4514     want to keep the redshift value with the resulting calculation.
4515
4516‘-Y’
4517‘--usedvelocity’
4518     The velocity (in km/s) that was used in this run.  The conversion
4519     from redshift will be done with the more general and accurate
4520     relativistic equation of $1+z=\sqrt{(c+v)/(c-v)}$, not the
4521     simplified $z\approx v/c$.
4522
4523‘-G’
4524‘--agenow’
4525     The current age of the universe (given the input parameters) in Ga
4526     (Giga annum, or billion years).
4527
4528‘-C’
4529‘--criticaldensitynow’
4530     The current critical density (given the input parameters) in grams
4531     per centimeter-cube ($g/cm^3$).
4532
4533‘-d’
4534‘--properdistance’
4535     The proper distance (at current time) to object at the given
4536     redshift in Megaparsecs (Mpc).  See *note Distance on a 2D curved
4537     space:: for a description of the proper distance.
4538
4539‘-A’
4540‘--angulardimdist’
4541     The angular diameter distance to object at given redshift in
4542     Megaparsecs (Mpc).
4543
4544‘-s’
4545‘--arcsectandist’
4546     The tangential distance covered by 1 arc-seconds at the given
4547     redshift in kiloparsecs (Kpc).  This can be useful when trying to
4548     estimate the resolution or pixel scale of an instrument (usually in
4549     units of arc-seconds) at a given redshift.
4550
4551‘-L’
4552‘--luminositydist’
4553     The luminosity distance to object at given redshift in Megaparsecs
4554     (Mpc).
4555
4556‘-u’
4557‘--distancemodulus’
4558     The distance modulus at given redshift.
4559
4560‘-a’
4561‘--absmagconv’
4562     The conversion factor (addition) to absolute magnitude.  Note that
4563     this is practically the distance modulus added with
4564     $-2.5\log{(1+z)}$ for the desired redshift based on the input
4565     parameters.  Once the apparent magnitude and redshift of an object
4566     is known, this value may be added with the apparent magnitude to
4567     give the object’s absolute magnitude.
4568
4569‘-g’
4570‘--age’
4571     Age of the universe at given redshift in Ga (Giga annum, or billion
4572     years).
4573
4574‘-b’
4575‘--lookbacktime’
4576     The look-back time to given redshift in Ga (Giga annum, or billion
4577     years).  The look-back time at a given redshift is defined as the
4578     current age of the universe (‘--agenow’) subtracted by the age of
4579     the universe at the given redshift.
4580
4581‘-c’
4582‘--criticaldensity’
4583     The critical density at given redshift in grams per centimeter-cube
4584     ($g/cm^3$).
4585
4586‘-v’
4587‘--onlyvolume’
4588     The comoving volume in Megaparsecs cube (Mpc$^3$) until the desired
4589     redshift based on the input parameters.
4590
4591
4592File: gnuastro.info,  Node: CosmicCalculator spectral line calculations,  Prev: CosmicCalculator basic cosmology calculations,  Up: Invoking astcosmiccal
4593
45949.1.3.3 CosmicCalculator spectral line calculations
4595...................................................
4596
4597At different redshifts, observed spectral lines are shifted compared to
4598their rest frame wavelengths with this simple relation:
4599$\lambda_{obs}=\lambda_{rest}(1+z)$.  Although this relation is very
4600simple and can be done for one line in the head (or a simple
4601calculator!), it slowly becomes tiring when dealing with a lot of lines
4602or redshifts, or some precision is necessary.  The options in this
4603section are thus provided to greatly simplify usage of this simple
4604equation, and also helping by storing a list of pre-defined spectral
4605line wavelengths.
4606
4607   For example if you want to know the wavelength of the $H\alpha$ line
4608(at 6562.8 Angstroms in rest frame), when $Ly\alpha$ is at 8000
4609Angstroms, you can call CosmicCalculator like the first example below.
4610And if you want the wavelength of all pre-defined spectral lines at this
4611redshift, you can use the second command.
4612
4613     $ astcosmiccal --obsline=lyalpha,8000 --lineatz=halpha
4614     $ astcosmiccal --obsline=lyalpha,8000 --listlinesatz
4615
4616   Bellow you can see the printed/output calculations of
4617CosmicCalculator that are related to spectral lines.  Note that
4618‘--obsline’ is an input parameter, so its discussed (with the full list
4619of known lines) in *note CosmicCalculator input options::.
4620
4621‘--listlines’
4622     List the pre-defined rest frame spectral line wavelengths and their
4623     names on standard output, then abort CosmicCalculator.  When this
4624     option is given, other operations on the command-line will be
4625     ignored.  This is convenient when you forget the specific name of
4626     the spectral line used within Gnuastro, or when you forget the
4627     exact wavelength of a certain line.
4628
4629     These names can be used with the options that deal with spectral
4630     lines, for example ‘--obsline’ and ‘--lineatz’ (*note
4631     CosmicCalculator basic cosmology calculations::).
4632
4633     The format of the output list is a two-column table, with
4634     Gnuastro’s text table format (see *note Gnuastro text table
4635     format::).  Therefore, if you are only looking for lines in a
4636     specific range, you can pipe the output into Gnuastro’s table
4637     program and use its ‘--range’ option on the ‘wavelength’ (first)
4638     column.  For example, if you only want to see the lines between
4639     4000 and 6000 Angstroms, you can run this command:
4640
4641          $ astcosmiccal --listlines \
4642                         | asttable --range=wavelength,4000,6000
4643
4644     And if you want to use the list later and have it as a table in a
4645     file, you can easily add the ‘--output’ (or ‘-o’) option to the
4646     ‘asttable’ command, and specify the filename, for example
4647     ‘--output=lines.fits’ or ‘--output=lines.txt’.
4648
4649‘--listlinesatz’
4650     Similar to ‘--listlines’ (above), but the printed wavelength is not
4651     in the rest frame, but redshifted to the given redshift.  Recall
4652     that the redshift can be specified by ‘--redshift’ directly or by
4653     ‘--obsline’, see *note CosmicCalculator input options::.
4654
4655‘-i STR/FLT4656‘--lineatz=STR/FLT4657     The wavelength of the specified line at the redshift given to
4658     CosmicCalculator.  The line can be specified either by its name or
4659     directly as a number (its wavelength).  To get the list of
4660     pre-defined names for the lines and their wavelength, you can use
4661     the ‘--listlines’ option, see *note CosmicCalculator input
4662     options::.  In the former case (when a name is given), the returned
4663     number is in units of Angstroms.  In the latter (when a number is
4664     given), the returned value is the same units of the input number
4665     (assuming its a wavelength).
4666
4667
4668File: gnuastro.info,  Node: Installed scripts,  Next: Library,  Prev: High-level calculations,  Up: Top
4669
467010 Installed scripts
4671********************
4672
4673Gnuastro’s programs (introduced in previous chapters) are designed to be
4674highly modular and thus contain lower-level operations on the data.
4675However, in many contexts, certain higher-level are also shared between
4676many contexts.  For example a sequence of calls to multiple Gnuastro
4677programs, or a special way of running a program and treating the output.
4678To facilitate such higher-level data analysis, Gnuastro also installs
4679some scripts on your system with the (‘astscript-’) prefix (in contrast
4680to the other programs that only have the ‘ast’ prefix).
4681
4682   Like all of Gnuastro’s source code, these scripts are also heavily
4683commented.  They are written in portable shell scripts (command-line
4684environments), which doesn’t need compilation.  Therefore, if you open
4685the installed scripts in a text editor, you can actually read them(1).
4686For example with this command (just replace ‘nano’ with your favorite
4687text editor, like ‘emacs’ or ‘vim’):
4688
4689     $ nano $(which astscript-NAME)
4690
4691   Shell scripting is the same language that you use when typing on the
4692command-line.  Therefore shell scripting is much more widely known and
4693used compared to C (the language of other Gnuastro programs).  Because
4694Gnuastro’s installed scripts do higher-level operations, customizing
4695these scripts for a special project will be more common than the
4696programs.
4697
4698   These scripts also accept options and are in many ways similar to the
4699programs (see *note Common options::) with some minor differences:
4700
4701   • Currently they don’t accept configuration files themselves.
4702     However, the configuration files of the Gnuastro programs they call
4703     are indeed parsed and used by those programs.
4704
4705     As a result, they don’t have the following options:
4706     ‘--checkconfig’, ‘--config’, ‘--lastconfig’, ‘--onlyversion’,
4707     ‘--printparams’, ‘--setdirconf’ and ‘--setusrconf’.
4708
4709   • They don’t directly allocate any memory, so there is no
4710     ‘--minmapsize’.
4711
4712   • They don’t have an independent ‘--usage’ option: when called with
4713     ‘--usage’, they just recommend running ‘--help’.
4714
4715   • The output of ‘--help’ is not configurable like the programs (see
4716     *note --help::).
4717
4718   • The scripts will commonly use your installed shell and other basic
4719     command-line tools (for example AWK or SED). Different systems have
4720     different versions and implementations of these basic tools (for
4721     example GNU/Linux systems use GNU Bash, GNU AWK and GNU SED which
4722     are far more advanced and up to date then the minimalist AWK and
4723     SED of most other systems).  Therefore, unexpected errors in these
4724     tools might come up when you run these scripts on non-GNU/Linux
4725     operating systems.  If you do confront such strange errors, please
4726     submit a bug report so we fix it as soon as possible (see *note
4727     Report a bug::).
4728
4729* Menu:
4730
4731* Sort FITS files by night::    Sort many files by date.
4732* Generate radial profile::     Radial profile of an object in an image.
4733* SAO DS9 region files from table::  Create ds9 region file from a table.
4734
4735   ---------- Footnotes ----------
4736
4737   (1) Gnuastro’s installed programs (those only starting with ‘ast’)
4738aren’t human-readable.  They are written in C and need to be compiled
4739before execution.  Compilation optimizes the steps into the low-level
4740hardware CPU instructions/language to improve efficiency.  Because
4741compiled programs don’t need an interpreter like Bash on every run, they
4742are much faster and more independent than scripts.  To read the source
4743code of the programs, look into the ‘bin/progname’ directory of
4744Gnuastro’s source (*note Downloading the source::).  If you would like
4745to read more about why C was chosen for the programs, please see *note
4746Why C::.
4747
4748
4749File: gnuastro.info,  Node: Sort FITS files by night,  Next: Generate radial profile,  Prev: Installed scripts,  Up: Installed scripts
4750
475110.1 Sort FITS files by night
4752=============================
4753
4754FITS images usually contain (several) keywords for preserving important
4755dates.  In particular, for lower-level data, this is usually the
4756observation date and time (for example, stored in the ‘DATE-OBS’ keyword
4757value).  When analyzing observed datasets, many calibration steps (like
4758the dark, bias or flat-field), are commonly calculated on a
4759per-observing-night basis.
4760
4761   However, the FITS standard’s date format (‘YYYY-MM-DDThh:mm:ss.ddd’)
4762is based on the western (Gregorian) calendar.  Dates that are stored in
4763this format are complicated for automatic processing: a night starts in
4764the final hours of one calendar day, and extends to the early hours of
4765the next calendar day.  As a result, to identify datasets from one
4766night, we commonly need to search for two dates.  However calendar
4767peculiarities can make this identification very difficult.  For example
4768when an observation is done on the night separating two months (like the
4769night starting on March 31st and going into April 1st), or two years
4770(like the night starting on December 31st 2018 and going into January
47711st, 2019).  To account for such situations, it is necessary to keep
4772track of how many days are in a month, and leap years, etc.
4773
4774   Gnuastro’s ‘astscript-sort-by-night’ script is created to help in
4775such important scenarios.  It uses *note Fits:: to convert the FITS date
4776format into the Unix epoch time (number of seconds since 00:00:00 of
4777January 1st, 1970), using the ‘--datetosec’ option.  The Unix epoch time
4778is a single number (integer, if not given in sub-second precision),
4779enabling easy comparison and sorting of dates after January 1st, 1970.
4780
4781   You can use this script as a basis for making a much more highly
4782customized sorting script.  Here are some examples
4783
4784   • If you need to copy the files, but only need a single extension
4785     (not the whole file), you can add a step just before the making of
4786     the symbolic links, or copies, and change it to only copy a certain
4787     extension of the FITS file using the Fits program’s ‘--copy’
4788     option, see *note HDU information and manipulation::.
4789
4790   • If you need to classify the files with finer detail (for example
4791     the purpose of the dataset), you can add a step just before the
4792     making of the symbolic links, or copies, to specify a file-name
4793     prefix based on other certain keyword values in the files.  For
4794     example when the FITS files have a keyword to specify if the
4795     dataset is a science, bias, or flat-field image.  You can read it
4796     and to add a ‘sci-’, ‘bias-’, or ‘flat-’ to the created file (after
4797     the ‘--prefix’) automatically.
4798
4799     For example, let’s assume the observing mode is stored in the
4800     hypothetical ‘MODE’ keyword, which can have three values of
4801     ‘BIAS-IMAGE’, ‘SCIENCE-IMAGE’ and ‘FLAT-EXP’.  With the step below,
4802     you can generate a mode-prefix, and add it to the generated
4803     link/copy names (just correct the filename and extension of the
4804     first line to the script’s variables):
4805
4806          modepref=$(astfits infile.fits -h1 \
4807                             | sed -e"s/'/ /g" \
4808                             | awk '$1=="MODE"{ \
4809                                 if($3=="BIAS-IMAGE") print "bias-"; \
4810                                 else if($3=="SCIENCE-IMAGE") print "sci-"; \
4811                                 else if($3==FLAT-EXP) print "flat-"; \
4812                                 else print $3, "NOT recognized"; exit 1}')
4813
4814     Here is a description of it.  We first use ‘astfits’ to print all
4815     the keywords in extension ‘1’ of ‘infile.fits’.  In the FITS
4816     standard, string values (that we are assuming here) are placed in
4817     single quotes (<'>) which are annoying in this context/use-case.
4818     Therefore, we pipe the output of ‘astfits’ into ‘sed’ to remove all
4819     such quotes (substituting them with a blank space).  The result is
4820     then piped to AWK for giving us the final mode-prefix: with
4821     ‘$1=="MODE"’, we ask AWK to only consider the line where the first
4822     column is ‘MODE’.  There is an equal sign between the key name and
4823     value, so the value is the third column (‘$3’ in AWK). We thus use
4824     a simple ‘if-else’ structure to look into this value and print our
4825     custom prefix based on it.  The output of AWK is then stored in the
4826     ‘modepref’ shell variable which you can add to the link/copy name.
4827
4828     With the solution above, the increment of the file counter for each
4829     night will be independent of the mode.  If you want the counter to
4830     be mode-dependent, you can add a different counter for each mode
4831     and use that counter instead of the generic counter for each night
4832     (based on the value of ‘modepref’).  But we’ll leave the
4833     implementation of this step to you as an exercise.
4834
4835* Menu:
4836
4837* Invoking astscript-sort-by-night::  Inputs and outputs to this script.
4838
4839
4840File: gnuastro.info,  Node: Invoking astscript-sort-by-night,  Prev: Sort FITS files by night,  Up: Sort FITS files by night
4841
484210.1.1 Invoking astscript-sort-by-night
4843---------------------------------------
4844
4845This installed script will read a FITS date formatted value from the
4846given keyword, and classify the input FITS files into individual nights.
4847For more on installed scripts please see (see *note Installed
4848scripts::).  This script can be used with the following general
4849template:
4850
4851     $ astscript-sort-by-night [OPTION...] FITS-files
4852
4853One line examples:
4854
4855     ## Use the DATE-OBS keyword
4856     $ astscript-sort-by-night --key=DATE-OBS /path/to/data/*.fits
4857
4858     ## Make links to the input files with the `img-' prefix
4859     $ astscript-sort-by-night --link --prefix=img- /path/to/data/*.fits
4860
4861   This script will look into a HDU/extension (‘--hdu’) for a keyword
4862(‘--key’) in the given FITS files and interpret the value as a date.
4863The inputs will be separated by "night"s (11:00a.m to next day’s
486410:59:59a.m, spanning two calendar days, exact hour can be set with
4865‘--hour’).
4866
4867   The default output is a list of all the input files along with the
4868following two columns: night number and file number in that night
4869(sorted by time).  With ‘--link’ a symbolic link will be made (one for
4870each input) that contains the night number, and number of file in that
4871night (sorted by time), see the description of ‘--link’ for more.  When
4872‘--copy’ is used instead of a link, a copy of the inputs will be made
4873instead of symbolic link.
4874
4875   Below you can see one example where all the ‘target-*.fits’ files in
4876the ‘data’ directory should be separated by observing night according to
4877the ‘DATE-OBS’ keyword value in their second extension (number ‘1’,
4878recall that HDU counting starts from 0).  You can see the output after
4879the ‘ls’ command.
4880
4881     $ astscript-sort-by-night -pimg- -h1 -kDATE-OBS data/target-*.fits
4882     $ ls
4883     img-n1-1.fits img-n1-2.fits img-n2-1.fits ...
4884
4885   The outputs can be placed in a different (already existing) directory
4886by including that directory’s name in the ‘--prefix’ value, for example
4887‘--prefix=sorted/img-’ will put them all under the ‘sorted’ directory.
4888
4889   This script can be configured like all Gnuastro’s programs (through
4890command-line options, see *note Common options::), with some minor
4891differences that are described in *note Installed scripts::.  The
4892particular options to this script are listed below:
4893
4894‘-h STR’
4895‘--hdu=STR’
4896     The HDU/extension to use in all the given FITS files.  All of the
4897     given FITS files must have this extension.
4898
4899‘-k STR’
4900‘--key=STR’
4901     The keyword name that contains the FITS date format to
4902     classify/sort by.
4903
4904‘-H FLT’
4905‘--hour=FLT’
4906     The hour that defines the next “night”.  By default, all times
4907     before 11:00a.m are considered to belong to the previous calendar
4908     night.  If a sub-hour value is necessary, it should be given in
4909     units of hours, for example ‘--hour=9.5’ corresponds to 9:30a.m.
4910
4911     *Dealing with time zones:* The time that is recorded in ‘--key’ may
4912     be in UTC (Universal Time Coordinate).  However, the organization
4913     of the images taken during the night depends on the local time.  It
4914     is possible to take this into account by setting the ‘--hour’
4915     option to the local time in UTC.
4916
4917     For example, consider a set of images taken in Auckland (New
4918     Zealand, UTC+12) during different nights.  If you want to classify
4919     these images by night, you have to know at which time (in UTC time)
4920     the Sun rises (or any other separator/definition of a different
4921     night).  For example if your observing night finishes before
4922     9:00a.m in Auckland, you can use ‘--hour=21’.  Because in Auckland
4923     the local time of 9:00 corresponds to 21:00 UTC.
4924
4925‘-l’
4926‘--link’
4927     Create a symbolic link for each input FITS file.  This option
4928     cannot be used with ‘--copy’.  The link will have a standard name
4929     in the following format (variable parts are written in ‘CAPITAL’
4930     letters and described after it):
4931
4932          PnN-I.fits
4933
4934     ‘P’
4935          This is the value given to ‘--prefix’.  By default, its value
4936          is ‘./’ (to store the links in the directory this script was
4937          run in).  See the description of ‘--prefix’ for more.
4938     ‘N’
4939          This is the night-counter: starting from 1.  ‘N’ is just
4940          incremented by 1 for the next night, no matter how many nights
4941          (without any dataset) there are between two subsequent
4942          observing nights (its just an identifier for each night which
4943          you can easily map to different calendar nights).
4944     ‘I’
4945          File counter in that night, sorted by time.
4946
4947‘-c’
4948‘--copy’
4949     Make a copy of each input FITS file with the standard naming
4950     convention described in ‘--link’.  With this option, instead of
4951     making a link, a copy is made.  This option cannot be used with
4952     ‘--link’.
4953
4954‘-p STR’
4955‘--prefix=STR’
4956     Prefix to append before the night-identifier of each newly created
4957     link or copy.  This option is thus only relevant with the ‘--copy’
4958     or ‘--link’ options.  See the description of ‘--link’ for how its
4959     used.  For example, with ‘--prefix=img-’, all the created file
4960     names in the current directory will start with ‘img-’, making
4961     outputs like ‘img-n1-1.fits’ or ‘img-n3-42.fits’.
4962
4963     ‘--prefix’ can also be used to store the links/copies in another
4964     directory relative to the directory this script is being run (it
4965     must already exist).  For example
4966     ‘--prefix=/path/to/processing/img-’ will put all the links/copies
4967     in the ‘/path/to/processing’ directory, and the files (in that
4968     directory) will all start with ‘img-’.
4969
4970‘--stdintimeout=INT’
4971     Number of micro-seconds to wait for standard input within this
4972     script.  This does not correspond to general inputs into the
4973     script, inputs to the script should always be given as a file.
4974     However, within the script, pipes are often used to pass the output
4975     of one program to another.  The value given to this option will be
4976     passed to those internal pipes.  When running this script, if you
4977     confront an error, saying “No input!”, you should be able to fix it
4978     by giving a larger number to this option (the default value is
4979     10000000 micro-seconds or 10 seconds).
4980
4981
4982File: gnuastro.info,  Node: Generate radial profile,  Next: SAO DS9 region files from table,  Prev: Sort FITS files by night,  Up: Installed scripts
4983
498410.2 Generate radial profile
4985============================
4986
4987The 1 dimensional radial profile of an object is an important parameter
4988in many aspects of astronomical image processing.  For example, you want
4989to study how the light of a galaxy is distributed as a function of the
4990radial distance from the center.  In other cases, the radial profile of
4991a star can show the PSF (see *note PSF::).  Gnuastro’s
4992‘astscript-radial-profile’ script is created to obtain such radial
4993profiles for one object within an image.  This script uses *note
4994MakeProfiles:: to generate elliptical apertures with the values equal to
4995the distance from the center of the object and *note MakeCatalog:: for
4996measuring the values over the apertures.
4997
4998* Menu:
4999
5000* Invoking astscript-radial-profile::  How to call astscript-radial-profile
5001
5002
5003File: gnuastro.info,  Node: Invoking astscript-radial-profile,  Prev: Generate radial profile,  Up: Generate radial profile
5004
500510.2.1 Invoking astscript-radial-profile
5006----------------------------------------
5007
5008This installed script will measure the radial profile of an object
5009within an image.  For more on installed scripts please see (see *note
5010Installed scripts::).  This script can be used with the following
5011general template:
5012
5013     $ astscript-radial-profile [OPTION...] FITS-file
5014
5015Examples:
5016
5017     ## Generate the radial profile with default options (assuming the
5018     ## object is in the center of the image, and using the mean).
5019     $ astscript-radial-profile image.fits
5020
5021     ## Generate the radial profile centered at x=44 and y=37 (in pixels),
5022     ## up to  a radial distance of 19 pixels, use the mean value.
5023     $ astscript-radial-profile image.fits --center=44,37 --rmax=19
5024
5025     ## Generate the radial profile centered at x=44 and y=37 (in pixels),
5026     ## up to a radial distance of 100 pixels, compute sigma clipped
5027     ## mean and standard deviation (sigclip-mean and sigclip-std) using
5028     ## 5 sigma and 0.1 tolerance (default is 3 sigma and 0.2 tolerance).
5029     $ astscript-radial-profile image.fits --center=44,37 --rmax=100 \
5030                                --sigmaclip=5,0.1 \
5031                                --measure=sigclip-mean,sigclip-std
5032
5033     ## Generate the radial profile centered at RA=20.53751695,
5034     ## DEC=0.9454292263, up to a radial distance of 88 pixels,
5035     ## axis ratio equal to 0.32, and position angle of 148 deg.
5036     ## Name the output table as `radial-profile.fits'
5037     $ astscript-radial-profile image.fits --mode=wcs \
5038                                --center=20.53751695,0.9454292263 \
5039                                --rmax=88 --axisratio=0.32 \
5040                                --positionangle=148 -oradial-profile.fits
5041
5042     ## Generate the radial profile centered at RA=40.062675270971,
5043     ## DEC=-8.1511992735126, up to a radial distance of 20 pixels,
5044     ## and calculate the SNR using the INPUT-NO-SKY and SKY-STD
5045     ## extensions of the NoiseChisel output file.
5046     $ astscript-radial-profile image_detected.fits -hINPUT-NO-SKY \
5047                                --mode=wcs --measure=sn \
5048                                --center=40.062675270971,-8.1511992735126 \
5049                                --rmax=20 --stdhdu=SKY_STD
5050
5051     ## Generate the radial profile centered at RA=40.062675270971,
5052     ## DEC=-8.1511992735126, up to a radial distance of 20 pixels,
5053     ## and compute the SNR with a fixed value for std, std=10.
5054     $ astscript-radial-profile image.fits -h1 --mode=wcs --rmax=20 \
5055                                --center=40.062675270971,-8.1511992735126 \
5056                                --measure=sn --instd=10
5057
5058     ## Generate the radial profile centered at X=1201, Y=1201 pixels, up
5059     ## to a radial distance of 20 pixels and compute the median and the
5060     ## SNR using the first extension of sky-std.fits as the dataset for std
5061     ## values.
5062     $ astscript-radial-profile image.fits -h1 --mode=img --rmax=20 \
5063                                --center=1201,1201 --measure=median,sn \
5064                                --instd=sky-std.fits
5065
5066   This installed script will read a FITS image and will use it as the
5067basis for constructing the radial profile.  The output radial profile is
5068a table (FITS or plain-text) containing the radial distance from the
5069center in the first row and the specified measurements in the other
5070columns (mean, median, sigclip-mean, sigclip-median, etc.).
5071
5072   To measure the radial profile, this script needs to generate
5073temporary files.  All these temporary files will be created within the
5074directory given to the ‘--tmpdir’ option.  When ‘--tmpdir’ is not
5075called, a temporary directory (with a name based on the inputs) will be
5076created in the running directory.  If the directory doesn’t exist at
5077run-time, this script will create it.  After the output is created, this
5078script will delete the directory by default, unless you call the
5079‘--keeptmp’ option.
5080
5081   With the default options, the script will generate a circular radial
5082profile using the mean value and centered at the center of the image.
5083In order to have more flexibility, several options are available to
5084configure for the desired radial profile.  In this sense, you can change
5085the center position, the maximum radius, the axis ratio and the position
5086angle (elliptical apertures are considered), the operator for obtaining
5087the profiles, and others (described below).
5088
5089*Debug your profile:* to debug your results, especially close to the
5090center of your object, you can see the radial distance associated to
5091every pixel in your input.  To do this, use ‘--keeptmp’ to keep the
5092temporary files, and compare ‘crop.fits’ (crop of your input image
5093centered on your desired coordinate) with ‘apertures.fits’ (radial
5094distance of each pixel).
5095
5096*Finding properties of your elliptical target: * you want to measure the
5097radial profile of a galaxy, but don’t know its exact location, position
5098angle or axis ratio.  To obtain these values, you can use *note
5099NoiseChisel:: to detect signal in the image, feed it to *note Segment::
5100to do basic segmentation, then use *note MakeCatalog:: to measure the
5101center (‘--x’ and ‘--y’ in MakeCatalog), axis ratio (‘--axisratio’) and
5102position angle (‘--positionangle’).
5103
5104*Masking other sources:* The image of an astronomical object will
5105usually have many other sources with your main target.  A crude solution
5106is to use sigma-clipped measurements for the profile.  However,
5107sigma-clipped measurements can easily be biased when the number of
5108sources at each radial distance increases at larger distances.
5109Therefore a robust solution is to mask all other detections within the
5110image.  You can use *note NoiseChisel:: and *note Segment:: to detect
5111and segment the sources, then set all pixels that don’t belong to your
5112target to blank using *note Arithmetic:: (in particular, its ‘where’
5113operator).
5114
5115‘-h STR’
5116‘--hdu=STR’
5117     The HDU/extension of the input image to use.
5118
5119‘-o STR’
5120‘--output=STR’
5121     Filename of measured radial profile.  It can be either a FITS
5122     table, or plain-text table (determined from your given file name
5123     suffix).
5124
5125‘-c FLT[,FLT[,...]]’
5126‘--center=FLT[,FLT[,...]]’
5127     The central position of the radial profile.  This option is used
5128     for placing the center of the profiles.  This parameter is used in
5129     *note Crop:: to center an crop the region.  The positions along
5130     each dimension must be separated by a comma (<,>) and fractions are
5131     also acceptable.  The number of values given to this option must be
5132     the same as the dimensions of the input dataset.  The units of the
5133     coordinates are read based on the value to the ‘--mode’ option, see
5134     below.
5135
5136‘-O STR’
5137‘--mode=STR’
5138     Interpret the center position of the object (values given to
5139     ‘--center’) in image or WCS coordinates.  This option thus accepts
5140     only two values: ‘img’ or ‘wcs’.  By default, it is ‘--mode=img’.
5141
5142‘-R FLT’
5143‘--rmax=FLT’
5144     Maximum radius for the radial profile (in pixels).  By default, the
5145     radial profile will be computed up to a radial distance equal to
5146     the maximum radius that fits into the image (assuming circular
5147     shape).
5148
5149‘-Q FLT’
5150‘--axisratio=FLT’
5151     The axis ratio of the apertures (minor axis divided by the major
5152     axis in a 2D ellipse).  By default (when this option isn’t given),
5153     the radial profile will be circular (axis ratio of 1).  This
5154     parameter is used as the option ‘--qcol’ in the generation of the
5155     apertures with ‘astmkprof’.
5156
5157‘-p FLT’
5158‘--positionangle=FLT’
5159     The position angle (in degrees) of the profiles relative to the
5160     first FITS axis (horizontal when viewed in SAO DS9).  By default,
5161     it is ‘--pangle=0’, which means that the semi-major axis of the
5162     profiles will be parallel to the first FITS axis.
5163
5164‘-m STR’
5165‘--measure=STR’
5166     The operator for measuring the values over each radial distance.
5167     The values given to this option will be directly passed to *note
5168     MakeCatalog::.  As a consequence, all MakeCatalog measurements like
5169     the magnitude, magnitude error, median, mean, signal-to-noise ratio
5170     (S/N), std, surface brightness, sigclip-mean, sigclip-number, etc.
5171     can be used here.  For a full list of MakeCatalog’s measurements,
5172     please run ‘astmkcatalog --help’.  Multiple values can be given to
5173     this option, each separated by a comma.  This option can also be
5174     called multiple times.
5175
5176     Some measurements by MakeCatalog require a per-pixel sky standard
5177     deviation (for example magnitude error or S/N). Therefore when
5178     asking for such measurements, use the ‘--instd’ option (described
5179     below) to specify the per-pixel sky standard deviation over each
5180     pixel.  For other measurements like the magnitude or surface
5181     brightness, MakeCatalog will need a Zeropoint, which you can set
5182     with the ‘--zeropoint’ option.
5183
5184     For example, by setting ‘--measure=mean,sigclip-mean
5185     --measure=median’, the mean, sigma-clipped mean and median values
5186     will be computed.  The output radial profile will have 4 columns in
5187     this order: radial distance, mean, sigma-clipped and median.  By
5188     default (when this option isn’t given), the mean of all pixels at
5189     each radial position will be computed.
5190
5191‘-s FLT,FLT’
5192‘--sigmaclip=FLT,FLT’
5193     Sigma clipping parameters: only relevant if sigma-clipping
5194     operators are requested by ‘--measure’.  For more on
5195     sigma-clipping, see *note Sigma clipping::.  If given, the value to
5196     this option is directly passed to the ‘--sigmaclip’ option of *note
5197     MakeCatalog::, see *note MakeCatalog inputs and basic settings::.
5198     By default (when this option isn’t given), the default values
5199     within MakeCatalog will be used.  To see the default value of this
5200     option in MakeCatalog, you can run this command:
5201
5202          $ astmkcatalog -P | grep " sigmaclip "
5203
5204‘-z FLT’
5205‘--zeropoint=FLT’
5206     The Zeropoint of the input dataset.  This is necessary when you
5207     request measurements like Magnitude, or Surface brightness.
5208
5209‘-v INT’
5210‘--oversample=INT’
5211     Oversample the input dataset to the fraction given to this option.
5212     Therefore if you set ‘--rmax=20’ for example and ‘--oversample=5’,
5213     your output will have 100 rows (without ‘--oversample’ it will only
5214     have 20 rows).  Unless the object is heavily undersampled (the
5215     pixels are larger than the actual object), this method provides a
5216     much more accurate result and there are sufficient number of pixels
5217     to get the profile accurately.
5218
5219‘-u INT’
5220‘--undersample=INT’
5221     Undersample the input dataset by the number given to this option.
5222     This option is for considering larger apertures than the original
5223     pixel size (aperture size is equal to 1 pixel).  For example, if a
5224     radial profile computed by default has 100 different radii
5225     (apertures of 1 pixel width), by considering ‘--undersample=2’ the
5226     radial profile will be computed over apertures of 2 pixels, so the
5227     final radial profile will have 50 different radii.  This option is
5228     good to measure over a larger number of pixels to improve the
5229     measurement.
5230
5231‘-i FLT/STR5232‘--instd=FLT/STR5233     Sky standard deviation as a single number (FLT) or as the filename
5234     (STR) containing the image with the std value for each pixel (the
5235     HDU within the file should be given to the ‘--stdhdu’ option
5236     mentioned below).  This is only necessary when the requested
5237     measurement (value given to ‘--measure’) by MakeCatalog needs the
5238     Standard deviation (for example the signal-to-noise ratio or
5239     magnitude error).  If your measurements don’t require a standard
5240     deviation, it is best to ignore this option (because it will slow
5241     down the script).
5242
5243‘-d INT/STR5244‘--stdhdu=INT/STR5245     HDU/extension of the sky standard deviation image specified with
5246     ‘--instd’.
5247
5248‘-t STR’
5249‘--tmpdir=STR’
5250     Several intermediate files are necessary to obtain the radial
5251     profile.  All of these temporal files are saved into a temporal
5252     directory.  With this option, you can directly specify this
5253     directory.  By default (when this option isn’t called), it will be
5254     built in the running directory and given an input-based name.  If
5255     the directory doesn’t exist at run-time, this script will create
5256     it.  Once the radial profile has been obtained, this directory is
5257     removed.  You can disable the deletion of the temporary directory
5258     with the ‘--keeptmp’ option.
5259
5260‘-k’
5261‘--keeptmp’
5262     Don’t delete the temporary directory (see description of ‘--tmpdir’
5263     above).  This option is useful for debugging.  For example, to
5264     check that the profiles generated for obtaining the radial profile
5265     have the desired center, shape and orientation.
5266
5267
5268File: gnuastro.info,  Node: SAO DS9 region files from table,  Prev: Generate radial profile,  Up: Installed scripts
5269
527010.3 SAO DS9 region files from table
5271====================================
5272
5273Once your desired catalog (containing the positions of some objects) is
5274created (for example with *note MakeCatalog::, *note Match::, or *note
5275Table::) it often happens that you want to see your selected objects on
5276an image for a feeling of the spatial properties of your objects.  For
5277example you want to see their positions relative to each other.
5278
5279   In this section we describe a simple installed script that is
5280provided within Gnuastro for converting your given columns to an SAO DS9
5281region file to help in this process.  SAO DS9(1) is one of the most
5282common FITS image visualization tools in astronomy and is free software.
5283
5284* Menu:
5285
5286* Invoking astscript-ds9-region::  How to call astscript-ds9-region
5287
5288   ---------- Footnotes ----------
5289
5290   (1) <http://ds9.si.edu>
5291
5292
5293File: gnuastro.info,  Node: Invoking astscript-ds9-region,  Prev: SAO DS9 region files from table,  Up: SAO DS9 region files from table
5294
529510.3.1 Invoking astscript-ds9-region
5296------------------------------------
5297
5298This installed script will read two positional columns within an input
5299table and generate an SAO DS9 region file to visualize the position of
5300the given objects over an image.  For more on installed scripts please
5301see (see *note Installed scripts::).  This script can be used with the
5302following general template:
5303
5304     ## Use the RA and DEC columns of 'table.fits' for the region file.
5305     $ astscript-ds9-region table.fits --column=RA,DEC \
5306                            --output=ds9.reg
5307
5308     ## Select objects with a magnitude between 18 to 20, and generate the
5309     ## region file directly (through a pipe), each region with radius of
5310     ## 0.5 arcseconds.
5311     $ asttable table.fits --range=MAG,18:20 --column=RA,DEC \
5312           | astscript-ds9-region --column=1,2 --radius=0.5
5313
5314     ## With the first command, select objects with a magnitude of 25 to 26
5315     ## as red regions in 'bright.reg'. With the second command, select
5316     ## objects with a magnitude between 28 to 29 as a green region and
5317     ## show both.
5318     $ asttable cat.fits --range=MAG_F160W,25:26 -cRA,DEC \
5319           | ./astscript-ds9-region -c1,2 --color=red -obright.reg
5320     $ asttable cat.fits --range=MAG_F160W,28:29 -cRA,DEC \
5321           | ./astscript-ds9-region -c1,2 --color=green \
5322                         --command="ds9 image.fits -regions bright.reg"
5323
5324   The input can either be passed as a named file, or from standard
5325input (a pipe).  Only the ‘--column’ option is mandatory (to specify the
5326input table columns): two columns from the input table must be
5327specified, either by name (recommended) or number.  You can optionally
5328also specify the region’s radius, width and color of the regions with
5329the ‘--radius’, ‘--width’ and ‘--color’ options, otherwise default
5330values will be used for these (described under each option).
5331
5332   The created region file will be written into the file name given to
5333‘--output’.  When ‘--output’ isn’t called, the default name of ‘ds9.reg5334will be used (in the running directory).  If the file exists before
5335calling this script, it will be overwritten, unless you pass the
5336‘--dontdelete’ option.  Optionally you can also use the ‘--command’
5337option to give the full command that should be run to execute SAO DS9
5338(see example above and description below).  In this mode, the created
5339region file will be deleted once DS9 is closed (unless you pass the
5340‘--dontdelete’ option).  A full description of each option is given
5341below.
5342
5343‘-h INT/STR5344‘--hdu INT/STR5345     The HDU of the input table when a named FITS file is given as
5346     input.  The HDU (or extension) can be either a name or number
5347     (counting from zero).  For more on this option, see *note Input
5348     output options::.
5349
5350‘-c STR,STR’
5351‘--column=STR,STR’
5352     Identifiers of the two positional columns to use in the DS9 region
5353     file from the table.  They can either be in WCS (RA and Dec) or
5354     image (pixel) coordinates.  The mode can be specified with the
5355     ‘--mode’ option, described below.
5356
5357‘-n STR’
5358‘--namecol=STR’
5359     The column containing the name (or label) of each region.  The type
5360     of the column (numeric or a character-based string) is irrelevant:
5361     you can use both types of columns as a name or label for the
5362     region.  This feature is useful when you need to recognize each
5363     region with a certain ID or property (for example magnitude or
5364     redshift).
5365
5366‘-m wcs|img’
5367‘--mode=wcs|org’
5368     The coordinate system of the positional columns (can be either
5369     ‘--mode=wcs’ and ‘--mode=img’).  In the WCS mode, the values within
5370     the columns are interpreted to be RA and Dec.  In the image mode,
5371     they are interpreted to be pixel X and Y positions.  This option
5372     also affects the interpretation of the value given to ‘--radius’.
5373     When this option isn’t explicitly given, the columns are assumed to
5374     be in WCS mode.
5375
5376‘-C STR’
5377‘--color=STR’
5378     The color to use for created regions.  These will be directly
5379     interpreted by SAO DS9 when it wants to open the region file so it
5380     must be recognizable by SAO DS9.  As of SAO DS9 8.2, the recognized
5381     color names are ‘black’, ‘white’, ‘red’, ‘green’, ‘blue’, ‘cyan’,
5382     ‘magenta’ and ‘yellow’.  The default color (when this option is not
5383     called) is ‘green’
5384
5385‘-w INT’
5386‘--width=INT’
5387     The line width of the regions.  These will be directly interpreted
5388     by SAO DS9 when it wants to open the region file so it must be
5389     recognizable by SAO DS9.  The default value is ‘1’.
5390
5391‘-r FLT’
5392‘--radius=FLT’
5393     The radius of all the regions.  In WCS mode, the radius is assumed
5394     to be in arc-seconds, in image mode, it is in pixel units.  If this
5395     option is not explicitly given, in WCS mode the default radius is 1
5396     arc-seconds and in image mode it is 3 pixels.
5397
5398‘--dontdelete’
5399     If the output file name exists, abort the program and don’t
5400     over-write the contents of the file.  This option is thus good if
5401     you want to avoid accidentally writing over an important file.
5402     Also, don’t delete the created region file when ‘--command’ is
5403     given (by default, when ‘--command’ is given, the created region
5404     file will be deleted after SAO DS9 closes).
5405
5406‘-o STR’
5407‘--output=STR’
5408     Write the created SAO DS9 region file into the name given to this
5409     option.  If not explicitly given on the command-line, a default
5410     name of ‘ds9.reg’ will be used.  If the file already exists, it
5411     will be over-written, you can avoid the deletion (or over-writing)
5412     of an existing file with the ‘--dontdelete’.
5413
5414‘--command="STR"’
5415     After creating the region file, run the string given to this option
5416     as a command-line command.  The SAO DS9 region command will be
5417     appended to the end of the given command.  Because the command will
5418     mostly likely contain white-space characters it is recommended to
5419     put the given string in double quotations.
5420
5421     For example, let’s assume ‘--command="ds9 image.fits -zscale"’.
5422     After making the region file (assuming it is called ‘ds9.reg’), the
5423     following command will be executed:
5424
5425          ds9 image.fits -zscale -regions ds9.reg
5426
5427     You can customize all aspects of SAO DS9 with its command-line
5428     options, therefore the value of this option can be as long and
5429     complicated as you like.  For example if you also want the image to
5430     fit into the window, this option will be: ‘--command="ds9
5431     image.fits -zscale -zoom to fit"’.  You can see the SAO DS9
5432     command-line descriptions by clicking on the “Help” menu and
5433     selecting “Reference Manual”.  In the opened window, click on
5434     “Command Line Options”.
5435
5436
5437File: gnuastro.info,  Node: Library,  Next: Developing,  Prev: Installed scripts,  Up: Top
5438
543911 Library
5440**********
5441
5442Each program in Gnuastro that was discussed in the prior chapters (or
5443any program in general) is a collection of functions that is compiled
5444into one executable file which can communicate directly with the outside
5445world.  The outside world in this context is the operating system.  By
5446communication, we mean that control is directly passed to a program from
5447the operating system with a (possible) set of inputs and after it is
5448finished, the program will pass control back to the operating system.
5449For programs written in C and C++, the unique ‘main’ function is in
5450charge of this communication.
5451
5452   Similar to a program, a library is also a collection of functions
5453that is compiled into one executable file.  However, unlike programs,
5454libraries don’t have a ‘main’ function.  Therefore they can’t
5455communicate directly with the outside world.  This gives you the chance
5456to write your own ‘main’ function and call library functions from within
5457it.  After compiling your program into a binary executable, you just
5458have to _link_ it to the library and you are ready to run (execute) your
5459program.  In this way, you can use Gnuastro at a much lower-level, and
5460in combination with other libraries on your system, you can
5461significantly boost your creativity.
5462
5463   This chapter starts with a basic introduction to libraries and how
5464you can use them in *note Review of library fundamentals::.  The
5465separate functions in the Gnuastro library are then introduced
5466(classified by context) in *note Gnuastro library::.  If you end up
5467routinely using a fixed set of library functions, with a well-defined
5468input and output, it will be much more beneficial if you define a
5469program for the job.  Therefore, in its *note Version controlled
5470source::, Gnuastro comes with the *note The TEMPLATE program:: to easily
5471define your own programs(s).
5472
5473* Menu:
5474
5475* Review of library fundamentals::  Guide on libraries and linking.
5476* BuildProgram::                Link and run source files with this library.
5477* Gnuastro library::            Description of all library functions.
5478* Library demo programs::       Demonstration for using libraries.
5479
5480
5481File: gnuastro.info,  Node: Review of library fundamentals,  Next: BuildProgram,  Prev: Library,  Up: Library
5482
548311.1 Review of library fundamentals
5484===================================
5485
5486Gnuastro’s libraries are written in the C programming language.  In
5487*note Why C::, we have thoroughly discussed the reasons behind this
5488choice.  C was actually created to write Unix, thus understanding the
5489way C works can greatly help in effectively using programs and libraries
5490in all Unix-like operating systems.  Therefore, in the following
5491subsections some important aspects of C, as it relates to libraries (and
5492thus programs that depend on them) on Unix are reviewed.  First we will
5493discuss header files in *note Headers:: and then go onto *note
5494Linking::.  This section finishes with *note Summary and example on
5495libraries::.  If you are already familiar with these concepts, please
5496skip this section and go directly to *note Gnuastro library::.
5497
5498   In theory, a full operating system (or any software) can be written
5499as one function.  Such a software would not need any headers or linking
5500(that are discussed in the subsections below).  However, writing that
5501single function and maintaining it (adding new features, fixing bugs,
5502documentation, etc) would be a programmer or scientist’s worst
5503nightmare!  Furthermore, all the hard work that went into creating it
5504cannot be reused in other software: every other programmer or scientist
5505would have to re-invent the wheel.  The ultimate purpose behind
5506libraries (which come with headers and have to be linked) is to address
5507this problem and increase modularity: “the degree to which a system’s
5508components may be separated and recombined” (from Wikipedia).  The more
5509modular the source code of a program or library, the easier maintaining
5510it will be, and all the hard work that went into creating it can be
5511reused for a wider range of problems.
5512
5513* Menu:
5514
5515* Headers::                     Header files included in source.
5516* Linking::                     Linking the compiled source files into one.
5517* Summary and example on libraries::  A summary and example on using libraries.
5518
5519
5520File: gnuastro.info,  Node: Headers,  Next: Linking,  Prev: Review of library fundamentals,  Up: Review of library fundamentals
5521
552211.1.1 Headers
5523--------------
5524
5525C source code is read from top to bottom in the source file, therefore
5526program components (for example variables, data structures and
5527functions) should all be _defined_ or _declared_ closer to the top of
5528the source file: before they are used.  _Defining_ something in C or C++
5529is jargon for providing its full details.  _Declaring_ it, on the
5530other-hand, is jargon for only providing the minimum information needed
5531for the compiler to pass it temporarily and fill in the detailed
5532definition later.
5533
5534   For a function, the _declaration_ only contains the inputs and their
5535data-types along with the output’s type(1).  The _definition_ adds to
5536the declaration by including the exact details of what operations are
5537done to the inputs to generate the output.  As an example, take this
5538simple summation function:
5539
5540     double
5541     sum(double a, double b)
5542     {
5543       return a + b;
5544     }
5545What you see above is the _definition_ of this function: it shows you
5546(and the compiler) exactly what it does to the two ‘double’ type inputs
5547and that the output also has a ‘double’ type.  Note that a function’s
5548internal operations are rarely so simple and short, it can be
5549arbitrarily long and complicated.  This unreasonably short and simple
5550function was chosen here for ease of reading.  The declaration for this
5551function is:
5552
5553     double
5554     sum(double a, double b);
5555
5556You can think of a function’s declaration as a building’s address in the
5557city, and the definition as the building’s complete blueprints.  When
5558the compiler confronts a call to a function during its processing, it
5559doesn’t need to know anything about how the inputs are processed to
5560generate the output.  Just as the postman doesn’t need to know the inner
5561structure of a building when delivering the mail.  The declaration
5562(address) is enough.  Therefore by _declaring_ the functions once at the
5563start of the source files, we don’t have to worry about _defining_ them
5564after they are used.
5565
5566   Even for a simple real-world operation (not a simple summation like
5567above!), you will soon need many functions (for example, some for
5568reading/preparing the inputs, some for the processing, and some for
5569preparing the output).  Although it is technically possible, managing
5570all the necessary functions in one file is not easy and is contrary to
5571the modularity principle (see *note Review of library fundamentals::),
5572for example the functions for preparing the input can be usable in your
5573other projects with a different processing.  Therefore, as we will see
5574later (in *note Linking::), the functions don’t necessarily need to be
5575defined in the source file where they are used.  As long as their
5576definitions are ultimately linked to the final executable, everything
5577will be fine.  For now, it is just important to remember that the
5578functions that are called within one source file must be declared within
5579the source file (declarations are mandatory), but not necessarily
5580defined there.
5581
5582   In the spirit of modularity, it is common to define contextually
5583similar functions in one source file.  For example, in Gnuastro,
5584functions that calculate the median, mean and other statistical
5585functions are defined in ‘lib/statistics.c’, while functions that deal
5586directly with FITS files are defined in ‘lib/fits.c’.
5587
5588   Keeping the definition of similar functions in a separate file
5589greatly helps their management and modularity, but this fact alone
5590doesn’t make things much easier for the caller’s source code: recall
5591that while definitions are optional, declarations are mandatory.  So if
5592this was all, the caller would have to manually copy and paste
5593(_include_) all the declarations from the various source files into the
5594file they are working on now.  To address this problem, programmers have
5595adopted the header file convention: the header file of a source code
5596contains all the declarations that a caller would need to be able to use
5597any of its functions.  For example, in Gnuastro, ‘lib/statistics.c5598(file containing function definitions) comes with
5599lib/gnuastro/statistics.h’ (only containing function declarations).
5600
5601   The discussion above was mainly focused on functions, however, there
5602are many more programming constructs such as pre-processor macros and
5603data structures.  Like functions, they also need to be known to the
5604compiler when it confronts a call to them.  So the header file also
5605contains their definitions or declarations when they are necessary for
5606the functions.
5607
5608   Pre-processor macros (or macros for short) are replaced with their
5609defined value by the pre-processor before compilation.  Conventionally
5610they are written only in capital letters to be easily recognized.  It is
5611just important to understand that the compiler doesn’t see the macros,
5612it sees their fixed values.  So when a header specifies macros you can
5613do your programming without worrying about the actual values.  The
5614standard C types (for example ‘int’, or ‘float’) are very low-level and
5615basic.  We can collect multiple C types into a _structure_ for a
5616higher-level way to keep and pass-along data.  See *note Generic data
5617container:: for some examples of macros and data structures.
5618
5619   The contents in the header need to be _include_d into the caller’s
5620source code with a special pre-processor command: ‘#include
5621<path/to/header.h>’.  As the name suggests, the _pre-processor_ goes
5622through the source code prior to the processor (or compiler).  One of
5623its jobs is to include, or merge, the contents of files that are
5624mentioned with this directive in the source code.  Therefore the
5625compiler sees a single entity containing the contents of the main file
5626and all the included files.  This allows you to include many (sometimes
5627thousands of) declarations into your code with only one line.  Since the
5628headers are also installed with the library into your system, you don’t
5629even need to keep a copy of them for each separate program, making
5630things even more convenient.
5631
5632   Try opening some of the ‘.c’ files in Gnuastro’s ‘lib/’ directory
5633with a text editor to check out the include directives at the start of
5634the file (after the copyright notice).  Let’s take ‘lib/fits.c’ as an
5635example.  You will notice that Gnuastro’s header files (like
5636gnuastro/fits.h’) are indeed within this directory (the ‘fits.h’ file
5637is in the ‘gnuastro/’ directory).  You will notice that files like
5638stdio.h’, or ‘string.h’ are not in this directory (or anywhere within
5639Gnuastro).
5640
5641   On most systems the basic C header files (like ‘stdio.h’ and
5642string.h’ mentioned above) are located in ‘/usr/include/’(2).  Your
5643compiler is configured to automatically search that directory (and
5644possibly others), so you don’t have to explicitly mention these
5645directories.  Go ahead, look into the ‘/usr/include’ directory and find
5646stdio.h’ for example.  When the necessary header files are not in those
5647specific libraries, the pre-processor can also search in places other
5648than the current directory.  You can specify those directories with this
5649pre-processor option(3):
5650
5651‘-I DIR’
5652     “Add the directory ‘DIR’ to the list of directories to be searched
5653     for header files.  Directories named by ’-I’ are searched before
5654     the standard system include directories.  If the directory ‘DIR’ is
5655     a standard system include directory, the option is ignored to
5656     ensure that the default search order for system directories and the
5657     special treatment of system headers are not defeated...” (quoted
5658     from the GNU Compiler Collection manual).  Note that the space
5659     between <I> and the directory is optional and commonly not used.
5660
5661   If the pre-processor can’t find the included files, it will abort
5662with an error.  In fact a common error when building programs that
5663depend on a library is that the compiler doesn’t not know where a
5664library’s header is (see *note Known issues::).  So you have to manually
5665tell the compiler where to look for the library’s headers with the ‘-I’
5666option.  For a small software with one or two source files, this can be
5667done manually (see *note Summary and example on libraries::).  However,
5668to enhance modularity, Gnuastro (and most other bin/libraries) contain
5669many source files, so the compiler is invoked many times(4).  This makes
5670manual addition or modification of this option practically impossible.
5671
5672   To solve this problem, in the GNU build system, there are
5673conventional environment variables for the various kinds of compiler
5674options (or flags).  These environment variables are used in every call
5675to the compiler (they can be empty).  The environment variable used for
5676the C Pre-Processor (or CPP) is ‘CPPFLAGS’.  By giving ‘CPPFLAGS’ a
5677value once, you can be sure that each call to the compiler will be
5678affected.  See *note Known issues:: for an example of how to set this
5679variable at configure time.
5680
5681   As described in *note Installation directory::, you can select the
5682top installation directory of a software using the GNU build system,
5683when you ‘./configure’ it.  All the separate components will be put in
5684their separate sub-directory under that, for example the programs,
5685compiled libraries and library headers will go into ‘$prefix/bin5686(replace ‘$prefix’ with a directory), ‘$prefix/lib’, and
5687‘$prefix/include’ respectively.  For enhanced modularity, libraries that
5688contain diverse collections of functions (like GSL, WCSLIB, and
5689Gnuastro), put their header files in a sub-directory unique to
5690themselves.  For example all Gnuastro’s header files are installed in
5691‘$prefix/include/gnuastro’.  In your source code, you need to keep the
5692library’s sub-directory when including the headers from such libraries,
5693for example ‘#include <gnuastro/fits.h>’(5).  Not all libraries need to
5694follow this convention, for example CFITSIO only has one header
5695(‘fitsio.h’) which is directly installed in ‘$prefix/include’.
5696
5697   ---------- Footnotes ----------
5698
5699   (1) Recall that in C, functions only have one output.
5700
5701   (2) The ‘include/’ directory name is taken from the pre-processor’s
5702‘#include’ directive, which is also the motivation behind the ‘I’ in the
5703‘-I’ option to the pre-processor.
5704
5705   (3) Try running Gnuastro’s ‘make’ and find the directories given to
5706the compiler with the ‘-I’ option.
5707
5708   (4) Nearly every command you see being executed after running ‘make’
5709is one call to the compiler.
5710
5711   (5) the top ‘$prefix/include’ directory is usually known to the
5712compiler
5713
5714
5715File: gnuastro.info,  Node: Linking,  Next: Summary and example on libraries,  Prev: Headers,  Up: Review of library fundamentals
5716
571711.1.2 Linking
5718--------------
5719
5720To enhance modularity, similar functions are defined in one source file
5721(with a ‘.c’ suffix, see *note Headers:: for more).  After running
5722‘make’, each human-readable, ‘.c’ file is translated (or compiled) into
5723a computer-readable “object” file (ending with ‘.o’).  Note that object
5724files are also created when building programs, they aren’t particular to
5725libraries.  Try opening Gnuastro’s ‘lib/’ and ‘bin/progname/’
5726directories after running ‘make’ to see these object files(1).
5727Afterwards, the object files are _linked_ together to create an
5728executable program or a library.
5729
5730   The object files contain the full definition of the functions in the
5731respective ‘.c’ file along with a list of any other function (or
5732generally “symbol”) that is referenced there.  To get a list of those
5733functions you can use the ‘nm’ program which is part of GNU Binutils.
5734For example from the top Gnuastro directory, run:
5735
5736     $ nm bin/arithmetic/arithmetic.o
5737
5738This will print a list of all the functions (more generally, ‘symbols’)
5739that were called within ‘bin/arithmetic/arithmetic.c’ along with some
5740further information (for example a ‘T’ in the second column shows that
5741this function is actually defined here, ‘U’ says that it is undefined
5742here).  Try opening the ‘.c’ file to check some of these functions for
5743your self.  Run ‘info nm’ for more information.
5744
5745   To recap, the _compiler_ created the separate object files mentioned
5746above for each ‘.c’ file.  The _linker_ will then combine all the
5747symbols of the various object files (and libraries) into one program or
5748library.  In the case of Arithmetic (a program) the contents of the
5749object files in ‘bin/arithmetic/’ are copied (and re-ordered) into one
5750final executable file which we can run from the operating system.
5751
5752   There are two ways to _link_ all the necessary symbols: static and
5753dynamic/shared.  When the symbols (computer-readable function
5754definitions in most cases) are copied into the output, it is called
5755_static_ linking.  When the symbols are kept in their original file and
5756only a reference to them is kept in the executable, it is called
5757_dynamic_, or _shared_ linking.
5758
5759   Let’s have a closer look at the executable to understand this better:
5760we’ll assume you have built Gnuastro without any customization and
5761installed Gnuastro into the default ‘/usr/local/’ directory (see *note
5762Installation directory::).  If you tried the ‘nm’ command on one of
5763Arithmetic’s object files above, then with the command below you can
5764confirm that all the functions that were defined in the object file
5765above (had a ‘T’ in the second column) are also defined in the
5766‘astarithmetic’ executable:
5767
5768     $ nm /usr/local/bin/astarithmetic
5769
5770These symbols/function have been statically linked (copied) in the final
5771executable.  But you will notice that there are still many undefined
5772symbols in the executable (those with a ‘U’ in the second column).  One
5773class of such functions are Gnuastro’s own library functions that start
5774with ‘‘gal_’’:
5775
5776     $ nm /usr/local/bin/astarithmetic | grep gal_
5777
5778   These undefined symbols (functions) are present in another file and
5779will be linked to the Arithmetic program every time you run it.
5780Therefore they are known as dynamically _linked_ libraries (2).  As we
5781saw above, static linking is done when the executable is being built.
5782However, when a program is dynamically linked to a library, at
5783build-time, the library’s symbols are only checked with the available
5784libraries: they are not actually copied into the program’s executable.
5785Every time you run the program, the (dynamic) linker will be activated
5786and will try to link the program to the installed library before the
5787program starts.
5788
5789   If you want all the libraries to be statically linked to the
5790executables, you have to tell Libtool (which Gnuastro uses for the
5791linking) to disable shared libraries at configure time(3):
5792
5793     $ configure --disable-shared
5794
5795Try configuring Gnuastro with the command above, then build and install
5796it (as described in *note Quick start::).  Afterwards, check the ‘gal_’
5797symbols in the installed Arithmetic executable like before.  You will
5798see that they are actually copied this time (have a ‘T’ in the second
5799column).  If the second column doesn’t convince you, look at the
5800executable file size with the following command:
5801
5802     $ ls -lh /usr/local/bin/astarithmetic
5803
5804It should be around 4.2 Megabytes with this static linking.  If you
5805configure and build Gnuastro again with shared libraries enabled (which
5806is the default), you will notice that it is roughly 100 Kilobytes!
5807
5808   This huge difference would have been very significant in the old
5809days, but with the roughly Terabyte storage drives commonly in use
5810today, it is negligible.  Fortunately, output file size is not the only
5811benefit of dynamic linking: since it links to the libraries at run-time
5812(rather than build-time), you don’t have to re-build a higher-level
5813program or library when an update comes for one of the lower-level
5814libraries it depends on.  You just install the new low-level library and
5815it will automatically be used/linked next time in the programs that use
5816it.  To be fair, this also creates a few complications(4):
5817
5818   • Reproducibility: Even though your high-level tool has the same
5819     version as before, with the updated library, you might not get the
5820     same results.
5821   • Broken links: if some functions have been changed or removed in the
5822     updated library, then the linker will abort with an error at
5823     run-time.  Therefore you need to re-build your higher-level program
5824     or library.
5825
5826   To see a list of all the shared libraries that are needed for a
5827program or a shared library to run, you can use GNU C library’s ‘ldd’(5)
5828program, for example:
5829
5830     $ ldd /usr/local/bin/astarithmetic
5831
5832   Library file names (in their installation directory) start with a
5833‘lib’ and their ending (suffix) shows if they are static (‘.a’) or
5834dynamic (‘.so’), as described below.  The name of the library is in the
5835middle of these two, for example ‘libgsl.a’ or ‘libgnuastro.a’ (GSL and
5836Gnuastro’s static libraries), and ‘libgsl.so.23.0.0’ or
5837libgnuastro.so.4.0.0’ (GSL and Gnuastro’s shared library, the numbers
5838may be different).
5839
5840   • A static library is known as an archive file and has the ‘.a’
5841     suffix.  A static library is not an executable file.
5842
5843   • A shared library ends with the ‘.so.X.Y.Z’ suffix and is
5844     executable.  The three numbers in the suffix, describe the version
5845     of the shared library.  Shared library versions are defined to
5846     allow multiple versions of a shared library simultaneously on a
5847     system and to help detect possible updates in the library and
5848     programs that depend on it by the linker.
5849
5850     It is very important to mention that this version number is
5851     different from the software version number (see *note Version
5852     numbering::), so do not confuse the two.  See the “Library
5853     interface versions” chapter of GNU Libtool for more.
5854
5855     For each shared library, we also have two symbolic links ending
5856     with ‘.so.X’ and ‘.so’.  They are automatically set by the
5857     installer, but you can change them (point them to another version
5858     of the library) when you have multiple versions of a library on
5859     your system.
5860
5861   Libraries that are built with GNU Libtool (including Gnuastro and its
5862dependencies), build both static and dynamic libraries by default and
5863install them in ‘prefix/lib/’ directory (for more on ‘prefix’, see *note
5864Installation directory::).  In this way, programs depending on the
5865libraries can link with them however they prefer.  See the contents of
5866/usr/local/lib’ with the command below to see both the static and
5867shared libraries available there, along with their executable nature and
5868the symbolic links:
5869
5870     $ ls -l /usr/local/lib/
5871
5872   To link with a library, the linker needs to know where to find the
5873library.  _At compilation time_, these locations can be passed to the
5874linker with two separate options (see *note Summary and example on
5875libraries:: for an example) as described below.  You can see these
5876options and their usage in practice while building Gnuastro (after
5877running ‘make’):
5878
5879‘-L DIR’
5880     Will tell the linker to look into ‘DIR’ for the libraries.  For
5881     example ‘-L/usr/local/lib’, or ‘-L/home/yourname/.local/lib’.  You
5882     can make multiple calls to this option, so the linker looks into
5883     several directories at compilation time.  Note that the space
5884     between <L> and the directory is optional and commonly ignored
5885     (written as ‘-LDIR’).
5886
5887‘-lLIBRARY’
5888     Specify the unique library identifier/name (not containing
5889     directory or shared/dynamic nature) to be linked with the
5890     executable.  As discussed above, library file names have fixed
5891     parts which must not be given to this option.  So ‘-lgsl’ will
5892     guide the linker to either look for ‘libgsl.a’ or ‘libgsl.so5893     (depending on the type of linking it is suppose to do).  You can
5894     link many libraries by repeated calls to this option.
5895
5896     *Very important: * The place of this option on the compiler’s
5897     command matters.  This is often a source of confusion for
5898     beginners, so let’s assume you have asked the linker to link with
5899     library A using this option.  As soon as the linker confronts this
5900     option, it looks into the list of the undefined symbols it has
5901     found until that point and does a search in library A for any of
5902     those symbols.  If any pending undefined symbol is found in library
5903     A, it is used.  After the search in undefined symbols is complete,
5904     the contents of library A are completely discarded from the
5905     linker’s memory.  Therefore, if a later object file or library uses
5906     an unlinked symbol in library A, the linker will abort after it has
5907     finished its search in all the input libraries or object files.
5908
5909     As an example, Gnuastro’s ‘gal_fits_img_read’ function depends on
5910     the ‘fits_read_pix’ function of CFITSIO (specified with
5911     ‘-lcfitsio’, which in turn depends on the cURL library, called with
5912     ‘-lcurl’).  So the proper way to link something that uses this
5913     function is ‘-lgnuastro -lcfitsio -lcurl’.  If instead, you give:
5914     ‘-lcfitsio -lgnuastro’ the linker will complain and abort.  To
5915     avoid such linking complexities when using Gnuastro’s library, we
5916     recommend using *note BuildProgram::.
5917
5918   If you have compiled and linked your program with a dynamic library,
5919then the dynamic linker also needs to know the location of the libraries
5920after building the program: _every time_ the program is run afterwards.
5921Therefore, it may happen that you don’t get any errors when
5922compiling/linking a program, but are unable to run your program because
5923of a failure to find a library.  This happens because the dynamic linker
5924hasn’t found the dynamic library _at run time_.
5925
5926   To find the dynamic libraries at run-time, the linker looks into the
5927paths, or directories, in the ‘LD_LIBRARY_PATH’ environment variable.
5928For a discussion on environment variables, especially search paths like
5929‘LD_LIBRARY_PATH’, and how you can add new directories to them, see
5930*note Installation directory::.
5931
5932   ---------- Footnotes ----------
5933
5934   (1) Gnuastro uses GNU Libtool for portable library creation.  Libtool
5935will also make a ‘.lo’ file for each ‘.c’ file when building libraries
5936(‘.lo’ files are human-readable).
5937
5938   (2) Do not confuse dynamically _linked_ libraries with dynamically
5939_loaded_ libraries.  The former (that is discussed here) are only loaded
5940once at the program startup.  However, the latter can be loaded anytime
5941during the program’s execution, they are also known as plugins.
5942
5943   (3) Libtool is very common and is commonly used.  Therefore, you can
5944use this option to configure on most programs using the GNU build system
5945if you want static linking.
5946
5947   (4) Both of these can be avoided by joining the mailing lists of the
5948lower-level libraries and checking the changes in newer versions before
5949installing them.  Updates that result in such behaviors are generally
5950heavily emphasized in the release notes.
5951
5952   (5) If your operating system is not using the GNU C library, you
5953might need another tool.
5954
5955
5956File: gnuastro.info,  Node: Summary and example on libraries,  Prev: Linking,  Up: Review of library fundamentals
5957
595811.1.3 Summary and example on libraries
5959---------------------------------------
5960
5961After the mostly abstract discussions of *note Headers:: and *note
5962Linking::, we’ll give a small tutorial here.  But before that, let’s
5963recall the general steps of how your source code is prepared, compiled
5964and linked to the libraries it depends on so you can run it:
5965
5966  1. The *pre-processor* includes the header (‘.h’) files into the
5967     function definition (‘.c’) files, expands pre-processor macros.
5968     Generally the pre-processor prepares the human-readable source for
5969     compilation (reviewed in *note Headers::).
5970
5971  2. The *compiler* will translate (compile) the human-readable contents
5972     of each source (merged ‘.c’ and the ‘.h’ files, or generally the
5973     output of the pre-processor) into the computer-readable code of
5974     ‘.o’ files.
5975
5976  3. The *linker* will link the called function definitions from various
5977     compiled files to create one unified object.  When the unified
5978     product has a ‘main’ function, this function is the product’s only
5979     entry point, enabling the operating system or user to directly
5980     interact with it, so the product is a program.  When the product
5981     doesn’t have a ‘main’ function, the linker’s product is a library
5982     and its exported functions can be linked to other executables (it
5983     has many entry points).
5984
5985   The GNU Compiler Collection (or GCC for short) will do all three
5986steps.  So as a first example, from Gnuastro’s source, go to
5987tests/lib/’.  This directory contains the library tests, you can use
5988these as some simple tutorials.  For this demonstration, we will compile
5989and run the ‘arraymanip.c’.  This small program will call Gnuastro
5990library for some simple operations on an array (open it and have a
5991look).  To compile this program, run this command inside the directory
5992containing it.
5993
5994     $ gcc arraymanip.c -lgnuastro -lm -o arraymanip
5995
5996The two ‘-lgnuastro’ and ‘-lm’ options (in this order) tell GCC to first
5997link with the Gnuastro library and then with C’s math library.  The ‘-o’
5998option is used to specify the name of the output executable, without it
5999the output file name will be ‘a.out’ (on most OSs), independent of your
6000input file name(s).
6001
6002   If your top Gnuastro installation directory (let’s call it ‘$prefix’,
6003see *note Installation directory::) is not recognized by GCC, you will
6004get pre-processor errors for unknown header files.  Once you fix it, you
6005will get linker errors for undefined functions.  To fix both, you should
6006run GCC as follows: additionally telling it which directories it can
6007find Gnuastro’s headers and compiled library (see *note Headers:: and
6008*note Linking::):
6009
6010     $ gcc -I$prefix/include -L$prefix/lib arraymanip.c -lgnuastro -lm     \
6011           -o arraymanip
6012
6013This single command has done all the pre-processor, compilation and
6014linker operations.  Therefore no intermediate files (object files in
6015particular) were created, only a single output executable was created.
6016You are now ready to run the program with:
6017
6018     $ ./arraymanip
6019
6020   The Gnuastro functions called by this program only needed to be
6021linked with the C math library.  But if your program needs WCS
6022coordinate transformations, needs to read a FITS file, needs special
6023math operations (which include its linear algebra operations), or you
6024want it to run on multiple CPU threads, you also need to add these
6025libraries in the call to GCC: ‘-lgnuastro -lwcs -lcfitsio -lgsl
6026-lgslcblas -pthread -lm’.  In *note Gnuastro library::, where each
6027function is documented, it is mentioned which libraries (if any) must
6028also be linked when you call a function.  If you feel all these linkings
6029can be confusing, please consider Gnuastro’s *note BuildProgram::
6030program.
6031
6032