core
/
eccodes
mirror of https://github.com/ecmwf/eccodes.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

ECC-1818: GRIB Geoiterator issues for Lambert azimuthal equal area

This commit is contained in:
Shahram Najm 2024-05-02 10:25:11 +00:00
parent e06656c673
commit c04bc4fa03
1 changed files with 11 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
src/grib_iterator_class_lambert_azimuthal_equal_area.cc
View File

@ -209,7 +209,11 @@ static int init_oblate(grib_handle* h,
sinphi_ = sin(standardParallelInRadians); /* (P->phi0); */
Q__sinb1 = pj_qsfn(sinphi_, e, one_es) / Q__qp;
Q__cosb1 = sqrt(1.0 - Q__sinb1 * Q__sinb1);
if (Q__cosb1 == 0) {
Q__dd = 1.0;
} else {
Q__dd = cos(standardParallelInRadians) / (sqrt(1. - es * sinphi_ * sinphi_) * Q__rq * Q__cosb1);
}
Q__ymf = (Q__xmf = Q__rq) / Q__dd;
Q__xmf *= Q__dd;
@ -253,7 +257,12 @@ static int init_oblate(grib_handle* h,
xy_y *= Q__dd;
rho = hypot(xy_x, xy_y);
Assert(rho >= EPS10); /* TODO(masn): check */
sCe = 2. * asin(.5 * rho / Q__rq);
const double asin_arg = (0.5 * rho / Q__rq);
if (asin_arg < -1.0 || asin_arg > 1.0) {
grib_context_log(h->context, GRIB_LOG_ERROR, "Invalid value: arcsin argument=%g", asin_arg);
return GRIB_GEOCALCULUS_PROBLEM;
}
sCe = 2. * asin(asin_arg);
cCe = cos(sCe);
sCe = sin(sCe);
xy_x *= sCe;