~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to vector/v.mapcalc/vector.c

Tags: 7.0.0~rc1+ds1-1~exp1
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
#define _ISOC99_SOURCE          /* to get isfinite() */
2
 
#include <stdio.h>
3
 
#include <stdlib.h>
4
 
#include <string.h>
5
 
#include <math.h>
6
 
#include <grass/gis.h>
7
 
#include "list.h"
8
 
#include "mapcalc.h"
9
 
#include "vector.h"
10
 
 
11
 
typedef struct Vecfunc
12
 
{
13
 
    char *name;
14
 
    void *func;
15
 
    char *proto;
16
 
} VECFUNC;
17
 
 
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"},
43
 
    {NULL, NULL, NULL}
44
 
};
45
 
 
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);
51
 
 
52
 
double nanval;
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 };
57
 
 
58
 
void init_vec(void);
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);
85
 
 
86
 
 
87
 
void init_vec(void)
88
 
{
89
 
    SYMBOL *sym;
90
 
    int i;
91
 
 
92
 
    for (i = 0; vf[i].name; i++) {
93
 
        sym = putsym(vf[i].name);
94
 
        switch (vf[i].proto[0]) {
95
 
        case 'p':
96
 
            sym->type = st_pfunc;
97
 
            sym->rettype = st_pnt;
98
 
            break;
99
 
        case 'd':
100
 
            sym->type = st_nfunc;
101
 
            sym->rettype = st_num;
102
 
            break;
103
 
        }
104
 
        sym->itype = sym->type;
105
 
        sym->v.p = vf[i].func;
106
 
        sym->proto = vf[i].proto + 2;
107
 
    }
108
 
 
109
 
    /* add some handy constants */
110
 
    sym = putsym("pnt_o");
111
 
    sym->type = sym->itype = st_pnt;
112
 
    sym->v.p = &pnt_o;
113
 
 
114
 
    sym = putsym("pnt_i");
115
 
    sym->type = sym->itype = st_pnt;
116
 
    sym->v.p = &pnt_i;
117
 
 
118
 
    sym = putsym("pnt_j");
119
 
    sym->type = sym->itype = st_pnt;
120
 
    sym->v.p = &pnt_j;
121
 
 
122
 
    sym = putsym("pnt_k");
123
 
    sym->type = sym->itype = st_pnt;
124
 
    sym->v.p = &pnt_k;
125
 
 
126
 
    /* initialize NaN */
127
 
    nanval = sqrt(-1);
128
 
}
129
 
 
130
 
void printvec(SYMBOL * sym)
131
 
{
132
 
    VECTOR *v;
133
 
 
134
 
    v = (VECTOR *) sym->v.p;
135
 
 
136
 
    fprintf(stdout, "\t(");
137
 
    if (!isfinite(v->x))
138
 
        fprintf(stdout, "??.??");
139
 
    else if (v->x == (int)v->x)
140
 
        fprintf(stdout, "%d", (int)v->x);
141
 
    else
142
 
        fprintf(stdout, "%g", v->x);
143
 
    fprintf(stdout, ", ");
144
 
    if (!isfinite(v->y))
145
 
        fprintf(stdout, "??.??");
146
 
    else if (v->y == (int)v->y)
147
 
        fprintf(stdout, "%d", (int)v->y);
148
 
    else
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);
154
 
        else
155
 
            fprintf(stdout, "%g", v->z);
156
 
    }
157
 
    fprintf(stdout, ")\n");
158
 
}
159
 
 
160
 
void showvec(SYMBOL * sym)
161
 
{
162
 
 
163
 
    VECTOR *v;
164
 
 
165
 
    v = (VECTOR *) sym->v.p;
166
 
    printvec(sym);
167
 
 
168
 
    if (v && --v->refcnt > 0)
169
 
        sym->v.p = NULL;
170
 
    freesym(sym);
171
 
}
172
 
 
173
 
void setpnt(SYMBOL * var, SYMBOL * pnt)
174
 
{
175
 
    SYMBOL *sym;
176
 
 
177
 
    if (var->name) {
178
 
        sym = getsym(var->name);
179
 
        if (sym) {
180
 
            if (--((VECTOR *) sym->v.p)->refcnt < 1)
181
 
                G_free(sym->v.p);
182
 
            /*
183
 
             * If refcnt(pnt) == 1, this was anonymous, else it's used
184
 
             * somewhere else. Must we dup then?
185
 
             */
186
 
            sym->v.p = pnt->v.p;
187
 
        }
188
 
    }
189
 
 
190
 
    if (--((VECTOR *) var->v.p)->refcnt < 1)
191
 
        G_free(var->v.p);
192
 
    var->v.p = NULL;
193
 
    freesym(var);
194
 
 
195
 
    printvec(pnt);
196
 
    pnt->v.p = NULL;
197
 
    freesym(pnt);
198
 
}
199
 
 
200
 
