better with more bendingness costs

[moebius2.git] / common.c

diff --git a/common.c b/common.c
index 146e5d54ef44a3f6fd9fc3f9879d1bc5d1fc7c34..bfe11b7dc23851e0e427f45533a86d9df2f0b112 100644 (file)
--- a/common.c
+++ b/common.c
@@ -4,15 +4,16 @@
 
 #include "common.h"
 
-double hypotD1(const double pq[D3]) {
-  gsl_vector v;
-
-  v.size= D3;
-  v.stride= 1;
-  v.vector.data= pq;
-  /* owner and block ought not to be used */
+double magnD(const double pq[D3]) {
+  return sqrt(magnD2(pq));
+}
 
-  return gsl_blas_snrm2(&v);
+double magnD2(const double pq[D3]) {
+  int k;
+  double d2= 0;
+  
+  K d2= ffsqa(pq[k], d2);
+  return d2;
 }
 
 double hypotD(const double p[D3], const double q[D3]) {
@@ -20,22 +21,58 @@ double hypotD(const double p[D3], const double q[D3]) {
   double pq[D3];
 
   K pq[k]= p[k] - q[k];
-  return hypotD1(pq);
+  return magnD(pq);
 }
 
 double hypotD2(const double p[D3], const double q[D3]) {
+  int k;
   double d2= 0;
+  
   K d2= ffsqa(p[k] - q[k], d2);
   return d2;
 }
 
 double hypotD2plus(const double p[D3], const double q[D3], double d2) {
+  int k;
   K d2= ffsqa(p[k] - q[k], d2);
   return d2;
 }
 
 void xprod(double r[D3], const double a[D3], const double b[D3]) {
-  r[0]= a[1]*b[2] - a[2]*b[1]);
-  r[1]= a[2]*b[0] - a[0]*b[2]);
-  r[2]= a[0]*b[1] - a[1]*b[0]);
+  r[0]= a[1]*b[2] - a[2]*b[1];
+  r[1]= a[2]*b[0] - a[0]*b[2];
+  r[2]= a[0]*b[1] - a[1]*b[0];
+}
+
+void normalise(double v[D3], double one, double absepsilon) {
+  int k;
+  double multby= one/(magnD(v) + absepsilon);
+  K v[k] *= multby;
 }
+
+void xprod_norm(double r[D3], const double a[D3], const double b[D3],
+               double one, double absepsilon) {
+  xprod(r,a,b);
+  normalise(r, absepsilon, one);
+}
+
+double dotprod(const double a[D3], const double b[D3]) {
+  int k;
+  double result= 0;
+  K result += a[k] * b[k];
+  return result;
+}
+
+void libdie(const char *lib, int l, const char *str) {
+  fprintf(stderr,"%s library call failed, line %d: %s\n", lib, l, str);
+}
+
+void gsldie(int l, const char *what, int status) {
+  fprintf(stderr,"gsl function failed, line %d: %s: %s\n",
+         l, what, gsl_strerror(status));
+  exit(-1);
+}
+
+void diee(const char *what) { perror(what); exit(16); }
+void fail(const char *emsg) { fputs(emsg,stderr); exit(12); }
+void flushoutput(void) { if (ferror(stdout)||fflush(stdout)) diee("stdout"); }