====== Zajęcia 1 ======
<code matlab>n = input('Podaj dlugosc ');
disp(n); % tylko wartosc zmiennej
display(n); % nazwa zmiennej i wartosc

% alokujemy pamiec tworzac tablice skladajaca sie z samych zer
% n,1 = wymiary
% podajac jeden argument = nxn
% mozna wiecej dodawac wymiarow po przecinku
tab = zeros(n,1);

% tab = 1:1:n stworzy wektor od 1 z krokiem 1 do n

 %i ma sie zmieniac od 1, z krokiem 1, do n (wykona sie dla n)
for i = 1:1:n
    
    % tablice sa indeksowane od 1
    tab(i)=i;
end

% brak { }. Po wszystkim zawsze end

% if
% else
% end

% wypisujemy tablice
disp(tab);

i = 1
while i <= n
    tab(i) = i;
    i = i + 1;
end

% jezeli na koncu nie dodamy srednika wypisze na ekranie np efekt

% tablice od -pi z krokiem 0.1 do pi
x = -pi:0.1:pi;
y = sin(x);
plot(x, y);

figure %otwiera nowe okno
plot(y);

% mnozymy y x y
%y = y * y; % nie mozna mnozyc przes siebie dwoch takich samych wektorow
% y = y * y.' % transpozycja [wynikiem jest skalar] wektor poziomy razy pionowy
y = y.' * y % transpozycja .' [wynikiem jest tablica] wektor pionowy razy poziomy
surf(y);

% Ax = b
% x =   A^{-1} * b
% x =   A \ b       <=> A^{-1} * b
%       b / A       <=> b * A^{-1}</code>
====== Zajęcia 2 ======
===== Trójkąt Pascala =====
<code matlab>n = input('podaj wymiar: ');
T = cell(n, 1);

for i = 1:1:n
    T{i} = ones(1, i);
    
    for j = 2:1:i-1
        T{i}(j) = T{i-1}(j-1) + T{i-1}(j);
    end
end

for i = 1:1:n
    disp(T{i});
end</code>
===== Całka metodą prostokątów =====
<code matlab>s = 0.1;
x = -1:s:1;
y = x.^2;

c = 0;

% end = length(y)
%for i = 1 : end - 1
for i = 1 : length(y) - 1
    c = c + y(i) * s;
end

display(c)

% lub
d = sum(s * y(1: end - 1))</code>
===== Całka metodą trapezów =====
<code matlab>s = 0.1;
x = -1:s:1;
y = x.^2;

% lub
c = sum( (y(1:end-1) + y(2:end)) * s/2 );

display(c)

% lub
c = trapz(x, y);

display(c)</code>

===== Całka metodą Monte Carlo =====
<code matlab>Pp = 2;
Nf = 0;
N = 1000000;
x = 2 * rand(N, 1) - 1;
y = rand(N, 1);

for i = 1:1:N
    if y(i) <= x(i)^2
        Nf = Nf + 1;
    end
end

c = (Nf / N) * Pp;

display(c)

% lub
w = x.^2 > y;
%nnz = number of non zero elements
Nf = nnz(w);

c = (Nf / N) * Pp;

display(c)</code>
====== Zajęcia 3 ======
===== Macierze i wektory =====
<code matlab>w = rand(1,10);
% co trzeci element wektora zamieniamy na odwrotny
w(1:3:end) = -w(1:3:end);
disp(w);

% odwroc elementy mniejsze od 0.5
w = rand(1,10);
w(w < 0.5) = -w(w < 0.5);
disp(w);

% & - dla tablic, wektorów
% && - dla skalarow
% odwracamy elementy wieksze od 0.3 i mniejsze od 0.6
w = rand(1,10);
w(w > 0.3 & w < 0.6) = -w(w > 0.3 & w < .6)
disp(w);

%w = w.*.2; % mnozenie tablicowy, kazdy element mnozymy * 0.2 (0 mozna
%pominac) | 0.2 = .2</code>

<code matlab>% randi - liczba calkowita od 0,1
w = randi(10,1,5);

%mean - srednia
% srednia bez wartosci najwiekszej i najmniejszej
%sr = mean(w(w > min(w) & w < max(w)));
%disp(sr);

