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
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <math/gsurfit.h>
$if (datatype == r)
include "gsurfitdef.h"
$else
include "dgsurfitdef.h"
$endif
# GSINIT -- Procedure to initialize the surface descriptor.
$if (datatype == r)
procedure gsinit (sf, surface_type, xorder, yorder, xterms, xmin, xmax,
ymin, ymax)
$else
procedure dgsinit (sf, surface_type, xorder, yorder, xterms, xmin, xmax,
ymin, ymax)
$endif
pointer sf # surface descriptor
int surface_type # type of surface to be fitted
int xorder # x order of surface to be fit
int yorder # y order of surface to be fit
int xterms # presence of cross terms
PIXEL xmin # minimum value of x
PIXEL xmax # maximum value of x
PIXEL ymin # minimum value of y
PIXEL ymax # maximum value of y
int order
errchk malloc, calloc
begin
if (xorder < 1 || yorder < 1)
call error (0, "GSINIT: Illegal order.")
if (xmax <= xmin)
call error (0, "GSINIT: xmax <= xmin.")
if (ymax <= ymin)
call error (0, "GSINIT: ymax <= ymin.")
# allocate space for the gsurve descriptor
call calloc (sf, LEN_GSSTRUCT, TY_STRUCT)
# specify the surface-type dependent parameters
switch (surface_type) {
case GS_CHEBYSHEV, GS_LEGENDRE:
GS_XORDER(sf) = xorder
GS_YORDER(sf) = yorder
GS_NXCOEFF(sf) = xorder
GS_NYCOEFF(sf) = yorder
GS_XTERMS(sf) = xterms
switch (xterms) {
case GS_XNONE:
GS_NCOEFF(sf) = xorder + yorder - 1
case GS_XHALF:
order = min (xorder, yorder)
GS_NCOEFF(sf) = xorder * yorder - order * (order - 1) / 2
default:
GS_NCOEFF(sf) = xorder * yorder
}
GS_XRANGE(sf) = 2. / (xmax - xmin)
GS_XMAXMIN(sf) = - (xmax + xmin) / 2.
GS_YRANGE(sf) = 2. / (ymax - ymin)
GS_YMAXMIN(sf) = - (ymax + ymin) / 2.
case GS_POLYNOMIAL:
GS_XORDER(sf) = xorder
GS_YORDER(sf) = yorder
GS_NXCOEFF(sf) = xorder
GS_NYCOEFF(sf) = yorder
GS_XTERMS(sf) = xterms
switch (xterms) {
case GS_XNONE:
GS_NCOEFF(sf) = xorder + yorder - 1
case GS_XHALF:
order = min (xorder, yorder)
GS_NCOEFF(sf) = xorder * yorder - order * (order - 1) / 2
default:
GS_NCOEFF(sf) = xorder * yorder
}
GS_XRANGE(sf) = 1.0
GS_XMAXMIN(sf) = 0.0
GS_YRANGE(sf) = 1.0
GS_YMAXMIN(sf) = 0.0
default:
call error (0, "GSINIT: Unknown surface type.")
}
# set remaining parameters
GS_TYPE(sf) = surface_type
GS_XREF(sf) = INDEF
GS_YREF(sf) = INDEF
GS_ZREF(sf) = INDEF
GS_XMIN(sf) = xmin
GS_XMAX(sf) = xmax
GS_YMAX(sf) = ymax
GS_YMIN(sf) = ymin
# allocate space for the matrix and vectors
switch (surface_type ) {
case GS_LEGENDRE, GS_CHEBYSHEV, GS_POLYNOMIAL:
$if (datatype == r)
call calloc (GS_MATRIX(sf), GS_NCOEFF(sf) ** 2, TY_REAL)
call calloc (GS_CHOFAC(sf), GS_NCOEFF(sf) ** 2, TY_REAL)
call calloc (GS_VECTOR(sf), GS_NCOEFF(sf), TY_REAL)
call calloc (GS_COEFF(sf), GS_NCOEFF(sf), TY_REAL)
$else
call calloc (GS_MATRIX(sf), GS_NCOEFF(sf) ** 2, TY_DOUBLE)
call calloc (GS_CHOFAC(sf), GS_NCOEFF(sf) ** 2, TY_DOUBLE)
call calloc (GS_VECTOR(sf), GS_NCOEFF(sf), TY_DOUBLE)
call calloc (GS_COEFF(sf), GS_NCOEFF(sf), TY_DOUBLE)
$endif
default:
call error (0, "GSINIT: Unknown surface type.")
}
# initialize pointer to basis functions to null
GS_XBASIS(sf) = NULL
GS_YBASIS(sf) = NULL
GS_WZ(sf) = NULL
# set data points counter
GS_NPTS(sf) = 0
end
|
