gaby

function[ybp]=JLCSplines(x,y,ypL,ypR)
// Function to calculate clamped cubic splines
//  Detailed explanation goes here
//   x independient variable
//   y dependient variable
//   ypL = dy/dx in the left side i.e. x = xmin
//   ypR = dy/dx in the right side, i.e. x = xmax
//   ybp = d ( dy/dx )/dx = Second derivative

[m,nL]=size(x);
    if m>nL then
        ngg=m
    else
        ngg=nL
    end

N=ngg-1;

for i=1:N
    h(i)=x(i+1)-x(i);
end

a(N+1)=h(N)/6.0;
b(N+1)=h(N)/2.0-h(N)/6.0;
c(N+1)=0.0;
b(1)=h(1)/2.0-h(1)/6.0;
c(1)=h(1)/6.0;
a(1)=0.0;

for i=1:N-1
    a(i+1)=h(i)/6.0;
    b(i+1)=(h(i)+h(i+1))/2.0-(h(i)+h(i+1))/6.0;
    c(i+1)=h(i+1)/6.0;
    r(i+1)=(y(i+2)-y(i+1))/h(i+1)-(y(i+1)-y(i))/h(i);
end

r(1)=(y(2)-y(1))/h(1)-ypL;
r(N+1)=ypR-(y(N+1)-y(N))/h(N);
aa=a';bb=b';cc=c';rr=r';
ybp=JLTridiag(aa,bb,cc,rr);

yp = (ypL+ypR)/2;

kappa=-ybp/(1+(yp)^2)^(3/2);

R = 1/kappa;

endfunction