Filling in Degrees in a Circle (Astronomy)

[Gerard flanagan]

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 * 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

 

[Lie]

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.

--
            Wayne Watson (Watson Adventures, Prop., Nevada City, CA)

              (121.015 Deg. W, 39.262 Deg. N) GMT-8 hr std. time)
               Obz Site:  39° 15' 7" N, 121° 2' 32" W, 2700 feet

                     Web Page: <www.speckledwithstars.net/>

Linear interpolation is simple, any mathematician would know. It's
just like this:

hA = height at point A
hB = height at point B
hN = height at point N (what we're finding)
l = horizontal distance between A and B
n = horizontal position measured from A

hN = hA + ((n / l) * (hB - hA))

so:

def interpolate(n, A, B):
# A, B is two-tuple of (angle, height) of
# the nearest known points surrounding n
oA, hA = A
oB, hB = B
l = oB - oA
return hA + ((n / l) * (hB - hA))