SHTns  3.4.6
SHT_example.c

A Fortran example program that performs backward and forward Spherical Harmonic Transforms using SHTns. Compile using:

make SHT_example
/*
* Copyright (c) 2010-2021 Centre National de la Recherche Scientifique.
* written by Nathanael Schaeffer (CNRS, ISTerre, Grenoble, France).
*
* nathanael.schaeffer@univ-grenoble-alpes.fr
*
* This software is governed by the CeCILL license under French law and
* abiding by the rules of distribution of free software. You can use,
* modify and/or redistribute the software under the terms of the CeCILL
* license as circulated by CEA, CNRS and INRIA at the following URL
* "http://www.cecill.info".
*
* The fact that you are presently reading this means that you have had
* knowledge of the CeCILL license and that you accept its terms.
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <math.h>
#include <shtns.h>
void write_vect(char *fn, double *vec, int N);
void write_mx(char *fn, double *mx, int N1, int N2);
int main()
{
shtns_cfg shtns; // handle to a sht transform configuration
long int NLM;
complex double *Slm, *Tlm; // spherical harmonics coefficients (l,m space): complex numbers.
double *Sh, *Th; // real space : theta,phi
long int i,im,lm;
double t;
const int lmax = 5; // maximum degree of spherical harmonics
const int mmax = 3; // maximum order of spherical harmonics
const int mres = 1; // periodicity in phi (1 for full-sphere, 2 for half the sphere, 3 for 1/3, etc...)
const int nlat = 64; // number of points in the latitude direction (constraint: nlat >= lmax+1)
const int nphi = 16; // number of points in the longitude direction (constraint: nphi >= 2*mmax+1)
shtns_verbose(1); // displays informations during initialization.
shtns_use_threads(0); // enable multi-threaded transforms (if supported).
shtns = shtns_init( sht_gauss, lmax, mmax, mres, nlat, nphi );
// shtns = shtns_create(lmax, mmax, mres, sht_orthonormal | SHT_REAL_NORM);
// shtns_set_grid(shtns, sht_gauss, 0.0, nlat, nphi);
NLM = shtns->nlm;
// Memory allocation : the use of shtns_malloc is recommended for proper alignement,
// or to use 'pinned' memory for faster GPU transfers (if CUDA is used).
// Use shtns_free() to free the memory when no more needed.
// allocate spatial fields.
Sh = (double *) shtns_malloc( NSPAT_ALLOC(shtns) * sizeof(double));
Th = (double *) shtns_malloc( NSPAT_ALLOC(shtns) * sizeof(double));
// allocate SH representations.
Slm = (complex double *) shtns_malloc( NLM * sizeof(complex double));
Tlm = (complex double *) shtns_malloc( NLM * sizeof(complex double));
// SH_to_spat
LM_LOOP(shtns, Slm[lm]=0.0; Tlm[lm] = 0.0; ) /* this is the same as :
for (lm=0; lm<shtns->nlm; lm++) {
Slm[lm] = 0.0; Tlm[lm] = 0.0;
} */
// Slm[LM(shtns, 1,1)] = sh11_st(shtns); // access to SH coefficient
Slm[LM(shtns, 2,0)] = 1.0;
// Slm[LiM(shtns, 1,0)] = sh10_ct(shtns);
// Slm[LiM(shtns, 0,0)] = 0.5*sh00_1(shtns);
SH_to_spat(shtns, Slm,Sh);
SHtor_to_spat(shtns, Slm,Th,Sh);
write_vect("ylm",(double *) Slm,NLM*2);
write_mx("spat",Sh,nphi,nlat);
// compute value of SH expansion at a given physical point.
double t2;
SHqst_to_point(shtns, Tlm, Tlm, Slm, shtns->ct[nlat/3], 2.*M_PI/(mres*nphi),&t2,&t2,&t);
printf("check if SH_to_point coincides with SH_to_spat : %f = %f\n",t,Sh[nlat/3]);
printf("ct*st = %f\n",shtns->ct[nlat/3]*shtns->st[nlat/3]);
// check non-linear behaviour
for (im=0;im<nphi;im++) {
for (i=0;i<nlat;i++) {
Sh[im*nlat+i] *= Sh[im*nlat+i];
}
}
spat_to_SH(shtns, Sh, Tlm); // /!\ WARNING! this destroys the spatial data in Sh
write_vect("ylm_nl",(double *) Tlm, NLM*2);
// vector transform
LM_LOOP(shtns, Slm[lm]=0.0; Tlm[lm] = 0.0; )
SHsphtor_to_spat(shtns, Slm,Tlm, Sh,Th); // vector transform
write_mx("spatt",Sh,nphi,nlat);
write_mx("spatp",Th,nphi,nlat);
// spat_to_SH
for (im=0;im<nphi;im++) {
for (i=0;i<nlat;i++) {
Sh[im*nlat+i] = shtns->ct[i]; // use cos(theta) array
}
}
spat_to_SH(shtns, Sh,Slm); // /!\ WARNING! this destroys the spatial data in Sh
write_vect("ylm_v",(double *) Slm,NLM*2);
}
void write_vect(char *fn, double *vec, int N)
{
FILE *fp;
int i;
fp = fopen(fn,"w");
for (i=0;i<N;i++) {
fprintf(fp,"%.6g ",vec[i]);
}
fclose(fp);
}
void write_mx(char *fn, double *mx, int N1, int N2)
{
FILE *fp;
int i,j;
fp = fopen(fn,"w");
for (i=0;i<N1;i++) {
for(j=0;j<N2;j++) {
fprintf(fp,"%.6g ",mx[i*N2+j]);
}
fprintf(fp,"\n");
}
fclose(fp);
}
void * shtns_malloc(size_t bytes)
alloc appropriate memory (pinned for gpu, aligned for avx, ...). Use shtns_free to free it.
Definition: sht_init.c:1279
int shtns_use_threads(int num_threads)
Enables multi-thread transform using OpenMP with num_threads (if available). Returns number of thread...
Definition: sht_init.c:1533
shtns_cfg shtns_init(enum shtns_type flags, int lmax, int mmax, int mres, int nlat, int nphi)
Definition: sht_init.c:1514
void SHqst_to_point(shtns_cfg, cplx *Qlm, cplx *Slm, cplx *Tlm, double cost, double phi, double *vr, double *vt, double *vp)
Evaluate vector SH representation Qlm at physical point defined by cost = cos(theta) and phi.
Definition: sht_func.c:596
void SHtor_to_spat(shtns_cfg, cplx *Tlm, double *Vt, double *Vp)
transform toroidal spherical harmonic coefficients Tlm to the spatial theta and phi components (Vt,...
Definition: sht_com.c:50
void SHsphtor_to_spat(shtns_cfg, cplx *Slm, cplx *Tlm, double *Vt, double *Vp)
transform spheroidal-toroidal spherical harmonic coefficients (Slm,Tlm) to the spatial theta and phi ...
Definition: sht_com.c:34
void SH_to_spat(shtns_cfg shtns, cplx *Qlm, double *Vr)
transform the spherical harmonic coefficients Qlm into its spatial representation Vr.
Definition: sht_com.c:26
shtns.h is the definition file for SHTns : include this file in your source code to use SHTns.
#define NSPAT_ALLOC(shtns)
total number of 'doubles' required for a spatial field (includes FFTW reserved space).
Definition: shtns.h:114
void shtns_verbose(int)
controls output during initialization: 0=no output (default), 1=some output, 2=degug (if compiled in)
Definition: sht_init.c:52
#define LM(shtns, l, m)
LM(shtns, l,m) : macro returning array index for given l and m, corresponding to config shtns.
Definition: shtns.h:100
#define LM_LOOP(shtns, action)
LM_LOOP( shtns, action ) : macro that performs "action" for every (l,m), with lm set,...
Definition: shtns.h:103
@ sht_gauss
use gaussian grid and quadrature. Allows on-the-fly or matrix-based algorithms.
Definition: shtns.h:53
Definition: shtns.h:73
const double *const ct
cos(theta) array (size nlat)
Definition: shtns.h:84
const double *const st
sin(theta) array (size nlat)
Definition: shtns.h:85
const unsigned int nlm
total number of (l,m) spherical harmonics components.
Definition: shtns.h:74