source: sasmodels/sasmodels/guyou.py @ b297ba9

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

lint

  • Property mode set to 100644
File size: 11.9 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, with degrees rather than radians
34"""
35Convert between latitude-longitude and Guyou map coordinates.
36"""
37
38from __future__ import division, print_function
39
40import numpy as np
41from numpy import sqrt, pi, tan, cos, sin, sign, radians, degrees
42from numpy import sinh, arctan as atan
43
44# scipy version of special functions
45from scipy.special import ellipj as ellipticJ, ellipkinc as ellipticF
46
47_ = """
48# mpmath version of special functions
49import mpmath as mp
50sn, cn, dn = [mp.ellipfun(v) for v in 'sn', 'cn', 'dn']
51
52def ellipticJi(u, v, m):
53    z = u+v*1j
54    return sn(z, m), cn(z, m), dn(z, m)
55
56def ellipticJ(u, m):
57    s, c, d = sn(u, m), cn(u, m), dn(u, m)
58    phi = mp.asin(s)
59    return s, c, d, phi
60
61def ellipticFi(phi, psi, m):
62    z = phi + psi*1j
63    return mp.ellipf(z, m)
64
65def ellipticF(phi, m):
66    return mp.ellipf(phi, m)
67"""
68
69def ellipticJi(u, v, m):
70    """Returns [sn, cn, dn](u + iv|m)."""
71    scalar = np.isscalar(u) and np.isscalar(v) and np.isscalar(m)
72    u, v, m = np.broadcast_arrays(u, v, m)
73    result = np.empty_like([u, u, u], 'D')
74    real = (v == 0)
75    imag = (u == 0)
76    mixed = ~(real|imag)
77    result[:, real] = _ellipticJi_real(u[real], m[real])
78    result[:, imag] = _ellipticJi_imag(v[imag], m[imag])
79    result[:, mixed] = _ellipticJi(u[mixed], v[mixed], m[mixed])
80    return result[0, :] if scalar else result
81
82def _ellipticJi_real(u, m):
83    sn, cn, dn, phi = ellipticJ(u, m)
84    return sn, cn, dn
85
86def _ellipticJi_imag(v, m):
87    sn, cn, dn, phi = ellipticJ(v, 1-m)
88    return 1j*sn/cn, 1/cn, dn/cn
89
90def _ellipticJi(u, v, m):
91    # Ignoring special cases for now.
92    #   u=0: (1j*b[0]/b[1], 1/b[1], b[2]/b[1])
93    #   v=0: (a[0], a[1], a[2])
94    a = ellipticJ(u, m)
95    b = ellipticJ(v, 1 - m)
96    c = b[1]**2 + m * (a[0] * b[0])**2
97    return [
98        (a[0] * b[2] / c) + 1j*(a[1] * a[2] * b[0] * b[1] / c),
99        (a[1] * b[1] / c) + 1j*(-a[0] * a[2] * b[0] * b[2] / c),
100        (a[2] * b[1] * b[2] / c) + 1j*(-m * a[0] * a[1] * b[0] / c),
101    ]
102
103def ellipticFi(phi, psi, m):
104    """Returns F(phi+ipsi|m). See Abramowitz and Stegun, 17.4.11."""
105    if np.any(phi == 0):
106        scalar = np.isscalar(phi) and np.isscalar(psi) and np.isscalar(m)
107        phi, psi, m = np.broadcast_arrays(phi, psi, m)
108        result = np.empty_like(phi, 'D')
109        index = (phi == 0)
110        result[index] = ellipticF(atan(sinh(abs(phi[index]))),
111                                  1-m[index]) * sign(psi[index])
112        result[~index] = ellipticFi(phi[~index], psi[~index], m[~index])
113        return result.reshape(1)[0] if scalar else result
114
115    r = abs(phi)
116    i = abs(psi)
117    sinhpsi2 = sinh(i)**2
118    cscphi2 = 1 / sin(r)**2
119    cotphi2 = 1 / tan(r)**2
120    b = -(cotphi2 + m * (sinhpsi2 * cscphi2) - 1 + m)
121    c = (m - 1) * cotphi2
122    cotlambda2 = (-b + sqrt(b * b - 4 * c)) / 2
123    re = ellipticF(atan(1 / sqrt(cotlambda2)), m) * sign(phi)
124    im = ellipticF(atan(sqrt(np.maximum(0, (cotlambda2 / cotphi2 - 1) / m))),
125                   1 - m) * sign(psi)
126    return re + 1j*im
127
128SQRT2 = sqrt(2)
129
130# [PAK] renamed k_ => cos_u, k => sin_u, k*k => sinsq_u to avoid k,K confusion
131# cos_u = 0.171572875253809902396622551580603842860656249246103853646...
132# sinsq_u = 0.970562748477140585620264690516376942836062504523376878120...
133# K = 3.165103454447431823666270142140819753058976299237578486994...
134def guyou(lam, phi):
135    """Transform from (latitude, longitude) to point (x, y)"""
136    # [PAK] wrap into [-pi/2, pi/2] radians
137    x, y = np.asarray(lam), np.asarray(phi)
138    xn, x = divmod(x+90, 180)
139    yn, y = divmod(y+90, 180)
140    xn, lam = xn*180, radians(x-90)
141    yn, phi = yn*180, radians(y-90)
142
143    # Compute constant K
144    cos_u = (SQRT2 - 1) / (SQRT2 + 1)
145    sinsq_u = 1 - cos_u**2
146    K = ellipticF(pi/2, sinsq_u)
147
148    # [PAK] simplify expressions, using the fact that f = -1
149    # Note: exp(f log(x)) => 1/x,  cos(f x) => cos(x), sin(f x) => -sin(x)
150    r = 1/(tan(pi/4 + abs(phi)/2) * sqrt(cos_u))
151    at = atan(r * (cos(lam) - 1j*sin(lam)))
152    t = ellipticFi(at.real, at.imag, sinsq_u)
153    x, y = (-t.imag, sign(phi + (phi == 0))*(0.5 * K - t.real))
154
155    # [PAK] convert to degrees, and return to original tile
156    return degrees(x)+xn, degrees(y)+yn
157
158def guyou_invert(x, y):
159    """Transform from point (x, y) on plot to (latitude, longitude)"""
160    # [PAK] wrap into [-pi/2, pi/2] radians
161    x, y = np.asarray(x), np.asarray(y)
162    xn, x = divmod(x+90, 180)
163    yn, y = divmod(y+90, 180)
164    xn, x = xn*180, radians(x-90)
165    yn, y = yn*180, radians(y-90)
166
167    # compute constant K
168    cos_u = (SQRT2 - 1) / (SQRT2 + 1)
169    sinsq_u = 1 - cos_u**2
170    K = ellipticF(pi/2, sinsq_u)
171
172    # [PAK] simplify expressions, using the fact that f = -1
173    j = ellipticJi(K/2 - y, -x, sinsq_u)
174    tn = j[0]/j[1]  # j[0], j[1] are complex
175    # Note: -atan2(im(x)/re(x)) => angle(x)
176    lam = -np.angle(tn)
177    # Note: exp(0.5/f log(a re(x)^2 + a im(x)^2)) => 1/(sqrt(a) |x|)
178    phi = 2*atan(1/sqrt(cos_u)/abs(tn)) - pi/2
179
180    # [PAK] convert to degrees, and return to original tile
181    return degrees(lam)+xn, degrees(phi)+yn
182
183def plot_grid():
184    """Plot the latitude-longitude grid for Guyou transform"""
185    import matplotlib.pyplot as plt
186    from numpy import linspace
187    lat_line = linspace(-90, 90, 400)
188    long_line = linspace(-90, 90, 400)
189
190    #scale = 1
191    limit, step, scale = 90, 10, 2
192    plt.subplot(211)
193    for lat in range(-limit, limit+1, step):
194        x, y = guyou(scale*lat_line, scale*lat)
195        plt.plot(x, y, 'g')
196
197    for longitude in range(-limit, limit+1, step):
198        x, y = guyou(scale*longitude, scale*long_line)
199        plt.plot(x, y, 'b')
200    #plt.xlabel('longitude')
201    plt.ylabel('latitude')
202    plt.title('forward transform')
203
204    plt.subplot(212)
205    for lat in range(-limit, limit+1, step):
206        x, y = guyou_invert(scale*lat_line, scale*lat)
207        plt.plot(x, y, 'g')
208
209    for longitude in range(-limit, limit+1, step):
210        x, y = guyou_invert(scale*longitude, scale*long_line)
211        plt.plot(x, y, 'b')
212    plt.xlabel('longitude')
213    plt.ylabel('latitude')
214    plt.title('inverse transform')
215
216def main():
217    """Show the Guyou transformation"""
218    plot_grid()
219    import matplotlib.pyplot as plt
220    plt.show()
221
222if __name__ == "__main__":
223    main()
224
225_ = """
226// Javascript source for elliptic functions
227//
228// Returns [sn, cn, dn](u + iv|m).
229export function ellipticJi(u, v, m) {
230  var a, b, c;
231  if (!u) {
232    b = ellipticJ(v, 1 - m);
233    return [
234      [0, b[0] / b[1]],
235      [1 / b[1], 0],
236      [b[2] / b[1], 0]
237    ];
238  }
239  a = ellipticJ(u, m);
240  if (!v) return [[a[0], 0], [a[1], 0], [a[2], 0]];
241  b = ellipticJ(v, 1 - m);
242  c = b[1] * b[1] + m * a[0] * a[0] * b[0] * b[0];
243  return [
244    [a[0] * b[2] / c, a[1] * a[2] * b[0] * b[1] / c],
245    [a[1] * b[1] / c, -a[0] * a[2] * b[0] * b[2] / c],
246    [a[2] * b[1] * b[2] / c, -m * a[0] * a[1] * b[0] / c]
247  ];
248}
249
250// Returns [sn, cn, dn, ph](u|m).
251export function ellipticJ(u, m) {
252  var ai, b, phi, t, twon;
253  if (m < epsilon) {
254    t = sin(u);
255    b = cos(u);
256    ai = m * (u - t * b) / 4;
257    return [
258      t - ai * b,
259      b + ai * t,
260      1 - m * t * t / 2,
261      u - ai
262    ];
263  }
264  if (m >= 1 - epsilon) {
265    ai = (1 - m) / 4;
266    b = cosh(u);
267    t = tanh(u);
268    phi = 1 / b;
269    twon = b * sinh(u);
270    return [
271      t + ai * (twon - u) / (b * b),
272      phi - ai * t * phi * (twon - u),
273      phi + ai * t * phi * (twon + u),
274      2 * atan(exp(u)) - halfPi + ai * (twon - u) / b
275    ];
276  }
277
278  var a = [1, 0, 0, 0, 0, 0, 0, 0, 0],
279      c = [sqrt(m), 0, 0, 0, 0, 0, 0, 0, 0],
280      i = 0;
281  b = sqrt(1 - m);
282  twon = 1;
283
284  while (abs(c[i] / a[i]) > epsilon && i < 8) {
285    ai = a[i++];
286    c[i] = (ai - b) / 2;
287    a[i] = (ai + b) / 2;
288    b = sqrt(ai * b);
289    twon *= 2;
290  }
291
292  phi = twon * a[i] * u;
293  do {
294    t = c[i] * sin(b = phi) / a[i];
295    phi = (asin(t) + phi) / 2;
296  } while (--i);
297
298  // [PAK] Cephes uses dn = sqrt(1 - m*sin^2 phi) rather than cos(phi)/cos(phi-b)
299  // DLMF says the second version is unstable near x = K.
300  return [sin(phi), t = cos(phi), t / cos(phi - b), phi];
301}
302
303// calculate F(phi+ipsi|m).
304// see Abramowitz and Stegun, 17.4.11.
305export function ellipticFi(phi, psi, m) {
306  var r = abs(phi),
307      i = abs(psi),
308      sinhpsi = sinh(i);
309  if (r) {
310    var cscphi = 1 / sin(r),
311        cotphi2 = 1 / (tan(r) * tan(r)),
312        b = -(cotphi2 + m * (sinhpsi * sinhpsi * cscphi * cscphi) - 1 + m),
313        c = (m - 1) * cotphi2,
314        cotlambda2 = (-b + sqrt(b * b - 4 * c)) / 2;
315    return [
316      ellipticF(atan(1 / sqrt(cotlambda2)), m) * sign(phi),
317      ellipticF(atan(sqrt((cotlambda2 / cotphi2 - 1) / m)), 1 - m) * sign(psi)
318    ];
319  }
320  return [
321    0,
322    ellipticF(atan(sinhpsi), 1 - m) * sign(psi)
323  ];
324}
325
326// Calculate F(phi|m) where m = k^2 = sin(alpha)^2.
327// See Abramowitz and Stegun, 17.6.7.
328export function ellipticF(phi, m) {
329  if (!m) return phi;
330  if (m === 1) return log(tan(phi / 2 + quarterPi));
331  var a = 1,
332      b = sqrt(1 - m),
333      c = sqrt(m);
334  for (var i = 0; abs(c) > epsilon; i++) {
335    if (phi % pi) {
336      var dPhi = atan(b * tan(phi) / a);
337      if (dPhi < 0) dPhi += pi;
338      phi += dPhi + ~~(phi / pi) * pi;
339    } else phi += phi;
340    c = (a + b) / 2;
341    b = sqrt(a * b);
342    c = ((a = c) - b) / 2;
343  }
344  return phi / (pow(2, i) * a);
345
346
347export function guyouRaw(lambda, phi) {
348  var k_ = (sqrt2 - 1) / (sqrt2 + 1),
349      k = sqrt(1 - k_ * k_),
350      K = ellipticF(halfPi, k * k),
351      f = -1,
352      psi = log(tan(pi / 4 + abs(phi) / 2)),
353      r = exp(f * psi) / sqrt(k_),
354      at = guyouComplexAtan(r * cos(f * lambda), r * sin(f * lambda)),
355      t = ellipticFi(at[0], at[1], k * k);
356  return [-t[1], (phi >= 0 ? 1 : -1) * (0.5 * K - t[0])];
357}
358
359function guyouComplexAtan(x, y) {
360  var x2 = x * x,
361      y_1 = y + 1,
362      t = 1 - x2 - y * y;
363  return [
364   0.5 * ((x >= 0 ? halfPi : -halfPi) - atan2(t, 2 * x)),
365    -0.25 * log(t * t + 4 * x2) +0.5 * log(y_1 * y_1 + x2)
366  ];
367}
368
369function guyouComplexDivide(a, b) {
370  var denominator = b[0] * b[0] + b[1] * b[1];
371  return [
372    (a[0] * b[0] + a[1] * b[1]) / denominator,
373    (a[1] * b[0] - a[0] * b[1]) / denominator
374  ];
375}
376
377guyouRaw.invert = function(x, y) {
378  var k_ = (sqrt2 - 1) / (sqrt2 + 1),
379      k = sqrt(1 - k_ * k_),
380      K = ellipticF(halfPi, k * k),
381      f = -1,
382      j = ellipticJi(0.5 * K - y, -x, k * k),
383      tn = guyouComplexDivide(j[0], j[1]),
384      lambda = atan2(tn[1], tn[0]) / f;
385  return [
386    lambda,
387    2 * atan(exp(0.5 / f * log(k_ * tn[0] * tn[0] + k_ * tn[1] * tn[1]))) - halfPi
388  ];
389};
390"""
Note: See TracBrowser for help on using the repository browser.