SYMBOL *mkpnt(double x, double y, double z)
201
 
{
202
 
    SYMBOL *pnt;
203
 
    VECTOR *vec;
204
 
 
205
 
    vec = (VECTOR *) listitem(sizeof(VECTOR));
206
 
    vec->x = x;
207
 
    vec->y = y;
208
 
    vec->z = z;
209
 
    vec->refcnt = 1;
210
 
 
211
 
    pnt = (SYMBOL *) listitem(sizeof(SYMBOL));
212
 
    pnt->type = pnt->itype = st_pnt;
213
 
    pnt->v.p = vec;
214
 
 
215
 
    return pnt;
216
 
}
217
 
 
218
 
SYMBOL *mkpntvar(SYMBOL * var, SYMBOL * pnt)
219
 
{
220
 
    var->type = var->itype = st_pnt;
221
 
    var->name = var->v.p;
222
 
    var->v.p = pnt->v.p;
223
 
    pnt->v.p = NULL;
224
 
    freesym(pnt);
225
 
 
226
 
    symtab = (SYMBOL *) listadd((LIST *) symtab, (LIST *) var, cmpsymsym);
227
 
 
228
 
    printvec(var);
229
 
 
230
 
    return var;
231
 
}
232
 
 
233
 
SYMBOL *pntfunc(SYMBOL * func, SYMBOL * arglist)
234
 
{
235
 
    SYMBOL *sym, *sptr;
236
 
    VECTOR *res = NULL;
237
 
    int argc = -1, i;
238
 
 
239
 
    sym = (SYMBOL *) listitem(sizeof(SYMBOL));
240
 
    sym->type = sym->itype = st_pnt;
241
 
 
242
 
    if (!func || !func->v.p || func->type != st_pfunc) {
243
 
        parseerror = 1;
244
 
        G_warning(_("Can't call bad function"));
245
 
    }
246
 
    else
247
 
        argc = listcnt((LIST *) arglist);
248
 
 
249
 
    for (i = 0, sptr = arglist; sptr; sptr = sptr->next, i++) {
250
 
        if (func->proto[i] == 'r')
251
 
            i++;
252
 
        if (func->proto[i] == 'p') {
253
 
            if (sptr->itype != st_pnt || !sptr->v.p)
254
 
                argc = -1;
255
 
        }
256
 
        else if (func->proto[i] == 'd' && sptr->itype != st_num)
257
 
            argc = -1;
258
 
    }
259
 
 
260
 
    res = (VECTOR *) listitem(sizeof(VECTOR));
261
 
 
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,
270
 
                                         arglist->next->v.d);
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,
275
 
                                         arglist->next->v.p);
276
 
    else if (argc == 3 && !strcmp(func->proto, "ppp"))
277
 
        res = (*(p_func_ppp) func->v.p) (arglist->v.p,
278
 
                                         arglist->next->v.p,
279
 
                                         arglist->next->next->v.p);
280
 
    else {
281
 
        G_warning(_("Bad arguments to pointfunc %s"), func->name);
282
 
        parseerror = 1;
283
 
        sym = (SYMBOL *) listdelall((LIST *) sym, freesym);
284
 
        if (res)
285
 
            G_free(res);
286
 
        return NULL;
287
 
    }
288
 
    sym->v.p = res;
289
 
 
290
 
    listdelall((LIST *) arglist, freesym);
291
 
 
292
 
    return sym;
293
 
}
294
 
 
295
 
SYMBOL *pntop(int op, SYMBOL * pnt1, SYMBOL * pnt2)
296
 
{
297
 
    SYMBOL *func, *arglist, *res = NULL;
298
 
    char buf[32];
299
 
 
300
 
    sprintf(buf, "pnt_op_func_%c", op);
301
 
 
302
 
    func = getsym(buf);
303
 
    if (!func) {
304
 
        G_warning(_("No function defined to perform ``point %c point''"), op);
305
 
        parseerror = 1;
306
 
    }
307
 
    else if (func->rettype == st_pnt) {
308
 
        res = (SYMBOL *) listitem(sizeof(SYMBOL));
309
 
        symcpy(res, func);
310
 
        res->next = NULL;
311
 
        func = res;
312
 
        arglist = (SYMBOL *) listapp(NULL, (LIST *) pnt1);
313
 
        arglist = (SYMBOL *) listapp((LIST *) arglist, (LIST *) pnt2);
314
 
 
315
 
        res = pntfunc(func, arglist);
316
 
    }
317
 
 
318
 
    return res;
319
 
}
320
 
 
321
 
