Skip to content

Commit 0e434c0

Browse files
simoneoteliesimoneotelie
committed
Analytical and numerical path length works. Assume we only use drift for analytical currently.
1 parent c7243e6 commit 0e434c0

File tree

4 files changed

+36
-20
lines changed

4 files changed

+36
-20
lines changed

‎xcoll/geometry/c_init/c_init.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,7 @@
4949
5050
#ifndef XC_GEOM_ROOT_GRID_POINTS
5151
#define XC_GEOM_ROOT_GRID_POINTS {XC_GEOM_ROOT_GRID_POINTS}
52+
#endif
5253
5354
#ifndef XC_GEOM_SIMPSON_SUBINTERVALS
5455
#define XC_GEOM_SIMPSON_SUBINTERVALS {XC_GEOM_SIMPSON_SUBINTERVALS}

‎xcoll/geometry/crossings/find_root.h‎

Lines changed: 23 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,15 @@
1010
#include <stdio.h>
1111

1212
/*gpufun*/
13-
void LineSegment_crossing_drift(FindRoot finder, LineSegment seg, double s0, double x0, double xm){
13+
void DriftTrajectory_analytical_solution_l(FindRoot finder, DriftTrajectory traj, double s, double x){
14+
if (DriftTrajectory_get_cos_t0(traj) - DriftTrajectory_get_cos_t0(traj) < XC_GEOM_EPSILON){
15+
FindRoot_set_solution_l(finder, FindRoot_get_num_solutions(finder), (x - DriftTrajectory_get_x0(traj)) / DriftTrajectory_get_sin_t0(traj));
16+
} else {
17+
FindRoot_set_solution_l(finder, FindRoot_get_num_solutions(finder), (s - DriftTrajectory_get_s0(traj)) / DriftTrajectory_get_cos_t0(traj));}
18+
}
19+
20+
/*gpufun*/
21+
void LineSegment_crossing_drift(FindRoot finder, LineSegment seg, DriftTrajectory traj, double s0, double x0, double xm){
1422
// Get segment data
1523
double s1 = LineSegment_get_s1(seg);
1624
double x1 = LineSegment_get_x1(seg);
@@ -24,9 +32,12 @@ void LineSegment_crossing_drift(FindRoot finder, LineSegment seg, double s0, dou
2432
// TODO: This is situational; we should return s1 if <get_s_first and current_s if after current_s
2533
// For now we hit twice (because we go nor IN nor OUT)
2634
FindRoot_set_solution_t(finder, FindRoot_get_num_solutions(finder), 0.0);
35+
DriftTrajectory_analytical_solution_l(finder, traj, LineSegment_func_s(seg, 0.0), LineSegment_func_x(seg, 0.0));
2736
FindRoot_set_converged(finder, FindRoot_get_num_solutions(finder), 1);
2837
FindRoot_set_num_solutions(finder, FindRoot_get_num_solutions(finder)+1);
38+
2939
FindRoot_set_solution_t(finder, FindRoot_get_num_solutions(finder), 1.0);
40+
DriftTrajectory_analytical_solution_l(finder, traj, LineSegment_func_s(seg, 1.0), LineSegment_func_x(seg, 1.0));
3041
FindRoot_set_converged(finder, FindRoot_get_num_solutions(finder), 1);
3142
FindRoot_set_num_solutions(finder, FindRoot_get_num_solutions(finder)+1);
3243
} else {
@@ -37,13 +48,14 @@ void LineSegment_crossing_drift(FindRoot finder, LineSegment seg, double s0, dou
3748
double t = (x0 - x1 - (s0 - s1)*xm) / denom;
3849
if (t >= 0 && t <= 1){
3950
FindRoot_set_solution_t(finder, FindRoot_get_num_solutions(finder), t);
51+
DriftTrajectory_analytical_solution_l(finder, traj, LineSegment_func_s(seg, t), LineSegment_func_x(seg, t));
4052
FindRoot_set_converged(finder, FindRoot_get_num_solutions(finder), 1);
4153
FindRoot_set_num_solutions(finder, FindRoot_get_num_solutions(finder)+1);
4254
}
4355
}
4456
}
4557
/*gpufun*/
46-
void HalfOpenLineSegment_crossing_drift(FindRoot finder, HalfOpenLineSegment seg, double s0, double x0, double xm){
58+
void HalfOpenLineSegment_crossing_drift(FindRoot finder, HalfOpenLineSegment seg, DriftTrajectory traj, double s0, double x0, double xm){
4759
// Get segment data
4860
double s1 = HalfOpenLineSegment_get_s1(seg);
4961
double x1 = HalfOpenLineSegment_get_x1(seg);
@@ -54,6 +66,7 @@ void HalfOpenLineSegment_crossing_drift(FindRoot finder, HalfOpenLineSegment seg
5466
// Trajectory is parallel to the segment
5567
if (fabs((x0 - x1)/(s0 - s1) - xm) < XC_GEOM_EPSILON){
5668
FindRoot_set_solution_t(finder, FindRoot_get_num_solutions(finder), 0.0);
69+
DriftTrajectory_analytical_solution_l(finder, traj, HalfOpenLineSegment_func_s(seg, 0.0), HalfOpenLineSegment_func_x(seg, 0.0));
5770
FindRoot_set_converged(finder, FindRoot_get_num_solutions(finder), 1);
5871
FindRoot_set_num_solutions(finder, FindRoot_get_num_solutions(finder)+1);
5972
// We do not add t=1 as it is a half-open segment
@@ -65,6 +78,7 @@ void HalfOpenLineSegment_crossing_drift(FindRoot finder, HalfOpenLineSegment seg
6578
double t = (x0 - x1 - (s0 - s1)*xm) / denom;
6679
if (t >= 0){ // We do not check for t<=1 as it is a half-open segment
6780
FindRoot_set_solution_t(finder, FindRoot_get_num_solutions(finder), t);
81+
DriftTrajectory_analytical_solution_l(finder, traj, HalfOpenLineSegment_func_s(seg, t), HalfOpenLineSegment_func_x(seg, t));
6882
FindRoot_set_converged(finder, FindRoot_get_num_solutions(finder), 1);
6983
FindRoot_set_num_solutions(finder, FindRoot_get_num_solutions(finder)+1);
7084
}
@@ -285,6 +299,8 @@ void slice_before_newton(FindRoot finder, LocalSegment seg, LocalTrajectory traj
285299
return;
286300
}
287301
num = FindRoot_get_num_solutions(finder);
302+
printf("Guess found at i,j: %d, %d (nest level %d). t = %f, l = %f\n", i, j, nest_level, t + 0.5 * t_step, l + 0.5 * l_step);
303+
fflush(stdout);
288304
FindRoot_set_guess_t(finder, num, t + 0.5 * t_step);
289305
FindRoot_set_guess_l(finder, num, l + 0.5 * l_step);
290306
FindRoot_set_num_solutions(finder, num + 1);
@@ -326,11 +342,11 @@ void FindRoot_find_crossing(FindRoot finder, LocalSegment seg, LocalTrajectory t
326342
double x0 = DriftTrajectory_get_x0((DriftTrajectory) LocalTrajectory_member(traj));
327343
switch (LocalSegment_typeid(seg)){
328344
case LocalSegment_LineSegment_t:
329-
LineSegment_crossing_drift(finder, (LineSegment) LocalSegment_member(seg), s0, x0, xp);
345+
LineSegment_crossing_drift(finder, (LineSegment) LocalSegment_member(seg), (DriftTrajectory) LocalTrajectory_member(traj), s0, x0, xp);
330346
return;
331347
break;
332348
case LocalSegment_HalfOpenLineSegment_t:
333-
return HalfOpenLineSegment_crossing_drift(finder, (HalfOpenLineSegment) LocalSegment_member(seg), s0, x0, xp);
349+
return HalfOpenLineSegment_crossing_drift(finder, (HalfOpenLineSegment) LocalSegment_member(seg), (DriftTrajectory) LocalTrajectory_member(traj), s0, x0, xp);
334350
break;
335351
case LocalSegment_BezierSegment_t:
336352
return BezierSegment_crossing_drift(finder, (BezierSegment) LocalSegment_member(seg), s0, x0, xp);
@@ -339,9 +355,9 @@ void FindRoot_find_crossing(FindRoot finder, LocalSegment seg, LocalTrajectory t
339355
return find_crossing_approximate(finder, seg, traj);
340356
}
341357
break;
342-
default:
343-
printf("Using approximate crossing method.\n");
344-
return find_crossing_approximate(finder, seg, traj);
358+
default:
359+
printf("Using approximate crossing method.\n");
360+
return find_crossing_approximate(finder, seg, traj);
345361
}
346362
}
347363
#endif /* XCOLL_GEOM_FIND_ROOT_H */

‎xcoll/geometry/crossings/find_root.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,7 @@ class FindRoot(xo.Struct):
4141
c_name="FindRoot_find_path_length",
4242
args=[
4343
xo.Arg(xo.ThisClass, name="finder"),
44+
xo.Arg(LocalSegment, name="seg"),
4445
xo.Arg(LocalTrajectory, name="traj"),
4546
], ret=None)}
4647
_needs_compilation = True

‎xcoll/geometry/crossings/path_length.h‎

Lines changed: 11 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -54,18 +54,16 @@ void find_path_length_analytic(FindRoot finder, LocalTrajectory traj){
5454
// only for drift for now
5555
}
5656
/*gpufun*/
57-
void FindRoot_find_path_length(FindRoot finder, LocalTrajectory traj){
58-
// TODO: we should find a find to classify if we are actually INSIDE or NOT. Drift outside dont need path length.
59-
switch (LocalTrajectory_typeid(traj)){
60-
case LocalTrajectory_DriftTrajectory_t:
61-
// TODO: Check with Frederik. I assume we always use the first solution as the second solution might not
62-
// even be relevant (say out - in again, then we actually never go in again. IN-OUT -> we change to mcs before out)
63-
//find_path_length_analytic(finder, traj);
64-
return;
65-
break;
66-
default:
67-
printf("Using approximate crossing method.\n");
68-
return simpson(finder, traj, XC_GEOM_SIMPSON_SUBINTERVALS);
69-
}
57+
void FindRoot_find_path_length(FindRoot finder, LocalSegment seg, LocalTrajectory traj){
58+
if (LocalSegment_typeid(seg) == LocalSegment_BezierSegment_t){
59+
return simpson(finder, traj, XC_GEOM_SIMPSON_SUBINTERVALS);
60+
}
61+
// TODO: Check with Frederik. I assume we always use the first solution as the second solution might not
62+
// even be relevant (say out - in again, then we actually never go in again. IN-OUT -> we change to mcs before out)
63+
if (LocalTrajectory_typeid(traj) == LocalTrajectory_DriftTrajectory_t){
64+
return find_path_length_analytic(finder, traj);
65+
} else {
66+
return simpson(finder, traj, XC_GEOM_SIMPSON_SUBINTERVALS);
67+
}
7068
}
7169
#endif /* XCOLL_GEOM_SIMPSON_H */

0 commit comments

Comments
 (0)