Merge pull request #188 from taiga-project/186-release-quadratic-bezier

Add Bezier

Makefile

@@ -33,7 +33,7 @@ taiga_renate_fast.exe: src/main.cu | $(BIN)
 
 t: test
 
-test: $(OBJ)/tests.o  $(OBJ)/test_bspline.o  $(OBJ)/test_solver.o $(OBJ)/test_basic_functions.o $(OBJ)/basic_functions.o | $(BIN)
+test: $(OBJ)/tests.o  $(OBJ)/test_bspline.o  $(OBJ)/test_solver.o $(OBJ)/test_basic_functions.o $(OBJ)/basic_functions.o $(OBJ)/test_bezier.o | $(BIN)
 	$(GCC) $(TEST_FLAGS) -Isrc -Itests $^ -lm -o $(BIN)/test.exe
 
 $(OBJ)/%.o: tests/%.c $(OBJ)

include/core/maths/maths.cuh

@@ -5,4 +5,7 @@ __device__ double cross(double *u, double *v, int index);
 __device__ double interpolate(double y1, double y2, double x, double x1, double x2);
 __device__ double interpolate_from_vector(double *x_vector, double *y_vector, long length, double x_value);
 
+__device__ double solve_quadratic(double a, double b, double c);
+__device__ void interpolate_bezier(double X_prev[6], double X[6], double detector[4], double timestep);
+
 #endif //MATHS_CUH

include/interface/parameter_reader.h

@@ -7,5 +7,6 @@ int runnumber_reader(ShotProp *shot, RunProp *run);
 int parameter_reader(BeamProp *beam, ShotProp *shot, RunProp *run);
 void set_solver(RunProp *run, char* solver);
 void set_field_interpolation_method(RunProp *run, char* method);
+void set_detect_interpolation_method(RunProp *run, char* method);
 
 #endif //PARAMETER_READER_H

include/utils/prop.h

