Skip to content
GitLab
Commit 1190cca0 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

cgribexlib update.

parent 9d386a23
Hide whitespace changes
Inline Side-by-side
src/cgribexlib.c
View file @ 1190cca0
/* Automatically generated by m214003 at 2019-09-25, do not edit */
/* Automatically generated by m214003 at 2019-10-29, do not edit */
/* CGRIBEXLIB_VERSION="1.9.4" */
......@@ -6644,31 +6644,26 @@ double TEMPLATE(calculate_pfactor,T)(const T *spectralField, long fieldTruncatio
{
/*long n_vals = ((fieldTruncation+1)*(fieldTruncation+2));*/
long loop, index, m, n = 0;
double pFactor, zeps = 1.0e-15;
double 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;
double weightedSumOverX = 0.0, weightedSumOverY = 0.0, sumOfWeights = 0.0;
double numerator = 0.0, denominator = 0.0;
/*
// Setup the weights
*/
range = (double) (ismax - ismin +1);
double range = (double) (ismax - ismin +1);
weights = (double*) Malloc(((size_t)ismax+1)*sizeof(double));
double *weights = (double*) Malloc(((size_t)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(((size_t)ismax+1)*sizeof(double));
double *norms = (double*) Malloc(((size_t)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++ )
......@@ -6683,9 +6678,8 @@ double TEMPLATE(calculate_pfactor,T)(const T *spectralField, long fieldTruncatio
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++ ) {
......@@ -6698,49 +6692,40 @@ double TEMPLATE(calculate_pfactor,T)(const T *spectralField, long fieldTruncatio
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).
*/
// 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];
double x = log( (double) (loop*(loop+1)) );
double y = log( norms[loop] );
weightedSumOverX += x * weights[loop];
weightedSumOverY += y * weights[loop];
sumOfWeights = sumOfWeights + weights[loop];
}
weightedSumOverX = weightedSumOverX / sumOfWeights;
weightedSumOverY = weightedSumOverY / sumOfWeights;
weightedSumOverX /= sumOfWeights;
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));
double x = log( (double)(loop*(loop+1)) );
double y = log( norms[loop] );
numerator += weights[loop] * (y-weightedSumOverY) * (x-weightedSumOverX);
denominator += weights[loop] * ((x-weightedSumOverX) * (x-weightedSumOverX));
}
slope = numerator / denominator;
double slope = numerator / denominator;
Free(weights);
Free(norms);
pFactor = -slope;
double pFactor = -slope;
if( pFactor < -9999.9 ) pFactor = -9999.9;
if( pFactor > 9999.9 ) pFactor = 9999.9;
......@@ -6863,14 +6848,6 @@ void TEMPLATE(gather_complex,T)(T *fpdata, size_t pcStart, size_t trunc, size_t
static void TEMPLATE(scm0,T)(T *pdl, T *pdr, T *pfl, T *pfr, int klg)
{
/* System generated locals */
double r_1;
/* Local variables */
int jl;
double zfac, zeps, zbeta;
double zalpha;
/* **** SCM0 - Apply SCM0 limiter to derivative estimates. */
/* output: */
/* pdl = the limited derivative at the left edge of the interval */
......@@ -6884,15 +6861,16 @@ static void TEMPLATE(scm0,T)(T *pdl, T *pdr, T *pfl, T *pfr, int klg)
/* define constants */
zeps = 1.0e-12;
zfac = (1.0 - zeps) * 3.0;
double zeps = 1.0e-12;
double zfac = (1.0 - zeps) * 3.0;
for ( jl = 0; jl < klg; ++jl )
for ( int jl = 0; jl < klg; ++jl )
{
double r_1;
if ( (r_1 = pfr[jl] - pfl[jl], fabs(r_1)) > zeps )
{
zalpha = pdl[jl] / (pfr[jl] - pfl[jl]);
zbeta = pdr[jl] / (pfr[jl] - pfl[jl]);
double zalpha = pdl[jl] / (pfr[jl] - pfl[jl]);
double zbeta = pdr[jl] / (pfr[jl] - pfl[jl]);
if ( zalpha <= 0.0 ) pdl[jl] = 0.0;
if ( zbeta <= 0.0 ) pdr[jl] = 0.0;
if ( zalpha > zfac ) pdl[jl] = (T)(zfac * (pfr[jl] - pfl[jl]));
......@@ -7047,7 +7025,7 @@ C -----------------------------------------------------------------
/* Get the current array position(minus 1) from the weight - */
/* note the implicit truncation. */
ip = (int) zwt;
/* Adjust the weight to range (0.0 to 1.0) */
zwt -= ip;
......@@ -7394,31 +7372,26 @@ double TEMPLATE(calculate_pfactor,T)(const T *spectralField, long fieldTruncatio
{
/*long n_vals = ((fieldTruncation+1)*(fieldTruncation+2));*/
long loop, index, m, n = 0;
double pFactor, zeps = 1.0e-15;
double 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;
double weightedSumOverX = 0.0, weightedSumOverY = 0.0, sumOfWeights = 0.0;
double numerator = 0.0, denominator = 0.0;
/*
// Setup the weights
*/
range = (double) (ismax - ismin +1);
double range = (double) (ismax - ismin +1);
weights = (double*) Malloc(((size_t)ismax+1)*sizeof(double));
double *weights = (double*) Malloc(((size_t)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(((size_t)ismax+1)*sizeof(double));
double *norms = (double*) Malloc(((size_t)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++ )
......@@ -7433,9 +7406,8 @@ double TEMPLATE(calculate_pfactor,T)(const T *spectralField, long fieldTruncatio
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++ ) {
......@@ -7448,49 +7420,40 @@ double TEMPLATE(calculate_pfactor,T)(const T *spectralField, long fieldTruncatio
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).
*/
// 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];
double x = log( (double) (loop*(loop+1)) );
double y = log( norms[loop] );
weightedSumOverX += x * weights[loop];
weightedSumOverY += y * weights[loop];
sumOfWeights = sumOfWeights + weights[loop];
}
weightedSumOverX = weightedSumOverX / sumOfWeights;
weightedSumOverY = weightedSumOverY / sumOfWeights;
weightedSumOverX /= sumOfWeights;
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));
double x = log( (double)(loop*(loop+1)) );
double y = log( norms[loop] );
numerator += weights[loop] * (y-weightedSumOverY) * (x-weightedSumOverX);
denominator += weights[loop] * ((x-weightedSumOverX) * (x-weightedSumOverX));
}
slope = numerator / denominator;
double slope = numerator / denominator;
Free(weights);
Free(norms);
pFactor = -slope;
double pFactor = -slope;
if( pFactor < -9999.9 ) pFactor = -9999.9;
if( pFactor > 9999.9 ) pFactor = 9999.9;
......@@ -7613,14 +7576,6 @@ void TEMPLATE(gather_complex,T)(T *fpdata, size_t pcStart, size_t trunc, size_t
static void TEMPLATE(scm0,T)(T *pdl, T *pdr, T *pfl, T *pfr, int klg)
{
/* System generated locals */
double r_1;
/* Local variables */
int jl;
double zfac, zeps, zbeta;
double zalpha;
/* **** SCM0 - Apply SCM0 limiter to derivative estimates. */
/* output: */
/* pdl = the limited derivative at the left edge of the interval */
......@@ -7634,15 +7589,16 @@ static void TEMPLATE(scm0,T)(T *pdl, T *pdr, T *pfl, T *pfr, int klg)
/* define constants */
zeps = 1.0e-12;
zfac = (1.0 - zeps) * 3.0;
double zeps = 1.0e-12;
double zfac = (1.0 - zeps) * 3.0;
for ( jl = 0; jl < klg; ++jl )
for ( int jl = 0; jl < klg; ++jl )
{
double r_1;
if ( (r_1 = pfr[jl] - pfl[jl], fabs(r_1)) > zeps )
{
zalpha = pdl[jl] / (pfr[jl] - pfl[jl]);
zbeta = pdr[jl] / (pfr[jl] - pfl[jl]);
double zalpha = pdl[jl] / (pfr[jl] - pfl[jl]);
double zbeta = pdr[jl] / (pfr[jl] - pfl[jl]);
if ( zalpha <= 0.0 ) pdl[jl] = 0.0;
if ( zbeta <= 0.0 ) pdr[jl] = 0.0;
if ( zalpha > zfac ) pdl[jl] = (T)(zfac * (pfr[jl] - pfl[jl]));
......@@ -7797,7 +7753,7 @@ C -----------------------------------------------------------------
/* Get the current array position(minus 1) from the weight - */
/* note the implicit truncation. */
ip = (int) zwt;
/* Adjust the weight to range (0.0 to 1.0) */
zwt -= ip;
......
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