Skip to content

GitLab

Commit 103f655f authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Moved isGaussLat() to grid.c.

parent 38ef22b4
Hide whitespace changes
Inline Side-by-side
src/gaussgrid.c
View file @ 103f655f
......@@ -3,6 +3,7 @@
#endif
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
......@@ -10,7 +11,6 @@
#define M_SQRT2 1.41421356237309504880168872420969808
#endif
#include "dmemory.h"
#include "cdi_int.h"
......@@ -18,40 +18,38 @@ 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 */
// 1.0 Newton iteration step
double zdlx = pdx;
double zdlk = 0.0;
if ( kodd == 0 ) zdlk = 0.5*pfn[0];
double zdlxn = 0.0;
double zdlldn = 0.0;
size_t ik = 1;
if ( kflag == 0 )
if (kflag == 0)
{
for ( size_t jn = 2-kodd; jn <= kn; jn += 2 )
double zdlk = 0.5*pfn[0];
for (size_t jn = 2-kodd; jn <= kn; jn += 2)
{
/* normalised ordinary Legendre polynomial == \overbar{p_n}^0 */
// 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) */
// normalised derivative == d/d\theta(\overbar{p_n}^0)
zdlldn = zdlldn - pfn[ik]*(double)(jn)*sin((double)(jn)*zdlx);
ik++;
}
/* Newton method */
// Newton method
double zdlmod = -(zdlk/zdlldn);
zdlxn = zdlx + zdlmod;
double zdlxn = zdlx + zdlmod;
*pdxn = zdlxn;
*pxmod = zdlmod;
}
/* 2.0 Compute weights */
// 2.0 Compute weights
if ( kflag == 1 )
{
for ( size_t jn = 2-kodd; jn <= kn; jn += 2 )
for (size_t jn = 2-kodd; jn <= kn; jn += 2)
{
/* normalised derivative */
// normalised derivative
zdlldn = zdlldn - pfn[ik]*(double)(jn)*sin((double)(jn)*zdlx);
ik++;
}
......@@ -64,11 +62,11 @@ void cpledn(size_t kn, size_t kodd, double *pfn, double pdx, int kflag,
static
void gawl(double *pfn, double *pl, double *pw, size_t kn)
{
double pmod = 0;
double zw = 0;
double zdlxn = 0;
double pmod = 0.0;
double zw = 0.0;
double zdlxn = 0.0;
/* 1.0 Initizialization */
// 1.0 Initizialization
int iflag = 0;
int itemax = 20;
......@@ -77,9 +75,9 @@ void gawl(double *pfn, double *pl, double *pw, size_t kn)
double zdlx = *pl;
/* 2.0 Newton iteration */
// 2.0 Newton iteration
for ( int jter = 1; jter <= itemax+1; jter++ )
for (int jter = 1; jter <= itemax+1; jter++)
{
cpledn(kn, iodd, pfn, zdlx, iflag, &zw, &zdlxn, &pmod);
zdlx = zdlxn;
......@@ -102,11 +100,11 @@ void gauaw(size_t kn, double *restrict pl, double *restrict pw)
* 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));
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++ )
for (size_t jn = 1; jn <= kn; jn++)
{
double zfnn = zfn[0];
for (size_t jgl = 1; jgl <= jn; jgl++)
......@@ -117,7 +115,7 @@ void gauaw(size_t kn, double *restrict pl, double *restrict pw)
zfn[jn*(kn+1)+jn] = zfnn;
size_t iodd = jn % 2;
for ( size_t jgl = 2; jgl <= jn-iodd; jgl += 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)));
......@@ -125,43 +123,40 @@ void gauaw(size_t kn, double *restrict pl, double *restrict pw)
}
/* 2.0 Gaussian latitudes and weights */
// 2.0 Gaussian latitudes and weights
size_t iodd = kn % 2;
size_t ik = iodd;
for ( size_t jgl = iodd; jgl <= kn; jgl += 2 )
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 */
// 2.1 Find first approximation of the roots of the Legendre polynomial of degree kn
size_t ins2 = kn/2+(kn % 2);
double z;
for ( size_t jgl = 1; jgl <= ins2; jgl++ )
for (size_t jgl = 1; jgl <= ins2; jgl++)
{
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 */
// 2.2 Computes roots and weights for transformed theta
for ( size_t jgl = ins2; jgl >= 1 ; jgl-- )
for (size_t jgl = ins2; jgl >= 1 ; jgl--)
{
size_t jglm1 = jgl-1;
gawl(zfnlat, &(pl[jglm1]), &(pw[jglm1]), kn);
}
/* convert to physical latitude */
// convert to physical latitude
for ( size_t jgl = 0; jgl < ins2; jgl++ )
{
pl[jgl] = cos(pl[jgl]);
}
for (size_t jgl = 0; jgl < ins2; jgl++) pl[jgl] = cos(pl[jgl]);
for ( size_t jgl = 1; jgl <= kn/2; jgl++ )
for (size_t jgl = 1; jgl <= kn/2; jgl++)
{
size_t jglm1 = jgl-1;
size_t isym = kn-jgl;
......@@ -169,8 +164,8 @@ void gauaw(size_t kn, double *restrict pl, double *restrict pw)
pw[isym] = pw[jglm1];
}
Free(zfnlat);
Free(zfn);
free(zfnlat);
free(zfn);
return;
}
......@@ -178,71 +173,9 @@ void gauaw(size_t kn, double *restrict pl, double *restrict pw)
void gaussianLatitudes(double *latitudes, double *weights, size_t nlat)
{
//gauaw_old(latitudes, weights, nlat);
gauaw(nlat, latitudes, weights);
}
bool isGaussGrid(size_t ysize, double yinc, const double *yvals)
{
bool lgauss = false;
if ( IS_EQUAL(yinc, 0) && ysize > 2 ) // check if gaussian
{
size_t i;
double *yv = (double *) Malloc(ysize*sizeof(double));
double *yw = (double *) Malloc(ysize*sizeof(double));
gaussianLatitudes(yv, yw, ysize);
Free(yw);
for ( i = 0; i < ysize; i++ )
yv[i] = asin(yv[i])/M_PI*180.0;
for ( i = 0; i < ysize; i++ )
if ( fabs(yv[i] - yvals[i]) > ((yv[0] - yv[1])/500) ) break;
if ( i == ysize ) lgauss = true;
// check S->N
if ( lgauss == false )
{
for ( i = 0; i < ysize; i++ )
if ( fabs(yv[i] - yvals[ysize-i-1]) > ((yv[0] - yv[1])/500) ) break;
if ( i == ysize ) lgauss = true;
}
Free(yv);
}
return lgauss;
}
/*
#define NGL 48
int main (int rgc, char *argv[])
{
int ngl = NGL;
double plo[NGL], pwo[NGL];
double pl[NGL], pw[NGL];
int i;
gauaw(ngl, pl, pw);
for (i = 0; i < ngl; i++)
{
pl[i] = asin(pl[i])/M_PI*180.0;
plo[i] = asin(plo[i])/M_PI*180.0;
}
for (i = 0; i < ngl; i++)
{
fprintf(stderr, "%4d%25.18f%25.18f%25.18f%25.18f\n", i+1, pl[i], pw[i], pl[i]-plo[i], pw[i]-pwo[i]);
}
return 0;
}
*/
/*
* Local Variables:
* c-file-style: "Java"
......
src/grid.c
View file @ 103f655f
......@@ -4051,6 +4051,42 @@ void gridInqUUID(int gridID, unsigned char uuid[CDI_UUID_SIZE])
}
bool isGaussLat(size_t ysize, double yinc, const double *yvals)
{
bool is_gauss_grid = false;
if (IS_EQUAL(yinc, 0) && ysize > 2) // check if gaussian
{
size_t i;
double *yv = (double *) malloc(ysize*sizeof(double));
double *yw = (double *) malloc(ysize*sizeof(double));
gaussianLatitudes(yv, yw, ysize);
free(yw);
for (i = 0; i < ysize; i++)
yv[i] = asin(yv[i])/M_PI*180.0;
for (i = 0; i < ysize; i++)
if (fabs(yv[i] - yvals[i]) > ((yv[0] - yv[1]) / 500.0)) break;
if (i == ysize) is_gauss_grid = true;
// check S->N
if (is_gauss_grid == false)
{
for (i = 0; i < ysize; i++)
if (fabs(yv[i] - yvals[ysize-i-1]) > ((yv[0] - yv[1]) / 500.0)) break;
if (i == ysize) is_gauss_grid = true;
}
free(yv);
}
return is_gauss_grid;
}
void cdiGridGetIndexList(unsigned ngrids, int * gridIndexList)
{
reshGetResHListOfType(ngrids, gridIndexList, &gridOps);
......
src/grid.h
View file @ 103f655f
......@@ -154,7 +154,7 @@ int gridVerifyGribParamLCC(double missval, double *lon_0, double *lat_0, double
int gridVerifyGribParamSTERE(double missval, double *lon_0, double *lat_ts, double *lat_0,
double *a, double *xval_0, double *yval_0, double *x_0, double *y_0);
bool isGaussGrid(size_t ysize, double yinc, const double *yvals);
bool isGaussLat(size_t ysize, double yinc, const double *yvals);
#endif
/*
......
src/stream_cdf_i.c
View file @ 103f655f
......@@ -2190,7 +2190,7 @@ void cdf_check_gridtype(int *gridtype, bool islon, bool islat, size_t xsize, siz
break;
}
}
if ( ysize < 10000 && isGaussGrid(ysize, yinc, grid->y.vals) )
if ( ysize < 10000 && isGaussLat(ysize, yinc, grid->y.vals) )
{
*gridtype = GRID_GAUSSIAN;
grid->np = (int)(ysize/2);
......
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