Marjeta Krajnc: Interpolation with spatial rational Pythagorean-hodograph curves of class 4 - notebook with examples

Functions

ClearAll[Norma]
(* Euclidean norm of a vector. *)
Norma[x_]:=(Sqrt[x.x])

ClearAll[bernstein];
(* Bernstein basis polynomials. *)
bernstein[x_,n_,i_]:=Module[{},
If[x==0 && i==0,Return[1]];
If[x==1 && i==n,Return[1]];
Return[Binomial[n,i]*x^(i)*(1-x)^(n-i)]
]

ClearAll[Rotacija]
(* Rotation matrix from a quaternion Q. *)
Rotacija[Q_]:=Module[{U,u0},
U={{Q[[1]]^2+Q[[2]]^2-Q[[3]]^2-Q[[4]]^2,2(Q[[2]]Q[[3]]-Q[[1]]Q[[4]]),2(Q[[2]]Q[[4]]+Q[[1]]Q[[3]])},{2(Q[[2]]Q[[3]]+Q[[1]]Q[[4]]),Q[[1]]^2-Q[[2]]^2+Q[[3]]^2-Q[[4]]^2,2(Q[[3]]Q[[4]]-Q[[1]]Q[[2]])},{2(Q[[2]]Q[[4]]-Q[[1]]Q[[3]]),2(Q[[3]]Q[[4]]+Q[[1]]Q[[2]]),Q[[1]]^2-Q[[2]]^2-Q[[3]]^2+Q[[4]]^2}};
u0=Q[[1]]^2+Q[[2]]^2+Q[[3]]^2+Q[[4]]^2;
Return[{U,u0,1/u0 U}]
]

ClearAll[DolociKvat]
(* quaternion Q from a rotation matrix. *)
DolociKvat[M_]:=Module[{QQ1,QQ2,QQ3,QQ4,QQ,q},
QQ1={1+M[[1,1]]+M[[2,2]]+M[[3,3]],M[[3,2]]-M[[2,3]],M[[1,3]]-M[[3,1]],M[[2,1]]-M[[1,2]]};
QQ2={M[[3,2]]-M[[2,3]],1+M[[1,1]]-M[[2,2]]-M[[3,3]],M[[2,1]]+M[[1,2]],M[[1,3]]+M[[3,1]]};
QQ3={M[[1,3]]-M[[3,1]],M[[2,1]]-M[[1,2]],1-M[[1,1]]+M[[2,2]]-M[[3,3]],M[[3,2]]+M[[2,3]]};
QQ4={M[[2,1]]-M[[1,2]],M[[1,3]]+M[[3,1]],M[[3,2]]+M[[2,3]],1-M[[1,1]]-M[[2,2]]+M[[3,3]]};
QQ=Select[{QQ1,QQ2,QQ3,QQ4},#≠{0,0,0,0}&];
q=QQ[[1]];
Return[q];
]

ClearAll[KvaternionProd,KvaternionConj]
(* Quaternion product. *)
KvaternionProd[X_,Y_]:=Module[{x0,x,y0,y},
x0=X[[1]];
x=Drop[X,1];
y0=Y[[1]];
y=Drop[Y,1];
Return[Flatten[{x0 y0-x.y,x0 y+y0 x + Cross[x,y]}]];
]
(* Quaternion conjugate. *)
KvaternionConj[X_]:=Flatten[{X[[1]],-Drop[X,1]}]

ClearAll[KvaternionNorm]
(* Quaternion norm. *)
KvaternionNorm[X_]:=((Sqrt[KvaternionProd[X,KvaternionConj[X]]]))[[1]]

ClearAll[OPStar]
(* Operation star between two quaternions. *)
OPStar[A_,B_]:=1/2(KvaternionProd[A,KvaternionProd[{0,1,0,0},KvaternionConj[B]]]+KvaternionProd[B,KvaternionProd[{0,1,0,0},KvaternionConj[A]]])
ClearAll[OPStar2]
OPStar2[A_,B_]:=1/2(KvaternionProd[A,KvaternionProd[{0,0,1,0},KvaternionConj[B]]]+KvaternionProd[B,KvaternionProd[{0,0,1,0},KvaternionConj[A]]])
ClearAll[OPStar3]
OPStar3[A_,B_]:=1/2(KvaternionProd[A,KvaternionProd[{0,0,0,1},KvaternionConj[B]]]+KvaternionProd[B,KvaternionProd[{0,0,0,1},KvaternionConj[A]]])
ClearAll[QToR]
QToR[Q_,p_]:=Simplify[KvaternionProd[Q,KvaternionProd[p,KvaternionConj[Q]]]];


ClearAll[SolveEQStar]
(* Solution of a star equation. *)
SolveEQStar[a_,phi_]:=Module[{A,nA,X},
A=Prepend[a,0];
nA=KvaternionNorm[A];
X=(Sqrt[nA])KvaternionProd[((A/nA+{0,1,0,0})/(KvaternionNorm[A/nA+{0,1,0,0}])),{Cos[phi],Sin[phi],0,0}];
Return[X];
]

ClearAll[WedgeProd]
(* Wedge product between three vectors in R^4. *)
WedgeProd[v1_,v2_,v3_]:=Module[{A},
A=Transpose[{v1,v2,v3}];
Return[{-Det[{A[[2]],A[[3]],A[[4]]}],Det[{A[[1]],A[[3]],A[[4]]}],Det[{A[[1]],A[[4]],A[[2]]}],Det[{A[[1]],A[[2]],A[[3]]}]}]
]


ClearAll[IzracunERF]
(* Computation of ER frame from the quaternion polynomial. *)
IzracunERF[AA_,t_]:=Module[{uu,vv,pp,qq,du,dv,dp,dq,WW1,WW2,WW3,WI,W1,W2,W3,EE1,EE2,EE3,UN},
uu=AA[[1]];
vv=AA[[2]];
pp=AA[[3]];
qq=AA[[4]];
du=D[uu,t];
dv=D[vv,t];
dp=D[pp,t];
dq=D[qq,t];
(* Compute Wi *)
WW1=2(uu dv -du vv -pp dq + dp qq);
WW2=2(uu dp - du pp + vv dq-dv qq);
WW3=2(uu dq - du qq-vv dp+dv pp);
WI=uu^2+vv^2+pp^2+qq^2;
W1=WW1/WI;
W2=WW2/WI;
W3=WW3/WI;
(* compute ERF frame *)
EE1={uu^2+vv^2-pp^2-qq^2,2(vv pp + uu qq),2(vv qq -uu pp)};
EE2={2(vv pp - uu qq),uu^2-vv^2+pp^2-qq^2,2(pp qq + uu vv)};
EE3={2(vv qq + uu pp),2(pp qq- uu vv), uu^2-vv^2-pp^2+qq^2};
UN=WW2 EE2 + WW3 EE3;
Return[{WW1,WW2,WW3,WI,EE1,EE2,EE3,UN}]
]


ClearAll[DeCasteljau]
(* De Casteljau algorithm. *)
DeCasteljau[KT_,t_]:=Module[{M,n,r,i,j},
n=Length[KT]-1;
M=Table[Table[0,{i,0,n}],{j,0,n}];
For[i=0,i≤n,i++,
M[[i+1,1]]=KT[[i+1]];
];
For[r=1,r≤n,r++,
For[i=0,i≤n-r,i++,
M[[i+1,r+1]]=(1-t)M[[i+1,r]]+t M[[i+2,r]];
];
];
Return[M]
]

ClearAll[PrehodnaM]
(* Connection matrix between a Bernstein and standard polynomial basis. *)
PrehodnaM[n_]:=Module[{p,t,i},
p=Table[CoefficientList[bernstein[t,n,i],t],{i,0,n}];
Return[p]
]

ClearAll[BezKoef]
(* Coefficients in a Bernstein basis from the coefficients in a standard basis. *)
BezKoef[n_,koef_]:=Module[{IA,i,koef1},
koef1=koef;
For[i=Length[koef],i≤n,i++,
koef1=Append[koef1,0]
];
IA=Inverse[PrehodnaM[n]];
Return[koef1.IA]
]



ClearAll[FactorPSI]
(* Factor Psi from a dual form. *)
FactorPSI[f_,g_,h_,t_]:=D[f/g Det[{h,D[h,{t,2}],D[h,{t,3}]}]/Det[{h,D[h,t],D[h,{t,2}]}]^2+D[(D[f/g,t]1/Det[{h,D[h,t],D[h,{t,2}]}]),t],t]+f/g Det[{D[h,t],D[h,{t,2}],D[h,{t,3}]}]/Det[{h,D[h,t],D[h,{t,2}]}]^2


ClearAll[Frene]
(* Frene frame. *)
Frene[r_,t_]:=Module[{tt,bb,nn},
tt=D[r,t]/Norma[D[r,t]];
bb=Cross[D[r,t],D[r,{t,2}]]/Norma[Cross[D[r,t],D[r,{t,2}]]];
nn=Cross[bb,tt];
Return[{tt,nn,bb}]
]

ClearAll[Simpson]
(* Simpson integration rule. *)
Simpson[f_,t_,a_,b_,m_]:=Module[{x,y1,y2,y3,i},
x=Table[a+i*(b-a)/(2m),{i,0,2m}];
y1=4*Plus@@(Table[f/.{t->x[[2i]]},{i,1,m}]);
y2=2*Plus@@(Table[f/.{t->x[[2i+1]]},{i,1,m-1}]);
y3=(f/.{t->a})+(f/.{t->b});
Return[(b-a)/(6m)*(y1+y2+y3)]
]

ClearAll[Ukriv]
(* Curvature. *)
Ukriv[R_,t_]:=Norma[Cross[D[R,t],D[R,{t,2}]]]/Norma[D[R,t]]^3;

ClearAll[Red1];
(* Function computes the approximation order of two consecutive measurements. *)
Red1[napaka1_,h_,napaka2_,k_]:=(Log[napaka1]-Log[napaka2])/(Log[h]-Log[k])

ClearAll[RedAproksimacije];
(* Function computes an approximation order. *)
RedAproksimacije[TabH_,TabN_]:=Module[{tabela,i},
If[Length[TabH]≠Length[TabN],Return[{}]];
tabela={{"h","napaka","eksponent padanja"}};
tabela=Append[tabela,{TabH[[1]],TabN[[1]],"/"}];
For[i=2,i≤Length[TabH],i++,
tabela=Append[tabela,{TabH[[i]],TabN[[i]],Red1[TabN[[i-1]],TabH[[i-1]],TabN[[i]],TabH[[i]]]}];
];
Return[tabela];
]


ClearAll[Red1];
(* Function computes the approximation order of two consecutive measurements. *)
Red1[napaka1_,h_,napaka2_,k_]:=(Log[napaka1]-Log[napaka2])/(Log[h]-Log[k])

ClearAll[RedAproksimacije];
(* Function computes an approximation order. *)
RedAproksimacije[TabH_,TabN_]:=Module[{tabela,i},
If[Length[TabH]≠Length[TabN],Return[{}]];
tabela={{"h","error","decay exponent"}};
tabela=Append[tabela,{TabH[[1]],TabN[[1]],"/"}];
For[i=2,i≤Length[TabH],i++,
tabela=Append[tabela,{TabH[[i]],TabN[[i]],Red1[TabN[[i-1]],TabH[[i-1]],TabN[[i]],TabH[[i]]]}];
];
Return[tabela];
]

ClearAll[RedAproksimacijeZaZlepek];
(* Function computes an approximation order. *)
RedAproksimacijeZaZlepek[TabM_,TabN_]:=Module[{tabela,i},
If[Length[TabM]≠Length[TabN],Return[{}]];
tabela={{"m","error","decay exponent"}};
tabela=Append[tabela,{TabM[[1]],TabN[[1]],"/"}];
For[i=2,i≤Length[TabM],i++,
tabela=Append[tabela,{TabM[[i]],TabN[[i]],Red1[TabN[[i-1]],1/TabM[[i-1]],TabN[[i]],1/TabM[[i]]]}];
];
Return[tabela];
]

(* ------ functions for rational G1 interpolant of class 4 ---------------- *)

ClearAll[IntervalForPhi]
(* Function returns the admissible interval for phi based on the sign of ksi - signKsi ∈ {-1,1} *)
IntervalForPhi[P1_,d0_,d1_,signKsi_]:=Module[{CC,d,P,N1,N2,inter,N3,N3a,N3b},
{(Subscript[d, 0,1]),(Subscript[d, 0,2]),(Subscript[d, 0,3])}=d0;
{(Subscript[d, 1,1]),(Subscript[d, 1,2]),(Subscript[d, 1,3])}=d1;
{(Subscript[P, 1]),(Subscript[P, 2]),(Subscript[P, 3])}=P1;
{(Subscript[CC, 1]),(Subscript[CC, 2]),(Subscript[CC, 3])}={(Subscript[P, 3]) (Subscript[d, 1,2])-(Subscript[P, 2]) (Subscript[d, 1,3]),(Subscript[P, 2]) (Subscript[d, 1,2])+(Subscript[P, 3]) (Subscript[d, 1,3]),(Subscript[P, 1]) (-1+(Subscript[d, 1,1]^2))+(Subscript[d, 1,1]) ((Subscript[P, 2]) (Subscript[d, 1,2])+(Subscript[P, 3]) (Subscript[d, 1,3]))};
N1=ArcCot[(Subscript[CC, 2])/(Subscript[CC, 1])];
If[N1<0,N1=N1+Pi];
N2=ArcCot[ -(Subscript[CC, 3]) /(Subscript[CC, 1])];
If[N2<0,N2=N2+Pi];
N3a=ArcCot[(((Subscript[CC, 2])-(Subscript[CC, 3])+(Sqrt[(Subscript[CC, 2]^2)+(Subscript[CC, 3]^2)-2 (Subscript[CC, 2]) (Subscript[CC, 3]) (Subscript[d, 1,1])+(Subscript[CC, 1]^2) (1-(Subscript[d, 1,1]^2))]))/((Subscript[CC, 1]) (1-(Subscript[d, 1,1]))))];
N3b=ArcCot[(((Subscript[CC, 2])-(Subscript[CC, 3])-(Sqrt[(Subscript[CC, 2]^2)+(Subscript[CC, 3]^2)-2 (Subscript[CC, 2]) (Subscript[CC, 3]) (Subscript[d, 1,1])+(Subscript[CC, 1]^2) (1-(Subscript[d, 1,1]^2))]))/((Subscript[CC, 1]) (1-(Subscript[d, 1,1]))))];
If[N3a<0,N3a=N3a+Pi];
If[N3b<0,N3b=N3b+Pi];
N3=Sort[{N3a,N3b},#1<#2&];
(* Final intervals where we can choose φ from: *)
If[(signKsi (Subscript[CC, 1]))>0,inter={Max[N3],Pi},inter={0,Min[N3]}];
Return[{N1,N2,N3,inter}]
]

ClearAll[ConstC1C2C3]
(* Constants C1, C2, C3. *)
ConstC1C2C3[P1_,d1_]:=Module[{d,P},
{(Subscript[d, 1,1]),(Subscript[d, 1,2]),(Subscript[d, 1,3])}=d1;
{(Subscript[P, 1]),(Subscript[P, 2]),(Subscript[P, 3])}=P1;
Return[{(Subscript[P, 3]) (Subscript[d, 1,2])-(Subscript[P, 2]) (Subscript[d, 1,3]),(Subscript[P, 2]) (Subscript[d, 1,2])+(Subscript[P, 3]) (Subscript[d, 1,3]),(Subscript[P, 1]) (-1+(Subscript[d, 1,1]^2))+(Subscript[d, 1,1]) ((Subscript[P, 2]) (Subscript[d, 1,2])+(Subscript[P, 3]) (Subscript[d, 1,3]))}]
]

ClearAll[FunDelta]
(* Functions (Subscript[δ, 1]), (Subscript[δ, 2]), (Subscript[δ, 3]). *)
FunDelta[φ_,P1_,d0_,d1_]:=Module[{d,P,CC,del1,del2,del3},
{(Subscript[d, 0,1]),(Subscript[d, 0,2]),(Subscript[d, 0,3])}=d0;
{(Subscript[d, 1,1]),(Subscript[d, 1,2]),(Subscript[d, 1,3])}=d1;
{(Subscript[P, 1]),(Subscript[P, 2]),(Subscript[P, 3])}=P1;
{(Subscript[CC, 1]),(Subscript[CC, 2]),(Subscript[CC, 3])}={(Subscript[P, 3]) (Subscript[d, 1,2])-(Subscript[P, 2]) (Subscript[d, 1,3]),(Subscript[P, 2]) (Subscript[d, 1,2])+(Subscript[P, 3]) (Subscript[d, 1,3]),(Subscript[P, 1]) (-1+(Subscript[d, 1,1]^2))+(Subscript[d, 1,1]) ((Subscript[P, 2]) (Subscript[d, 1,2])+(Subscript[P, 3]) (Subscript[d, 1,3]))};
del1=-Cos[φ] (Subscript[CC, 1])+Sin[φ] (Subscript[CC, 2]);
del2=Cos[φ] (Subscript[CC, 1])+Sin[φ] (Subscript[CC, 3]);
del3=(Subscript[CC, 2]) (Subscript[CC, 3])+Sin[2 φ] ((Subscript[CC, 1]) (Subscript[CC, 2])-(Subscript[CC, 1]) (Subscript[CC, 3]))+Cos[2 φ] (-(Subscript[CC, 1]^2)-(Subscript[CC, 2]) (Subscript[CC, 3]))+(Subscript[CC, 1]^2) (Subscript[d, 1,1]);
Return[{del1,del2,del3}]
]

ClearAll[IntervalForKsi]
(* Function returns the admissible interval for ξ>0 based od the chosen φ. *)
IntervalForKsi[P1_,d0_,d1_,phi_]:=Module[{CC,d,P,del1,del2,del4,int},
{(Subscript[d, 0,1]),(Subscript[d, 0,2]),(Subscript[d, 0,3])}=d0;
{(Subscript[d, 1,1]),(Subscript[d, 1,2]),(Subscript[d, 1,3])}=d1;
{(Subscript[P, 1]),(Subscript[P, 2]),(Subscript[P, 3])}=P1;
{(Subscript[CC, 1]),(Subscript[CC, 2]),(Subscript[CC, 3])}={(Subscript[P, 3]) (Subscript[d, 1,2])-(Subscript[P, 2]) (Subscript[d, 1,3]),(Subscript[P, 2]) (Subscript[d, 1,2])+(Subscript[P, 3]) (Subscript[d, 1,3]),(Subscript[P, 1]) (-1+(Subscript[d, 1,1]^2))+(Subscript[d, 1,1]) ((Subscript[P, 2]) (Subscript[d, 1,2])+(Subscript[P, 3]) (Subscript[d, 1,3]))};
del1=Sin[phi] (Subscript[CC, 2])-Cos[phi] (Subscript[CC, 1]);
del2=(Cos[phi] (Subscript[CC, 1])+Sin[phi] ((Subscript[CC, 2]) (Subscript[d, 1,1])+(Subscript[P, 1]) (-1+(Subscript[d, 1,1]^2))));
del4=(-2del1 del2)/((Subscript[CC, 1])^2(1+(Subscript[d, 1,1])));
int={1/(Sqrt[2 del4-2]),(Sqrt[((1+(Sqrt[del4]))/(4 del4-4))])};
Return[int];

]

ClearAll[G1InterpolantSol]
(* Function returns two solutions for A0,A1,f3,f4,lambda. *)
G1InterpolantSol[P1_,d0_,d1_,phi_,xi_]:=Module[{A0,A1,lam,φ,ξ,d,λ,f,P,f3,f4,resLam},
ξ=xi;
φ=phi;
{(Subscript[d, 0,1]),(Subscript[d, 0,2]),(Subscript[d, 0,3])}=d0;
{(Subscript[d, 1,1]),(Subscript[d, 1,2]),(Subscript[d, 1,3])}=d1;
{(Subscript[P, 1]),(Subscript[P, 2]),(Subscript[P, 3])}=P1;
A0={1,0,0,0};
A1={((λ Cos[φ] (Sqrt[1+(Subscript[d, 1,1])]))/(Sqrt[2])),-((λ Sin[φ] (Sqrt[1+(Subscript[d, 1,1])]))/(Sqrt[2])),((λ (-Sin[φ] (Subscript[d, 1,2])-Cos[φ] (Subscript[d, 1,3])))/(Sqrt[2+2 (Subscript[d, 1,1])])),((λ (Cos[φ] (Subscript[d, 1,2])-Sin[φ] (Subscript[d, 1,3])))/(Sqrt[2+2 (Subscript[d, 1,1])]))};
{f3,f4}={-(1/ξ)4 λ (λ ((Subscript[P, 3]) (Subscript[d, 1,2])-(Subscript[P, 2]) (Subscript[d, 1,3]))+2 λ (ξ^2) ((Subscript[P, 3]) (Subscript[d, 1,2])-(Subscript[P, 2]) (Subscript[d, 1,3]))+((2 (Sqrt[2]) (ξ^2) (Cos[φ] ((Subscript[P, 3]) (Subscript[d, 1,2])-(Subscript[P, 2]) (Subscript[d, 1,3]))-Sin[φ] ((Subscript[P, 2]) (Subscript[d, 1,2])+(Subscript[P, 3]) (Subscript[d, 1,3]))))/(Sqrt[1+(Subscript[d, 1,1])]))),(1/((ξ^2) (Sqrt[1+(Subscript[d, 1,1])])))λ ((2 (Sqrt[2]) (ξ^2) (-1+4 ξ) Cos[φ]+λ (-1+4 ξ (1+ξ (-1+2 ξ))) (Sqrt[1+(Subscript[d, 1,1])])) ((Subscript[P, 3]) (Subscript[d, 1,2])-(Subscript[P, 2]) (Subscript[d, 1,3]))+2 (Sqrt[2]) (1-4 ξ) (ξ^2) Sin[φ] ((Subscript[P, 2]) (Subscript[d, 1,2])+(Subscript[P, 3]) (Subscript[d, 1,3])))};
resLam={{λ->(-(Subscript[P, 3]) (Sqrt[1+(Subscript[d, 1,1])]) (Subscript[d, 1,2])-4 (ξ^2) (Subscript[P, 3]) (Sqrt[1+(Subscript[d, 1,1])]) (Subscript[d, 1,2])+(Subscript[P, 2]) (Sqrt[1+(Subscript[d, 1,1])]) (Subscript[d, 1,3])+4 (ξ^2) (Subscript[P, 2]) (Sqrt[1+(Subscript[d, 1,1])]) (Subscript[d, 1,3])+(1/(Sqrt[2]))(√(2 (((1+4 (ξ^2)))^2) (1+(Subscript[d, 1,1])) ((((Subscript[P, 3]) (Subscript[d, 1,2])-(Subscript[P, 2]) (Subscript[d, 1,3])))^2)+64 (ξ^4) ((Subscript[P, 2]) (Sin[φ] (Subscript[d, 1,2])+Cos[φ] (Subscript[d, 1,3]))+(Subscript[P, 3]) (-Cos[φ] (Subscript[d, 1,2])+Sin[φ] (Subscript[d, 1,3]))) (Sin[φ] (Subscript[P, 1]) (-1+(Subscript[d, 1,1]^2))+(Subscript[P, 2]) (Sin[φ] (Subscript[d, 1,1]) (Subscript[d, 1,2])-Cos[φ] (Subscript[d, 1,3]))+(Subscript[P, 3]) (Cos[φ] (Subscript[d, 1,2])+Sin[φ] (Subscript[d, 1,1]) (Subscript[d, 1,3]))))))/(4 (Sqrt[2]) (ξ^2) (Sin[φ] (Subscript[P, 1]) (-1+(Subscript[d, 1,1]^2))+(Subscript[P, 2]) (Sin[φ] (Subscript[d, 1,1]) (Subscript[d, 1,2])-Cos[φ] (Subscript[d, 1,3]))+(Subscript[P, 3]) (Cos[φ] (Subscript[d, 1,2])+Sin[φ] (Subscript[d, 1,1]) (Subscript[d, 1,3]))))},{λ->-(((Subscript[P, 3]) (Sqrt[1+(Subscript[d, 1,1])]) (Subscript[d, 1,2])+4 (ξ^2) (Subscript[P, 3]) (Sqrt[1+(Subscript[d, 1,1])]) (Subscript[d, 1,2])-(Subscript[P, 2]) (Sqrt[1+(Subscript[d, 1,1])]) (Subscript[d, 1,3])-4 (ξ^2) (Subscript[P, 2]) (Sqrt[1+(Subscript[d, 1,1])]) (Subscript[d, 1,3])+(1/(Sqrt[2]))(√(2 (((1+4 (ξ^2)))^2) (1+(Subscript[d, 1,1])) ((((Subscript[P, 3]) (Subscript[d, 1,2])-(Subscript[P, 2]) (Subscript[d, 1,3])))^2)+64 (ξ^4) ((Subscript[P, 2]) (Sin[φ] (Subscript[d, 1,2])+Cos[φ] (Subscript[d, 1,3]))+(Subscript[P, 3]) (-Cos[φ] (Subscript[d, 1,2])+Sin[φ] (Subscript[d, 1,3]))) (Sin[φ] (Subscript[P, 1]) (-1+(Subscript[d, 1,1]^2))+(Subscript[P, 2]) (Sin[φ] (Subscript[d, 1,1]) (Subscript[d, 1,2])-Cos[φ] (Subscript[d, 1,3]))+(Subscript[P, 3]) (Cos[φ] (Subscript[d, 1,2])+Sin[φ] (Subscript[d, 1,1]) (Subscript[d, 1,3]))))))/(4 (Sqrt[2]) (ξ^2) (Sin[φ] (Subscript[P, 1]) (-1+(Subscript[d, 1,1]^2))+(Subscript[P, 2]) (Sin[φ] (Subscript[d, 1,1]) (Subscript[d, 1,2])-Cos[φ] (Subscript[d, 1,3]))+(Subscript[P, 3]) (Cos[φ] (Subscript[d, 1,2])+Sin[φ] (Subscript[d, 1,1]) (Subscript[d, 1,3])))))}};
Return[{{A0,A1,f3,f4,λ/.resLam[[1]]}/.resLam[[1]],{A0,A1,f3,f4,λ/.resLam[[2]]}/.resLam[[2]]}]

]

ClearAll[G1Interpolant]
(* Function returns the chosen G1 interpolant. For num==1 we obtain (SubPlus[r]) and for num=2 we obtain (SubMinus[r]). *)
G1Interpolant[t_,P1_,d0_,d1_,phi_,xi_,num_]:=Module[{a,A1,f3,f4,resLam,ξ,φ,UR,L,DL,D2L,PP,RR,i},
ξ=xi;
φ=phi;
{{(Subscript[a, 0,0]),(Subscript[a, 0,1]),(Subscript[a, 0,2]),(Subscript[a, 0,3])},{(Subscript[a, 1,0]),(Subscript[a, 1,1]),(Subscript[a, 1,2]),(Subscript[a, 1,3])},f3,f4,resLam}=G1InterpolantSol[P1,d0,d1,phi,xi][[num]];
UR={8 (-(Subscript[a, 0,2]) (Subscript[a, 1,0])-(Subscript[a, 0,3]) (Subscript[a, 1,1])+(Subscript[a, 0,0]) (Subscript[a, 1,2])+(Subscript[a, 0,1]) (Subscript[a, 1,3])) ((-(-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,1])+t (t+2 ξ-2 t ξ) (Subscript[a, 1,1])) (-(-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,2])+t (t+2 ξ-2 t ξ) (Subscript[a, 1,2]))-(-(-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,0])+t (t+2 ξ-2 t ξ) (Subscript[a, 1,0])) (-(-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,3])+t (t+2 ξ-2 t ξ) (Subscript[a, 1,3])))+8 (-(Subscript[a, 0,3]) (Subscript[a, 1,0])+(Subscript[a, 0,2]) (Subscript[a, 1,1])-(Subscript[a, 0,1]) (Subscript[a, 1,2])+(Subscript[a, 0,0]) (Subscript[a, 1,3])) ((-(-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,0])+t (t+2 ξ-2 t ξ) (Subscript[a, 1,0])) (-(-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,2])+t (t+2 ξ-2 t ξ) (Subscript[a, 1,2]))+(-(-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,1])+t (t+2 ξ-2 t ξ) (Subscript[a, 1,1])) (-(-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,3])+t (t+2 ξ-2 t ξ) (Subscript[a, 1,3]))),4 (-(Subscript[a, 0,2]) (Subscript[a, 1,0])-(Subscript[a, 0,3]) (Subscript[a, 1,1])+(Subscript[a, 0,0]) (Subscript[a, 1,2])+(Subscript[a, 0,1]) (Subscript[a, 1,3])) (((((-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,0])-t (t+2 ξ-2 t ξ) (Subscript[a, 1,0])))^2)-((((-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,1])-t (t+2 ξ-2 t ξ) (Subscript[a, 1,1])))^2)+((((-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,2])-t (t+2 ξ-2 t ξ) (Subscript[a, 1,2])))^2)-((((-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,3])-t (t+2 ξ-2 t ξ) (Subscript[a, 1,3])))^2))+8 (-(Subscript[a, 0,3]) (Subscript[a, 1,0])+(Subscript[a, 0,2]) (Subscript[a, 1,1])-(Subscript[a, 0,1]) (Subscript[a, 1,2])+(Subscript[a, 0,0]) (Subscript[a, 1,3])) (-(-(-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,0])+t (t+2 ξ-2 t ξ) (Subscript[a, 1,0])) (-(-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,1])+t (t+2 ξ-2 t ξ) (Subscript[a, 1,1]))+(-(-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,2])+t (t+2 ξ-2 t ξ) (Subscript[a, 1,2])) (-(-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,3])+t (t+2 ξ-2 t ξ) (Subscript[a, 1,3]))),4 (-(Subscript[a, 0,3]) (Subscript[a, 1,0])+(Subscript[a, 0,2]) (Subscript[a, 1,1])-(Subscript[a, 0,1]) (Subscript[a, 1,2])+(Subscript[a, 0,0]) (Subscript[a, 1,3])) (((((-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,0])-t (t+2 ξ-2 t ξ) (Subscript[a, 1,0])))^2)-((((-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,1])-t (t+2 ξ-2 t ξ) (Subscript[a, 1,1])))^2)-((((-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,2])-t (t+2 ξ-2 t ξ) (Subscript[a, 1,2])))^2)+((((-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,3])-t (t+2 ξ-2 t ξ) (Subscript[a, 1,3])))^2))+8 (-(Subscript[a, 0,2]) (Subscript[a, 1,0])-(Subscript[a, 0,3]) (Subscript[a, 1,1])+(Subscript[a, 0,0]) (Subscript[a, 1,2])+(Subscript[a, 0,1]) (Subscript[a, 1,3])) ((-(-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,0])+t (t+2 ξ-2 t ξ) (Subscript[a, 1,0])) (-(-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,1])+t (t+2 ξ-2 t ξ) (Subscript[a, 1,1]))+(-(-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,2])+t (t+2 ξ-2 t ξ) (Subscript[a, 1,2])) (-(-1+t) (1+t (-1+2 ξ)) (Subscript[a, 0,3])+t (t+2 ξ-2 t ξ) (Subscript[a, 1,3])))};
L={-(f3 t^3+f4 t^4) , UR[[1]], UR[[2]], UR[[3]]};
DL=D[L,t];
D2L=D[L,{t,2}];
PP=WedgeProd[L,DL,D2L];
(*RR=Table[PP[[i]]/PP[[1]],{i,2,4}];*)
Return[PP]
]

