Weddles

function IW = weddles(a, b, n)

    f = @(x) sqrt(1 - 0.162 * (sin(x) ^ 2));
    h = (b - a) / n;

    sum = f(a) + f(b);

    for i = 1 : n - 1
        md = mod(i, 6);
        func_val = f(a + i * h);

        if md == 1 || md == 5
            sum = sum + 5 * func_val;
        elseif md == 2 || md == 4
            sum = sum + func_val;
        elseif md == 3
            sum = sum + 6 * func_val;
        else
            sum = sum + 2 * func_val;
        end

    end

    IW = (3 * sum * h) / 10;
end