aboutsummaryrefslogtreecommitdiff
path: root/pkg/xtools/ranges/rgbinr.x
blob: 81fb9f70315a0828d32d7d1e77064455b38f58cf (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.

include	<pkg/rg.h>

# RG_BIN -- Average or median of data.
#
# The ranges are broken up into subranges of at most abs (nbin) points.  The
# subranges are averaged if nbin > 1 and medianed if nbin < 1.
# The output array must be large enough to contain the desired points.
# If the ranges are merged then the input and output arrays may be the same.

procedure rg_binr (rg, nbin, in, nin, out, nout)

pointer	rg				# Ranges
int	nbin				# Maximum points in average or median
real	in[nin]				# Input array
int	nin				# Number of input points
real	out[ARB]			# Output array
int	nout				# Number of output points

int	i, j, k, n, npts, ntemp

real	asumr(), amedr()

errchk	rg_packr

begin
	# Error check the range pointer.

	if (rg == NULL)
	    call error (0, "Range descriptor undefined")

	# If the bin size is exactly one then move the selected input points
	# to the output array.

	if (abs (nbin) == 1) {
	    call rg_packr (rg, in, out)
	    return
	}

	# Determine the subranges and take the median or average.

	npts = abs (nbin)
	ntemp = 0

	do i = 1, RG_NRGS(rg) {
	    if (RG_X1(rg, i) > RG_X2(rg, i)) {
		j = min (nin, RG_X1(rg, i))
		k = max (1, RG_X2(rg, i))
		while (j >= k) {
		    n = max (0, min (npts, j - k + 1))
		    k = k - n
		    ntemp = ntemp + 1
		    if (nbin > 0)
		        out[ntemp] = asumr (in[k + 1], n) / n
		    else
			out[ntemp] = amedr (in[k+1], n)
		}
	    } else {
		j = max (1, RG_X1(rg, i))
		k = min (nin, RG_X2(rg, i))
	        while (j <= k) {
		    n = max (0, min (npts, k - j + 1))
		    ntemp = ntemp + 1
		    if (nbin > 0)
		        out[ntemp] = asumr (in[j], n) / n
		    else
			out[ntemp] = amedr (in[j], n)
		    j = j + n
		}
	    }
	}

	nout = ntemp
end