% comparing continuum and collision models by density output

T = .5;             % total time
m = 20;             % number of blocks in collision model
    j = 1:m;
    
dx = 1/m;           % grid width for continuum model
dt = 1/125;         % time step for both models
mu = dt/dx;
c=(2/3);

% initial distribution of blocks, velocity for collision model
xst=[0:dx:1-dx];
x(j)= (1/c)*(xst(j)+dx/2); 
u(j) = sin((2*pi*c*x(j)));

    
[uh,p,rho,xh,x1h,th]=ice_simplex(dx,mu,T);
[pos,vel,den,te]=comp_icemodel(m,T,dt,x,u);

% translating variable rho into concentration/density
conc = 1./(1+rho);

% comparing initial densities to check
figure(1)
hold on
plot((x1h/c),conc(1,:),'b.-')
plot(pos(1,:),den(1,:),'g.')
title('initial densities')
xlabel('position')
ylabel('density')
axis([0 1.5 0 1])
legend('continuum model','collision model')


% comparing densities at point immediately at first collision
[q1,q11] = size(uh);
[q2,q22] = size(pos);

    % find that point on the continuum model
    for k=1:q1
        if any(p(k,:))
            contpoint=k;
            break;
        end
    end  
    % find that point on the collision model
    for k=1:q2     
        if den(k,(ceil(m/2))) > 1-exp(-15)
             collpoint=k;
             break;
        end
    end     

    % roughly integrate along the velocities in continuum model to find the
    % positions of grid points at that first collision timestep
    mpos = xh/c + (sum(uh(1:contpoint,:),1).*dt);

    % now we get the positions halfway between, in order to approximate the
    % grid points where rho is evaluated by the continuum model
    for k = 1:m-1
        mmid(k) = mpos(k)+((mpos(k+1)-mpos(k))/2);
    end
    mmid(m) = mpos(m)+(((1/c)-mpos(m))/2);

figure(2)
plot(mmid,conc(contpoint,:),'b.-')
hold on
plot(pos(collpoint,:),den(collpoint,:),'g.-')
title('densities immediately at first collision')
xlabel('position')
ylabel('density')
legend('continuum model','collision model')


% comparing densities at end of T
    % roughly integrate along the velocities to find the final positions of
    % grid points where velocity was computed
    fpos = xh/c + (sum(uh,1).*dt);

    % now we get the positions halfway between, in order to approximate the
    % grid points where rho is evaluated by the continuum model
    for k = 1:m-1
        fmid(k) = fpos(k)+((fpos(k+1)-fpos(k))/2);
    end
    fmid(m) = fpos(m)+(((1/c)-fpos(m))/2);

figure(3)
plot(fmid,conc(q1,:),'b.-')
hold on
plot(pos(q2,:),den(q2,:),'g.-')
title('densities after time T')
xlabel('position')
ylabel('density')
legend('continuum model','collision model')

% integrating densities to recover mass, a rough Reimann sum
% continuum model
    % the initial width of each bar is always dx
    % the width of each bar after time T, each point imagined to be in the
    % middle of a bar, except at end points, which stretch to 0 and 1/c
    for k = 2:m-1
        contbarf(k) = (fmid(k+1)-fmid(k-1))/2;
    end
    contbarf(1) = ((fmid(2)-fmid(1))/2) + fmid(1);
    contbarf(m) = (1/c)-fmid(m) + ((fmid(m)-fmid(m-1))/2);

    contdensumi = sum(conc(1,:)*(dx/c))
    contdensumf = sum(conc(q1,:).*contbarf)

% collision model
    % the initial width of each bar --- is always dx
    % the width of each bar after time T, each point imagined to be in the
    % middle of a bar, except at end points, which stretch to 0 and 1/c
    for k = 2:m-1
        collbarf(k) = (pos(q2,k+1)-pos(q2,k-1))/2;
    end
    collbarf(1) = ((pos(q2,2)-pos(q2,1))/2) + pos(q2,1);
    collbarf(m) = (1/c)-pos(q2,m) + ((pos(q2,m)-pos(q2,m-1))/2);

    colldensumi = sum(den(1,:)*(dx/c))
    colldensumf = sum(den(q2,:).*collbarf)