Skip to content
GitLab
Commit 0241edb6 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

cgribexlib update

parent 4df0b395
Hide whitespace changes
Inline Side-by-side
src/cgribex.h
View file @ 0241edb6
......@@ -243,7 +243,8 @@ int gribVersion(unsigned char *buffer, size_t buffersize);
int grib_info_for_grads(off_t recpos, long recsize, unsigned char *gribbuffer, int *intnum, float *fltnum, off_t *bignum);
double calculate_pfactor(const double* spectralField, long fieldTruncation, long subsetTruncation);
double calculate_pfactor_float(const float* spectralField, long fieldTruncation, long subsetTruncation);
double calculate_pfactor_double(const double* spectralField, long fieldTruncation, long subsetTruncation);
#endif /* _CGRIBEX_H */
src/cgribexlib.c
View file @ 0241edb6
/* Automatically generated by m214003 at 2015-10-30, do not edit */
/* Automatically generated by m214003 at 2015-11-11, do not edit */
/* CGRIBEXLIB_VERSION="1.7.4" */
......@@ -7027,113 +7027,6 @@ int gribUnzip(unsigned char *dbuf, long dbufsize, unsigned char *sbuf, long sbu
#include <math.h>
/* calculate_pfactor: source code from grib_api-1.8.0 */
double calculate_pfactor(const double* spectralField, long fieldTruncation, long subsetTruncation)
{
/*long n_vals = ((fieldTruncation+1)*(fieldTruncation+2));*/
long loop, index, m, n = 0;
double pFactor, zeps = 1.0e-15;
long ismin = (subsetTruncation+1), ismax = (fieldTruncation+1);
double* weights, range, * norms;
double weightedSumOverX = 0.0, weightedSumOverY = 0.0, sumOfWeights = 0.0, x, y;
double numerator = 0.0, denominator = 0.0, slope;
/*
// Setup the weights
*/
range = (double) (ismax - ismin +1);
weights = (double*) Malloc((ismax+1)*sizeof(double));
for( loop = ismin; loop <= ismax; loop++ )
weights[loop] = range / (double) (loop-ismin+1);
/*
// Compute norms
// Handle values 2 at a time (real and imaginary parts).
*/
norms = (double*) Malloc((ismax+1)*sizeof(double));
for( loop = 0; loop < ismax+1; loop++ ) norms[loop] = 0.0;
/*
// Form norms for the rows which contain part of the unscaled subset.
*/
index = -2;
for( m = 0; m < subsetTruncation; m++ )
for( n = m; n <= fieldTruncation; n++ ) {
index += 2;
if( n >= subsetTruncation ) {
double tval = spectralField[index];
tval=tval<0?-tval:tval;
norms[n] = norms[n] > tval ? norms[n] : tval;
tval = spectralField[index+1];
tval=tval<0?-tval:tval;
norms[n] = norms[n] > tval ? norms[n] : tval;
}
}
/*
// Form norms for the rows which do not contain part of the unscaled subset.
*/
for( m = subsetTruncation; m <= fieldTruncation; m++ )
for( n = m; n <= fieldTruncation; n++ ) {
double tval = spectralField[index];
index += 2;
tval=tval<0?-tval:tval;
norms[n] = norms[n] > tval ? norms[n] : tval;
tval = spectralField[index+1];
tval=tval<0?-tval:tval;
norms[n] = norms[n] > tval ? norms[n] : tval;
}
/*
// Ensure the norms have a value which is not too small in case of
// problems with math functions (e.g. LOG).
*/
for( loop = ismin; loop <= ismax; loop++ ) {
norms[n] = norms[n] > zeps ? norms[n] : zeps;
if( IS_EQUAL(norms[n], zeps) ) weights[n] = 100.0 * zeps;
}
/*
// Do linear fit to find the slope
*/
for( loop = ismin; loop <= ismax; loop++ ) {
x = log( (double) (loop*(loop+1)) );
y = log( norms[loop] );
weightedSumOverX = weightedSumOverX + x * weights[loop];
weightedSumOverY = weightedSumOverY + y * weights[loop];
sumOfWeights = sumOfWeights + weights[loop];
}
weightedSumOverX = weightedSumOverX / sumOfWeights;
weightedSumOverY = weightedSumOverY / sumOfWeights;
/*
// Perform a least square fit for the equation
*/
for( loop = ismin; loop <= ismax; loop++ ) {
x = log( (double)(loop*(loop+1)) );
y = log( norms[loop] );
numerator =
numerator + weights[loop] * (y-weightedSumOverY) * (x-weightedSumOverX);
denominator =
denominator + weights[loop] * ((x-weightedSumOverX) * (x-weightedSumOverX));
}
slope = numerator / denominator;
Free(weights);
Free(norms);
pFactor = -slope;
if( pFactor < -9999.9 ) pFactor = -9999.9;
if( pFactor > 9999.9 ) pFactor = 9999.9;
return pFactor;
}
static
int rowina2(double *p, int ko, int ki, double *pw,
......@@ -7474,6 +7367,114 @@ L900:
#define T double
#ifdef T
/* calculate_pfactor: source code from grib_api-1.8.0 */
double TEMPLATE(calculate_pfactor,T)(const T *spectralField, long fieldTruncation, long subsetTruncation)
{
/*long n_vals = ((fieldTruncation+1)*(fieldTruncation+2));*/
long loop, index, m, n = 0;
double pFactor, zeps = 1.0e-15;
long ismin = (subsetTruncation+1), ismax = (fieldTruncation+1);
double* weights, range, * norms;
double weightedSumOverX = 0.0, weightedSumOverY = 0.0, sumOfWeights = 0.0, x, y;
double numerator = 0.0, denominator = 0.0, slope;
/*
// Setup the weights
*/
range = (double) (ismax - ismin +1);
weights = (double*) Malloc((ismax+1)*sizeof(double));
for( loop = ismin; loop <= ismax; loop++ )
weights[loop] = range / (double) (loop-ismin+1);
/*
// Compute norms
// Handle values 2 at a time (real and imaginary parts).
*/
norms = (double*) Malloc((ismax+1)*sizeof(double));
for( loop = 0; loop < ismax+1; loop++ ) norms[loop] = 0.0;
/*
// Form norms for the rows which contain part of the unscaled subset.
*/
index = -2;
for( m = 0; m < subsetTruncation; m++ )
for( n = m; n <= fieldTruncation; n++ ) {
index += 2;
if( n >= subsetTruncation ) {
double tval = spectralField[index];
tval=tval<0?-tval:tval;
norms[n] = norms[n] > tval ? norms[n] : tval;
tval = spectralField[index+1];
tval=tval<0?-tval:tval;
norms[n] = norms[n] > tval ? norms[n] : tval;
}
}
/*
// Form norms for the rows which do not contain part of the unscaled subset.
*/
for( m = subsetTruncation; m <= fieldTruncation; m++ )
for( n = m; n <= fieldTruncation; n++ ) {
double tval = spectralField[index];
index += 2;
tval=tval<0?-tval:tval;
norms[n] = norms[n] > tval ? norms[n] : tval;
tval = spectralField[index+1];
tval=tval<0?-tval:tval;
norms[n] = norms[n] > tval ? norms[n] : tval;
}
/*
// Ensure the norms have a value which is not too small in case of
// problems with math functions (e.g. LOG).
*/
for( loop = ismin; loop <= ismax; loop++ ) {
norms[n] = norms[n] > zeps ? norms[n] : zeps;
if( IS_EQUAL(norms[n], zeps) ) weights[n] = 100.0 * zeps;
}
/*
// Do linear fit to find the slope
*/
for( loop = ismin; loop <= ismax; loop++ ) {
x = log( (double) (loop*(loop+1)) );
y = log( norms[loop] );
weightedSumOverX = weightedSumOverX + x * weights[loop];
weightedSumOverY = weightedSumOverY + y * weights[loop];
sumOfWeights = sumOfWeights + weights[loop];
}
weightedSumOverX = weightedSumOverX / sumOfWeights;
weightedSumOverY = weightedSumOverY / sumOfWeights;
/*
// Perform a least square fit for the equation
*/
for( loop = ismin; loop <= ismax; loop++ ) {
x = log( (double)(loop*(loop+1)) );
y = log( norms[loop] );
numerator =
numerator + weights[loop] * (y-weightedSumOverY) * (x-weightedSumOverX);
denominator =
denominator + weights[loop] * ((x-weightedSumOverX) * (x-weightedSumOverX));
}
slope = numerator / denominator;
Free(weights);
Free(norms);
pFactor = -slope;
if( pFactor < -9999.9 ) pFactor = -9999.9;
if( pFactor > 9999.9 ) pFactor = 9999.9;
return pFactor;
}
void TEMPLATE(scale_complex,T)(T *fpdata, int pcStart, int pcScale, int trunc, int inv)
{
double power;
......@@ -8137,6 +8138,114 @@ L900:
#define T float
#ifdef T
/* calculate_pfactor: source code from grib_api-1.8.0 */
double TEMPLATE(calculate_pfactor,T)(const T *spectralField, long fieldTruncation, long subsetTruncation)
{
/*long n_vals = ((fieldTruncation+1)*(fieldTruncation+2));*/
long loop, index, m, n = 0;
double pFactor, zeps = 1.0e-15;
long ismin = (subsetTruncation+1), ismax = (fieldTruncation+1);
double* weights, range, * norms;
double weightedSumOverX = 0.0, weightedSumOverY = 0.0, sumOfWeights = 0.0, x, y;
double numerator = 0.0, denominator = 0.0, slope;
/*
// Setup the weights
*/
range = (double) (ismax - ismin +1);
weights = (double*) Malloc((ismax+1)*sizeof(double));
for( loop = ismin; loop <= ismax; loop++ )
weights[loop] = range / (double) (loop-ismin+1);
/*
// Compute norms
// Handle values 2 at a time (real and imaginary parts).
*/
norms = (double*) Malloc((ismax+1)*sizeof(double));
for( loop = 0; loop < ismax+1; loop++ ) norms[loop] = 0.0;
/*
// Form norms for the rows which contain part of the unscaled subset.
*/
index = -2;
for( m = 0; m < subsetTruncation; m++ )
for( n = m; n <= fieldTruncation; n++ ) {
index += 2;
if( n >= subsetTruncation ) {
double tval = spectralField[index];
tval=tval<0?-tval:tval;
norms[n] = norms[n] > tval ? norms[n] : tval;
tval = spectralField[index+1];
tval=tval<0?-tval:tval;
norms[n] = norms[n] > tval ? norms[n] : tval;
}
}
/*
// Form norms for the rows which do not contain part of the unscaled subset.
*/
for( m = subsetTruncation; m <= fieldTruncation; m++ )
for( n = m; n <= fieldTruncation; n++ ) {
double tval = spectralField[index];
index += 2;
tval=tval<0?-tval:tval;
norms[n] = norms[n] > tval ? norms[n] : tval;
tval = spectralField[index+1];
tval=tval<0?-tval:tval;
norms[n] = norms[n] > tval ? norms[n] : tval;
}
/*
// Ensure the norms have a value which is not too small in case of
// problems with math functions (e.g. LOG).
*/
for( loop = ismin; loop <= ismax; loop++ ) {
norms[n] = norms[n] > zeps ? norms[n] : zeps;
if( IS_EQUAL(norms[n], zeps) ) weights[n] = 100.0 * zeps;
}
/*
// Do linear fit to find the slope
*/
for( loop = ismin; loop <= ismax; loop++ ) {
x = log( (double) (loop*(loop+1)) );
y = log( norms[loop] );
weightedSumOverX = weightedSumOverX + x * weights[loop];
weightedSumOverY = weightedSumOverY + y * weights[loop];
sumOfWeights = sumOfWeights + weights[loop];
}
weightedSumOverX = weightedSumOverX / sumOfWeights;
weightedSumOverY = weightedSumOverY / sumOfWeights;
/*
// Perform a least square fit for the equation
*/
for( loop = ismin; loop <= ismax; loop++ ) {
x = log( (double)(loop*(loop+1)) );
y = log( norms[loop] );
numerator =
numerator + weights[loop] * (y-weightedSumOverY) * (x-weightedSumOverX);
denominator =
denominator + weights[loop] * ((x-weightedSumOverX) * (x-weightedSumOverX));
}
slope = numerator / denominator;
Free(weights);
Free(norms);
pFactor = -slope;
if( pFactor < -9999.9 ) pFactor = -9999.9;
if( pFactor > 9999.9 ) pFactor = 9999.9;
return pFactor;
}
void TEMPLATE(scale_complex,T)(T *fpdata, int pcStart, int pcScale, int trunc, int inv)
{
double power;
......@@ -14086,7 +14195,7 @@ void encode_dummy(void)
(void) encode_array_unrolled_double(0, 0, 0, NULL, NULL, 0, 0, NULL);
(void) encode_array_unrolled_float(0, 0, 0, NULL, NULL, 0, 0, NULL);
}
static const char grb_libvers[] = "1.7.4" " of ""Oct 30 2015"" ""16:20:41";
static const char grb_libvers[] = "1.7.4" " of ""Nov 11 2015"" ""14:13:54";
const char *
cgribexLibraryVersion(void)
{
......
src/stream_cgribex.c
View file @ 0241edb6
......@@ -2205,7 +2205,10 @@ size_t cgribexEncode(int memtype, int varID, int levelID, int vlistID, int gridI
if ( isec4[2] == 128 && isec4[3] == 64 )
{
isec4[16] = (int) (1000*calculate_pfactor(data, ISEC2_PentaJ, isec4[17]));
if ( memtype == MEMTYPE_FLOAT )
isec4[16] = (int) (1000*calculate_pfactor_float((const float*) data, ISEC2_PentaJ, isec4[17]));
else
isec4[16] = (int) (1000*calculate_pfactor_double((const double*) data, ISEC2_PentaJ, isec4[17]));
if ( isec4[16] < -10000 ) isec4[16] = -10000;
if ( isec4[16] > 10000 ) isec4[16] = 10000;
}
......
Supports Markdown
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

Imprint | Privacy Policy