Akifumi Takahashi / CubicSpline

Fork of TRP105F_Spline by Akifumi Takahashi

Files at this revision

API Documentation at this revision

Comitter:
aktk
Date:
Sun May 29 12:15:19 2016 +0000
Parent:
11:a279e31d8499
Child:
13:9a51747773af
Commit message:
value of y dont modified except in ival:0. In ival:1~_Sample_Num-2 y = getY(x) all has the value of the biggest value of y in ival:0.

Changed in this revision

CubicSpline.cpp Show annotated file Show diff for this revision Revisions of this file
--- a/CubicSpline.cpp	Sun May 29 09:31:44 2016 +0000
+++ b/CubicSpline.cpp	Sun May 29 12:15:19 2016 +0000
@@ -1,4 +1,9 @@
 #define DEBUG
+#define VERSION_C
+#define DEBUG_MAKE_MODEL
+//#define DEBUG_SOLVE
+#define DEBUG_GETX "DEBUG_GETX\n"
+//#define DEBUG_GETY "DEBUG_GETY\n"
 #include "CubicSpline.h"
 
 //  To get voltage of TRP105F
@@ -189,8 +194,6 @@
                    +pow(_Sample_Set[i].y - _Sample_Set[i-1].y, 2));
 }
 
-#define VERSION_C
-#define DEBUG_MAKE_MODEL
 //
 //  Function to define _u_spline, specific constants of spline.
 //
@@ -326,7 +329,6 @@
 //  Function to solve a cubic polinomial
 //  by using Gardano-Tartaglia formula
 //
-#define DEBUG_SOLVE
 void CubicSpline2d::_solve_cubic_f(
     std::complex<double>* arg_t,
     const double  arg_C[4],
@@ -378,8 +380,8 @@
         v = std::complex<double>(-q,-sqrt(-D));
         //u = pow(u, 1/3);
         //v = pow(v, 1/3);
-        u = std::exp(std::log(u)/3);
-        v = std::exp(std::log(u)/3);
+        u = std::exp(std::log(u)/std::complex<double>(3.0,0.0));
+        v = std::exp(std::log(u)/std::complex<double>(3.0,0.0));
     }
     //One real root and two complex root
     else {
@@ -411,7 +413,6 @@
 #endif
 }
 
-#define DEBUG_GETX "DEBUG_GETX\n"
 double  CubicSpline2d::getX(double arg_y)
 {
     double x;
@@ -450,15 +451,15 @@
                 } else*/ if (_Sample_Set[ival].t <= t_sol[i].real() && t_sol[i].real() < _Sample_Set[ival+1].t) {
                     t_real.push_back(t_sol[i].real());
                     t_ival.push_back(ival);
+#ifdef DEBUG_GETX
+                    g_Serial_Signal.printf("(t, i) = (%f, %d)\n", t_real[t_real.size() - 1], ival);
+#endif
                 }
-#ifdef DEBUG_GETX
-                g_Serial_Signal.printf("(t, i) = (%f, %d)\n", t_real[t_real.size() - 1], ival);
-#endif
             }
         }
     }
 
-    if(t_real.size > 0) {
+    if(t_real.size() > 0) {
         the_t = t_real[0];
         the_i = t_ival[0];
         //if t's size is bigger than 1
@@ -469,9 +470,12 @@
             }
         }
     } else {
+#ifdef DEBUG_GETX
+        g_Serial_Signal.printf("LastPoint\n");
+#endif
         the_t = _Last_Point.t;
-        for (int i = 0; i < _Sumple_Num - 1; i++)
-            if(_Sumple_Set[i].t <= the_t && the_t <= _Sample_Set[i+1])
+        for (int i = 0; i < _Sample_Num - 1; i++)
+            if(_Sample_Set[i].t <= the_t && the_t <= _Sample_Set[i+1].t)
                 the_i = i;
     }
     for(int i = 0; i < 4; i++) C[i] = _C_x[i][the_i];
@@ -479,13 +483,15 @@
 #ifdef DEBUG_GETX
     g_Serial_Signal.printf("(the_t, the_i):= (%f , %d)\n",the_t, the_i);
 #endif
-
+    _Last_Point = (Vxyt) {
+        x, arg_y, the_t
+    };
     return x;
 }
 
