Commit b6608a78 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Added gaussian_latitudes.c.

parent b1ac136b
Pipeline #5358 passed with stages
in 15 minutes and 47 seconds
......@@ -48,6 +48,7 @@
#include "cdo_options.h"
#include "mpmo_color.h"
#include "cdi_lockedIO.h"
#include "gaussian_latitudes.h"
int scan_par_obsolete(char *namelist, const char *name, int def);
int scan_par(int verbose, char *namelist, const char *name, int def);
......@@ -1018,7 +1019,7 @@ after_defineGrid(const AfterControl &globs, struct Variable *vars)
const size_t nlat = globs.Latitudes;
std::vector<double> lats(nlat), latw(nlat);
gaussianLatitudes(lats.data(), latw.data(), nlat);
gaussian_latitudes(nlat, lats.data(), latw.data());
for (size_t j = 0; j < nlat; ++j) lats[j] = 180. / M_PI * std::asin(lats[j]);
gridDefYvals(gaussGridID, lats.data());
}
......
......@@ -23,6 +23,7 @@
#include "process_int.h"
#include "cdo_cdi_wrapper.h"
#include "printinfo.h"
#include "gaussian_latitudes.h"
extern "C"
{
......@@ -75,37 +76,6 @@ rev_vals(double *vals, int n)
}
}
static bool
y_is_gauss(Varray<double> &gridyvals, int ysize)
{
bool lgauss = false;
int i;
if (ysize > 2)
{
Varray<double> yvals(ysize);
Varray<double> yw(ysize);
gaussianLatitudes(yvals.data(), yw.data(), (size_t) ysize);
for (i = 0; i < ysize; i++) yvals[i] = std::asin(yvals[i]) / M_PI * 180.0;
for (i = 0; i < ysize; i++)
if (std::fabs(yvals[i] - gridyvals[i]) > ((yvals[0] - yvals[1]) / 500)) break;
if (i == ysize) lgauss = true;
// check S->N
if (lgauss == false)
{
for (i = 0; i < ysize; i++)
if (std::fabs(yvals[i] - gridyvals[ysize - i - 1]) > ((yvals[0] - yvals[1]) / 500)) break;
if (i == ysize) lgauss = true;
}
}
return lgauss;
}
static int
define_grid(dsets_t *pfi)
{
......@@ -120,7 +90,7 @@ define_grid(dsets_t *pfi)
if (pfi->yrflg) rev_vals(yvals.data(), ny);
const bool lgauss = (pfi->linear[1] == 0) ? y_is_gauss(yvals, (size_t) ny) : false;
const bool lgauss = (pfi->linear[1] == 0) ? is_gaussian_latitudes((size_t) ny, yvals.data()) : false;
const int gridtype = lgauss ? GRID_GAUSSIAN : GRID_LONLAT;
......
......@@ -102,6 +102,8 @@ libcdo_la_SOURCES = after_dvtrans.cc \
functs.h \
eof_mode.cc \
eof_mode.h \
gaussian_latitudes.c \
gaussian_latitudes.h \
getMemorySize.c \
getRSS.c \
grid_area.cc \
......
......@@ -25,8 +25,7 @@
#include "cimdOmp.h"
#include "constants.h"
#include "array.h"
extern "C" void gaussianLatitudes(double *latitudes, double *weights, size_t nlat);
#include "gaussian_latitudes.h"
static void
jspleg1(double *pleg, double plat, long ktrunc, double *work)
......@@ -285,7 +284,7 @@ after_legini_full(long ntr, long nlat, double *restrict poli, double *restrict p
long dimsp = (ntr + 1) * (ntr + 2); // used in omp loop
Varray<double> gmu(nlat), gwt(nlat);
gaussianLatitudes(gmu.data(), gwt.data(), nlat);
gaussian_latitudes(nlat, gmu.data(), gwt.data());
#ifndef _OPENMP
Varray<double> pnm(dimsp), hnm(dimsp), ztemp1(waves << 1), ztemp2(waves << 1);
......@@ -343,7 +342,7 @@ after_legini(long ntr, long nlat, double *restrict poli, double *restrict pold,
#endif
Varray<double> gmu(nlat), gwt(nlat);
gaussianLatitudes(gmu.data(), gwt.data(), nlat);
gaussian_latitudes(nlat, gmu.data(), gwt.data());
for (long jgl = 0; jgl < nlat; ++jgl) gwt[jgl] *= 0.5;
for (long jgl = 0; jgl < nlat; ++jgl) coslat[jgl] = std::sqrt(1.0 - gmu[jgl] * gmu[jgl]);
......
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <float.h>
#include <math.h>
#ifndef M_SQRT2
#define M_SQRT2 1.41421356237309504880168872420969808
#endif
static
void cpledn(size_t kn, size_t kodd, double *pfn, double pdx, int kflag,
double *pw, double *pdxn, double *pxmod)
{
// 1.0 Newton iteration step
double zdlx = pdx;
double zdlldn = 0.0;
size_t ik = 1;
if (kflag == 0)
{
double zdlk = 0.5*pfn[0];
for (size_t jn = 2-kodd; jn <= kn; jn += 2)
{
// normalised ordinary Legendre polynomial == \overbar{p_n}^0
zdlk = zdlk + pfn[ik]*cos((double)(jn)*zdlx);
// normalised derivative == d/d\theta(\overbar{p_n}^0)
zdlldn = zdlldn - pfn[ik]*(double)(jn)*sin((double)(jn)*zdlx);
ik++;
}
// Newton method
double zdlmod = -(zdlk/zdlldn);
double zdlxn = zdlx + zdlmod;
*pdxn = zdlxn;
*pxmod = zdlmod;
}
// 2.0 Compute weights
if ( kflag == 1 )
{
for (size_t jn = 2-kodd; jn <= kn; jn += 2)
{
// normalised derivative
zdlldn = zdlldn - pfn[ik]*(double)(jn)*sin((double)(jn)*zdlx);
ik++;
}
*pw = (double)(2*kn+1)/(zdlldn*zdlldn);
}
return;
}
static
void gawl(double *pfn, double *pl, double *pw, size_t kn)
{
double pmod = 0.0;
double zw = 0.0;
double zdlxn = 0.0;
// 1.0 Initizialization
int iflag = 0;
int itemax = 20;
size_t iodd = (kn % 2);
double zdlx = *pl;
// 2.0 Newton iteration
for (int jter = 1; jter <= itemax+1; jter++)
{
cpledn(kn, iodd, pfn, zdlx, iflag, &zw, &zdlxn, &pmod);
zdlx = zdlxn;
if (iflag == 1) break;
if (fabs(pmod) <= DBL_EPSILON*1000.0) iflag = 1;
}
*pl = zdlxn;
*pw = zw;
return;
}
static
void gauaw(size_t kn, double *restrict pl, double *restrict pw)
{
/*
* 1.0 Initialize Fourier coefficients for ordinary Legendre polynomials
*
* Belousov, Swarztrauber, and ECHAM use zfn(0,0) = sqrt(2)
* IFS normalisation chosen to be 0.5*Integral(Pnm**2) = 1 (zfn(0,0) = 2.0)
*/
double *zfn = (double *) malloc((kn+1) * (kn+1) * sizeof(double));
double *zfnlat = (double *) malloc((kn/2+1+1)*sizeof(double));
zfn[0] = M_SQRT2;
for (size_t jn = 1; jn <= kn; jn++)
{
double zfnn = zfn[0];
for (size_t jgl = 1; jgl <= jn; jgl++)
{
zfnn *= sqrt(1.0-0.25/((double)(jgl*jgl)));
}
zfn[jn*(kn+1)+jn] = zfnn;
size_t iodd = jn % 2;
for (size_t jgl = 2; jgl <= jn-iodd; jgl += 2)
{
zfn[jn*(kn+1)+jn-jgl] = zfn[jn*(kn+1)+jn-jgl+2]
* ((double)((jgl-1)*(2*jn-jgl+2))) / ((double)(jgl*(2*jn-jgl+1)));
}
}
// 2.0 Gaussian latitudes and weights
size_t iodd = kn % 2;
size_t ik = iodd;
for (size_t jgl = iodd; jgl <= kn; jgl += 2)
{
zfnlat[ik] = zfn[kn*(kn+1)+jgl];
ik++;
}
// 2.1 Find first approximation of the roots of the Legendre polynomial of degree kn
size_t ins2 = kn/2+(kn % 2);
for (size_t jgl = 1; jgl <= ins2; jgl++)
{
double z = ((double)(4*jgl-1)) * M_PI / ((double)(4*kn+2));
pl[jgl-1] = z + 1.0 / (tan(z)*((double)(8*kn*kn)));
}
// 2.2 Computes roots and weights for transformed theta
for (size_t jgl = ins2; jgl >= 1 ; jgl--)
{
size_t jglm1 = jgl-1;
gawl(zfnlat, &(pl[jglm1]), &(pw[jglm1]), kn);
}
// convert to physical latitude
for (size_t jgl = 0; jgl < ins2; jgl++) pl[jgl] = cos(pl[jgl]);
for (size_t jgl = 1; jgl <= kn/2; jgl++)
{
size_t jglm1 = jgl-1;
size_t isym = kn-jgl;
pl[isym] = -pl[jglm1];
pw[isym] = pw[jglm1];
}
free(zfnlat);
free(zfn);
return;
}
void gaussian_latitudes(size_t nlats, double *latitudes, double *weights)
{
gauaw(nlats, latitudes, weights);
}
bool is_gaussian_latitudes(size_t nlats, const double *latitudes)
{
bool is_gauss_lats = false;
if (nlats > 2) // check if gaussian
{
size_t i;
double *yv = (double *) malloc(nlats*sizeof(double));
double *yw = (double *) malloc(nlats*sizeof(double));
gaussian_latitudes(nlats, yv, yw);
free(yw);
for (i = 0; i < nlats; i++)
yv[i] = asin(yv[i]) / M_PI * 180.0;
for (i = 0; i < nlats; i++)
if (fabs(yv[i] - latitudes[i]) > ((yv[0] - yv[1]) / 500.0)) break;
if (i == nlats) is_gauss_lats = true;
// check S->N
if (is_gauss_lats == false)
{
for (i = 0; i < nlats; i++)
if (fabs(yv[i] - latitudes[nlats-i-1]) > ((yv[0] - yv[1]) / 500.0)) break;
if (i == nlats) is_gauss_lats = true;
}
free(yv);
}
return is_gauss_lats;
}
#ifndef GAUSSIAN_LATITUDES_H_
#define GAUSSIAN_LATITUDES_H_
#include <stdio.h>
#ifdef __cplusplus
extern "C" {
#endif
void gaussian_latitudes(size_t nlats, double *latitudes, double *weights);
bool is_gaussian_latitudes(size_t nlats, const double *latitudes);
#ifdef __cplusplus
}
#endif
#endif /* GAUSSIAN_LATITUDES_H_ */
......@@ -21,6 +21,7 @@
#include <mpim_grid.h>
#include "griddes.h"
#include "util_string.h"
#include "gaussian_latitudes.h"
size_t genIcosphereCoords(int subdivisions, bool lbounds, std::vector<double> &xvals, std::vector<double> &yvals,
std::vector<double> &xbounds, std::vector<double> &ybounds);
......@@ -238,7 +239,7 @@ gaussianLatitudesInDegrees(std::vector<double> &lats, std::vector<double> &lat_b
// lat_bounds(nlat+1)
std::vector<double> latw(nlat), latw_cumsum(nlat);
gaussianLatitudes(lats.data(), latw.data(), nlat);
gaussian_latitudes(nlat, lats.data(), latw.data());
for (size_t j = 0; j < nlat; ++j) lats[j] = RAD2DEG * std::asin(lats[j]);
......
......@@ -23,6 +23,7 @@
#include <mpim_grid.h>
#include "cdo_output.h"
#include "compare.h"
#include "gaussian_latitudes.h"
static void
skip_nondigit_lines(FILE *gfp)
......@@ -156,21 +157,7 @@ grid_read_pingo(FILE *gfp)
return gridID;
}
bool lgauss = false;
if (nlat > 2) // check if gaussian
{
std::vector<double> yvals(grid.ysize), yw(grid.ysize);
gaussianLatitudes(yvals.data(), yw.data(), grid.ysize);
for (i = 0; i < (int) grid.ysize; i++) yvals[i] = std::asin(yvals[i]) * RAD2DEG;
for (i = 0; i < (int) grid.ysize; i++)
if (std::fabs(yvals[i] - grid.yvals[i]) > ((yvals[0] - yvals[1]) / 500)) break;
if (i == (int) grid.ysize) lgauss = true;
}
grid.type = lgauss ? GRID_GAUSSIAN : GRID_LONLAT;
grid.type = is_gaussian_latitudes(nlat, grid.yvals.data()) ? GRID_GAUSSIAN : GRID_LONLAT;
}
if (grid.type != CDI_UNDEFID) gridID = gridDefine(grid);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment