% model for ice flow, inelastic collisions, implicit relaxation

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

tolr = exp(-25);    % tolerance

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

spc = .001;         % 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
    w = wv(x,t,j,n);
     
    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));
            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
%             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)
            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);
            cbs(lxr) = 1;         
        end
        if cbs(j) == 0;
            u(j) = (u(j)+alpha*dt*w(j))/(1+alpha*dt);
        end
        cbs(j) = 0;
        t(n+1) = t(n) + tcs(1);
    else
        x(j) = x(j) + u(j)*dt;
%            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
        u(j) = u(j) - alpha*dt*(u(j)-w(j));
        t(n+1) = t(n) + dt;
    end
    n=n+1;
    pos(n,:)=x;
    vel(n,:)=u;
%     den(n,:)=d;
end