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