function prestot=gasv(nprint,ntimes,nparticles,vmax)
% schrijf zelf commentaar
x=unifrnd(0,1,[1 nparticles]);
y=unifrnd(0,1,[1 nparticles]);
vx=unifrnd(-vmax,vmax,[1 nparticles]);
vy=unifrnd(-vmax,vmax,[1 nparticles]);
pres=zeros(5,1);
prestot=zeros(ntimes/nprint,5);
iprint=0;

for i=1:1:ntimes
    newx=step(x,vx,.1);
    newy=step(y,vy,.1);
    for j=1:1:nparticles
        if (~inbox(newx(j),newy(j),0,0,1,1))
[x(j),y(j),newx(j),newy(j),vx(j),vy(j),pres]= ...
correct(x(j),y(j),newx(j),newy(j),vx(j),vy(j),pres,0,0,1,1);
        end
    end
    x=newx;
    y=newy;
    if (~mod(i,nprint))
        iprint=iprint+1;
        prestot(iprint,:)=pres(:);
        prestot(iprint,5)=mean((abs(pres(1:4)))); %calculate mean pressure over 4 walls
        pres=zeros(5,1);
    end
end
        plot(prestot(:,5));
        xlabel('Time')
        ylabel('Pressure(Pa)')
        title('Pressure')


function q=step(p,v,deltat)
q=p+deltat*v;

function out=inbox(x,y,llx,lly,urx,ury)
if (x>llx&&x<urx&&y>lly&&y<ury)
    out=1;
else
    out=0;
end
function [x,y,newx,newy,vx,vy,pres]= ...
correct(x,y,newx,newy,vx,vy,pres,llx,lly,urx,ury)
while (~inbox(newx,newy,llx,lly,urx,ury))
    if (newx>urx)
        f=(newy-y)*(urx-x)/(newx-x);
        if(y+f>lly)&&(y+f<ury)
            y=y+f;
            x=urx;
            newx=urx-(newx-urx);
            vx=-vx;
            pres(1)=pres(1)+2*vx;
        end
    end
    if (newx<llx)
        f=(newy-y)*(llx-x)/(newx-x);
        if(y+f>lly)&&(y+f<ury)
            y=y+f;
            x=llx;
            newx=llx-(newx-llx);
            vx=-vx;
            pres(2)=pres(2)+2*vx;
        end
    end
    if (newy>ury)
        f=(newx-x)*(ury-y)/(newy-y);
        if(x+f>llx)&&(x+f<urx)
            x=x+f;
            y=ury;
            newy=ury-(newy-ury);
            vy=-vy;
            pres(3)=pres(3)+2*vy;
        end
    end
    if (newy<lly)
        f=(newx-x)*(lly-y)/(newy-y);
        if(x+f>llx)&&(x+f<urx)
            x=x+f;
            y=lly;
            newy=lly-(newy-lly);
            vy=-vy;
            pres(4)=pres(4)+2*vy;
        end
    end
end
