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

Moved isGaussianLatitudes to gaussian_latitudes.c

parent d10a502f
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <float.h>
#include <math.h>
......@@ -163,7 +164,43 @@ void gauaw(size_t kn, double *restrict pl, double *restrict pw)
}
void gaussianLatitudes(double *latitudes, double *weights, size_t nlat)
void gaussianLatitudes(size_t nlats, double *latitudes, double *weights)
{
gauaw(nlat, latitudes, weights);
gauaw(nlats, latitudes, weights);
}
bool isGaussianLatitudes(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));
gaussianLatitudes(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;
}
......@@ -359,9 +359,9 @@ void gridGenXvals(int xsize, double xfirst, double xlast, double xinc, double *r
static
void calc_gaussgrid(double *restrict yvals, size_t ysize, double yfirst, double ylast)
{
double *restrict yw = (double *) Malloc(ysize * sizeof(double));
gaussianLatitudes(yvals, yw, ysize);
Free(yw);
double *restrict yw = (double *) malloc(ysize * sizeof(double));
gaussianLatitudes(ysize, yvals, yw);
free(yw);
for (size_t i = 0; i < ysize; i++ )
yvals[i] = asin(yvals[i])/M_PI*180.0;
......@@ -4051,42 +4051,6 @@ 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);
......
......@@ -154,8 +154,8 @@ 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 isGaussLat(size_t ysize, double yinc, const double *yvals);
void gaussianLatitudes(double *latitudes, double *weights, size_t nlat);
bool isGaussianLatitudes(size_t nlats, const double *latitudes);
void gaussianLatitudes(size_t nlats, double *latitudes, double *weights);
#endif
/*
......
......@@ -2190,7 +2190,7 @@ void cdf_check_gridtype(int *gridtype, bool islon, bool islat, size_t xsize, siz
break;
}
}
if ( ysize < 10000 && isGaussLat(ysize, yinc, grid->y.vals) )
if ( ysize < 10000 && IS_EQUAL(yinc, 0.0) && isGaussianLatitudes(ysize, 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