1 subroutine SSTAT1(XTAB, NX, STATS, IHIST, NCELLS, X1, X2) 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5C>> 1997-04-25 SSTAT1 Krogh Simplified code (WVS suggestion) 6C>> 1994-11-11 SSTAT1 Krogh Declared all vars. 7C>> 1994-10-20 SSTAT1 Krogh Changes to use M77CON 8C>> 1989-10-20 SSTAT1 CLL 9C>> 1987-05-01 SSTAT1 Lawson Initial code. 10c--S replaces "?": ?STAT1, ?STAT2 11c 12c This subr computes basic statistics for X, storing them is 13c STATS() as follows: 14c 15c STATS(1) = Total count 16c STATS(2) = Min 17c STATS(3) = Max 18c STATS(4) = Mean 19c STATS(5) = Standard deviation 20c 21c This subr also accumulates counts in IHIST() to develop a 22c histogram of values of X. 23c 24c The data to be treated is given in XTAB(1:NX). If the 25c value of STATS(1) on entry is positive , say = COUNT, it is 26c assumed that COUNT data values have been processed previously 27c and results from that processing are in IHIST() and STATS(). 28c These results will be updated to reflect the additional set of 29c NX values. 30c 31c Alternatively, if STATS(1) is zero, the initial contents of 32c STATS(2:5) and IHIST() will be ignored and results will be 33c computed just for the present data set XTAB(1:NX). 34c 35c The user must specify the range and resolution of the histogram 36c by setting X1, X2, and NCELLS. The end cells, IHIST(1) and 37c IHIST(NCELLS) will be used to count occurences of X less than X1 38c or greater than X2 respectively. 39c The cells IHIST(2) through IHIST(NCELLS-1) will 40c be used to count occurences of X in NCELLS-2 equal-length 41c subintervals of [X1, X2]. 42c 43c Define h = (X2 - X1)/(NCELLS-2). X-intervals will be 44c associated with elements of IHIST() as follows. 45c 46c X interval Counting cell 47c 48c (-Infinity, X1) IHIST(1) 49c [X1+(i-2)*h, X1+(i-1)*h) IHIST(i), i = 2,...,NCELLS-2 50c [X2-h, X2] IHIST(NCELLS-1) 51c (X2, Infinity) IHIST(NCELLS) 52c 53c After use of this subroutine, the user can call 54c SSTAT2, to produce a printer-plot of the histogram and print the 55c statistics. 56c 57c Remark: It is more efficient to call this subroutine one 58c time giving it N points rather than calling it N times giving it 59c one point each time, but the results will be the same to within 60c arithmetic limitations either way. 61c ------------------------------------------------------------------ 62c Subroutine arguments 63c 64c XTAB() [in] Array of NX values whose statistics are to be 65c computed. 66c NX [in] Number of values given in XTAB(). 67c Require NX .ge. 1. 68c STATS() [inout] Array of length 5 into which statistics are or 69c will be stored. Initial value of STATS(1) must be positive 70c if IHIST() and STATS() contain prior results that are to be 71c updated. Otherwise the initial value of STATS(1) must be 72c zero. 73c IHIST() [inout] Integer array of length at least NCELLS into 74c which counts will be accumulated. 75c NCELLS [in] Total number of classification regions. 76c Require NCELLS .ge. 3. 77c X1,X2 [in] Lower and upper boundaries, respectively defining 78c the range of y values to be classified into NCELLS-2 equal 79c intervals. Require X1 < X2. 80c ------------------------------------------------------------------ 81c The value of FAC is not critical. It should be greater than 82c one. The program does less computation each time the test 83c (abs(DELTA) .lt. TEST) is satisfied. It will be true more 84c frequently if FAC is larger. There is probably not much advantage 85c in setting FAC larger than 4, so 64 is probably more than ample. 86c ------------------------------------------------------------------ 87c C. L. Lawson and S. Y. Chiu, JPL, Apr 1987. 88C 1989-10-20 CLL Moved integer declaration earlier to avoid warning 89c msg from Cray compiler. 90c ------------------------------------------------------------------ 91 integer NCELLS, NX 92 integer IHIST(NCELLS) 93 integer I, INDEX, J 94 real COUNT, DELTA, FAC, PREV 95 real SCALE, RSCALE, SCLNEW, STATS(5), SUMSCL 96 real TEMP, TEST, X, X1, X2, XMAX, XMEAN, XMIN 97 real XTAB(NX) 98 parameter(FAC = 64.0E0) 99c ------------------------------------------------------------------ 100 if(NX .lt. 1) return 101 COUNT = STATS(1) 102 if(COUNT .eq. 0.E0) then 103 do 10 I=1,NCELLS 104 IHIST(I) = 0 105 10 continue 106c 107c Begin: Special for first point 108c 109 PREV = 0.E0 110 X = XTAB(1) 111 XMIN = X 112 XMAX = X 113 XMEAN = 0.E0 114 TEST = -1.0E0 115 SUMSCL = 0.E0 116c End: Special for first point 117c 118 else 119c Here when COUNT is > zero on entry. 120 XMIN = STATS(2) 121 XMAX = STATS(3) 122 XMEAN = STATS(4) 123c 124 if(STATS(5) .eq. 0.E0) then 125 TEST = -1.0E0 126 SUMSCL = 0.E0 127 else 128c 129c STATS(5) contains the old value of Sigma. Since it is 130c nonzero (positive) here, COUNT must be at least 2. 131c 132 SCALE = STATS(5) 133 RSCALE = 1.0E0/SCALE 134 TEST = FAC * SCALE 135 SUMSCL = COUNT - 1.0E0 136 endif 137 endif 138c 139 do 30 J = 1, NX 140 PREV = COUNT 141 COUNT = COUNT + 1.0E0 142 X = XTAB(J) 143 XMIN = min(X, XMIN) 144 XMAX = max(X, XMAX) 145c . Begin: Tally in histogram. 146 if(X .lt. X1) then 147 IHIST(1) = IHIST(1) + 1 148 elseif(X .gt. X2) then 149 IHIST(NCELLS) = IHIST(NCELLS) + 1 150 else 151c Following stmt converts integer to float. 152 TEMP = NCELLS-2 153c Following stmt converts float to integer. 154 INDEX = TEMP*(X-X1)/(X2-X1) 155 IHIST(INDEX + 2) = IHIST(INDEX + 2) + 1 156 endif 157c . End: Tally in histogram. 158c 159 DELTA = X - XMEAN 160c 161c Expect abs(DELTA) .le. TEST most of the time. 162c 163 if(abs(DELTA) .gt. TEST) then 164 if( DELTA .eq. 0.E0 ) go to 20 165c 166c Here abs(DELTA) .gt. TEST and DELTA .ne. 0.E0 167c Must compute new SCALE, RSCALE and TEST 168c and update SUMSCL if it is nonzero. 169c 170 SCLNEW = abs(DELTA) 171 RSCALE = 1.0E0 / SCLNEW 172 TEST = FAC * SCLNEW 173 if(SUMSCL .ne. 0.E0) SUMSCL = SUMSCL * (SCALE * RSCALE)**2 174 SCALE = SCLNEW 175 endif 176 XMEAN = XMEAN + DELTA / COUNT 177 SUMSCL = SUMSCL + (PREV/COUNT) * (DELTA*RSCALE)**2 178 20 continue 179 30 continue 180c 181 STATS(1) = COUNT 182 STATS(2) = XMIN 183 STATS(3) = XMAX 184 STATS(4) = XMEAN 185 if(PREV .eq. 0.E0) then 186 STATS(5) = 0.E0 187 else 188 STATS(5) = SCALE * sqrt( SUMSCL / PREV ) 189 endif 190 return 191 end 192