% srednia z zabezpieczeniem odrzucenia np. dwoch maksymalnych elementow
%sr = (sum(w) - min(w) - max(w)) / length(w) - 2;
%disp(sr);

disp(w);

% wyswietl ile jest liczb parzystych, a ile nieparzystych
a = [ length(w(mod(w,2)==0)), nnz(mod(w,2))];
disp(a);

% randn - liczba losowa o rozkladzie normalnym o wartosci oczekiwanej od 0,1
w = randn(1,10);

% dodanie wszystkich ujemnych, dodanie wszystkich dodatnich
a = [ sum(w(w<0)), sum(w(w>0))];
disp(a);
</code>
<code matlab>% macierz 3x5
n = randi(10,3,5);
disp(n);

% dodalismy numery wierszy do macierzy
%n = [(1:3).' n];
%disp(n);

% dodalismy numery kolumn do macierzy
%n = [0:5; n];
%n = [0:size(n,2)-1 ; n];
%disp(n);

% ; - nowy wiersz w macierzy
% [ a,b ] = macierz pozioma
% [ a ; b ] = macierz pionowa

%n = [0:5 ; [(1:3).' n]];
% lub 
%n = [0:size(n,2) ; [(1:size(n,1)).' n]];
%disp(n);

% usuwamy kolumne macierzy w ktorej wystepuje wartosc mniejsza od 3
%n(:,i) - i-ta kolumna macierzy
%disp(n);

%n(:, min(n) < 3) = [];
% lub
%n(:, ~(min(n) > 3)) = [];
% lub
%n = n(:, min(n) > 3);

% usuwamy wiersz macierzy w ktorej wystepuje wartosc mniejse od 6
%n(min(n.').' < 3, :) = [];
%disp(n);

% pomijamy drugi argument, szukaj w wierszach(2) lub w kolumnach(1) min(n, [], 2)
%n(min(n, [], 2) < 3, :) = [];
%disp(n);

% znajdz wszystkie wartosci minimalne i zamien na NaN;
% n(min(n, [], 2) == min(min(n)), min(n) == min(min(n))) = NaN;
% disp(n);</code>
===== Sortowanie =====
<code matlab>% macierz 3x5
n = randi(10,3,5);
disp(n);

%[a b] = sort(n)
% a = posortowana macierz
% b = kolejnosc jak on to zmienia

% sortuje tylko pierwsza kolumne (i zwraca do a, n jest po staremu)
[ a b ] = sort(n(:, 1));
% lub ~pomija argument [tylko w nowych wersjach]
%[ ~, b ] = sort(n(:, 1));

% sortujemy tylko pierwsza kolumne
%n(:,1) = n(b,1);

%sortujemy wszystkie kolumny zgodnie z kluczem pierwszej kolumny
n = n(b,:);

% sortujemy pierwszy wiersz
[a b ] = sort(n(1,:))
% sortujemy kolumny wg klucza z poprzedniego sortowania
n = n(:, b)

disp(n);</code>
====== Zajęcia 4 ======
===== Funkcja całki z wielomianu =====
==== calka.m ====
<code matlab>
function wy = calka(w, a, b)
    if (nargin == 1)
        wy = polyint(w);
    elseif (nargin == 2)
        wy = polyint(w, a);
    else
        tmp = polyint(w);
        wy = polyval(tmp, b) - polyval(tmp, a);
    end
end
</code>
==== main.m ====
<code matlab>
p = [ 1 0 0 ];
%x^4 + 5x^3 - 2x^2 + 0x + 3

%polyint(p); %calka z wielomianu
%polyint(p, k); %warunek couchiego drugi artument, wartosc wyrazu wolnego

n = calka(p, 0, 1);
disp(n);
</code>
===== Miejsca zerowe i ekstrema na wykresie =====
<code matlab>
% zamknij wszystkie otwarte wykresy
close all;

%p = [ 1 -2 -1 2 4 ];
p = conv(conv(conv([1 1], [1 -1]), [1, 0]), [1 2]);
x = -2:0.01:2;

y = polyval(p, x);
dp = polyder(p); % wielomian ktory jest pochodna wielomianu p

%plot(x, y);
%hold on;
%plot(r, zeros(size(r)), 'rx');
% mozna tak jak wyzej lub jeden plot na dole

%roots - funkcja zwraca pierwiastki wielomianu
r = roots(p);
g = roots(dp);
%plot(r, zeros(size(r)), 'rx', x, y);

% polyval(p,g) - wartosc wielomianu p w punktach g
plot(g, polyval(p, g), 'X', r, zeros(size(r)), '*', x, y);
</code>
===== Ekstrema funkcji =====
<code matlab>
% zamknij wszystkie otwarte wykresy
close all;

% szukamy elementow ktore nie sa ekstremami, ich druga pochodna musi byc = 0
p = [ 1 0 0 0 ];
x = -2:0.01:2;

y = polyval(p, x);
z = roots(p);

dp = polyder(p);
rdp = roots(dp);

dp2 = polyder(dp);
rdp2 = roots(dp2);

for i=1:length(rdp2)
   rdp( rdp == rdp2(i) ) = []; 
end

disp(rdp);
disp(rdp2);

plot(rdp, polyval(p, rdp), 'X', z, zeros(size(z)), 'r*', x, y);
</code>
===== Praca domowa =====
<code matlab>
p1 = [1 0 0];
p2 = 1;

% funkcja kwadratowa
% funkcja stala

% napisac program ktory dostaje dwa wieomiany p1 i p2 i zwraca wspolrzedne
% punktow przeciecia

% roznica dwoch wielomianow, i liczymy miejsca zerowe tej roznicy
% wielomianow
</code>
====== Zajęcia 5 ======
Kolokwium
====== Zajęcia 6 ======
===== GUI =====

<code matlab>

function pushbutton1_Callback(hObject, eventdata, handles)
   wiersze = str2double(get(handles.edit1,'String'));
   kolumny = str2double(get(handles.edit2,'String'));
   A = ones(wiersze, kolumny, 3);
   
   n = str2double(get(handles.edit3,'String'));
   for i = 1:1:n
    x = randi(wiersze);
    y = randi(kolumny);
    rgb = rand(1,3);
    A(x,y,:) = rgb;
   end
   set(handles.figure1,'UserData', A);
   imshow(A);

% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
    set(handles.pushbutton2,'UserData', 1);
    A = get(handles.figure1,'UserData');
    B = get(handles.figure1,'UserData');
    n = str2double(get(handles.edit1,'String'));
    m = str2double(get(handles.edit2,'String'));
    todo = true;
    while todo && get(handles.pushbutton2,'UserData')==1
        todo = false;
        for i = 2:1:n-1
           for j = 2:1:m-1
               if(all(B(i,j,:)==1))
                   if(~all(A(i-1,j,:)==1))
                       B(i,j,:) = A(i-1,j,:) ;
                       todo = true;
                   end
                   if(~all(A(i,j-1,:)==1))
                       B(i,j,:) = A(i,j-1,:);
                       todo = true;
                   end
                   if(~all(A(i,j+1,:)==1))
                       B(i,j,:) = A(i,j+1,:);
                       todo = true;
                   end
                   if(~all(A(i+1,j,:)==1))
                       B(i,j,:) = A(i+1,j,:);
                       todo = true;
                   end
               end
           end
        end
        set(handles.figure1,'UserData', A);
        A = B;
        imshow(A);
        pause(0.03);
    end
    
    set(handles.pushbutton2,'UserData', 0);

</code>

====== Zajęcia 7 ======
===== fun_celu =====
<code matlab>
function y=fun_celu(X,r, a, b)
    N = size(X,1)/2;
    X = reshape(X, N, 2);
    suma = 0;
    for i = 1:N
        if X(i,1)<=a+r || X(i,1)>=b-r || X(i,2)<=a+r || X(i,2)>=b-r
            suma = suma +1000;
        end
        for j = i+1:N
            odl = norm(X(i,:)-X(j,:));
            if odl<2*r
                suma = suma +(2*r - odl);
            end
        end
    end
    y = suma;
end
</code>

===== main =====
<code matlab>

a = 0;
b = 10;
N = 10;
r = 1.5;
k = 0:0.01:2*pi;
X = rand(N,2)*(b-a-2*r)+a+r;

for i = 1:N
    x = r * cos(k) + X(i,1);
    y = r * sin(k) + X(i,2);
    fill(x,y,'r','FaceAlpha',0.5);
    hold on;
end
set(gca,'Xlim',[a,b]);
set(gca,'Ylim',[a,b]);
axis equal;

X = reshape(X, 2*N, 1);
[x_opt y_opt] = fminsearch(@(x)fun_celu(x,r,a,b), X);

X = reshape(x_opt, N, 2);
figure;
for i = 1:N
    x = r * cos(k) + X(i,1);
    y = r * sin(k) + X(i,2);
    fill(x,y,'g','FaceAlpha',0.5);
    hold on;
end
set(gca,'Xlim',[a,b]);
set(gca,'Ylim',[a,b]);
axis equal;


</code>
====== Zajęcia 8 ======
===== 1 =====
Obiekt zawieszony na sprężynie
<code matlab>m = 1;
b = 0.4;
k = 0.5;

dx = @(t, x)[x(2);...
    -(k*x(1)+b*x(2))/m];
    
%x(0) = 1
%x'(0) = 0
[t,x] = ode45(dx, 0:0.1:100, [1 0]);

plot(t, x(:,1));</code>
===== 2 =====
==== fun_celu ====
<code matlab>
function ts = fun_celu(m, X)
    dt = 0.05;
    dx = @(t, x)[x(2);-(X(1)*x(1)+X(2)*x(2))/m];

    [t,x] = ode45(dx, 0:0.1:100, [1 0]);
    pol = x(:,1);
    
    for i=size(pol):-1:0
        if pol(i) > dt || pol(i) < -dt
            ts = t(i);
            break;
        end
    end
end
</code>
==== main ====
<code matlab>
m = 1;
b = 0.4;
k = 0.5;

X = [k;b];

[x_opt y_opt] = fmincon(@(X)fun_celu(m, X),[0.5;0.4],[],[],[],[],[0;0],[5;5]);

display(y_opt);
display(x_opt);

%znalezc czas stabilizacji zeby k,b -> ts -> 0
%do optymalizacji potrzebujemy funkcji celu
</code>
====== Zajęcia 9 ======
===== Ciężarek z zajęć nr 8 =====
<code matlab>simulink

plot(x);</code>
{{:studia:ciezarek.png?400|}}
==Zawieszenie samochodu==
{{:studia:zawieszenie.png?400|}}
====== Zajęcia 10 ======
===== Wariant 1 =====
==== main ====
<code matlab>
n = 10;
m = 1;
b1 = 0.1;
b2 = 0.5;
F = 1;

%a1 = x1;
%a2 = x1';
%a3 = x2;
%a4 = x2';

da = @(t, a)[a(2); (sila(t)-b1*a(2)-b2*(a(2)-a(4)))/n;...
             a(4); b2*(a(2)-a(4))/m];
 
%                  [0 10] - od 0 do 10 sekund, dla takiego czasu dostajemy rozwiazanie         
%                          [0 - poczatkowe polozenie wozka, 0 - pocztkowa predkosc wozka, 0 - poczatkowe polozenie klocka, 0 - poczatkowa predkosc klocka]
[t, a] = ode45(da, [0 10], [0 0 0 0]);
%                          [a1 a2 a3 a4]

plot(t, a);</code>
==== sila ====
<code matlab>
function [ f ] = sila( t )
    if t < 5
        f = 1;
    else
        f = -1;
    end
end
</code>
===== Wariant 2 =====
==== main ====
<code matlab>

M = 10;
m = 1;
b1 = 0.5;
b2 = 0.5;
T = 10;
F = 1;



[F, odl]= fmincon(@FC,2*rand(3,1)-1,[],[],[],[],[-2;-2;-2],[2;2;2]);


 display(F);
 display(odl);
 
 da = @(t,a)[a(2);(sila(t,F)-b1*a(2)-b2*(a(2)-a(4)))/M ;...
     a(4);b2*(a(2)-a(4))/m];
 
 [t,a] = ode45(da,[0 30], [0 0 0 0]);

 display(a);
 plot(t,a);

</code>

==== sila ====
<code matlab>


function f=sila(t,F)
    if t<10
        f=F(1);
    elseif t<20 && t>=10
        f=F(2);
    else
        f=F(3);
    end
end
</code>

==== fc ====
<code matlab>
function fc = FC(F)
M = 10;
m = 1;
b1 = 0.5;
b2 = 0.5;

da = @(t,a)[a(2);(sila(t,F)-b1*a(2)-b2*(a(2)-a(4)))/M ;...
    a(4);b2*(a(2)-a(4))/m];
[t,a] = ode45(da,[0 30], [0 0 0 0]);

fc = (a(end,1)-5)^2+ 10*a(end,2)^2; %ograniczenia na przejechana droge a(end,1)-5)^2 i na to zeby predkosc na koncu byla 0

if max(abs(a(:,1)-a(:,1))) >0.1     % ograniczenie zeby klocek nie spadl
    fc = 1000;
end
end

</code>
====== Zajęcia 11 ======
===== Calka =====
==== main ====
<code matlab>p = [ 1 0 0 ];

% liczenie analityczne
% calka               , 1 - przedzial calkowania (od 0 do 1)
c = polyval(polyint(p), 1) - polyint(polyder(p), 0);
disp(c);

%liczenie calki numerycznie
% c1 = CALKA[0 1] x^2 dx
c1 = quad(@(x)polyval(p,x), 0, 1);
disp(c1);

%liczenie calki numerycznie ktorej funkcja podcalkowa jest w fun.m
p = [ 1 0 ];
L = quad(@(x)fun(p, x), 0, 1);
disp(L);</code>
==== fun.m ====
Całka służąca do liczenia długości odcinka
<code matlab>function y  = fun( p, x )
    dp = polyder(p);
    y = sqrt(1 + polyval(dp, x).^2);  
end
</code>
===== Calka 2 =====
<code matlab>
x = 0:10;
y = rand(size(x));
arr = zeros(length(x)-1, 2);

% zwraca wspolczynniki wielomianu
for i = 1:length(x)-1
    %dopasowanie wielomianu
    p = polyfit(x, y, i);
    arr(i, 1) = quad(@(x)fun(p, x), 0, 1);
    arr(i, 2) = i;
    
    figure(i);
    plot(p);
end

bar(arr(:, 1));
</code>
===== Funkcja sklejalna =====
==== main ====
<code matlab>%dopasowanie funkcji sklejanej
% s = spline(x, y);
% y = ppval(s, x);

x = 0:10;
y = x.^2;

s = spline(x ,y);
sum = 0;
for i = 1:s.pieces
    p = s.coefs(i, :);   
    
    sum = sum + quad(@(x)fun(p, x), 0, 1);
end

disp(sum);</code>
==== fun.m ====
Całka służąca do liczenia długości odcinka
<code matlab>function y  = fun( p, x )
    dp = polyder(p);
    y = sqrt(1 + polyval(dp, x).^2);  
end
</code>
====== Zajęcia 12 ======
<code matlab>T = 1;
fs = 8192;
t = 0:1/fs:T;
n = length(t);

y1 = 2 * sin(2 * pi * t * 200);

y2 = sin (2 * pi * t * 600);

plot(t(1:200), y1(1:200), t(1:200), y2(1:200));

y = y1 + y2;
figure;
plot(t(1:200), y(1:200));

% widmo sygnalu
% fft - fast fourier transform
w = abs(fft(y))/n;
% f - wektor czestotoliwosci
f = (0 : n/2) * fs / n;

bar (f, w(1:length(f)));

% filtrowanie sygnalu
% 1szy arg - rzad filtru
% 2gi  arg - czestotoliwosc odciecia (musi byc z zakresu [0 1]
% 3gi  arg - 'low' - przepuszczaj czestotliwosci ponizej
%            'high' - przepuszczaj czestotliwosci powyzej
[l,m] = butter(12, 2 * 400 / fs, 'low');
y3 = filter(l, m, y);
figure;
plot(t(1:200), y3(1:200), t(1:200), y1(1:200));

% rysowanie charakterystyki filtra
[h, theta] = freqz(l, m);
figure;
plot(theta * fs/(2 * pi), abs(h));

% dodamy szum do y1
y1 = y1 + randn(size(y1));
figure;
plot(t(1:200), y1(1:200));

% widmo sygnalu zaszumionego
s = abs(fft(y1))/n;
% find(s > 0.5)
figure;
plot(f, s(1:length(f)));

% probujemy odfiltrowac szum
% dwa filtry (dolnoprzepustowy [wszystko ponizej 220] + gornoprzepustowy [wszystko powyzej 180])
% [l1, m1] = butter(12, 2 * 180 / h, 'high');
% [l2, m2] = butter(12, 2 * 220 / h, 'low');
% y4 = filter(l1, m1, y1);
% y5 = filter(l2, m2, y4);
% figure;
% plot(t(1:200), y5(1:200));

% symulacja klawiatury telefonicznej (klawisz 5)
f1 = 1336;
f2 = 770;
y = sin(2 * pi * t * f1) + sin ( 2 * pi * t * f2);
sound(y * fs);</code>
====== Zajęcia 13 ======
===== Korelacja =====
<code matlab>T = 1;
fs = 8192;
t = 0:1/fs:T;
y1 = sin(2 * pi * t * 200);
y2 = sin (2 * pi * t * 500);
y = y1 + y2;

%wspolczynnik korelacji
% c = corr(y.', y.');
% disp(c);
% 
% x = 0:0.01:2*pi;
% y1 = sin(x);
% y2 = -sin(x);
% c1 = corr(y1.', y2.');
% disp(c1);
% 
% x = -1:0.01:1;
% y = x.^2;
% c2 = corr(x.', y.');
% disp(c2);

% for i=1:200
%    c(i) = corr(y(1:1000).', y(i:i+999).');
% end
% plot((0:199)/fs, c);</code>
==Dokładność aproksymacji wielomianowej==
<code matlab>% x = -1:0.1:1;
% y = x.^2 + 0.1 * randn(size(x));
% plot(x, y, 'r+');
% p = polyfit(x,y,2);
% 
% x1 = -1:0.01:1;
% y1 = polyval(p, x1);
% hold on;
% plot (x1, y1);

x = -1:0.1:1;
y = x.^2 + 0.1 * randn(size(x));
plot(x, y, 'r+');
[p S] = polyfit(x,y,2);

x1 = -1:0.01:1;
y1 = polyval(p, x1);
hold on;
plot (x1, y1);
[Y D] = polyconf(p, x, S);
plot(x, Y-D, 'r--');
plot(x, Y+D, 'g--');</code>

===== Funkcje statystyczne =====
<code matlab>x = [1 1 1 1 1 1 100];

%minimum
v_min = min(x);
v_max = max(x);
%srednia
v_mean = mean(x);
%srednia z x ale odrzucamy 20% pomiarow
%(10 najmniejszych wartosci, 20 najwiekszych)
v_trimmean = trimmean(x, 20);

%srednia geometryczna
v_geomean = geomean(x);

%srednia harmoniczna
v_harmmean = harmmean(x);

%mediana
v_median = median(x);

%moda (najczesciej wystepujaca wartosc)
v_mode = mode(x);

%zakres
v_range = range(x);

%odchylenie standardowe
v_std = std(x);

%wariancja
v_var = var(x);

% roznica pomiedzy 1szym a 3tym kwantylem x
v_iqr = iqr(x);</code>
===== Spadek z wysokości =====
<code matlab>%spadek z wysokosci
close all;
g = 9.81;
h0 = [100 200 300 400 500];
n = length(h0);
m = 20;
t = zeros(m, n);
v = zeros(m, n);

for i=1:n
   for j=1:m
      t(j,i) = sqrt(2 * h0(i) / g ) + 0.1 * randn();
      v(j,i) = sqrt(2 * h0(i) * g ) + randn();
   end
end

% wykres ramka-wasy
boxplot(t, h0);
figure;
boxplot(v, h0);
   
h1 = [];
for i = 1:n
    h1 = [h1; ones(m,1) * h0(i)];
end
t1 = reshape(t, n * m, 1);
v1 = reshape(v, n * m, 1);

figure;
gscatter(v1, t1, h1, 'rbgky', 'xo+sx');

t2 = t(:, 1);
v2 = v(:, 1);
figure;
gscatter(v2, t2);

% wyznaczenie najbardziej odstajacych, wyznacza kontury wewnatrz ktorych sa
% wszystkie punkty
k = convhull(v2, t2);
hold on;
plot(v2(k), t2(k));
hold off;

figure;

gscatter(v2, t2);
% wygenerowanie trojkatow
T = delaunay(v2, t2);
hold on;
triplot(T, v2, t2);
hold off;

[m1 s1] = normfit(v2);
</code>