File: legendre.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 (217 lines) | stat: -rw-r--r-- 4,269 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
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
# File rvsao/Util/legendre.x
# April 13, 1994
# By Doug Mink, Harvard-Smithsonian Center for Astrophysics
# After John Tonry (11/19/82) and Guillermo Torres (Jan/1989)

# Copyright(c) 1994 Smithsonian Astrophysical Observatory
# You may do anything you like with this file except remove this copyright.
# The Smithsonian Astrophysical Observatory makes no representations about
# the suitability of this software for any purpose.  It is provided "as is"
# without express or implied warranty.

# All sorts of routines to fit and evaluate legendre polynomials, and to
# convert coefficients to those of plain polynomials

procedure fitlpo (npt,x,y,scale,ncoeff,coeff)

# x and y are data arrays with npt points to fit.
# scale receives xMin and xMax, and returns a and b.
# where z = a * (x - b)
# a polynomial of degree ncoeff-1 is fit and the coefficients are
# returned in coeff.

int	npt
real	x[npt]
real	y[npt]
real	scale[ARB]
int	ncoeff
double	coeff[ncoeff]

double	cov[8*8], vec[8]
double	a, b, z, lpi, det
int	i,j,n,maxfit, ij, ji, irst[2,8]

begin
	maxfit = 8

	a = 2.0d0 / (scale[2] - scale[1])
	b = 0.5d0 * (scale[2] + scale[1])

	do j = 1, ncoeff {
	    vec[J] = 0
	    }
	do i = 1, ncoeff*ncoeff {
	    cov[i] = 0
	    }

	do n = 1, npt {
	    z = a * (double (x[n]) - b)

	    do j = 1, ncoeff {
		coeff[j] = lpi (j,z)
		vec[j] = vec[j] + coeff[j] * y[n]
		do i = 1, j {
		    ji = (j - 1) * ncoeff + i
		    cov[ji] = cov[ji] + coeff[i] * coeff[j]
		   }
		}
	    }

	do j = 2, ncoeff {
	    do i = 1, j-1 {
		ij = (i - 1) * ncoeff + j
		ji = (j - 1) * ncoeff + i
		cov[ij] = cov[ji]
		}
	    }

	call invert (ncoeff,cov,irst,det)

	if (det == 0) {
	    call printf ("fitlpoly: singular covariance matrix\n")
	    return	
	    }

	do j = 1, ncoeff {
	    coeff[j] = 0
	    do i = 1, ncoeff {
		coeff[j] = coeff[j] + cov[(j-1)*ncoeff+i] * vec[i]
		}
	    }

	scale[1] = a
	scale[2] = b

	return
end


real procedure lpoly (x,scale,ncoeff,coeff)

real	x
double	coeff[ARB]
int	ncoeff
real	scale[ARB]

double	temp, z, lpi()
int	i

begin
	z = scale[1] * (x - scale[2])
	temp = 0
	do i = 1, ncoeff {
	    temp = temp + coeff[i] * lpi (i,z)
	    }
	lpoly = temp
	return
end


# legendre polynomial of order i-1

double procedure lpi (i,z)

int	i
double	z

begin
	switch (i) {
	    case 1:
		lpi = 1.d0
	    case 2:
		lpi = z
	    case 3:
		lpi =  -.5d0 + 1.5d0*z*z
	    case 4:
		lpi = z*(-1.5d0 + 2.5d0*z*z)
	    case 5:
		lpi = .375d0 + z*z*(-3.75d0 + 4.375d0*z*z)
	    case 6:
		lpi = z*(1.875d0 + z*z*(-8.75d0 + 7.875*z*z))
	    case 7:
		lpi = -.3125d0 + z*z*(6.5625d0 + z*z*(-19.6875d0 +
		      14.4375d0*z*z))
	    case 8:
		lpi = z*(2.1875d0 + z*z*(19.6875d0 + z*z*(-43.3125d0 +
		      26.8125d0*z*z)))
	    default:
	    }

	return
end


procedure polyc (nc,scale,coeff,pcoeff)

# routine to convert coefficients of a 1-d legendre polynomial to 
# coefficients of a plain polynomial
# The nc coefficients are stored with the low order coefficient first

int	nc
real	scale[2]
double	coeff[ARB]
double	pcoeff[ARB]

double	sc[2], temp
int	i,k,ibc()

# Legendre polynomials 0 - 7

double lcoeff[8,8]
data lcoeff/ 1d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,
	     0d0,1d0,0d0,0d0,0d0,0d0,0d0,0d0,
	   -.5d0,0d0,1.5d0,0d0,0d0,0d0,0d0,0d0,
	     0d0,-1.5d0,0d0,2.5d0,0d0,0d0,0d0,0d0,
	  .375d0,0d0,-3.75d0,0d0,4.375d0,0d0,0d0,0d0,
	     0d0,1.875d0,0d0,-8.75d0,0d0,7.875,0d0,0d0,
	-.3125d0,0d0,6.5625d0,0d0,-19.6875d0,0d0,14.4375d0,0d0,
	     0d0,2.1875d0,0d0,19.6875d0,0d0,-43.3125d0,0d0,26.8125d0/

begin
	sc[1] = double (scale[1])
	sc[2] = double (scale[2])
	do k = 0, nc-1 {
	    temp = 0
	    do i = k, nc-1 {
	    	temp = temp + coeff[i+1] * lcoeff[k+1,i+1]
		}
	    pcoeff[k+1] = temp
	    }

	do k = 0, nc-1 {
	    temp = 0
	    do i = k, nc-1 {
	        if (sc(2) != 0)
	    	    temp = temp + pcoeff[i+1] * sc[1]**i * 
     			   (-sc[2])**(i-k) * double (ibc(k,i))
	        else
	    	    temp = temp + pcoeff[i+1] * sc[1]**i
		}
	    pcoeff[k+1] = temp
	    }

	return
end


#  return the (m n) binomial coefficient

int procedure ibc (m, n)

int	m,n

int	iret,i

begin
	iret = 1
	do i = m+1, n {
	    iret = iret * i
	    }
	do i = 2, n-m {
	    iret = iret / i
	    }
	return iret
end
# Sep 17 1991	New program

# Apr 13 1994	Add working array for invert