ClearAll[G1InterpolantFactor]
(* Function returns (Subscript[ψ, 0(Subscript[,,  ])])(Subscript[ψ, 1, ])numerator(ψ), theta. *)
G1InterpolantFactor[t_,P1_,d0_,d1_,phi_,xi_,num_]:=Module[{d,P,theta,ξ,φ,a,f3,f4,resLam,bezK,facN,fac,fac0,fac1,i},
ξ=xi;
φ=phi;
{(Subscript[d, 0,1]),(Subscript[d, 0,2]),(Subscript[d, 0,3])}=d0;
{(Subscript[d, 1,1]),(Subscript[d, 1,2]),(Subscript[d, 1,3])}=d1;
{(Subscript[P, 1]),(Subscript[P, 2]),(Subscript[P, 3])}=P1;
{{(Subscript[a, 0,0]),(Subscript[a, 0,1]),(Subscript[a, 0,2]),(Subscript[a, 0,3])},{(Subscript[a, 1,0]),(Subscript[a, 1,1]),(Subscript[a, 1,2]),(Subscript[a, 1,3])},f3,f4,resLam}=G1InterpolantSol[P1,d0,d1,phi,xi][[num]];
theta=t (1-2 ξ)+ξ+(t^2) (-1+2 ξ);
bezK={6 (ξ^2) f3,3 ξ ((-1+4 ξ) f3+4 ξ f4),-3 ((1-4 ξ+2 (ξ^2)) f3-4 ξ f4)}/(-64 ((Subscript[a, 0,1]^3) (Subscript[a, 1,0]) ((Subscript[a, 1,2]^2)+(Subscript[a, 1,3]^2))-((Subscript[a, 0,0]) (Subscript[a, 1,1])+(Subscript[a, 0,3]) (Subscript[a, 1,2])-(Subscript[a, 0,2]) (Subscript[a, 1,3])) ((Subscript[a, 0,2]^2) ((Subscript[a, 1,0]^2)+(Subscript[a, 1,1]^2))+(Subscript[a, 0,3]^2) ((Subscript[a, 1,0]^2)+(Subscript[a, 1,1]^2))-2 (Subscript[a, 0,0]) (Subscript[a, 0,3]) ((Subscript[a, 1,1]) (Subscript[a, 1,2])+(Subscript[a, 1,0]) (Subscript[a, 1,3]))+2 (Subscript[a, 0,0]) (Subscript[a, 0,2]) (-(Subscript[a, 1,0]) (Subscript[a, 1,2])+(Subscript[a, 1,1]) (Subscript[a, 1,3]))+(Subscript[a, 0,0]^2) ((Subscript[a, 1,2]^2)+(Subscript[a, 1,3]^2)))+(Subscript[a, 0,1]^2) (-(Subscript[a, 0,0]) (Subscript[a, 1,1]) ((Subscript[a, 1,2]^2)+(Subscript[a, 1,3]^2))+(Subscript[a, 0,2]) (-2 (Subscript[a, 1,0]) (Subscript[a, 1,1]) (Subscript[a, 1,2])-2 (Subscript[a, 1,0]^2) (Subscript[a, 1,3])+(Subscript[a, 1,2]^2) (Subscript[a, 1,3])+(Subscript[a, 1,3]^3))+(Subscript[a, 0,3]) (2 (Subscript[a, 1,0]^2) (Subscript[a, 1,2])-2 (Subscript[a, 1,0]) (Subscript[a, 1,1]) (Subscript[a, 1,3])-(Subscript[a, 1,2]) ((Subscript[a, 1,2]^2)+(Subscript[a, 1,3]^2))))+(Subscript[a, 0,1]) (-2 (Subscript[a, 0,0]) (Subscript[a, 0,3]) (2 (Subscript[a, 1,0]) (Subscript[a, 1,1]) (Subscript[a, 1,2])+(Subscript[a, 1,0]^2) (Subscript[a, 1,3])-(Subscript[a, 1,1]^2) (Subscript[a, 1,3]))+(Subscript[a, 0,3]^2) ((Subscript[a, 1,0]^3)+(Subscript[a, 1,0]) ((Subscript[a, 1,1]^2)-2 (Subscript[a, 1,2]^2))+2 (Subscript[a, 1,1]) (Subscript[a, 1,2]) (Subscript[a, 1,3]))+(Subscript[a, 0,0]^2) (Subscript[a, 1,0]) ((Subscript[a, 1,2]^2)+(Subscript[a, 1,3]^2))+(Subscript[a, 0,2]^2) ((Subscript[a, 1,0]^3)-2 (Subscript[a, 1,1]) (Subscript[a, 1,2]) (Subscript[a, 1,3])+(Subscript[a, 1,0]) ((Subscript[a, 1,1]^2)-2 (Subscript[a, 1,3]^2)))+(Subscript[a, 0,2]) ((Subscript[a, 0,0]) (-2 (Subscript[a, 1,0]^2) (Subscript[a, 1,2])+2 (Subscript[a, 1,1]^2) (Subscript[a, 1,2])+4 (Subscript[a, 1,0]) (Subscript[a, 1,1]) (Subscript[a, 1,3]))+2 (Subscript[a, 0,3]) (2 (Subscript[a, 1,0]) (Subscript[a, 1,2]) (Subscript[a, 1,3])+(Subscript[a, 1,1]) ((Subscript[a, 1,2]^2)-(Subscript[a, 1,3]^2)))))));
facN=Plus@@Table[bezK[[i+1]]bernstein[t,2,i],{i,0,2}];
fac=facN/theta^4;
fac0=bezK[[1]]/ξ^4;
fac1=bezK[[3]]/ξ^4;
Return[{fac0,fac1,facN,theta}]
]

