1 /*
2 * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
3 * Applied Mathematics, Norway.
4 *
5 * Contact information: E-mail: tor.dokken@sintef.no
6 * SINTEF ICT, Department of Applied Mathematics,
7 * P.O. Box 124 Blindern,
8 * 0314 Oslo, Norway.
9 *
10 * This file is part of SISL.
11 *
12 * SISL is free software: you can redistribute it and/or modify
13 * it under the terms of the GNU Affero General Public License as
14 * published by the Free Software Foundation, either version 3 of the
15 * License, or (at your option) any later version.
16 *
17 * SISL is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Affero General Public License for more details.
21 *
22 * You should have received a copy of the GNU Affero General Public
23 * License along with SISL. If not, see
24 * <http://www.gnu.org/licenses/>.
25 *
26 * In accordance with Section 7(b) of the GNU Affero General Public
27 * License, a covered work must retain the producer line in every data
28 * file that is created or manipulated using SISL.
29 *
30 * Other Usage
31 * You can be released from the requirements of the license by purchasing
32 * a commercial license. Buying such a license is mandatory as soon as you
33 * develop commercial activities involving the SISL library without
34 * disclosing the source code of your own applications.
35 *
36 * This file may be used in accordance with the terms contained in a
37 * written agreement between you and SINTEF ICT.
38 */
39
40 #include "sisl-copyright.h"
41
42 /*
43 *
44 * $Id: s1615.c,v 1.2 2001-03-19 15:58:52 afr Exp $
45 *
46 */
47
48
49 #define S1615
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
s1615(double epoint[],int inbpnt,int idim,int eptyp[],int * jstat)54 void s1615(double epoint[], int inbpnt, int idim, int eptyp[],int *jstat)
55 #else
56 void s1615(epoint, inbpnt, idim, eptyp, jstat)
57 double epoint[];
58 int inbpnt;
59 int idim;
60 int eptyp[];
61 int *jstat;
62 #endif
63 /*
64 *************************************************************************
65 *
66 * PURPOSE: To test if the points lie on two branches of a hyperbola
67 *
68 * INPUT:
69 * Epoint - Points/tangents describing the conic.
70 * Inbpnt - No. of points/tangents in the epoint array.
71 * Idim - The dimension of the space in which the points lie.
72 * Eptyp - Type indicator for points/tangents :
73 * 1 - Ordinary point.
74 * 2 - Knuckle point. (Is treated as an ordinary point.)
75 * 3 - Tangent to next point.
76 * 4 - Tangent to prior point.
77 *
78 * Output:
79 * Jstat - status messages:
80 * < 0 : Error.
81 * = 0 : Ok.
82 * > 0 : Warning.
83 *
84 * Method:
85 * If 3D points, the subroutine assumes that the third coordinate
86 * is zero, and that the first entry in Epoint is a point.
87 * Derivative vectors and vectors between the points are extracted
88 * and put into a temporary array. Then the cross product between
89 * adjacent ones of these are calculated. If both positive and negative
90 * cross products are found, the points lie on different branches
91 * of a hyperbola.
92 *-
93 * Calls: s6err.
94 *
95 * Written by: A.M. Ytrehus, SI Oslo, Oct.91.
96 * After FORTRAN, (P1615), written by: T. Dokken SI.
97 *****************************************************************
98 */
99 {
100 double svector[8];
101 double *spoint = SISL_NULL;
102 int ki, kk, kp, kki;
103 int ktyp;
104 double tcross;
105 int kant = 4;
106 int kpluss = 0;
107 int kneg = 0;
108 int kdim = 2;
109 int kstore = 0;
110 int kpos = 0;
111
112 *jstat = 0;
113
114
115 /* Allocate local array. */
116
117 spoint = newarray (kdim * inbpnt, DOUBLE);
118 if (spoint == SISL_NULL)
119 goto err101;
120
121
122 /* Kant is the no. of segments between the points (that is inbpnt-1).
123 We can have max. 5 points. */
124
125 if (inbpnt < 5)
126 kant = inbpnt - 1;
127
128
129 /* If the no. of segments is two or less, the points cannot
130 belong to two different conic branches. (Ok, go out). */
131
132 if (kant < 3)
133 goto out;
134
135
136 /* Store the positions of points and tangents (next/prior point
137 minus/plus tangent/vector) in the working array spoit. */
138
139
140 for (ki = 0; ki < inbpnt; ki++)
141 {
142 ktyp = eptyp[ki];
143 kki = idim * ki;
144
145 if (ktyp == 1 || ktyp == 2)
146 {
147 /* Store this point directly. */
148
149 kk = kdim * kstore;
150
151 for (kp = 0; kp < kdim; kp++)
152 spoint[kk + kp] = epoint[kki + kp];
153
154 kstore++;
155 }
156 else if (ktyp == 3)
157 {
158 /* Tangent to next point. Store the coordinates of
159 next point minus this one. */
160
161 kk = kdim * kstore;
162
163 for (kp = 0; kp < kdim; kp++)
164 spoint[kk + kp] = epoint[kki + idim + kp] - epoint[kki + kp];
165
166 kstore++;
167
168 }
169 else if (ktyp == 4)
170 {
171 /* Tangent to prior point. Store the coordinates of
172 last point plus this one. */
173
174 kk = kdim * kstore;
175
176 for (kp = 0; kp < kdim; kp++)
177 spoint[kk + kp] = epoint[kki - idim + kp] + epoint[kki + kp];
178
179 kstore++;
180 }
181 }
182
183
184 /* Run through the points/derivatives and put difference
185 vectors into the svector-array. */
186
187
188 kstore = 0;
189
190 for (ki = 1; ki < inbpnt; ki++)
191 {
192 kk = kdim * kstore;
193 kki = kdim * ki;
194
195 for (kp = 0; kp < kdim; kp++)
196 svector[kk + kp] = spoint[kki + kp] - spoint[kki - kdim + kp];
197
198 kstore++;
199 }
200
201
202 /* If the cross product of adjacent vectors all have the same
203 direction, then the points all lie on the same conic branch. */
204
205 for (ki = 0; ki < kant - 1; ki++)
206 {
207 kk = kdim * ki;
208
209 tcross = svector[kk] * svector[kk + 3]
210 - svector[kk + 1] * svector[kk + 2];
211
212 if (tcross > 0)
213 kpluss++;
214 if (tcross < 0)
215 kneg++;
216 }
217
218 /* If both kpluss and kneg are nonzero, then the points
219 lie on different branches. */
220
221 if (kpluss > 0 && kneg > 0)
222 *jstat = 1;
223
224 goto out;
225
226 /* Allocation error. */
227
228 err101:
229 *jstat = -101;
230 s6err ("s1615", *jstat, kpos);
231 goto out;
232
233 out:
234
235 /* Free local arrays. */
236
237 if (spoint != SISL_NULL)
238 freearray (spoint);
239
240 return;
241 }
242