% collision model for ice flow, inelastic collisions, produces "density"

function [pos,vel,den,t] = comp_icemodel(m,T,dt,x,u);
% m = number of points
% T = time length
% dt = timestep size
% x = initial ice positions
% u = initial ice velocity

tolr = exp(-25);    % tolerance

j = 1:m;            % m-length vectors
k = 1:m-1;

spc = 1/m;         % minimum spacing between blocks (block width)
cbs = zeros(1,m);   % on/off to mark blocks collided within timestep

% density
d(1)= spc/(x(2)-x(1));
d(m)= spc/(x(m)-x(m-1));
for r=2:m-1
    d(r) = 2*spc/(x(r+1)-x(r-1));
end

    n = 1; 
    t(n) = 0;
    pos(1,:)=x;
    vel(1,:)=u;
    den(1,:)=d;

while t(n) < T
    dif(k) =  x(k+1) + u(k+1)*dt - (x(k) + u(k)*dt) - spc;

    if any(dif < -tolr)
        for y = 1:m-1
            if u(y)-u(y+1) <= tolr
                tcvec(y) = 100*dt;
            else
                tcvec(y) = (x(y+1)-x(y)-spc)/(u(y)-u(y+1));     % tcvec stores time till collision between blocks
            end
        end
        [tcs,i] = sort(tcvec);                 % sorts tcvec, gives sorted tc's & their indices
        ind = i(find(abs(tcs-tcs(1))<tolr));   % returns lowest tcs' block indices
        
        holdx = x;
        x(j) = u(j)*tcs(1) + x(j);          % all the blocks moved until earliest collision time
            
        % density calculation
            d(1)= spc/(x(2)-x(1));
            d(m)= spc/(x(m)-x(m-1));
            for r=2:m-1
                d(r) = 2*spc/(x(r+1)-x(r-1));
            end
       
        for r=1:length(ind)     % this loop changes the velocities to reflect collisions
            q = ind(r);
            mar = 1;
            mal = 1;
            rx = q+1;
            lx = q;
            while rx < m
                rcon = x(rx+1) - x(rx) - spc + holdx(rx+1) - holdx(rx) - spc;
                if abs(rcon) < tolr
                    mar = mar + 1;
                    rx = rx+1;
                else
                    break;
                end
            end
            while lx > 1
                lcon = x(lx) - x(lx-1) - spc + holdx(lx) - holdx(lx-1) - spc;
                if abs(lcon) < tolr
                    mal = mal + 1;
                    lx = lx-1;
                else
                    break;
                end
            end                     
            lxr = lx:rx;
            holdx(lxr) = x(lxr);
            u(lxr) = (mal*u(q) + mar*u(q+1))/(mal+mar);     
        end
        t(n+1) = t(n) + tcs(1);
    else
        x(j) = x(j) + u(j)*dt;
        
        % density calculation
           d(1)= spc/(x(2)-x(1));
           d(m)= spc/(x(m)-x(m-1));
            for r=2:m-1
                d(r) = 2*spc/(x(r+1)-x(r-1));
            end
        t(n+1) = t(n) + dt;
    end
    n=n+1;
    pos(n,:)=x;
    vel(n,:)=u;
    den(n,:)=d;
end