@@ -35,6 +35,7 @@ typedef struct RunProp{
     int solver;
     int init_source;
     int field_interpolation_method;
+    int detect_interpolation_method;
     bool is_electric_field_on;
     bool is_magnetic_field_perturbation;
     bool is_ionisation_on;
@@ -80,6 +81,7 @@ typedef struct TaigaCommons{
     double timestep;
     int solver;
     int field_interpolation_method;
+    int detect_interpolation_method;
     int *grid_size;
     double *spline_rgrid;
     double *spline_zgrid;

include/utils/taiga_constants.h

@@ -26,6 +26,9 @@
     #define CUBIC_SPLINE 3
     #define CUBIC_BSPLINE 13
 
+    #define LINEAR_INTERPOLATION 0
+    #define BEZIER_INTERPOLATION 1
+
     #define SPLINE_INDEX_ERROR -1
 
     #define SERVICE_VAR_LENGTH 10

src/core/detection.cu

@@ -1,6 +1,7 @@
 #include "core/maths/maths.cuh"
 
-__device__ int calculate_detection_position(double *X, double *X_prev, double *detector_geometry, double timestep){
+__device__ int calculate_detection_position(double *X, double *X_prev, double *detector_geometry, double timestep,
+                                            int method){
     double detector_R   = detector_geometry[0];
     double detector_Z   = detector_geometry[1];
     double detector_T   = detector_geometry[2];
@@ -19,7 +20,16 @@ __device__ int calculate_detection_position(double *X, double *X_prev, double *d
         (X_prev[2]-detector_T)*sin_beta;
 
     if((detector_distance*detector_distance_prev<=0) && (X[3]>0)){
-        for (int i=0; i<6; ++i) X[i] = interpolate(X_prev[i], X[i], 0, detector_distance_prev, detector_distance);
+        if (method == BEZIER_INTERPOLATION) {
+            double detector_plane[4] = {cos_alpha, sin_alpha * cos_beta, sin_beta,
+                                        -(cos_alpha * detector_R +
+                                        sin_alpha * cos_beta * detector_Z +
+                                        sin_beta * detector_T)};
+            interpolate_bezier(X_prev, X, detector_plane, timestep);
+        }else{
+            for (int i = 0; i < 6; ++i)
+                X[i] = interpolate(X_prev[i], X[i], 0, detector_distance_prev, detector_distance);
+        }
         X[TIME_OF_FLIGHT_ID] += interpolate(0, 1, 0, detector_distance_prev, detector_distance) * timestep;
         return CALCULATION_FINISHED;
     }

src/core/maths/maths.cu

@@ -20,4 +20,51 @@ __device__ double interpolate_from_vector(double *x_vector, double *y_vector, lo
     long i=0;
     for (i=0; (i<length-1) && (x_vector[i]>x_value); ++i);
     return y_vector[i+1] - (y_vector[i+1]-y_vector[i])*(x_value-x_vector[i])/(x_vector[i+1]-x_vector[i]);
+}
+
+// Bezier
+
+__device__ double solve_quadratic(double a, double b, double c) {
+    double D = b * b - 4.0 * a * c;
+    if (D >= 0) {
+        double sqrtD__2a = sqrt(D) / 2.0 / a;
+        double mb__2a = -b / 2.0 / a;
+        double t = mb__2a + sqrtD__2a;
+        if ( t>= 0 && t <= 1) {
+            return t;
+        }
+        t = mb__2a - sqrtD__2a;
+        if ( t>= 0 && t <= 1) {
+            return t;
+        }
+    }
+    return OUT_OF_RANGE;
+}
+
+__device__ void interpolate_bezier(double X_prev[6], double X[6], double D[4], double timestep) {
+    double P[3];
+
+    for (int i = 0 ; i < 3; ++i) {
+        P[i] = X_prev[i] + (1.0 / 3.0) * X_prev[i+3] * timestep;
+    }
+
+    double a = D[0] * (X_prev[0] - 2.0 * P[0] + X[0]) +
+               D[1] * (X_prev[1] - 2.0 * P[1] + X[1]) +
+               D[2] * (X_prev[2] - 2.0 * P[2] + X[2]);
+    double b = 2.0 * D[0] * (P[0] - X_prev[0]) +
+               2.0 * D[1] * (P[1] - X_prev[1]) +
+               2.0 * D[2] * (P[2] - X_prev[2]);
+    double c = D[0] * X_prev[0] + D[1] * X_prev[1] + D[2] * X_prev[2] + D[3];
+
+    double t = solve_quadratic(a, b, c);
+
+    if (t == OUT_OF_RANGE) {
+        return;
+    }
+
+    for (int i = 0 ; i < 3; ++i) {
+        X[i] = (1.0 - t) * (1.0 - t) * X_prev[i] +
+                2.0 * (1.0 - t) * t * P[i] +
+                t * t * X[i];
+    }
 }

src/core/traj.cu

@@ -91,7 +91,8 @@ __device__ int calculate_trajectory(TaigaCommons *c, double X[X_SIZE], int detce
 
         (*update_intensity)(psi_n, X, c, local_ts_index, local_ts_psi, timestep);
 
-        detcellid = calculate_detection_position(X, X_prev, c->detector_geometry, timestep);
+        detcellid = calculate_detection_position(X, X_prev, c->detector_geometry, timestep,
+                                                 c->detect_interpolation_method);
     }
 
     return detcellid;

src/init/device/init.cu

@@ -61,6 +61,7 @@ void init_device_structs(BeamProp beam, ShotProp shot, RunProp run, TaigaGlobals
     shared_common->timestep              = run.timestep;
     shared_common->solver                = run.solver;
     shared_common->field_interpolation_method = run.field_interpolation_method;
+    shared_common->detect_interpolation_method = run.detect_interpolation_method;
     shared_common->is_ionisation_on      = run.is_ionisation_on;
     shared_common->ionisation_energy     = get_ionisation_energy(beam.species, beam.charge);
 }

src/interface/parameter_reader.c

@@ -49,6 +49,9 @@ void set_taiga_parameter(char* par_name, char* par_value, BeamProp *beam, ShotPr
     else if (!strcmp(par_name, "field_interpolation_method")){
         set_field_interpolation_method(run, clean_string(par_value));
     }
+    else if (!strcmp(par_name, "detect_interpolation_method")){
+        set_detect_interpolation_method(run, clean_string(par_value));
+    }
     else if (!strcmp(par_name, "secondary_ionisation"))     run->is_ionisation_on = (bool)par_value_d;
     else if (!strcmp(par_name, "electric_field_value") || !strcmp(par_name, "matlab"))
         ;
@@ -125,3 +128,12 @@ void set_field_interpolation_method(RunProp *run, char* method){
         printf("Warning: Unvalid interpolation method: %s\nSet to bicubic spline [spline].", method);
     }
 }
+
+void set_detect_interpolation_method(RunProp *run, char* method){
+    string_to_lowercase(method);
+    if (!strcmp(method, "bezier")){
+        run->detect_interpolation_method = BEZIER_INTERPOLATION;
+    }else{
+        run->detect_interpolation_method = LINEAR_INTERPOLATION;
+    }
+}

src/utils/prop.c

@@ -31,6 +31,7 @@ void init_run_prop(RunProp *run){
     run->solver = SOLVER_RK45;
     run->init_source = READ_COORDINATES;
     run->field_interpolation_method = CUBIC_SPLINE;
+    run->detect_interpolation_method = LINEAR_INTERPOLATION;
     run->is_electric_field_on = false;
     run->is_magnetic_field_perturbation = false;
     run->is_ionisation_on = false;

tests/test_bezier.c

@@ -0,0 +1,32 @@
+#include "test_bezier.h"
+#include "taiga_test.h"
+#include "utils/taiga_constants.h"
+#define __device__ ;
+#include "core/maths/maths.cu"
+
+double get_bezier(double X_prev[6], double X[6], double detector[4], double timestep, int index){
+    interpolate_bezier(X_prev, X, detector, timestep);
+    return X[index];
+}
+
+int test_bezier() {
+    TAIGA_INIT_TEST("BEZIER");
+
+    double X_prev[6] = {-1, 1, 0, sqrt(2), sqrt(2), 0};
+    double X[6] = {2, 0, 0, 0, -1, 0};
+    double D[4] = {1, 0, 0, -1};
+    TAIGA_ASSERT_ALMOST_EQ(1, get_bezier(X_prev, X, D, 1, 0), "x @ x=1 vs x^2+y^2=4");
+
+    double X2_prev[6] = {0, sqrt(3), 0, sqrt(3), 1, 0};
+    double X2[6] = {2, sqrt(3), 0, sqrt(3), -1, 0};
+    TAIGA_ASSERT_ALMOST_EQ(1, get_bezier(X2_prev, X2, D, 1, 0), "x @ x=1 vs (x-1)^2+y^2=4");
+
+    double X3_prev[6] = {1, 2, 0, 1, 1, 0};
+    double X3[6] = {2, 1, 0, 1, -1, 0};
+    double D3[4] = {1.0, 1.0, 0, -1.0};
+    TAIGA_ASSERT_ALMOST_EQ(2, get_bezier(X3_prev, X3, D3, 1, 0), "x @ x=1 vs (x-1)^2+y^2=4");
+    TAIGA_ASSERT_ALMOST_EQ(1, get_bezier(X3_prev, X3, D3, 1, 1), "y @ x=1 vs (x-1)^2+y^2=4");
+
+
+    return TAIGA_ASSERT_SUMMARY();
+}

tests/test_bezier.h

@@ -0,0 +1,6 @@
+#ifndef TEST_INTERPOLATE_H
+#define TEST_INTERPOLATE_H
+
+int test_bezier();
+
+#endif //TEST_INTERPOLATE_H

tests/test_solver.c

@@ -10,7 +10,6 @@
 
 #include "utils/taiga_constants.h"
 #include "utils/prop.c"
-#include "core/maths/maths.cu"
 #include "core/solvers/rk4.cu"
 #include "core/solvers/runge_kutta_nystrom.cu"
 #include "core/solvers/yoshida.cu"

tests/tests.c

@@ -4,12 +4,13 @@
 #include "test_basic_functions.h"
 #include "test_solver.h"
 #include "test_bspline.h"
+#include "test_bezier.h"
 
 int main(){
-    return test_concat() +
+    return /*test_concat() +
     test_interpolate() +
     test_to_string() +
     test_bspline() +
-    test_solver();
+    test_solver() +*/
+    test_bezier();
 }
-