findcurve: compiles

[moebius3.git] / findcurve.c

diff --git a/findcurve.c b/findcurve.c
index 6ba51ef486f4b3819955b5c1af3686af1dd65e95..169116a826df4969cd9cfdab8c1ae9ad0c037829 100644 (file)
--- a/findcurve.c
+++ b/findcurve.c
@@ -1,13 +1,25 @@
 /*
  */
 
+#include <math.h>
+
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_sf.h>
 #include <gsl/gsl_multiroots.h>
 
 #include "symbolic.c"
 
 #define X(i) gsl_vector_get(x,i)
 #define J_END_COL(i) \
-  for (j=0; j<N; j++) gsl_matrix_set(f,i,j,J_COL[j]);
+  for (j=0; j<N; j++) gsl_matrix_set(J,i,j,J_COL[j]);
+
+static inline _Bool IS_SMALL(double v) {
+  return v < 1E6;
+}
+
+static inline double sinc(double x) {
+  return gsl_sf_sinc(x / M_PI);  
+}
 
 static int cb_f(const gsl_vector *x, void *params, gsl_vector *f) {
   double *target = params;
@@ -16,7 +28,7 @@ static int cb_f(const gsl_vector *x, void *params, gsl_vector *f) {
   int i;
   X_EXTRACT;
   F_POPULATE;
-  for (i=0; i<N; i++) gsl_vector_set(f,i, F(i) - target[i]);
+  for (i=0; i<N; i++) gsl_vector_set(f,i, F[i] - target[i]);
   return 0;
 }
 
@@ -25,15 +37,17 @@ static int cb_df(const gsl_vector *x, void *params, gsl_matrix *J) {
   int j;
   X_EXTRACT;
   J_POPULATE;
+  return 0;
 }
 
 static int cb_fdf(const gsl_vector *x, void *params,
                  gsl_vector *f, gsl_matrix *J) {
   cb_f(x,params,f);
-  cb_df(x,params,df);
+  cb_df(x,params,J);
+  return 0;
 }
 
-static double *scan1double(void) {
+static double scan1double(void) {
   double v;
   int r;
 
@@ -44,8 +58,9 @@ static double *scan1double(void) {
   return v;
 }
 
-static int main(void) {
+int main(void) {
   double target[N], epsilon;
+  int i;
 
   gsl_multiroot_function_fdf function;
   function.f = cb_f;
@@ -57,7 +72,7 @@ static int main(void) {
   gsl_vector *startpoint = gsl_vector_alloc(N);
 
   gsl_multiroot_fdfsolver *solver
-    = gsl_multiroot_fdfsolver_alloc(gsl_multiroot_fdfsolver_hybrid, N);
+    = gsl_multiroot_fdfsolver_alloc(gsl_multiroot_fdfsolver_hybridsj, N);
 
   for (;;) {
     /* 2N+1 doubles: target, initial guess, epsilon for residual */
@@ -67,8 +82,9 @@ static int main(void) {
       gsl_vector_set(startpoint,i,scan1double());
     epsilon = scan1double();
 
-    gsl_multiroot_fdfsolver_set(solver, function, startpoint);
+    gsl_multiroot_fdfsolver_set(solver, &function, startpoint);
 
+    int r;
     for (;;) {
       r = gsl_multiroot_fdfsolver_iterate(solver);
       if (r) break;
@@ -78,7 +94,7 @@ static int main(void) {
     }
 
     if (r) {
-      fprintf(stderr,"ERROR %s\n",gsl_strerror(gsl_errno));
+      fprintf(stderr,"ERROR %s\n",gsl_strerror(r));
       exit(-1);
     }
   }