Akifumi Takahashi
This lib is considered to be used as a sensor's calibration program. Calibration with Spline Interpolation might be useful in the case that you want some model expressing relationship such like between a value of physical quantity and your sensor's voltage, but you cannot estimate a model such as liner, square, cubic polynomial, or sine curve. This makes (Parametric) Cubic Spline Polynomial Model (Coefficients of the polynomial) from some sample plots(e.g. sets of (value, voltage)). The inverse function (x,y)->(y,x) has been implemented so as to get analog data (not stepping or leveled data).
Fork of
TRP105F_Spline
by 
Akifumi Takahashi
Diff: CubicSpline.cpp
- Revision:
- 9:1903c6f8d5a9
- Child:
- 10:607a68db6303
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/CubicSpline.cpp Sun May 29 09:14:54 2016 +0000
@@ -0,0 +1,786 @@
+#define DEBUG
+#include "CubicSpline.h"
+
+// To get voltage of TRP105F
+AnalogIn g_Sensor_Voltage(p16);
+// To get sample distance via seral com
+Serial g_Serial_Signal(USBTX, USBRX);
+
+LocalFileSystem local("local"); // マウントポイントを定義（ディレクトリパスになる）
+// for debug
+#ifdef DEBUG
+DigitalOut led1(LED1);
+DigitalOut led2(LED2);
+DigitalOut led3(LED3);
+DigitalOut led4(LED4);
+#endif
+
+CubicSpline2d::CubicSpline2d()
+ :_useType(AsMODULE)
+{
+ _Sample_Num = 5;
+ _Sample_Set = (Vxyt *)malloc(_Sample_Num * sizeof(Vxyt));
+ _Last_Point = (Vxyt) {
+ 0,0,0
+ };
+ for(int i = 0; i < 4; i++) {
+ _C_x[i]= (double*)malloc((_Sample_Num - 1)* sizeof(double));;
+ _C_y[i]= (double*)malloc((_Sample_Num - 1)* sizeof(double));;
+ }
+ //calibrateSensor();
+}
+
+CubicSpline2d::CubicSpline2d(
+ unsigned int arg_num
+)
+ :_useType(AsMODULE)
+{
+ _Sample_Num = arg_num;
+ _Sample_Set = (Vxyt *)malloc(_Sample_Num * sizeof(Vxyt));
+ _Last_Point = (Vxyt) {
+ 0,0,0
+ };
+ for(int i = 0; i < 4; i++) {
+ _C_x[i]= (double*)malloc((_Sample_Num - 1)* sizeof(double));;
+ _C_y[i]= (double*)malloc((_Sample_Num - 1)* sizeof(double));;
+ }
+ //calibrateSensor();
+}
+
+CubicSpline2d::CubicSpline2d(
+ unsigned int arg_num,
+ UseType arg_useType
+)
+ :_useType(arg_useType)
+{
+ _Sample_Num = arg_num;
+ _Sample_Set = (Vxyt *)malloc(_Sample_Num * sizeof(Vxyt));
+ _Last_Point = (Vxyt) {
+ 0,0,0
+ };
+ for(int i = 0; i < 4; i++) {
+ _C_x[i]= (double*)malloc((_Sample_Num - 1)* sizeof(double));;
+ _C_y[i]= (double*)malloc((_Sample_Num - 1)* sizeof(double));;
+ }
+ //calibrateSensor();
+}
+
+CubicSpline2d::~CubicSpline2d()
+{
+ free(_Sample_Set);
+ //free(_u_param);
+ for(int i = 0; i < 4; i++) {
+ free(_C_x[i]);
+ free(_C_y[i]);
+ }
+}
+
+void CubicSpline2d::_sampleData()
+{
+ int tmp;
+ char sig;
+ Vxyt tmp_set;
+
+ int floatflag = 0;
+
+ // For evry set,
+ // 1, get dst data via serai com,
+ // 2, get vol data,
+ // and then do same for next index set.
+ for(int i = 0; i < _Sample_Num; i++) {
+ if(_useType == AsDEBUG) {
+ //
+ // Recieve a Distance datus and store it into member
+ //
+ g_Serial_Signal.printf("X:");
+ _Sample_Set[i].x = 0;
+ do {
+ sig = g_Serial_Signal.getc();
+ if('0' <= sig && sig <= '9') {
+ if(floatflag == 0) {
+ _Sample_Set[i].x = 10 * _Sample_Set[i].x + sig - 48;
+ } else {
+ _Sample_Set[i].x = _Sample_Set[i].x + (sig - 48) * pow((double)10, (double)- floatflag);
+ floatflag++;
+ }
+ g_Serial_Signal.putc(char(sig));
+ } else if(sig == '.') {
+ if(floatflag == 0) {
+ floatflag = 1;
+ g_Serial_Signal.putc(char(sig));
+ }
+ } else if(sig == 0x08) {
+ _Sample_Set[i].x = 0;
+ g_Serial_Signal.printf("[canseled!]");
+ g_Serial_Signal.putc('\n');
+ g_Serial_Signal.putc('>');
+ }
+ } while (!(sig == 0x0a || sig == 0x0d));
+ floatflag = 0;
+ g_Serial_Signal.putc('\n');
+ g_Serial_Signal.printf("x:%f|",_Sample_Set[i].x);
+ //
+ // Recieve a Voltage datus and store it into member
+ //
+ // LOW PASS FILTERED
+ // Get 10 data and store mean as a sample.
+ // After get one original sample, system waits for 0.1 sec,
+ // thus it takes 1 sec evry sampling.
+ _Sample_Set[i].y = 0;
+ for(int j = 0; j < 10; j++) {
+ tmp_set.y = g_Sensor_Voltage.read();
+#ifdef DEBUG
+ g_Serial_Signal.printf("%f,",tmp_set.y);
+#endif
+ _Sample_Set[i].y += (tmp_set.y / 10);
+ wait(0.1);
+ }
+#ifdef DEBUG
+ g_Serial_Signal.printf("(%f)\n",_Sample_Set[i].y);
+#endif
+ }
+
+ // if the input data is over the bound, it is calibrated
+ if (_Sample_Set[i].x < 0)
+ _Sample_Set[i].x = 0;
+ }
+ //
+ // Sort set data array in x-Ascending order
+ //
+ tmp = 0;
+ for( int i = 0 ; i < _Sample_Num; i++) {
+ for(int j = _Sample_Num - 1; i < j ; j--) {
+ // use dst as index for dst range [2,20]
+ if (_Sample_Set[i].x > _Sample_Set[j].x) {
+ tmp_set.x = _Sample_Set[i].x;
+ tmp_set.y = _Sample_Set[i].y;
+ _Sample_Set[i].x = _Sample_Set[j].x;
+ _Sample_Set[i].y = _Sample_Set[j].y;
+ _Sample_Set[j].x = tmp_set.x;
+ _Sample_Set[j].y = tmp_set.y;
+ }
+ // if a same dst has been input, calcurate mean.
+ else if (_Sample_Set[i].x == _Sample_Set[j].x) {
+ tmp_set.y = (_Sample_Set[i].y + _Sample_Set[j].y)/2;
+ _Sample_Set[i].y = tmp_set.y;
+ for(int k = j; k < _Sample_Num - 1; k++)
+ _Sample_Set[k] = _Sample_Set[k+1];
+ tmp++;
+ }
+ }
+ }
+#ifdef DEBUG
+ g_Serial_Signal.printf(" _Sample_num: %d\n", _Sample_Num );
+ g_Serial_Signal.printf("-) tmp: %d\n", tmp );
+#endif
+ // substruct tmp from number of sample.
+ _Sample_Num -= tmp;
+#ifdef DEBUG
+ g_Serial_Signal.printf("-----------------\n");
+ g_Serial_Signal.printf(" _Sample_num: %d\n", _Sample_Num );
+#endif
+
+ // generate t which is parameter related to x,y
+ _Sample_Set[0].t = 0;
+ for(int i = 1; i < _Sample_Num; i++)
+ _Sample_Set[i].t =
+ _Sample_Set[i-1].t
+ + sqrt(pow(_Sample_Set[i].x - _Sample_Set[i-1].x, 2)
+ +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.
+//
+void CubicSpline2d::_makeModel(const double* arg_sampled_t, const double* arg_sampled_ft, double* arg_C[4], const unsigned int arg_num)
+{
+ // arg_t : t; The variable of f(t)
+ // arg_ft: f(t); The cubic poliminal in Interval-j.
+ // arg_C[i]: Ci; The coefficient of t^i of f(t) that defines Spline Model Poliminal f(t).
+ // arg_num: j in [0,_Sample_Num-1]; The number of interval.
+ // f(t)j = C3j*t^3 + C2j*t^2 + C1j*t + C0j
+ //
+ // N: max of index <=> (_Sample_Num - 1)
+ //
+ // u[i] === d^2/dx^2(Spline f)[i]
+ // i:[0,N]
+ // u[0] = u[N] = 0
+#if defined (VERSION_C)
+ double *u = (double*)malloc((arg_num ) * sizeof(double));
+#elif defined (VERSION_Cpp)
+ double *u = new double[arg_num];
+#elif defined (VERSION_Cpp11)
+ std::array<double,arg_num> u;
+#endif
+ //
+ // h[i] = x[i+1] - x[i]
+ // i:[0,N-1]; num of elm: N<=>_Sample_Num - 1
+ double *h = (double*)malloc((arg_num - 1) * sizeof(double));
+ //
+ // v[i] = 6*((y[i+2]-y[i+1])/h[i+1] + (y[i+1]-y[i])/h[i])
+ // i:[0,N-2]
+ double *v = (double*)malloc((arg_num - 2) * sizeof(double));
+ //
+ // temporary array whose num of elm equals v array
+ double *w = (double*)malloc((arg_num - 2) * sizeof(double));
+ //
+ // [ 2(h[0]+h[1]) , h[1] , O ] [u[1] ] [v[0] ]
+ // [ h[1] , 2(h[1]+h[2]) , h[2] ] [u[2] ] [v[1] ]
+ // [ ... ] * [... ] = [... ]
+ // [ h[j] , 2(h[j]+h[j+1]) , h[j+1] ] [u[j+1]] [v[j] ]
+ // [ ... ] [ ... ] [ ... ]
+ // [ h[N-3] , 2(h[N-3]+h[N-2]), h[N-2] ] [u[j+1]] [v[j] ]
+ // [ O h[N-2] , 2(h[N-2]+h[N-1]) ] [u[N-1]] [v[N-2]]
+ //
+ // For LU decomposition
+ double *Upper = (double*)malloc((arg_num - 2) * sizeof(double));
+ double *Lower = (double*)malloc((arg_num - 2) * sizeof(double));
+#ifdef DEBUG_MAKE_MODEL
+ _printOutDataCouple(arg_sampled_t, arg_sampled_ft, arg_num, "\nargment set\n");
+#endif
+
+ for(int i = 0; i < arg_num - 1; i++)
+ h[i] = (double)(arg_sampled_t[i + 1] - arg_sampled_t[i]);
+
+ for(int i = 0; i < arg_num - 2; i++)
+ v[i] = 6 * (
+ ((double)(arg_sampled_ft[i + 2] - arg_sampled_ft[i + 1])) / h[i + 1]
+ -
+ ((double)(arg_sampled_ft[i + 1] - arg_sampled_ft[i])) / h[i]
+ );
+
+ //
+ // LU decomposition
+ //
+ Upper[0] = 2 * (h[0] + h[1]);
+ Lower[0] = 0;
+ for (int i = 1; i < arg_num - 2; i++) {
+ Lower[i] = h[i] / Upper[i - 1];
+ Upper[i] = 2 * (h[i] + h[i + 1]) - Lower[i] * h[i];
+ }
+ //
+ // forward substitution
+ //
+ w[0] = v[0];
+ for (int i = 1; i < arg_num - 2; i ++) {
+ w[i] = v[i] - Lower[i] * w[i-1];
+ }
+ //
+ // backward substitution
+ //
+ u[arg_num - 2] = w[arg_num - 3] / Upper[arg_num - 3];
+ for(int i = arg_num - 3; i > 0; i--) {
+ u[i] = (w[(i - 1)] - h[(i)] * u[(i) + 1]) / Upper[(i - 1)];
+ }
+ // _u_spline[i] === d^2/dx^2(Spline f)[i]
+ u[0] = u[arg_num - 1] = 0.0;
+#ifdef DEBUG_MAKE_MODEL
+ _printOutData(h, arg_num - 1, "h");
+ _printOutData(v, arg_num - 2, "v");
+ _printOutData(w, arg_num - 2, "w");
+ _printOutData(Upper, arg_num - 2, "Upper");
+ _printOutData(Lower, arg_num - 2, "Lower");
+ _printOutData(u, arg_num , "u");
+#endif
+
+ for(int ival = 0; ival < arg_num - 1; ival++) {
+ arg_C[3][ival] = (u[ival + 1] - u[ival]) / 6.0 / (arg_sampled_t[ival + 1] - arg_sampled_t[ival]);
+ arg_C[2][ival] = (u[ival]) / 2.0;
+ arg_C[1][ival] = (arg_sampled_ft[ival + 1] - arg_sampled_ft[ival]) / (arg_sampled_t[ival + 1] - arg_sampled_t[ival])
+ -
+ (arg_sampled_t[ival + 1] - arg_sampled_t[ival]) * (u[ival + 1] + 2.0 * u[ival]) / 6.0;
+ arg_C[0][ival] = (arg_sampled_ft[ival]);
+ }
+#ifdef DEBUG_MAKE_MODEL
+ for(int ival = 0; ival < arg_num - 1; ival++) {
+ for(int i = 0; i < 4; i++)
+ g_Serial_Signal.printf("C[%d][%d]: %f\n", i, ival, arg_C[i][ival]);
+ }
+#endif
+
+ free(h);
+ free(u);
+ free(v);
+ free(w);
+ free(Upper);
+ free(Lower);
+}
+void CubicSpline2d::_makeModel(const double* arg_t, const double* arg_ft, double* arg_C[4])
+{
+ _makeModel(arg_t, arg_ft, arg_C, _Sample_Num);
+}
+//
+// Fuction to return the value of Cubic polynomial f(t)
+//
+double CubicSpline2d::_cubic_f(const double arg_t, const double arg_C[4])
+{
+ double ft; //the value of Spline f(t).
+
+ ft = arg_C[3] * pow(arg_t, 3) + arg_C[2] * pow(arg_t, 2) + arg_C[1] * arg_t + arg_C[0];
+
+ return ft;
+}
+//
+// 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],
+ const double arg_ft)
+{
+#ifdef DEBUG_SOLVE
+ for(int i = 0; i < 4; i++)
+ g_Serial_Signal.printf("C%d: %f\n", i, arg_C[i]);
+#endif
+
+ double c[3];
+ //f(t) = arg_ft/arg_C[3]
+ // = t^3 + c[2]*t^2 + c[1]*t + c[0].
+ for(int i = 0; i < 3; i++) {
+ c[i] = arg_C[i] / arg_C[3];
+ }
+ //modify the formula
+ //t^3 + c[2]*t^2 + c[1]*t + (c[0] - ft) = 0.
+ c[0] -= arg_ft / arg_C[3];
+#ifdef DEBUG_SOLVE
+ for(int i = 0; i < 3; i++)
+ g_Serial_Signal.printf("c%d: %f\n", i, c[i]);
+#endif
+
+ //The values defined from coefficients of the formula
+ //that identify solutions
+ double p,q,d;
+ p = ( -pow(c[2], 2) + 3 * c[1]) / 9;
+ q = (2 * pow(c[2], 3) - 9 * c[2] * c[1] + 27 * c[0]) / 54;
+ d = - c[2] / 3;
+
+ //Discriminant section
+ double D;
+ D = pow(p, 3) + pow(q, 2);
+#ifdef DEBUG_SOLVE
+ g_Serial_Signal.printf("p: %f\n", p);
+ g_Serial_Signal.printf("q: %f\n", q);
+ g_Serial_Signal.printf("d: %f\n", d);
+ g_Serial_Signal.printf("D: %f\n", D);
+#endif
+
+ //The values defined from p and q
+ //that idetify solutions
+ std::complex<double> u,v;
+
+ //Real root only
+ if(D <= 0) {
+ u = std::complex<double>(-q, sqrt(-D));
+ v = std::complex<double>(-q,-sqrt(-D));
+ u = pow(u, 1/3);
+ v = pow(v, 1/3);
+ }
+ //One real root and two complex root
+ else {
+ u = std::complex<double>(-q+sqrt(D),0.0);
+ v = std::complex<double>(-q-sqrt(D),0.0);
+ u = std::complex<double>(cbrt(u.real()), 0.0);
+ v = std::complex<double>(cbrt(v.real()), 0.0);
+ }
+#ifdef DEBUG_SOLVE
+ g_Serial_Signal.printf("u: %f + (%f)i\n", u.real(), u.imag());
+ g_Serial_Signal.printf("v: %f + (%f)i\n", v.real(), v.imag());
+#endif
+
+ //Cubic root of 1
+ std::complex<double> omega[3]= {
+ std::complex<double>( 1.0, 0.0),
+ std::complex<double>(-1/2, sqrt(3.0)/2),
+ std::complex<double>(-1/2,-sqrt(3.0)/2)
+ };
+
+ //Solution of the formula
+ arg_t[0] = omega[0] * u + omega[0] * v + d;
+ arg_t[1] = omega[1] * u + omega[2] * v + d;
+ arg_t[2] = omega[2] * u + omega[1] * v + d;
+
+#ifdef DEBUG_SOLVE
+ for(int i = 0; i < 3; i++)
+ g_Serial_Signal.printf("t%d: %f + (%f)i\n", i, arg_t[i].real(), arg_t[i].imag() );
+#endif
+}
+
+#define DEBUG_GETX "DEBUG_GETX\n"
+double CubicSpline2d::getX(double arg_y)
+{
+ double x;
+ double C[4];
+ double the_t;
+ int the_i;
+ std::complex<double>t_sol[3];
+ std::vector<double> t_real;
+ std::vector<int> t_ival;
+
+#ifdef DEBUG_GETX
+ g_Serial_Signal.printf(DEBUG_GETX);
+#endif
+ // For the every Intervals of Spline,
+ //it solves the polynomial defined by C[i] of the interval,
+ //checks the solutions are real number,
+ //and ckecks the solutions are in the interval.
+ // And if not-excluded solutions are more than one,
+ //it trys to find which one is more nearest to last point.
+ for(int ival = 0; ival < _Sample_Num - 1; ival++) {
+ for(int i = 0; i < 4; i++) C[i] = _C_y[i][ival];
+ _solve_cubic_f(t_sol, C, arg_y);
+#ifdef DEBUG_GETX
+ 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) {
+ 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) {
+ 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
+ }
+ }
+ }
+
+ the_t = t_real[0];
+ the_i = t_ival[0];
+ //if t's size is bigger than 1
+ for(int i = 1; i < t_real.size(); i++) {
+ if(std::abs(t_real[i] - _Last_Point.t) < std::abs(the_t - _Last_Point.t)) {
+ the_t = t_real[i];
+ the_i = t_ival[i];
+ }
+ }
+ for(int i = 0; i < 4; i++) C[i] = _C_x[i][the_i];
+ x = _cubic_f(the_t - _Sample_Set[the_i].t, C);
+#ifdef DEBUG_GETX
+ g_Serial_Signal.printf("(the_t, the_i):= (%f , %d)\n",the_t, the_i);
+#endif
+
+ return x;
+}
+
+#define DEBUG_GETY "DEBUG_GETY\n"
+double CubicSpline2d::getY(double arg_x)
+{
+ double y;
+ double C[4];
+ double the_t;
+ int the_i;
+ std::complex<double>t_sol[3];
+ std::vector<double> t_real;
+ std::vector<int> t_ival;
+
+#ifdef DEBUG_GETY
+ g_Serial_Signal.printf(DEBUG_GETY);
+#endif
+ // For the every Intervals of Spline,
+ //it solves the polynomial defined by C[i] of the interval,
+ //checks the solutions are real number,
+ //and ckecks the solutions are in the interval.
+ // And if not-excluded solutions are more than one,
+ //it trys to find which one is more nearest to last point.
+ 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);
+ 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) {
+ 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) {
+ t_real.push_back(t_sol[i].real());
+ t_ival.push_back(ival);
+ }
+ }
+ }
+
+
+ the_t = t_real[0];
+ the_i = t_ival[0];
+ //if t's size is bigger than 1
+ for(int i = 1; i < t_real.size(); i++) {
+ if(std::abs(t_real[i] - _Last_Point.t) < std::abs(the_t - _Last_Point.t)) {
+ the_t = t_real[i];
+ the_i = t_ival[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);
+#endif
+
+ return y;
+}
+
+
+void CubicSpline2d::calibrateSensor()
+{
+ double t[_Sample_Num];
+ double ft[_Sample_Num];
+
+ _sampleData();
+ _Last_Point = _Sample_Set[0];
+
+ for(int i = 0; i < _Sample_Num; i++) {
+ t[i] = _Sample_Set[i].t;
+ 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;
+ }
+ _makeModel(t,ft,_C_y);
+
+}
+
+void CubicSpline2d::saveSetting()
+{
+ FILE *fp;
+
+ fp = fopen("/local/savedata.log", "wb");
+
+ // Save _Sample_Num
+ fwrite(&_Sample_Num, sizeof(unsigned int), 1, fp);
+ fputc(0x3b, fp);
+ // Save _Sample_Set
+ for(int i = 0; i < _Sample_Num; i++) {
+ fwrite(&_Sample_Set[i].x, sizeof(double), 1, fp);
+ fputc(0x2c, fp);
+ fwrite(&_Sample_Set[i].y, sizeof(double), 1, fp);
+ fputc(0x2c, fp);
+ fwrite(&_Sample_Set[i].t, sizeof(double), 1, fp);
+ fputc(0x3b, fp);
+ }
+ // Save _C_x
+ for(int i = 0; i < _Sample_Num - 1; i++) {
+ for(int j = 0; j < 4; j++) {
+ fwrite(&_C_x[j][i], sizeof(double), 1, fp);
+ fputc((j != 3)? 0x2c : 0x3b, fp);
+ }
+ }
+ // Save _C_y
+ for(int i = 0; i < _Sample_Num - 1; i++) {
+ for(int j = 0; j < 4; j++) {
+ fwrite(&_C_y[j][i], sizeof(double), 1, fp);
+ fputc((j != 3)? 0x2c : 0x3b, fp);
+ }
+ }
+
+ fclose(fp);
+
+}
+
+void CubicSpline2d::saveSetting(
+ const char *filename
+)
+{
+ FILE *fp;
+ char *filepath;
+ int fnnum = 0;
+
+ while (filename[fnnum] != 0) fnnum++;
+ filepath = (char *)malloc((fnnum + 8) * sizeof(char)); // "/local/" are 7 char and \0 is 1 char.
+
+ sprintf(filepath, "/local/%s", filename);
+ fp = fopen(filepath, "wb");
+
+ // Save _Sample_Num
+ fwrite(&_Sample_Num, sizeof(unsigned int), 1, fp);
+ fputc(0x3b, fp);
+ // Save _Sample_Set
+ for(int i = 0; i < _Sample_Num; i++) {
+ fwrite(&_Sample_Set[i].x, sizeof(double), 1, fp);
+ fputc(0x2c, fp);
+ fwrite(&_Sample_Set[i].y, sizeof(double), 1, fp);
+ fputc(0x2c, fp);
+ fwrite(&_Sample_Set[i].t, sizeof(double), 1, fp);
+ fputc(0x3b, fp);
+ }
+ // Save _C_x
+ for(int i = 0; i < _Sample_Num - 1; i++) {
+ for(int j = 0; j < 4; j++) {
+ fwrite(&_C_x[j][i], sizeof(double), 1, fp);
+ fputc((j != 3)? 0x2c : 0x3b, fp);
+ }
+ }
+ // Save _C_y
+ for(int i = 0; i < _Sample_Num - 1; i++) {
+ for(int j = 0; j < 4; j++) {
+ fwrite(&_C_y[j][i], sizeof(double), 1, fp);
+ fputc((j != 3)? 0x2c : 0x3b, fp);
+ }
+ }
+
+ fclose(fp);
+ free(filepath);
+}
+
+void CubicSpline2d::loadSetting()
+{
+ FILE *fp;
+ char tmp;
+
+ //sprintf(filepath, "/local/%s", filename);
+ //fp = fopen(filepath, "rb");
+ fp = fopen("/local/savedata.log", "rb");
+
+ // Load _Sample_Num
+ fread(&_Sample_Num, sizeof(unsigned int), 1, fp);
+ fread(&tmp, sizeof(char), 1, fp);
+
+ // Load _Sample_Set
+ for(int i = 0; i < _Sample_Num; i++) {
+ fread(&_Sample_Set[i].x, sizeof(double), 1, fp);
+ fread(&tmp, sizeof(char),1,fp);
+ fread(&_Sample_Set[i].y, sizeof(double), 1, fp);
+ fread(&tmp, sizeof(char),1,fp);
+ fread(&_Sample_Set[i].t, sizeof(double), 1, fp);
+ fread(&tmp, sizeof(char),1,fp);
+ }
+
+ // Load _C_x
+ for(int i = 0; i < _Sample_Num - 1; i++) {
+ for(int j = 0; j < 4; j++) {
+ fread(&_C_x[j][i], sizeof(double), 1, fp);
+ fread(&tmp, sizeof(char),1,fp);
+ }
+ }
+ // Load _C_y
+ for(int i = 0; i < _Sample_Num - 1; i++) {
+ for(int j = 0; j < 4; j++) {
+ fread(&_C_y[j][i], sizeof(double), 1, fp);
+ fread(&tmp, sizeof(char),1,fp);
+ }
+ }
+ fclose(fp);
+}
+
+
+void CubicSpline2d::loadSetting(
+ const char *filename
+)
+{
+ FILE *fp;
+ char *filepath;
+ char tmp;
+ int fnnum = 0;
+
+ while (filename[fnnum] != 0) fnnum++;
+ filepath = (char *)malloc((fnnum + 8) * sizeof(char)); // "/local/" are 7 char and \0 is 1 char.
+
+ sprintf(filepath, "/local/%s", filename);
+ fp = fopen(filepath, "rb");
+
+ // Load _Sample_Num
+ fread(&_Sample_Num, sizeof(unsigned int), 1, fp);
+ fread(&tmp, sizeof(char), 1, fp);
+
+ // Load _Sample_Set
+ for(int i = 0; i < _Sample_Num; i++) {
+ fread(&_Sample_Set[i].x, sizeof(double), 1, fp);
+ fread(&tmp, sizeof(char),1,fp);
+ fread(&_Sample_Set[i].y, sizeof(double), 1, fp);
+ fread(&tmp, sizeof(char),1,fp);
+ fread(&_Sample_Set[i].t, sizeof(double), 1, fp);
+ fread(&tmp, sizeof(char),1,fp);
+ }
+
+ // Load _C_x
+ for(int i = 0; i < _Sample_Num - 1; i++) {
+ for(int j = 0; j < 4; j++) {
+ fread(&_C_x[j][i], sizeof(double), 1, fp);
+ fread(&tmp, sizeof(char),1,fp);
+ }
+ }
+
+ // Load _C_y
+ for(int i = 0; i < _Sample_Num - 1; i++) {
+ for(int j = 0; j < 4; j++) {
+ fread(&_C_y[j][i], sizeof(double), 1, fp);
+ fread(&tmp, sizeof(char),1,fp);
+ }
+ }
+ fclose(fp);
+ free(filepath);
+}
+
+void CubicSpline2d::printOutData()
+{
+ FILE *fp;
+ 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, "\nSample:dst, vol\n");
+ for(int i = 0; i < _Sample_Num; i++) {
+ fprintf(fp, "%f,%f,(%f)\n", _Sample_Set[i].x, _Sample_Set[i].y, _Sample_Set[i].t);
+ }
+ fclose(fp);
+
+}
+void CubicSpline2d::_printOutData(const double *arg, const int num, const char* name)
+{
+ FILE *fp;
+
+ fp = fopen("/local/varlog.txt", "a"); // open file in add mode
+ fprintf(fp, "%10s\n", name);
+ for(int i = 0; i < num; i++) {
+ fprintf(fp, "%.2f, ", arg[i]);
+ }
+ fprintf(fp, "\n");
+ fclose(fp);
+}
+void CubicSpline2d::_printOutDataCouple(const double *arg1, const double *arg2, const int num, const char* name)
+{
+ FILE *fp;
+
+ fp = fopen("/local/varlog.txt", "a"); // open file in add mode
+ fprintf(fp, "%10s\n", name);
+ for(int i = 0; i < num; i++) {
+ fprintf(fp, "(%.2f, %.2f)\n", arg1[i], arg2[i]);
+ }
+ fprintf(fp, "\n");
+ fclose(fp);
+}
+void CubicSpline2d::_printOutData(const Vxyt *arg, int num, const char* name)
+{
+ FILE *fp;
+
+ fp = fopen("/local/varlog.txt", "a"); // open file in add mode
+ fprintf(fp, "%10s\n", name);
+ for(int i = 0; i < num; i++) {
+ fprintf(fp, "%f, ", arg[i].y);
+ }
+ fprintf(fp, "\n");
+ fclose(fp);
+}