1
#define _ISOC99_SOURCE /* to get isfinite() */
11
typedef struct Vecfunc
18
static VECFUNC vf[] = {
19
{"v_copy", v_copy, "p=rp"},
20
{"v_add", v_add, "p=rpp"},
21
{"pnt_op_func_+", v_add, "p=rpp"},
22
{"v_sub", v_sub, "p=rpp"},
23
{"pnt_op_func_-", v_sub, "p=rpp"},
24
{"v_abs", v_abs, "p=rp"},
25
{"v_neg", v_neg, "p=rp"},
26
{"pnt_op_func__", v_neg, "p=rp"},
27
{"v_mul", v_mul, "p=rpd"},
28
{"pnt_op_func_*", v_mul, "p=rpd"},
29
{"v_div", v_div, "p=rpd"},
30
{"pnt_op_func_/", v_div, "p=rpd"},
31
{"v_unit", v_unit, "p=rp"},
32
{"v_cross", v_cross, "p=rpp"},
33
{"pnt_op_func_^", v_cross, "p=rpp"},
34
{"v_val", v_val, "d=p"},
35
{"v_dot", v_dot, "d=pp"},
36
{"pnt_op_func_%", v_dot, "d=pp"},
37
{"v_area", v_area, "d=pp"},
38
{"v_eq", v_eq, "d=pp"},
39
{"v_eq_epsilon", v_eq_epsilon, "d=ppp"},
40
{"v_isortho", v_isortho, "d=pp"},
41
{"v_ispara", v_ispara, "d=pp"},
42
{"v_isacute", v_isacute, "d=pp"},
46
typedef VECTOR *(*p_func) (void);
47
typedef VECTOR *(*p_func_p) (void *p0);
48
typedef VECTOR *(*p_func_pp) (void *p0, void *p1);
49
typedef VECTOR *(*p_func_ppp) (void *p0, void *p1, void *p2);
50
typedef VECTOR *(*p_func_ppd) (void *p0, void *p1, double d);
53
static VECTOR pnt_o = { NULL, 0, 0, 0, 1 };
54
static VECTOR pnt_i = { NULL, 1, 0, 0, 1 };
55
static VECTOR pnt_j = { NULL, 0, 1, 0, 1 };
56
static VECTOR pnt_k = { NULL, 0, 0, 1, 1 };
59
void printvec(SYMBOL * sym);
60
void showvec(SYMBOL * sym);
61
void setpnt(SYMBOL * var, SYMBOL * pnt);
62
SYMBOL *mkpnt(double x, double y, double z);
63
SYMBOL *mkpntvar(SYMBOL * var, SYMBOL * pnt);
64
SYMBOL *pntfunc(SYMBOL * func, SYMBOL * arglist);
65
SYMBOL *pntop(int op, SYMBOL * pnt1, SYMBOL * pnt2);
66
SYMBOL *pntapp(SYMBOL * head, SYMBOL * elt);
67
VECTOR *v_copy(VECTOR * p, VECTOR * p1);
68
VECTOR *v_add(VECTOR * p, VECTOR * p1, VECTOR * p2);
69
VECTOR *v_sub(VECTOR * p, VECTOR * p1, VECTOR * p2);
70
VECTOR *v_abs(VECTOR * p, VECTOR * p1);
71
VECTOR *v_neg(VECTOR * p, VECTOR * p1);
72
static inline int _is_zero(double r);
73
double v_eq(VECTOR * p1, VECTOR * p2);
74
double v_eq_epsilon(VECTOR * p1, VECTOR * p2, VECTOR * e);
75
VECTOR *v_mul(VECTOR * p, VECTOR * p1, double d);
76
VECTOR *v_div(VECTOR * p, VECTOR * p1, double d);
77
double v_val(VECTOR * p);
78
VECTOR *v_unit(VECTOR * p, VECTOR * p1);
79
double v_dot(VECTOR * p1, VECTOR * p2);
80
VECTOR *v_cross(VECTOR * p, VECTOR * p1, VECTOR * p2);
81
double v_isortho(VECTOR * p1, VECTOR * p2);
82
double v_ispara(VECTOR * p1, VECTOR * p2);
83
double v_isacute(VECTOR * p1, VECTOR * p2);
84
double v_area(VECTOR * p1, VECTOR * p2);
92
for (i = 0; vf[i].name; i++) {
93
sym = putsym(vf[i].name);
94
switch (vf[i].proto[0]) {
97
sym->rettype = st_pnt;
100
sym->type = st_nfunc;
101
sym->rettype = st_num;
104
sym->itype = sym->type;
105
sym->v.p = vf[i].func;
106
sym->proto = vf[i].proto + 2;
109
/* add some handy constants */
110
sym = putsym("pnt_o");
111
sym->type = sym->itype = st_pnt;
114
sym = putsym("pnt_i");
115
sym->type = sym->itype = st_pnt;
118
sym = putsym("pnt_j");
119
sym->type = sym->itype = st_pnt;
122
sym = putsym("pnt_k");
123
sym->type = sym->itype = st_pnt;
130
void printvec(SYMBOL * sym)
134
v = (VECTOR *) sym->v.p;
136
fprintf(stdout, "\t(");
138
fprintf(stdout, "??.??");
139
else if (v->x == (int)v->x)
140
fprintf(stdout, "%d", (int)v->x);
142
fprintf(stdout, "%g", v->x);
143
fprintf(stdout, ", ");
145
fprintf(stdout, "??.??");
146
else if (v->y == (int)v->y)
147
fprintf(stdout, "%d", (int)v->y);
149
fprintf(stdout, "%g", v->y);
150
if (isfinite(v->z)) {
151
fprintf(stdout, ", ");
152
if (v->z == (int)v->z)
153
fprintf(stdout, "%d", (int)v->z);
155
fprintf(stdout, "%g", v->z);
157
fprintf(stdout, ")\n");
160
void showvec(SYMBOL * sym)
165
v = (VECTOR *) sym->v.p;
168
if (v && --v->refcnt > 0)
173
void setpnt(SYMBOL * var, SYMBOL * pnt)
178
sym = getsym(var->name);
180
if (--((VECTOR *) sym->v.p)->refcnt < 1)
183
* If refcnt(pnt) == 1, this was anonymous, else it's used
184
* somewhere else. Must we dup then?
190
if (--((VECTOR *) var->v.p)->refcnt < 1)
200
SYMBOL *mkpnt(double x, double y, double z)
205
vec = (VECTOR *) listitem(sizeof(VECTOR));
211
pnt = (SYMBOL *) listitem(sizeof(SYMBOL));
212
pnt->type = pnt->itype = st_pnt;
218
SYMBOL *mkpntvar(SYMBOL * var, SYMBOL * pnt)
220
var->type = var->itype = st_pnt;
221
var->name = var->v.p;
226
symtab = (SYMBOL *) listadd((LIST *) symtab, (LIST *) var, cmpsymsym);
233
SYMBOL *pntfunc(SYMBOL * func, SYMBOL * arglist)
239
sym = (SYMBOL *) listitem(sizeof(SYMBOL));
240
sym->type = sym->itype = st_pnt;
242
if (!func || !func->v.p || func->type != st_pfunc) {
244
G_warning(_("Can't call bad function"));
247
argc = listcnt((LIST *) arglist);
249
for (i = 0, sptr = arglist; sptr; sptr = sptr->next, i++) {
250
if (func->proto[i] == 'r')
252
if (func->proto[i] == 'p') {
253
if (sptr->itype != st_pnt || !sptr->v.p)
256
else if (func->proto[i] == 'd' && sptr->itype != st_num)
260
res = (VECTOR *) listitem(sizeof(VECTOR));
262
if (argc == 0 && (!func->proto || !*func->proto))
263
res = (*(p_func) func->v.p) ();
264
else if (argc == 1 && !strcmp(func->proto, "p"))
265
res = (*(p_func_p) func->v.p) (arglist->v.p);
266
else if (argc == 1 && !strcmp(func->proto, "rp"))
267
res = (*(p_func_pp) func->v.p) (res, arglist->v.p);
268
else if (argc == 2 && !strcmp(func->proto, "rpd"))
269
res = (*(p_func_ppd) func->v.p) (res, arglist->v.p,
271
else if (argc == 2 && !strcmp(func->proto, "pp"))
272
res = (*(p_func_pp) func->v.p) (arglist->v.p, arglist->next->v.p);
273
else if (argc == 2 && !strcmp(func->proto, "rpp"))
274
res = (*(p_func_ppp) func->v.p) (res, arglist->v.p,
276
else if (argc == 3 && !strcmp(func->proto, "ppp"))
277
res = (*(p_func_ppp) func->v.p) (arglist->v.p,
279
arglist->next->next->v.p);
281
G_warning(_("Bad arguments to pointfunc %s"), func->name);
283
sym = (SYMBOL *) listdelall((LIST *) sym, freesym);
290
listdelall((LIST *) arglist, freesym);
295
SYMBOL *pntop(int op, SYMBOL * pnt1, SYMBOL * pnt2)
297
SYMBOL *func, *arglist, *res = NULL;
300
sprintf(buf, "pnt_op_func_%c", op);
304
G_warning(_("No function defined to perform ``point %c point''"), op);
307
else if (func->rettype == st_pnt) {
308
res = (SYMBOL *) listitem(sizeof(SYMBOL));
312
arglist = (SYMBOL *) listapp(NULL, (LIST *) pnt1);
313
arglist = (SYMBOL *) listapp((LIST *) arglist, (LIST *) pnt2);
315
res = pntfunc(func, arglist);
321
SYMBOL *pntapp(SYMBOL * head, SYMBOL * elt)
323
return (SYMBOL *) listapp((LIST *) head, (LIST *) elt);
327
* Utility function to copy a point: p = p1;
328
* The dimension (2D/3D) depends on p1. Note, that copying a constant
329
* will always yield 3D.
331
VECTOR *v_copy(VECTOR * p, VECTOR * p1)
342
* Result is 2D if at least one of p1 or p2 is 2D.
344
VECTOR *v_add(VECTOR * p, VECTOR * p1, VECTOR * p2)
346
p->x = p1->x + p2->x;
347
p->y = p1->y + p2->y;
349
* resist the tentation setting p->z to nanval and then testing for
350
* dimension, as p might be the same as p1 or p2
352
if (!isnan(p1->z) && !isnan(p2->z))
353
p->z = p1->z + p2->z;
361
* Vector substraction
362
* Result is 2D if at least one of p1 or p2 is 2D.
364
VECTOR *v_sub(VECTOR * p, VECTOR * p1, VECTOR * p2)
366
p->x = p1->x - p2->x;
367
p->y = p1->y - p2->y;
368
if (!isnan(p1->z) && !isnan(p2->z))
369
p->z = p1->z + p2->z;
377
* Utility function to make all coordinates positive
379
VECTOR *v_abs(VECTOR * p, VECTOR * p1)
392
* Utility function to negate all coordinates
394
VECTOR *v_neg(VECTOR * p, VECTOR * p1)
407
* Utility function to compare two doubles for equality without epsilon
408
* This is not really true, as we consider NaN to be zero.
410
static inline int _is_zero(double r)
412
if ((isfinite(r) && r == 0.0) || isnan(r))
418
* Test for equality of two points. No epsion applied
420
double v_eq(VECTOR * p1, VECTOR * p2)
425
if (!isnan(p1->z) && !isnan(p2->z))
430
if (_is_zero(p.x) && _is_zero(p.y) && (dim == 2 || _is_zero(p.z)))
436
* Test for equality of two points by a given epsilon
437
* epsilon is supposed to have positive values only.
439
double v_eq_epsilon(VECTOR * p1, VECTOR * p2, VECTOR * e)
444
if (!isnan(p1->z) && !isnan(p2->z))
449
if (p.x < e->x && p.y < e->y && (dim == 2 || (p.z < e->z)))
455
* Multiply a vector by a scalar
457
VECTOR *v_mul(VECTOR * p, VECTOR * p1, double d)
470
* Divide a vector by a scalar
472
VECTOR *v_div(VECTOR * p, VECTOR * p1, double d)
474
if (!isfinite(d) || d == 0.0) {
489
* Compute the value of a vector
491
double v_val(VECTOR * p)
493
return sqrt(p->x * p->x + p->y * p->y +
494
((isnan(p->z)) ? 0 : p->z * p->z));
498
* The only way to get a value of zero is that P1 is the origin.
499
* The unit vector of the origin doesn't exist, but we return the
502
VECTOR *v_unit(VECTOR * p, VECTOR * p1)
504
double val = v_val(p1);
507
return v_copy(p, &pnt_o);
520
* Compute the dot product of P1 and P2
522
double v_dot(VECTOR * p1, VECTOR * p2)
526
if (!isnan(p1->z) && !isnan(p2->z))
528
return p1->x * p2->x + p1->y * p2->y + ((dim == 2) ? 0 : p1->z * p2->z);
532
* Compute the cross product of P1 and P2
533
* Return (0,0) for 2D
536
VECTOR *v_cross(VECTOR * p, VECTOR * p1, VECTOR * p2)
540
if (!isnan(p1->z) && !isnan(p2->z)) {
541
p0.x = p1->y * p2->z + p1->z * p2->y;
542
p0.y = p1->z * p2->x + p1->x * p2->z;
543
p0.z = p1->x * p2->y + p1->y * p2->x;
555
* Decide if vector P1 is ortogonal to vector P2
556
* Should test if either P1 or P2 are (0,0,0);
558
double v_isortho(VECTOR * p1, VECTOR * p2)
560
return v_dot(p1, p2) == 0;
564
* Decide if P1 and P2 are parallel. If they are but have a diferent
565
* direction, -1 is returned.
567
double v_ispara(VECTOR * p1, VECTOR * p2)
569
double dot, val, dif;
575
dif = fabs(dot - val);
579
dif = fabs(dot + val);
587
* Decide if P1 and P2 have and angle alpha which: 0 < alpha < 90.0
589
double v_isacute(VECTOR * p1, VECTOR * p2)
591
return v_dot(p1, p2) > 0;
595
* Return the area spanned by the two vectors P1 and P2
598
double v_area(VECTOR * p1, VECTOR * p2)
602
return 0.5 * v_val(v_cross(&p, p1, p2));