---

% ЛБ1, часть 1, Вариант 5 (без локальных функций)
clear; clc; close all;

%% ================= ЗАДАЧА 1 =================
x  = [10.2 10.4 10.6 10.8 11.0];   % Sol
y  = [90 20 0 45 84];              % T, град C
xq = 10.5;                         % марсианский полдень

%% 1.1 Геометрическое представление табулированной функции
figure(1);
subplot(2,2,1); plot(x, y, 'ro-', 'LineWidth', 1.5, 'MarkerSize', 8); grid on;
xlabel('Sol'); ylabel('T, ^\circC'); title('plot(x,y)');
subplot(2,2,2); bar(x, y, 0.6); grid on;
xlabel('Sol'); ylabel('T, ^\circC'); title('bar(x,y)');
subplot(2,2,3); line(x, y, 'LineWidth', 1.5, 'Color', 'b'); grid on;
xlabel('Sol'); ylabel('T, ^\circC'); title('line(x,y)');
subplot(2,2,4); stem(x, y, 'filled', 'LineWidth', 1.5); grid on;
xlabel('Sol'); ylabel('T, ^\circC'); title('stem(x,y)');

%% 1.2 Интерполяционный полином Ньютона
n = length(x);
F = zeros(n);
F(:,1) = y(:);
for j = 2:n
    for i = 1:(n-j+1)
        F(i,j) = (F(i+1,j-1) - F(i,j-1)) / (x(i+j-1) - x(i));
    end
end
a = F(1,:);                        % разделённые разности
fprintf('Разделённые разности (коэф. Ньютона):\n'); disp(a);

% Полином Ньютона в символьном виде (выражение для отчёта)
% (нужен Symbolic Math Toolbox; если его нет — закомментируй этот блок)
syms s
Ps = a(1); w = 1;
for k = 2:n
    w = w*(s - x(k-1));
    Ps = Ps + a(k)*w;
end
Ps = expand(Ps);
fprintf('Полином Ньютона (развёрнутый вид):\n'); disp(vpa(Ps, 6));

% Численный расчёт полинома Ньютона на плотной сетке (без local function)
Xq = linspace(min(x), max(x), 201);
Ynew = zeros(size(Xq));
for idx = 1:numel(Xq)
    P = a(1); w = 1;
    for k = 2:n
        w = w*(Xq(idx) - x(k-1));
        P = P + a(k)*w;
    end
    Ynew(idx) = P;
end

% Значение в марсианский полдень Sol = 10.5
P = a(1); w = 1;
for k = 2:n
    w = w*(xq - x(k-1));
    P = P + a(k)*w;
end
Tnew = P;
fprintf('T(10.5) по полиному Ньютона = %.4f град C\n', Tnew);

%% 1.3 interp1: все способы интерполяции
Ynear = interp1(x, y, Xq, 'nearest');
Ylin  = interp1(x, y, Xq, 'linear');
Yspl  = interp1(x, y, Xq, 'spline');
Ypc   = interp1(x, y, Xq, 'pchip');

Tnear = interp1(x, y, xq, 'nearest');
Tlin  = interp1(x, y, xq, 'linear');
Tspl  = interp1(x, y, xq, 'spline');
Tpc   = interp1(x, y, xq, 'pchip');

fprintf('T(10.5): nearest = %.4f | linear = %.4f | spline = %.4f | pchip = %.4f\n', ...
        Tnear, Tlin, Tspl, Tpc);

figure(2);
plot(x, y, 'ro', 'MarkerSize', 9, 'LineWidth', 2); hold on;
plot(Xq, Ynew,  'k-',  'LineWidth', 1.8);
plot(Xq, Ynear, 'c--', 'LineWidth', 1.2);
plot(Xq, Ylin,  'b-.', 'LineWidth', 1.2);
plot(Xq, Yspl,  'g-',  'LineWidth', 1.2);
plot(Xq, Ypc,   'm:',  'LineWidth', 1.5);
plot(xq, Tnew, 'kp', 'MarkerSize', 13, 'LineWidth', 2);
grid on; xlabel('Sol'); ylabel('T, ^\circC');
title('Вариант 5: сравнение методов интерполяции');
legend('Табличные данные', 'Ньютон', 'nearest', 'linear', 'spline', ...
       'pchip', 'T(10.5)', 'Location', 'best');

%% ================= ЗАДАЧА 2 (двумерная интерполяция) =================
% Матрица с рисунка (строки сверху вниз)
Ztop = [121  99  52  40  59  123;
        114  97  49  36  42   78;
         69  52  34  39  75   80;
         20  15  12  10  11   15;
         29  25  20  21  25   29;
        136  98  71  65  52   38];
Z = flipud(Ztop);                  % строка 1 матрицы Matlab <-> минимальный y

h  = 2.5;                          % единичный отрезок сетки, см
nx = size(Z, 2); ny = size(Z, 1);
[Xg, Yg] = meshgrid((0:nx-1)*h, (0:ny-1)*h);

figure(3);
subplot(2,3,1);
surf(Xg, Yg, Z); shading interp; colorbar;
xlabel('x, см'); ylabel('y, см'); title('Табулированная функция');

% Более мелкая сетка для интерполяции
step = h/5;
[Xi, Yi] = meshgrid(0:step:(nx-1)*h, 0:step:(ny-1)*h);
ZN = interp2(Xg, Yg, Z, Xi, Yi, 'nearest');
ZL = interp2(Xg, Yg, Z, Xi, Yi, 'linear');
ZC = interp2(Xg, Yg, Z, Xi, Yi, 'cubic');
ZS = interp2(Xg, Yg, Z, Xi, Yi, 'spline');

subplot(2,3,2); surf(Xi, Yi, ZN); shading interp; colorbar; title('nearest');
subplot(2,3,3); surf(Xi, Yi, ZL); shading interp; colorbar; title('bilinear');
subplot(2,3,4); surf(Xi, Yi, ZC); shading interp; colorbar; title('bicubic');
subplot(2,3,5); surf(Xi, Yi, ZS); shading interp; colorbar; title('spline');
subplot(2,3,6); meshc(Xi, Yi, ZS); colorbar; title('Линии уровня (meshc)');
xlabel('x, см'); ylabel('y, см');