-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathtoycython.pyx
133 lines (116 loc) · 4.75 KB
/
toycython.pyx
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
# cython: profile=False
#
# 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;
cimport numpy;
ctypedef numpy.float_t FLOAT_t
ctypedef numpy.uint8_t BOOL_t
def circleangleloop(numpy.ndarray[FLOAT_t, ndim=1] r, double p, double 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.
cdef numpy.ndarray[FLOAT_t, ndim=1] answer = numpy.zeros_like(r)
cdef double pminusz = p-z
cdef double pplusz = p+z
cdef double zsquared = z*z
cdef double psquared = p*p
cdef double ri
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(numpy.ndarray[FLOAT_t, ndim=1] r, double p, double 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
"""
cdef numpy.ndarray[BOOL_t, ndim=1, cast=True] inside = (r < p-z)
cdef numpy.ndarray[BOOL_t, ndim=1, cast=True] intersect = (r < p+z) & (z < r+p) & numpy.logical_not(inside)
cdef numpy.ndarray[FLOAT_t, ndim=1] 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(numpy.ndarray[FLOAT_t, ndim=1] r, double p, double 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.
cdef numpy.ndarray[FLOAT_t, ndim=1] 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;