1# -*- coding: ISO-8859-1 -*- 2# 3# Copyright (C) 2011 Michael Schindler <m-schindler@users.sourceforge.net> 4# 5# This file is part of PyX (https://pyx-project.org/). 6# 7# PyX is free software; you can redistribute it and/or modify 8# it under the terms of the GNU General Public License as published by 9# the Free Software Foundation; either version 2 of the License, or 10# (at your option) any later version. 11# 12# PyX is distributed in the hope that it will be useful, 13# but WITHOUT ANY WARRANTY; without even the implied warranty of 14# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15# GNU General Public License for more details. 16# 17# You should have received a copy of the GNU General Public License 18# along with PyX; if not, write to the Free Software 19# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA 20 21from math import atan2, sin, cos, sqrt, pi 22 23################################################################################ 24# Internal functions of MetaPost 25# 26# This file re-implements some of the functionality of MetaPost 27# (http://tug.org/metapost). MetaPost was developed by John D. Hobby and 28# others. The code of Metapost is in the public domain, which we understand as 29# an implicit permission to reuse the code here (see the comment at 30# http://www.gnu.org/licenses/license-list.html) 31# 32# This file is based on the MetaPost version distributed by TeXLive: 33# svn://tug.org/texlive/trunk/Build/source/texk/web2c/mplibdir revision 22737 # 34# (2011-05-31) 35################################################################################ 36 37# from mplib.h: 38mp_endpoint = 0 39mp_explicit = 1 40mp_given = 2 41mp_curl = 3 42mp_open = 4 43mp_end_cycle = 5 44 45# from mpmath.c: 46unity = 1.0 47two = 2.0 48fraction_half = 0.5 49fraction_one = 1.0 50fraction_three = 3.0 51one_eighty_deg = pi 52three_sixty_deg = 2*pi 53 54def mp_make_choices(knots, epsilon): # <<< 55 """Implements mp_make_choices from metapost (mp.c)""" 56 # 334: If consecutive knots are equal, join them explicitly 57 p = knots 58 while True: 59 q = p.next 60 if p.rtype > mp_explicit and (p.x_pt-q.x_pt)**2 + (p.y_pt-q.y_pt)**2 < epsilon**2: 61 p.rtype = mp_explicit 62 if p.ltype == mp_open: 63 p.ltype = mp_curl 64 p.set_left_curl(unity) 65 q.ltype = mp_explicit 66 if q.rtype == mp_open: 67 q.rtype = mp_curl 68 q.set_right_curl(unity) 69 p.rx_pt = p.x_pt 70 q.lx_pt = p.x_pt 71 p.ry_pt = p.y_pt 72 q.ly_pt = p.y_pt 73 p = q 74 if p is knots: 75 break 76 77 # 335: 78 # If there are no breakpoints, it is necessary to compute the direction angles around an entire cycle. 79 # In this case the mp left type of the first node is temporarily changed to end cycle. 80 # Find the first breakpoint, h, on the path 81 # insert an artificial breakpoint if the path is an unbroken cycle 82 h = knots 83 while True: 84 if h.ltype != mp_open or h.rtype != mp_open: 85 break 86 h = h.next 87 if h is knots: 88 h.ltype = mp_end_cycle 89 break 90 91 p = h 92 while True: 93 # 336: 94 # Fill in the control points between p and the next breakpoint, then advance p to that breakpoint 95 q = p.next 96 if p.rtype >= mp_given: 97 while q.ltype == mp_open and q.rtype == mp_open: 98 q = q.next 99 # the breakpoints are now p and q 100 101 # 346: 102 # Calculate the turning angles psi_k and the distances d(k, k+1) 103 # set n to the length of the path 104 k = 0 105 s = p 106 n = knots.linked_len() 107 delta_x, delta_y, delta, psi = [], [], [], [None] 108 while True: 109 t = s.next 110 assert len(delta_x) == k 111 delta_x.append(t.x_pt - s.x_pt) 112 delta_y.append(t.y_pt - s.y_pt) 113 delta.append(mp_pyth_add(delta_x[k], delta_y[k])) 114 if k > 0: 115 sine = mp_make_fraction(delta_y[k-1], delta[k-1]) 116 cosine = mp_make_fraction(delta_x[k-1], delta[k-1]) 117 psi.append(mp_n_arg( 118 mp_take_fraction(delta_x[k], cosine) + mp_take_fraction(delta_y[k], sine), 119 mp_take_fraction(delta_y[k], cosine) - mp_take_fraction(delta_x[k], sine))) 120 k += 1 121 s = t 122 if s == q: 123 n = k 124 if k >= n and s.ltype != mp_end_cycle: 125 break 126 if k == n: 127 psi.append(0) 128 else: 129 # for closed paths: 130 psi.append(psi[1]) 131 132 # 347: Remove open types at the breakpoints 133 if q.ltype == mp_open: 134 delx_pt = q.rx_pt - q.x_pt 135 dely_pt = q.ry_pt - q.y_pt 136 if delx_pt**2 + dely_pt**2 < epsilon**2: 137 # use curl if the controls are not usable for giving an angle 138 q.ltype = mp_curl 139 q.set_left_curl(unity) 140 else: 141 q.ltype = mp_given 142 q.set_left_given(mp_n_arg(delx_pt, dely_pt)) 143 144 if p.rtype == mp_open and p.ltype == mp_explicit: 145 delx_pt = p.x_pt - p.lx_pt 146 dely_pt = p.y_pt - p.ly_pt 147 if delx_pt**2 + dely_pt**2 < epsilon**2: 148 p.rtype = mp_curl 149 p.set_right_curl(unity) 150 else: 151 p.rtype = mp_given 152 p.set_right_given(mp_n_arg(delx_pt, dely_pt)) 153 154 # call the internal solving routine 155 mp_solve_choices(p, q, n, delta_x, delta_y, delta, psi) 156 157 elif p.rtype == mp_endpoint: 158 # 337: Give reasonable values for the unused control points between p and q 159 p.rx_pt = p.x_pt 160 p.ry_pt = p.y_pt 161 q.lx_pt = q.x_pt 162 q.ly_pt = q.y_pt 163 164 p = q 165 if p is h: 166 break 167# >>> 168def mp_solve_choices(p, q, n, delta_x, delta_y, delta, psi): # <<< 169 """Implements mp_solve_choices form metapost (mp.c)""" 170 uu = [None]*(len(delta)+1) # relations between adjacent angles ("matrix" entries) 171 ww = [None]*len(uu) # additional matrix entries for the cyclic case 172 vv = [None]*len(uu) # angles ("rhs" entries) 173 theta = [None]*len(uu) # solution of the linear system of equations 174 # 348: 175 # the "matrix" is in tridiagonal form, the solution is obtained by Gaussian elimination. 176 # uu and ww are of type "fraction", vv and theta are of type "angle" 177 # k is the current knot number 178 # r, s, t registers for list traversal 179 k = 0 180 s = p 181 r = 0 182 while True: 183 t = s.next 184 if k == 0: # <<< 185 # 354: 186 # Get the linear equations started 187 # or return with the control points in place, if linear equations needn't be solved 188 189 if s.rtype == mp_given: # <<< 190 if t.ltype == mp_given: 191 # 372: Reduce to simple case of two givens and return 192 aa = mp_n_arg(delta_x[0], delta_y[0]) 193 ct, st = mp_n_sin_cos(p.right_given() - aa) 194 cf, sf = mp_n_sin_cos(q.left_given() - aa) 195 mp_set_controls(p, q, delta_x[0], delta_y[0], st, ct, -sf, cf) 196 return 197 else: 198 # 362: 199 vv[0] = s.right_given() - mp_n_arg(delta_x[0], delta_y[0]) 200 vv[0] = reduce_angle(vv[0]) 201 uu[0] = 0 202 ww[0] = 0 203 # >>> 204 elif s.rtype == mp_curl: # <<< 205 if t.ltype == mp_curl: 206 # 373: (mp.pdf) Reduce to simple case of straight line and return 207 p.rtype = mp_explicit 208 q.ltype = mp_explicit 209 lt = abs(q.left_tension()) 210 rt = abs(p.right_tension()) 211 212 ff = mp_make_fraction(unity, 3.0*rt) 213 p.rx_pt = p.x_pt + mp_take_fraction(delta_x[0], ff) 214 p.ry_pt = p.y_pt + mp_take_fraction(delta_y[0], ff) 215 216 ff = mp_make_fraction(unity, 3.0*lt) 217 q.lx_pt = q.x_pt - mp_take_fraction(delta_x[0], ff) 218 q.ly_pt = q.y_pt - mp_take_fraction(delta_y[0], ff) 219 return 220 221 else: # t.ltype != mp_curl 222 # 363: 223 cc = s.right_curl() 224 lt = abs(t.left_tension()) 225 rt = abs(s.right_tension()) 226 uu[0] = mp_curl_ratio(cc, rt, lt) 227 vv[0] = -mp_take_fraction(psi[1], uu[0]) 228 ww[0] = 0 229 # >>> 230 elif s.rtype == mp_open: # <<< 231 uu[0] = 0 232 vv[0] = 0 233 ww[0] = fraction_one 234 # >>> 235 # end of 354 >>> 236 else: # k > 0 <<< 237 238 if s.ltype == mp_end_cycle or s.ltype == mp_open: # <<< 239 # 356: Set up equation to match mock curvatures at z_k; 240 # then finish loop with theta_n adjusted to equal theta_0, if a 241 # cycle has ended 242 243 # 357: Calculate the values 244 # aa = Ak/Bk, bb = Dk/Ck, dd = (3-alpha_{k-1})d(k,k+1), 245 # ee = (3-beta_{k+1})d(k-1,k), cc=(Bk-uk-Ak)/Bk 246 aa = mp_make_fraction(unity, 3.0*abs(r.right_tension()) - unity) 247 dd = mp_take_fraction(delta[k], 248 fraction_three - mp_make_fraction(unity, abs(r.right_tension()))) 249 bb = mp_make_fraction(unity, 3*abs(t.left_tension()) - unity) 250 ee = mp_take_fraction(delta[k-1], 251 fraction_three - mp_make_fraction(unity, abs(t.left_tension()))) 252 cc = fraction_one - mp_take_fraction(uu[k-1], aa) 253 254 # 358: Calculate the ratio ff = Ck/(Ck + Bk - uk-1Ak) 255 dd = mp_take_fraction(dd, cc) 256 lt = abs(s.left_tension()) 257 rt = abs(s.right_tension()) 258 if lt < rt: 259 dd *= (lt/rt)**2 260 elif lt > rt: 261 ee *= (rt/lt)**2 262 ff = mp_make_fraction(ee, ee + dd) 263 264 uu[k] = mp_take_fraction(ff, bb) 265 266 # 359: Calculate the values of vk and wk 267 acc = -mp_take_fraction(psi[k+1], uu[k]) 268 if r.rtype == mp_curl: 269 ww[k] = 0 270 vv[k] = acc - mp_take_fraction(psi[1], fraction_one - ff) 271 else: 272 ff = mp_make_fraction(fraction_one - ff, cc) 273 acc = acc - mp_take_fraction(psi[k], ff) 274 ff = mp_take_fraction(ff, aa) 275 vv[k] = acc - mp_take_fraction(vv[k-1], ff) 276 ww[k] = -mp_take_fraction(ww[k-1], ff) 277 278 if s.ltype == mp_end_cycle: 279 # 360: Adjust theta_n to equal theta_0 and finish loop 280 281 aa = 0 282 bb = fraction_one 283 while True: 284 k -= 1 285 if k == 0: 286 k = n 287 aa = vv[k] - mp_take_fraction(aa, uu[k]) 288 bb = ww[k] - mp_take_fraction(bb, uu[k]) 289 if k == n: 290 break 291 aa = mp_make_fraction(aa, fraction_one - bb) 292 theta[n] = aa 293 vv[0] = aa 294 for k in range(1, n): 295 vv[k] = vv[k] + mp_take_fraction(aa, ww[k]) 296 break 297 # >>> 298 elif s.ltype == mp_curl: # <<< 299 # 364: 300 cc = s.left_curl() 301 lt = abs(s.left_tension()) 302 rt = abs(r.right_tension()) 303 ff = mp_curl_ratio(cc, lt, rt) 304 theta[n] = -mp_make_fraction(mp_take_fraction(vv[n-1], ff), 305 fraction_one - mp_take_fraction(ff, uu[n-1])) 306 break 307 # >>> 308 elif s.ltype == mp_given: # <<< 309 # 361: 310 theta[n] = s.left_given() - mp_n_arg(delta_x[n-1], delta_y[n-1]) 311 theta[n] = reduce_angle(theta[n]) 312 break 313 # >>> 314 315 # end of k == 0, k != 0 >>> 316 317 r = s 318 s = t 319 k += 1 320 321 # 367: 322 # Finish choosing angles and assigning control points 323 for k in range(n-1, -1, -1): 324 theta[k] = vv[k] - mp_take_fraction(theta[k+1], uu[k]) 325 s = p 326 k = 0 327 while True: 328 t = s.next 329 ct, st = mp_n_sin_cos(theta[k]) 330 cf, sf = mp_n_sin_cos(-psi[k+1]-theta[k+1]) 331 mp_set_controls(s, t, delta_x[k], delta_y[k], st, ct, sf, cf) 332 k += 1 333 s = t 334 if k == n: 335 break 336# >>> 337def mp_n_arg(x, y): # <<< 338 return atan2(y, x) 339# >>> 340def mp_n_sin_cos(z): # <<< 341 """Given an integer z that is 2**20 times an angle theta in degrees, the 342 purpose of n sin cos(z) is to set x = r cos theta and y = r sin theta 343 (approximately), for some rather large number r. The maximum of x and y 344 will be between 2**28 and 2**30, so that there will be hardly any loss of 345 accuracy. Then x and y are divided by r.""" 346 # 67: mpmath.pdf 347 return cos(z), sin(z) 348# >>> 349def mp_set_controls(p, q, delta_x, delta_y, st, ct, sf, cf): # <<< 350 """The set controls routine actually puts the control points into a pair of 351 consecutive nodes p and q. Global variables are used to record the values 352 of sin(theta), cos(theta), sin(phi), and cos(phi) needed in this 353 calculation. 354 355 See mp.pdf, item 370""" 356 lt = abs(q.left_tension()) 357 rt = abs(p.right_tension()) 358 rr = mp_velocity(st, ct, sf, cf, rt) 359 ss = mp_velocity(sf, cf, st, ct, lt) 360 if p.right_tension() < 0 or q.left_tension() < 0: 361 # 371: Decrease the velocities, if necessary, to stay inside the bounding triangle 362 # this is the only place where the sign of the tension counts 363 if (st >= 0 and sf >= 0) or (st <= 0 and sf <= 0): 364 sine = mp_take_fraction(abs(st), cf) + mp_take_fraction(abs(sf), ct) # sin(theta+phi) 365 if sine > 0: 366 #sine = mp_take_fraction(sine, fraction_one + unity) # safety factor 367 sine *= 1.00024414062 # safety factor 368 if p.right_tension() < 0: 369 if mp_ab_vs_cd(abs(sf), fraction_one, rr, sine) < 0: 370 rr = mp_make_fraction(abs(sf), sine) 371 if q.left_tension() < 0: 372 if mp_ab_vs_cd(abs(st), fraction_one, ss, sine) < 0: 373 ss = mp_make_fraction(abs(st), sine) 374 375 p.rx_pt = p.x_pt + mp_take_fraction(mp_take_fraction(delta_x, ct) - mp_take_fraction(delta_y, st), rr) 376 p.ry_pt = p.y_pt + mp_take_fraction(mp_take_fraction(delta_y, ct) + mp_take_fraction(delta_x, st), rr) 377 q.lx_pt = q.x_pt - mp_take_fraction(mp_take_fraction(delta_x, cf) + mp_take_fraction(delta_y, sf), ss) 378 q.ly_pt = q.y_pt - mp_take_fraction(mp_take_fraction(delta_y, cf) - mp_take_fraction(delta_x, sf), ss) 379 p.rtype = mp_explicit 380 q.ltype = mp_explicit 381# >>> 382def mp_make_fraction(p, q): # <<< 383 # 17: mpmath.pdf 384 """The make fraction routine produces the fraction equivalent of p/q, given 385 integers p and q; it computes the integer f = floor(2**28 p/q + 1/2), when 386 p and q are positive. 387 388 In machine language this would simply be (2**28*p)div q""" 389 return p/q 390# >>> 391def mp_take_fraction(q, f): # <<< 392 # 20: mpmath.pdf 393 """The dual of make fraction is take fraction, which multiplies a given 394 integer q by a fraction f. When the operands are positive, it computes 395 p = floor(q*f/2**28 + 1/2), a symmetric function of q and f.""" 396 return q*f 397# >>> 398def mp_pyth_add(a, b): # <<< 399 # 44: mpmath.pdf 400 """Pythagorean addition sqrt(a**2 + b**2) is implemented by an elegant 401 iterative scheme due to Cleve Moler and Donald Morrison [IBM Journal of 402 Research and Development 27 (1983), 577-581]. It modifies a and b in such a 403 way that their Pythagorean sum remains invariant, while the smaller 404 argument decreases.""" 405 return sqrt(a*a + b*b) 406# >>> 407def mp_curl_ratio(gamma, a_tension, b_tension): # <<< 408 """The curl ratio subroutine has three arguments, which our previous 409 notation encourages us to call gamma, 1/alpha, and 1/beta. It is a somewhat 410 tedious program to calculate 411 [(3-alpha)alpha^2 gamma + beta^3] / [alpha^3 gamma + (3-beta)beta^2], 412 with the result reduced to 4 if it exceeds 4. (This reduction of curl is 413 necessary only if the curl and tension are both large.) The values of alpha 414 and beta will be at most 4/3. 415 416 See mp.pdf (items 365, 366).""" 417 alpha = 1.0/a_tension 418 beta = 1.0/b_tension 419 return min(4.0, ((3.0-alpha)*alpha**2*gamma + beta**3) / 420 (alpha**3*gamma + (3.0-beta)*beta**2)) 421# >>> 422def mp_ab_vs_cd(a, b, c, d): # <<< 423 """Tests rigorously if ab is greater than, equal to, or less than cd, given 424 integers (a, b, c, d). In most cases a quick decision is reached. The 425 result is +1, 0, or -1 in the three respective cases. 426 See mpmath.pdf (item 33).""" 427 if a*b == c*d: 428 return 0 429 if a*b > c*d: 430 return 1 431 return -1 432# >>> 433def mp_velocity(st, ct, sf, cf, t): # <<< 434 """Metapost's standard velocity subroutine for cubic Bezier curves. 435 See mpmath.pdf (item 30) and mp.pdf (item 339).""" 436 return min(4.0, (2.0 + sqrt(2)*(st-sf/16.0)*(sf-st/16.0)*(ct-cf)) / 437 (1.5*t*(2+(sqrt(5)-1)*ct + (3-sqrt(5))*cf))) 438# >>> 439def reduce_angle(A): # <<< 440 """A macro in mp.c""" 441 if abs(A) > one_eighty_deg: 442 if A > 0: 443 A -= three_sixty_deg 444 else: 445 A += three_sixty_deg 446 return A 447# >>> 448 449# vim:foldmethod=marker:foldmarker=<<<,>>> 450