## File: polfit.x

package info (click to toggle)
iraf-rvsao 2.8.3-1
 `123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131` ``````# 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 ``````