00001 #ident "$Header: /home/cvsuser/cvsroot/JSOC/proj/lev1.5_hmi/libs/lev15/fstats2.c,v 1.2 2011/11/23 21:54:23 couvidat Exp $"
00002
00003 #include <stdlib.h>
00004 #include <math.h>
00005 #include <float.h>
00006 #include <string.h>
00007
00008 #define NBINS 65536
00009 #define OK 0
00010 #define TOO_FEW_GOOD_POINTS -2
00011
00012 int fstats2(int n, float arr[], double *min, double *max, double *medn,
00013 double *mean, double *sig, double *skew, double *kurt, int *ngood)
00014 {
00015 int i, bin, nv = 0;
00016 float fmin = FLT_MAX, fmax = -FLT_MAX;
00017 double s = 0.0, s2 = 0.0, s3 = 0.0, s4 = 0.0, avg, var;
00018 static int hist[NBINS];
00019 double delta;
00020
00021 *min = *max = *medn = *mean = *sig = *skew = *kurt = __builtin_nan("");
00022
00023 for (i = 0; i < n; ++i) {
00024 float t = arr[i];
00025 if (isnanf(t))
00026 continue;
00027 if (fmin > t)
00028 fmin = t;
00029 if (fmax < t)
00030 fmax = t;
00031 s += t;
00032 ++nv;
00033 }
00034 *ngood = nv;
00035
00036 if (nv < 2) {
00037 if (nv == 1)
00038 *min = *max = *medn = *mean = s;
00039 return TOO_FEW_GOOD_POINTS;
00040 }
00041
00042 if (fmin == fmax) {
00043 *min = *max = *medn = *mean = fmin;
00044 *sig = *skew = *kurt = 0;
00045 return OK;
00046 }
00047
00048 delta = (fmax - fmin) / NBINS;
00049 memset(hist, 0, sizeof(int)*NBINS);
00050 avg = s / nv;
00051 for (i = 0; i < n; ++i) {
00052 double tmp = arr[i];
00053 double t = tmp - avg;
00054 double tt;
00055 if (isnan(tmp))
00056 continue;
00057 tt = t * t;
00058 s2 += tt;
00059 tt *= t;
00060 s3 += tt;
00061 tt *= t;
00062 s4 += tt;
00063 bin = (tmp - fmin) / delta;
00064 if (bin >= NBINS)
00065 bin = NBINS - 1;
00066 ++hist[bin];
00067 }
00068
00069 bin = 0;
00070 n = hist[bin];
00071 while (2*n < nv)
00072 n += hist[++bin];
00073 *medn = fmin + delta * (bin + 0.5);
00074
00075 *min = fmin;
00076 *max = fmax;
00077 *mean = avg;
00078 var = s2 / (nv - 1);
00079 *sig = sqrt(var);
00080 *skew = s3 / (var * *sig * nv);
00081 *kurt = s4 / (var * var * nv) - 3.0;
00082
00083 return OK;
00084 }
