File: polfit.x

package info (click to toggle)
iraf-rvsao 2.8.3-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 16,456 kB
  • sloc: ansic: 963; lisp: 651; fortran: 397; makefile: 27
file content (131 lines) | stat: -rw-r--r-- 2,807 bytes parent folder | download
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
125
126
127
128
129
130
131
# File rvsao/Util/polfit.x
# October 31, 1991
# By Doug Mink, Harvard-Smithsonian Center for Astrophysics
# After Bevington, page 141

# Polynomial least squares fitting program, almost identical to the
# one in Bevington, "Data Reduction and Error Analysis for the
# Physical Sciences," page 141.  I changed the argument list and
# removed the weighting:  y = a(1) + a(2)*x + a(3)*x**2 + a(3)*x**3 + . . .

procedure polfit (x,y,npts,nterms,a,chisqr)

double	x[ARB]		# Array of independent variable points
double	y[ARB]		# Array of dependent variable points
int	npts		# Number of data points to fit
int	nterms		# Number of parameters to fit
double	a[ARB]		# Vector containing current fit values
double	chisqr

double xterm,yterm,delta,chisq,determ()
double sumx[19],sumy[10],array[10,10],xi,yi
int	i,j,k,l,n,nmax

begin
	nmax = (2 * nterms) - 1 
	call aclrd (sumx,nmax)
	call aclrd (sumy,nterms)

# accumulate weighted sums
	chisq = 0.d0
	do i = 1, npts {
	    xi = x[i] 
	    yi = y[i] 
	    xterm = 1.d0
	    do n = 1, nmax {
		sumx[n] = sumx[n] + xterm
		xterm = xterm * xi
		}
	    yterm = yi 
	    do n = 1, nterms {
		sumy[n] = sumy[n] + yterm
		yterm = yterm * xi
		}
	    chisq = chisq + yi*yi
	    }

# construct matrices and calculate coeffients
	do j = 1, nterms {
	    do k = 1, nterms {
		n = j + k - 1 
		array[j,k] = sumx[n]
		}
	    }
	delta =  determ (array,nterms)
	if (delta == 0.d0) {
	    chisqr = 0.d0
	    call aclrd (a, nterms)
	    return
	    }

	do l = 1, nterms {
	    do j = 1, nterms {
		do k = 1, nterms {
		    n = j + k - 1 
		    array[j,k] = sumx[n]
		    }
		array[j,l] = sumy[j]
		}
	    a[l] = determ (array,nterms) / delta
	    }

# calculate chi square
	do j = 1, nterms {
	    chisq = chisq - (2.d0 * a[j] * sumy[j])
	    do k = 1, nterms {
		n = j + k - 1
		chisq = chisq + (a[j] * a[k] * sumx[n])
		}
	    }
	chisqr = chisq / (npts - nterms)

	return
end


#--- calculate the determinant of a square matrix
#    this subprogram destroys the input matrix array
#    from bevington, page 294.

double procedure determ (array,norder)

double array[10,10]	# Input matrix array
int norder		# Order of determinant (degree of matrix)

double	save,det
int	i,j,k,k1

begin
	det = 1.d0
	do k = 1, norder {

	# Interchange columns if diagnoal element is zero
	    if (array[k,k] == 0) {
		do j = k, norder {
		    if (array[k,j] == 0) {
			det = 0.d0
			return det
			}
		    }
		do i = k, norder {
		    save = array[i,j] 
		    array[i,j] = array[i,k]
		    array[i,k] = save 
		    }
		det = -det
		}

	# Subtract row k from lower rows to get diagonal matrix 
	    det = det * array[k,k]
	    if (k < norder) {
		k1 = k+1
		do i = k1, norder {
		    do j = k1, norder {
			array[i,j] = array[i,j]-array[i,k]*array[k,j]/array[k,k]
			}
		    }
		}
	    }

	return det
end