source: sasmodels/explore/guyou.py @ 9149238

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 9149238 was f475f90, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

further guyou code cleanup

  • Property mode set to 100644
File size: 11.4 KB
Line 
1# Copyright 2013-2016 Mike Bostock
2# All rights reserved.
3#
4# Redistribution and use in source and binary forms, with or without modification,
5# are permitted provided that the following conditions are met:
6#
7# * Redistributions of source code must retain the above copyright notice, this
8#   list of conditions and the following disclaimer.
9#
10# * Redistributions in binary form must reproduce the above copyright notice,
11#   this list of conditions and the following disclaimer in the documentation
12#   and/or other materials provided with the distribution.
13#
14# * Neither the name of the author nor the names of contributors may be used to
15#   endorse or promote products derived from this software without specific prior
16#   written permission.
17#
18# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
19# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
20# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
22# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
23# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
24# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
25# ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
26# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
27# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28#
29# https://github.com/d3/d3-geo-projection
30# commit fd2886555e46b35163c7b898d43c7d1bcebbba7c 2016-07-02
31#
32# 2017-11-01 Paul Kienzle
33# * converted to python, using degrees rather than radians
34from __future__ import division, print_function
35
36import numpy as np
37from numpy import sqrt, pi, tan, cos, sin, log, exp, arctan2 as atan2, sign, radians, degrees
38
39_ = """
40# mpmath version of special functions
41import mpmath as mp
42sn, cn, dn = [mp.ellipfun(v) for v in 'sn', 'cn', 'dn']
43
44def ellipticJi(u, v, m):
45    z = u+v*1j
46    return sn(z, m), cn(z, m), dn(z, m)
47
48def ellipticJ(u, m):
49    s, c, d = sn(u, m), cn(u, m), dn(u, m)
50    phi = mp.asin(s)
51    return s, c, d, phi
52
53def ellipticFi(phi, psi, m):
54    z = phi + psi*1j
55    return mp.ellipf(z, m)
56
57def ellipticF(phi, m):
58    return mp.ellipf(phi, m)
59"""
60
61# scipy version of special functions
62from scipy.special import ellipj as ellipticJ, ellipkinc as ellipticF
63from numpy import sinh, sign, arctan as atan
64
65def ellipticJi(u, v, m):
66    scalar = np.isscalar(u) and np.isscalar(v) and np.isscalar(m)
67    u, v, m = np.broadcast_arrays(u, v, m)
68    result = np.empty_like([u,u,u], 'D')
69    real = v==0
70    imag = u==0
71    mixed = ~(real|imag)
72    result[:, real] = _ellipticJi_real(u[real], m[real])
73    result[:, imag] = _ellipticJi_imag(v[imag], m[imag])
74    result[:, mixed] = _ellipticJi(u[mixed], v[mixed], m[mixed])
75    return result[0,:] if scalar else result
76
77def _ellipticJi_real(u, m):
78    sn, cn, dn, phi = ellipticJ(u, m)
79    return sn, cn, dn
80
81def _ellipticJi_imag(v, m):
82    sn, cn, dn, phi = ellipticJ(v, 1-m)
83    return 1j*sn/cn, 1/cn, dn/cn
84
85def _ellipticJi(u, v, m):
86    # Ignoring special cases for now.
87    #   u=0: (1j*b[0]/b[1], 1/b[1], b[2]/b[1])
88    #   v=0: (a[0], a[1], a[2])
89    a = ellipticJ(u, m)
90    b = ellipticJ(v, 1 - m)
91    c = b[1]**2 + m * (a[0] * b[0])**2
92    return [
93        (a[0] * b[2] / c) + 1j*(a[1] * a[2] * b[0] * b[1] / c),
94        (a[1] * b[1] / c) + 1j*(-a[0] * a[2] * b[0] * b[2] / c),
95        (a[2] * b[1] * b[2] / c) + 1j*(-m * a[0] * a[1] * b[0] / c),
96    ]
97
98# calculate F(phi+ipsi|m).
99# see Abramowitz and Stegun, 17.4.11.
100def ellipticFi(phi, psi, m):
101    if np.any(phi == 0):
102        scalar = np.isscalar(phi) and np.isscalar(psi) and np.isscalar(m)
103        phi, psi, m = np.broadcast_arrays(phi, psi, m)
104        result = np.empty_like(phi, 'D')
105        index = (phi == 0)
106        result[index] = ellipticF(atan(sinh(abs(phi[index]))), 1-m[index]) * sign(psi[index])
107        result[~index] = ellipticFi(phi[~index], psi[~index], m[~index])
108        return result.reshape(1)[0] if scalar else result
109
110    r = abs(phi)
111    i = abs(psi)
112    sinhpsi2 = sinh(i)**2
113    cscphi2 = 1 / sin(r)**2
114    cotphi2 = 1 / tan(r)**2
115    b = -(cotphi2 + m * (sinhpsi2 * cscphi2) - 1 + m)
116    c = (m - 1) * cotphi2
117    cotlambda2 = (-b + sqrt(b * b - 4 * c)) / 2
118    re = ellipticF(atan(1 / sqrt(cotlambda2)), m) * sign(phi)
119    im = ellipticF(atan(sqrt(np.maximum(0,(cotlambda2 / cotphi2 - 1) / m))), 1 - m) * sign(psi)
120    return re + 1j*im
121
122sqrt2 = sqrt(2)
123
124# [PAK] renamed k_ => cos_u, k => sin_u, k*k => sinsq_u to avoid k,K confusion
125# cos_u = 0.171572875253809902396622551580603842860656249246103853646...
126# sinsq_u = 0.970562748477140585620264690516376942836062504523376878120...
127# K = 3.165103454447431823666270142140819753058976299237578486994...
128def guyou(lam, phi):
129    # [PAK] wrap into [-pi/2, pi/2] radians
130    x, y = np.asarray(lam), np.asarray(phi)
131    xn, x = divmod(x+90, 180)
132    yn, y = divmod(y+90, 180)
133    xn, lam = xn*180, radians(x-90)
134    yn, phi = yn*180, radians(y-90)
135
136    # Compute constant K
137    cos_u = (sqrt2 - 1) / (sqrt2 + 1)
138    sinsq_u = 1 - cos_u**2
139    K = ellipticF(pi/2, sinsq_u)
140
141    # [PAK] simplify expressions, using the fact that f = -1
142    # Note: exp(f log(x)) => 1/x,  cos(f x) => cos(x), sin(f x) => -sin(x)
143    r = 1/(tan(pi/4 + abs(phi)/2) * sqrt(cos_u))
144    at = atan(r * (cos(lam) - 1j*sin(lam)))
145    t = ellipticFi(at.real, at.imag, sinsq_u)
146    x, y = (-t.imag, sign(phi + (phi==0))*(0.5 * K - t.real))
147
148    # [PAK] convert to degrees, and return to original tile
149    return degrees(x)+xn, degrees(y)+yn
150
151def guyou_invert(x, y):
152    # [PAK] wrap into [-pi/2, pi/2] radians
153    x, y = np.asarray(x), np.asarray(y)
154    xn, x = divmod(x+90, 180)
155    yn, y = divmod(y+90, 180)
156    xn, x = xn*180, radians(x-90)
157    yn, y = yn*180, radians(y-90)
158
159    # compute constant K
160    cos_u = (sqrt2 - 1) / (sqrt2 + 1)
161    sinsq_u = 1 - cos_u**2
162    K = ellipticF(pi/2, sinsq_u)
163
164    # [PAK] simplify expressions, using the fact that f = -1
165    j = ellipticJi(K/2 - y, -x, sinsq_u)
166    tn = j[0]/j[1]  # j[0], j[1] are complex
167    # Note: -atan2(im(x)/re(x)) => angle(x)
168    lam = -np.angle(tn)
169    # Note: exp(0.5/f log(a re(x)^2 + a im(x)^2)) => 1/(sqrt(a) |x|)
170    phi = 2*atan(1/sqrt(cos_u)/abs(tn)) - pi/2
171
172    # [PAK] convert to degrees, and return to original tile
173    return degrees(lam)+xn, degrees(phi)+yn
174
175def plot_grid():
176    import matplotlib.pyplot as plt
177    from numpy import linspace
178    lat_line = linspace(-90, 90, 400)
179    long_line = linspace(-90, 90, 400)
180
181    #scale = 1
182    limit, step, scale = 90, 10, 2
183    plt.subplot(211)
184    for lat in range(-limit, limit+1, step):
185        x, y = guyou(scale*lat_line, scale*lat)
186        plt.plot(x, y, 'g')
187
188    for long in range(-limit, limit+1, step):
189        x, y = guyou(scale*long, scale*long_line)
190        plt.plot(x, y, 'b')
191    #plt.xlabel('longitude')
192    plt.ylabel('latitude')
193
194    plt.subplot(212)
195    for lat in range(-limit, limit+1, step):
196        x, y = guyou_invert(scale*lat_line, scale*lat)
197        plt.plot(x, y, 'g')
198
199    for long in range(-limit, limit+1, step):
200        x, y = guyou_invert(scale*long, scale*long_line)
201        plt.plot(x, y, 'b')
202    plt.xlabel('longitude')
203    plt.ylabel('latitude')
204
205if __name__ == "__main__":
206    plot_grid()
207    import matplotlib.pyplot as plt; plt.show()
208
209
210_ = """
211// Javascript source for elliptic functions
212//
213// Returns [sn, cn, dn](u + iv|m).
214export function ellipticJi(u, v, m) {
215  var a, b, c;
216  if (!u) {
217    b = ellipticJ(v, 1 - m);
218    return [
219      [0, b[0] / b[1]],
220      [1 / b[1], 0],
221      [b[2] / b[1], 0]
222    ];
223  }
224  a = ellipticJ(u, m);
225  if (!v) return [[a[0], 0], [a[1], 0], [a[2], 0]];
226  b = ellipticJ(v, 1 - m);
227  c = b[1] * b[1] + m * a[0] * a[0] * b[0] * b[0];
228  return [
229    [a[0] * b[2] / c, a[1] * a[2] * b[0] * b[1] / c],
230    [a[1] * b[1] / c, -a[0] * a[2] * b[0] * b[2] / c],
231    [a[2] * b[1] * b[2] / c, -m * a[0] * a[1] * b[0] / c]
232  ];
233}
234
235// Returns [sn, cn, dn, ph](u|m).
236export function ellipticJ(u, m) {
237  var ai, b, phi, t, twon;
238  if (m < epsilon) {
239    t = sin(u);
240    b = cos(u);
241    ai = m * (u - t * b) / 4;
242    return [
243      t - ai * b,
244      b + ai * t,
245      1 - m * t * t / 2,
246      u - ai
247    ];
248  }
249  if (m >= 1 - epsilon) {
250    ai = (1 - m) / 4;
251    b = cosh(u);
252    t = tanh(u);
253    phi = 1 / b;
254    twon = b * sinh(u);
255    return [
256      t + ai * (twon - u) / (b * b),
257      phi - ai * t * phi * (twon - u),
258      phi + ai * t * phi * (twon + u),
259      2 * atan(exp(u)) - halfPi + ai * (twon - u) / b
260    ];
261  }
262
263  var a = [1, 0, 0, 0, 0, 0, 0, 0, 0],
264      c = [sqrt(m), 0, 0, 0, 0, 0, 0, 0, 0],
265      i = 0;
266  b = sqrt(1 - m);
267  twon = 1;
268
269  while (abs(c[i] / a[i]) > epsilon && i < 8) {
270    ai = a[i++];
271    c[i] = (ai - b) / 2;
272    a[i] = (ai + b) / 2;
273    b = sqrt(ai * b);
274    twon *= 2;
275  }
276
277  phi = twon * a[i] * u;
278  do {
279    t = c[i] * sin(b = phi) / a[i];
280    phi = (asin(t) + phi) / 2;
281  } while (--i);
282
283  // [PAK] Cephes uses dn = sqrt(1 - m*sin^2 phi) rather than cos(phi)/cos(phi-b)
284  // DLMF says the second version is unstable near x = K.
285  return [sin(phi), t = cos(phi), t / cos(phi - b), phi];
286}
287
288// calculate F(phi+ipsi|m).
289// see Abramowitz and Stegun, 17.4.11.
290export function ellipticFi(phi, psi, m) {
291  var r = abs(phi),
292      i = abs(psi),
293      sinhpsi = sinh(i);
294  if (r) {
295    var cscphi = 1 / sin(r),
296        cotphi2 = 1 / (tan(r) * tan(r)),
297        b = -(cotphi2 + m * (sinhpsi * sinhpsi * cscphi * cscphi) - 1 + m),
298        c = (m - 1) * cotphi2,
299        cotlambda2 = (-b + sqrt(b * b - 4 * c)) / 2;
300    return [
301      ellipticF(atan(1 / sqrt(cotlambda2)), m) * sign(phi),
302      ellipticF(atan(sqrt((cotlambda2 / cotphi2 - 1) / m)), 1 - m) * sign(psi)
303    ];
304  }
305  return [
306    0,
307    ellipticF(atan(sinhpsi), 1 - m) * sign(psi)
308  ];
309}
310
311// Calculate F(phi|m) where m = k^2 = sin(alpha)^2.
312// See Abramowitz and Stegun, 17.6.7.
313export function ellipticF(phi, m) {
314  if (!m) return phi;
315  if (m === 1) return log(tan(phi / 2 + quarterPi));
316  var a = 1,
317      b = sqrt(1 - m),
318      c = sqrt(m);
319  for (var i = 0; abs(c) > epsilon; i++) {
320    if (phi % pi) {
321      var dPhi = atan(b * tan(phi) / a);
322      if (dPhi < 0) dPhi += pi;
323      phi += dPhi + ~~(phi / pi) * pi;
324    } else phi += phi;
325    c = (a + b) / 2;
326    b = sqrt(a * b);
327    c = ((a = c) - b) / 2;
328  }
329  return phi / (pow(2, i) * a);
330
331
332export function guyouRaw(lambda, phi) {
333  var k_ = (sqrt2 - 1) / (sqrt2 + 1),
334      k = sqrt(1 - k_ * k_),
335      K = ellipticF(halfPi, k * k),
336      f = -1,
337      psi = log(tan(pi / 4 + abs(phi) / 2)),
338      r = exp(f * psi) / sqrt(k_),
339      at = guyouComplexAtan(r * cos(f * lambda), r * sin(f * lambda)),
340      t = ellipticFi(at[0], at[1], k * k);
341  return [-t[1], (phi >= 0 ? 1 : -1) * (0.5 * K - t[0])];
342}
343
344function guyouComplexAtan(x, y) {
345  var x2 = x * x,
346      y_1 = y + 1,
347      t = 1 - x2 - y * y;
348  return [
349   0.5 * ((x >= 0 ? halfPi : -halfPi) - atan2(t, 2 * x)),
350    -0.25 * log(t * t + 4 * x2) +0.5 * log(y_1 * y_1 + x2)
351  ];
352}
353
354function guyouComplexDivide(a, b) {
355  var denominator = b[0] * b[0] + b[1] * b[1];
356  return [
357    (a[0] * b[0] + a[1] * b[1]) / denominator,
358    (a[1] * b[0] - a[0] * b[1]) / denominator
359  ];
360}
361
362guyouRaw.invert = function(x, y) {
363  var k_ = (sqrt2 - 1) / (sqrt2 + 1),
364      k = sqrt(1 - k_ * k_),
365      K = ellipticF(halfPi, k * k),
366      f = -1,
367      j = ellipticJi(0.5 * K - y, -x, k * k),
368      tn = guyouComplexDivide(j[0], j[1]),
369      lambda = atan2(tn[1], tn[0]) / f;
370  return [
371    lambda,
372    2 * atan(exp(0.5 / f * log(k_ * tn[0] * tn[0] + k_ * tn[1] * tn[1]))) - halfPi
373  ];
374};
375"""
Note: See TracBrowser for help on using the repository browser.