chiark / gitweb /
wip new energy functions
[moebius2.git] / common.c
1 /*
2  * Generally useful stuff.
3  */
4
5 #include "common.h"
6
7 double magnD(const double pq[D3]) {
8   gsl_vector v;
9
10   v.size= D3;
11   v.stride= 1;
12   v.data= (double*)pq;
13   /* owner and block ought not to be used */
14
15   return gsl_blas_dnrm2(&v);
16 }
17
18 double magnD2(const double pq[D3]) {
19   int k;
20   double d2= 0;
21   
22   K d2= ffsqa(pq[k], d2);
23   return d2;
24 }
25
26 double hypotD(const double p[D3], const double q[D3]) {
27   int k;
28   double pq[D3];
29
30   K pq[k]= p[k] - q[k];
31   return magnD(pq);
32 }
33
34 double hypotD2(const double p[D3], const double q[D3]) {
35   int k;
36   double d2= 0;
37   
38   K d2= ffsqa(p[k] - q[k], d2);
39   return d2;
40 }
41
42 double hypotD2plus(const double p[D3], const double q[D3], double d2) {
43   int k;
44   K d2= ffsqa(p[k] - q[k], d2);
45   return d2;
46 }
47
48 void xprod(double r[D3], const double a[D3], const double b[D3]) {
49   r[0]= a[1]*b[2] - a[2]*b[1];
50   r[1]= a[2]*b[0] - a[0]*b[2];
51   r[2]= a[0]*b[1] - a[1]*b[0];
52 }
53
54 void xprod_norm(double r[D3], const double a[D3], const double b[D3],
55                 double absepsilon, double one) {
56   int k;
57   xprod(r,a,b);
58   double multby= one/(magnD(r) + absepsilon);
59   K r[k] *= multby;
60 }
61
62 double dotprod(const double a[D3], const double b[D3]) {
63   int k;
64   double result= 0;
65   K result += a[k] * b[k];
66   return result;
67 }
68
69 void libdie(const char *lib, int l, const char *str) {
70   fprintf(stderr,"%s library call failed, line %d: %s\n", lib, l, str);
71 }
72
73 void gsldie(int l, const char *what, int status) {
74   fprintf(stderr,"gsl function failed, line %d: %s: %s\n",
75           l, what, gsl_strerror(status));
76   exit(-1);
77 }
78
79 void diee(const char *what) { perror(what); exit(16); }
80 void fail(const char *emsg) { fputs(emsg,stderr); exit(12); }
81 void flushoutput(void) { if (fflush(stdout)||ferror(stdout)) diee("stdout"); }