00001 #ident "$Header: /home/cvsuser/cvsroot/JSOC/proj/lev1.5_hmi/libs/lev15/dstats2.c,v 1.2 2011/11/23 21:54:33 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 dstats2(int n, double 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 double dmin = DBL_MAX, dmax = -DBL_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 double t = arr[i];
00025 if (isnan(t))
00026 continue;
00027 if (dmin > t)
00028 dmin = t;
00029 if (dmax < t)
00030 dmax = 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 (dmin == dmax) {
00043 *min = *max = *medn = *mean = dmin;
00044 *sig = *skew = *kurt = 0;
00045 return OK;
00046 }
00047
00048 delta = (dmax - dmin) / 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 - dmin) / 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 = dmin + delta * (bin + 0.5);
00074
00075 *min = dmin;
00076 *max = dmax;
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 }
