TSpline derivative continuity issue

Hello,

I am trying to use the TSpline(3/5) class to parametrize some data points and then use the resulting TSpline object to retrieve the higher oder derivatives.
Reading the (somewhat poor) documentation of the TSpline(3/5) classes, I thought that the spline is build from data at knots requiring that the spline derivatives are continuous up to 2nd order for a TSpline3 and up to 4th order for a TSpline5. I am probably missing something or made a bug in the following simple macro which obviously give non continuous derivatives already at 1st order.
Thanks for any suggestion.

//_______________________________________________________________________________________________________
void GetTaylorTermsAtX(TSpline3 *spl, double x, double *SExpansion){
  // SExpansion array of size 4 (value of function and 3 derivatives at point x
  
  int i = spl->FindX(x);
  double Xi;
  double A0, A1, A2, A3;//  An = d^n spl(x) / dx^n  / n!
  
  spl->GetCoeff(i, Xi, A0, A1, A2, A3);
 
  double p = x - Xi;
  SExpansion[0] = A0 + A1 * p + A2 * p * p + A3 * p * p * p;
  SExpansion[1] = A1 + 2 * A2 * p + 3 * A3 * p * p;
  SExpansion[2] = (2 * A2 + 6 * A3 * p) / 2.;
  SExpansion[3] = A3;

  return;
}// end of function GTaylorTermsAtX
//_______________________________________________________________________________________________________

//_______________________________________________________________________________________________________
void GetTaylorTermsAtX(TSpline5 *spl, double x, double *SExpansion){
  // SExpansion array of size 6 (value of function and 5 derivatives at point x
  
  int i = spl->FindX(x);
  double Xi;
  double A0, A1, A2, A3, A4, A5;//  An = d^n spl(x) / dx^n  / n!
  
  spl->GetCoeff(i, Xi, A0, A1, A2, A3, A4, A5);
  
  double p = x - Xi;
  SExpansion[0] = A0 + A1 * p + A2 * p * p + A3 * p * p * p + A4 * p * p * p * p + A5 * p * p * p * p * p;
  SExpansion[1] = A1 + 2 * A2 * p + 3 * A3 * p * p + 4 * A4 * p * p * p + 5 * A5 * p * p * p * p;
  SExpansion[2] = (2 * A2 + 6 * A3 * p + 12 * A4 * p * p + 20 * A5 * p * p * p ) / 2.;
  SExpansion[3] = (6 * A3 + 24 * A4 * p + 60 * A5 * p * p) / 6.;
  SExpansion[4] = (24 * A4 + 120 * A5 * p) / 24.;
  SExpansion[5] = A5;// (120 / 120)

  return;
}// end of function GTaylorTermsAtX
//_______________________________________________________________________________________________________