SYMBOL *pntapp(SYMBOL * head, SYMBOL * elt)
322
 
{
323
 
    return (SYMBOL *) listapp((LIST *) head, (LIST *) elt);
324
 
}
325
 
 
326
 
/*
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.
330
 
 */
331
 
VECTOR *v_copy(VECTOR * p, VECTOR * p1)
332
 
{
333
 
    p->x = p1->x;
334
 
    p->y = p1->y;
335
 
    p->z = p1->z;
336
 
 
337
 
    return p;
338
 
}
339
 
 
340
 
/*
341
 
 * Vector addition
342
 
 * Result is 2D if at least one of p1 or p2 is 2D.
343
 
 */
344
 
VECTOR *v_add(VECTOR * p, VECTOR * p1, VECTOR * p2)
345
 
{
346
 
    p->x = p1->x + p2->x;
347
 
    p->y = p1->y + p2->y;
348
 
    /*
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
351
 
     */
352
 
    if (!isnan(p1->z) && !isnan(p2->z))
353
 
        p->z = p1->z + p2->z;
354
 
    else
355
 
        p->z = nanval;
356
 
 
357
 
    return p;
358
 
}
359
 
 
360
 
/*
361
 
 * Vector substraction
362
 
 * Result is 2D if at least one of p1 or p2 is 2D.
363
 
 */
364
 
VECTOR *v_sub(VECTOR * p, VECTOR * p1, VECTOR * p2)
365
 
{
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;
370
 
    else
371
 
        p->z = nanval;
372
 
 
373
 
    return p;
374
 
}
375
 
 
376
 
/*
377
 
 * Utility function to make all coordinates positive
378
 
 */
379
 
VECTOR *v_abs(VECTOR * p, VECTOR * p1)
380
 
{
381
 
    p->x = fabs(p1->x);
382
 
    p->y = fabs(p1->y);
383
 
    if (!isnan(p1->z))
384
 
        p->z = fabs(p1->z);
385
 
    else
386
 
        p->z = nanval;
387
 
 
388
 
    return p;
389
 
}
390
 
 
391
 
/*
392
 
 * Utility function to negate all coordinates
393
 
 */
394
 
VECTOR *v_neg(VECTOR * p, VECTOR * p1)
395
 
{
396
 
    p->x = -p1->x;
397
 
    p->y = -p1->y;
398
 
    if (!isnan(p1->z))
399
 
        p->z = -p1->z;
400
 
    else
401
 
        p->z = nanval;
402
 
 
403
 
    return p;
404
 
}
405
 
 
406
 
/*
407
 
 * Utility function to compare two doubles for equality without epsilon
408
 
 * This is not really true, as we consider NaN to be zero.
409
 
 */
410
 
static inline int _is_zero(double r)
411
 
{
412
 
    if ((isfinite(r) && r == 0.0) || isnan(r))
413
 
        return 1;
414
 
    return 0;
415
 
}
416
 
 
417
 
/*
418
 
 * Test for equality of two points. No epsion applied
419
 
 */
420
 
double v_eq(VECTOR * p1, VECTOR * p2)
421
 
{
422
 
    VECTOR p;
423
 
    int dim = 2;
424
 
 
425
 
    if (!isnan(p1->z) && !isnan(p2->z))
426
 
        dim = 3;
427
 
 
428
 
    v_sub(&p, p1, p2);
429
 
    v_abs(&p, &p);
430
 
    if (_is_zero(p.x) && _is_zero(p.y) && (dim == 2 || _is_zero(p.z)))
431
 
        return 1;
432
 
    return 0;
433
 
}
434
 
 
435
 
/*
436
 
 * Test for equality of two points by a given epsilon
437
 
 * epsilon is supposed to have positive values only.
438
 
 */
439
 
double v_eq_epsilon(VECTOR * p1, VECTOR * p2, VECTOR * e)
440
 
{
441
 
    VECTOR p;
442
 
    int dim = 2;
443
 
 
444
 
    if (!isnan(p1->z) && !isnan(p2->z))
445
 
        dim = 3;
446
 
 
447
 
    v_sub(&p, p1, p2);
448
 
    v_abs(&p, &p);
449
 
    if (p.x < e->x && p.y < e->y && (dim == 2 || (p.z < e->z)))
450
 
        return 1;
451
 
    return 0;
452
 
}
453
 
 
454
 
/*
455
 
 * Multiply a vector by a scalar
456
 
 */
457
 
VECTOR *v_mul(VECTOR * p, VECTOR * p1, double d)
458
 
