-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathtoypython.py
129 lines (112 loc) · 4.29 KB
/
toypython.py
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
#
# coding=utf8
# Copyright 2013, 2014 Bence Béky
#
# This file is part of Spotrod.
#
# Spotrod is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# Spotrod is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with Spotrod. If not, see <http://www.gnu.org/licenses/>.
#
import numpy;
def circleangleloop(r, p, z):
"""circleangleloop(r, p, z)
Calculate half central angle of the arc of circle of radius r
(which concentrically spans the inside of the star during integration)
that is inside a circle of radius p (planet)
with separation of centers z.
This is a zeroth order homogeneous function, that is,
circleangle(alpha*r, alpha*p, alpha*z) = circleangle(r, p, z).
This version uses a loop over r.
Input:
r one dimensional numpy array
p scalar
z scalar
They should all be non-negative, but there is no other restriction.
Output:
circleangle one dimensional numpy array, same size as r
"""
# If the circle arc of radius r is disjoint from the circular disk
# of radius p, then the angle is zero.
answer = numpy.zeros_like(r);
pminusz = p-z;
pplusz = p+z;
zsquared = z*z;
psquared = p*p;
for i in xrange(r.shape[0]):
ri = r[i];
# If the planet entirely covers the circle, the half central angle is pi.
if (ri <= pminusz):
answer[i] = numpy.pi;
# If the triangle inequalities hold between z, r, and p,
# then we have partial overlap.
# If alpha is the half central angle in the triangle with sides r, p, and z,
# with p opposite the angle, then p^2 = r^2 + z^2 - 2 rz cos(alpha)
elif (ri < pplusz) & (ri > -pminusz):
answer[i] = numpy.arccos((ri*ri+zsquared-psquared)/(2*z*ri));
return answer;
def circleanglemask(r, p, z):
"""circleanglemask(r, p, z)
Calculate half central angle of the arc of circle of radius r
(which concentrically spans the inside of the star during integration)
that is inside a circle of radius p (planet)
with separation of centers z.
This is a zeroth order homogeneous function, that is,
circleangle(alpha*r, alpha*p, alpha*z) = circleangle(r, p, z).
This version uses masked array operations.
Input:
r one dimensional numpy array
p scalar
z scalar
They should all be non-negative, but there is no other restriction.
Output:
circleangle one dimensional numpy array, same size as r
"""
inside = (r < p-z);
intersect = (r < p+z) & (z < r+p) & numpy.logical_not(inside);
answer = numpy.zeros_like(r);
answer[inside] = numpy.pi;
answer[intersect] = numpy.arccos((numpy.power(r[intersect],2)+z*z-p*p)/(2*z*r[intersect]));
return answer;
def circleanglesorted(r, p, z):
"""circleanglesorted(r, p, z)
Calculate half central angle of the arc of circle of radius r
(which concentrically spans the inside of the star during integration)
that is inside a circle of radius p (planet)
with separation of centers z.
This is a zeroth order homogeneous function, that is,
circleangle(alpha*r, alpha*p, alpha*z) = circleangle(r, p, z).
This version uses a binary search on the sorted r.
Input:
r one dimensional numpy array, must be increasing
p scalar
z scalar
They should all be non-negative, but there is no other restriction.
Output:
circleangle one dimensional numpy array, same size as r
"""
# If the circle arc of radius r is disjoint from the circular disk
# of radius p, then the angle is zero.
answer = numpy.empty_like(r);
if (p > z):
# Planet covers center of star.
a, b = numpy.searchsorted(r, [p-z, p+z], side="right");
answer[:a] = numpy.pi;
answer[a:b] = numpy.arccos((r[a:b]*r[a:b]+z*z-p*p)/(2*z*r[a:b]));
answer[b:] = 0.0;
else:
# Planet does not cover center of star.
a, b = numpy.searchsorted(r, [z-p, z+p], side="right");
answer[:a] = 0.0;
answer[a:b] = numpy.arccos((r[a:b]*r[a:b]+z*z-p*p)/(2*z*r[a:b]));
answer[b:] = 0.0;
return answer;