CHarm#
This chapter demonstrates how to perform some basic tasks with
CHarm. For a detailed description of the functions to be used, see the
CHarm (C API) chapter. The source codes of the examples can be found in the
cookbook/c directory.
A single header file of the entire library is available as:
charm/charmf.hfor single precision,charm/charm.hfor double precision, andcharm/charmq.hfor quadruple precision.
You may combine various precisions in a single program, provided that you follow the rules from the CHarm in single and quadruple precision chapter below. It is not recommended though to combine various releases (e.g., single precision from v0.0.0 and double precision form v0.1.0).
The chapters that follow show how to work with CHarm in double, single and quadruple precision, respectively.
CHarm in double precision#
In this section, we assume you have successfully installed CHarm (see Installing) in double precision with OpenMP enabled into the default installation paths, and that you know how to compile a C code and link external libraries. A brief example of the compilation under Linux with GCC is provided in the Compilation on Linux section below.
Compilation on Linux#
Assuming your working directory is cookbook/c, an example compilation
with the static library might look like this:
gcc -fopenmp shcs.c -l:libcharm.a -lfftw3 -lfftw3_omp -lm
Compilation using the shared library:
gcc -fopenmp -Wl,-rpath -Wl,/usr/local/lib shcs.c \
-lcharm -lfftw3 -lfftw3_omp -lm
After a successful compilation, you are ready to execute the compiled binary as
./a.out
You should see the following output:
C( 9, 0) = 2.7671430085300001e-08
S( 9, 0) = 0.0000000000000000e+00
C( 9, 4) = -9.0017922533600003e-09
S( 9, 4) = 1.9466677947499999e-08
C( 9, 9) = -4.7747538613200003e-08
S( 9, 9) = 9.6641284771399995e-08
Degree variance for harmonic degree 0 = 1.0000000000000000e+00
Degree variance for harmonic degree 4 = 2.5183167531865940e-12
Degree variance for harmonic degree 10 = 1.2631494966213168e-13
Difference degree variance for harmonic degree 0 = 0.0000000000000000e+00
Difference degree variance for harmonic degree 4 = 0.0000000000000000e+00
Difference degree variance for harmonic degree 10 = 0.0000000000000000e+00
Great, all done!
Note
Various approaches exist to compile the code using static/shared libraries. Use whatever works best for you.
Working with spherical harmonic coefficients#
This example shows how to read/write spherical harmonic coefficients from/to files and how to compute (difference) degree variances.
#include <stdio.h>
#include <stdlib.h>
#include <charm/charm.h>
int main(void)
{
/* INPUTS */
/* ===================================================================== */
/* Define the path to an input "gfc" text file with spherical harmonic
* coefficients. For details on the structure of the "gfc" file, see the
* description of the "charm_shc_read_gfc" function in the "charm_shc"
* module. */
char shcs_in_file[] = "../../data/input/EGM96-degree10.gfc";
/* Maximum harmonic degree to initialize, read and write spherical harmonic
* coefficients */
unsigned long nmax = 10;
/* Define the path to an output binary file with spherical harmonic
* coefficients. For details on the structure of the output binary file,
* see the description of the "charm_shc_write_bin" function in the
* "charm_shc" module. */
char shcs_out_file[] = "../../data/output/EGM96-degree10.shcs";
/* ===================================================================== */
/* ===================================================================== */
/* Read the spherical harmonic coefficients */
/* --------------------------------------------------------------------- */
/* Initialize a "charm_shc" structure up to degree "nmax". All
* coefficients will be initialized to zero and both the scaling constant
* "shcs->mu" and the radius of the reference sphere "shcs->r" will be set
* to "1.0". These values will later be overwritten when reading the input
* file. In case of failure, returned is "NULL", so always check the
* returned pointer for the "NULL" value! */
charm_shc *shcs = charm_shc_calloc(nmax, 1.0, 1.0);
if (shcs == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_shc\" structure");
exit(CHARM_FAILURE);
}
/* Many CHarm functions take a "charm_err" structure as their last
* parameter. This is to allow the called function to report an error (if
* encountered) back to the caller, including some useful information, such
* as the error message, error code, function name and so on. The error
* structure is defined in the "charm_err" module and can be initialized
* like this: */
charm_err *err = charm_err_init();
if (err == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_err\" structure.\n");
exit(CHARM_FAILURE);
}
/* Now let's read spherical harmonic coefficients, the scaling constant and
* the radius of the reference sphere from the input "gfc" file. We are
* reading a static "gfc" model, so the "epoch" parameter is set to "NULL".
* */
charm_shc_read_gfc(shcs_in_file, nmax, NULL, shcs, err);
/* At this point, "shcs" should store the loaded spherical harmonic
* coefficients. However, before we continue in our program, we should
* check whether the previous function call was successful or not. To this
* end, we shall look into the "err" structure by calling an error handler.
* The first parameter of the error handler is the error structure itself,
* "err". The other parameter says that if there is indeed an error
* message in "err", the program prints details on the error and
* subsequently terminates. If you want to print the error but do not want
* to terminate your program, replace "1" by "0" (further details in the
* "charm_err" module). After you call the error handler, you can reuse
* the same "err" structure in your program, since the function resets the
* error structure to the default empty values.
*
* It is not absolutely necessary to call the error handler, but it is
* really really really recommended to always do so. */
charm_err_handler(err, 1);
/* --------------------------------------------------------------------- */
/* Now print some more or less randomly chosen spherical harmonic
* coefficients */
/* --------------------------------------------------------------------- */
/* Let's start with zonal coefficients */
unsigned long n = 9; /* Harmonic degree */
unsigned long m = 0; /* Harmonic order */
printf("C(%3lu,%3lu) = %0.16e\n", n, m, shcs->c[m][n - m]);
printf("S(%3lu,%3lu) = %0.16e\n", n, m, shcs->s[m][n - m]);
/* Now some tesseral coefficients */
n = 9;
m = 4;
printf("C(%3lu,%3lu) = %0.16e\n", n, m, shcs->c[m][n - m]);
printf("S(%3lu,%3lu) = %0.16e\n", n, m, shcs->s[m][n - m]);
/* And finally some sectorial coefficients */
n = 9;
m = 9;
printf("C(%3lu,%3lu) = %0.16e\n", n, m, shcs->c[m][n - m]);
printf("S(%3lu,%3lu) = %0.16e\n", n, m, shcs->s[m][n - m]);
/* --------------------------------------------------------------------- */
/* Now let's save the coefficients to a binary file and then read them back
* to another structure for spherical harmonic cofficients. Just for the
* fun... */
/* --------------------------------------------------------------------- */
/* Write the coefficients. */
charm_shc_write_bin(shcs, nmax, shcs_out_file, err);
/* Again, we should call the error handler to see whether the previous call
* was successful or not. */
charm_err_handler(err, 1);
/* And now read back the coefficients to a new "charm_shc" structure called
* "shcs2" */
charm_shc *shcs2 = charm_shc_calloc(nmax, 1.0, 1.0);
if (shcs2 == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_shc\" structure.");
exit(CHARM_FAILURE);
}
charm_shc_read_bin(shcs_out_file, nmax, shcs2, err);
charm_err_handler(err, 1);
/* --------------------------------------------------------------------- */
/* Compute degree variances from "shcs" */
/* --------------------------------------------------------------------- */
/* Compute degree variances of the input signal */
double *dv = (double *)malloc((nmax + 1) * sizeof(double));
if (dv == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
charm_shc_dv(shcs, nmax, dv, err);
charm_err_handler(err, 1);
/* Print some degree variances */
n = 0;
printf("\n\n");
printf("Degree variance for harmonic degree %lu = %0.16e\n", n, dv[n]);
n = 4;
printf("Degree variance for harmonic degree %lu = %0.16e\n", n, dv[n]);
n = 10;
printf("Degree variance for harmonic degree %lu = %0.16e\n", n, dv[n]);
/* --------------------------------------------------------------------- */
/* Now check whether "shcs" and "shcs2" contain the same coefficients by
* computing difference degree variances */
/* --------------------------------------------------------------------- */
double *ddv = (double *)malloc((nmax + 1) * sizeof(double));
if (ddv == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
charm_shc_ddv(shcs, shcs2, nmax, ddv, err);
charm_err_handler(err, 1);
/* Print some difference degree variances */
n = 0;
printf("\n\n");
printf("Difference degree variance for harmonic degree %lu = %0.16e\n", n,
ddv[n]);
n = 4;
printf("Difference degree variance for harmonic degree %lu = %0.16e\n", n,
ddv[n]);
n = 10;
printf("Difference degree variance for harmonic degree %lu = %0.16e\n", n,
ddv[n]);
/* --------------------------------------------------------------------- */
/* --------------------------------------------------------------------- */
/* Remember: always free the memory associated with the CHarm structures
* via the special CHarm "charm_*_free*" functions. Calling the usual
* "free" function will not deallocate the memory properly and will lead to
* memory leaks. */
charm_err_free(err);
charm_shc_free(shcs);
charm_shc_free(shcs2);
free(dv), free(ddv);
/* --------------------------------------------------------------------- */
printf("\nGreat, all done!\n");
exit(CHARM_SUCCESS);
/* ===================================================================== */
}
Spherical harmonic synthesis and analysis#
This section shows several examples.
A closed-loop test of analysis and synthesis with point data values. First, spherical harmonic coefficients are loaded from a text file. Then, they are used to synthesize a signal, from which a new set of coefficients is finally computed by surface spherical harmonic analysis. The two coefficients sets are then finally compared (should be equal, that is, the difference degree amplitudes should be negligibly small).
A solid spherical harmonic synthesis of point values at scattered points.
A solid spherical harmonic synthesis of point values at a custom grid of points.
A solid spherical harmonic synthesis of block-mean values at scattered cells.
A solid spherical harmonic synthesis of block mean values at a grid of cells.
Surface spherical harmonic analysis with block-mean values in cells.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <charm/charm.h>
int main(void)
{
/* INPUTS */
/* ===================================================================== */
/* Define the path to an input text file with spherical harmonic
* coefficients. For details on the structure of the text file, see the
* description of the "charm_shc_read_mtx" function in the "charm_shc"
* module. */
char shcs_in_file[] = "../../data/input/EGM96-degree10-mtx.txt";
/* Maximum harmonic degree of coefficients to read. The same degree is
* used later in the harmonic synthesis and harmonic analysis. */
unsigned long nmax = 10;
/* ===================================================================== */
/* ===================================================================== */
printf("Closed-loop experiment -- Spherical harmonic synthesis and "
"analysis\n");
printf("===========================\n");
/* Initialize a "charm_shc" structure */
charm_shc *shcs = charm_shc_calloc(nmax, 1.0, 1.0);
if (shcs == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_shc\" structure.\n");
exit(CHARM_FAILURE);
}
/* Initialize a "charm_err" structure */
charm_err *err = charm_err_init();
if (err == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_err\" structure.\n");
exit(CHARM_FAILURE);
}
/* Read spherical harmonic coefficients from the input text file. */
charm_shc_read_mtx(shcs_in_file, nmax, shcs, err);
charm_err_handler(err, 1);
/* Now let's say we do not want to use the zero-degree term. This can be
* achieved simply by setting the respective "C00" coefficient to zero. */
shcs->c[0][0] = 0.0;
/* Compute the Gauss--Legendre grid for a given "nmax" on a sphere that is
* "1000.0" metres above the reference sphere of spherical harmonic
* coefficients. */
charm_point *grd_pnt = charm_crd_point_gl(nmax, shcs->r + 1000.0);
if (grd_pnt == NULL)
{
fprintf(stderr, "Failed to compute the Gauss--Legendre grid.\n");
exit(CHARM_FAILURE);
}
/* Allocate the memory for the synthesized signal */
double *f = (double *)malloc(grd_pnt->npoint * sizeof(double));
if (f == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
/* Perform the synthesis. Since the Gauss--Legendre grid resides "1000.0"
* metres above the reference sphere of the coefficients, performed is
* solid spherical harmonic synthesis. */
charm_shs_point(grd_pnt, shcs, nmax, f, err);
charm_err_handler(err, 1);
/* Print some synthesized values */
printf("Print some synthesized values of the signal...\n");
size_t i = 0;
size_t j = 0;
printf("f(%zu, %zu) = %0.16e\n", i, j, f[i * grd_pnt->nlon + j]);
i = grd_pnt->nlat / 2;
j = grd_pnt->nlon / 2;
printf("f(%zu, %zu) = %0.16e\n", i, j, f[i * grd_pnt->nlon + j]);
/* Initialize a new structure of spherical harmonic coefficients for
* coefficients to be computed by the harmonic analysis. Used are the same
* scalling constants as in "shcs". */
charm_shc *shcs2 = charm_shc_calloc(nmax, shcs->mu, shcs->r);
if (shcs2 == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_shc\" structure.\n");
exit(CHARM_FAILURE);
}
/* Now use the synthesized signal and compute back its spherical harmonic
* coefficients by harmonic analysis. The output coefficients in "shcs2"
* should be the same as the input ones "shcs". Note that the signal "f"
* resides on a sphere that is "1000.0" metres above the reference sphere.
* Therefore, the coefficients are automatically properly rescaled to the
* desired sphere defined by "shcs2->r". */
charm_sha_point(grd_pnt, f, nmax, shcs2, err);
charm_err_handler(err, 1);
/* Now check whether "shcs" and "shcs2" are the same by computing their
* difference degree amplitudes */
double *dda = (double *)malloc((nmax + 1) * sizeof(double));
if (dda == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
charm_shc_dda(shcs, shcs2, nmax, dda, err);
charm_err_handler(err, 1);
/* Print some difference degree amplitudes */
unsigned long n = 0;
printf("\n\n");
printf("Now print the difference degree amplitudes. ");
printf("These should be very small, say, at the order of 1e-18 "
"or less...\n");
printf("Difference degree amplitude for harmonic degree %lu = %0.16e\n",
n, dda[n]);
n = 4;
printf("Difference degree amplitude for harmonic degree %lu = %0.16e\n",
n, dda[n]);
n = 10;
printf("Difference degree amplitude for harmonic degree %lu = %0.16e\n",
n, dda[n]);
/* We will not need some of the structures and arrays anymore, so let's
* free them. */
charm_crd_point_free(grd_pnt);
charm_shc_free(shcs2);
free(f);
free(dda);
printf("===========================\n\n");
/* ===================================================================== */
/* ===================================================================== */
printf("\n");
printf("Solid spherical harmonic synthesis at scattered points\n");
printf("===========================\n");
/* Let's create a "charm_point" structure to store 3 scattered points. */
charm_point *sctr_pnt = charm_crd_point_calloc(CHARM_CRD_POINT_SCATTERED,
3, 3);
if (sctr_pnt == NULL)
{
fprintf(stderr,
"Failed to initialize the \"charm_point\" structure.\n");
exit(CHARM_FAILURE);
}
/* Now let's feed "sctr_pnt" with some computation points (angular values
* must be provided in radians). */
sctr_pnt->lat[0] = 0.1;
sctr_pnt->lat[1] = 0.436231;
sctr_pnt->lat[2] = -0.9651;
sctr_pnt->lon[0] = 0.0;
sctr_pnt->lon[1] = 1.53434;
sctr_pnt->lon[2] = 4.2316;
sctr_pnt->r[0] = shcs->r + 1000.0;
sctr_pnt->r[1] = shcs->r + 2000.0;
sctr_pnt->r[2] = shcs->r + 3000.0;
/* Allocate memory to store the synthesized signal */
f = (double *)malloc(sctr_pnt->npoint * sizeof(double));
if (f == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
/* Do the synthesis */
charm_shs_point(sctr_pnt, shcs, nmax, f, err);
charm_err_handler(err, 1);
/* It's really this easy! */
/* Print the synthesized values */
for (size_t i = 0; i < sctr_pnt->npoint; i++)
printf("%0.16e\n", f[i]);
charm_crd_point_free(sctr_pnt);
free(f);
printf("===========================\n\n");
/* ===================================================================== */
/* ===================================================================== */
printf("\n");
printf("Solid spherical harmonic synthesis at a custom grid of points\n");
printf("===========================\n");
/* Initialize a charm_point structure to hold a "5 x 10" point grid. */
grd_pnt = charm_crd_point_calloc(CHARM_CRD_POINT_GRID, 5, 10);
if (grd_pnt == NULL)
{
fprintf(stderr,
"Failed to initialize the \"charm_point\" structure.\n");
exit(CHARM_FAILURE);
}
/* Define some grid points */
for (size_t i = 0; i < grd_pnt->nlat; i++)
{
grd_pnt->lat[i] = M_PI_2 - ((double)i / (double)grd_pnt->nlat) * M_PI;
grd_pnt->r[i] = shcs->r + (double)i;
}
for (size_t j = 0; j < grd_pnt->nlon; j++)
{
grd_pnt->lon[j] = ((double)j / (double)grd_pnt->nlon) * (2.0 * M_PI);
}
/* Initialize an array to store the synthesized signal */
f = (double *)malloc(grd_pnt->npoint * sizeof(double));
if (f == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
/* Do the synthesis */
charm_shs_point(grd_pnt, shcs, nmax, f, err);
charm_err_handler(err, 1);
/* Print the synthesized values */
for (size_t i = 0; i < grd_pnt->nlat; i++)
for (size_t j = 0; j < grd_pnt->nlon; j++)
printf("%0.16e\n", f[i * grd_pnt->nlon + j]);
charm_crd_point_free(grd_pnt);
free(f);
printf("===========================\n");
/* ===================================================================== */
/* ===================================================================== */
printf("\n");
printf("Solid spherical harmonic synthesis at scattered cells\n");
printf("===========================\n");
/* Initialize a "charm_cell" structure to hold 3 scattered cells. */
charm_cell *sctr_cell = charm_crd_cell_calloc(CHARM_CRD_CELL_SCATTERED,
3, 3);
if (sctr_cell == NULL)
{
fprintf(stderr, "Failed to initialize scattered cells.\n");
exit(CHARM_FAILURE);
}
/* Define more or less randomly some scattered cells */
/* --------------------------------------------------------------------- */
/* The maximum latitude of the first cell */
sctr_cell->latmax[0] = 0.323413435;
/* The minimum latitude of the first cell */
sctr_cell->latmin[0] = sctr_cell->latmax[0] - 0.234323;
/* The maximum latitude of the second cell */
sctr_cell->latmax[1] = -0.90234320952;
/* The minimum latitude of the second cell */
sctr_cell->latmin[1] = sctr_cell->latmax[1] - 0.4456;
/* And so on... */
sctr_cell->latmax[2] = 0.0;
sctr_cell->latmin[2] = sctr_cell->latmax[2] - M_PI_2;
/* And now the same, but with the longitudes... */
sctr_cell->lonmin[0] = 0.123456789;
sctr_cell->lonmax[0] = sctr_cell->lonmin[0] + 1.3235;
sctr_cell->lonmin[1] = 4.3445234;
sctr_cell->lonmax[1] = sctr_cell->lonmin[1] + 0.01;
sctr_cell->lonmin[2] = 0.0;
sctr_cell->lonmax[2] = sctr_cell->lonmin[2] + 2.0 * M_PI;
/* Finally, we define spherical radii of the scattered cells. Importantly,
* the spherical radius is *constants* over the cells (but may vary from
* cell to cell). */
sctr_cell->r[0] = shcs->r;
sctr_cell->r[1] = shcs->r + 1000.0;
sctr_cell->r[2] = shcs->r + 2000.0;
/* --------------------------------------------------------------------- */
/* Initialize an array to store the synthesized signal */
f = (double *)malloc(sctr_cell->ncell * sizeof(double));
if (f == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
/* Do the synthesis */
charm_shs_cell(sctr_cell, shcs, nmax, f, err);
charm_err_handler(err, 1);
/* Print the synthesized values */
for (size_t i = 0; i < sctr_cell->ncell; i++)
printf("%0.16e\n", f[i]);
charm_crd_cell_free(sctr_cell);
free(f);
printf("===========================\n\n");
/* ===================================================================== */
/* ===================================================================== */
printf("\n");
printf("Solid spherical harmonic synthesis at a grid "
"of cells\n");
printf("===========================\n");
/* Initialize a "charm_cell" structure to hold a grid of cells with "15"
* cells in the latitudinal direction and "30" cells in the longitudinal
* direction. */
charm_cell *grd_cell = charm_crd_cell_calloc(CHARM_CRD_CELL_GRID, 15, 30);
/* Define some grid cells */
for (size_t i = 0; i < grd_cell->nlat; i++)
{
/* The maximum cell latitudes */
grd_cell->latmax[i] = M_PI_2 -
((double)i / (double)grd_cell->nlat) * M_PI;
/* The minimum cell latitudes */
grd_cell->latmin[i] = M_PI_2 -
((double)(i + 1) / (double)grd_cell->nlat) *
M_PI;
/* The spherical radii (may vary with the latitude, but we use
* a constant value here, because of the example that follows after
* this one). */
grd_cell->r[i] = shcs->r + 1000.0;
}
for (size_t j = 0; j < grd_cell->nlon; j++)
{
/* The minimum cell longitudes */
grd_cell->lonmin[j] = ((double)j / (double)grd_cell->nlon) *
(2.0 * M_PI);
/* The maximum cell longitudes */
grd_cell->lonmax[j] = ((double)(j + 1) / (double)grd_cell->nlon) *
(2.0 * M_PI);
}
/* Initialize an array to store the synthesized signal */
f = (double *)malloc(grd_cell->ncell * sizeof(double));
if (f == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
/* Do the synthesis */
charm_shs_cell(grd_cell, shcs, nmax, f, err);
charm_err_handler(err, 1);
/* Print some of the synthesized values */
/* --------------------------------------------------------------------- */
i = 0;
j = 10;
printf("f(%zu, %zu) = %0.16e\n", i, j, f[i * grd_cell->nlon + j]);
i = 4;
j = 3;
printf("f(%zu, %zu) = %0.16e\n", i, j, f[i * grd_cell->nlon + j]);
/* --------------------------------------------------------------------- */
printf("===========================\n\n");
/* ===================================================================== */
/* ===================================================================== */
printf("\n");
printf("Spherical harmonic analysis of block-mean values in cells\n");
printf("===========================\n");
/* In this example, we use the "grd_cell" structure and the signal "f" from
* the previous example. */
/* Initialize a new structure of spherical harmonic coefficients for
* coefficients to be recovered from the harmonic analysis. Used are the
* same scalling constants as in "shcs". */
shcs2 = charm_shc_calloc(nmax, shcs->mu, shcs->r);
if (shcs2 == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_shc\" structure.\n");
exit(CHARM_FAILURE);
}
/* Do the analysis using the method of the approximate quadrature */
charm_sha_cell(grd_cell, f, nmax, CHARM_SHA_CELL_AQ, shcs2, err);
charm_err_handler(err, 1);
/* Print some of the computed coefficients. Note that the harmonic
* analysis with block-mean values in cells in *not* exact, hence the
* coefficients will not be equal to the input ones. */
/* --------------------------------------------------------------------- */
n = 2;
unsigned long m = 0;
printf("C(%3lu,%3lu) = %0.16e\n", n, m, shcs2->c[m][n - m]);
printf("S(%3lu,%3lu) = %0.16e\n", n, m, shcs2->s[m][n - m]);
n = 9;
m = 0;
printf("C(%3lu,%3lu) = %0.16e\n", n, m, shcs2->c[m][n - m]);
printf("S(%3lu,%3lu) = %0.16e\n", n, m, shcs2->s[m][n - m]);
n = 9;
m = 4;
printf("C(%3lu,%3lu) = %0.16e\n", n, m, shcs2->c[m][n - m]);
printf("S(%3lu,%3lu) = %0.16e\n", n, m, shcs2->s[m][n - m]);
n = 9;
m = 9;
printf("C(%3lu,%3lu) = %0.16e\n", n, m, shcs2->c[m][n - m]);
printf("S(%3lu,%3lu) = %0.16e\n", n, m, shcs2->s[m][n - m]);
/* --------------------------------------------------------------------- */
charm_crd_cell_free(grd_cell);
charm_shc_free(shcs2);
free(f);
printf("===========================\n\n");
/* ===================================================================== */
/* Free the heap memory */
/* ===================================================================== */
charm_err_free(err);
charm_shc_free(shcs);
printf("\nGreat, all done!\n");
exit(CHARM_SUCCESS);
/* ===================================================================== */
}
First- and second-order gradients (gravitational vector and tensor)#
This example shows how to compute the full gravitational vector and the full gravitational tensor in the Local north-oriented reference frame.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <charm/charm.h>
int main(void)
{
/* INPUTS */
/* ===================================================================== */
/* Path to the input file with spherical harmonic coefficients in the "gfc"
* format. */
char shcs_file[] = "../../data/input/EGM96-degree10.gfc";
/* ===================================================================== */
/* ===================================================================== */
/* Initialize a "charm_err" structure */
charm_err *err = charm_err_init();
if (err == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_err\" structure.\n");
exit(CHARM_FAILURE);
}
/* First, let's get the maximum harmonic degree stored in "shcs_file".
* This is especially useful to read all coefficients in "shcs_file"
* without knowing its maximum harmonic degree a priori. */
unsigned long nmax = charm_shc_read_gfc(shcs_file, CHARM_SHC_NMAX_MODEL,
NULL, NULL, err);
charm_err_handler(err, 1);
/* Allocate memory for the spherical harmonic coefficients up to degree
* "nmax" */
charm_shc *shcs = charm_shc_malloc(nmax, 1.0, 1.0);
if (shcs == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_shc\" structure.\n");
exit(CHARM_FAILURE);
}
/* Now read all coefficients in "shcs_file" */
charm_shc_read_gfc(shcs_file, nmax, NULL, shcs, err);
charm_err_handler(err, 1);
/* Compute the Gauss--Legendre point grid for an "nmax" given by the
* maximum degree stored in "shcs" and for a radius equal to the reference
* sphere, to which the spherical harmonic coefficients are scaled to. We
* intentionally use here the Gauss--Legendre grid instead of the
* Driscoll--Healy grids in order to avoid the inaccuracies due to the
* singularities at the poles (see the documentation). */
charm_point *grd = charm_crd_point_gl(shcs->nmax, shcs->r);
if (grd == NULL)
{
fprintf(stderr, "Failed to compute the Gauss--Legendre grid.\n");
exit(CHARM_FAILURE);
}
/* Total number of gravitational vector elements */
size_t np = 3;
/* Allocate memory for the gravitational vector elements */
double **f = (double **)malloc(np * sizeof(double *));
if (f == NULL)
{
fprintf(stderr, "malloc failure\n");
exit(CHARM_FAILURE);
}
for (size_t i = 0; i < np; i++)
{
f[i] = (double *)malloc(grd->npoint * sizeof(double));
if (f[i] == NULL)
{
fprintf(stderr, "malloc failure\n");
exit(CHARM_FAILURE);
}
}
/* Compute the full first-order gradient (vector) */
charm_shs_point_grad1(grd, shcs, nmax, f, err);
charm_err_handler(err, 1);
/* Allocate memory for the magnitude of the gravitational vector */
double *g = (double *)malloc(grd->npoint * sizeof(double));
if (g == NULL)
{
fprintf(stderr, "malloc failure\n");
exit(CHARM_FAILURE);
}
/* Compute the magnitude of the gravitational acceleration (no contribution
* due to the centrifugal force is considered here) */
for (size_t i = 0; i < grd->npoint; i++)
g[i] = sqrt(f[0][i] * f[0][i] + f[1][i] * f[1][i] + f[2][i] * f[2][i]);
/* Free the memory associated with the gravitational vector */
for (size_t i = 0; i < np; i++)
free(f[i]);
free(f);
/* Allocate memory for the gravitational tensor elements */
np = 6; /* There are six tensor elements */
f = (double **)malloc(np * sizeof(double *));
if (f == NULL)
{
fprintf(stderr, "malloc failure\n");
exit(CHARM_FAILURE);
}
for (size_t i = 0; i < np; i++)
{
f[i] = (double *)malloc(grd->npoint * sizeof(double));
if (f[i] == NULL)
{
fprintf(stderr, "malloc failure\n");
exit(CHARM_FAILURE);
}
}
/* Now compute the full second-order gradient (tensor) */
charm_shs_point_grad2(grd, shcs, nmax, f, err);
charm_err_handler(err, 1);
/* Check whether "fxx + fyy + fzz" is zero within numerical errors */
double trace_error = fabs(f[0][0] + f[3][0] + f[5][0]); /* Trace at the
* first grid
* point of "grd"
* */
double tmp;
for (size_t i = 1; i < grd->npoint; i++)
tmp = fabs(f[0][i] + f[3][i] + f[5][i]);
if (tmp > trace_error)
trace_error = tmp;
printf("The largest error of the gravitational tensor trace is "
"%0.16e s^-2.\n", trace_error);
/* Free the heap memory */
charm_shc_free(shcs);
charm_crd_point_free(grd);
charm_err_free(err);
for (size_t i = 0; i < np; i++)
free(f[i]);
free(f);
free(g);
printf("Great, all done!\n");
/* ===================================================================== */
exit(CHARM_SUCCESS);
}
Fourier coefficients of Legendre functions#
#include <stdio.h>
#include <stdlib.h>
#include <charm/charm.h>
int main(void)
{
/* Maximum harmonic degree to compute the Fourier coefficients of Legendre
* functions */
unsigned long nmax = 500;
/* Initialize a structure to store the Fourier coefficients of Legendre
* functions */
charm_pnmj *pnmj = charm_leg_pnmj_calloc(nmax, CHARM_LEG_PMNJ);
if (pnmj == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_pnmj\" "
"structure.\n");
exit(CHARM_FAILURE);
}
/* Initialize an "charm_err" structure */
charm_err *err = charm_err_init();
if (err == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_err\" structure.\n");
exit(CHARM_FAILURE);
}
/* Compute the Fourier coefficients */
charm_leg_pnmj_coeffs(pnmj, nmax, err);
charm_err_handler(err, 1);
/* Print some Fourier coefficients */
/* Harmonic degree */
unsigned long n = 123;
/* Harmonic order */
unsigned long m = 23;
/* Wave-number-related variable */
unsigned long j = 12;
/* Wave-number (computed from "j"; for details, see the "charm_leg"
* module)*/
unsigned long k = charm_leg_pnmj_j2k(n, j);
printf("Fourier coefficient for degree %lu, order %lu and "
"wave-number %lu = %0.16e\n", n, m, k, pnmj->pnmj[m][n - m][j]);
n = 360; m = 358; j = 101; k = charm_leg_pnmj_j2k(n, j);
printf("Fourier coefficient for degree %lu, order %lu and "
"wave-number %lu = %0.16e\n", n, m, k, pnmj->pnmj[m][n - m][j]);
/* Free the heap memory */
charm_leg_pnmj_free(pnmj);
charm_err_free(err);
printf("\nGreat, all done!\n");
exit(CHARM_SUCCESS);
}
Integrals#
Finally, let’s compute an integral of a product of two spherical harmonics over a restricted domain on the unit sphere. Then, we do the same, but with the product of Legendre functions only.
#include <stdio.h>
#include <stdlib.h>
#include <charm/charm.h>
int main(void)
{
/* Type of the first spherical harmonic function, its harmonic degree and
* order, respectively */
_Bool i1 = 0; /* The "cos" spherical harmonic function */
unsigned long n1 = 287;
unsigned long m1 = 122;
/* Type of the second spherical harmonic function, its harmonic degree and
* order, respectively */
_Bool i2 = 1; /* The "sin" spherical harmonic function */
unsigned long n2 = 34;
unsigned long m2 = 9;
/* Minimum and maximum co-latitudes of the integration domain */
double cltmin = 0.1;
double cltmax = 0.9;
/* Minimum and maximum longitudes of the integration domain */
double lonmin = 0.5;
double lonmax = 0.6;
/* Initialize the Fourier coefficients of Legendre functions */
unsigned long nmax = CHARM_MAX(n1, n2);
charm_pnmj *pnmj = charm_leg_pnmj_calloc(nmax, CHARM_LEG_PMNJ);
if (pnmj == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_pnmj\" "
"structure.\n");
exit(CHARM_FAILURE);
}
/* Initialize an "charm_err" structure */
charm_err *err = charm_err_init();
if (err == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_err\" structure.\n");
exit(CHARM_FAILURE);
}
/* Compute the Fourier coefficients */
charm_leg_pnmj_coeffs(pnmj, nmax, err);
charm_err_handler(err, 1);
/* Compute the integral of a product of two spherical harmonics */
double iy = charm_integ_yi1n1m1yi2n2m2(cltmin, cltmax, lonmin, lonmax,
i1, n1, m1,
i2, n2, m2, pnmj, err);
charm_err_handler(err, 1);
/* Print the value of the integral */
printf("The integral of the product of two spherical harmonics "
"i1 = %d, n1 = %lu, m1 = %lu, i2 = %d, n2 = %lu, m2 = %lu "
"is %0.16e\n", i1, n1, m1, i2, n2, m2, iy);
/* Now compute the integral of a product of two Legendre functions over a
* restricted domain */
double ip = charm_integ_pn1m1pn2m2(cltmin, cltmax, n1, m1, n2, m2, pnmj,
err);
charm_err_handler(err, 1);
printf("The integral of the product of two Legendre functions "
"n1 = %lu, m1 = %lu, n2 = %lu, m2 = %lu "
"is %0.16e\n", n1, m1, n2, m2, ip);
/* Free the heap memory */
charm_err_free(err);
charm_leg_pnmj_free(pnmj);
exit(CHARM_SUCCESS);
}
CHarm in single and quadruple precision#
A few simple rules need to be obeyed in single and quadruple precision.
Make sure you have compiled CHarm in single or quadruple precision (see Installing), whichever you need.
Replace every occurrence of the
charm_*prefix withcharmf_*for single precision andcharmq_*for quadruple precision. This applies to your own code, including the header files, and the compilation of your program.If the prefix is written in capital letters, do not alter it. The
CHARM_*symbols are common across all precisions of CHarm.In your code, replace every
doublevariable entering CHarm functions byfloatfor single precision or by__float128for quadruple precision.To compile your code in quadruple precision, link the
libquadmathlibrary, that is, add-lquadmathbefore-lm(see Compilation on Linux).When compiling your program, use the
fftwf*orfftwq*prefix to link FFTW in single or quadruple precision, respectively.
Example compilation in single precision#
Assuming your working directory is cookbook/c, an example compilation
with the static library might look like this:
gcc -fopenmp shcsf.c -l:libcharmf.a -lfftw3f -lfftw3f_omp -lm
Compilation with the shared library:
gcc -fopenmp -Wl,-rpath -Wl,/usr/local/lib shcsf.c \
-lcharmf -lfftw3f -lfftw3f_omp -lm
Example compilation in quadruple precision#
Static library:
gcc -fopenmp shcsq.c -l:libcharmq.a -lfftw3q -lfftw3q_omp \
-lquadmath -lm
Shared library:
gcc -fopenmp -Wl,-rpath -Wl,/usr/local/lib shcsq.c \
-lcharmq -lfftw3q -lfftw3q_omp -lquadmath -lm