all clear
format short g

% x = [341 1506 10481];
% y = [359.26 317.6 273.71];
x = [341 3643 7126 10481];
y = [359.26 295.93 281.48 273.71];
% x = [341 661 1381 3828 10481];
% y = [359.26 339.26 319.26 294.82 273.71];
% x = [341 590 1071 2060 4451 10481];
% y = [359.26 342.59 325.93 309.26 291.48 273.71];
% x = [341 546 908 1574 2858 5466	10481];
% y = [359.26	344.82 330.37 315.93 301.48 287.04 273.71];

% for i = 1:2,
%      Q1(i) = (x(1)-x(i+1))/(y(1)-y(i+1));
% end
%Q2(i) = (x(2)-x(3))/(Q1(1)-Q1(2));


for i = 1:3,
    Q1(i) = (x(1)-x(i+1))/(y(1)-y(i+1));
end
for i = 1:2,
    Q2(i) = (x(2)-x(i+2))/(Q1(1)-Q1(i+1));
end
Q3(1) = (x(3)-x(4))/(Q2(1)-Q2(2));


% for i = 1:4,
%     Q1(i) = (x(1)-x(i+1))/(y(1)-y(i+1));
% end
% for i = 1:3,
%     Q2(i) = (x(2)-x(i+2))/(Q1(1)-Q1(i+1));
% end
% for i = 1:2,
%     Q3(i) = (x(3)-x(i+3))/(Q2(1)-Q2(i+1));
% end
% Q4(1) = (x(4)-x(5))/(Q3(1)-Q3(2));


% for i = 1:5,
%     Q1(i) = (x(1)-x(i+1))/(y(1)-y(i+1));
% end
% for i = 1:4,
%     Q2(i) = (x(2)-x(i+2))/(Q1(1)-Q1(i+1));
% end
% for i = 1:3,
%     Q3(i) = (x(3)-x(i+3))/(Q2(1)-Q2(i+1));
% end
% for i = 1:2,
%     Q4(i) = (x(4)-x(i+4))/(Q3(1)-Q3(i+1));
% end
% Q5(1) = (x(5)-x(6))/(Q4(1)-Q4(2));


% for i = 1:6,
%     Q1(i) = (x(1)-x(i+1))/(y(1)-y(i+1));
% end
% for i = 1:5,
%     Q2(i) = (x(2)-x(i+2))/(Q1(1)-Q1(i+1));
% end
% for i = 1:4,
%     Q3(i) = (x(3)-x(i+3))/(Q2(1)-Q2(i+1));
% end
% for i = 1:3,
%     Q4(i) = (x(4)-x(i+4))/(Q3(1)-Q3(i+1));
% end
% for i = 1:2,
%     Q5(i) = (x(5)-x(i+5))/(Q4(1)-Q4(i+1));
% end
% Q6(1) = (x(6)-x(7))/(Q5(1)-Q5(2));

%--------------------------------------------------------------------------
r = [341
354
366
379
393
407
422
438
454
471
489
507
526
546
567
590
612
636
661
688
715
744
774
805
838
872
908
946
986
1027
1071
1116
1164
1215
1267
1323
1381
1442
1506
1574
1645
1719
1799
1881
1968
2060
2157
2258
2366
2479
2598
2725
2858
2998
3146
3303
3468
3643
3828
4023
4232
4451
4684
4929
5190
5466
5759
6070
6400
6753
7126
7521
7941
8388
8863
9369
9908
10481];

t = [359.26
358.15
357.04
355.93
354.82
353.71
352.59
351.48
350.37
349.26
348.15
347.04
345.93
344.82
343.71
342.59
341.48
340.37
339.26
338.15
337.04
335.93
334.82
333.71
332.59
331.48
330.37
329.26
328.15
327.04
325.93
324.82
323.71
322.59
321.48
320.37
319.26
318.15
317.04
315.93
314.82
313.71
312.59
311.48
310.37
309.26
308.15
307.04
305.93
304.82
303.71
302.59
301.48
300.37
299.26
298.15
297.04
295.93
294.82
293.71
292.59
291.48
290.37
289.26
288.15
287.04
285.93
284.82
283.71
282.59
281.48
280.37
279.26
278.15
277.04
275.93
274.82
273.71];
%--------------------------------------------------------------------------

xv = r;

% p1 = (Q1(1).*y(1) + (xv-x(1)));
% q1 = Q1(1);
% p2 = (Q2(1).*p1 + (xv-x(2)).*y(1));
% q2 = (Q2(1).*q1 + (xv-x(2)));
% fx = p2./q2;

p1 = (Q1(1).*y(1) + (xv-x(1)));
q1 = Q1(1);
p2 = (Q2(1).*p1 + (xv-x(2)).*y(1));
q2 = (Q2(1).*q1 + (xv-x(2)));
p3 = (Q3(1).*p2 + (xv-x(3)).*p1);
q3 = (Q3(1).*q2 + (xv-x(3)).*q1);
fx = p3./q3;

yv = 1./(1.4e-3 + 2.37e-4*log(xv) + 9.9e-8*log(xv).^3); %Thermistor
% p1 = (Q1(1).*y(1) + (xv-x(1)));
% q1 = Q1(1);
% p2 = (Q2(1).*p1 + (xv-x(2)).*y(1));
% q2 = (Q2(1).*q1 + (xv-x(2)));
% p3 = (Q3(1).*p2 + (xv-x(3)).*p1);
% q3 = (Q3(1).*q2 + (xv-x(3)).*q1);
% p4 = (Q4(1).*p3 + (xv-x(4)).*p2);
% q4 = (Q4(1).*q3 + (xv-x(4)).*q2);
% fx = p4./q4;

% p1 = (Q1(1).*y(1) + (xv-x(1)));
% q1 = Q1(1);
% p2 = (Q2(1).*p1 + (xv-x(2)).*y(1));
% q2 = (Q2(1).*q1 + (xv-x(2)));
% p3 = (Q3(1).*p2 + (xv-x(3)).*p1);
% q3 = (Q3(1).*q2 + (xv-x(3)).*q1);
% p4 = (Q4(1).*p3 + (xv-x(4)).*p2);
% q4 = (Q4(1).*q3 + (xv-x(4)).*q2);
% p5 = (Q5(1).*p4 + (xv-x(5)).*p3);
% q5 = (Q5(1).*q4 + (xv-x(5)).*q3);
% fx = p5./q5;

% p1 = (Q1(1).*y(1) + (xv-x(1)));
% q1 = Q1(1);
% p2 = (Q2(1).*p1 + (xv-x(2)).*y(1));
% q2 = (Q2(1).*q1 + (xv-x(2)));
% p3 = (Q3(1).*p2 + (xv-x(3)).*p1);
% q3 = (Q3(1).*q2 + (xv-x(3)).*q1);
% p4 = (Q4(1).*p3 + (xv-x(4)).*p2);
% q4 = (Q4(1).*q3 + (xv-x(4)).*q2);
% p5 = (Q5(1).*p4 + (xv-x(5)).*p3);
% q5 = (Q5(1).*q4 + (xv-x(5)).*q3);
% p6 = (Q6(1).*p5 + (xv-x(6)).*p4);
% q6 = (Q6(1).*q5 + (xv-x(6)).*q4);
% fx = p6./q6;

err = 100*(yv-t)./t;
%plot(x,y,'o',xv,fx);
plot(xv,err);

x = [3 4 5 6 7 8 9];
y = [1.1981 1.1443 0.08299 0.032542 0.0207 0.0103 0.007];