function [RotZ]=calcRotZ(theta)
    RotZ = [cos(theta),-sin(theta),0,0; sin(theta),cos(theta),0,0; 0,0,1,0; 0,0,0,1]
endfunction

function [TransZ]=calcTransZ(d)
    TransZ = [1,0,0,0; 0,1,0,0; 0,0,1,d; 0,0,0,1]
endfunction

function [TransX]=calcTransX(a)
    TransX = [1,0,0,a; 0,1,0,0; 0,0,1,0; 0,0,0,1]
endfunction

function [RotX]=calcRotX(alfa)
    RotX = [1,0,0,0; 0,cos(alfa),-sin(alfa),0; 0,sin(alfa), cos(alfa),0; 0,0,0,1]
endfunction

//O1
theta = (0 * %pi / 180) //zmienny kąt
RotZ = calcRotZ(theta)
TransZ = calcTransZ(290)
RotX = calcRotX((%pi/2))
O1 = RotZ * TransZ * RotX

for i = -11:11

    //O2
    theta = ((i*10) * %pi / 180) //zmienny kąt
    RotZ = calcRotZ(theta)
    TransX = calcTransX(270)
    O2 = RotZ * TransX
    disp(i);

    for j = -11:7

    //O3
        theta = ((j*10) * %pi / 180) //zmienny kąt
        RotZ = calcRotZ(theta)
    TransX = calcTransX(270)
    TransX2 = calcTransX(70)
    RotX = calcRotX((%pi/2))
    TransZ = calcTransZ(302)

        O3 = RotZ * TransX * TransX2 * RotX * TransZ

    for j = -12:12

        //O4
        theta = (0 * %pi / 180) //zmienny kąt
        RotZ = calcRotZ(theta)
        TransX = calcTransX(72)
        RotX = calcRotX((%pi/2))
        O4 = RotZ * TransX * RotX

        T = O1 * O2 * O3 * O4
        //disp(T)

        plot(T(1,4), T(3,4), '.')

    end
end