|
rpm-build |
8267b0 |
/*
|
|
rpm-build |
8267b0 |
ppedit - A pattern plate editor for Spiro splines.
|
|
rpm-build |
8267b0 |
Copyright (C) 2007 Raph Levien
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
This program is free software; you can redistribute it and/or
|
|
rpm-build |
8267b0 |
modify it under the terms of the GNU General Public License
|
|
rpm-build |
8267b0 |
as published by the Free Software Foundation; either version 3
|
|
rpm-build |
8267b0 |
of the License, or (at your option) any later version.
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
This program is distributed in the hope that it will be useful,
|
|
rpm-build |
8267b0 |
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
rpm-build |
8267b0 |
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
rpm-build |
8267b0 |
GNU General Public License for more details.
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
You should have received a copy of the GNU General Public License
|
|
rpm-build |
8267b0 |
along with this program; if not, write to the Free Software
|
|
rpm-build |
8267b0 |
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
|
|
rpm-build |
8267b0 |
02110-1301, USA.
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
*/
|
|
rpm-build |
8267b0 |
/* C implementation of third-order polynomial spirals. */
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#include <math.h>
|
|
rpm-build |
8267b0 |
#include <stdlib.h>
|
|
rpm-build |
8267b0 |
#include <string.h>
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#include "bezctx_intf.h"
|
|
rpm-build |
8267b0 |
#include "spiro.h"
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#include "spiro-config.h"
|
|
rpm-build |
8267b0 |
#ifdef VERBOSE
|
|
rpm-build |
8267b0 |
#include <stdio.h>
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
typedef struct {
|
|
rpm-build |
8267b0 |
double a[11]; /* band-diagonal matrix */
|
|
rpm-build |
8267b0 |
double al[5]; /* lower part of band-diagonal decomposition */
|
|
rpm-build |
8267b0 |
} bandmat;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#ifndef M_PI
|
|
rpm-build |
8267b0 |
#define M_PI 3.14159265358979323846 /* pi */
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#ifndef N_IS
|
|
rpm-build |
8267b0 |
//int n = 4;
|
|
rpm-build |
8267b0 |
#define N_IS 4
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#ifndef ORDER
|
|
rpm-build |
8267b0 |
#define ORDER 12
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
/* Integrate polynomial spiral curve over range -.5 .. .5. */
|
|
rpm-build |
8267b0 |
static void
|
|
rpm-build |
8267b0 |
integrate_spiro(const double ks[4], double xy[2], int n)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
#if 0
|
|
rpm-build |
8267b0 |
int n = 1024;
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
double th1 = ks[0];
|
|
rpm-build |
8267b0 |
double th2 = .5 * ks[1];
|
|
rpm-build |
8267b0 |
double th3 = (1./6) * ks[2];
|
|
rpm-build |
8267b0 |
double th4 = (1./24) * ks[3];
|
|
rpm-build |
8267b0 |
double x, y;
|
|
rpm-build |
8267b0 |
double ds = 1. / n;
|
|
rpm-build |
8267b0 |
double ds2 = ds * ds;
|
|
rpm-build |
8267b0 |
double ds3 = ds2 * ds;
|
|
rpm-build |
8267b0 |
double k0 = ks[0] * ds;
|
|
rpm-build |
8267b0 |
double k1 = ks[1] * ds;
|
|
rpm-build |
8267b0 |
double k2 = ks[2] * ds;
|
|
rpm-build |
8267b0 |
double k3 = ks[3] * ds;
|
|
rpm-build |
8267b0 |
int i;
|
|
rpm-build |
8267b0 |
double s = .5 * ds - .5;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
x = 0;
|
|
rpm-build |
8267b0 |
y = 0;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
for (i = 0; i < n; i++) {
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#if ORDER > 2
|
|
rpm-build |
8267b0 |
double u, v;
|
|
rpm-build |
8267b0 |
double km0, km1, km2, km3;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
if (n == 1) {
|
|
rpm-build |
8267b0 |
km0 = k0;
|
|
rpm-build |
8267b0 |
km1 = k1 * ds;
|
|
rpm-build |
8267b0 |
km2 = k2 * ds2;
|
|
rpm-build |
8267b0 |
} else {
|
|
rpm-build |
8267b0 |
km0 = (((1./6) * k3 * s + .5 * k2) * s + k1) * s + k0;
|
|
rpm-build |
8267b0 |
km1 = ((.5 * k3 * s + k2) * s + k1) * ds;
|
|
rpm-build |
8267b0 |
km2 = (k3 * s + k2) * ds2;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
km3 = k3 * ds3;
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#if ORDER == 4
|
|
rpm-build |
8267b0 |
double km0_2 = km0 * km0;
|
|
rpm-build |
8267b0 |
u = 24 - km0_2;
|
|
rpm-build |
8267b0 |
v = km1;
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#if ORDER == 6
|
|
rpm-build |
8267b0 |
double km0_2 = km0 * km0;
|
|
rpm-build |
8267b0 |
double km0_4 = km0_2 * km0_2;
|
|
rpm-build |
8267b0 |
u = 24 - km0_2 + (km0_4 - 4 * km0 * km2 - 3 * km1 * km1) * (1./80);
|
|
rpm-build |
8267b0 |
v = km1 + (km3 - 6 * km0_2 * km1) * (1./80);
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#if ORDER == 8
|
|
rpm-build |
8267b0 |
double t1_1 = km0;
|
|
rpm-build |
8267b0 |
double t1_2 = .5 * km1;
|
|
rpm-build |
8267b0 |
double t1_3 = (1./6) * km2;
|
|
rpm-build |
8267b0 |
double t1_4 = (1./24) * km3;
|
|
rpm-build |
8267b0 |
double t2_2 = t1_1 * t1_1;
|
|
rpm-build |
8267b0 |
double t2_3 = 2 * (t1_1 * t1_2);
|
|
rpm-build |
8267b0 |
double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
|
|
rpm-build |
8267b0 |
double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
|
|
rpm-build |
8267b0 |
double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
|
|
rpm-build |
8267b0 |
double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
|
|
rpm-build |
8267b0 |
double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
|
|
rpm-build |
8267b0 |
double t4_4 = t2_2 * t2_2;
|
|
rpm-build |
8267b0 |
double t4_5 = 2 * (t2_2 * t2_3);
|
|
rpm-build |
8267b0 |
double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
|
|
rpm-build |
8267b0 |
double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
|
|
rpm-build |
8267b0 |
double t6_6 = t4_4 * t2_2;
|
|
rpm-build |
8267b0 |
u = 1;
|
|
rpm-build |
8267b0 |
v = 0;
|
|
rpm-build |
8267b0 |
v += (1./12) * t1_2 + (1./80) * t1_4;
|
|
rpm-build |
8267b0 |
u -= (1./24) * t2_2 + (1./160) * t2_4 + (1./896) * t2_6;
|
|
rpm-build |
8267b0 |
v -= (1./480) * t3_4 + (1./2688) * t3_6;
|
|
rpm-build |
8267b0 |
u += (1./1920) * t4_4 + (1./10752) * t4_6;
|
|
rpm-build |
8267b0 |
v += (1./53760) * t5_6;
|
|
rpm-build |
8267b0 |
u -= (1./322560) * t6_6;
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#if ORDER == 10
|
|
rpm-build |
8267b0 |
double t1_1 = km0;
|
|
rpm-build |
8267b0 |
double t1_2 = .5 * km1;
|
|
rpm-build |
8267b0 |
double t1_3 = (1./6) * km2;
|
|
rpm-build |
8267b0 |
double t1_4 = (1./24) * km3;
|
|
rpm-build |
8267b0 |
double t2_2 = t1_1 * t1_1;
|
|
rpm-build |
8267b0 |
double t2_3 = 2 * (t1_1 * t1_2);
|
|
rpm-build |
8267b0 |
double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
|
|
rpm-build |
8267b0 |
double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
|
|
rpm-build |
8267b0 |
double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
|
|
rpm-build |
8267b0 |
double t2_7 = 2 * (t1_3 * t1_4);
|
|
rpm-build |
8267b0 |
double t2_8 = t1_4 * t1_4;
|
|
rpm-build |
8267b0 |
double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
|
|
rpm-build |
8267b0 |
double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
|
|
rpm-build |
8267b0 |
double t3_8 = t2_4 * t1_4 + t2_5 * t1_3 + t2_6 * t1_2 + t2_7 * t1_1;
|
|
rpm-build |
8267b0 |
double t4_4 = t2_2 * t2_2;
|
|
rpm-build |
8267b0 |
double t4_5 = 2 * (t2_2 * t2_3);
|
|
rpm-build |
8267b0 |
double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
|
|
rpm-build |
8267b0 |
double t4_7 = 2 * (t2_2 * t2_5 + t2_3 * t2_4);
|
|
rpm-build |
8267b0 |
double t4_8 = 2 * (t2_2 * t2_6 + t2_3 * t2_5) + t2_4 * t2_4;
|
|
rpm-build |
8267b0 |
double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
|
|
rpm-build |
8267b0 |
double t5_8 = t4_4 * t1_4 + t4_5 * t1_3 + t4_6 * t1_2 + t4_7 * t1_1;
|
|
rpm-build |
8267b0 |
double t6_6 = t4_4 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_7 = t4_4 * t2_3 + t4_5 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2;
|
|
rpm-build |
8267b0 |
double t7_8 = t6_6 * t1_2 + t6_7 * t1_1;
|
|
rpm-build |
8267b0 |
double t8_8 = t6_6 * t2_2;
|
|
rpm-build |
8267b0 |
u = 1;
|
|
rpm-build |
8267b0 |
v = 0;
|
|
rpm-build |
8267b0 |
v += (1./12) * t1_2 + (1./80) * t1_4;
|
|
rpm-build |
8267b0 |
u -= (1./24) * t2_2 + (1./160) * t2_4 + (1./896) * t2_6 + (1./4608) * t2_8;
|
|
rpm-build |
8267b0 |
v -= (1./480) * t3_4 + (1./2688) * t3_6 + (1./13824) * t3_8;
|
|
rpm-build |
8267b0 |
u += (1./1920) * t4_4 + (1./10752) * t4_6 + (1./55296) * t4_8;
|
|
rpm-build |
8267b0 |
v += (1./53760) * t5_6 + (1./276480) * t5_8;
|
|
rpm-build |
8267b0 |
u -= (1./322560) * t6_6 + (1./1.65888e+06) * t6_8;
|
|
rpm-build |
8267b0 |
v -= (1./1.16122e+07) * t7_8;
|
|
rpm-build |
8267b0 |
u += (1./9.28973e+07) * t8_8;
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#if ORDER == 12
|
|
rpm-build |
8267b0 |
double t1_1 = km0;
|
|
rpm-build |
8267b0 |
double t1_2 = .5 * km1;
|
|
rpm-build |
8267b0 |
double t1_3 = (1./6) * km2;
|
|
rpm-build |
8267b0 |
double t1_4 = (1./24) * km3;
|
|
rpm-build |
8267b0 |
double t2_2 = t1_1 * t1_1;
|
|
rpm-build |
8267b0 |
double t2_3 = 2 * (t1_1 * t1_2);
|
|
rpm-build |
8267b0 |
double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
|
|
rpm-build |
8267b0 |
double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
|
|
rpm-build |
8267b0 |
double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
|
|
rpm-build |
8267b0 |
double t2_7 = 2 * (t1_3 * t1_4);
|
|
rpm-build |
8267b0 |
double t2_8 = t1_4 * t1_4;
|
|
rpm-build |
8267b0 |
double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
|
|
rpm-build |
8267b0 |
double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
|
|
rpm-build |
8267b0 |
double t3_8 = t2_4 * t1_4 + t2_5 * t1_3 + t2_6 * t1_2 + t2_7 * t1_1;
|
|
rpm-build |
8267b0 |
double t3_10 = t2_6 * t1_4 + t2_7 * t1_3 + t2_8 * t1_2;
|
|
rpm-build |
8267b0 |
double t4_4 = t2_2 * t2_2;
|
|
rpm-build |
8267b0 |
double t4_5 = 2 * (t2_2 * t2_3);
|
|
rpm-build |
8267b0 |
double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
|
|
rpm-build |
8267b0 |
double t4_7 = 2 * (t2_2 * t2_5 + t2_3 * t2_4);
|
|
rpm-build |
8267b0 |
double t4_8 = 2 * (t2_2 * t2_6 + t2_3 * t2_5) + t2_4 * t2_4;
|
|
rpm-build |
8267b0 |
double t4_9 = 2 * (t2_2 * t2_7 + t2_3 * t2_6 + t2_4 * t2_5);
|
|
rpm-build |
8267b0 |
double t4_10 = 2 * (t2_2 * t2_8 + t2_3 * t2_7 + t2_4 * t2_6) + t2_5 * t2_5;
|
|
rpm-build |
8267b0 |
double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
|
|
rpm-build |
8267b0 |
double t5_8 = t4_4 * t1_4 + t4_5 * t1_3 + t4_6 * t1_2 + t4_7 * t1_1;
|
|
rpm-build |
8267b0 |
double t5_10 = t4_6 * t1_4 + t4_7 * t1_3 + t4_8 * t1_2 + t4_9 * t1_1;
|
|
rpm-build |
8267b0 |
double t6_6 = t4_4 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_7 = t4_4 * t2_3 + t4_5 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_9 = t4_4 * t2_5 + t4_5 * t2_4 + t4_6 * t2_3 + t4_7 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_10 = t4_4 * t2_6 + t4_5 * t2_5 + t4_6 * t2_4 + t4_7 * t2_3 + t4_8 * t2_2;
|
|
rpm-build |
8267b0 |
double t7_8 = t6_6 * t1_2 + t6_7 * t1_1;
|
|
rpm-build |
8267b0 |
double t7_10 = t6_6 * t1_4 + t6_7 * t1_3 + t6_8 * t1_2 + t6_9 * t1_1;
|
|
rpm-build |
8267b0 |
double t8_8 = t6_6 * t2_2;
|
|
rpm-build |
8267b0 |
double t8_9 = t6_6 * t2_3 + t6_7 * t2_2;
|
|
rpm-build |
8267b0 |
double t8_10 = t6_6 * t2_4 + t6_7 * t2_3 + t6_8 * t2_2;
|
|
rpm-build |
8267b0 |
double t9_10 = t8_8 * t1_2 + t8_9 * t1_1;
|
|
rpm-build |
8267b0 |
double t10_10 = t8_8 * t2_2;
|
|
rpm-build |
8267b0 |
u = 1;
|
|
rpm-build |
8267b0 |
v = 0;
|
|
rpm-build |
8267b0 |
v += (1./12) * t1_2 + (1./80) * t1_4;
|
|
rpm-build |
8267b0 |
u -= (1./24) * t2_2 + (1./160) * t2_4 + (1./896) * t2_6 + (1./4608) * t2_8;
|
|
rpm-build |
8267b0 |
v -= (1./480) * t3_4 + (1./2688) * t3_6 + (1./13824) * t3_8 + (1./67584) * t3_10;
|
|
rpm-build |
8267b0 |
u += (1./1920) * t4_4 + (1./10752) * t4_6 + (1./55296) * t4_8 + (1./270336) * t4_10;
|
|
rpm-build |
8267b0 |
v += (1./53760) * t5_6 + (1./276480) * t5_8 + (1./1.35168e+06) * t5_10;
|
|
rpm-build |
8267b0 |
u -= (1./322560) * t6_6 + (1./1.65888e+06) * t6_8 + (1./8.11008e+06) * t6_10;
|
|
rpm-build |
8267b0 |
v -= (1./1.16122e+07) * t7_8 + (1./5.67706e+07) * t7_10;
|
|
rpm-build |
8267b0 |
u += (1./9.28973e+07) * t8_8 + (1./4.54164e+08) * t8_10;
|
|
rpm-build |
8267b0 |
v += (1./4.08748e+09) * t9_10;
|
|
rpm-build |
8267b0 |
u -= (1./4.08748e+10) * t10_10;
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#if ORDER == 14
|
|
rpm-build |
8267b0 |
double t1_1 = km0;
|
|
rpm-build |
8267b0 |
double t1_2 = .5 * km1;
|
|
rpm-build |
8267b0 |
double t1_3 = (1./6) * km2;
|
|
rpm-build |
8267b0 |
double t1_4 = (1./24) * km3;
|
|
rpm-build |
8267b0 |
double t2_2 = t1_1 * t1_1;
|
|
rpm-build |
8267b0 |
double t2_3 = 2 * (t1_1 * t1_2);
|
|
rpm-build |
8267b0 |
double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
|
|
rpm-build |
8267b0 |
double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
|
|
rpm-build |
8267b0 |
double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
|
|
rpm-build |
8267b0 |
double t2_7 = 2 * (t1_3 * t1_4);
|
|
rpm-build |
8267b0 |
double t2_8 = t1_4 * t1_4;
|
|
rpm-build |
8267b0 |
double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
|
|
rpm-build |
8267b0 |
double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
|
|
rpm-build |
8267b0 |
double t3_8 = t2_4 * t1_4 + t2_5 * t1_3 + t2_6 * t1_2 + t2_7 * t1_1;
|
|
rpm-build |
8267b0 |
double t3_10 = t2_6 * t1_4 + t2_7 * t1_3 + t2_8 * t1_2;
|
|
rpm-build |
8267b0 |
double t3_12 = t2_8 * t1_4;
|
|
rpm-build |
8267b0 |
double t4_4 = t2_2 * t2_2;
|
|
rpm-build |
8267b0 |
double t4_5 = 2 * (t2_2 * t2_3);
|
|
rpm-build |
8267b0 |
double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
|
|
rpm-build |
8267b0 |
double t4_7 = 2 * (t2_2 * t2_5 + t2_3 * t2_4);
|
|
rpm-build |
8267b0 |
double t4_8 = 2 * (t2_2 * t2_6 + t2_3 * t2_5) + t2_4 * t2_4;
|
|
rpm-build |
8267b0 |
double t4_9 = 2 * (t2_2 * t2_7 + t2_3 * t2_6 + t2_4 * t2_5);
|
|
rpm-build |
8267b0 |
double t4_10 = 2 * (t2_2 * t2_8 + t2_3 * t2_7 + t2_4 * t2_6) + t2_5 * t2_5;
|
|
rpm-build |
8267b0 |
double t4_11 = 2 * (t2_3 * t2_8 + t2_4 * t2_7 + t2_5 * t2_6);
|
|
rpm-build |
8267b0 |
double t4_12 = 2 * (t2_4 * t2_8 + t2_5 * t2_7) + t2_6 * t2_6;
|
|
rpm-build |
8267b0 |
double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
|
|
rpm-build |
8267b0 |
double t5_8 = t4_4 * t1_4 + t4_5 * t1_3 + t4_6 * t1_2 + t4_7 * t1_1;
|
|
rpm-build |
8267b0 |
double t5_10 = t4_6 * t1_4 + t4_7 * t1_3 + t4_8 * t1_2 + t4_9 * t1_1;
|
|
rpm-build |
8267b0 |
double t5_12 = t4_8 * t1_4 + t4_9 * t1_3 + t4_10 * t1_2 + t4_11 * t1_1;
|
|
rpm-build |
8267b0 |
double t6_6 = t4_4 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_7 = t4_4 * t2_3 + t4_5 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_9 = t4_4 * t2_5 + t4_5 * t2_4 + t4_6 * t2_3 + t4_7 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_10 = t4_4 * t2_6 + t4_5 * t2_5 + t4_6 * t2_4 + t4_7 * t2_3 + t4_8 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_11 = t4_4 * t2_7 + t4_5 * t2_6 + t4_6 * t2_5 + t4_7 * t2_4 + t4_8 * t2_3 + t4_9 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_12 = t4_4 * t2_8 + t4_5 * t2_7 + t4_6 * t2_6 + t4_7 * t2_5 + t4_8 * t2_4 + t4_9 * t2_3 + t4_10 * t2_2;
|
|
rpm-build |
8267b0 |
double t7_8 = t6_6 * t1_2 + t6_7 * t1_1;
|
|
rpm-build |
8267b0 |
double t7_10 = t6_6 * t1_4 + t6_7 * t1_3 + t6_8 * t1_2 + t6_9 * t1_1;
|
|
rpm-build |
8267b0 |
double t7_12 = t6_8 * t1_4 + t6_9 * t1_3 + t6_10 * t1_2 + t6_11 * t1_1;
|
|
rpm-build |
8267b0 |
double t8_8 = t6_6 * t2_2;
|
|
rpm-build |
8267b0 |
double t8_9 = t6_6 * t2_3 + t6_7 * t2_2;
|
|
rpm-build |
8267b0 |
double t8_10 = t6_6 * t2_4 + t6_7 * t2_3 + t6_8 * t2_2;
|
|
rpm-build |
8267b0 |
double t8_11 = t6_6 * t2_5 + t6_7 * t2_4 + t6_8 * t2_3 + t6_9 * t2_2;
|
|
rpm-build |
8267b0 |
double t8_12 = t6_6 * t2_6 + t6_7 * t2_5 + t6_8 * t2_4 + t6_9 * t2_3 + t6_10 * t2_2;
|
|
rpm-build |
8267b0 |
double t9_10 = t8_8 * t1_2 + t8_9 * t1_1;
|
|
rpm-build |
8267b0 |
double t9_12 = t8_8 * t1_4 + t8_9 * t1_3 + t8_10 * t1_2 + t8_11 * t1_1;
|
|
rpm-build |
8267b0 |
double t10_10 = t8_8 * t2_2;
|
|
rpm-build |
8267b0 |
double t10_11 = t8_8 * t2_3 + t8_9 * t2_2;
|
|
rpm-build |
8267b0 |
double t10_12 = t8_8 * t2_4 + t8_9 * t2_3 + t8_10 * t2_2;
|
|
rpm-build |
8267b0 |
double t11_12 = t10_10 * t1_2 + t10_11 * t1_1;
|
|
rpm-build |
8267b0 |
double t12_12 = t10_10 * t2_2;
|
|
rpm-build |
8267b0 |
u = 1;
|
|
rpm-build |
8267b0 |
v = 0;
|
|
rpm-build |
8267b0 |
v += (1./12) * t1_2 + (1./80) * t1_4;
|
|
rpm-build |
8267b0 |
u -= (1./24) * t2_2 + (1./160) * t2_4 + (1./896) * t2_6 + (1./4608) * t2_8;
|
|
rpm-build |
8267b0 |
v -= (1./480) * t3_4 + (1./2688) * t3_6 + (1./13824) * t3_8 + (1./67584) * t3_10 + (1./319488) * t3_12;
|
|
rpm-build |
8267b0 |
u += (1./1920) * t4_4 + (1./10752) * t4_6 + (1./55296) * t4_8 + (1./270336) * t4_10 + (1./1.27795e+06) * t4_12;
|
|
rpm-build |
8267b0 |
v += (1./53760) * t5_6 + (1./276480) * t5_8 + (1./1.35168e+06) * t5_10 + (1./6.38976e+06) * t5_12;
|
|
rpm-build |
8267b0 |
u -= (1./322560) * t6_6 + (1./1.65888e+06) * t6_8 + (1./8.11008e+06) * t6_10 + (1./3.83386e+07) * t6_12;
|
|
rpm-build |
8267b0 |
v -= (1./1.16122e+07) * t7_8 + (1./5.67706e+07) * t7_10 + (1./2.6837e+08) * t7_12;
|
|
rpm-build |
8267b0 |
u += (1./9.28973e+07) * t8_8 + (1./4.54164e+08) * t8_10 + (1./2.14696e+09) * t8_12;
|
|
rpm-build |
8267b0 |
v += (1./4.08748e+09) * t9_10 + (1./1.93226e+10) * t9_12;
|
|
rpm-build |
8267b0 |
u -= (1./4.08748e+10) * t10_10 + (1./1.93226e+11) * t10_12;
|
|
rpm-build |
8267b0 |
v -= (1./2.12549e+12) * t11_12;
|
|
rpm-build |
8267b0 |
u += (1./2.55059e+13) * t12_12;
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#if ORDER == 16
|
|
rpm-build |
8267b0 |
double t1_1 = km0;
|
|
rpm-build |
8267b0 |
double t1_2 = .5 * km1;
|
|
rpm-build |
8267b0 |
double t1_3 = (1./6) * km2;
|
|
rpm-build |
8267b0 |
double t1_4 = (1./24) * km3;
|
|
rpm-build |
8267b0 |
double t2_2 = t1_1 * t1_1;
|
|
rpm-build |
8267b0 |
double t2_3 = 2 * (t1_1 * t1_2);
|
|
rpm-build |
8267b0 |
double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
|
|
rpm-build |
8267b0 |
double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
|
|
rpm-build |
8267b0 |
double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
|
|
rpm-build |
8267b0 |
double t2_7 = 2 * (t1_3 * t1_4);
|
|
rpm-build |
8267b0 |
double t2_8 = t1_4 * t1_4;
|
|
rpm-build |
8267b0 |
double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
|
|
rpm-build |
8267b0 |
double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
|
|
rpm-build |
8267b0 |
double t3_8 = t2_4 * t1_4 + t2_5 * t1_3 + t2_6 * t1_2 + t2_7 * t1_1;
|
|
rpm-build |
8267b0 |
double t3_10 = t2_6 * t1_4 + t2_7 * t1_3 + t2_8 * t1_2;
|
|
rpm-build |
8267b0 |
double t3_12 = t2_8 * t1_4;
|
|
rpm-build |
8267b0 |
double t4_4 = t2_2 * t2_2;
|
|
rpm-build |
8267b0 |
double t4_5 = 2 * (t2_2 * t2_3);
|
|
rpm-build |
8267b0 |
double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
|
|
rpm-build |
8267b0 |
double t4_7 = 2 * (t2_2 * t2_5 + t2_3 * t2_4);
|
|
rpm-build |
8267b0 |
double t4_8 = 2 * (t2_2 * t2_6 + t2_3 * t2_5) + t2_4 * t2_4;
|
|
rpm-build |
8267b0 |
double t4_9 = 2 * (t2_2 * t2_7 + t2_3 * t2_6 + t2_4 * t2_5);
|
|
rpm-build |
8267b0 |
double t4_10 = 2 * (t2_2 * t2_8 + t2_3 * t2_7 + t2_4 * t2_6) + t2_5 * t2_5;
|
|
rpm-build |
8267b0 |
double t4_11 = 2 * (t2_3 * t2_8 + t2_4 * t2_7 + t2_5 * t2_6);
|
|
rpm-build |
8267b0 |
double t4_12 = 2 * (t2_4 * t2_8 + t2_5 * t2_7) + t2_6 * t2_6;
|
|
rpm-build |
8267b0 |
double t4_13 = 2 * (t2_5 * t2_8 + t2_6 * t2_7);
|
|
rpm-build |
8267b0 |
double t4_14 = 2 * (t2_6 * t2_8) + t2_7 * t2_7;
|
|
rpm-build |
8267b0 |
double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
|
|
rpm-build |
8267b0 |
double t5_8 = t4_4 * t1_4 + t4_5 * t1_3 + t4_6 * t1_2 + t4_7 * t1_1;
|
|
rpm-build |
8267b0 |
double t5_10 = t4_6 * t1_4 + t4_7 * t1_3 + t4_8 * t1_2 + t4_9 * t1_1;
|
|
rpm-build |
8267b0 |
double t5_12 = t4_8 * t1_4 + t4_9 * t1_3 + t4_10 * t1_2 + t4_11 * t1_1;
|
|
rpm-build |
8267b0 |
double t5_14 = t4_10 * t1_4 + t4_11 * t1_3 + t4_12 * t1_2 + t4_13 * t1_1;
|
|
rpm-build |
8267b0 |
double t6_6 = t4_4 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_7 = t4_4 * t2_3 + t4_5 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_9 = t4_4 * t2_5 + t4_5 * t2_4 + t4_6 * t2_3 + t4_7 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_10 = t4_4 * t2_6 + t4_5 * t2_5 + t4_6 * t2_4 + t4_7 * t2_3 + t4_8 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_11 = t4_4 * t2_7 + t4_5 * t2_6 + t4_6 * t2_5 + t4_7 * t2_4 + t4_8 * t2_3 + t4_9 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_12 = t4_4 * t2_8 + t4_5 * t2_7 + t4_6 * t2_6 + t4_7 * t2_5 + t4_8 * t2_4 + t4_9 * t2_3 + t4_10 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_13 = t4_5 * t2_8 + t4_6 * t2_7 + t4_7 * t2_6 + t4_8 * t2_5 + t4_9 * t2_4 + t4_10 * t2_3 + t4_11 * t2_2;
|
|
rpm-build |
8267b0 |
double t6_14 = t4_6 * t2_8 + t4_7 * t2_7 + t4_8 * t2_6 + t4_9 * t2_5 + t4_10 * t2_4 + t4_11 * t2_3 + t4_12 * t2_2;
|
|
rpm-build |
8267b0 |
double t7_8 = t6_6 * t1_2 + t6_7 * t1_1;
|
|
rpm-build |
8267b0 |
double t7_10 = t6_6 * t1_4 + t6_7 * t1_3 + t6_8 * t1_2 + t6_9 * t1_1;
|
|
rpm-build |
8267b0 |
double t7_12 = t6_8 * t1_4 + t6_9 * t1_3 + t6_10 * t1_2 + t6_11 * t1_1;
|
|
rpm-build |
8267b0 |
double t7_14 = t6_10 * t1_4 + t6_11 * t1_3 + t6_12 * t1_2 + t6_13 * t1_1;
|
|
rpm-build |
8267b0 |
double t8_8 = t6_6 * t2_2;
|
|
rpm-build |
8267b0 |
double t8_9 = t6_6 * t2_3 + t6_7 * t2_2;
|
|
rpm-build |
8267b0 |
double t8_10 = t6_6 * t2_4 + t6_7 * t2_3 + t6_8 * t2_2;
|
|
rpm-build |
8267b0 |
double t8_11 = t6_6 * t2_5 + t6_7 * t2_4 + t6_8 * t2_3 + t6_9 * t2_2;
|
|
rpm-build |
8267b0 |
double t8_12 = t6_6 * t2_6 + t6_7 * t2_5 + t6_8 * t2_4 + t6_9 * t2_3 + t6_10 * t2_2;
|
|
rpm-build |
8267b0 |
double t8_13 = t6_6 * t2_7 + t6_7 * t2_6 + t6_8 * t2_5 + t6_9 * t2_4 + t6_10 * t2_3 + t6_11 * t2_2;
|
|
rpm-build |
8267b0 |
double t8_14 = t6_6 * t2_8 + t6_7 * t2_7 + t6_8 * t2_6 + t6_9 * t2_5 + t6_10 * t2_4 + t6_11 * t2_3 + t6_12 * t2_2;
|
|
rpm-build |
8267b0 |
double t9_10 = t8_8 * t1_2 + t8_9 * t1_1;
|
|
rpm-build |
8267b0 |
double t9_12 = t8_8 * t1_4 + t8_9 * t1_3 + t8_10 * t1_2 + t8_11 * t1_1;
|
|
rpm-build |
8267b0 |
double t9_14 = t8_10 * t1_4 + t8_11 * t1_3 + t8_12 * t1_2 + t8_13 * t1_1;
|
|
rpm-build |
8267b0 |
double t10_10 = t8_8 * t2_2;
|
|
rpm-build |
8267b0 |
double t10_11 = t8_8 * t2_3 + t8_9 * t2_2;
|
|
rpm-build |
8267b0 |
double t10_12 = t8_8 * t2_4 + t8_9 * t2_3 + t8_10 * t2_2;
|
|
rpm-build |
8267b0 |
double t10_13 = t8_8 * t2_5 + t8_9 * t2_4 + t8_10 * t2_3 + t8_11 * t2_2;
|
|
rpm-build |
8267b0 |
double t10_14 = t8_8 * t2_6 + t8_9 * t2_5 + t8_10 * t2_4 + t8_11 * t2_3 + t8_12 * t2_2;
|
|
rpm-build |
8267b0 |
double t11_12 = t10_10 * t1_2 + t10_11 * t1_1;
|
|
rpm-build |
8267b0 |
double t11_14 = t10_10 * t1_4 + t10_11 * t1_3 + t10_12 * t1_2 + t10_13 * t1_1;
|
|
rpm-build |
8267b0 |
double t12_12 = t10_10 * t2_2;
|
|
rpm-build |
8267b0 |
double t12_13 = t10_10 * t2_3 + t10_11 * t2_2;
|
|
rpm-build |
8267b0 |
double t12_14 = t10_10 * t2_4 + t10_11 * t2_3 + t10_12 * t2_2;
|
|
rpm-build |
8267b0 |
double t13_14 = t12_12 * t1_2 + t12_13 * t1_1;
|
|
rpm-build |
8267b0 |
double t14_14 = t12_12 * t2_2;
|
|
rpm-build |
8267b0 |
u = 1;
|
|
rpm-build |
8267b0 |
u -= 1./24 * t2_2 + 1./160 * t2_4 + 1./896 * t2_6 + 1./4608 * t2_8;
|
|
rpm-build |
8267b0 |
u += 1./1920 * t4_4 + 1./10752 * t4_6 + 1./55296 * t4_8 + 1./270336 * t4_10 + 1./1277952 * t4_12 + 1./5898240 * t4_14;
|
|
rpm-build |
8267b0 |
u -= 1./322560 * t6_6 + 1./1658880 * t6_8 + 1./8110080 * t6_10 + 1./38338560 * t6_12 + 1./176947200 * t6_14;
|
|
rpm-build |
8267b0 |
u += 1./92897280 * t8_8 + 1./454164480 * t8_10 + 4.6577500191e-10 * t8_12 + 1.0091791708e-10 * t8_14;
|
|
rpm-build |
8267b0 |
u -= 2.4464949595e-11 * t10_10 + 5.1752777990e-12 * t10_12 + 1.1213101898e-12 * t10_14;
|
|
rpm-build |
8267b0 |
u += 3.9206649992e-14 * t12_12 + 8.4947741650e-15 * t12_14;
|
|
rpm-build |
8267b0 |
u -= 4.6674583324e-17 * t14_14;
|
|
rpm-build |
8267b0 |
v = 0;
|
|
rpm-build |
8267b0 |
v += 1./12 * t1_2 + 1./80 * t1_4;
|
|
rpm-build |
8267b0 |
v -= 1./480 * t3_4 + 1./2688 * t3_6 + 1./13824 * t3_8 + 1./67584 * t3_10 + 1./319488 * t3_12;
|
|
rpm-build |
8267b0 |
v += 1./53760 * t5_6 + 1./276480 * t5_8 + 1./1351680 * t5_10 + 1./6389760 * t5_12 + 1./29491200 * t5_14;
|
|
rpm-build |
8267b0 |
v -= 1./11612160 * t7_8 + 1./56770560 * t7_10 + 1./268369920 * t7_12 + 8.0734333664e-10 * t7_14;
|
|
rpm-build |
8267b0 |
v += 2.4464949595e-10 * t9_10 + 5.1752777990e-11 * t9_12 + 1.1213101898e-11 * t9_14;
|
|
rpm-build |
8267b0 |
v -= 4.7047979991e-13 * t11_12 + 1.0193728998e-13 * t11_14;
|
|
rpm-build |
8267b0 |
v += 6.5344416654e-16 * t13_14;
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
if (n == 1) {
|
|
rpm-build |
8267b0 |
#if ORDER == 2
|
|
rpm-build |
8267b0 |
x = 1;
|
|
rpm-build |
8267b0 |
y = 0;
|
|
rpm-build |
8267b0 |
#else
|
|
rpm-build |
8267b0 |
x = u;
|
|
rpm-build |
8267b0 |
y = v;
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
} else {
|
|
rpm-build |
8267b0 |
double th = (((th4 * s + th3) * s + th2) * s + th1) * s;
|
|
rpm-build |
8267b0 |
double cth = cos(th);
|
|
rpm-build |
8267b0 |
double sth = sin(th);
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#if ORDER == 2
|
|
rpm-build |
8267b0 |
x += cth;
|
|
rpm-build |
8267b0 |
y += sth;
|
|
rpm-build |
8267b0 |
#else
|
|
rpm-build |
8267b0 |
x += cth * u - sth * v;
|
|
rpm-build |
8267b0 |
y += cth * v + sth * u;
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
s += ds;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#if ORDER == 4 || ORDER == 6
|
|
rpm-build |
8267b0 |
xy[0] = x * (1./24 * ds);
|
|
rpm-build |
8267b0 |
xy[1] = y * (1./24 * ds);
|
|
rpm-build |
8267b0 |
#else
|
|
rpm-build |
8267b0 |
xy[0] = x * ds;
|
|
rpm-build |
8267b0 |
xy[1] = y * ds;
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
static double
|
|
rpm-build |
8267b0 |
compute_ends(const double ks[4], double ends[2][4], double seg_ch)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
double xy[2];
|
|
rpm-build |
8267b0 |
double ch, th;
|
|
rpm-build |
8267b0 |
double l, l2, l3;
|
|
rpm-build |
8267b0 |
double th_even, th_odd;
|
|
rpm-build |
8267b0 |
double k0_even, k0_odd;
|
|
rpm-build |
8267b0 |
double k1_even, k1_odd;
|
|
rpm-build |
8267b0 |
double k2_even, k2_odd;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
integrate_spiro(ks, xy, N_IS);
|
|
rpm-build |
8267b0 |
ch = hypot(xy[0], xy[1]);
|
|
rpm-build |
8267b0 |
th = atan2(xy[1], xy[0]);
|
|
rpm-build |
8267b0 |
l = ch / seg_ch;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
th_even = .5 * ks[0] + (1./48) * ks[2];
|
|
rpm-build |
8267b0 |
th_odd = .125 * ks[1] + (1./384) * ks[3] - th;
|
|
rpm-build |
8267b0 |
ends[0][0] = th_even - th_odd;
|
|
rpm-build |
8267b0 |
ends[1][0] = th_even + th_odd;
|
|
rpm-build |
8267b0 |
k0_even = l * (ks[0] + .125 * ks[2]);
|
|
rpm-build |
8267b0 |
k0_odd = l * (.5 * ks[1] + (1./48) * ks[3]);
|
|
rpm-build |
8267b0 |
ends[0][1] = k0_even - k0_odd;
|
|
rpm-build |
8267b0 |
ends[1][1] = k0_even + k0_odd;
|
|
rpm-build |
8267b0 |
l2 = l * l;
|
|
rpm-build |
8267b0 |
k1_even = l2 * (ks[1] + .125 * ks[3]);
|
|
rpm-build |
8267b0 |
k1_odd = l2 * .5 * ks[2];
|
|
rpm-build |
8267b0 |
ends[0][2] = k1_even - k1_odd;
|
|
rpm-build |
8267b0 |
ends[1][2] = k1_even + k1_odd;
|
|
rpm-build |
8267b0 |
l3 = l2 * l;
|
|
rpm-build |
8267b0 |
k2_even = l3 * ks[2];
|
|
rpm-build |
8267b0 |
k2_odd = l3 * .5 * ks[3];
|
|
rpm-build |
8267b0 |
ends[0][3] = k2_even - k2_odd;
|
|
rpm-build |
8267b0 |
ends[1][3] = k2_even + k2_odd;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
return l;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
static void
|
|
rpm-build |
8267b0 |
compute_pderivs(const spiro_seg *s, double ends[2][4], double derivs[4][2][4],
|
|
rpm-build |
8267b0 |
int jinc)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
double recip_d = 2e6;
|
|
rpm-build |
8267b0 |
double delta = 1./ recip_d;
|
|
rpm-build |
8267b0 |
double try_ks[4];
|
|
rpm-build |
8267b0 |
double try_ends[2][4];
|
|
rpm-build |
8267b0 |
int i, j, k;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
compute_ends(s->ks, ends, s->seg_ch);
|
|
rpm-build |
8267b0 |
for (i = 0; i < jinc; i++) {
|
|
rpm-build |
8267b0 |
for (j = 0; j < 4; j++)
|
|
rpm-build |
8267b0 |
try_ks[j] = s->ks[j];
|
|
rpm-build |
8267b0 |
try_ks[i] += delta;
|
|
rpm-build |
8267b0 |
compute_ends(try_ks, try_ends, s->seg_ch);
|
|
rpm-build |
8267b0 |
for (k = 0; k < 2; k++)
|
|
rpm-build |
8267b0 |
for (j = 0; j < 4; j++)
|
|
rpm-build |
8267b0 |
derivs[j][k][i] = recip_d * (try_ends[k][j] - ends[k][j]);
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
static double
|
|
rpm-build |
8267b0 |
mod_2pi(double th)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
double u = th / (2 * M_PI);
|
|
rpm-build |
8267b0 |
return 2 * M_PI * (u - floor(u + 0.5));
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
static spiro_seg *
|
|
rpm-build |
8267b0 |
setup_path(const spiro_cp *src, int n)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
int n_seg = src[0].ty == '{' ? n - 1 : n;
|
|
rpm-build |
8267b0 |
spiro_seg *r = (spiro_seg *)malloc((n_seg + 1) * sizeof(spiro_seg));
|
|
rpm-build |
8267b0 |
if ( r==NULL ) return 0;
|
|
rpm-build |
8267b0 |
int i;
|
|
rpm-build |
8267b0 |
int ilast;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
for (i = 0; i < n_seg; i++) {
|
|
rpm-build |
8267b0 |
r[i].x = src[i].x;
|
|
rpm-build |
8267b0 |
r[i].y = src[i].y;
|
|
rpm-build |
8267b0 |
r[i].ty = src[i].ty;
|
|
rpm-build |
8267b0 |
r[i].ks[0] = 0.;
|
|
rpm-build |
8267b0 |
r[i].ks[1] = 0.;
|
|
rpm-build |
8267b0 |
r[i].ks[2] = 0.;
|
|
rpm-build |
8267b0 |
r[i].ks[3] = 0.;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
r[n_seg].x = src[n_seg % n].x;
|
|
rpm-build |
8267b0 |
r[n_seg].y = src[n_seg % n].y;
|
|
rpm-build |
8267b0 |
r[n_seg].ty = src[n_seg % n].ty;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#ifdef CHECK_INPUT_FINITENESS
|
|
rpm-build |
8267b0 |
/* Verify that input values are within realistic limits */
|
|
rpm-build |
8267b0 |
for (i = 0; i < n; i++) {
|
|
rpm-build |
8267b0 |
if (IS_FINITE(r[i].x)==0 || IS_FINITE(r[i].y)==0) {
|
|
rpm-build |
8267b0 |
#ifdef VERBOSE
|
|
rpm-build |
8267b0 |
fprintf(stderr, "ERROR: LibSpiro: #%d={'%c',%g,%g} is not finite.\n", \
|
|
rpm-build |
8267b0 |
i, src[i].ty, r[i].x, r[i].y);
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
free(r);
|
|
rpm-build |
8267b0 |
return 0;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
for (i = 0; i < n_seg; i++) {
|
|
rpm-build |
8267b0 |
double dx = r[i + 1].x - r[i].x;
|
|
rpm-build |
8267b0 |
double dy = r[i + 1].y - r[i].y;
|
|
rpm-build |
8267b0 |
#ifndef CHECK_INPUT_FINITENESS
|
|
rpm-build |
8267b0 |
r[i].seg_ch = hypot(dx, dy);
|
|
rpm-build |
8267b0 |
#else
|
|
rpm-build |
8267b0 |
if (IS_FINITE(dx)==0 || IS_FINITE(dy)==0 || \
|
|
rpm-build |
8267b0 |
IS_FINITE((r[i].seg_ch = hypot(dx, dy)))==0) {
|
|
rpm-build |
8267b0 |
#ifdef VERBOSE
|
|
rpm-build |
8267b0 |
fprintf(stderr, "ERROR: LibSpiro: #%d={'%c',%g,%g} hypot error.\n", \
|
|
rpm-build |
8267b0 |
i, src[i].ty, r[i].x, r[i].y);
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
free(r);
|
|
rpm-build |
8267b0 |
return 0;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
r[i].seg_th = atan2(dy, dx);
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
ilast = n_seg - 1;
|
|
rpm-build |
8267b0 |
for (i = 0; i < n_seg; i++) {
|
|
rpm-build |
8267b0 |
if (r[i].ty == '{' || r[i].ty == '}' || r[i].ty == 'v')
|
|
rpm-build |
8267b0 |
r[i].bend_th = 0.;
|
|
rpm-build |
8267b0 |
else
|
|
rpm-build |
8267b0 |
r[i].bend_th = mod_2pi(r[i].seg_th - r[ilast].seg_th);
|
|
rpm-build |
8267b0 |
ilast = i;
|
|
rpm-build |
8267b0 |
#ifdef VERBOSE
|
|
rpm-build |
8267b0 |
printf("input #%d={'%c',%g,%g}, hypot=%g, atan2=%g, bend_th=%g\n", \
|
|
rpm-build |
8267b0 |
i, src[i].ty, r[i].x, r[i].y, r[i]. seg_th, r[i].seg_th, r[i].bend_th);
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
#ifdef VERBOSE
|
|
rpm-build |
8267b0 |
if (n_seg < n)
|
|
rpm-build |
8267b0 |
printf("input #%d={'%c',%g,%g}\n", i, src[i].ty, r[i].x, r[i].y);
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
return r;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
static void
|
|
rpm-build |
8267b0 |
bandec11(bandmat *m, int *perm, int n)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
int i, j, k;
|
|
rpm-build |
8267b0 |
int l;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
/* pack top triangle to the left. */
|
|
rpm-build |
8267b0 |
for (i = 0; i < 5; i++) {
|
|
rpm-build |
8267b0 |
for (j = 0; j < i + 6; j++)
|
|
rpm-build |
8267b0 |
m[i].a[j] = m[i].a[j + 5 - i];
|
|
rpm-build |
8267b0 |
for (; j < 11; j++)
|
|
rpm-build |
8267b0 |
m[i].a[j] = 0.;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
l = 5;
|
|
rpm-build |
8267b0 |
for (k = 0; k < n; k++) {
|
|
rpm-build |
8267b0 |
int pivot = k;
|
|
rpm-build |
8267b0 |
double pivot_val = m[k].a[0];
|
|
rpm-build |
8267b0 |
double pivot_scale;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
l = l < n ? l + 1 : n;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
for (j = k + 1; j < l; j++)
|
|
rpm-build |
8267b0 |
if (fabs(m[j].a[0]) > fabs(pivot_val)) {
|
|
rpm-build |
8267b0 |
pivot_val = m[j].a[0];
|
|
rpm-build |
8267b0 |
pivot = j;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
perm[k] = pivot;
|
|
rpm-build |
8267b0 |
if (pivot != k) {
|
|
rpm-build |
8267b0 |
for (j = 0; j < 11; j++) {
|
|
rpm-build |
8267b0 |
double tmp = m[k].a[j];
|
|
rpm-build |
8267b0 |
m[k].a[j] = m[pivot].a[j];
|
|
rpm-build |
8267b0 |
m[pivot].a[j] = tmp;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
if (fabs(pivot_val) < 1e-12) pivot_val = 1e-12;
|
|
rpm-build |
8267b0 |
pivot_scale = 1. / pivot_val;
|
|
rpm-build |
8267b0 |
for (i = k + 1; i < l; i++) {
|
|
rpm-build |
8267b0 |
double x = m[i].a[0] * pivot_scale;
|
|
rpm-build |
8267b0 |
m[k].al[i - k - 1] = x;
|
|
rpm-build |
8267b0 |
for (j = 1; j < 11; j++)
|
|
rpm-build |
8267b0 |
m[i].a[j - 1] = m[i].a[j] - x * m[k].a[j];
|
|
rpm-build |
8267b0 |
m[i].a[10] = 0.;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
static void
|
|
rpm-build |
8267b0 |
banbks11(const bandmat *m, const int *perm, double *v, int n)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
int i, k, l;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
/* forward substitution */
|
|
rpm-build |
8267b0 |
l = 5;
|
|
rpm-build |
8267b0 |
for (k = 0; k < n; k++) {
|
|
rpm-build |
8267b0 |
i = perm[k];
|
|
rpm-build |
8267b0 |
if (i != k) {
|
|
rpm-build |
8267b0 |
double tmp = v[k];
|
|
rpm-build |
8267b0 |
v[k] = v[i];
|
|
rpm-build |
8267b0 |
v[i] = tmp;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
if (l < n) l++;
|
|
rpm-build |
8267b0 |
for (i = k + 1; i < l; i++)
|
|
rpm-build |
8267b0 |
v[i] -= m[k].al[i - k - 1] * v[k];
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
/* back substitution */
|
|
rpm-build |
8267b0 |
l = 1;
|
|
rpm-build |
8267b0 |
for (i = n - 1; i >= 0; i--) {
|
|
rpm-build |
8267b0 |
double x = v[i];
|
|
rpm-build |
8267b0 |
for (k = 1; k < l; k++)
|
|
rpm-build |
8267b0 |
x -= m[i].a[k] * v[k + i];
|
|
rpm-build |
8267b0 |
v[i] = x / m[i].a[0];
|
|
rpm-build |
8267b0 |
if (l < 11) l++;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
static int compute_jinc(char ty0, char ty1)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
if (ty0 == 'o' || ty1 == 'o' ||
|
|
rpm-build |
8267b0 |
ty0 == ']' || ty1 == '[')
|
|
rpm-build |
8267b0 |
return 4;
|
|
rpm-build |
8267b0 |
else if (ty0 == 'c' && ty1 == 'c')
|
|
rpm-build |
8267b0 |
return 2;
|
|
rpm-build |
8267b0 |
else if (((ty0 == '{' || ty0 == 'v' || ty0 == '[') && ty1 == 'c') ||
|
|
rpm-build |
8267b0 |
(ty0 == 'c' && (ty1 == '}' || ty1 == 'v' || ty1 == ']')))
|
|
rpm-build |
8267b0 |
return 1;
|
|
rpm-build |
8267b0 |
else
|
|
rpm-build |
8267b0 |
return 0;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
static int count_vec(const spiro_seg *s, int nseg)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
int i;
|
|
rpm-build |
8267b0 |
int n = 0;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
for (i = 0; i < nseg; i++)
|
|
rpm-build |
8267b0 |
n += compute_jinc(s[i].ty, s[i + 1].ty);
|
|
rpm-build |
8267b0 |
return n;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
static void
|
|
rpm-build |
8267b0 |
add_mat_line(bandmat *m, double *v,
|
|
rpm-build |
8267b0 |
double derivs[4], double x, double y, int j, int jj, int jinc,
|
|
rpm-build |
8267b0 |
int nmat)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
int k;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
if (jj >= 0) {
|
|
rpm-build |
8267b0 |
int joff = (j + 5 - jj + nmat) % nmat;
|
|
rpm-build |
8267b0 |
if (nmat < 6) {
|
|
rpm-build |
8267b0 |
joff = j + 5 - jj;
|
|
rpm-build |
8267b0 |
} else if (nmat == 6) {
|
|
rpm-build |
8267b0 |
joff = 2 + (j + 3 - jj + nmat) % nmat;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
#ifdef VERBOSE
|
|
rpm-build |
8267b0 |
printf("add_mat_line j=%d jj=%d jinc=%d nmat=%d joff=%d\n", j, jj, jinc, nmat, joff);
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
v[jj] += x;
|
|
rpm-build |
8267b0 |
for (k = 0; k < jinc; k++)
|
|
rpm-build |
8267b0 |
m[jj].a[joff + k] += y * derivs[k];
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
static double
|
|
rpm-build |
8267b0 |
spiro_iter(spiro_seg *s, bandmat *m, int *perm, double *v, int n, int nmat)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
int cyclic = s[0].ty != '{' && s[0].ty != 'v';
|
|
rpm-build |
8267b0 |
int i, j, jj;
|
|
rpm-build |
8267b0 |
double norm;
|
|
rpm-build |
8267b0 |
int n_invert;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
for (i = 0; i < nmat; i++) {
|
|
rpm-build |
8267b0 |
v[i] = 0.;
|
|
rpm-build |
8267b0 |
for (j = 0; j < 11; j++)
|
|
rpm-build |
8267b0 |
m[i].a[j] = 0.;
|
|
rpm-build |
8267b0 |
for (j = 0; j < 5; j++)
|
|
rpm-build |
8267b0 |
m[i].al[j] = 0.;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
j = 0;
|
|
rpm-build |
8267b0 |
if (s[0].ty == 'o')
|
|
rpm-build |
8267b0 |
jj = nmat - 2;
|
|
rpm-build |
8267b0 |
else if (s[0].ty == 'c')
|
|
rpm-build |
8267b0 |
jj = nmat - 1;
|
|
rpm-build |
8267b0 |
else
|
|
rpm-build |
8267b0 |
jj = 0;
|
|
rpm-build |
8267b0 |
for (i = 0; i < n; i++) {
|
|
rpm-build |
8267b0 |
char ty0 = s[i].ty;
|
|
rpm-build |
8267b0 |
char ty1 = s[i + 1].ty;
|
|
rpm-build |
8267b0 |
int jinc = compute_jinc(ty0, ty1);
|
|
rpm-build |
8267b0 |
double th = s[i].bend_th;
|
|
rpm-build |
8267b0 |
double ends[2][4];
|
|
rpm-build |
8267b0 |
double derivs[4][2][4];
|
|
rpm-build |
8267b0 |
int jthl = -1, jk0l = -1, jk1l = -1, jk2l = -1;
|
|
rpm-build |
8267b0 |
int jthr = -1, jk0r = -1, jk1r = -1, jk2r = -1;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
compute_pderivs(&s[i], ends, derivs, jinc);
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
/* constraints crossing left */
|
|
rpm-build |
8267b0 |
if (ty0 == 'o' || ty0 == 'c' || ty0 == '[' || ty0 == ']') {
|
|
rpm-build |
8267b0 |
jthl = jj++;
|
|
rpm-build |
8267b0 |
jj %= nmat;
|
|
rpm-build |
8267b0 |
jk0l = jj++;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
if (ty0 == 'o') {
|
|
rpm-build |
8267b0 |
jj %= nmat;
|
|
rpm-build |
8267b0 |
jk1l = jj++;
|
|
rpm-build |
8267b0 |
jk2l = jj++;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
/* constraints on left */
|
|
rpm-build |
8267b0 |
if ((ty0 == '[' || ty0 == 'v' || ty0 == '{' || ty0 == 'c') &&
|
|
rpm-build |
8267b0 |
jinc == 4) {
|
|
rpm-build |
8267b0 |
if (ty0 != 'c')
|
|
rpm-build |
8267b0 |
jk1l = jj++;
|
|
rpm-build |
8267b0 |
jk2l = jj++;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
/* constraints on right */
|
|
rpm-build |
8267b0 |
if ((ty1 == ']' || ty1 == 'v' || ty1 == '}' || ty1 == 'c') &&
|
|
rpm-build |
8267b0 |
jinc == 4) {
|
|
rpm-build |
8267b0 |
if (ty1 != 'c')
|
|
rpm-build |
8267b0 |
jk1r = jj++;
|
|
rpm-build |
8267b0 |
jk2r = jj++;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
/* constraints crossing right */
|
|
rpm-build |
8267b0 |
if (ty1 == 'o' || ty1 == 'c' || ty1 == '[' || ty1 == ']') {
|
|
rpm-build |
8267b0 |
jthr = jj;
|
|
rpm-build |
8267b0 |
jk0r = (jj + 1) % nmat;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
if (ty1 == 'o') {
|
|
rpm-build |
8267b0 |
jk1r = (jj + 2) % nmat;
|
|
rpm-build |
8267b0 |
jk2r = (jj + 3) % nmat;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
add_mat_line(m, v, derivs[0][0], th - ends[0][0], 1, j, jthl, jinc, nmat);
|
|
rpm-build |
8267b0 |
add_mat_line(m, v, derivs[1][0], ends[0][1], -1, j, jk0l, jinc, nmat);
|
|
rpm-build |
8267b0 |
add_mat_line(m, v, derivs[2][0], ends[0][2], -1, j, jk1l, jinc, nmat);
|
|
rpm-build |
8267b0 |
add_mat_line(m, v, derivs[3][0], ends[0][3], -1, j, jk2l, jinc, nmat);
|
|
rpm-build |
8267b0 |
add_mat_line(m, v, derivs[0][1], -ends[1][0], 1, j, jthr, jinc, nmat);
|
|
rpm-build |
8267b0 |
add_mat_line(m, v, derivs[1][1], -ends[1][1], 1, j, jk0r, jinc, nmat);
|
|
rpm-build |
8267b0 |
add_mat_line(m, v, derivs[2][1], -ends[1][2], 1, j, jk1r, jinc, nmat);
|
|
rpm-build |
8267b0 |
add_mat_line(m, v, derivs[3][1], -ends[1][3], 1, j, jk2r, jinc, nmat);
|
|
rpm-build |
8267b0 |
if (jthl >= 0)
|
|
rpm-build |
8267b0 |
v[jthl] = mod_2pi(v[jthl]);
|
|
rpm-build |
8267b0 |
if (jthr >= 0)
|
|
rpm-build |
8267b0 |
v[jthr] = mod_2pi(v[jthr]);
|
|
rpm-build |
8267b0 |
j += jinc;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
if (cyclic) {
|
|
rpm-build |
8267b0 |
memcpy(m + nmat, m, sizeof(bandmat) * nmat);
|
|
rpm-build |
8267b0 |
memcpy(m + 2 * nmat, m, sizeof(bandmat) * nmat);
|
|
rpm-build |
8267b0 |
memcpy(v + nmat, v, sizeof(double) * nmat);
|
|
rpm-build |
8267b0 |
memcpy(v + 2 * nmat, v, sizeof(double) * nmat);
|
|
rpm-build |
8267b0 |
n_invert = 3 * nmat;
|
|
rpm-build |
8267b0 |
j = nmat;
|
|
rpm-build |
8267b0 |
} else {
|
|
rpm-build |
8267b0 |
n_invert = nmat;
|
|
rpm-build |
8267b0 |
j = 0;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#ifdef VERBOSE
|
|
rpm-build |
8267b0 |
for (i = 0; i < n; i++) {
|
|
rpm-build |
8267b0 |
int k;
|
|
rpm-build |
8267b0 |
for (k = 0; k < 11; k++)
|
|
rpm-build |
8267b0 |
printf(" %2.4f", m[i].a[k]);
|
|
rpm-build |
8267b0 |
printf(": %2.4f\n", v[i]);
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
printf("---\n");
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
bandec11(m, perm, n_invert);
|
|
rpm-build |
8267b0 |
banbks11(m, perm, v, n_invert);
|
|
rpm-build |
8267b0 |
norm = 0.;
|
|
rpm-build |
8267b0 |
for (i = 0; i < n; i++) {
|
|
rpm-build |
8267b0 |
int jinc = compute_jinc(s[i].ty, s[i + 1].ty);
|
|
rpm-build |
8267b0 |
int k;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
for (k = 0; k < jinc; k++) {
|
|
rpm-build |
8267b0 |
double dk = v[j++];
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
#ifdef VERBOSE
|
|
rpm-build |
8267b0 |
printf("s[%d].ks[%d] += %f\n", i, k, dk);
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
s[i].ks[k] += dk;
|
|
rpm-build |
8267b0 |
norm += dk * dk;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
s[i].ks[0] = 2.0*mod_2pi(s[i].ks[0]/2.0);
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
return norm;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
static int
|
|
rpm-build |
8267b0 |
check_finiteness(spiro_seg * segs, int num_segs)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
/* Check if all values are "finite", return true=0, else return fail=-1 */
|
|
rpm-build |
8267b0 |
int i, j;
|
|
rpm-build |
8267b0 |
for (i = 0; i < num_segs; ++i)
|
|
rpm-build |
8267b0 |
for (j = 0; j < 4; ++j)
|
|
rpm-build |
8267b0 |
if ( IS_FINITE( segs[i].ks[j])==0 ) return -1;
|
|
rpm-build |
8267b0 |
return 0;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
static int
|
|
rpm-build |
8267b0 |
solve_spiro(spiro_seg *s, int nseg)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
bandmat *m;
|
|
rpm-build |
8267b0 |
double *v;
|
|
rpm-build |
8267b0 |
int *perm;
|
|
rpm-build |
8267b0 |
int nmat = count_vec(s, nseg);
|
|
rpm-build |
8267b0 |
int n_alloc = nmat;
|
|
rpm-build |
8267b0 |
double norm;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
if (nmat == 0)
|
|
rpm-build |
8267b0 |
return 1; // just means no convergence problems
|
|
rpm-build |
8267b0 |
if (s[0].ty != '{' && s[0].ty != 'v')
|
|
rpm-build |
8267b0 |
n_alloc *= 3;
|
|
rpm-build |
8267b0 |
if (n_alloc < 5)
|
|
rpm-build |
8267b0 |
n_alloc = 5;
|
|
rpm-build |
8267b0 |
m = (bandmat *)malloc(sizeof(bandmat) * n_alloc);
|
|
rpm-build |
8267b0 |
v = (double *)malloc(sizeof(double) * n_alloc);
|
|
rpm-build |
8267b0 |
perm = (int *)malloc(sizeof(int) * n_alloc);
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
int i = 0, converged = 0; // not solved (yet)
|
|
rpm-build |
8267b0 |
if ( m!=NULL && v!=NULL && perm!=NULL ) {
|
|
rpm-build |
8267b0 |
while (i++ < 60) {
|
|
rpm-build |
8267b0 |
norm = spiro_iter(s, m, perm, v, nseg, nmat);
|
|
rpm-build |
8267b0 |
#ifdef VERBOSE
|
|
rpm-build |
8267b0 |
printf("iteration #%d, %% norm = %g\n", i, norm);
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
if (check_finiteness(s, nseg)) break;
|
|
rpm-build |
8267b0 |
if (norm < 1e-12) { converged = 1; break; }
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
#ifdef VERBOSE
|
|
rpm-build |
8267b0 |
if (converged==0)
|
|
rpm-build |
8267b0 |
fprintf(stderr, "ERROR: LibSpiro: failed to converge after %d attempts to converge.\n", i);
|
|
rpm-build |
8267b0 |
} else {
|
|
rpm-build |
8267b0 |
fprintf(stderr, "ERROR: LibSpiro: failed to alloc memory.\n");
|
|
rpm-build |
8267b0 |
#endif
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
free(m);
|
|
rpm-build |
8267b0 |
free(v);
|
|
rpm-build |
8267b0 |
free(perm);
|
|
rpm-build |
8267b0 |
return converged;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
static void
|
|
rpm-build |
8267b0 |
spiro_seg_to_bpath(const double ks[4],
|
|
rpm-build |
8267b0 |
double x0, double y0, double x1, double y1,
|
|
rpm-build |
8267b0 |
bezctx *bc, int depth)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
double bend = fabs(ks[0]) + fabs(.5 * ks[1]) + fabs(.125 * ks[2]) +
|
|
rpm-build |
8267b0 |
fabs((1./48) * ks[3]);
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
if (bend <= 1e-8) {
|
|
rpm-build |
8267b0 |
bezctx_lineto(bc, x1, y1);
|
|
rpm-build |
8267b0 |
} else {
|
|
rpm-build |
8267b0 |
double seg_ch = hypot(x1 - x0, y1 - y0);
|
|
rpm-build |
8267b0 |
double seg_th = atan2(y1 - y0, x1 - x0);
|
|
rpm-build |
8267b0 |
double xy[2];
|
|
rpm-build |
8267b0 |
double ch, th;
|
|
rpm-build |
8267b0 |
double scale, rot;
|
|
rpm-build |
8267b0 |
double th_even, th_odd;
|
|
rpm-build |
8267b0 |
double ul, vl;
|
|
rpm-build |
8267b0 |
double ur, vr;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
integrate_spiro(ks, xy, N_IS);
|
|
rpm-build |
8267b0 |
ch = hypot(xy[0], xy[1]);
|
|
rpm-build |
8267b0 |
th = atan2(xy[1], xy[0]);
|
|
rpm-build |
8267b0 |
scale = seg_ch / ch;
|
|
rpm-build |
8267b0 |
rot = seg_th - th;
|
|
rpm-build |
8267b0 |
if (depth > 5 || bend < 1.) {
|
|
rpm-build |
8267b0 |
th_even = (1./384) * ks[3] + (1./8) * ks[1] + rot;
|
|
rpm-build |
8267b0 |
th_odd = (1./48) * ks[2] + .5 * ks[0];
|
|
rpm-build |
8267b0 |
ul = (scale * (1./3)) * cos(th_even - th_odd);
|
|
rpm-build |
8267b0 |
vl = (scale * (1./3)) * sin(th_even - th_odd);
|
|
rpm-build |
8267b0 |
ur = (scale * (1./3)) * cos(th_even + th_odd);
|
|
rpm-build |
8267b0 |
vr = (scale * (1./3)) * sin(th_even + th_odd);
|
|
rpm-build |
8267b0 |
bezctx_curveto(bc, x0 + ul, y0 + vl, x1 - ur, y1 - vr, x1, y1);
|
|
rpm-build |
8267b0 |
} else {
|
|
rpm-build |
8267b0 |
/* subdivide */
|
|
rpm-build |
8267b0 |
double ksub[4];
|
|
rpm-build |
8267b0 |
double thsub;
|
|
rpm-build |
8267b0 |
double xysub[2];
|
|
rpm-build |
8267b0 |
double xmid, ymid;
|
|
rpm-build |
8267b0 |
double cth, sth;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
ksub[0] = .5 * ks[0] - .125 * ks[1] + (1./64) * ks[2] - (1./768) * ks[3];
|
|
rpm-build |
8267b0 |
ksub[1] = .25 * ks[1] - (1./16) * ks[2] + (1./128) * ks[3];
|
|
rpm-build |
8267b0 |
ksub[2] = .125 * ks[2] - (1./32) * ks[3];
|
|
rpm-build |
8267b0 |
ksub[3] = (1./16) * ks[3];
|
|
rpm-build |
8267b0 |
thsub = rot - .25 * ks[0] + (1./32) * ks[1] - (1./384) * ks[2] + (1./6144) * ks[3];
|
|
rpm-build |
8267b0 |
cth = .5 * scale * cos(thsub);
|
|
rpm-build |
8267b0 |
sth = .5 * scale * sin(thsub);
|
|
rpm-build |
8267b0 |
integrate_spiro(ksub, xysub, N_IS);
|
|
rpm-build |
8267b0 |
xmid = x0 + cth * xysub[0] - sth * xysub[1];
|
|
rpm-build |
8267b0 |
ymid = y0 + cth * xysub[1] + sth * xysub[0];
|
|
rpm-build |
8267b0 |
spiro_seg_to_bpath(ksub, x0, y0, xmid, ymid, bc, depth + 1);
|
|
rpm-build |
8267b0 |
ksub[0] += .25 * ks[1] + (1./384) * ks[3];
|
|
rpm-build |
8267b0 |
ksub[1] += .125 * ks[2];
|
|
rpm-build |
8267b0 |
ksub[2] += (1./16) * ks[3];
|
|
rpm-build |
8267b0 |
spiro_seg_to_bpath(ksub, xmid, ymid, x1, y1, bc, depth + 1);
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
spiro_seg *
|
|
rpm-build |
8267b0 |
run_spiro(const spiro_cp *src, int n)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
if (src==NULL || n <= 0) return 0;
|
|
rpm-build |
8267b0 |
int nseg = src[0].ty == '{' ? n - 1 : n;
|
|
rpm-build |
8267b0 |
spiro_seg *s = setup_path(src, n);
|
|
rpm-build |
8267b0 |
if (s) {
|
|
rpm-build |
8267b0 |
int converged = 1 ; // this value is for when nseg == 1; else actual value determined below
|
|
rpm-build |
8267b0 |
if (nseg > 1) converged = solve_spiro(s, nseg);
|
|
rpm-build |
8267b0 |
if (converged) return s;
|
|
rpm-build |
8267b0 |
free(s);
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
return 0;
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
void
|
|
rpm-build |
8267b0 |
free_spiro(spiro_seg *s)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
if (s) free(s);
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
void
|
|
rpm-build |
8267b0 |
spiro_to_bpath(const spiro_seg *s, int n, bezctx *bc)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
if (s==NULL || n <= 0 || bc==NULL) return;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
int i;
|
|
rpm-build |
8267b0 |
int nsegs = s[n - 1].ty == '}' ? n - 1 : n;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
for (i = 0; i < nsegs; i++) {
|
|
rpm-build |
8267b0 |
double x0 = s[i].x;
|
|
rpm-build |
8267b0 |
double y0 = s[i].y;
|
|
rpm-build |
8267b0 |
double x1 = s[i + 1].x;
|
|
rpm-build |
8267b0 |
double y1 = s[i + 1].y;
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
if (i == 0)
|
|
rpm-build |
8267b0 |
bezctx_moveto(bc, x0, y0, s[0].ty == '{');
|
|
rpm-build |
8267b0 |
bezctx_mark_knot(bc, i);
|
|
rpm-build |
8267b0 |
spiro_seg_to_bpath(s[i].ks, x0, y0, x1, y1, bc, 0);
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
double
|
|
rpm-build |
8267b0 |
get_knot_th(const spiro_seg *s, int i)
|
|
rpm-build |
8267b0 |
{
|
|
rpm-build |
8267b0 |
double ends[2][4];
|
|
rpm-build |
8267b0 |
|
|
rpm-build |
8267b0 |
if (i == 0) {
|
|
rpm-build |
8267b0 |
compute_ends(s[i].ks, ends, s[i].seg_ch);
|
|
rpm-build |
8267b0 |
return s[i].seg_th - ends[0][0];
|
|
rpm-build |
8267b0 |
} else {
|
|
rpm-build |
8267b0 |
compute_ends(s[i - 1].ks, ends, s[i - 1].seg_ch);
|
|
rpm-build |
8267b0 |
return s[i - 1].seg_th + ends[1][0];
|
|
rpm-build |
8267b0 |
}
|
|
rpm-build |
8267b0 |
}
|