ClearAll[G1InterpolantLength]
(* Function returns (Subscript[ψ, 0(Subscript[,,  ])])(Subscript[ψ, 1, ])numerator(ψ), theta, length of r. *)
G1InterpolantLength[P1_,d0_,d1_,phi_,xi_]:=Module[{ξ,φ,d,P},
ξ=xi;
φ=phi;
{(Subscript[d, 0,1]),(Subscript[d, 0,2]),(Subscript[d, 0,3])}=d0;
{(Subscript[d, 1,1]),(Subscript[d, 1,2]),(Subscript[d, 1,3])}=d1;
{(Subscript[P, 1]),(Subscript[P, 2]),(Subscript[P, 3])}=P1;
Return[(((Subscript[P, 1]) (1+(Subscript[d, 1,1]))+(Subscript[P, 3]) (-Cot[φ] (Subscript[d, 1,2])+(Subscript[d, 1,3]))+(Subscript[P, 2]) ((Subscript[d, 1,2])+Cot[φ] (Subscript[d, 1,3])))/(1+(Subscript[d, 1,1])))];
]

ClearAll[G1InterpolantGen]
(* Function returns class 4 G1 interpolant for general data P0,P1,d0,d1. Free angle is determined so that the length of the interpolant is equal to scalFac*|P1-P0|. Parameter Izbxi determines the value of xi. If Izbxi=0, we choose the left boundary value + 10^(-10), if Izbxi=1/2 we choose the middle value in the admissible interval, else we choose the boundary value. 'num' determines which solution (r_+ or r_-) is chosen. *)
G1InterpolantGen[t_,P0_,P1_,d0_,d1_,scalFac_,Izbxi_,num_,disp_]:=Module[{QQR,QQRi,P0t,P1t,d0t,d1t,φ,ξ,pom,dolz,izbFi,intKsi,PP1,RR1},
(* rotation *)
If[Norma[d0]≠1, d0=d0/Norma[d0]];
If[Norma[d1]≠1, d1=d1/Norma[d1]];
QQR=SolveEQStar[d0,0];
QQRi=KvaternionConj[QQR]/Norm[QQR];
(* transformed data *)
P0t={0,0,0};
P1t=Rotacija[QQRi][[3]].(P1-P0);
d0t={1,0,0};
d1t=Rotacija[QQRi][[3]].d1;
(* length of the interpolant *)
dolz=G1InterpolantLength[P1t,d0t, d1t,φ,ξ];
pom=FindRoot[(dolz-scalFac(Norma[P1-P0]))==0,{φ,Pi/2}];
izbFi=φ/.pom[[1]];
intKsi=IntervalForKsi[P1t,d0t,d1t,izbFi];
If[Izbxi==0,PP1=Expand[G1Interpolant[t,P1t,d0t,d1t,izbFi,intKsi[[1]]+10^(-10),num]]];
If[Izbxi==1/2,PP1=Expand[G1Interpolant[t,P1t,d0t,d1t,izbFi,(intKsi[[1]]+intKsi[[2]])/2,num]]];
If[(Izbxi≠0)&&(Izbxi!=1/2),PP1=Expand[G1Interpolant[t,P1t,d0t,d1t,izbFi,intKsi[[2]]-10^(-10),num]]];
RR1=Rotacija[QQR][[3]].(Table[PP1[[i]]/PP1[[1]],{i,2,4}])+P0;
If[disp==1,Print["Check the solution: ",Norm[{P0-RR1/.{t->0},P1-RR1/.{t->1},d0-D[RR1,t]/Norma[D[RR1,t]]/.{t->0},d1-D[RR1,t]/Norma[D[RR1,t]]/.{t->1}}]]];
Return[RR1];
]

ClearAll[BendEnergy]
(* Function computes the bend energy.  *)
BendEnergy[r_,t_,a_,b_,WP_]:=Module[{u1,u2,ii},
u1=Norm[Cross[D[r,t],D[r,{t,2}]]];
u2=Norm[D[r,t]];
ii=NIntegrate[u1^2/u2^5,{t,a,b},WorkingPrecision->WP];
Return[ii]
]

ClearAll[BendEnergyApprox]
(* Function computes the approximated bend energy.  *)
BendEnergyApprox[r_,t_,a_,b_,WP_]:=Module[{u1,u2,ii},
u1=Norm[D[r,{t,2}]];
ii=NIntegrate[u1^2,{t,a,b},WorkingPrecision->WP];
Return[ii]
]

ClearAll[PlotRegionBetweenTwoCurves]
PlotRegionBetweenTwoCurves[f1_,f2_,t_,a_,b_,h_]:=
Polygon[Join[Table[{x,f1/.{t->x}},{x,a,b,h}],Reverse[Table[{x,f2/.{t->x}},{x,a,b,h}]]]];

Derivation and properties of class 4 rational PH curves

General derivation

Print["Quadratic quaternion polynomial A: "]
ClearAll[a,ξ]
KA={{(Subscript[a, 0,0]),(Subscript[a, 0,1]),(Subscript[a, 0,2]),(Subscript[a, 0,3])},(Subscript[ξ, 0 ]){(Subscript[a, 0,0]),(Subscript[a, 0,1]),(Subscript[a, 0,2]),(Subscript[a, 0,3])}+ (Subscript[ξ, 1 ]){(Subscript[a, 1,0]),(Subscript[a, 1,1]),(Subscript[a, 1,2]),(Subscript[a, 1,3])},{(Subscript[a, 1,0]),(Subscript[a, 1,1]),(Subscript[a, 1,2]),(Subscript[a, 1,3])}};
AA=Plus@@Table[KA[[i+1]]bernstein[t,2,i],{i,0,2}]

Print["Compute the ER frame vectors h:=e1=(SuperStar[AiA]),e2=(SuperStar[AjA]),e3=(SuperStar[AkA]) and a binormal v=h×h':"]
{WW1,WW2,WW3,WI,EE1,EE2,EE3,UN}=Simplify[IzracunERF[AA,t]];
Print["h = ", Simplify[EE1]]
Print["v = ", Simplify[UN]]

theta=Collect[((t^2) (Subscript[ξ, 0])+(-1+t) (-t+(-1+t) (Subscript[ξ, 1]))),t];
Print["Polynomial θ(t) = ", theta]

Print["Reduced binormal (Subscript[v, R]):"]
UR=Simplify[Simplify[PolynomialQuotient[WW2,theta,t]]EE2+Simplify[PolynomialQuotient[WW3,theta,t]]EE3]
Simplify[UR-Table[Simplify[PolynomialQuotient[UN[[i]],theta,t]],{i,1,3}]]

Print["Dual form: "]
ClearAll[f]
ff=Plus@@Table[(Subscript[f, i]) t^i,{i,0,4}];
L={ff , UR[[1]], UR[[2]], UR[[3]]};
DL=D[L,t];
D2L=D[L,{t,2}];
(* Wedge product *)
{D0,D1,D2,D3}=WedgeProd[L,DL,D2L];
Print["Rational PH curve:"]
RR={D1/D0,D2/D0,D3/D0}

Reparameterization

Print["Quadratic quaternion polynomial A: "]
ClearAll[a,ξ]
KA={{(Subscript[a, 0,0]),(Subscript[a, 0,1]),(Subscript[a, 0,2]),(Subscript[a, 0,3])},(Subscript[ξ, 0 ]){(Subscript[a, 0,0]),(Subscript[a, 0,1]),(Subscript[a, 0,2]),(Subscript[a, 0,3])}+ (Subscript[ξ, 1 ]){(Subscript[a, 1,0]),(Subscript[a, 1,1]),(Subscript[a, 1,2]),(Subscript[a, 1,3])},{(Subscript[a, 1,0]),(Subscript[a, 1,1]),(Subscript[a, 1,2]),(Subscript[a, 1,3])}}
AA=Plus@@Table[KA[[i+1]]bernstein[t,2,i],{i,0,2}]

Print["Compute the ER frame vectors h:=e1=(SuperStar[AiA]),e2=(SuperStar[AjA]),e3=(SuperStar[AkA]) and a binormal v=h×h':"]
{WW1,WW2,WW3,WI,EE1,EE2,EE3,UN}=Simplify[IzracunERF[AA,t]];
hh=Simplify[EE1];
vv=Simplify[UN];
Print["h = ",hh]
Print["v = ", vv]
Print["θ(t) = "]
theta=Collect[((t^2) (Subscript[ξ, 0])+(-1+t) (-t+(-1+t) (Subscript[ξ, 1]))),t]

Print["Reduced binormal (Subscript[v, R]):"]
UR=Simplify[Simplify[PolynomialQuotient[WW2,theta,t]]EE2+Simplify[PolynomialQuotient[WW3,theta,t]]EE3]
Simplify[UR-Table[Simplify[PolynomialQuotient[vv[[i]],theta,t]],{i,1,3}]]

Print["Dual form: "]
ClearAll[f]
ff=Plus@@Table[(Subscript[f, i]) t^i,{i,0,4}];
L={-ff , UR[[1]], UR[[2]], UR[[3]]};
DL=D[L,t];
D2L=D[L,{t,2}];
(* Wedge product *)
{D0,D1,D2,D3}=WedgeProd[L,DL,D2L];
Print["Rational PH curve:"]
RR={D1/D0,D2/D0,D3/D0}

Print["Let us reparameterize the quaternion polynomial with a reparameterization t->ρ(1-u)+u: "]
rep={t->ρ(1-u)+u}

Print["Control points after the reparameterization:"]
KArep=Simplify[BezKoef[2,CoefficientList[#,u]]]&/@(AA/.rep)

Print["Let us introduce new quaternion C: "]
pomC=Flatten[Solve[{(Subscript[c, 0]),(Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3])}==Transpose[KArep][[1]],{(Subscript[c, 0]),(Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3])}]]//Simplify
pomA0=First[Simplify[Solve[Transpose[KArep][[1]]=={(Subscript[c, 0]),(Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3])},KA[[1]]]]]

Print["Determine ρ so that the weights will be equal: "]
pom3=CoefficientList[Transpose[KArep][[2,1]]/.pomA0,{(Subscript[a, 1,0]),(Subscript[c, 0])}];
Solve[pom3[[1,2]]==pom3[[2,1]],{ρ}]
resRho=Last[Solve[pom3[[1,2]]==pom3[[2,1]],{ρ}]]

Print["Connection between A0 and C at a particularly chosen reparameterization: "]
pomCa=FullSimplify[pomC/.resRho]
pomA0a=Simplify[pomA0/.resRho]

Print["Let's check: "]
Simplify[KArep/.pomA0/.resRho]//MatrixForm

Print["Let us introduce new ξ:"]
pomKi={ξξ->-(((Subscript[ξ, 0]) (-1+(Subscript[ξ, 0])+(Subscript[ξ, 1])) )/(1+2 (Subscript[ξ, 0]^2)-2 (Subscript[ξ, 0]) (1+(Subscript[ξ, 1]))))}
pomKi1=Flatten[Simplify[Solve[-(((Subscript[ξ, 0]) (-1+(Subscript[ξ, 0])+(Subscript[ξ, 1])) )/(1+2 (Subscript[ξ, 0]^2)-2 (Subscript[ξ, 0]) (1+(Subscript[ξ, 1]))))==ξξ,(Subscript[ξ, 1])]]]
Simplify[KArep/.pomA0/.resRho/.pomKi1]//MatrixForm
Print["The weights are equal!"]

Print["Connection between A0 and C at a particularly chosen reparameterization and new ξξ: "]
pomCb=FullSimplify[pomCa/.pomKi1]
pomA0b=Simplify[pomA0a/.pomKi1]

Print["Derived reparameterization: "]
repN=rep/.resRho
repNInv=Simplify[First[Solve[t==(t/.repN),u]]]
odvRep=Simplify[D[(t/.repN),u]];

Print["New quaternion polynomial: "]
KA1={{(Subscript[c, 0]),(Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3])},ξξ{(Subscript[c, 0]),(Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3])}+ ξξ{(Subscript[a, 1,0]),(Subscript[a, 1,1]),(Subscript[a, 1,2]),(Subscript[a, 1,3])},{(Subscript[a, 1,0]),(Subscript[a, 1,1]),(Subscript[a, 1,2]),(Subscript[a, 1,3])}}
AA1=Plus@@Table[KA1[[i+1]]bernstein[u,2,i],{i,0,2}]

Print["Let's check: "]
Simplify[AA1/.pomKi]-Simplify[AA/.repN/.pomA0a]