-#define DEBUG_GETY "DEBUG_GETY\n"
 double  CubicSpline2d::getY(double arg_x)
 {
+
     double y;
     double C[4];
     double the_t;
@@ -506,24 +512,31 @@
     for(int ival = 0; ival < _Sample_Num - 1; ival++) {
         for(int i = 0; i < 4; i++) C[i] = _C_x[i][ival];
         _solve_cubic_f(t_sol, C, arg_x);
+#ifdef DEBUG_GETY
+        g_Serial_Signal.printf("interval:%d \t %f < t < %f\n", ival, _Sample_Set[ival].t, _Sample_Set[ival+1].t);
+#endif
         for(int i = 0; i < 3; i++) {
             //  regarding only real solution
             //  acuracy (error range) is supposed +-10E-3 here(groundless)
             if(std::abs(t_sol[i].imag()) < 0.000001) {
-                /*  */ if (ival == 0               && t_sol[i].real() < _Sample_Set[ival].t) {
+                /*   if (ival == 0               && t_sol[i].real() < _Sample_Set[ival].t) {
                     t_real.push_back(_Sample_Set[ival].t);
                     t_ival.push_back(ival);
                 } else if (ival == _Sample_Num - 2 && _Sample_Set[ival + 1].t <= t_sol[i].real()) {
                     t_real.push_back(_Sample_Set[ival + 1].t);
                     t_ival.push_back(ival);
-                } else if (_Sample_Set[ival].t <= t_sol[i].real() && t_sol[i].real() < _Sample_Set[ival+1].t) {
+                } else*/ if (_Sample_Set[ival].t <= t_sol[i].real() && t_sol[i].real() < _Sample_Set[ival+1].t) {
                     t_real.push_back(t_sol[i].real());
                     t_ival.push_back(ival);
+#ifdef DEBUG_GETY
+                    g_Serial_Signal.printf("(t, i) = (%f, %d)\n", t_real[t_real.size() - 1], ival);
+#endif
                 }
             }
         }
+    }
 
-
+    if(t_real.size() > 0) {
         the_t = t_real[0];
         the_i = t_ival[0];
         //if t's size is bigger than 1
@@ -533,14 +546,23 @@
                 the_i = t_ival[i];
             }
         }
+    } else {
+#ifdef DEBUG_GETY
+        g_Serial_Signal.printf("LastPoint\n");
+#endif
+        the_t = _Last_Point.t;
+        for (int i = 0; i < _Sample_Num - 1; i++)
+            if(_Sample_Set[i].t <= the_t && the_t <= _Sample_Set[i+1].t)
+                the_i = i;
     }
-
     for(int i = 0; i < 4; i++) C[i] = _C_y[i][the_i];
     y = _cubic_f(the_t - _Sample_Set[the_i].t, C);
 #ifdef DEBUG_GETY
-    g_Serial_Signal.printf("(the_t, the_i):= (%f , %d)i\n",the_t, the_i);
+    g_Serial_Signal.printf("(the_t, the_i):= (%f , %d)\n",the_t, the_i);
 #endif
-
+    _Last_Point = (Vxyt) {
+        y, arg_x, the_t
+    };
     return y;
 }
 
@@ -558,6 +580,7 @@
         ft[i]= _Sample_Set[i].x;
     }
     _makeModel(t,ft,_C_x);
+    
     for(int i = 0; i < _Sample_Num; i++) {
         ft[i]= _Sample_Set[i].y;
     }
@@ -739,15 +762,19 @@
 void CubicSpline2d::printOutData()
 {
     FILE *fp;
+    double x,y;
     double d = (_Sample_Set[_Sample_Num - 1].x - _Sample_Set[0].x) / 100;
 
     fp = fopen("/local/log.txt", "w");  // open file in writing mode
 
-    fprintf(fp, "x, y\n");
-    for(int ival = 0; ival < _Sample_Num; ival++) {
-        for(double x = _Sample_Set[ival].x; x < _Sample_Set[ival+1].x; x += d) {
-            fprintf(fp, "%f,%f\n", x, getY(x));
+    fprintf(fp, "x, y, (t)\n");
+    for(int ival = 0; ival < _Sample_Num - 1; ival++) {
+        fprintf(fp, "ival: %d \n", ival);
+        for(x = _Sample_Set[ival].x; x < _Sample_Set[ival + 1].x; x += d) {
+            y = getY(x);
+            fprintf(fp, "%f,%f,(%f)\n", x, y, sqrt(x*x + y*y));
         }
+        fprintf(fp, "ival: %d \n", ival);
     }
 
     fprintf(fp, "\nSample:dst, vol\n");