harj55.m

function harj55()

    % Ladataan data
    data = load('data51.txt');

    % Datan koko, sarakkeet (N) ovat x-koordinaatit, rivit (M) y-koordinaatit
    [M,N] = size(data);

    Imax = max(max(data));   % Suurin alkio
    Ifiltteri = 0.05*Imax;   % Suodatusraja
    
    % Muutetaan alkiot jotka pienempiä kuin suodatusraja nollaksi
    for i = 1:N
        for j = 1:M
            if data(j,i) <= Ifiltteri
                data(j,i) = 0;
            end  
        end
    end

    % Edellä oleva silmukka voidaan myös korvata yhdellä käskyllä:
    %data(data<=Ifiltteri) = 0;

    Isum = sum(sum(data)); % Kaikkien alkioiden summa

    % Pisteiden koordinaatit millimetrinä ja milliradiaaneina
    x = linspace(-29.5,30,N);
    y = linspace(-87.3,89.7,M);

    % Lasketaan odotusarvot
    x2 = 0;
    y2 = 0;
    xy = 0;
    for i = 1:N
        for j = 1:M
            x2 = x2 + data(j,i) * x(i)^2;
            y2 = y2 + data(j,i) * y(j)^2;
            xy = xy + data(j,i) * x(i)*y(j);
        end
    end
    x2 = x2/Isum;
    y2 = y2/Isum;
    xy = xy/Isum;

    % Silmukat voidaan korvata myös meshgrid-rakenteella:
    %[xx,yy] = meshgrid(x,y);
    %x2 = sum(sum(data.*xx.^2))/Isum;
    %y2 = sum(sum(data.*yy.^2))/Isum;
    %xy = sum(sum(data.*xx.*yy))/Isum;
    
    epsilon = sqrt(y2*x2-xy^2);
    fprintf( 'Emittanssi on %g mm mrad\n', epsilon )
end