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 normalise(double v[D3], double one, double absepsilon) {
  55   int k;
  56   double multby= one/(magnD(v) + absepsilon);
  57   K v[k] *= multby;
  58 }
  59
  60 void xprod_norm(double r[D3], const double a[D3], const double b[D3],
  61                 double one, double absepsilon) {
  62   xprod(r,a,b);
  63   normalise(r, absepsilon, one);
  64 }
  65
  66 double dotprod(const double a[D3], const double b[D3]) {
  67   int k;
  68   double result= 0;
  69   K result += a[k] * b[k];
  70   return result;
  71 }
  72
  73 void libdie(const char *lib, int l, const char *str) {
  74   fprintf(stderr,"%s library call failed, line %d: %s\n", lib, l, str);
  75 }
  76
  77 void gsldie(int l, const char *what, int status) {
  78   fprintf(stderr,"gsl function failed, line %d: %s: %s\n",
  79           l, what, gsl_strerror(status));
  80   exit(-1);
  81 }
  82
  83 void diee(const char *what) { perror(what); exit(16); }
  84 void fail(const char *emsg) { fputs(emsg,stderr); exit(12); }
  85 void flushoutput(void) { if (fflush(stdout)||ferror(stdout)) diee("stdout"); }