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-2015 Centre National de la Recherche Scientifique.
 * written by Nathanael Schaeffer (CNRS, ISTerre, Grenoble, France).
 * 
 * nathanael.schaeffer@ujf-grenoble.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 <fftw3.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 lmax,mmax,nlat,nphi,mres, 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;

        lmax = 5;       nlat = 32;
        mmax = 3;       nphi = 10;
        mres = 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 fftw_malloc is required because we need proper 16-byte alignement.
// allocate spatial fields.
        Sh = (double *) fftw_malloc( NSPAT_ALLOC(shtns) * sizeof(double));
        Th = (double *) fftw_malloc( NSPAT_ALLOC(shtns) * sizeof(double));

// allocate SH representations.
        Slm = (complex double *) fftw_malloc( NLM * sizeof(complex double));
        Tlm = (complex double *) fftw_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);
        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);
        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);
}