Filling in Degrees in a Circle (Astronomy)

[Gerard flanagan]

W. eWatson wrote:
> The other night I surveyed a site for astronomical use by measuring the 
> altitude (0-90 degrees above the horizon) and az (azimuth, 0 degrees 
> north clockwise around the site to 360 degrees, almost north again) of 
> obstacles, trees. My purpose was to feed this profile of obstacles 
> (trees) to an astronomy program that would then account for not sighting 
> objects below the trees.
> 
> When I got around to entering them into the program by a file, I found 
> it required the alt at 360 azimuth points in order from 0 to 360 (same 
> as 0). Instead I have about 25 points, and expected the program to be 
> able to do simple linear interpolation between those.
> 
> Is there some simple operational device in Python that would allow me to 
> create an array (vector) of 360 points from my data by interpolating 
> between azimuth points when necessary? All my data I rounded to the 
> nearest integer. Maybe there's an interpolation operator?
> 
> As an example, supposed I had made 3 observations: (0,0) (180,45) and 
> (360,0). I would want some thing like (note the slope of the line from 0 
> to 179 is 45/180 or 0.25):
> alt: 0, 0.25, 0.50, 0.75, ... 44.75, 45.0
> az : 0, 1,    2,    3,              180
> 
> Of course, I don't need the az.
> 

Using one of the Python maths packages (scipy, sage, ...) is no doubt 
better, but out of interest here is some first-principles interpolation:

8<-----------------------------------------------------

class NewtonInterpolatingPolynomial(object):

     def __init__(self):
         self._domain = []
         self._codomain = []
         self._diffs = None
         self._coeffs = []

     def add_data_point(self, x, y):
         self._domain.append(x)
         self._codomain.append(y)
         if self._diffs is None:
             self._diffs = {}
         else:
             degree = len(self._domain) - 2
             _x = self._domain[-2]
             _y = self._codomain[-2]
             self._diffs[(degree, degree+1)] = (y - _y) / (x - _x)
             indices = range(degree+2)
             for t in ( tuple(indices[i:]) for i in 
reversed(indices[:-2]) ):
                 denominator = self._domain[t[0]] - self._domain[t[-1]]
                 k, _k = self._diffs[t[1:]], self._diffs[t[:-1]]
                 self._diffs[t] = (_k - k) / denominator
             self._coeffs.append(self._diffs[tuple(indices)])

     def __str__(self):
         N = len(self._domain)
         if not N:
             return ''
         parts = [str(self._codomain[0])]
         multipliers = [''.join(('(X - ',str(C),')')) for C in 
self._domain[:-1]]
         for i, k in enumerate(self._coeffs):
             parts.append('*'.join([str(k)] + multipliers[:i+1]))
         return ' + '.join(parts)

     def interpolate(self, gamma):
         #return eval(str(self).replace('X', str(gamma)))
         ret = self._codomain[0]
         K = 1
         multipliers = [gamma-C for C in self._domain[:-1]]
         for i, k in enumerate(multipliers):
             K *= k
             ret += self._coeffs[i] * K
         return ret

     def __call__(self, x):
         return self.interpolate(x)


rawdata = '''
0    18
18   18
27   16
34   20
48   20
59  28
72  32
'''

data = [map(float, line.split() ) for line in rawdata.splitlines() if line]

newton = NewtonInterpolatingPolynomial()

for x, y in data:
     newton.add_data_point(x, y)
print newton
for P in range(80):
     print P, '        -> ', newton(P)

18.0 + 0.0*(X - 0.0) + -0.0082304526749*(X - 0.0)*(X - 18.0) + 
0.00170098903759*(X - 0.0)*(X - 18.0)*(X - 27.0) + -8.87803681143e-05*(X 
- 0.0)*(X - 18.0)*(X - 27.0)*(X - 34.0) + 3.29057245545e-06*(X - 0.0)*(X 
- 18.0)*(X - 27.0)*(X - 34.0)*(X - 48.0) + -8.98633510787e-08*(X - 
0.0)*(X - 18.0)*(X - 27.0)*(X - 34.0)*(X - 48.0)*(X - 59.0)
0         ->  18.0
1         ->  26.0156268043
2         ->  31.8038369501
3         ->  35.719116702
4         ->  38.0797664434
5         ->  39.1701395406
6         ->  39.2428165062
7         ->  38.5207144606
8         ->  37.1991318912
9         ->  35.4477287116
10         ->  33.4124416172
11         ->  31.2173347409
12         ->  28.9663856061
13         ->  26.7452063785
14         ->  24.6227004166
15         ->  22.6526541194
16         ->  20.8752640745
17         ->  19.3185995022
18         ->  18.0
19         ->  16.9274085845
20         ->  16.1006400316
21         ->  15.5125845157
22         ->  15.1503465468
23         ->  14.996319206
24         ->  15.0291936797
25         ->  15.2249040921
26         ->  15.5575076355
27         ->  16.0
28         ->  16.525066101
29         ->  17.1057661048
30         ->  17.7161567532
31         ->  18.3318479864
32         ->  18.9304948638
33         ->  19.4922247833
34         ->  20.0
35         ->  20.4399154416
36         ->  20.8014318237
37         ->  21.0775440624
38         ->  21.2648849859
39         ->  21.3637643445
40         ->  21.3781431185
41         ->  21.3155431251
42         ->  21.1868919232
43         ->  21.0063030167
44         ->  20.7907913565
45         ->  20.5599241405
46         ->  20.3354069122
47         ->  20.1406049574
48         ->  20.0
49         ->  19.9385821951
50         ->  19.9811774214
51         ->  20.1517098718
52         ->  20.472399942
53         ->  20.962897418
54         ->  21.6393499612
55         ->  22.5134068932
56         ->  23.5911582776
57         ->  24.872009301
58         ->  26.3474899523
59         ->  28.0
60         ->  29.8014892685
61         ->  31.712073212
62         ->  33.6785837876
63         ->  35.6330556265
64         ->  37.4911475029
65         ->  39.1504991026
66         ->  40.4890230889
67         ->  41.3631324673
68         ->  41.6059032486
69         ->  41.0251724103
70         ->  39.4015711567
71         ->  36.4864934765
72         ->  32.0
73         ->  25.6286571539
74         ->  17.0233116147
75         ->  5.79680006015
76         ->  -8.47840578009
77         ->  -26.2726187763
78         ->  -48.1014207514
79         ->  -74.5282115362