AUTO_TOPLEVELS := $(foreach m,$(USING_AUTOS),$(shell $(PLAY)/toplevel-find $m))
+symbolic.c: symbolic.py
+ ./symbolic.py -q >$@.tmp && mv -f $@.tmp $@
+
+findcurve.o: findcurve.c symbolic.c
+
autoincs: $(AUTO_INCS)
scads: $(addsuffix .auto.scad, $(AUTO_TOPLEVELS))
stls: $(addsuffix .auto.stl, $(AUTO_TOPLEVELS))
--- /dev/null
+/*
+ */
+
+#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]);
+
+static int cb_f(const gsl_vector *x, void *params, gsl_vector *f) {
+ double *target = params;
+ double F[6];
+
+ int i;
+ X_EXTRACT;
+ F_POPULATE;
+ for (i=0; i<N; i++) gsl_vector_set(f,i, F(i) - target[i]);
+ return 0;
+}
+
+static int cb_df(const gsl_vector *x, void *params, gsl_matrix *J) {
+ double J_COL[N];
+ int j;
+ X_EXTRACT;
+ J_POPULATE;
+}
+
+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);
+}
+
+static double *scan1double(void) {
+ double v;
+ int r;
+
+ r = scanf("%lf",&v);
+ if (ferror(stdin)) { perror("stdin"); exit(-1); }
+ if (feof(stdin)) exit(0);
+ if (r!=1) { fputs("bad input\n",stderr); exit(-1); }
+ return v;
+}
+
+static int main(void) {
+ double target[N], epsilon;
+
+ gsl_multiroot_function_fdf function;
+ function.f = cb_f;
+ function.df = cb_df;
+ function.fdf = cb_fdf;
+ function.n = N;
+ function.params = target;
+
+ gsl_vector *startpoint = gsl_vector_alloc(N);
+
+ gsl_multiroot_fdfsolver *solver
+ = gsl_multiroot_fdfsolver_alloc(gsl_multiroot_fdfsolver_hybrid, N);
+
+ for (;;) {
+ /* 2N+1 doubles: target, initial guess, epsilon for residual */
+ for (i=0; i<N; i++)
+ target[i] = scan1double();
+ for (i=0; i<N; i++)
+ gsl_vector_set(startpoint,i,scan1double());
+ epsilon = scan1double();
+
+ gsl_multiroot_fdfsolver_set(solver, function, startpoint);
+
+ for (;;) {
+ r = gsl_multiroot_fdfsolver_iterate(solver);
+ if (r) break;
+
+ r = gsl_multiroot_test_residual(solver->f, epsilon);
+ if (r != GSL_CONTINUE) break;
+ }
+
+ if (r) {
+ fprintf(stderr,"ERROR %s\n",gsl_strerror(gsl_errno));
+ exit(-1);
+ }
+ }
+}
dbg('j')
j = cse_prep_cprint(j, 'jtmp')
for ix in range(0, j.cols):
- cprint(ourccode(j.col(ix), 'J(%d)' % ix))
+ cprint(ourccode(j.col(ix), 'J_COL'))
cprint('J_END_COL(%d)' % ix)
else:
small = smalls[0]
)
cprint('} /* %s small */' % small)
+def gen_misc():
+ cprintraw('// AUTOGENERATED - DO NOT EDIT\n')
+ cprintraw('#define N %d\n' % len(params))
+
def gen_x_extract():
- cprint('#define X_EXTRACT \\')
+ cprint('#define X_EXTRACT')
for ix in range(0, len(params)):
- cprint('double %s = X(%d)' % (params[ix], ix))
+ cprint('double %s = X(%d);' % (params[ix], ix))
cprintraw()
def gen_f_populate():
gen_diff(result_dirnscaled, (sh*sh*la, th*th*la))
cprintraw('')
-gen_x_extract()
-gen_f_populate()
-gen_j_populate()
+def gen_C():
+ gen_misc()
+ gen_x_extract()
+ gen_f_populate()
+ gen_j_populate()
+
+gen_C()