rate_libary.m

Incorporating Prior Knowledge in Cubic Spline Approximation - Application to the Identification of R

%RATE_LIBARY

X = [];
Y = [];
b = [];


if libarycommand == 1,
  %Calculate b1,..b4 based on Ca Cb Cc vs. t
  %inputs: t(),c()
  %output: b()
  %work-variables: X(),Y(),n,i,dt,intca,intcb,intcc
  b = zeros(4,1);
  n = length(t)-1;
  % A
  j = 1;
  intca = 0; intcb = 0; intcc = 0;
  Y(1:n,1) = c(2:n+1,j)-c(1,j);
  for i=1:n,
    dt = t(i+1)-t(i);
    intca = intca + (c(i+1,1)+c(i,1))/2*dt;
    intcb = intcb + (c(i+1,2)+c(i,2))/2*dt;
    intcc = intcc + (c(i+1,3)+c(i,3))/2*dt;
    X(i,1) = -intca;
    X(i,2) = +intcb;
  end
  % C 
  j = 3;
  intca = 0; intcb = 0; intcc = 0;
  Y(n+1:2*n,1) = c(2:n+1,j)-c(1,j);
  for i=1:n,
    dt = t(i+1)-t(i);
    intca = intca + (c(i+1,1)+c(i,1))/2*dt;
    intcb = intcb + (c(i+1,2)+c(i,2))/2*dt;
    intcc = intcc + (c(i+1,3)+c(i,3))/2*dt;
    X(n+i,3) = +intcb;
    X(n+i,4) = -intcc;
  end
  % B 
  j = 2;
  intca = 0; intcb = 0; intcc = 0;
  Y(2*n+1:3*n,1) = c(2:n+1,j)-c(1,j);
  for i=1:n,
    dt = t(i+1)-t(i);
    intca = intca + (c(i+1,1)+c(i,1))/2*dt;
    intcb = intcb + (c(i+1,2)+c(i,2))/2*dt;
    intcc = intcc + (c(i+1,3)+c(i,3))/2*dt;
    X(2*n+i,1) = +intca;
    X(2*n+i,2) = -intcb;
    X(2*n+i,3) = -intcb;
    X(2*n+i,4) = +intcc;
  end
  %result
  b = X \ Y;
  return
end


if libarycommand == 2,
  %Calculate b1,..b4 based on Ca Cb vs. t
  %inputs: t(),c()
  %output: b()
  %work-variables: X(),Y(),n,i,dt,intca,intcb,intcc
  b = zeros(4,1);
  n = length(t)-1;
  % A
  j = 1;
  intca = 0; intcb = 0; intcc = 0;
  Y(1:n,1) = c(2:n+1,j)-c(1,j);
  for i=1:n,
    dt = t(i+1)-t(i);
    intca = intca + (c(i+1,1)+c(i,1))/2*dt;
    intcb = intcb + (c(i+1,2)+c(i,2))/2*dt;
    intcc = intcc + (c(i+1,3)+c(i,3))/2*dt;
    X(i,1) = -intca;
    X(i,2) = +intcb;
  end
  % C 
  j = 3;
  intca = 0; intcb = 0; intcc = 0;
  Y(n+1:2*n,1) = c(2:n+1,j)-c(1,j);
  for i=1:n,
    dt = t(i+1)-t(i);
    intca = intca + (c(i+1,1)+c(i,1))/2*dt;
    intcb = intcb + (c(i+1,2)+c(i,2))/2*dt;
    intcc = intcc + (c(i+1,3)+c(i,3))/2*dt;
    X(n+i,3) = +intcb;
    X(n+i,4) = -intcc;
  end
  %result
  b = X \ Y;
  return
end

if libarycommand == 3,
  %Calculate b1,..b4 based on Ca Cb Cc vs. t
  %inputs: t(),c(),ic()
  %output: b()
  %work-variables: X(),Y(),n,i,dt,intca,intcb,intcc
  b = zeros(4,1);
  n = length(t)-1;
  % A
  j = 1;
  Y(1:n,1) = c(2:n+1,j)-c(1,j);
  X(1:n,1) = -ic(2:n+1,1);
  X(1:n,2) = +ic(2:n+1,2);
  % C 
  j = 3;
  intca = 0; intcb = 0; intcc = 0;
  Y(n+1:2*n,1) = c(2:n+1,j)-c(1,j);
  X(n+1:2*n,3) = +ic(2:n+1,2);
  X(n+1:2*n,4) = -ic(2:n+1,3);
  % B 
  j = 2;
  Y(2*n+1:3*n,1) = c(2:n+1,j)-c(1,j);
  X(2*n+1:3*n,1) = +ic(2:n+1,1);
  X(2*n+1:3*n,2) = -ic(2:n+1,2);
  X(2*n+1:3*n,3) = -ic(2:n+1,2);
  X(2*n+1:3*n,4) = +ic(2:n+1,3);
  %result
  b = X \ Y;
  return
end


if libarycommand == 4,
  %Calculate b1,..b4 based on Ca Cb vs. t
  %inputs: t(),c()
  %output: b()
  %work-variables: X(),Y(),n,i,dt,intca,intcb,intcc
  b = zeros(4,1);
  n = length(t)-1;
  % A
  j = 1;
  Y(1:n,1) = c(2:n+1,j)-c(1,j);
  X(1:n,1) = -ic(2:n+1,1);
  X(1:n,2) = +ic(2:n+1,2);
  % C 
  j = 3;
  intca = 0; intcb = 0; intcc = 0;
  Y(n+1:2*n,1) = c(2:n+1,j)-c(1,j);
  X(n+1:2*n,3) = +ic(2:n+1,2);
  X(n+1:2*n,4) = -ic(2:n+1,3);
  %result
  b = X \ Y;
  return
end