Print["Compute the ER frame vectors and a binormal for a reparameterized A:"]
{WW11,WW21,WW31,WI1,EE11,EE21,EE31,UN1}=Simplify[IzracunERF[AA1,u]];
Print["Reparameterized θ:"]
theta1=Simplify[theta/.repN/.pomKi1]
Print["Reduced binormal:"]
hh1=Cross[EE11,D[EE11,u]];
UR1=Simplify[PolynomialQuotient[#,theta1,t]&/@hh1];

Print["New f:"]
ClearAll[fu]
resF=Simplify[First[Solve[Simplify[CoefficientList[Plus@@Table[(Subscript[f, i]) t^i,{i,0,4}]/.repN,u]]==Table[(Subscript[fu, i]),{i,0,4}],Table[(Subscript[fu, i]),{i,0,4}]]]/.pomKi1]

(* Dual form *)
CleaAll[fu,u]
ff1=Plus@@Table[(Subscript[fu, i]) u^i,{i,0,4}];
(*Simplify[ff-Simplify[ff1/.repNInv/.resF/.pomKi]]*)
L1={-ff1 odvRep , UR1[[1]], UR1[[2]], UR1[[3]]};
DL1=D[L1,u];
D2L1=D[L1,{u,2}];
(* Wedge produkt *)
{D01,D11,D21,D31}=WedgeProd[L1,DL1,D2L1];
Print["Curve:"]
RR1={D11/D01,D21/D01,D31/D01};

(* --- example --- *)

Print["Choesen data: "]
pod={(Subscript[a, 0,0])->1,(Subscript[a, 0,1])->0,(Subscript[a, 0,2])->0,(Subscript[a, 0,3])->0,(Subscript[a, 1,0])->2,(Subscript[a, 1,1])->-1,(Subscript[a, 1,2])->3,(Subscript[a, 1,3])->-2,(Subscript[ξ, 0])->1/5,(Subscript[ξ, 1])->-1/12,(Subscript[f, 0])->0,(Subscript[f, 1])->0,(Subscript[f, 2])->0,(Subscript[f, 3])->2,(Subscript[f, 4])->1}
Print["Curve: "]
kriv1=Simplify[RR/.pod]

Print["New data:"]
podc=pomC/.resRho/.pod
podKi=pomKi/.pod
podF=resF/.pod/.podKi

Print["Curve from the new data: "]
kriv2=Simplify[Simplify[RR1/.podc/.pod/.podKi/.podF]/.(pomKi/.pod)]

Print["Let's check if the new curve is a reparameterized old curve: "]
Simplify[Simplify[kriv1/.repN/.pod]-kriv2]

Show[{ParametricPlot3D[kriv1,{t,-5,0},PlotStyle->RGBColor[1,0,0]],ParametricPlot3D[kriv2,{u,-5,-1},PlotRange->All]}]

Properties

Print["Quadratic quaternion polynomial A after a reparameterization and rotation: "]
ClearAll[a,ξ]
KA={{(Subscript[a, 0,0]),(Subscript[a, 0,1]),(Subscript[a, 0,2]),(Subscript[a, 0,3])},(Subscript[ξ, 0 ]){(Subscript[a, 0,0]),(Subscript[a, 0,1]),(Subscript[a, 0,2]),(Subscript[a, 0,3])}+ (Subscript[ξ, 1 ]){(Subscript[a, 1,0]),(Subscript[a, 1,1]),(Subscript[a, 1,2]),(Subscript[a, 1,3])},{(Subscript[a, 1,0]),(Subscript[a, 1,1]),(Subscript[a, 1,2]),(Subscript[a, 1,3])}}/.{(Subscript[ξ, 0 ])->ξ,(Subscript[ξ, 1 ])->ξ,(Subscript[a, 0,2])->0,(Subscript[a, 0,3])->0}
AA=Plus@@Table[KA[[i+1]]bernstein[t,2,i],{i,0,2}]

Print["Compute the ER frame vectors h:=e1=(SuperStar[AiA]),e2=(SuperStar[AjA]),e3=(SuperStar[AkA]) and a binormal v=h×h':"]
{WW1,WW2,WW3,WI,EE1,EE2,EE3,UN}=Simplify[IzracunERF[AA,t]];
hh=Simplify[EE1];
vv=Simplify[UN];
Print["h = ",hh]
Print["v = ", vv]
Print["θ(t) = "]
theta=Collect[((t^2) ξ+(-1+t) (-t+(-1+t) ξ)),t]

Print["Reduced binormal (Subscript[v, R]):"]
UR=Simplify[Simplify[PolynomialQuotient[WW2,theta,t]]EE2+Simplify[PolynomialQuotient[WW3,theta,t]]EE3]
Simplify[UR-Table[Simplify[PolynomialQuotient[vv[[i]],theta,t]],{i,1,3}]]

Print["Dual form: "]
ClearAll[f]
ff=Plus@@Table[(Subscript[f, i]) t^i,{i,3,4}];
L={ff , UR[[1]], UR[[2]], UR[[3]]};
DL=D[L,t];
D2L=D[L,{t,2}];
(* Wedge product *)
{D0,D1,D2,D3}=WedgeProd[L,DL,D2L];
Print["Rational PH curve:"]
RR={D1/D0,D2/D0,D3/D0};

Print["Curve's denominator:"]
imen=Simplify[D0]

Print["Factor psi: "]
fac=Simplify[FactorPSI[-ff theta,1,EE1,t]]

Print["Poles: "]
poles=Simplify[Solve[theta==0,t]]

Print["Zeros of θ correspond to points on A where a tangent line passes through the origin: "]
koef=Simplify[AA[[1]]/D[AA[[1]],t]]
Simplify[Simplify[AA-koef D[AA,t]]/.poles]

Print["Irregular points:"]
Simplify[Solve[fac==0,t]]

(* -------- The curve can't have nodes --------------- *)

Print["Suppose that at t=t1 and t=t2 the curve has a double point. Then the following equations must hold true: "]
eq=Simplify[RR/.{t->t1}]-Simplify[RR/.{t->t2}];

Print["Solution: "]
Solve[{eq[[2]],eq[[3]]}=={0,0},{(Subscript[f, 3]),(Subscript[f, 4])}]
RR/.{(Subscript[f, 3])->0,(Subscript[f, 4])->0}
Print["This solution gives r=0."]

Print["Let check this using Groebner Basis: "]
G=Simplify[GroebnerBasis[eq,{(Subscript[f, 3]),(Subscript[f, 4]),t1}]];
Length[G]

Print["Element G(9) implies (Subscript[f, 4])=0 or ξ=0: "]
G[[9]]

Print["Suppose first that (Subscript[f, 4])=0. Then the element G(14) equals "]
Simplify[G[[14]]/.{(Subscript[f, 4])->0}]
Print["It is zero if (Subscript[f, 3])=0, t1=t2 or"]
res=Simplify[Solve[Simplify[G[[14]]/.{(Subscript[f, 4])->0}]==0,ξ]]
Print["The element G(18): "]
Simplify[G[[18]]/.{(Subscript[f, 4])->0}/.res[[1]]]
Simplify[G[[18]]/.{(Subscript[f, 4])->0}/.res[[2]]]
Simplify[G[[18]]/.{(Subscript[f, 4])->0}/.res[[3]]]
Simplify[G[[18]]/.{(Subscript[f, 4])->0}/.res[[4]]]
Print["This implies (Subscript[f, 3])=0. "]

Print["Suppose now that ξ=0. Then the Groebner basis simplifies to "]
Simplify[G/.{ξ->0}]
Print["This implies also (Subscript[f, 3])=0. But the curve in this case reduces to a point: "]
Simplify[RR/.{ξ->0,(Subscript[f, 3])->0}]

(* ----------------------------------------- *)
Print["Compute length of a curve."]
der=fac Norma[AA]^2;
ClearAll[c]
(* indefinite integral *)
Ninteg=Plus@@Table[(Subscript[c, i])t^i,{i,1,6}]/theta^3;
NDinteg=Simplify[D[Ninteg,t]];
polp6=der*theta^4;
r1=Simplify[CoefficientList[Together[polp6],t]];
r2=CoefficientList[Together[NDinteg*theta^4],t];
(* Solution using matrices *)
(*Print["Matrix 7×6: "]*)
SS=Transpose[Table[r2/.Join[Table[(Subscript[c, i])->0,{i,1,k-1}],{(Subscript[c, k])->1},Table[(Subscript[c, i])->0,{i,k+1,6}]],{k,1,6}]];
Simplify[(SS.Table[(Subscript[c, i]),{i,1,6}])-r2];
(*Print["Submatrix, where the last equation is left out: "]*)
SS6=Drop[SS,{7}];
(*Print["Its determinant: "]*)
Det[SS6];
(*Print["Unknown koefficents c_i, i=1,...,6"]*)
resC6=Simplify[Inverse[Drop[SS,{7}]].Drop[r1,{7}]];
(*Print["Check the last equation: "]*)
Simplify[Take[SS,{7}].resC6-Take[r1,{7}]];
resC6a=Table[(Subscript[c, i])->resC6[[i]],{i,1,6}];
arcLength=Ninteg/.resC6a;
Print["Length of the curve for t in [0,1]:"]
arcLength1=Simplify[arcLength/.{t->1}]

Print["Coefficients c_i expressed with given cc_i: "]
(Simplify[Inverse[SS6].Table[(Subscript[cc, i]),{i,0,5}]])
resC61=(Simplify[Inverse[SS6].Table[(Subscript[cc, i]),{i,0,5}]])/.Table[(Subscript[cc, i])->r1[[i+1]],{i,0,6}];
resC6b=Table[(Subscript[c, i])->resC61[[i]],{i,1,6}];
Print["Lets check."]
Simplify[resC6-resC61]
(r1-r2)/.Simplify[Flatten[Simplify[Solve[(#==0)&/@Take[Table[(Subscript[cc, i]),{i,0,6}]-r2,{1,6}],Table[(Subscript[c, i]),{i,1,6}]]]]/.Table[(Subscript[cc, i])->r1[[i+1]],{i,0,6}]]//Simplify

Print["If ksi=0:"]
pom=Simplify[(SS/.{ξ->0}).Table[(Subscript[c, i]),{i,1,6}]-Simplify[r1/.{ξ->0}/.{(Subscript[a, 0,1])->0,(Subscript[a, 0,2])->0,(Subscript[a, 0,3])->0}]];
pom1=Flatten[Solve[(#==0)&/@{pom[[3]],pom[[5]],pom[[6]],pom[[7]]},Table[(Subscript[c, i]),{i,2,5}]]]/.{(Subscript[c, 6])->0,(Subscript[c, 1])->0}
Simplify[pom/.pom1]/.{(Subscript[c, 6])->0,(Subscript[c, 1])->0}

Example

Print["Quadratic quaternion polynomial: "]
ClearAll[a,ξ]
KA={{1,0,0,0},ξ{1,0,0,0}+ ξ{1/2,-1,2,1/4},{1/2,-1,2,1/4}}/.{ ξ->1};
AA=Plus@@Table[KA[[i+1]]bernstein[t,2,i],{i,0,2}]
Print["Compute the ER frame vectors h:=e1=(SuperStar[AiA]),e2=(SuperStar[AjA]),e3=(SuperStar[AkA]) and a binormal v=h×h':"]
{WW1,WW2,WW3,WI,EE1,EE2,EE3,UN}=Simplify[IzracunERF[AA,t]];
hh=Simplify[EE1];
vv=Simplify[UN];
Print["h = ",hh]
Print["v = ", vv]
Print["θ(t) = "]
theta=Collect[((t^2) ξ+(-1+t) (-t+(-1+t) ξ)),t]/.{ ξ->1}

Print["Reduced binormal: "]
UR=Simplify[Simplify[PolynomialQuotient[WW2,theta,t]]EE2+Simplify[PolynomialQuotient[WW3,theta,t]]EE3]
(* Dual form *)
ff=3 t^3+5 t^4
L={-ff , UR[[1]], UR[[2]], UR[[3]]};
DL=Simplify[D[L,t]];
D2L=Simplify[D[L,{t,2}]];
(* Wedge produkt *)
PP=WedgeProd[L,DL,D2L];
Print["Curve:"]
RR=Simplify[Table[PP[[i]]/PP[[1]],{i,2,4}]]

ParametricPlot3D[RR,{t,0,1},PlotRange->All,Axes->True,BoxStyle->GrayLevel[0.8],PlotStyle->GrayLevel[0]]

Print["Factor psi: "]
fac=Simplify[FactorPSI[ff theta,1,EE1,t]]
Print["Irregular points obtained at parameters: "]
Solve[fac==0,t]
%//N

Print["Compute length of a curve:"]
der=-Simplify[fac Norma[AA]^2];
pom=Simplify[Integrate[der,{t,0,s}]];
Print["Arc-length:"]
arcLength=Simplify[Expand[pom]]
Print["Length of the curve for t in [0,1]:"]
arcLength1=Simplify[arcLength/.{s->1}]

(G^1) interpolation - derivation of equations

Print["Data: "]
ClearAll[P,d]
P0={0,0,0}
P1={(Subscript[P, 1]),(Subscript[P, 2]),(Subscript[P, 3])}
d0={(Subscript[d, 0,1]),(Subscript[d, 0,2]),(Subscript[d, 0,3])}
d0={1,0,0}
d1={(Subscript[d, 1,1]),(Subscript[d, 1,2]),(Subscript[d, 1,3])}

Print["Quadratic quaternion polynomial: "]
ClearAll[a,ξ]
KA={{(Subscript[a, 0,0]),(Subscript[a, 0,1]),(Subscript[a, 0,2]),(Subscript[a, 0,3])},ξ{(Subscript[a, 0,0]),(Subscript[a, 0,1]),(Subscript[a, 0,2]),(Subscript[a, 0,3])}+ ξ{(Subscript[a, 1,0]),(Subscript[a, 1,1]),(Subscript[a, 1,2]),(Subscript[a, 1,3])},{(Subscript[a, 1,0]),(Subscript[a, 1,1]),(Subscript[a, 1,2]),(Subscript[a, 1,3])}};
AA=Plus@@Table[KA[[i+1]]bernstein[t,2,i],{i,0,2}]
Print["Compute the ER frame vectors h:=e1=(SuperStar[AiA]),e2=(SuperStar[AjA]),e3=(SuperStar[AkA]) and a binormal v=h×h':"]
{WW1,WW2,WW3,WI,EE1,EE2,EE3,UN}=Simplify[IzracunERF[AA,t]];
theta=Collect[(t+ξ-2 t ξ+(t^2) (-1+2 ξ)),t];
Print["Polynomial θ(t) = ", theta]

Print["Reduced binormal: "]
UR=Simplify[Simplify[PolynomialQuotient[WW2,theta,t]]EE2+Simplify[PolynomialQuotient[WW3,theta,t]]EE3];
(* Dual form *)
ClearAll[f]
ff=(Subscript[f, 3]) t^3+(Subscript[f, 4]) t^4;
L={-ff , UR[[1]], UR[[2]], UR[[3]]};
DL=Simplify[D[L,t]];
D2L=Simplify[D[L,{t,2}]];
(* Wedge produkt *)
PP=WedgeProd[L,DL,D2L];
Print["Curve:"]
RR=Table[PP[[i]]/PP[[1]],{i,2,4}];

Print["Factor psi: "]
fac=FactorPSI[ff theta,1,EE1,t];
Simplify[fac/.{(Subscript[a, 0,2])->0,(Subscript[a, 0,3])->0}]
Print["Factor psi * θ^4(t): "]
facBT=Simplify[Together[fac theta^4]];
Simplify[facBT/.{(Subscript[a, 0,2])->0,(Subscript[a, 0,3])->0}]
Print["Its numerator and denominator: "]
facNum=-Numerator[Together[facBT]];
facDenom=-Simplify[Denominator[Together[facBT]]];
Simplify[facNum/.{(Subscript[a, 0,2])->0,(Subscript[a, 0,3])->0}]
Simplify[facDenom/.{(Subscript[a, 0,2])->0,(Subscript[a, 0,3])->0}]
Print["Bezier coefficients of factor's numerator:"]
facBezCoef=Simplify[BezKoef[2,CoefficientList[facNum,t]]]

Print["Interpolation of tangent direction at t=0 -> determine A0 "]
ClearAll[φ]
pom=Simplify[Simplify[SolveEQStar[d0,(Subscript[φ, 0])]]];
resA0=Table[KA[[1,i]]->pom[[i]],{i,1,4}];
Print["A0: ",resA0]
normaA0=Simplify[Simplify[Norma[pom]]];
Print["||A0|| = ", normaA0]

Print["Interpolation of tangent direction at t=1 -> determine A2 "]
pom=Simplify[λ FullSimplify[SolveEQStar[d1,(Subscript[φ, 1])]]/.{d1[[2]]^2->1-d1[[1]]^2-d1[[3]]^2}];
resA2=Table[KA[[3,i]]->pom[[i]],{i,1,4}];
Print["A2: ", resA2]
normaA2=Simplify[Simplify[Norma[pom]]/.{d1[[1]]^2+d1[[2]]^2+d1[[3]]^2->1}];
Print["||A2|| = ",normaA2]

Print["ψ(0)= "]
fac0=Simplify[Simplify[fac/.{t->0}/.resA0/.resA2]]
Print["r'(0) = "]
Simplify[D[RR,t]/.{t->0}/.resA0/.resA2]
Print["ψ(0) (Subscript[A, 0]) i (Subscript[A, 0]^( * )) = "]
Simplify[fac0 Rest[OPStar[pom,pom]]]

Print["ψ(1)= "]
fac1=Simplify[Simplify[fac/.{t->1}/.resA0]/.resA2]
Print["ψ(1) (Subscript[A, 2]) i (Subscript[A, 2]^( * )) = "]
Simplify[fac1 Rest[OPStar[pom,pom]]/.{((Subscript[d, 1,2])^2)->1-((Subscript[d, 1,1])^2)-((Subscript[d, 1,3])^2)}]
Print["(Subscript[A, 2]) i (Subscript[A, 2]^( * )) = "]
Simplify[Rest[OPStar[pom,pom]]/.{((Subscript[d, 1,2])^2)->1-((Subscript[d, 1,1])^2)-((Subscript[d, 1,3])^2)}]

Print["Interpolation of a point at t=0 - fulfilled:"]
eqP0=Simplify[(RR/.{t->0})-P0]

Print["Interpolation of a point at t=1. Equations:"]
eP1=Simplify[Simplify[(RR/.{t->1})-P1]/.resA0/.resA2];
(* equations for point interpolation *)
Print["Denominator: "]
e0=Denominator[Together[eP1[[2]]]]Sin[(Subscript[φ, 0])-(Subscript[φ, 1])](Sqrt[2])/2
Print["Numerator:"]
{e1,e2,e3}=Simplify[Simplify[eP1*e0]/.{d1[[2]]^2->1-d1[[1]]^2-d1[[3]]^2}]

Print["New variables:"]
ClearAll[X,Y]
{-(Subscript[f, 3])+4 ((Subscript[f, 3])+(Subscript[f, 4])) ξ->-X,(-1+4 ξ-4 (ξ^2)+8 (ξ^3)) (Subscript[f, 3])+4 ξ (1+2 (ξ^2)) (Subscript[f, 4])->-Y}
nove1=Join[Flatten[Solve[{-(Subscript[f, 3])+4 ((Subscript[f, 3])+(Subscript[f, 4])) ξ==-X,(-1+4 ξ-4 (ξ^2)+8 (ξ^3)) (Subscript[f, 3])+4 ξ (1+2 (ξ^2)) (Subscript[f, 4])==-Y},{(Subscript[f, 3]),(Subscript[f, 4])}]],{(Subscript[φ, 0])->(Subscript[φ, 1])+φ}]//Simplify

Print["Equations in new variables: "]
e1n=Collect[FullSimplify[e1/.nove1],{X,Y}]
e2n=Collect[Simplify[e2/.nove1],{X,Y}]
e3n=Collect[Simplify[e3/.nove1],{X,Y}]

Print["From the second and the third equation we compute X and Y."]
resXY=Simplify[FullSimplify[Flatten[Simplify[Solve[{e2n,e3n}=={0,0},{X,Y}]]]]/.{d1[[2]]^2->1-d1[[1]]^2-d1[[3]]^2}]

Print["We introduce new constants to simplify the notation. "]
constCCinv={(Subscript[CC, 1])->((Subscript[P, 3]) (Subscript[d, 1,2])-(Subscript[P, 2]) (Subscript[d, 1,3])),(Subscript[CC, 2])->((Subscript[P, 2]) (Subscript[d, 1,2])+(Subscript[P, 3]) (Subscript[d, 1,3])),(Subscript[CC, 3])->(Subscript[P, 1]) (-1+(Subscript[d, 1,1]^2))+(Subscript[d, 1,1]) ((Subscript[P, 2]) (Subscript[d, 1,2])+(Subscript[P, 3]) (Subscript[d, 1,3]))}
constCC=Simplify[Flatten[Solve[{((Subscript[P, 3]) (Subscript[d, 1,2])-(Subscript[P, 2]) (Subscript[d, 1,3]))==(Subscript[CC, 1]),((Subscript[P, 2]) (Subscript[d, 1,2])+(Subscript[P, 3]) (Subscript[d, 1,3]))==(Subscript[CC, 2]),(Subscript[CC, 3])==(Subscript[P, 1]) (-1+(Subscript[d, 1,1]^2))+(Subscript[d, 1,1]) ((Subscript[P, 2]) (Subscript[d, 1,2])+(Subscript[P, 3]) (Subscript[d, 1,3]))},{(Subscript[P, 1]),(Subscript[P, 2]),(Subscript[P, 3])}]]/.{(Subscript[d, 1,2]^2)+(Subscript[d, 1,3]^2)->1-(Subscript[d, 1,1]^2)}]

Print["Solution X,Y:"]
pom=Simplify[{X,Y}/.resXY/.constCC];
resXY1=Simplify[{X->pom[[1]],Y->pom[[2]]}/.{(Subscript[d, 1,2]^2)+(Subscript[d, 1,3]^2)->1-(Subscript[d, 1,1]^2)}]
Print["Solution (Subscript[f, 3]),(Subscript[f, 4]):"]
resAB=FullSimplify[Drop[nove1,-1]/.resXY1]

Print["First equation after the substitution of X and Y:"]
e1nn=Simplify[Simplify[e1n/.resXY1]/(8 ξ  )]
Simplify[CoefficientList[e1nn,λ]]

Print["(Subscript[λ, 1]) expressed with (Subscript[λ, 1])^2 and a constant: "]
Lam1WithLam12=First[Simplify[Solve[(Simplify[CoefficientList[e1nn,λ]].{1,X,(λ^2)})=={0,0,0},X]]]

Print["Solution λ: "]
Simplify[Solve[e1nn==0,λ]];
resLam1=Simplify[Solve[e1nn==0,λ]/.constCCinv]

Admissible region - examples

Example 1 ((Subscript[C, 1])>0)

(*  example  data  *)
precision1=50;
P0p={0,0,0};
P1p={4,2,4};
d0p={1,0,0};
d1p={-6,3,2};
d1p=d1p/Norm[d1p];
Print["Data: P0 = ",P0p, ",  P1 = ", P1p , ",  d0 = ", d0p, ",  d1 = ", d1p]
G1=Graphics3D[{PointSize[0.02],Point[P0p],Point[P1p],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}],Dashing[0.01],Line[{P0p,P1p}]}];
Show[G1,Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Print["Constants (Subscript[C, 1]),(Subscript[C, 2]),(Subscript[C, 3]):"]
ConstC1C2C3[P1p,d1p]//N
Print["Norm of (Subscript[P, 1])-(Subscript[P, 0]):"]
N[Norm[P1p]]
Print["Length of the interpolant in dependence of phi: "]
dolz=G1InterpolantLength[P1p,d0p,d1p,φ,ξ];
N[Collect[dolz,Cot[φ]]]

Print["The admissible intervals for phi (for positive ξ):"]
intPhi=Last[IntervalForPhi[P1p,d0p,d1p,1]]
N[%]

Print["(Subscript[Z, m]), (Subscript[Z, M]):"]
{zm,zp}=N[IntervalForPhi[P1p,d0p,d1p,1][[3]]]

Print["The region of pairs (φ,ξ) for which the discriminant is positive: "]
g1=Show[{Graphics[{GrayLevel[0.8],Polygon[{{zm,-50},{zm,50},{zp,50},{zp,-50}}],PlotRegionBetweenTwoCurves[-IntervalForKsi[P1p,d0p,d1p,x][[2]],IntervalForKsi[P1p,d0p,d1p,x][[2]],x,0+0.0001,zm,1/100],PlotRegionBetweenTwoCurves[-IntervalForKsi[P1p,d0p,d1p,x][[2]],IntervalForKsi[P1p,d0p,d1p,x][[2]],x,zp+0.0001,Pi-0.0001,1/100]}],Graphics[{GrayLevel[0.7],Line[{{zm,-50},{zm,50}}],Line[{{zp,-50},{zp,50}}],GrayLevel[0.3],Point[{zm,0}],Point[{zp,0}],Text["(Subscript[Z, m])",{zm-0.06,-0.2}],Text["(Subscript[Z, M])",{zp-0.06,-0.2}]}]},Axes->True,PlotRange->{{0,Pi},{-2.5,2.5}},AxesLabel->{φ,ξ},AspectRatio->1,Ticks->{{0,Pi/4,Pi/2,3Pi/4,Pi},{0,1,2,3,4,-1,-2,-3,-4}},BaseStyle->{FontSize->16}]

Print["A subregion of admissible pairs (φ,ξ): "]
znak=Sign[ConstC1C2C3[P1p,d1p][[1]]];
Show[g1,Graphics[{GrayLevel[0.5],PlotRegionBetweenTwoCurves[znak IntervalForKsi[P1p,d0p,d1p,x][[1]],znak IntervalForKsi[P1p,d0p,d1p,x][[2]],x,zp+0.0001,Pi,1/100]}],Graphics[{GrayLevel[0.5],PlotRegionBetweenTwoCurves[-znak IntervalForKsi[P1p,d0p,d1p,x][[1]],-znak IntervalForKsi[P1p,d0p,d1p,x][[2]],x,0,zm-0.0001,1/100]}],Plot[znak IntervalForKsi[P1p,d0p,d1p,x][[1]],{x,zp,Pi-0.01},PlotRange->All,PlotStyle->{GrayLevel[0.3],Dashing[0.01]}],Plot[-znak IntervalForKsi[P1p,d0p,d1p,x][[1]],{x,0+0.01,zm},PlotRange->All,PlotStyle->{GrayLevel[0.3],Dashing[0.01]}],Plot[znak IntervalForKsi[P1p,d0p,d1p,x][[2]],{x,zp,Pi-0.01},PlotRange->All,PlotStyle->{GrayLevel[0.3]}],Plot[-znak IntervalForKsi[P1p,d0p,d1p,x][[2]],{x,0+0.01,zm},PlotRange->All,PlotStyle->{GrayLevel[0.3]}],AspectRatio->1,BaseStyle->{FontSize->16}]

Print["Table of interpolant's lengths for admissible φ:"]
N[Table[dolz,{φ,intPhi[[1]]+0.00001,intPhi[[2]]-0.00001,(intPhi[[2]]-intPhi[[1]])/50}]]

Show[{Plot[dolz,{φ,intPhi[[1]],intPhi[[2]]},PlotRange->{{-0.8,Pi},{-5,50}},AxesLabel->{φ,length},AspectRatio->1,Ticks->{{0,Pi/4,Pi/2,3Pi/4,Pi},{0,10,20,30,40,50}},BaseStyle->{FontSize->16}],Graphics[{GrayLevel[0.3],Point[{zm,0}],Point[{zp,0}],Text["(Subscript[Z, m])",{zm-0.1,-2}],Text["(Subscript[Z, M])",{zp-0.05,-2}],Text["\[LeftDoubleBracketingBar](Subscript[ΔP, 0])\[RightDoubleBracketingBar]",{-0.4,N[Norm[P1p]]}],Dashing[0.01],Line[{{0,Norm[P1p]},{Pi,Norm[P1p]}}],Line[{{zp,0},{zp,50}}]}]},BaseStyle->{FontSize->16},AspectRatio->1]

Example 2 ((Subscript[C, 1])<0)

precision1=50;
P0p={0,0,0};
P1p={4,2,4};
d0p={1,0,0};
d1p={2,-4,3};
d1p={-3,-4,3};
d1p={-6,-6,6};
d1p=d1p/Norm[d1p];
Print["Data: P0 = ",P0p, ",  P1 = ", P1p , ",  d0 = ", d0p, ",  d1 = ", d1p]
G1=Graphics3D[{PointSize[0.02],Point[P0p],Point[P1p],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}],Dashing[0.01],Line[{P0p,P1p}]}];
Show[G1,Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Print["Constants (Subscript[C, 1]),(Subscript[C, 2]),(Subscript[C, 3]):"]
ConstC1C2C3[P1p,d1p]//N
Print["Norm of (Subscript[P, 1])-(Subscript[P, 0]):"]
N[Norm[P1p]]
Print["Length of the interpolant in dependence of phi: "]
dolz=G1InterpolantLength[P1p,d0p,d1p,φ,ξ];
N[Collect[dolz,Cot[φ]]]

Print["The admissible intervals for phi (for positive ξ):"]
intPhi=Last[IntervalForPhi[P1p,d0p,d1p,1]];
N[%]

Print["(Subscript[Z, m]), (Subscript[Z, M]):"]
{zm,zp}=N[IntervalForPhi[P1p,d0p,d1p,1][[3]]]

Print["The region of pairs (φ,ξ) for which the discriminant is positive: "]
g1=Show[{Graphics[{GrayLevel[0.8],Polygon[{{zm,-50},{zm,50},{zp,50},{zp,-50}}],PlotRegionBetweenTwoCurves[-IntervalForKsi[P1p,d0p,d1p,x][[2]],IntervalForKsi[P1p,d0p,d1p,x][[2]],x,0+0.0001,zm,1/100],PlotRegionBetweenTwoCurves[-IntervalForKsi[P1p,d0p,d1p,x][[2]],IntervalForKsi[P1p,d0p,d1p,x][[2]],x,zp+0.0001,Pi-0.0001,1/100]}],Graphics[{GrayLevel[0.7],Line[{{zm,-50},{zm,50}}],Line[{{zp,-50},{zp,50}}],GrayLevel[0.3],Point[{zm,0}],Point[{zp,0}],Text["(Subscript[Z, m])",{zm+0.1,-0.2}],Text["(Subscript[Z, M])",{zp+0.1,-0.2}]}]},Axes->True,PlotRange->{{0,Pi},{-2.5,2.5}},AxesLabel->{φ,ξ},AspectRatio->1,Ticks->{{0,Pi/4,Pi/2,3Pi/4,Pi},{0,1,2,3,4,-1,-2,-3,-4}},BaseStyle->{FontSize->16}]

Print["A subregion of admissible pairs (φ,ξ): "]
znak=Sign[ConstC1C2C3[P1p,d1p][[1]]];
Show[g1,Graphics[{GrayLevel[0.5],PlotRegionBetweenTwoCurves[znak IntervalForKsi[P1p,d0p,d1p,x][[1]],znak IntervalForKsi[P1p,d0p,d1p,x][[2]],x,zp+0.0001,Pi,1/100]}],Graphics[{GrayLevel[0.5],PlotRegionBetweenTwoCurves[-znak IntervalForKsi[P1p,d0p,d1p,x][[1]],-znak IntervalForKsi[P1p,d0p,d1p,x][[2]],x,0,zm-0.0001,1/100]}],Plot[znak IntervalForKsi[P1p,d0p,d1p,x][[1]],{x,zp,Pi-0.01},PlotRange->All,PlotStyle->{GrayLevel[0.3],Dashing[0.01]}],Plot[-znak IntervalForKsi[P1p,d0p,d1p,x][[1]],{x,0+0.01,zm},PlotRange->All,PlotStyle->{GrayLevel[0.3],Dashing[0.01]}],Plot[znak IntervalForKsi[P1p,d0p,d1p,x][[2]],{x,zp,Pi-0.01},PlotRange->All,PlotStyle->{GrayLevel[0.3]}],Plot[-znak IntervalForKsi[P1p,d0p,d1p,x][[2]],{x,0+0.01,zm},PlotRange->All,PlotStyle->{GrayLevel[0.3]}],AspectRatio->1,BaseStyle->{FontSize->16}]

Print["Table of interpolant's lengths for admissible φ:"]
N[Table[dolz,{φ,intPhi[[1]]+0.00001,intPhi[[2]]-0.00001,(intPhi[[2]]-intPhi[[1]])/50}]]

Show[{Plot[dolz,{φ,intPhi[[1]],intPhi[[2]]},PlotRange->{{-0.7,Pi},{-5,50}},AxesLabel->{φ,length},AspectRatio->1,Ticks->{{0,Pi/4,Pi/2,3Pi/4,Pi},{0,10,20,30,40,50}},BaseStyle->{FontSize->16}],Graphics[{GrayLevel[0.3],Point[{zm,0}],Point[{zp,0}],Text["(Subscript[Z, m])",{zm+0.15,-2}],Text["(Subscript[Z, M])",{zp+0.1,-2}],Text["\[LeftDoubleBracketingBar](Subscript[ΔP, 0])\[RightDoubleBracketingBar]",{-0.4,N[Norm[P1p]]}],Dashing[0.01],Line[{{0,Norm[P1p]},{Pi,Norm[P1p]}}],Line[{{zm,0},{zm,50}}]}]},BaseStyle->{FontSize->16},AspectRatio->1]

(G^1) interpolation - examples with C1>0

Example 1

(*  example  data  *)
precision1=50;
P0p={0,0,0};
P1p={4,2,4};
d0p={1,0,0};
d1p={-6,3,2};
d1p=d1p/Norm[d1p];
Print["Data: P0 = ",P0p, ",  P1 = ", P1p , ",  d0 = ", d0p, ",  d1 = ", d1p]
G1=Graphics3D[{PointSize[0.02],Point[P0p],Point[P1p],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}],Dashing[0.01],Line[{P0p,P1p}]}];
Show[G1,Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Print["C1,C2,C3:"]
ConstC1C2C3[P1p,d1p]
Print["Norm of P1-P0:"]
N[Norm[P1p]]
Print["Length of the interpolant in dependence of phi: "]
dolz=G1InterpolantLength[P1p,d0p,d1p,φ,ξ];
N[Collect[dolz,Cot[φ]]]

Print["The admissible intervals for phi:"]
intPhi=Last[IntervalForPhi[P1p,d0p,d1p,1]]
N[%]

Print["Chosen phi: "]
izbFi=Pi/4
Print["Length of the interpolant for the chosen φ: "]
N[dolz/.{φ->izbFi},precision1]

Print["The admissible interval for xi based on the chosen phi:"]
intKsi=Simplify[IntervalForKsi[P1p,d0p,d1p,izbFi]]
N[%,8]

Print["Chosen ξ: "]
izbKsi=(intKsi[[1]]+intKsi[[2]])/2;
izbKsi=33/100

Print["First interpolant: "]
PP1=Expand[G1Interpolant[t,N[P1p,precision1],N[d0p,precision1],N[d1p,precision1],N[izbFi,precision1],N[izbKsi,precision1],1]];
RR1=Table[PP1[[i]]/PP1[[1]],{i,2,4}];
N[RR1]
G2a=ParametricPlot3D[RR1,{t,0,1},PlotStyle->GrayLevel[0.4]];
Show[{G1,G2a}]

Print["Second interpolant: "]
PP2=Expand[G1Interpolant[t,N[P1p,precision1],N[d0p,precision1],N[d1p,precision1],N[izbFi,precision1],N[izbKsi,precision1],2]];
RR2=Table[PP2[[i]]/PP2[[1]],{i,2,4}];
N[RR2]
G2b=ParametricPlot3D[RR2,{t,0,1},PlotStyle->GrayLevel[0.7]];
Show[{G1,G2b}]

Show[{G1,G2a,G2b},Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Example 1-negative ξ

(*  example  data  *)
precision1=50;
P0p={0,0,0};
P1p={4,2,4};
d0p={1,0,0};
d1p={-6,3,2};
d1p=d1p/Norm[d1p];
Print["Data: P0 = ",P0p, ",  P1 = ", P1p , ",  d0 = ", d0p, ",  d1 = ", d1p]
G1=Graphics3D[{PointSize[0.02],Point[P0p],Point[P1p],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}],Dashing[0.01],Line[{P0p,P1p}]}];
Show[G1,Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Print["C1,C2,C3:"]
ConstC1C2C3[P1p,d1p]
Print["Norm of P1-P0:"]
N[Norm[P1p]]
Print["Length of the interpolant in dependence of phi: "]
dolz=G1InterpolantLength[P1p,d0p,d1p,φ,ξ];
N[Collect[dolz,Cot[φ]]]

Print["The admissible intervals for phi (and negative ξ):"]
intPhi=Last[IntervalForPhi[P1p,d0p,d1p,-1]]
N[%]

Print["Chosen phi: "]
izbFi=3/10

Print["The admissible interval for xi based on the chosen phi:"]
intKsi=Sort[-Simplify[IntervalForKsi[P1p,d0p,d1p,izbFi]]]
N[%,8]

Print["Chosen ξ: "]
izbKsi=(intKsi[[1]]+intKsi[[2]])/2;
izbKsi=-105/100
N[%]

Print["First interpolant: "]
PP1=Expand[G1Interpolant[t,N[P1p,precision1],N[d0p,precision1],N[d1p,precision1],N[izbFi,precision1],N[izbKsi,precision1],1]];
RR1=Table[PP1[[i]]/PP1[[1]],{i,2,4}];
N[RR1]
G2a=ParametricPlot3D[RR1,{t,0,1},PlotStyle->GrayLevel[0.4]];
Show[{G1,G2a},Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

d0p-(D[RR1,t]/Norma[D[RR1,t]]/.{t->0})
d1p-(D[RR1,t]/Norma[D[RR1,t]]/.{t->1})

Print["Second interpolant: "]
PP2=Expand[G1Interpolant[t,N[P1p,precision1],N[d0p,precision1],N[d1p,precision1],N[izbFi,precision1],N[izbKsi,precision1],2]];
RR2=Table[PP2[[i]]/PP2[[1]],{i,2,4}];
N[RR2]
G2b=ParametricPlot3D[RR2,{t,0,1},PlotStyle->GrayLevel[0.7]];
Show[{G1,G2b},Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

d0p-(D[RR2,t]/Norma[D[RR2,t]]/.{t->0})
d1p-(D[RR2,t]/Norma[D[RR2,t]]/.{t->1})

Show[{G1,G2a,G2b},Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Example 1 - graphs of interpolants for different ξ at a fixed φ

(*  example  data  *)
precision1=40;
P0p=N[{0,0,0},precision1];
P1p=N[{4,2,4},precision1];
d0p=N[{1,0,0},precision1];
d1p=N[{-6,3,2},precision1];
d1p=d1p/Norm[d1p];

GData=Graphics3D[{PointSize[0.02],Point[P0p],Point[P1p],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}],Dashing[0.01],Line[{P0p,P1p}]}];
Show[GData,Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Print["Norm of P1-P0:"]
Norm[P1p]
Print["Length of the interpolant in dependence of phi: "]
dolz=G1InterpolantLength[P1p,d0p,d1p,φ,ξ]

Print["The admissible intervals for phi:"]
intPhi=Last[IntervalForPhi[P1p,d0p,d1p,1]]

Print["Chosen phi: "]
izbFi=Pi/4
dolz/.{φ->izbFi};

Print["The admissible interval for xi based on the chosen phi:"]
intKsi=Simplify[IntervalForKsi[P1p,d0p,d1p,izbFi]]

Print["Table of chosen values for xi: "]
tabKsi=Table[(1-i)intKsi[[1]]+i(intKsi[[2]]),{i,0,1,0.1}];
tabKsi=Table[intKsi[[1]]+i/9(intKsi[[2]]-intKsi[[1]]),{i,0,9}]
nn=Length[tabKsi];
tabKsi[[1]]=tabKsi[[1]]+10^(-7);
tabKsi[[nn]]=tabKsi[[nn]]-10^(-15);

Print["Table of first interpolants: "]
tabInt1={};
tabU={};
For[j=1,j≤Length[tabKsi],j++,
izb=tabKsi[[j]];
PP1=Expand[G1Interpolant[t,P1p,d0p,d1p,izbFi,izb,1]];
RR1=Table[PP1[[i]]/PP1[[1]],{i,2,4}];
tabInt1=Append[tabInt1,RR1];
tabU=Append[tabU,Ukriv[N[RR1],t]];
]
N[tabInt1];

G1=Show[{GData,Table[ParametricPlot3D[tabInt1[[i]],{t,0,1},PlotStyle->GrayLevel[0.3+(0.5/(1/nn-1))(i/nn-1)]],{i,1,nn}]},PlotRange->All,Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Print["Curvature: "]
G1U=Show[Table[Plot[tabU[[i]],{t,0,1},PlotStyle->GrayLevel[0.3+(0.5/(1/nn-1))(i/nn-1)]],{i,1,nn}],PlotRange->{All,{0,1}}]

Print["Bend energy:"]
tabBend1=Table[BendEnergy[tabInt1[[i]],t,0,1,20],{i,1,nn}]

Print["Approximated bend energy:"]
tabBendA1=Table[BendEnergyApprox[tabInt1[[i]],t,0,1,20],{i,1,nn}]

Print["Table of second interpolants: "]
tabInt2={};
tabU2={};
For[j=1,j≤Length[tabKsi],j++,
izb=tabKsi[[j]];
PP1=Expand[G1Interpolant[t,P1p,d0p,d1p,izbFi,izb,2]];
RR1=Table[PP1[[i]]/PP1[[1]],{i,2,4}];
tabInt2=Append[tabInt2,RR1];
tabU2=Append[tabU2,Ukriv[RR1,t]];
]
G2=Show[{GData,Table[ParametricPlot3D[tabInt2[[i]],{t,0,1},PlotStyle->GrayLevel[0.3+(0.5/(1/nn-1))(i/nn-1)]],{i,1,nn}]},PlotRange->All,Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]];

Print["Curvature: "]
G2U=Show[Table[Plot[tabU2[[i]],{t,0,1},PlotStyle->GrayLevel[0.3+(0.5/(1/nn-1))(i/nn-1)]],{i,1,nn}],PlotRange->{All,{0,1}}]

Print["Bend energy:"]
tabBend2=Table[BendEnergy[tabInt2[[i]],t,0,1,20],{i,1,nn}]

Print["Approximated bend energy:"]
tabBendA2=Table[BendEnergyApprox[tabInt2[[i]],t,0,1,20],{i,1,nn}]

Show[{G1,G2},Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Print["Interpolants at the boundary values for xi:"]
Show[{ParametricPlot3D[tabInt1[[1]],{t,0,1},PlotStyle->GrayLevel[0.6]],ParametricPlot3D[Last[tabInt1],{t,0,1},PlotStyle->GrayLevel[0]],ParametricPlot3D[tabInt2[[1]],{t,0,1},PlotStyle->GrayLevel[0.8]],ParametricPlot3D[Last[tabInt2],{t,0,1},PlotStyle->GrayLevel[0]],Graphics3D[{PointSize[0.02],Dashing[0.005],Point[P0p],Point[P1p],Line[{P0p,P1p}],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}]}]},PlotRange->All,Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Show[{G1U,G2U},Axes->{True,True},PlotRange->{{0,1},{0,1.5}}]

Transpose[Prepend[Transpose[N[{tabKsi,tabBend1,tabBend2,tabBendA1,tabBendA2}]],{"ξ","(SubPlus[E])","(SubMinus[E])","(SubPlus[EA])","(SubMinus[EA])"}]]//MatrixForm

Example 2 - graphs of interpolants for different φ

(*  example  data  *)
precision1=50;
P0p=N[{0,0,0},precision1];
P1p=N[{4,2,4},precision1];
d0p=N[{1,0,0},precision1];
d1p=N[{-6,3,2},precision1];
d1p=d1p/Norm[d1p];

GData=Graphics3D[{PointSize[0.02],Point[P0p],Point[P1p],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}],Dashing[0.01],Line[{P0p,P1p}]}];
Show[GData,Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Print["Norm of P1-P0:"]
Norm[P1p]
Print["Length of the interpolant in dependence of phi: "]
dolz=Simplify[G1InterpolantLength[P1p,d0p,d1p,φ,ξ]]

Print["The admissible intervals for phi:"]
intPhi=Last[IntervalForPhi[P1p,d0p,d1p,1]]

Print["Table of interpolants for different φ and ξ chosen at a boundary: "]
tabInt1={};
tabU={};
For[izbFi=intPhi[[1]]+10^(-10),izbFi≤intPhi[[2]]-10^(-10),izbFi=izbFi+0.1,
intKsi=Simplify[IntervalForKsi[P1p,d0p,d1p,izbFi]];
PP1=Expand[G1Interpolant[t,P1p,d0p,d1p,izbFi,intKsi[[2]]-10^(-15),1]];
RR1=Table[PP1[[i]]/PP1[[1]],{i,2,4}];
tabInt1=Append[tabInt1,RR1];
(*tabU=Append[tabU,Ukriv[RR1,t]];*)
]
nn=Length[tabInt1];
G1=Show[{Table[ParametricPlot3D[tabInt1[[i]],{t,0,1},PlotStyle->GrayLevel[0.3+0.4/(nn-1)(i-1)]],{i,1,nn}]},Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8],PlotRange->All]

Print["Table of interpolants for different φ and ξ chosen at a boundary: "]
tabPhi={0.5881,0.59,0.6,0.61,0.62,0.63,0.64,0.65,0.7,0.75,0.8,0.85,0.9,1}
tabInt1={};
tabU={};
For[j=1,j≤Length[tabPhi],j++,
izbFi=tabPhi[[j]];
intKsi=Simplify[IntervalForKsi[P1p,d0p,d1p,izbFi]];
PP1=Expand[G1Interpolant[t,P1p,d0p,d1p,izbFi,intKsi[[2]]-10^(-15),1]];
RR1=Table[PP1[[i]]/PP1[[1]],{i,2,4}];
tabInt1=Append[tabInt1,RR1];
(*tabU=Append[tabU,Ukriv[RR1,t]];*)
]
nn=Length[tabInt1];
G1=Show[{Table[ParametricPlot3D[tabInt1[[i]],{t,0,1},PlotStyle->GrayLevel[0.3+0.4/(nn-1)(i-1)]],{i,1,nn}],Graphics3D[{PointSize[0.02],Point[P0p],Point[P1p],Arrowheads[0.02],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}],Dashing[0.01],Line[{P0p,P1p}]}]},Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8],PlotRange->All]

Print["Table of interpolants for different φ and ξ chosen at a boundary: "]
tabPhi={1.2,1.4,1.6,1.8,2,2.5,3}
tabInt1={};
tabU={};
For[j=1,j≤Length[tabPhi],j++,
izbFi=tabPhi[[j]];
intKsi=Simplify[IntervalForKsi[P1p,d0p,d1p,izbFi]];
PP1=Expand[G1Interpolant[t,P1p,d0p,d1p,izbFi,intKsi[[2]]-10^(-15),1]];
RR1=Table[PP1[[i]]/PP1[[1]],{i,2,4}];
tabInt1=Append[tabInt1,RR1];
(*tabU=Append[tabU,Ukriv[RR1,t]];*)
]
nn=Length[tabInt1];
G1=Show[{Table[ParametricPlot3D[tabInt1[[i]],{t,0,1},PlotStyle->GrayLevel[0.5+0.3/(nn-1)(i-1)]],{i,1,nn}],Graphics3D[{PointSize[0.01],Point[P0p],Point[P1p],Arrowheads[0.01],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}],Dashing[0.005],Line[{P0p,P1p}]}]},Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8],PlotRange->All]

tabPhi=Join[{"φ:",0.5881,0.59,0.6,0.61,0.62,0.63,0.64,0.65,0.7,0.75,0.8,0.85,0.9,1},{1.2,1.4,1.6,1.8,2,2.5,3}];
tab={tabPhi,Join[{"L(φ):"},N[(dolz/.{φ->#})&/@Rest[tabPhi]]]};
Transpose[Transpose[tab]]//MatrixForm

(* Let us choose the length *)

pom=Solve[dolz-1.01(Norm[P1p-P0p])==0,φ]
izbFi=φ/.pom[[1]]

Print["Interpolant with a chosen length:"]
intKsi=Simplify[IntervalForKsi[P1p,d0p,d1p,izbFi]];
PP1=Expand[G1Interpolant[t,P1p,d0p,d1p,izbFi,intKsi[[2]]-10^(-10),1]];
RR1=Table[PP1[[i]]/PP1[[1]],{i,2,4}];

Show[{GData,ParametricPlot3D[RR1,{t,0,1}]}]

(G^1) interpolation - examples with C1<0

Example 1

(*  example  data  *)
precision1=50;
P0p={0,0,0};
P1p={4,2,4};
d0p={1,0,0};
d1p={6,-3,-2};
d1p=d1p/Norm[d1p];
Print["Data: P0 = ",P0p, ",  P1 = ", P1p , ",  d0 = ", d0p, ",  d1 = ", d1p]
G1=Graphics3D[{PointSize[0.02],Point[P0p],Point[P1p],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}],Dashing[0.01],Line[{P0p,P1p}]}];
Show[G1,Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Print["C1,C2,C3:"]
ConstC1C2C3[P1p,d1p]
Print["Norm of P1-P0:"]
N[Norm[P1p]]
Print["Length of the interpolant in dependence of phi: "]
dolz=Simplify[G1InterpolantLength[P1p,d0p,d1p,φ,ξ]]

Print["The admissible intervals for phi:"]
intPhi=Last[IntervalForPhi[P1p,d0p,d1p,1]]
N[%]

Print["Chosen phi: "]
izbFi=φ/.Flatten[Solve[dolz==10,φ]/.{C[1]->0}]
N[dolz/.{φ->izbFi},precision1]

Print["The admissible interval for xi based on the chosen phi:"]
intKsi=Simplify[IntervalForKsi[P1p,d0p,d1p,izbFi]]
N[%,7]

Print["Chosen ξ: "]
izbKsi=(intKsi[[1]]+intKsi[[2]])/2
N[izbKsi]

Print["First interpolant: "]
PP1=Expand[G1Interpolant[t,N[P1p,precision1],N[d0p,precision1],N[d1p,precision1],N[izbFi,precision1],N[izbKsi,precision1],1]];
RR1=Table[PP1[[i]]/PP1[[1]],{i,2,4}];
N[RR1]
G2a=ParametricPlot3D[RR1,{t,0,1},PlotStyle->GrayLevel[0.4]];
Show[{G1,G2a}]

Print["Second interpolant: "]
PP2=Expand[G1Interpolant[t,N[P1p,precision1],N[d0p,precision1],N[d1p,precision1],N[izbFi,precision1],N[izbKsi,precision1],2]];
RR2=Table[PP2[[i]]/PP2[[1]],{i,2,4}];
N[RR2]
G2b=ParametricPlot3D[RR2,{t,0,1},PlotStyle->GrayLevel[0.7]];
Show[{G1,G2b}]

Show[{G1,G2a,G2b},Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Example 2 - negative ξ

(*  example  data  *)
precision1=50;
P0p={0,0,0};
P1p={4,2,4};
d0p={1,0,0};
d1p={6,-3,-2};
d1p=d1p/Norm[d1p];
Print["Data: P0 = ",P0p, ",  P1 = ", P1p , ",  d0 = ", d0p, ",  d1 = ", d1p]
G1=Graphics3D[{PointSize[0.02],Point[P0p],Point[P1p],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}],Dashing[0.01],Line[{P0p,P1p}]}];
Show[G1,Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Print["C1,C2,C3:"]
ConstC1C2C3[P1p,d1p]
Print["Norm of P1-P0:"]
N[Norm[P1p]]
Print["Length of the interpolant in dependence of phi: "]
dolz=G1InterpolantLength[P1p,d0p,d1p,φ,ξ];
N[Collect[dolz,Cot[φ]]]

Print["The admissible intervals for phi:"]
intPhi=Last[IntervalForPhi[P1p,d0p,d1p,-1]]
N[%]

Print["Chosen phi: "]
izbFi=intPhi[[1]]-1/100;
izbFi=31/10

Print["The admissible interval for xi based on the chosen phi:"]
intKsi=Sort[-Simplify[IntervalForKsi[P1p,d0p,d1p,izbFi]]]
N[%,8]

Print["Chosen ξ: "]
izbKsi=(intKsi[[1]]+intKsi[[2]])/2
N[%]

Print["First interpolant: "]
PP1=Expand[G1Interpolant[t,N[P1p,precision1],N[d0p,precision1],N[d1p,precision1],N[izbFi,precision1],N[izbKsi,precision1],1]];
RR1=Table[PP1[[i]]/PP1[[1]],{i,2,4}];
N[RR1]
G2a=ParametricPlot3D[RR1,{t,0,1},PlotStyle->GrayLevel[0.4]];
Show[{G1,G2a},Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]
Print["Graph of psi*theta^4:"]
facTab=Simplify[G1InterpolantFactor[t,N[P1p,precision1],N[d0p,precision1],N[d1p,precision1],N[izbFi,precision1],N[izbKsi,precision1],1]];
Plot[facTab[[3]],{t,0,1}]

d0p-(D[RR1,t]/Norma[D[RR1,t]]/.{t->0})
d1p-(D[RR1,t]/Norma[D[RR1,t]]/.{t->1})

Print["Second interpolant: "]
PP2=Expand[G1Interpolant[t,N[P1p,precision1],N[d0p,precision1],N[d1p,precision1],N[izbFi,precision1],N[izbKsi,precision1],2]];
RR2=Table[PP2[[i]]/PP2[[1]],{i,2,4}];
N[RR2]
G2b=ParametricPlot3D[RR2,{t,0,1},PlotStyle->GrayLevel[0.7]];
Show[{G1,G2b},Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]
Print["Graph of psi*theta^4:"]
facTab=Simplify[G1InterpolantFactor[t,N[P1p,precision1],N[d0p,precision1],N[d1p,precision1],N[izbFi,precision1],N[izbKsi,precision1],2]];
Plot[facTab[[3]],{t,0,1}]

d0p-(D[RR2,t]/Norma[D[RR2,t]]/.{t->0})
d1p-(D[RR2,t]/Norma[D[RR2,t]]/.{t->1})

Show[{G1,G2a,G2b},Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Example 3 - graphs of interpolants for different ξ at a fixed φ

(*  example  data  *)
precision1=40;
P0p=N[{0,0,0},precision1];
P1p=N[{4,2,4},precision1];
d0p=N[{1,0,0},precision1];
d1p=N[{6,-3,-2},precision1];
d1p=d1p/Norm[d1p];

GData=Graphics3D[{PointSize[0.02],Point[P0p],Point[P1p],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}],Dashing[0.01],Line[{P0p,P1p}]}];
Show[GData,Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Print["Norm of P1-P0:"]
Norm[P1p]
Print["Length of the interpolant in dependence of phi: "]
dolz=Simplify[G1InterpolantLength[P1p,d0p,d1p,φ,ξ]]

Print["The admissible intervals for phi:"]
intPhi=Last[IntervalForPhi[P1p,d0p,d1p,1]]

Print["Chosen phi: "]
izbFi=φ/.Flatten[Solve[dolz==10,φ]/.{C[1]->0}]
Print["Length at a chosen phi: "]
N[dolz/.{φ->izbFi},precision1]

Print["The admissible interval for xi based on the chosen phi:"]
intKsi=Simplify[IntervalForKsi[P1p,d0p,d1p,izbFi]]

Print["Table of chosen values for xi: "]
tabKsi=Table[(1-i)intKsi[[1]]+i(intKsi[[2]]),{i,0,1,0.1}];
tabKsi=Table[intKsi[[1]]+i/9(intKsi[[2]]-intKsi[[1]]),{i,0,9}]
nn=Length[tabKsi];
tabKsi[[1]]=tabKsi[[1]]+10^(-7);
tabKsi[[nn]]=tabKsi[[nn]]-10^(-15);

Print["Table of first interpolants: "]
tabInt1={};
tabU={};
For[j=1,j≤Length[tabKsi],j++,
izb=tabKsi[[j]];
PP1=Expand[G1Interpolant[t,P1p,d0p,d1p,izbFi,izb,1]];
RR1=Table[PP1[[i]]/PP1[[1]],{i,2,4}];
tabInt1=Chop[Append[tabInt1,RR1]];
tabU=Append[tabU,Ukriv[N[RR1],t]];
]

G1=Show[{GData,Table[ParametricPlot3D[tabInt1[[i]],{t,0,1},PlotStyle->GrayLevel[0.3+(0.5/(1/nn-1))(i/nn-1)]],{i,1,nn}]},PlotRange->All,Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Print["Curvature: "]
G1U=Show[Table[Plot[tabU[[i]],{t,0,1},PlotStyle->GrayLevel[0.3+(0.5/(1/nn-1))(i/nn-1)]],{i,1,nn}],PlotRange->{All,All}]

Print["Table of second interpolants: "]
tabInt2={};
tabU2={};
For[j=1,j≤Length[tabKsi],j++,
izb=tabKsi[[j]];
PP1=Expand[G1Interpolant[t,P1p,d0p,d1p,izbFi,izb,2]];
RR1=Table[PP1[[i]]/PP1[[1]],{i,2,4}];
tabInt2=Chop[Append[tabInt2,RR1]];
tabU2=Append[tabU2,Ukriv[RR1,t]];
]
G2=Show[{GData,Table[ParametricPlot3D[tabInt2[[i]],{t,0,1},PlotStyle->{GrayLevel[0.3+(0.5/(1/nn-1))(i/nn-1)],Dashing[0.01]}],{i,1,nn}]},PlotRange->All,Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Print["Curvature: "]
G2U=Show[Table[Plot[tabU2[[i]],{t,0,1},PlotStyle->GrayLevel[0.3+(0.5/(1/nn-1))(i/nn-1)]],{i,1,nn}],PlotRange->{All,All}]

Show[{G1,G2},Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Print["Interpolants at the boundary values for xi:"]
Show[{ParametricPlot3D[tabInt1[[1]],{t,0,1},PlotStyle->GrayLevel[0.6]],ParametricPlot3D[Last[tabInt1],{t,0,1},PlotStyle->GrayLevel[0]],ParametricPlot3D[tabInt2[[1]],{t,0,1},PlotStyle->GrayLevel[0.8]],ParametricPlot3D[Last[tabInt2],{t,0,1},PlotStyle->GrayLevel[0]],Graphics3D[{PointSize[0.02],Dashing[0.005],Point[P0p],Point[P1p],Line[{P0p,P1p}],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}]}]},PlotRange->All,Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Example 4 - graphs of interpolants for different φ

(*  example  data  *)
precision1=50;
P0p=N[{0,0,0},precision1];
P1p=N[{4,2,4},precision1];
d0p=N[{1,0,0},precision1];
d1p=N[{6,-3,-2},precision1];
d1p=d1p/Norm[d1p];

GData=Graphics3D[{PointSize[0.02],Point[P0p],Point[P1p],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}],Dashing[0.01],Line[{P0p,P1p}]}];
Show[GData,Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8]]

Print["Norm of P1-P0:"]
Norm[P1p]
Print["Length of the interpolant in dependence of phi: "]
dolz=Simplify[G1InterpolantLength[P1p,d0p,d1p,φ,ξ]]

Print["The admissible intervals for phi:"]
intPhi=Last[IntervalForPhi[P1p,d0p,d1p,1]]

Print["Table of interpolants for different φ and ξ chosen at a boundary: "]
tabPhi=Table[i,{i,0.03,0.19,0.01}]
tabInt1={};
tabU={};
For[j=1,j≤Length[tabPhi],j++,
izbFi=tabPhi[[j]];
intKsi=Simplify[IntervalForKsi[P1p,d0p,d1p,izbFi]];
PP1=Expand[G1Interpolant[t,P1p,d0p,d1p,izbFi,intKsi[[2]]-10^(-15),1]];
RR1=Table[PP1[[i]]/PP1[[1]],{i,2,4}];
tabInt1=Append[tabInt1,RR1];
(*tabU=Append[tabU,Ukriv[RR1,t]];*)
]
nn=Length[tabInt1];
G1=Show[{Table[ParametricPlot3D[tabInt1[[i]],{t,0,1},PlotStyle->GrayLevel[0.3+0.4/(nn-1)(i-1)]],{i,1,nn}],Graphics3D[{PointSize[0.02],Point[P0p],Point[P1p],Arrowheads[0.02],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}],Dashing[0.01],Line[{P0p,P1p}]}]},Axes->True,ViewPoint->{0,-2,2},BoxStyle->GrayLevel[0.8],PlotRange->All]

(* Let us choose the length *)

pom=Solve[dolz-2(Norm[P1p-P0p])==0,φ]
izbFi=φ/.pom[[1]]

intKsi=Simplify[IntervalForKsi[P1p,d0p,d1p,izbFi]];
PP1=Expand[G1Interpolant[t,P1p,d0p,d1p,izbFi,intKsi[[2]]-10^(-10),1]];
RR1=Table[PP1[[i]]/PP1[[1]],{i,2,4}];

Show[{GData,ParametricPlot3D[RR1,{t,0,1}]}]

(G^1) interpolation - data in a general position

Example 1

(*  example  data  *)
precision1=50;
P0p={1,2,0};
P1p={4,5,5};
d0p={-2,2,1};
d0p=d0p/Norm[d0p];
d1p={3,-4,-5};
d1p={4,-16,-13};
d1p=d1p/Norm[d1p];
Print["Data: P0 = ",P0p, ",  P1 = ", P1p , ",  d0 = ", d0p, ",  d1 = ", d1p]
Gdata=Graphics3D[{PointSize[0.02],Point[P0p],Point[P1p],Arrowheads[0.025],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}],Dashing[0.01],Line[{P0p,P1p}]}];
Show[Gdata,Axes->True,ViewPoint->Bottom,BoxStyle->GrayLevel[0.8]]

Print["Transformed data:"]
QQR=Simplify[SolveEQStar[d0p,0]]; (* rotation *)
QQRi=KvaternionConj[QQR]/Norm[QQR];
P0t={0,0,0}
P1t=Rotacija[QQRi][[3]].(P1p-P0p)
d0t=Rotacija[QQRi][[3]].d0p
d1t=Rotacija[QQRi][[3]].d1p

Print["C1,C2,C3:"]
ConstC1C2C3[P1t,d1t]
Print["Norm of P1-P0:"]
N[Norm[P1t]]
Print["Length of the interpolant in dependence of phi: "]
dolz=G1InterpolantLength[P1t,d0t,d1t,φ,ξ];
N[Collect[dolz,Cot[φ]]]

Print["The admissible intervals for phi:"]
intPhi=Simplify[Last[IntervalForPhi[P1t,d0t,d1t,1]]]
N[%,7]

Print["Chosen length -> phi: "]
izbFi=φ/.N[Flatten[Solve[Simplify[dolz]==13/10Norma[P1p-P0p],φ]]/.{C[1]->1},precision1]

Print["The admissible interval for xi based on the chosen phi:"]
intKsi=Simplify[IntervalForKsi[P1t,d0t,d1t,izbFi]]

Print["Chosen ξ: "]
izbKsi=(intKsi[[1]]+intKsi[[2]])/2

Print["First interpolant: "]
PP1=Expand[G1Interpolant[t,N[P1t,precision1],N[d0t,precision1],N[d1t,precision1],N[izbFi,precision1],N[izbKsi,precision1],1]];
RR1=Table[PP1[[i]]/PP1[[1]],{i,2,4}];
Print["Transformed interpolant:"]
RR1=Chop[Simplify[Rotacija[QQRi][[3]].RR1+P0p]];
G1=ParametricPlot3D[RR1,{t,0,1},PlotStyle->GrayLevel[0.4]];
Show[{Gdata,G1},PlotRange->All,Axes->True,ViewPoint->Bottom,BoxStyle->GrayLevel[0.8]]

Print["Second interpolant: "]
PP2=Expand[G1Interpolant[t,N[P1t,precision1],N[d0t,precision1],N[d1t,precision1],N[izbFi,precision1],N[izbKsi,precision1],2]];
RR2=Table[PP2[[i]]/PP2[[1]],{i,2,4}];
Print["Transformed interpolant:"]
RR2=Chop[Simplify[Rotacija[QQRi][[3]].RR2+P0p]];
G2=ParametricPlot3D[RR2,{t,0,1},PlotStyle->GrayLevel[0.7]];
Show[{Gdata,G2},PlotRange->All,Axes->True,ViewPoint->Bottom,BoxStyle->GrayLevel[0.8]]

Show[{Gdata,G1,G2},PlotRange->All,Axes->True,ViewPoint->Bottom,BoxStyle->GrayLevel[0.8]]

Example 2

(*  example  data  *)
precision1=50;
P0p={1,1,-3};
P1p={3,-4,2};
d0p={-2,2,1};
d0p=d0p/Norm[d0p];
d1p={3,5,-4};
d1p={10,-1,-8};
d1p=d1p/Norm[d1p];
Print["Data: P0 = ",P0p, ",  P1 = ", P1p , ",  d0 = ", d0p, ",  d1 = ", d1p]
Gdata=Graphics3D[{PointSize[0.02],Point[P0p],Point[P1p],Arrowheads[0.025],Arrow[{P0p,P0p+d0p}],Arrow[{P1p,P1p+d1p}],Dashing[0.01],Line[{P0p,P1p}]}];
Show[Gdata,Axes->True,ViewPoint->Bottom,BoxStyle->GrayLevel[0.8]]

Print["Transformed data:"]
QQR=Simplify[SolveEQStar[d0p,0]]; (* rotation *)
QQRi=KvaternionConj[QQR]/Norm[QQR];
P0t={0,0,0}
P1t=Rotacija[QQRi][[3]].(P1p-P0p)
d0t=Rotacija[QQRi][[3]].d0p
d1t=Rotacija[QQRi][[3]].d1p

Print["C1,C2,C3:"]
ConstC1C2C3[P1t,d1t]
Print["Norm of P1-P0:"]
N[Norm[P1t]]
Print["Length of the interpolant in dependence of phi: "]
dolz=G1InterpolantLength[P1t,d0t,d1t,φ,ξ];
N[Collect[dolz,Cot[φ]]]

Print["The admissible intervals for phi:"]
intPhi=Simplify[Last[IntervalForPhi[P1t,d0t,d1t,1]]]
N[%,7]

Print["Chosen length -> phi: "]
izbFi=φ/.N[Flatten[Solve[Simplify[dolz]==13/10Norma[P1p-P0p],φ]]/.{C[1]->1},precision1]

Print["The admissible interval for xi based on the chosen phi:"]
intKsi=Simplify[IntervalForKsi[P1t,d0t,d1t,izbFi]]

Print["Chosen ξ: "]
izbKsi=(intKsi[[1]]+intKsi[[2]])/2

Print["First interpolant: "]
PP1=Expand[G1Interpolant[t,N[P1t,precision1],N[d0t,precision1],N[d1t,precision1],N[izbFi,precision1],N[izbKsi,precision1],1]];
RR1=Table[PP1[[i]]/PP1[[1]],{i,2,4}];
Print["Transformed interpolant:"]
RR1=Chop[Simplify[Rotacija[QQRi][[3]].RR1+P0p]];
G1=ParametricPlot3D[RR1,{t,0,1},PlotStyle->GrayLevel[0.4]];
Show[{Gdata,G1},PlotRange->All,Axes->True,ViewPoint->Bottom,BoxStyle->GrayLevel[0.8]]

Print["Second interpolant: "]
PP2=Expand[G1Interpolant[t,N[P1t,precision1],N[d0t,precision1],N[d1t,precision1],N[izbFi,precision1],N[izbKsi,precision1],2]];
RR2=Table[PP2[[i]]/PP2[[1]],{i,2,4}];
Print["Transformed interpolant:"]
RR2=Chop[Simplify[Rotacija[QQRi][[3]].RR2+P0p]];
G2=ParametricPlot3D[RR2,{t,0,1},PlotStyle->GrayLevel[0.7]];
Show[{Gdata,G2},PlotRange->All,Axes->True,ViewPoint->Bottom,BoxStyle->GrayLevel[0.8]]

Show[{Gdata,G1,G2},PlotRange->All,Axes->True,ViewPoint->Bottom,BoxStyle->GrayLevel[0.8]]

Spline examples

Example - spline 1

Print["Chosen curve: "]
curve={Cos[t]Log[4+t^2],Sin[t]Log[4+t^2],(Sqrt[4+(t+1)^2])}
curveDer=D[curve,t];
ParametricPlot3D[curve,{t,0,20},PlotStyle->GrayLevel[0],Axes->True,ViewPoint->Right,BoxStyle->GrayLevel[0.8]]

Print["The data are sampled from the chosen curve."]
izpis=0;
a=0;
b=6;
mm=6;
par=Table[a+i(b-a)/mm,{i,0,mm}];
Print["Parameters:",par];
data={};
For[i=1,i≤Length[par]-1,i++,
data=Append[data,N[{curve/.{t->par[[i]]},curve/.{t->par[[i+1]]},curveDer/Norma[curveDer]/.{t->par[[i]]},curveDer/Norma[curveDer]/.{t->par[[i+1]]}},50]];
];
GData= Show[Graphics3D[{PointSize[0.02],GrayLevel[0],Table[Point[{data[[i,1]],data[[i,2]]}],{i,1,Length[data]}],Arrowheads[0.02],Table[Arrow[{data[[i,1]],data[[i,1]]+data[[i,3]]}],{i,1,Length[data]}],Table[Arrow[{data[[i,2]],data[[i,2]]+data[[i,4]]}],{i,1,Length[data]}],Dashing[0.01],Table[Line[{data[[i,1]],data[[i,2]]}],{i,1,Length[data]}]}],PlotRange->All,Axes->True,ViewPoint->Right,BoxStyle->GrayLevel[0.8]]

Print["We determine the length of the interpolant on each segment by the formula (1+fac)||P_{i+1}-P_i|| for some chosen factor 'fac'."]

fac=1.01;
tabScal=Table[fac,{i,1,Length[par]-1}];
Print["Table of scalar factors: ", tabScal];
(*tabScal=Table[1.03,{i,1,Length[par]-1}]*)
Print["Interpolants at a computed scalar factor: "];
izbK=1; (* izbK=0->ξ=(Subscript[ξ, **]), izbK=1/2->ξ=((Subscript[ξ, **])+(SubStar[ξ]))/2, izbK=1->ξ=(SubStar[ξ]) *)
izbS=1; (* izbS=1->(SubPlus[r]), izbS=1->(SubMinus[r]) *)
tabRR={};
For[i=1,i≤Length[data],i++,
tabRR=Append[tabRR,Chop[G1InterpolantGen[t,data[[i,1]],data[[i,2]],data[[i,3]],data[[i,4]],tabScal[[i]],izbK,izbS,izpis]]];
];
G1=Show[Table[ParametricPlot3D[tabRR[[i]],{t,0,1},PlotStyle->{GrayLevel[0.3],Thickness[0.003]}],{i,1,Length[tabRR]}],PlotRange->All,Axes->True,ViewPoint->Right,BoxStyle->GrayLevel[0.8]]

fac=1.03;
tabScal=Table[fac,{i,1,Length[par]-1}];
Print["Table of scalar factors: ", tabScal];
(*tabScal=Table[1.03,{i,1,Length[par]-1}]*)
Print["Interpolants at a computed scalar factor: "];
izbK=1; (* izbK=0->ξ=(Subscript[ξ, **]), izbK=1/2->ξ=((Subscript[ξ, **])+(SubStar[ξ]))/2, izbK=1->ξ=(SubStar[ξ]) *)
izbS=1; (* izbS=1->(SubPlus[r]), izbS=1->(SubMinus[r]) *)
tabRR={};
For[i=1,i≤Length[data],i++,
tabRR=Append[tabRR,Chop[G1InterpolantGen[t,data[[i,1]],data[[i,2]],data[[i,3]],data[[i,4]],tabScal[[i]],izbK,izbS,izpis]]];
];
G2=Show[Table[ParametricPlot3D[tabRR[[i]],{t,0,1},PlotStyle->{GrayLevel[0.5],Thickness[0.003]}],{i,1,Length[tabRR]}],PlotRange->All,Axes->True,ViewPoint->Right,BoxStyle->GrayLevel[0.8]]

fac=1.05;
tabScal=Table[fac,{i,1,Length[par]-1}];
Print["Table of scalar factors: ", tabScal];
(*tabScal=Table[1.03,{i,1,Length[par]-1}]*)
Print["Interpolants at a computed scalar factor: "];
izbK=1; (* izbK=0->ξ=(Subscript[ξ, **]), izbK=1/2->ξ=((Subscript[ξ, **])+(SubStar[ξ]))/2, izbK=1->ξ=(SubStar[ξ]) *)
izbS=1; (* izbS=1->(SubPlus[r]), izbS=1->(SubMinus[r]) *)
tabRR={};
For[i=1,i≤Length[data],i++,
tabRR=Append[tabRR,Chop[G1InterpolantGen[t,data[[i,1]],data[[i,2]],data[[i,3]],data[[i,4]],tabScal[[i]],izbK,izbS,izpis]]];
];
G3=Show[Table[ParametricPlot3D[tabRR[[i]],{t,0,1},PlotStyle->{GrayLevel[0.7],Thickness[0.003]}],{i,1,Length[tabRR]}],PlotRange->All,Axes->True,ViewPoint->Right,BoxStyle->GrayLevel[0.8]]

fac=1.5;
tabScal=Table[fac,{i,1,Length[par]-1}];
Print["Table of scalar factors: ", tabScal];
(*tabScal=Table[1.03,{i,1,Length[par]-1}]*)
Print["Interpolants at a computed scalar factor: "];
izbK=1; (* izbK=0->ξ=(Subscript[ξ, **]), izbK=1/2->ξ=((Subscript[ξ, **])+(SubStar[ξ]))/2, izbK=1->ξ=(SubStar[ξ]) *)
izbS=1; (* izbS=1->(SubPlus[r]), izbS=1->(SubMinus[r]) *)
tabRR={};
For[i=1,i≤Length[data],i++,
tabRR=Append[tabRR,Chop[G1InterpolantGen[t,data[[i,1]],data[[i,2]],data[[i,3]],data[[i,4]],tabScal[[i]],izbK,izbS,izpis]]];
];
G4=Show[Table[ParametricPlot3D[tabRR[[i]],{t,0,1},PlotStyle->{GrayLevel[0.8],Thickness[0.003]}],{i,1,Length[tabRR]}],PlotRange->All,Axes->True,ViewPoint->Right,BoxStyle->GrayLevel[0.8]]

Show[{G1,G2,G3,G4,Graphics3D[{PointSize[0.02],GrayLevel[0],Table[Point[{data[[i,1]],data[[i,2]]}],{i,1,Length[data]}],Arrowheads[0.015],Table[Arrow[{data[[i,1]],data[[i,1]]+data[[i,3]]}],{i,1,Length[data]}],Table[Arrow[{data[[i,2]],data[[i,2]]+data[[i,4]]}],{i,1,Length[data]}]}]}]

(* --------------------------------------------------- *) 
(* Splines that interpolate also the length of the chosen curve for different number of segments.  *)

Print["The spline that interpolates also the length of the chosen curve. "](* Lengths computed from the curve *)
izpis=0;
a=0;
b=20;
mm=4;
par=Table[a+i(b-a)/mm,{i,0,mm}];
Print["Parameters:",par];
data={};
For[i=1,i≤Length[par]-1,i++,
data=Append[data,N[{curve/.{t->par[[i]]},curve/.{t->par[[i+1]]},curveDer/Norma[curveDer]/.{t->par[[i]]},curveDer/Norma[curveDer]/.{t->par[[i+1]]}},50]];
];GData= Show[{Graphics3D[{PointSize[0.025],GrayLevel[0],Table[Point[{data[[i,1]],data[[i,2]]}],{i,1,Length[data]}]}],ParametricPlot3D[curve,{t,par[[1]],Last[par]},PlotStyle->{GrayLevel[0]}]},PlotRange->All,Axes->True,ViewPoint->Right,BoxStyle->GrayLevel[0.8]];

tabScal=Table[NIntegrate[Norma[curveDer],{t,par[[i]],par[[i+1]]}]/Norma[data[[i,2]]-data[[i,1]]],{i,1,Length[par]-1}];
Print["Table of scalar factors: ", tabScal];
Print["Interpolants at a computed scalar factor: "];
izbK=1; (* izbK=0->ξ=(Subscript[ξ, **]), izbK=1/2->ξ=((Subscript[ξ, **])+(SubStar[ξ]))/2, izbK=1->ξ=(SubStar[ξ]) *)
izbS=1; (* izbS=1->(SubPlus[r]), izbS=1->(SubMinus[r]) *)
tabRR={};
For[i=1,i≤Length[data],i++,
tabRR=Append[tabRR,Chop[G1InterpolantGen[t,data[[i,1]],data[[i,2]],data[[i,3]],data[[i,4]],tabScal[[i]],izbK,izbS,izpis]]];
];
G1=Show[Table[ParametricPlot3D[tabRR[[i]],{t,0,1},PlotStyle->{GrayLevel[0.5]}],{i,1,Length[tabRR]}],PlotRange->All,Axes->True,ViewPoint->Right,BoxStyle->GrayLevel[0.8]];
Show[{G1,GData}]
TabNap={};
TabNapF={};
For[i=1,i≤Length[tabRR],i++,
vd0=Norma[D[tabRR[[i]],t]/.{t->0}]^2/((curveDer/.{t->par[[i]]}).(D[tabRR[[i]],t]/.{t->0}));
vd1=Norma[D[tabRR[[i]],t]/.{t->1}]^2/((curveDer/.{t->par[[i+1]]}).(D[tabRR[[i]],t]/.{t->1}));
rep=InterpolatingPolynomial[{{0,{par[[i]],vd0}},{1,{par[[i+1]],vd1}}},x];
napF=Norma[(curve/.{t->rep})-(tabRR[[i]]/.{t->x})];
TabNapF=Append[TabNapF,napF];
TabNap=Append[TabNap,NMaximize[{napF,0<=x≤1},x]];
];
Print["Error:"]
Max[First[#]&/@TabNap]

izpis=0;
a=0;
b=20;
mm=5;
par=Table[a+i(b-a)/mm,{i,0,mm}];
Print["Parameters:",par];
data={};
For[i=1,i≤Length[par]-1,i++,
data=Append[data,N[{curve/.{t->par[[i]]},curve/.{t->par[[i+1]]},curveDer/Norma[curveDer]/.{t->par[[i]]},curveDer/Norma[curveDer]/.{t->par[[i+1]]}},50]];
];GData= Show[{Graphics3D[{PointSize[0.025],GrayLevel[0],Table[Point[{data[[i,1]],data[[i,2]]}],{i,1,Length[data]}]}],ParametricPlot3D[curve,{t,par[[1]],Last[par]},PlotStyle->{GrayLevel[0]}]},PlotRange->All,Axes->True,ViewPoint->Right,BoxStyle->GrayLevel[0.8]];
tabScal=Table[NIntegrate[Norma[curveDer],{t,par[[i]],par[[i+1]]}]/Norma[data[[i,2]]-data[[i,1]]],{i,1,Length[par]-1}];
Print["Table of scalar factors: ", tabScal];
Print["Interpolants at a computed scalar factor: "];
izbK=1; (* izbK=0->ξ=(Subscript[ξ, **]), izbK=1/2->ξ=((Subscript[ξ, **])+(SubStar[ξ]))/2, izbK=1->ξ=(SubStar[ξ]) *)
izbS=1; (* izbS=1->(SubPlus[r]), izbS=1->(SubMinus[r]) *)
tabRR={};
For[i=1,i≤Length[data],i++,
tabRR=Append[tabRR,Chop[G1InterpolantGen[t,data[[i,1]],data[[i,2]],data[[i,3]],data[[i,4]],tabScal[[i]],izbK,izbS,izpis]]];
];
G1=Show[Table[ParametricPlot3D[tabRR[[i]],{t,0,1},PlotStyle->{GrayLevel[0.5]}],{i,1,Length[tabRR]}],PlotRange->All,Axes->True,ViewPoint->Right,BoxStyle->GrayLevel[0.8]];
Show[{G1,GData}]
TabNap={};
TabNapF={};
For[i=1,i≤Length[tabRR],i++,
vd0=Norma[D[tabRR[[i]],t]/.{t->0}]^2/((curveDer/.{t->par[[i]]}).(D[tabRR[[i]],t]/.{t->0}));
vd1=Norma[D[tabRR[[i]],t]/.{t->1}]^2/((curveDer/.{t->par[[i+1]]}).(D[tabRR[[i]],t]/.{t->1}));
rep=InterpolatingPolynomial[{{0,{par[[i]],vd0}},{1,{par[[i+1]],vd1}}},x];
napF=Norma[(curve/.{t->rep})-(tabRR[[i]]/.{t->x})];
TabNapF=Append[TabNapF,napF];
TabNap=Append[TabNap,NMaximize[{napF,0<=x≤1},x]];
];
Print["Error:"]
Max[First[#]&/@TabNap]

izpis=0;
a=0;
b=20;
mm=6;
par=Table[a+i(b-a)/mm,{i,0,mm}];
Print["Parameters:",par];
data={};
For[i=1,i≤Length[par]-1,i++,
data=Append[data,N[{curve/.{t->par[[i]]},curve/.{t->par[[i+1]]},curveDer/Norma[curveDer]/.{t->par[[i]]},curveDer/Norma[curveDer]/.{t->par[[i+1]]}},50]];
];GData= Show[{Graphics3D[{PointSize[0.025],GrayLevel[0],Table[Point[{data[[i,1]],data[[i,2]]}],{i,1,Length[data]}]}],ParametricPlot3D[curve,{t,par[[1]],Last[par]},PlotStyle->{GrayLevel[0]}]},PlotRange->All,Axes->True,ViewPoint->Right,BoxStyle->GrayLevel[0.8]];
tabScal=Table[NIntegrate[Norma[curveDer],{t,par[[i]],par[[i+1]]}]/Norma[data[[i,2]]-data[[i,1]]],{i,1,Length[par]-1}];
Print["Table of scalar factors: ", tabScal];
Print["Interpolants at a computed scalar factor: "];
izbK=1; (* izbK=0->ξ=(Subscript[ξ, **]), izbK=1/2->ξ=((Subscript[ξ, **])+(SubStar[ξ]))/2, izbK=1->ξ=(SubStar[ξ]) *)
izbS=1; (* izbS=1->(SubPlus[r]), izbS=1->(SubMinus[r]) *)
tabRR={};
For[i=1,i≤Length[data],i++,
tabRR=Append[tabRR,Chop[G1InterpolantGen[t,data[[i,1]],data[[i,2]],data[[i,3]],data[[i,4]],tabScal[[i]],izbK,izbS,izpis]]];
];
G1=Show[Table[ParametricPlot3D[tabRR[[i]],{t,0,1},PlotStyle->{GrayLevel[0.5]}],{i,1,Length[tabRR]}],PlotRange->All,Axes->True,ViewPoint->Right,BoxStyle->GrayLevel[0.8]];
Show[{G1,GData}]
TabNap={};
TabNapF={};
For[i=1,i≤Length[tabRR],i++,
vd0=Norma[D[tabRR[[i]],t]/.{t->0}]^2/((curveDer/.{t->par[[i]]}).(D[tabRR[[i]],t]/.{t->0}));
vd1=Norma[D[tabRR[[i]],t]/.{t->1}]^2/((curveDer/.{t->par[[i+1]]}).(D[tabRR[[i]],t]/.{t->1}));
rep=InterpolatingPolynomial[{{0,{par[[i]],vd0}},{1,{par[[i+1]],vd1}}},x];
napF=Norma[(curve/.{t->rep})-(tabRR[[i]]/.{t->x})];
TabNapF=Append[TabNapF,napF];
TabNap=Append[TabNap,NMaximize[{napF,0<=x≤1},x]];
];
Print["Error:"]
Max[First[#]&/@TabNap]

tabela={};
tabH={};
tabM={};
For[jj=20,jj≤50,jj=jj+5,
izpis=0;
a=0;
b=20;
mm=jj;
tabM=Append[tabM,mm];
tabH=Append[tabH,(b-a)/mm];
par=Table[a+i(b-a)/mm,{i,0,mm}];
Print["Parameters:",par];
data={};
For[i=1,i≤Length[par]-1,i++,
data=Append[data,N[{curve/.{t->par[[i]]},curve/.{t->par[[i+1]]},curveDer/Norma[curveDer]/.{t->par[[i]]},curveDer/Norma[curveDer]/.{t->par[[i+1]]}},50]];
];GData= Show[{Graphics3D[{PointSize[0.02],GrayLevel[0],Table[Point[{data[[i,1]],data[[i,2]]}],{i,1,Length[data]}]}],ParametricPlot3D[curve,{t,par[[1]],Last[par]},PlotStyle->{GrayLevel[0]}]},PlotRange->All,Axes->True,ViewPoint->Right,BoxStyle->GrayLevel[0.8]];
tabScal=Table[NIntegrate[Norma[curveDer],{t,par[[i]],par[[i+1]]}]/Norma[data[[i,2]]-data[[i,1]]],{i,1,Length[par]-1}];
Print["Table of scalar factors: ", tabScal];
Print["Interpolants at a computed scalar factor: "];
izbK=1; (* izbK=0->ξ=(Subscript[ξ, **]), izbK=1/2->ξ=((Subscript[ξ, **])+(SubStar[ξ]))/2, izbK=1->ξ=(SubStar[ξ]) *)
izbS=1; (* izbS=1->(SubPlus[r]), izbS=1->(SubMinus[r]) *)
tabRR={};
For[i=1,i≤Length[data],i++,
tabRR=Append[tabRR,Chop[G1InterpolantGen[t,data[[i,1]],data[[i,2]],data[[i,3]],data[[i,4]],tabScal[[i]],izbK,izbS,izpis]]];
];
(*G1=Show[Table[ParametricPlot3D[tabRR[[i]],{t,0,1},PlotStyle->{GrayLevel[0.5],Thickness[0.002]}],{i,1,Length[tabRR]}],PlotRange->All,Axes->True,ViewPoint->Right,BoxStyle->GrayLevel[0.8]];
Print[Show[{G1,GData}]];*)
TabNap={};
TabNapF={};
For[i=1,i≤Length[tabRR],i++,
vd0=Norma[D[tabRR[[i]],t]/.{t->0}]^2/((curveDer/.{t->par[[i]]}).(D[tabRR[[i]],t]/.{t->0}));
vd1=Norma[D[tabRR[[i]],t]/.{t->1}]^2/((curveDer/.{t->par[[i+1]]}).(D[tabRR[[i]],t]/.{t->1}));
rep=InterpolatingPolynomial[{{0,{par[[i]],vd0}},{1,{par[[i+1]],vd1}}},x];
napF=Norma[(curve/.{t->rep})-(tabRR[[i]]/.{t->x})];
TabNapF=Append[TabNapF,napF];
TabNap=Append[TabNap,NMaximize[{napF,0<=x≤1},x]];
];
tabela=Append[tabela,Max[First[#]&/@TabNap]];
Print["Error:", Max[First[#]&/@TabNap]];
]

Print["Approximation order:"]
RedAproksimacije[tabH,tabela]//MatrixForm

RedAproksimacijeZaZlepek[tabM,tabela]//MatrixForm

Example - spline 2

Print["Viviani's curve "]
curve=4{1+Cos[t],Sin[t],1Sin[1/2t]}
curveDer=D[curve,t];
a=-2Pi;
b=2Pi;
GData= Show[{ParametricPlot3D[curve,{t,a,b},PlotStyle->{GrayLevel[0]}]},PlotRange->All,Axes->True,ViewPoint->{2,0,2},BoxStyle->GrayLevel[0.8]]

izpis=0;
mm=4;
par=Table[a+(b-a)/mm i,{i,0,mm}];
Print["Parameters:",par];
data={};
For[i=1,i≤Length[par]-1,i++,
data=Append[data,N[{curve/.{t->par[[i]]},curve/.{t->par[[i+1]]},curveDer/Norma[curveDer]/.{t->par[[i]]},curveDer/Norma[curveDer]/.{t->par[[i+1]]}},50]];
];
tabScal=Table[NIntegrate[Norma[curveDer],{t,par[[i]],par[[i+1]]}]/Norma[data[[i,2]]-data[[i,1]]],{i,1,Length[par]-1}];
Print["Table of scalar factors: ", tabScal];
(*tabScal=Table[1.3,{i,1,Length[par]-1}]*)
Print["Interpolants at a computed scalar factor: "];
izbK=1; (* izbK=0->ξ=(Subscript[ξ, **]), izbK=1/2->ξ=((Subscript[ξ, **])+(SubStar[ξ]))/2, izbK=1->ξ=(SubStar[ξ]) *)
izbS=1; (* izbS=1->(SubPlus[r]), izbS=1->(SubMinus[r]) *)
tabRR={};
For[i=1,i≤Length[data],i++,
tabRR=Append[tabRR,Chop[G1InterpolantGen[t,data[[i,1]],data[[i,2]],data[[i,3]],data[[i,4]],tabScal[[i]],izbK,izbS,izpis]]];
];
GData= Show[{Graphics3D[{Text[Style["(Subscript[P, 0])=(Subscript[P, 2])=(Subscript[P, 4])",Bold,Italic,13],data[[1,1]]+{0.7,0,0}],Text[Style["(Subscript[P, 1])",Bold,Italic,13],data[[2,1]]+{0.7,0,0}],Text[Style["(Subscript[P, 3])",Bold,Italic,13],data[[4,1]]+{0.7,0,0}],PointSize[0.025],GrayLevel[0],Table[Point[{data[[i,1]],data[[i,2]]}],{i,1,Length[data]}]}],ParametricPlot3D[curve,{t,par[[1]],Last[par]},PlotStyle->{GrayLevel[0]}]},PlotRange->All,Axes->True,ViewPoint->Right,BoxStyle->GrayLevel[0.8]];
Show[{GData,Table[ParametricPlot3D[tabRR[[i]],{t,0,1},PlotStyle->{GrayLevel[0.6]}],{i,1,Length[tabRR]}]},PlotRange->All,Axes->True,ViewPoint->{2,0,2},BoxStyle->GrayLevel[0.8]]
TabNap={};
TabNapF={};
For[i=1,i≤Length[tabRR],i++,
vd0=Norma[D[tabRR[[i]],t]/.{t->0}]^2/((curveDer/.{t->par[[i]]}).(D[tabRR[[i]],t]/.{t->0}));
vd1=Norma[D[tabRR[[i]],t]/.{t->1}]^2/((curveDer/.{t->par[[i+1]]}).(D[tabRR[[i]],t]/.{t->1}));
rep=InterpolatingPolynomial[{{0,{par[[i]],vd0}},{1,{par[[i+1]],vd1}}},x];
napF=Norma[(curve/.{t->rep})-(tabRR[[i]]/.{t->x})];
TabNapF=Append[TabNapF,napF];
TabNap=Append[TabNap,NMaximize[{napF,0<=x≤1},x]];
];
Print["Error:"]
Max[First[#]&/@TabNap]

Example - spline 3

(* curve *)
curve={(Sqrt[1+Exp[t]]),(2+t^2),t Sin[3t]}
curveDer=D[curve,t];
a=-5;
b=5;
Show[{ParametricPlot3D[curve,{t,a,b},PlotStyle->{GrayLevel[0],Thickness[0.002]}]},PlotRange->All,Axes->True,ViewPoint->{-3,-2,4},BoxStyle->GrayLevel[0.8]]

izpis=0;
mm=16;
par=Table[a+(b-a)/mm i,{i,0,mm}];
Print["Parameters:",par];
data={};
For[i=1,i≤Length[par]-1,i++,
data=Append[data,N[{curve/.{t->par[[i]]},curve/.{t->par[[i+1]]},curveDer/Norma[curveDer]/.{t->par[[i]]},curveDer/Norma[curveDer]/.{t->par[[i+1]]}},50]];
];
tabScal=Table[NIntegrate[Norma[curveDer],{t,par[[i]],par[[i+1]]}]/Norma[data[[i,2]]-data[[i,1]]],{i,1,Length[par]-1}];
Print["Table of scalar factors: ", tabScal];
(*tabScal=Table[1.3,{i,1,Length[par]-1}]*)
Print["Interpolants at a computed scalar factor: "];
izbK=1; (* izbK=0->ξ=(Subscript[ξ, **]), izbK=1/2->ξ=((Subscript[ξ, **])+(SubStar[ξ]))/2, izbK=1->ξ=(SubStar[ξ]) *)
izbS=1; (* izbS=1->(SubPlus[r]), izbS=1->(SubMinus[r]) *)
tabRR={};
For[i=1,i≤Length[data],i++,
tabRR=Append[tabRR,Chop[G1InterpolantGen[t,data[[i,1]],data[[i,2]],data[[i,3]],data[[i,4]],tabScal[[i]],izbK,izbS,izpis]]];
];
GData= Show[{Graphics3D[{PointSize[0.015],GrayLevel[0],Table[Point[{data[[i,1]],data[[i,2]]}],{i,1,Length[data]}]}],ParametricPlot3D[curve,{t,par[[1]],Last[par]},PlotStyle->{GrayLevel[0],Thickness[0.002]}]},PlotRange->All,Axes->True,ViewPoint->Right,BoxStyle->GrayLevel[0.8]];
Show[{GData,Table[ParametricPlot3D[tabRR[[i]],{t,0,1},PlotStyle->{GrayLevel[0.6],Thickness[0.002]}],{i,1,Length[tabRR]}]},PlotRange->All,Axes->True,ViewPoint->{-3,-2,4},BoxStyle->GrayLevel[0.8]]
TabNap={};
TabNapF={};
For[i=1,i≤Length[tabRR],i++,
vd0=Norma[D[tabRR[[i]],t]/.{t->0}]^2/((curveDer/.{t->par[[i]]}).(D[tabRR[[i]],t]/.{t->0}));
vd1=Norma[D[tabRR[[i]],t]/.{t->1}]^2/((curveDer/.{t->par[[i+1]]}).(D[tabRR[[i]],t]/.{t->1}));
rep=InterpolatingPolynomial[{{0,{par[[i]],vd0}},{1,{par[[i+1]],vd1}}},x];
napF=Norma[(curve/.{t->rep})-(tabRR[[i]]/.{t->x})];
TabNapF=Append[TabNapF,napF];
TabNap=Append[TabNap,NMaximize[{napF,0<=x≤1},x]];
];
Print["Error:"]
Max[First[#]&/@TabNap]