/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:11 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "sstat1.h" #include /* PARAMETER translations */ #define FAC 64.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ sstat1( float xtab[], long nx, float stats[], long ihist[], long ncells, float x1, float x2) { long int i, index, j; float count, delta, prev, rscale, scale, sclnew, sumscl, temp, test, x, xmax, xmean, xmin; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Ihist = &ihist[0] - 1; float *const Stats = &stats[0] - 1; float *const Xtab = &xtab[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 1997-04-25 SSTAT1 Krogh Simplified code (WVS suggestion) *>> 1994-11-11 SSTAT1 Krogh Declared all vars. *>> 1994-10-20 SSTAT1 Krogh Changes to use M77CON *>> 1989-10-20 SSTAT1 CLL *>> 1987-05-01 SSTAT1 Lawson Initial code. *--S replaces "?": ?STAT1, ?STAT2 * * This subr computes basic statistics for X, storing them is * STATS() as follows: * * STATS(1) = Total count * STATS(2) = Min * STATS(3) = Max * STATS(4) = Mean * STATS(5) = Standard deviation * * This subr also accumulates counts in IHIST() to develop a * histogram of values of X. * * The data to be treated is given in XTAB(1:NX). If the * value of STATS(1) on entry is positive , say = COUNT, it is * assumed that COUNT data values have been processed previously * and results from that processing are in IHIST() and STATS(). * These results will be updated to reflect the additional set of * NX values. * * Alternatively, if STATS(1) is zero, the initial contents of * STATS(2:5) and IHIST() will be ignored and results will be * computed just for the present data set XTAB(1:NX). * * The user must specify the range and resolution of the histogram * by setting X1, X2, and NCELLS. The end cells, IHIST(1) and * IHIST(NCELLS) will be used to count occurences of X less than X1 * or greater than X2 respectively. * The cells IHIST(2) through IHIST(NCELLS-1) will * be used to count occurences of X in NCELLS-2 equal-length * subintervals of [X1, X2]. * * Define h = (X2 - X1)/(NCELLS-2). X-intervals will be * associated with elements of IHIST() as follows. * * X interval Counting cell * * (-Infinity, X1) IHIST(1) * [X1+(i-2)*h, X1+(i-1)*h) IHIST(i), i = 2,...,NCELLS-2 * [X2-h, X2] IHIST(NCELLS-1) * (X2, Infinity) IHIST(NCELLS) * * After use of this subroutine, the user can call * SSTAT2, to produce a printer-plot of the histogram and print the * statistics. * * Remark: It is more efficient to call this subroutine one * time giving it N points rather than calling it N times giving it * one point each time, but the results will be the same to within * arithmetic limitations either way. * ------------------------------------------------------------------ * Subroutine arguments * * XTAB() [in] Array of NX values whose statistics are to be * computed. * NX [in] Number of values given in XTAB(). * Require NX .ge. 1. * STATS() [inout] Array of length 5 into which statistics are or * will be stored. Initial value of STATS(1) must be positive * if IHIST() and STATS() contain prior results that are to be * updated. Otherwise the initial value of STATS(1) must be * zero. * IHIST() [inout] Integer array of length at least NCELLS into * which counts will be accumulated. * NCELLS [in] Total number of classification regions. * Require NCELLS .ge. 3. * X1,X2 [in] Lower and upper boundaries, respectively defining * the range of y values to be classified into NCELLS-2 equal * intervals. Require X1 < X2. * ------------------------------------------------------------------ * The value of FAC is not critical. It should be greater than * one. The program does less computation each time the test * (abs(DELTA) .lt. TEST) is satisfied. It will be true more * frequently if FAC is larger. There is probably not much advantage * in setting FAC larger than 4, so 64 is probably more than ample. * ------------------------------------------------------------------ * C. L. Lawson and S. Y. Chiu, JPL, Apr 1987. * 1989-10-20 CLL Moved integer declaration earlier to avoid warning * msg from Cray compiler. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ if (nx < 1) return; count = Stats[1]; if (count == 0.e0) { for (i = 1; i <= ncells; i++) { Ihist[i] = 0; } /* Begin: Special for first point * */ prev = 0.e0; x = Xtab[1]; xmin = x; xmax = x; xmean = 0.e0; test = -1.0e0; sumscl = 0.e0; /* End: Special for first point * */ } else { /* Here when COUNT is > zero on entry. */ xmin = Stats[2]; xmax = Stats[3]; xmean = Stats[4]; if (Stats[5] == 0.e0) { test = -1.0e0; sumscl = 0.e0; } else { /* STATS(5) contains the old value of Sigma. Since it is * nonzero (positive) here, COUNT must be at least 2. * */ scale = Stats[5]; rscale = 1.0e0/scale; test = FAC*scale; sumscl = count - 1.0e0; } } for (j = 1; j <= nx; j++) { prev = count; count += 1.0e0; x = Xtab[j]; xmin = fminf( x, xmin ); xmax = fmaxf( x, xmax ); /* . Begin: Tally in histogram. */ if (x < x1) { Ihist[1] += 1; } else if (x > x2) { Ihist[ncells] += 1; } else { /* Following stmt converts integer to float. */ temp = ncells - 2; /* Following stmt converts float to integer. */ index = temp*(x - x1)/(x2 - x1); Ihist[index + 2] += 1; } /* . End: Tally in histogram. * */ delta = x - xmean; /* Expect abs(DELTA) .le. TEST most of the time. * */ if (fabsf( delta ) > test) { if (delta == 0.e0) goto L_20; /* Here abs(DELTA) .gt. TEST and DELTA .ne. 0.E0 * Must compute new SCALE, RSCALE and TEST * and update SUMSCL if it is nonzero. * */ sclnew = fabsf( delta ); rscale = 1.0e0/sclnew; test = FAC*sclnew; if (sumscl != 0.e0) sumscl *= powif(scale*rscale,2); scale = sclnew; } xmean += delta/count; sumscl += (prev/count)*powif(delta*rscale,2); L_20: ; } Stats[1] = count; Stats[2] = xmin; Stats[3] = xmax; Stats[4] = xmean; if (prev == 0.e0) { Stats[5] = 0.e0; } else { Stats[5] = scale*sqrtf( sumscl/prev ); } return; } /* end of function */