harj54.m

% Lasketaan kappaleen hitausmomentti käyttäen hilan kokoa h.
% Aloitetaan arvosta h=1 ja tihennetään hilaa. Esitetään tulos
% graafisesti.
function harj54()
    h = 1
    i = 1;
    while h > 0.01
        hh(i) = h;
        J(i) = hitausmomentti( h );
        i = i+1;
        h = 0.8*h;
    end

    semilogx( hh, J )
end


% Alifunktio joka laskee hitausmomentin käyttäen hilan kokoa h.
function J = hitausmomentti( h )
    L = 30;
    P = 20;
    xc = L/2; % Keskipisteen
    yc = P/2; % koordinaatit
    R = 5;

    J = 0;    % Hitausmomentti alustetaan nollaksi

    rho = 1.3;
    dm = h^2*rho;
    
    % Käydään läpi neliskanttisen levyn pisteet
    for x = 0:h:L
        for y = 0:h:P
            % Mikäli ei olla reiän kohdalla summataan pisteen
            % kontribuutio hitausmomenttiin
            if sqrt((x-xc)^2+(y-yc)^2) >= R
                J = J + dm*(x^2+y^2);
            end
        end
    end

    fprintf( 'Hilakoko on %g, hitausmomentti on %g g/cm^2\n', h, J )
end