% Class 8
% Coded by NFR on 2.7.2019
% Numerical differentiation
%
close all
clc
% Warm up activity on splines
v = [0 5 12 30 100 200];
t = [0 100 200 300 400 500];
v_adj = [0 v 0];
tq = linspace(0,500,1000);
vq = spline(t,v_adj,tq);
plot(t,v,'ro',tq,vq,'b:')
%
% Example: 21.1 from the BOOK
F1 = @(x) exp(x);
xval = 2;
h = 0.1;
%O(h^2) first derivative (centered)
FD_h2 = (F1(xval+h)-F1(xval-h))/(2*h)
%O(h^4) first derivative (centered)
FD_h4 = (-F1(xval+2*h)+8*F1(xval+h)...
    -8*F1(xval-h)+F1(xval-2*h))/(12*h)
% You can do the same for 2nd and higher
% order derivatives...
%
% Protein Data Example
D = xlsread('sample Data.xlsx','I3:I290');
% Change RFU to protein concentration
C_vec = 0.0043.*D + 0.6626; %ug/ml protein
% Calculate the rate of protein production vs. time
% For first point - use the forward finite difference
h = 5; % min
FirstSlope = (C_vec(2)-C_vec(1))/h;
% Last slope - use the backwards finite diff
LastSlope = (C_vec(end)-C_vec(end-1))/h;
% Middle points
Slopes = (C_vec(3:end)-C_vec(1:end-2))./(2*h);
% To plot all the slopes, concantenate the vector
S_vec = [FirstSlope Slopes' LastSlope];
np = length(S_vec);
t_vec = (1:np)*h;

% Hmm.. that looks interesting...
% Dr. Reuel will check it
% 2.8.19 Dr. Reuel check, I see the issue! Plot the raw data
figure
subplot(1,2,1)
title('Raw Data and Fit')
plot(t_vec,C_vec,'k+')
xlabel('Time (min)')
ylabel('Protein Conc (ug/ml)')
hold on
% See how ragged the data looks?  This variance is common in real data.
% We fit is with a model (we'll cover in class 18):
 FT = fittype('K*(1-(tau1*exp(-(x-theta)/tau1)-tau2*exp(-(x-theta)/tau2))/(tau1-tau2))*heaviside(x-theta)+Rbak');
 [p5,gof5] = fit(t_vec',C_vec,FT,'StartPoint',[500 70 160 114 222]);
 MODEL = @(t) p5.K*(1-(p5.tau1.*exp(-(t-p5.theta)./p5.tau1)-p5.tau2.*exp(-(t-p5.theta)./p5.tau2))./(p5.tau1-p5.tau2)).*heaviside(t-p5.theta)+p5.Rbak;
 % Smooth model data
 SMD = MODEL(t_vec);
 plot(t_vec,SMD,'r-')
 % Find the rates again:
 FirstSlope = (SMD(2)-SMD(1))/h;
 % Last slope - use the backwards finite diff
 LastSlope = (SMD(end)-SMD(end-1))/h;
 % Middle points
 Slopes = (SMD(3:end)-SMD(1:end-2))./(2*h);
 % To plot all the slopes, concantenate the vector
 S_vec2 = [FirstSlope Slopes LastSlope];
 %
 subplot(1,2,2)
 title('First Derivative (rate)')
 plot(t_vec,S_vec,t_vec,S_vec2);
xlabel('Time (min)')
ylabel('Rate of Protein Production (ug/(ml*min))')
legend('From Raw Data','From Smooth Model')
 %
 %
% Numerical derivative for not equal spaced points
% Use La Grange interpolation
F2 = @(x) 2*x^4-6*x^3-12*x-8;
% Use equation given in notes
x0 = -.5;
x1 = 1;
x2 = 2;
x = 0;
Der_0 = F2(x0)*(2*x-x1-x2)/((x0-x1)*(x0-x2))...
    +F2(x1)*(2*x-x0-x2)/((x1-x0)*(x1-x2))...
    +F2(x2)*(2*x-x0-x1)/((x2-x0)*(x2-x1))