1 /*
2 # This file is part of the Astrometry.net suite.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
5 %module cutils
6
7 %{
8 #include <math.h>
9
deg2rad(double x)10 static double deg2rad(double x) {
11 return x * (M_PI / 180.0);
12 }
rad2deg(double x)13 static double rad2deg(double x) {
14 return x * (180.0 / M_PI);
15 }
16 %}
17 %inline %{
18
19 #if CONFUSE_EMACS
20 }
21 #endif
22
munu_to_prime(double mu,double nu,double color)23 PyObject* munu_to_prime(double mu, double nu, double color) {
24 return NULL;
25 // return (xp,yp)
26 /*
27 def munu_to_prime(self, mu, nu, color=0):
28 '''
29 mu = a + b * rowm + c * colm
30 nu = d + e * rowm + f * colm
31
32 So
33
34 [rowm; colm] = [b,c; e,f]^-1 * [mu-a; nu-d]
35
36 [b,c; e,f]^1 = [B,C; E,F] in the code below, so
37
38 [rowm; colm] = [B,C; E,F] * [mu-a; nu-d]
39
40 '''
41 a, b, c, d, e, f = self._get_abcdef()
42 #print 'mu,nu', mu, nu, 'a,d', a,d
43 determinant = b * f - c * e
44 #print 'det', determinant
45 B = f / determinant
46 C = -c / determinant
47 E = -e / determinant
48 F = b / determinant
49 #print 'B', B, 'mu-a', mu-a, 'C', C, 'nu-d', nu-d
50 #print 'E', E, 'mu-a', mu-a, 'F', F, 'nu-d', nu-d
51 mua = mu - a
52 # in field 6955, g3, 809 we see a~413
53 #if mua < -180.:
54 # mua += 360.
55 mua = 360. * (mua < -180.)
56 yprime = B * mua + C * (nu - d)
57 xprime = E * mua + F * (nu - d)
58 return xprime,yprime
59 */
60 }
61
prime_to_pixel(double xp,double yp)62 PyObject* prime_to_pixel(double xp, double yp) {
63 return NULL;
64 /*
65 def prime_to_pixel(self, xprime, yprime, color=0):
66 color0 = self._get_ricut()
67 g0, g1, g2, g3 = self._get_drow()
68 h0, h1, h2, h3 = self._get_dcol()
69 px, py, qx, qy = self._get_cscc()
70
71 # #$(%*&^(%$%*& bad documentation.
72 (px,py) = (py,px)
73 (qx,qy) = (qy,qx)
74
75 qx = qx * np.ones_like(xprime)
76 qy = qy * np.ones_like(yprime)
77 #print 'color', color.shape, 'px', px.shape, 'qx', qx.shape
78 xprime -= np.where(color < color0, px * color, qx)
79 yprime -= np.where(color < color0, py * color, qy)
80
81 # Now invert:
82 # yprime = y + g0 + g1 * x + g2 * x**2 + g3 * x**3
83 # xprime = x + h0 + h1 * x + h2 * x**2 + h3 * x**3
84 x = xprime - h0
85 # dumb-ass Newton's method
86 dx = 1.
87 # FIXME -- should just update the ones that aren't zero
88 # FIXME -- should put in some failsafe...
89 while max(np.abs(np.atleast_1d(dx))) > 1e-10:
90 xp = x + h0 + h1 * x + h2 * x**2 + h3 * x**3
91 dxpdx = 1 + h1 + h2 * 2*x + h3 * 3*x**2
92 dx = (xprime - xp) / dxpdx
93 #print 'Max Newton dx', max(abs(dx))
94 x += dx
95 y = yprime - (g0 + g1 * x + g2 * x**2 + g3 * x**3)
96 return (x, y)
97 */
98 }
99
100 /*
101 Convention: ra,dec in degrees
102 node,incl in radians
103 resulting mu,nu in degrees.
104 */
radec_to_munu(double ra,double dec,double node,double incl)105 PyObject* radec_to_munu(double ra, double dec, double node, double incl) {
106 ra = deg2rad(ra);
107 dec = deg2rad(dec);
108 /*
109 double mu = node + atan2(sin(ra - node) * cos(dec) * cos(incl)
110 + sin(dec) * sin(incl),
111 cos(ra - node) * cos(dec));
112 double nu = asin(-sin(ra - node) * cos(dec) * sin(incl)
113 + sin(dec) * cos(incl));
114 */
115 // This version is faster
116 const double sinramnode = sin(ra - node);
117 const double cosdec = cos(dec);
118 const double sindec = sin(dec);
119 const double sinincl = sin(incl);
120 const double cosincl = cos(incl);
121 double mu = node + atan2(sinramnode * cosdec * cosincl + sindec * sinincl,
122 cos(ra - node) * cosdec);
123 double nu = asin(-sinramnode * cosdec * sinincl + sindec * cosincl);
124
125 mu = rad2deg(mu);
126 nu = rad2deg(nu);
127 if (mu < 0)
128 mu += 360.0;
129 if (mu > 360)
130 mu -= 360.0;
131 return PyTuple_Pack(2, PyFloat_FromDouble(mu), PyFloat_FromDouble(nu));
132 }
133
getlistval(PyObject * cachedvals,int i)134 static double getlistval(PyObject* cachedvals, int i) {
135 PyObject* po = PyList_GET_ITEM(cachedvals, i);
136 assert(PyFloat_Check(po));
137 return PyFloat_AsDouble(po);
138 }
139
140 /*
141 Convention: ra,dec in degrees
142 node,incl in radians
143
144 ASSUMES color < riCut
145 ASSUMES color = 0
146 */
radec_to_pixel(double ra,double dec,PyObject * cachedvals)147 PyObject* radec_to_pixel(double ra, double dec, PyObject* cachedvals) {
148 //double node, double incl,
149 //double a, double b, double c,
150 //double d, double e, double f,
151 //double px, double py) {
152
153 assert(PyList_Check(cachedvals));
154 assert(PyList_Size(cachedvals) == 25);
155
156 double node = getlistval(cachedvals, 0);
157 double incl = getlistval(cachedvals, 1);
158 double a = getlistval(cachedvals, 2);
159 double d = getlistval(cachedvals, 5);
160 double B = getlistval(cachedvals, 8);
161 double C = getlistval(cachedvals, 9);
162 double E = getlistval(cachedvals, 10);
163 double F = getlistval(cachedvals, 11);
164 //double px = getlistval(cachedvals, 12);
165 //double py = getlistval(cachedvals, 13);
166 double g0 = getlistval(cachedvals, 16);
167 double g1 = getlistval(cachedvals, 17);
168 double g2 = getlistval(cachedvals, 18);
169 double g3 = getlistval(cachedvals, 19);
170 double h0 = getlistval(cachedvals, 20);
171 double h1 = getlistval(cachedvals, 21);
172 double h2 = getlistval(cachedvals, 22);
173 double h3 = getlistval(cachedvals, 23);
174
175 ra = deg2rad(ra);
176 dec = deg2rad(dec);
177
178 // from previous (radec_to_munu)
179 const double sinramnode = sin(ra - node);
180 const double cosdec = cos(dec);
181 const double sindec = sin(dec);
182 const double sinincl = sin(incl);
183 const double cosincl = cos(incl);
184 double mu = node + atan2(sinramnode * cosdec * cosincl + sindec * sinincl,
185 cos(ra - node) * cosdec);
186 double nu = asin(-sinramnode * cosdec * sinincl + sindec * cosincl);
187
188 mu = rad2deg(mu);
189 nu = rad2deg(nu);
190 if (mu > 360.0)
191 mu -= 360.0;
192 if (mu < 0)
193 mu += 360.0;
194
195 //printf("c: mu,nu %g,%g\n", mu, nu);
196
197 // munu_to_pixel
198 // munu_to_prime,
199 double mua = mu - a;
200 if (mua < -180)
201 mua += 360;
202 double xp = E * mua + F * (nu - d);
203 double yp = B * mua + C * (nu - d);
204
205 //printf("c: xp,yp %g,%g\n", xp,yp);
206 // prime_to_pixel
207 // color terms would go here...
208
209 // Now invert:
210 // yprime = y + g0 + g1 * x + g2 * x**2 + g3 * x**3
211 // xprime = x + h0 + h1 * x + h2 * x**2 + h3 * x**3
212 double x = xp - h0;
213 double dx = 1.;
214 double dxpdx;
215 double xpi;
216 do {
217 //xp = x + h0 + h1 * x + h2 * x*x + h3 * x*x*x;
218 xpi = x + h0 + x * (h1 + x * (h2 + x * h3));
219 dxpdx = 1. + h1 + x * (h2 * 2. + h3 * 3. * x);
220 dx = (xp - xpi) / dxpdx;
221 x += dx;
222 } while (dx > 1e-10);
223 double y = yp - (g0 + x * (g1 + x * (g2 + x * g3)));
224 return PyTuple_Pack(2, PyFloat_FromDouble(x), PyFloat_FromDouble(y));
225 }
226
227
228
229 %}
230
231