//_______________________________________________________________________________________________________
void test_spline(){
  
  const Int_t npts = 100;
  Double_t theta[npts], GE2[npts];
  theta[   0] =  2.44346e+00;  GE2[   0] =  9.05542e-01;
  theta[   1] =  2.44364e+00;  GE2[   1] =  9.05537e-01;
  theta[   2] =  2.44381e+00;  GE2[   2] =  9.05532e-01;
  theta[   3] =  2.44398e+00;  GE2[   3] =  9.05527e-01;
  theta[   4] =  2.44416e+00;  GE2[   4] =  9.05522e-01;
  theta[   5] =  2.44433e+00;  GE2[   5] =  9.05517e-01;
  theta[   6] =  2.44451e+00;  GE2[   6] =  9.05512e-01;
  theta[   7] =  2.44468e+00;  GE2[   7] =  9.05507e-01;
  theta[   8] =  2.44486e+00;  GE2[   8] =  9.05503e-01;
  theta[   9] =  2.44503e+00;  GE2[   9] =  9.05498e-01;
  theta[  10] =  2.44521e+00;  GE2[  10] =  9.05493e-01;
  theta[  11] =  2.44538e+00;  GE2[  11] =  9.05488e-01;
  theta[  12] =  2.44556e+00;  GE2[  12] =  9.05483e-01;
  theta[  13] =  2.44573e+00;  GE2[  13] =  9.05478e-01;
  theta[  14] =  2.44590e+00;  GE2[  14] =  9.05474e-01;
  theta[  15] =  2.44608e+00;  GE2[  15] =  9.05469e-01;
  theta[  16] =  2.44625e+00;  GE2[  16] =  9.05464e-01;
  theta[  17] =  2.44643e+00;  GE2[  17] =  9.05459e-01;
  theta[  18] =  2.44660e+00;  GE2[  18] =  9.05454e-01;
  theta[  19] =  2.44678e+00;  GE2[  19] =  9.05449e-01;
  theta[  20] =  2.44695e+00;  GE2[  20] =  9.05444e-01;
  theta[  21] =  2.44713e+00;  GE2[  21] =  9.05440e-01;
  theta[  22] =  2.44730e+00;  GE2[  22] =  9.05435e-01;
  theta[  23] =  2.44748e+00;  GE2[  23] =  9.05430e-01;
  theta[  24] =  2.44765e+00;  GE2[  24] =  9.05425e-01;
  theta[  25] =  2.44782e+00;  GE2[  25] =  9.05420e-01;
  theta[  26] =  2.44800e+00;  GE2[  26] =  9.05415e-01;
  theta[  27] =  2.44817e+00;  GE2[  27] =  9.05411e-01;
  theta[  28] =  2.44835e+00;  GE2[  28] =  9.05406e-01;
  theta[  29] =  2.44852e+00;  GE2[  29] =  9.05401e-01;
  theta[  30] =  2.44870e+00;  GE2[  30] =  9.05396e-01;
  theta[  31] =  2.44887e+00;  GE2[  31] =  9.05391e-01;
  theta[  32] =  2.44905e+00;  GE2[  32] =  9.05386e-01;
  theta[  33] =  2.44922e+00;  GE2[  33] =  9.05382e-01;
  theta[  34] =  2.44940e+00;  GE2[  34] =  9.05377e-01;
  theta[  35] =  2.44957e+00;  GE2[  35] =  9.05372e-01;
  theta[  36] =  2.44974e+00;  GE2[  36] =  9.05367e-01;
  theta[  37] =  2.44992e+00;  GE2[  37] =  9.05362e-01;
  theta[  38] =  2.45009e+00;  GE2[  38] =  9.05357e-01;
  theta[  39] =  2.45027e+00;  GE2[  39] =  9.05353e-01;
  theta[  40] =  2.45044e+00;  GE2[  40] =  9.05348e-01;
  theta[  41] =  2.45062e+00;  GE2[  41] =  9.05343e-01;
  theta[  42] =  2.45079e+00;  GE2[  42] =  9.05338e-01;
  theta[  43] =  2.45097e+00;  GE2[  43] =  9.05333e-01;
  theta[  44] =  2.45114e+00;  GE2[  44] =  9.05328e-01;
  theta[  45] =  2.45131e+00;  GE2[  45] =  9.05324e-01;
  theta[  46] =  2.45149e+00;  GE2[  46] =  9.05319e-01;
  theta[  47] =  2.45166e+00;  GE2[  47] =  9.05314e-01;
  theta[  48] =  2.45184e+00;  GE2[  48] =  9.05309e-01;
  theta[  49] =  2.45201e+00;  GE2[  49] =  9.05304e-01;
  theta[  50] =  2.45219e+00;  GE2[  50] =  9.05300e-01;
  theta[  51] =  2.45236e+00;  GE2[  51] =  9.05295e-01;
  theta[  52] =  2.45254e+00;  GE2[  52] =  9.05290e-01;
  theta[  53] =  2.45271e+00;  GE2[  53] =  9.05285e-01;
  theta[  54] =  2.45289e+00;  GE2[  54] =  9.05280e-01;
  theta[  55] =  2.45306e+00;  GE2[  55] =  9.05276e-01;
  theta[  56] =  2.45323e+00;  GE2[  56] =  9.05271e-01;
  theta[  57] =  2.45341e+00;  GE2[  57] =  9.05266e-01;
  theta[  58] =  2.45358e+00;  GE2[  58] =  9.05261e-01;
  theta[  59] =  2.45376e+00;  GE2[  59] =  9.05256e-01;
  theta[  60] =  2.45393e+00;  GE2[  60] =  9.05252e-01;
  theta[  61] =  2.45411e+00;  GE2[  61] =  9.05247e-01;
  theta[  62] =  2.45428e+00;  GE2[  62] =  9.05242e-01;
  theta[  63] =  2.45446e+00;  GE2[  63] =  9.05237e-01;
  theta[  64] =  2.45463e+00;  GE2[  64] =  9.05232e-01;
  theta[  65] =  2.45481e+00;  GE2[  65] =  9.05228e-01;
  theta[  66] =  2.45498e+00;  GE2[  66] =  9.05223e-01;
  theta[  67] =  2.45515e+00;  GE2[  67] =  9.05218e-01;
  theta[  68] =  2.45533e+00;  GE2[  68] =  9.05213e-01;
  theta[  69] =  2.45550e+00;  GE2[  69] =  9.05208e-01;
  theta[  70] =  2.45568e+00;  GE2[  70] =  9.05204e-01;
  theta[  71] =  2.45585e+00;  GE2[  71] =  9.05199e-01;
  theta[  72] =  2.45603e+00;  GE2[  72] =  9.05194e-01;
  theta[  73] =  2.45620e+00;  GE2[  73] =  9.05189e-01;
  theta[  74] =  2.45638e+00;  GE2[  74] =  9.05184e-01;
  theta[  75] =  2.45655e+00;  GE2[  75] =  9.05180e-01;
  theta[  76] =  2.45673e+00;  GE2[  76] =  9.05175e-01;
  theta[  77] =  2.45690e+00;  GE2[  77] =  9.05170e-01;
  theta[  78] =  2.45707e+00;  GE2[  78] =  9.05165e-01;
  theta[  79] =  2.45725e+00;  GE2[  79] =  9.05161e-01;
  theta[  80] =  2.45742e+00;  GE2[  80] =  9.05156e-01;
  theta[  81] =  2.45760e+00;  GE2[  81] =  9.05151e-01;
  theta[  82] =  2.45777e+00;  GE2[  82] =  9.05146e-01;
  theta[  83] =  2.45795e+00;  GE2[  83] =  9.05141e-01;
  theta[  84] =  2.45812e+00;  GE2[  84] =  9.05137e-01;
  theta[  85] =  2.45830e+00;  GE2[  85] =  9.05132e-01;
  theta[  86] =  2.45847e+00;  GE2[  86] =  9.05127e-01;
  theta[  87] =  2.45865e+00;  GE2[  87] =  9.05122e-01;
  theta[  88] =  2.45882e+00;  GE2[  88] =  9.05118e-01;
  theta[  89] =  2.45899e+00;  GE2[  89] =  9.05113e-01;
  theta[  90] =  2.45917e+00;  GE2[  90] =  9.05108e-01;
  theta[  91] =  2.45934e+00;  GE2[  91] =  9.05103e-01;
  theta[  92] =  2.45952e+00;  GE2[  92] =  9.05099e-01;
  theta[  93] =  2.45969e+00;  GE2[  93] =  9.05094e-01;
  theta[  94] =  2.45987e+00;  GE2[  94] =  9.05089e-01;
  theta[  95] =  2.46004e+00;  GE2[  95] =  9.05084e-01;
  theta[  96] =  2.46022e+00;  GE2[  96] =  9.05080e-01;
  theta[  97] =  2.46039e+00;  GE2[  97] =  9.05075e-01;
  theta[  98] =  2.46057e+00;  GE2[  98] =  9.05070e-01;
  theta[  99] =  2.46074e+00;  GE2[  99] =  9.05065e-01;

  const Int_t ngr = 4; 
  TSpline3 *GE2_spline = new TSpline3("GE2_spline", theta, GE2, npts);
  GE2_spline->SetLineColor(kBlue);
  GE2_spline->SetLineWidth(2);
  

  // ---------------------------------------------------------------- //
  // -------*-*-*-*-*------ Derivatives of GE2 -------*-*-*-*-*------ //
  // ---------------------------------------------------------------- //
  TCanvas *cnv_GE2 = new TCanvas("cnv_GE2", "cnv_GE2", 2500, 1500);
  cnv_GE2->Divide(ngr / 2, 2);
  Double_t leftmargin = 0.15, rightmargin = 0.15;

  Double_t SExpansion[ngr];
  TString gr_title;
  
  TGraph *gr_GE2_S[ngr];
  for (Int_t igr = 0; igr < ngr; ++igr){
    gr_GE2_S[igr] = new TGraph(npts);
  }
  
  for (Int_t ipt = 0; ipt < npts; ++ipt){
    
    GetTaylorTermsAtX(GE2_spline, theta[ipt], SExpansion);
     
    for (Int_t igr = 0; igr < ngr; ++igr){
      gr_GE2_S[igr]->SetPoint(ipt, theta[ipt], SExpansion[igr]);
    }
  }

  for (Int_t igr = 0; igr < ngr; ++igr){
    
    cnv_GE2->cd(igr+1);
    gPad->SetLeftMargin(leftmargin);
    gPad->SetRightMargin(rightmargin);
    
    gr_GE2_S[igr]->Draw("apl");
    gr_title.Form(";#theta [rd];G_{E}^{2(%d)}", igr);
    gr_GE2_S[igr]->SetTitle(gr_title.Data());
    
    if (igr == 0) GE2_spline->Draw("Lsame");
    
  }
  
}// end of test_spline function
//_______________________________________________________________________________________________________


_ROOT Version: 5.34 to 6.14
Platform: Ubuntu 18.04 LTS
Compiler: gcc 7.3


I think @moneta can most probably give some hints about this

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.