Aquila-Consortium/libsharp2
0
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Experimental: add a reduced Gauss-Legendre grid, which has

Browse source 
the minimum amount of pixels in phi direction (depending on theta)
that still allows SHTs at machine accuracy.
This commit is contained in:
Martin Reinecke 2013-01-11 13:42:44 +01:00
parent 0dd6e2a858
commit d6d008b4b3
1 changed files with 59 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

59
libsharp/sharp_testsuite.c
View file

@ -188,6 +188,32 @@ static void get_errors (dcmplx **alm, ptrdiff_t nalms, int ncomp, double *sqsum,
}
}
static int good_fft_size(int n)
{
if (n<=6) return n;
int bestfac=2*n;
for (int f2=1; f2<bestfac; f2*=2)
for (int f23=f2; f23<bestfac; f23*=3)
for (int f235=f23; f235<bestfac; f235*=5)
if (f235>=n) bestfac=f235;
return bestfac;
}
static int sharp_get_mlim (int lmax, int spin, double sth, double cth,
double ofs)
{
double b = -2*spin*fabs(cth);
double t1 = lmax*sth+ofs;
double c = spin*spin-t1*t1;
double discr = b*b-4*c;
if (discr<=0) return lmax;
double res=(-b+sqrt(discr))/2.;
if (res>lmax) res=lmax;
return (int)(res+0.5);
}
static void get_infos (const char *gname, int lmax, int mmax, int gpar1,
int gpar2, sharp_geom_info **ginfo, sharp_alm_info **ainfo)
{
@ -242,6 +268,39 @@ static void get_infos (const char *gname, int lmax, int mmax, int gpar1,
sharp_make_fejer1_geom_info (nlat, nlon, 0., 1, nlon, ginfo);
if (mytask==0) printf("Clenshaw-Curtis grid, nlat=%d, nlon=%d\n",nlat,nlon);
}
else if (strcmp(gname,"smallgauss")==0)
{
int nlat=gpar1, nlon=gpar2;
if (nlat<1) nlat=lmax+1;
if (nlon<1) nlon=2*mmax+1;
sharp_make_gauss_geom_info (nlat, nlon, 0., 1, nlon, ginfo);
ptrdiff_t npix_o=get_npix(*ginfo);
size_t ofs=0;
for (int i=0; i<(*ginfo)->npairs; ++i)
{
sharp_ringpair *pair=&((*ginfo)->pair[i]);
int pring=1+2*sharp_get_mlim(lmax,0,pair->r1.sth,pair->r1.cth,100.);
if (pring>nlon) pring=nlon;
pring=good_fft_size(pring);
pair->r1.nph=pring;
pair->r1.weight*=nlon*1./pring;
pair->r1.ofs=ofs;
ofs+=pring;
if (pair->r2.nph>0)
{
pair->r2.nph=pring;
pair->r2.weight*=nlon*1./pring;
pair->r2.ofs=ofs;
ofs+=pring;
}
}
if (mytask==0)
{
ptrdiff_t npix=get_npix(*ginfo);
printf("Small Gauss grid, nlat=%d, npix=%ld, savings=%.2f%%\n",
nlat,(long)npix,(npix_o-npix)*100./npix_o);
}
}
else
UTIL_FAIL("unknown grid geometry");