function yNew=reflective_boundary_conditions(y)
global nElectrons X Y Z
yRow=y';
xe=yRow(1:3:3*nElectrons);
ye=yRow(2:3:3*nElectrons);
ze=yRow(3:3:3*nElectrons);
vx=yRow(3*nElectrons+1:3:6*nElectrons);
vy=yRow(3*nElectrons+2:3:6*nElectrons);
vz=yRow(3*nElectrons+3:3:6*nElectrons);
for i=1:nElectrons
    if (xe(i)>X)
        xe(i)=X; %reflective boundary conditions
        vx(i)=-vx(i); %reflecting reverses velocity
    end
    if (xe(i)<0)
        xe(i)=0;
        vx(i)=-vx(i);
    end
    if (ye(i)>Y)
        ye(i)=Y;
        vy(i)=-vy(i);
    end
    if (ye(i)<0)
        ye(i)=0;
        vy(i)=-vy(i);
    end
    if (ze(i)>Z)
        ze(i)=Z;
        vz(i)=-vz(i);
    end
    if (ze(i)<0)
        ze(i)=0;
        vz(i)=-vz(i);
    end
end
r=[xe ye ze];
v=[vx vy vz];
for ir=0:nElectrons-1
    yr(3*ir+1:3*ir+3)=r(ir+1:nElectrons:nElectrons*3);
end
for iv=0:nElectrons-1
    yv(3*iv+1:3*iv+3)=v(iv+1:nElectrons:nElectrons*3);
yNewRow=[yr yv];
yNew=yNewRow';
end

