1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
| struct vortex_filament{
int n_seg;
double a;
double* t;
coord* c;
coord* tvec;
coord* nvec;
coord* bvec;
coord pcar;
};
void draw_space_curve(int n, coord *p){
for (int i = 0; i < n-1; i++){
glBegin(GL_LINES);
glVertex3f(p[i ].x, p[i ].y, p[i ].z);
glVertex3f(p[i+1].x, p[i+1].y, p[i+1].z);
glEnd();
}
}
void fd_derivative( int n, double delta_t0, coord shift, coord *X, coord *dX){
for (int i = 1; i < n-1; i++){
foreach_dimension()
dX[i].x = (X[i+1].x - X[i-1].x)/(2*delta_t0);
}
foreach_dimension(){
dX[0].x = (X[1].x - X[n-2].x + shift.x)/(2*delta_t0);
dX[n-1].x = (X[1].x - X[n-2].x + shift.x)/(2*delta_t0);
}
}
#include <gsl/gsl_spline.h>
#pragma autolink -lgsl -lgslcblas
coord gsl_interp1d( int n, double* t0, coord * P0, double tq){
coord Pq;
gsl_interp_accel *acc = gsl_interp_accel_alloc();
foreach_dimension(){
double *P_x;
P_x = malloc(sizeof(double)*n);
for (int i = 0; i < n; i++)
P_x[i] = P0[i].x;
gsl_spline *spline_x = gsl_spline_alloc(gsl_interp_cspline, n);
gsl_spline_init(spline_x, t0, P_x, n);
Pq.x = gsl_spline_eval (spline_x, tq, acc);
gsl_spline_free (spline_x);
free(P_x);
}
gsl_interp_accel_free (acc);
return Pq;
}
double frenet_projection (double pos, void *params){
struct vortex_filament *p = (struct vortex_filament *) params;
int n_seg = p->n_seg;
double* t = p->t;
coord* c = p->c;
coord* tvec = p->tvec;
coord* nvec = p->nvec;
coord* bvec = p->bvec;
coord pcar = p->pcar;
coord ccar, frenet[3];
ccar = gsl_interp1d( n_seg, t, c, pos);
frenet[0] = gsl_interp1d( n_seg, t, tvec, pos);
frenet[1] = gsl_interp1d( n_seg, t, nvec, pos);
frenet[2] = gsl_interp1d( n_seg, t, bvec, pos);
return vecdot(vecdiff(pcar, ccar), frenet[0]);
}
#include <gsl/gsl_roots.h>
#include <gsl/gsl_errno.h>
double frenet_projection_min ( int n_seg, double as, double* t0, coord* c, coord* tvec, coord* nvec, coord* bvec, coord pcar, double r){
struct vortex_filament params = {n_seg, as, t0, c, tvec, nvec, bvec, pcar};
int status, verbose = 0;
int iter = 0, max_iter = 100;
const gsl_root_fsolver_type *T;
gsl_root_fsolver *s;
double x_lo = r - 0.25, x_hi = r + 0.25;
gsl_function F;
F.function = &frenet_projection;
F.params = ¶ms;
T = gsl_root_fsolver_brent;
s = gsl_root_fsolver_alloc (T);
gsl_set_error_handler_off();
gsl_root_fsolver_set (s, &F, x_lo, x_hi);
if (verbose == 1) {
printf ("using %s method\n", gsl_root_fsolver_name (s));
printf ("%5s [%9s, %9s] %9s %10s %9s\n", "iter", "lower", "upper", "root", "err", "err(est)");
}
do {
iter++;
status = gsl_root_fsolver_iterate (s);
r = gsl_root_fsolver_root (s);
x_lo = gsl_root_fsolver_x_lower (s);
x_hi = gsl_root_fsolver_x_upper (s);
status = gsl_root_test_interval (x_lo, x_hi, 0, 1e-8);
if ((status == GSL_SUCCESS) && (verbose == 1)){
printf ("Converged:\n");
printf ("%5d [%.7f, %.7f] %.7f %.7f\n", iter, x_lo, x_hi, r, x_hi - x_lo);
}
}
while (status == GSL_CONTINUE && iter < max_iter);
gsl_root_fsolver_free (s);
return r;
}
|