{
459
 
    p->x = d * p1->x;
460
 
    p->y = d * p1->y;
461
 
    if (!isnan(p1->z))
462
 
        p->z = d * p1->z;
463
 
    else
464
 
        p->z = nanval;
465
 
 
466
 
    return p;
467
 
}
468
 
 
469
 
/*
470
 
 * Divide a vector by a scalar
471
 
 */
472
 
VECTOR *v_div(VECTOR * p, VECTOR * p1, double d)
473
 
{
474
 
    if (!isfinite(d) || d == 0.0) {
475
 
        parseerror = 1;
476
 
        return p;
477
 
    }
478
 
    p->x = p1->x / d;
479
 
    p->y = p1->y / d;
480
 
    if (!isnan(p1->z))
481
 
        p->z = p1->z / d;
482
 
    else
483
 
        p->z = nanval;
484
 
 
485
 
    return p;
486
 
}
487
 
 
488
 
/*
489
 
 * Compute the value of a vector
490
 
 */
491
 
double v_val(VECTOR * p)
492
 
{
493
 
    return sqrt(p->x * p->x + p->y * p->y +
494
 
                ((isnan(p->z)) ? 0 : p->z * p->z));
495
 
}
496
 
 
497
 
/*
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
500
 
 * origin.
501
 
 */
502
 
VECTOR *v_unit(VECTOR * p, VECTOR * p1)
503
 
{
504
 
    double val = v_val(p1);
505
 
 
506
 
    if (_is_zero(val))
507
 
        return v_copy(p, &pnt_o);
508
 
 
509
 
    p->x = p1->x / val;
510
 
    p->y = p1->y / val;
511
 
    if (!isnan(p1->z))
512
 
        p->z = p1->z / val;
513
 
    else
514
 
        p->z = nanval;
515
 
 
516
 
    return p;
517
 
}
518
 
 
519
 
/*
520
 
 * Compute the dot product of P1 and P2
521
 
 */
522
 
double v_dot(VECTOR * p1, VECTOR * p2)
523
 
{
524
 
    int dim = 2;
525
 
 
526
 
    if (!isnan(p1->z) && !isnan(p2->z))
527
 
        dim = 3;
528
 
    return p1->x * p2->x + p1->y * p2->y + ((dim == 2) ? 0 : p1->z * p2->z);
529
 
}
530
 
 
531
 
/*
532
 
 * Compute the cross product of P1 and P2
533
 
 * Return (0,0) for 2D
534
 
 */
535
 
 
536
 
VECTOR *v_cross(VECTOR * p, VECTOR * p1, VECTOR * p2)
537
 
{
538
 
    VECTOR p0;
539
 
 
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;
544
 
        v_copy(p, &p0);
545
 
    }
546
 
    else {
547
 
        p->x = p->y = 0;
548
 
        p->z = nanval;
549
 
    }
550
 
 
551
 
    return p;
552
 
}
553
 
 
554
 
/*
555
 
 * Decide if vector P1 is ortogonal to vector P2
556
 
 * Should test if either P1 or P2 are (0,0,0);
557
 
 */
558
 
double v_isortho(VECTOR * p1, VECTOR * p2)
559
 
{
560
 
    return v_dot(p1, p2) == 0;
561
 
}
562
 
 
563
 
/*
564
 
 * Decide if P1 and P2 are parallel. If they are but have a diferent
565
 
 * direction, -1 is returned.
566
 
 */
567
 
double v_ispara(VECTOR * p1, VECTOR * p2)
568
 
{
569
 
    double dot, val, dif;
570
 
 
571
 
    dot = v_dot(p1, p2);
572
 
    val = v_val(p1);
573
 
    val *= v_val(p2);
574
 
 
575
 
    dif = fabs(dot - val);
576
 
    if (_is_zero(dif))
577
 
        return 1;
578
 
 
579
 
    dif = fabs(dot + val);
580
 
    if (_is_zero(dif))
581
 
        return -1;
582
 
 
583
 
    return 0;
584
 
}
585
 
 
586
 
/*
587
 
 * Decide if P1 and P2 have and angle alpha which: 0 < alpha < 90.0
588
 
 */
589
 
double v_isacute(VECTOR * p1, VECTOR * p2)
590
 
{
591
 
    return v_dot(p1, p2) > 0;
592
 
}
593
 
 
594
 
/*
595
 
 * Return the area spanned by the two vectors P1 and P2
596
 
 * Works only in 3D.
597
 
 */
598
 
double v_area(VECTOR * p1, VECTOR * p2)
599
 
{
600
 
    VECTOR p;
601
 
 
602
 
    return 0.5 * v_val(v_cross(&p, p1, p2));
603
 
}