function varma1(a,b,w,v)
    %Viljum adeins loglegar staerdir
    if(b <= a)
        disp('Armur a verdur ad vera styttri en armur b');
        return
    elseif(w < 0)
        w = abs(w);
    elseif(v < 0)
        disp('R?mm?l brunaholfsins verdur ad vera staerra en 0');
        return
    end

    % Eiginleikar lofts
    cp = 1.006;
    cv = 0.7177;
    n = cp/cv;

    % Eiginleikar brunans
    H_comb = 47.8; %Brennsluvarmi oktans i kJ/g, reiknad ut fra 5460 kJ/mol
    m_air = 0.001225*(0.5 + v); %density*volume (i litrum)
    M_air = 0.029; %molmassi = kg/mol
    n_air = m_air/M_air; %molfjoldi lofts
    %C_8H_18 + 12.5O_2 = 47.8 kJ/g

    n_octane = n_air/12.5; %molfjoldi oktans
    M_octane = 114.23; %g/mol
    m_octane = n_octane*M_octane; %massi oktans
    Q = m_octane*H_comb*0.4*0.21; %Oxygen 21% of air, 40% efficiency

    % Reikna max horn alpha ? TDC/BDC me? P??ag?ras
    max_alpha = atan(a/sqrt(b^2-a^2));
    w1 = w;
    w2 = (max_alpha/pi)*w;

    % Breytingar
    dtheta = 2*pi/360;
    dV = 0.0005/180; %Notum m3 i stadinn fyrir l

    % Faersla stimpils
    dt = dtheta/w;
    dx = 2*a*(w/(2*pi))*dt;
    S = 2*a;

    % Horn i TDC (I upphafi)
    alpha = 0;
    beta = pi;
    dalpha = max_alpha/pi;
    dbeta = dtheta;

    % Astand 1

    W = 0;
    for i = 1:1
        for theta = linspace(0,2*pi,360)
            alpha = asin((sin(theta)*a)/b);
            phi = pi - alpha - theta; %Hornid milli a og b
            varmaplot1(a,b,beta,S,v) %Hreyfimynd
            if(theta-pi<0.01 && theta-pi>0)
                % ORKA BRUNA - ORKA INN I KERFI
                T2 = (Q+m_air*cv*T(length(T)))/(m_air*cv); %Lokahiti
                P2 = (P(length(P))/T(length(T)))*T2; %Lokathrystingur vid fast rummal
                P = [P P2];
                V = [V V(length(V))];
                T = [T T2];
            elseif(theta == 2*pi)
                P = [P P(1)];
                V = [V V(length(V))];
                T = [T T(1)];
            elseif(theta >= 0 && theta <= pi-dtheta)
                % Samthjoppun
                if(theta == 0)
                    P = [100];
                    V = [0.0005+(v/1000)];
                    R = 0.2871;
                    T = [(P(1)*V(1))/(m_air*R)];
                else
                    vt = V(length(V))-dV;
                    P = [P P(length(P))*(V(length(V))/vt)^n];
                    T = [T (P(length(P))*V(length(V)))/(m_air*R)];
                    V = [V vt];
                    S = S - dx;
                end
            elseif(theta >= pi && theta <= 2*pi)
                % Uthensla
                vt = V(length(V))+dV;
                P = [P P(length(P))*(V(length(V))/vt)^n];
                T = [T ((P(length(P))*vt))/(m_air*R)];
                V = [V vt];
                S = S + dx;
            end

            if(alpha == max_alpha || alpha == -max_alpha)
                dalpha = -dalpha;
            end
            alpha = alpha + dalpha;
            beta = beta - dbeta;
            W = W + abs(P(length(P))*(0.5e-3/a)*sin(phi));
        end
    end
    fprintf('Shaft work: %.2f kJ \n', W)
    close all

    plot(V,T)
    xlim([(v/1000)*0.5 1.2*max(V)])
    ylim([0 1.2*max(T)])
    xlabel('Volume (m3)')
    ylabel('Temperature (K)')
    title('TV graph')

    figure
    plot(V,P)
    xlim([(v/1000)*0.5 1.2*max(V)])
    ylim([0 1.2*max(P)])
    xlabel('Volume (m3)')
    ylabel('Pressure (kPa)')
    title('PV graph')

end