Commit 138092d9 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

added operator std1 and var1 for module Timstat

parent 09778229
......@@ -3,6 +3,19 @@
* using CDI library version 1.6.0
* Version 1.6.0 released
2013-01-23 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* New operator: timvar1 - Time variance [Divisor is (n-1)]
* New operator: timstd1 - Time standard deviation [Divisor is (n-1)]
* New operator: hourvar1 - Hourly variance [Divisor is (n-1)]
* New operator: hourstd1 - Hourly standard deviation [Divisor is (n-1)]
* New operator: dayvar1 - Daily variance [Divisor is (n-1)]
* New operator: daystd1 - Daily standard deviation [Divisor is (n-1)]
* New operator: monvar1 - Monthly variance [Divisor is (n-1)]
* New operator: monstd1 - Monthly standard deviation [Divisor is (n-1)]
* New operator: yearvar1 - Yearly variance [Divisor is (n-1)]
* New operator: yearstd1 - Yearly standard deviation [Divisor is (n-1)]
2013-01-23 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* New operator: ensstd1 - Ensemble standard deviation [Divisor is (n-1)]
......
......@@ -82,7 +82,7 @@ void *Ensstat(void *argument)
cdoOperatorAdd("ensmean", func_mean, 0, NULL);
cdoOperatorAdd("ensavg", func_avg, 0, NULL);
cdoOperatorAdd("ensstd", func_std, 0, NULL);
cdoOperatorAdd("ensstd", func_std1, 0, NULL);
cdoOperatorAdd("ensstd1", func_std1, 0, NULL);
cdoOperatorAdd("ensvar", func_var, 0, NULL);
cdoOperatorAdd("ensvar1", func_var1, 0, NULL);
cdoOperatorAdd("enspctl", func_pctl, 0, NULL);
......
......@@ -131,7 +131,8 @@ void *Fldstat(void *argument)
}
if ( operfunc == func_mean || operfunc == func_avg ||
operfunc == func_var || operfunc == func_std )
operfunc == func_var || operfunc == func_std ||
operfunc == func_var1 || operfunc == func_std1 )
needWeights = TRUE;
streamID1 = streamOpenRead(cdoStreamName(0));
......
......@@ -24,35 +24,45 @@
Timstat timmean Time mean
Timstat timavg Time average
Timstat timvar Time variance
Timstat timvar1 Time variance [Divisor is (n-1)]
Timstat timstd Time standard deviation
Timstat timstd1 Time standard deviation [Divisor is (n-1)]
Hourstat hourmin Hourly minimum
Hourstat hourmax Hourly maximum
Hourstat hoursum Hourly sum
Hourstat hourmean Hourly mean
Hourstat houravg Hourly average
Hourstat hourvar Hourly variance
Hourstat hourvar1 Hourly variance [Divisor is (n-1)]
Hourstat hourstd Hourly standard deviation
Hourstat hourstd1 Hourly standard deviation [Divisor is (n-1)]
Daystat daymin Daily minimum
Daystat daymax Daily maximum
Daystat daysum Daily sum
Daystat daymean Daily mean
Daystat dayavg Daily average
Daystat dayvar Daily variance
Daystat dayvar1 Daily variance [Divisor is (n-1)]
Daystat daystd Daily standard deviation
Daystat daystd1 Daily standard deviation [Divisor is (n-1)]
Monstat monmin Monthly minimum
Monstat monmax Monthly maximum
Monstat monsum Monthly sum
Monstat monmean Monthly mean
Monstat monavg Monthly average
Monstat monvar Monthly variance
Monstat monvar1 Monthly variance [Divisor is (n-1)]
Monstat monstd Monthly standard deviation
Monstat monstd1 Monthly standard deviation [Divisor is (n-1)]
Yearstat yearmin Yearly minimum
Yearstat yearmax Yearly maximum
Yearstat yearsum Yearly sum
Yearstat yearmean Yearly mean
Yearstat yearavg Yearly average
Yearstat yearvar Yearly variance
Yearstat yearvar1 Yearly variance [Divisor is (n-1)]
Yearstat yearstd Yearly standard deviation
Yearstat yearstd1 Yearly standard deviation [Divisor is (n-1)]
*/
......@@ -94,41 +104,51 @@ void *Timstat(void *argument)
cdoInitialize(argument);
cdoOperatorAdd("timmin", func_min, 31, NULL);
cdoOperatorAdd("timmax", func_max, 31, NULL);
cdoOperatorAdd("timsum", func_sum, 31, NULL);
cdoOperatorAdd("timmean", func_mean, 31, NULL);
cdoOperatorAdd("timavg", func_avg, 31, NULL);
cdoOperatorAdd("timvar", func_var, 31, NULL);
cdoOperatorAdd("timstd", func_std, 31, NULL);
cdoOperatorAdd("yearmin", func_min, 10, NULL);
cdoOperatorAdd("yearmax", func_max, 10, NULL);
cdoOperatorAdd("yearsum", func_sum, 10, NULL);
cdoOperatorAdd("yearmean", func_mean, 10, NULL);
cdoOperatorAdd("yearavg", func_avg, 10, NULL);
cdoOperatorAdd("yearvar", func_var, 10, NULL);
cdoOperatorAdd("yearstd", func_std, 10, NULL);
cdoOperatorAdd("monmin", func_min, 8, NULL);
cdoOperatorAdd("monmax", func_max, 8, NULL);
cdoOperatorAdd("monsum", func_sum, 8, NULL);
cdoOperatorAdd("monmean", func_mean, 8, NULL);
cdoOperatorAdd("monavg", func_avg, 8, NULL);
cdoOperatorAdd("monvar", func_var, 8, NULL);
cdoOperatorAdd("monstd", func_std, 8, NULL);
cdoOperatorAdd("daymin", func_min, 6, NULL);
cdoOperatorAdd("daymax", func_max, 6, NULL);
cdoOperatorAdd("daysum", func_sum, 6, NULL);
cdoOperatorAdd("daymean", func_mean, 6, NULL);
cdoOperatorAdd("dayavg", func_avg, 6, NULL);
cdoOperatorAdd("dayvar", func_var, 6, NULL);
cdoOperatorAdd("daystd", func_std, 6, NULL);
cdoOperatorAdd("hourmin", func_min, 4, NULL);
cdoOperatorAdd("hourmax", func_max, 4, NULL);
cdoOperatorAdd("hoursum", func_sum, 4, NULL);
cdoOperatorAdd("hourmean", func_mean, 4, NULL);
cdoOperatorAdd("houravg", func_avg, 4, NULL);
cdoOperatorAdd("hourvar", func_var, 4, NULL);
cdoOperatorAdd("hourstd", func_std, 4, NULL);
cdoOperatorAdd("timmin", func_min, 31, NULL);
cdoOperatorAdd("timmax", func_max, 31, NULL);
cdoOperatorAdd("timsum", func_sum, 31, NULL);
cdoOperatorAdd("timmean", func_mean, 31, NULL);
cdoOperatorAdd("timavg", func_avg, 31, NULL);
cdoOperatorAdd("timvar", func_var, 31, NULL);
cdoOperatorAdd("timvar1", func_var1, 31, NULL);
cdoOperatorAdd("timstd", func_std, 31, NULL);
cdoOperatorAdd("timstd1", func_std1, 31, NULL);
cdoOperatorAdd("yearmin", func_min, 10, NULL);
cdoOperatorAdd("yearmax", func_max, 10, NULL);
cdoOperatorAdd("yearsum", func_sum, 10, NULL);
cdoOperatorAdd("yearmean", func_mean, 10, NULL);
cdoOperatorAdd("yearavg", func_avg, 10, NULL);
cdoOperatorAdd("yearvar", func_var, 10, NULL);
cdoOperatorAdd("yearvar1", func_var1, 10, NULL);
cdoOperatorAdd("yearstd", func_std, 10, NULL);
cdoOperatorAdd("yearstd1", func_std1, 10, NULL);
cdoOperatorAdd("monmin", func_min, 8, NULL);
cdoOperatorAdd("monmax", func_max, 8, NULL);
cdoOperatorAdd("monsum", func_sum, 8, NULL);
cdoOperatorAdd("monmean", func_mean, 8, NULL);
cdoOperatorAdd("monavg", func_avg, 8, NULL);
cdoOperatorAdd("monvar", func_var, 8, NULL);
cdoOperatorAdd("monvar1", func_var1, 8, NULL);
cdoOperatorAdd("monstd", func_std, 8, NULL);
cdoOperatorAdd("monstd1", func_std1, 8, NULL);
cdoOperatorAdd("daymin", func_min, 6, NULL);
cdoOperatorAdd("daymax", func_max, 6, NULL);
cdoOperatorAdd("daysum", func_sum, 6, NULL);
cdoOperatorAdd("daymean", func_mean, 6, NULL);
cdoOperatorAdd("dayavg", func_avg, 6, NULL);
cdoOperatorAdd("dayvar", func_var, 6, NULL);
cdoOperatorAdd("dayvar1", func_var1, 6, NULL);
cdoOperatorAdd("daystd", func_std, 6, NULL);
cdoOperatorAdd("daystd1", func_std1, 6, NULL);
cdoOperatorAdd("hourmin", func_min, 4, NULL);
cdoOperatorAdd("hourmax", func_max, 4, NULL);
cdoOperatorAdd("hoursum", func_sum, 4, NULL);
cdoOperatorAdd("hourmean", func_mean, 4, NULL);
cdoOperatorAdd("houravg", func_avg, 4, NULL);
cdoOperatorAdd("hourvar", func_var, 4, NULL);
cdoOperatorAdd("hourvar1", func_var1, 4, NULL);
cdoOperatorAdd("hourstd", func_std, 4, NULL);
cdoOperatorAdd("hourstd1", func_std1, 4, NULL);
operatorID = cdoOperatorID();
operfunc = cdoOperatorF1(operatorID);
......@@ -204,7 +224,7 @@ void *Timstat(void *argument)
vars1 = field_malloc(vlistID1, FIELD_PTR);
samp1 = field_malloc(vlistID1, FIELD_NONE);
if ( operfunc == func_std || operfunc == func_var )
if ( operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1 )
vars2 = field_malloc(vlistID1, FIELD_PTR);
tsID = 0;
......@@ -288,7 +308,7 @@ void *Timstat(void *argument)
samp1[varID][levelID].ptr[i]++;
}
if ( operfunc == func_std || operfunc == func_var )
if ( operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1 )
{
farsumq(&vars2[varID][levelID], field);
farsum(&vars1[varID][levelID], field);
......@@ -300,7 +320,7 @@ void *Timstat(void *argument)
}
}
if ( nsets == 0 && (operfunc == func_std || operfunc == func_var) )
if ( nsets == 0 && (operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1) )
for ( varID = 0; varID < nvars; varID++ )
{
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
......@@ -330,27 +350,29 @@ void *Timstat(void *argument)
fardiv(&vars1[varID][levelID], samp1[varID][levelID]);
}
}
else if ( operfunc == func_std || operfunc == func_var )
else if ( operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1 )
for ( varID = 0; varID < nvars; varID++ )
{
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for ( levelID = 0; levelID < nlevel; levelID++ )
{
double divisor = 0;
if ( operfunc == func_std1 || operfunc == func_var1 ) divisor = 1;
if ( samp1[varID][levelID].ptr == NULL )
{
if ( operfunc == func_std )
farcstd(&vars1[varID][levelID], vars2[varID][levelID], 1.0/nsets);
if ( operfunc == func_std || operfunc == func_std1 )
farcstdx(&vars1[varID][levelID], vars2[varID][levelID], nsets, divisor);
else
farcvar(&vars1[varID][levelID], vars2[varID][levelID], 1.0/nsets);
farcvarx(&vars1[varID][levelID], vars2[varID][levelID], nsets, divisor);
}
else
{
farinv(&samp1[varID][levelID]);
if ( operfunc == func_std )
farstd(&vars1[varID][levelID], vars2[varID][levelID], samp1[varID][levelID]);
if ( operfunc == func_std || operfunc == func_std1 )
farstdx(&vars1[varID][levelID], vars2[varID][levelID], samp1[varID][levelID], divisor);
else
farvar(&vars1[varID][levelID], vars2[varID][levelID], samp1[varID][levelID]);
farvarx(&vars1[varID][levelID], vars2[varID][levelID], samp1[varID][levelID], divisor);
}
}
}
......@@ -441,7 +463,7 @@ void *Timstat(void *argument)
field_free(vars1, vlistID1);
field_free(samp1, vlistID1);
if ( operfunc == func_std || operfunc == func_var )
if ( operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1 )
field_free(vars2, vlistID1);
if ( cdoDiag ) streamClose(streamID3);
......
......@@ -76,7 +76,7 @@ double fldfun(field_t field, int function)
else if ( function == func_rank ) rval = fldrank(field);
else if ( function == func_roc ) rval = fldroc(field);
else cdoAbort("function %d not implemented!", function);
else cdoAbort("%s: function %d not implemented!", __func__, function);
return rval;
}
......
......@@ -151,8 +151,12 @@ void farmin(field_t *field1, field_t field2);
void farmax(field_t *field1, field_t field2);
void farvar(field_t *field1, field_t field2, field_t field3);
void farstd(field_t *field1, field_t field2, field_t field3);
void farvarx(field_t *field1, field_t field2, field_t field3, double divisor);
void farstdx(field_t *field1, field_t field2, field_t field3, double divisor);
void farcvar(field_t *field1, field_t field2, double rconst1);
void farcstd(field_t *field1, field_t field2, double rconst1);
void farcvarx(field_t *field1, field_t field2, double rconst1, double divisor);
void farcstdx(field_t *field1, field_t field2, double rconst1, double divisor);
void farmoq(field_t *field1, field_t field2);
void faratan2(field_t *field1, field_t field2);
......
......@@ -35,7 +35,7 @@ void farfun(field_t *field1, field_t field2, int function)
else if ( function == func_mul ) farmul(field1, field2);
else if ( function == func_div ) fardiv(field1, field2);
else if ( function == func_atan2 ) faratan2(field1, field2);
else cdoAbort("function %d not implemented!", function);
else cdoAbort("%s: function %d not implemented!", __func__, function);
}
static
......@@ -439,6 +439,8 @@ void farvar(field_t *field1, field_t field2, field_t field3)
int nmiss2 = field2.nmiss;
double missval2 = field2.missval;
double *array2 = field2.ptr;
int nmiss3 = field3.nmiss;
double missval3 = field3.missval;
double *array3 = field3.ptr;
len = gridInqSize(grid1);
......@@ -446,16 +448,23 @@ void farvar(field_t *field1, field_t field2, field_t field3)
if ( len != gridInqSize(grid2) )
cdoAbort("Fields have different gridsize (%s)", __func__);
if ( nmiss1 > 0 || nmiss2 > 0 )
if ( nmiss1 > 0 || nmiss2 > 0 || nmiss3 > 0 )
{
for ( i = 0; i < len; i++ )
if ( !DBL_IS_EQUAL(array1[i], missval1) && !DBL_IS_EQUAL(array2[i], missval2) )
{
array1[i] = array2[i]*array3[i] - (array1[i]*array3[i])*(array1[i]*array3[i]);
if ( array1[i] < 0 && array1[i] > -1.e-5 ) array1[i] = 0;
}
else
array1[i] = missval1;
{
// array3[i] = 0.111111;
if ( i < 10 )
printf("array123 %g %g %g", array1[i], array2[i], array3[i]);
if ( !DBL_IS_EQUAL(array1[i], missval1) && !DBL_IS_EQUAL(array2[i], missval2) && !DBL_IS_EQUAL(array3[i], missval3) )
{
array1[i] = array2[i]*array3[i] - (array1[i]*array3[i])*(array1[i]*array3[i]);
if ( array1[i] < 0 && array1[i] > -1.e-5 ) array1[i] = 0;
}
else
array1[i] = missval1;
if ( i < 10 )
printf(" %g\n", array1[i]);
}
}
else
{
......@@ -476,6 +485,40 @@ void farvar(field_t *field1, field_t field2, field_t field3)
}
void farvarx(field_t *field1, field_t field2, field_t field3, double divisor)
{
long i, len;
int grid1 = field1->grid;
double missval1 = field1->missval;
double *array1 = field1->ptr;
int grid2 = field2.grid;
double missval2 = field2.missval;
double *array2 = field2.ptr;
double *array3 = field3.ptr;
double temp;
len = gridInqSize(grid1);
if ( len != gridInqSize(grid2) )
cdoAbort("Fields have different gridsize (%s)", __func__);
for ( i = 0; i < len; i++ )
{
temp = DIV( MUL(array1[i], array1[i]), array3[i]);
array1[i] = DIV( SUB(array2[i], temp), array3[i]-divisor);
if ( array1[i] < 0 && array1[i] > -1.e-5 ) array1[i] = 0;
}
field1->nmiss = 0;
for ( i = 0; i < len; i++ )
if ( DBL_IS_EQUAL(array1[i], missval1) || array1[i] < 0 )
{
array1[i] = missval1;
field1->nmiss++;
}
}
void farstd(field_t *field1, field_t field2, field_t field3)
{
long i, len;
......@@ -505,6 +548,35 @@ void farstd(field_t *field1, field_t field2, field_t field3)
}
void farstdx(field_t *field1, field_t field2, field_t field3, double divisor)
{
long i, len;
int grid1 = field1->grid;
double missval1 = field1->missval;
double *array1 = field1->ptr;
int grid2 = field2.grid;
len = gridInqSize(grid1);
if ( len != gridInqSize(grid2) )
cdoAbort("Fields have different gridsize (%s)", __func__);
farvarx(field1, field2, field3, divisor);
field1->nmiss = 0;
for ( i = 0; i < len; i++ )
if ( DBL_IS_EQUAL(array1[i], missval1) || array1[i] < 0 )
{
array1[i] = missval1;
field1->nmiss++;
}
else
{
array1[i] = IS_NOT_EQUAL(array1[i], 0) ? sqrt(array1[i]) : 0;
}
}
void farcvar(field_t *field1, field_t field2, double rconst1)
{
long i, len;
......@@ -516,22 +588,27 @@ void farcvar(field_t *field1, field_t field2, double rconst1)
int nmiss2 = field2.nmiss;
double missval2 = field2.missval;
double *array2 = field2.ptr;
int nmiss3 = 0;
len = gridInqSize(grid1);
if ( len != gridInqSize(grid2) )
cdoAbort("Fields have different gridsize (%s)", __func__);
if ( nmiss1 > 0 || nmiss2 > 0 )
if ( DBL_IS_EQUAL(rconst1, missval1) ) nmiss3 = 1;
if ( nmiss1 > 0 || nmiss2 > 0 || nmiss3 > 0 )
{
for ( i = 0; i < len; i++ )
if ( !DBL_IS_EQUAL(array1[i], missval1) && !DBL_IS_EQUAL(array2[i], missval2) )
{
array1[i] = array2[i]*rconst1 - (array1[i]*rconst1)*(array1[i]*rconst1);
if ( array1[i] < 0 && array1[i] > -1.e-5 ) array1[i] = 0;
}
else
array1[i] = missval1;
{
if ( !DBL_IS_EQUAL(array1[i], missval1) && !DBL_IS_EQUAL(array2[i], missval2) && nmiss3 == 0 )
{
array1[i] = array2[i]*rconst1 - (array1[i]*rconst1)*(array1[i]*rconst1);
if ( array1[i] < 0 && array1[i] > -1.e-5 ) array1[i] = 0;
}
else
array1[i] = missval1;
}
}
else
{
......@@ -552,6 +629,39 @@ void farcvar(field_t *field1, field_t field2, double rconst1)
}
void farcvarx(field_t *field1, field_t field2, double rconst1, double divisor)
{
long i, len;
int grid1 = field1->grid;
double missval1 = field1->missval;
double *array1 = field1->ptr;
int grid2 = field2.grid;
double missval2 = field2.missval;
double *array2 = field2.ptr;
double temp;
len = gridInqSize(grid1);
if ( len != gridInqSize(grid2) )
cdoAbort("Fields have different gridsize (%s)", __func__);
for ( i = 0; i < len; i++ )
{
temp = DIV( MUL(array1[i], array1[i]), rconst1);
array1[i] = DIV( SUB(array2[i], temp), rconst1-divisor);
if ( array1[i] < 0 && array1[i] > -1.e-5 ) array1[i] = 0;
}
field1->nmiss = 0;
for ( i = 0; i < len; i++ )
if ( DBL_IS_EQUAL(array1[i], missval1) || array1[i] < 0 )
{
array1[i] = missval1;
field1->nmiss++;
}
}
void farcstd(field_t *field1, field_t field2, double rconst1)
{
long i, len;
......@@ -581,6 +691,35 @@ void farcstd(field_t *field1, field_t field2, double rconst1)
}
void farcstdx(field_t *field1, field_t field2, double rconst1, double divisor)
{
long i, len;
int grid1 = field1->grid;
double missval1 = field1->missval;
double *array1 = field1->ptr;
int grid2 = field2.grid;
len = gridInqSize(grid1);
if ( len != gridInqSize(grid2) )
cdoAbort("Fields have different gridsize (%s)", __func__);
farcvarx(field1, field2, rconst1, divisor);
field1->nmiss = 0;
for ( i = 0; i < len; i++ )
if ( DBL_IS_EQUAL(array1[i], missval1) || array1[i] < 0 )
{
array1[i] = missval1;
field1->nmiss++;
}
else
{
array1[i] = IS_NOT_EQUAL(array1[i], 0) ? sqrt(array1[i]) : 0;
}
}
void farmoq(field_t *field1, field_t field2)
{
long i, len;
......
......@@ -27,7 +27,7 @@ void farcfun(field_t *field, double rconst, int function)
else if ( function == func_mul ) farcmul(field, rconst);
else if ( function == func_div ) farcdiv(field, rconst);
else if ( function == func_mod ) farmod(field, rconst);
else cdoAbort("function %d not implemented!", function);
else cdoAbort("%s: function %d not implemented!", __func__, function);
}
void farcmul(field_t *field, double rconst)
......
......@@ -435,11 +435,11 @@ void *Maggraph(void *argument);
#define HourpctlOperators {"hourpctl"}
#define TimselpctlOperators {"timselpctl"}
#define TimselstatOperators {"timselmin", "timselmax", "timselsum", "timselmean", "timselavg", "timselvar", "timselstd"}
#define TimstatOperators {"timmin", "timmax", "timsum", "timmean", "timavg", "timvar", "timstd"}
#define YearstatOperators {"yearmin", "yearmax", "yearsum", "yearmean", "yearavg", "yearvar", "yearstd"}
#define MonstatOperators {"monmin", "monmax", "monsum", "monmean", "monavg", "monvar", "monstd"}
#define DaystatOperators {"daymin", "daymax", "daysum", "daymean", "dayavg", "dayvar", "daystd"}
#define HourstatOperators {"hourmin", "hourmax", "hoursum", "hourmean", "houravg", "hourvar", "hourstd"}
#define TimstatOperators {"timmin", "timmax", "timsum", "timmean", "timavg", "timvar", "timvar1", "timstd", "timstd1"}
#define YearstatOperators {"yearmin", "yearmax", "yearsum", "yearmean", "yearavg", "yearvar", "yearvar1", "yearstd", "yearstd1"}
#define MonstatOperators {"monmin", "monmax", "monsum", "monmean", "monavg", "monvar", "monvar1", "monstd", "monstd1"}
#define DaystatOperators {"daymin", "daymax", "daysum", "daymean", "dayavg", "dayvar", "dayvar1", "daystd", "daystd1"}
#define HourstatOperators {"hourmin", "hourmax", "hoursum", "hourmean", "houravg", "hourvar", "hourvar1", "hourstd", "hourstd1"}
#define TimcorOperators {"timcor"}
#define TimcovarOperators {"timcovar"}
#define Timstat3Operators {"meandiff2test", "varquot2test"}
......
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