clc, clear, close all;
data = readtable('ph_eq_data.xlsx', 'VariableNamingRule','preserve');
data = table2array(data);

set = [5e8; 5e9; 5e10; 5e11; 5e12; 5e13];
derivative = zeros(size(set,1),1);

for ci = 1:size(set,1)
    middle = ci * 2 + 1;
interpolatingdata = [data(middle - 2,:);data(middle - 1,:);data(middle,:);data(middle + 1,:)];

syms x;

y = 0;

for i = 1:size(interpolatingdata,1)
    p = sym(1);
    for j = 1:size(interpolatingdata,1)
        if i == j
            continue
        end
        p = p .* (x-interpolatingdata(j,1))/(interpolatingdata(i,1)-interpolatingdata(j,1));
    end
    y = y + p * interpolatingdata(i,2);
end

yp = diff(y);

ypm = matlabFunction(yp);

derivative(ci) = ypm(set(ci));
end

figure(1);
hold on;
plot(set, derivative)