aboutsummaryrefslogtreecommitdiff
path: root/math/curfit/cvinit.gx
blob: f3518dabcf9764994f882f0b85a48d2e25da1639 (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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.

include <math/curfit.h>
include	<mach.h>

$if (datatype == r)
include "curfitdef.h"
$else
include "dcurfitdef.h"
$endif

# CVINIT --  Procedure to initialize the curve  descriptor.

$if (datatype == r)
procedure cvinit (cv, curve_type, order, xmin, xmax)
$else
procedure dcvinit (cv, curve_type, order, xmin, xmax)
$endif

pointer	cv		# curve descriptor
int	curve_type	# type of curve to be fitted
int	order		# order of curve to be fitted, or in the case of the
			# spline the number of polynomial pieces to be fit
PIXEL	xmin		# minimum value of x
PIXEL	xmax		# maximum value of x

errchk	malloc, calloc

begin
	# check for bad parameters.
	cv = NULL
	if (order < 1)
	    call error (0, "CVINIT: Illegal order.")

	if (xmax <= xmin)
	    call error (0, "CVINIT: xmax <= xmin.")

	# allocate space for the curve descriptor
	call calloc (cv, LEN_CVSTRUCT, TY_STRUCT)

	# specify the curve-type dependent parameters
	switch (curve_type) {
	case CHEBYSHEV, LEGENDRE:
	    CV_ORDER(cv) = order
	    CV_NCOEFF(cv) = order
	    CV_RANGE(cv) = 2. / (xmax - xmin)
	    CV_MAXMIN(cv) = - (xmax + xmin) / 2.
	case SPLINE3:
	    CV_ORDER(cv) = SPLINE3_ORDER
	    CV_NCOEFF(cv) = order + SPLINE3_ORDER - 1
	    CV_NPIECES(cv) = order - 1
	    CV_SPACING(cv) = order / (xmax - xmin)
	case SPLINE1:
	    CV_ORDER(cv) = SPLINE1_ORDER
	    CV_NCOEFF(cv) = order + SPLINE1_ORDER - 1
	    CV_NPIECES(cv) = order - 1
	    CV_SPACING(cv) = order / (xmax - xmin)
	case USERFNC:
	    CV_ORDER(cv) = order
	    CV_NCOEFF(cv) = order
	    # Prevent abort for non-linear userfnc, where these values
	    #  may be arbitrary arguments to pass to external.
	    if ( abs(xmax-xmin) > EPSILON ) {
		CV_RANGE(cv) = 2. / (xmax - xmin)
	    } else {
		CV_RANGE(cv) = 0.
	    }
	    CV_MAXMIN(cv) = - (xmax + xmin) / 2.
	default:
	    call error (0, "CVINIT: Unknown curve type.")
	}

	# set remaining parameters
	CV_TYPE(cv) = curve_type
	CV_XMIN(cv) = xmin
	CV_XMAX(cv) = xmax

	# allocate space for the matrix and vectors
	call calloc (CV_XBASIS(cv), CV_ORDER(cv), TY_PIXEL)
	call calloc (CV_MATRIX(cv), CV_ORDER(cv)*CV_NCOEFF(cv), TY_PIXEL)
	call calloc (CV_CHOFAC(cv), CV_ORDER(cv)*CV_NCOEFF(cv), TY_PIXEL)
	call calloc (CV_VECTOR(cv), CV_NCOEFF(cv), TY_PIXEL)
	call calloc (CV_COEFF(cv), CV_NCOEFF(cv), TY_PIXEL)

	# initialize pointer to basis functions to null
	CV_BASIS(cv) = NULL
	CV_WY(cv) = NULL
	CV_LEFT(cv) = NULL

	# set null user function
	CV_USERFNC(cv) = NULL

	# set data points counter
	CV_NPTS(cv) = 0
end