%This script calculates the Taylor dispersivity due to flow through a wavy
%channel.  The lower wall is at y=0 and the upper wall is at
%h=1+e*cos(2*pi*x/w) where e is the amplitude of the waviness and w is the
%width of the unit cell.

e=.2;
w=5;

%We begin by introducing all of our particles:

n=20000;

nstart=round(n*(1+2*e));  %A bit more than n particles, as some are discarded.

x=w*rand(nstart,1); %The x locations
y=(1+e)*rand(nstart,1); %The y locations

%Now we need to determine the values which are inside our domain:

i=find(y<(1+e*cos(2*pi*x/w)));

x=x(i);
y=y(i);

%These are all the valid particles.  We really only want n of them:

x=x(1:n);
y=y(1:n);

%We can plot this up:
figure(1)
xp=w*[0:.01:1];
plot(w*[0,1],[0,0,],'k',xp,1+e*cos(2*pi*xp/w),'k',x,y,'ob')
xlabel('x')
ylabel('y')
axis('equal')
title('Initial distribution')
drawnow

%Now for the simulation:

tlim=0.3*w^2;
dt=0.001;

t=[dt:dt:tlim]'; %We have a column vector of times
varz=zeros(size(t));
zbar=zeros(size(t));

z=zeros(size(x));  %We start the particles at zero...

for j=1:length(t)
    %We have the random steps:
    dx=randn(n,1)*(2*dt)^.5;
    dy=randn(n,1)*(2*dt)^.5;
    
    xt=x+dx;
    yt=y+dy;
    
    %Now we check to see if particles have bounced out:
    
    ilow=find(yt<0); %the particles below the bottom wall
    ihigh=find(yt>(1+e*cos(2*pi*xt/w))); %the particles above the top wall
    
    yt(ilow)=abs(yt(ilow)); %Reflect those back!
    
    a=yt(ihigh)-(1+e*cos(2*pi*xt(ihigh)/w));
    theta=atan(-e*2*pi/w*sin(2*pi*xt(ihigh)/w));
    yt(ihigh)=(1+e*cos(2*pi*xt(ihigh)/w))-a.*cos(2*theta);
    xt(ihigh)=xt(ihigh)+a.*sin(2*theta);
    
    %Now for getting dz.  We use simpson's rule:
    h=1+e*cos(2*pi*x/w);
    ht=1+e*cos(2*pi*xt/w);
    
    %We calculate the midpoint of the jump:
    xm=(x+xt)/2;
    ym=(y+yt)/2;
    hm=1+e*cos(2*pi*xm/w);
    
    %and we get dz via Simpson's rule:
    dz=1/6*(6*y.*(h-y)+24*ym.*(hm-ym)+6*yt.*(ht-yt))*dt;
    
    %Now we deal with the change in z for particles bouncing off the walls:
    f=-y(ilow)./dy(ilow);
    dz(ilow)=0.5*dt*(f.*6.0.*y(ilow).*(h(ilow)-y(ilow))+(1-f)*6.0.*yt(ilow).*(ht(ilow)-yt(ilow)));
    
    f=(h(ihigh)-y(ihigh))./((h(ihigh)-y(ihigh))+a); %The value of a was calculated before!
    dz(ihigh)=0.5*dt*(f.*6.0.*y(ihigh).*(h(ihigh)-y(ihigh))+(1-f)*6.0.*yt(ihigh).*(ht(ihigh)-yt(ihigh)));
   
    %Move x values back into the unit cell:
    ileft=find(xt<0);
    xt(ileft)=w+xt(ileft);
    iright=find(xt>w);
    xt(iright)=xt(iright)-w;
    
    %We update z and get the statistics:
    z=z+dz;
    x=xt;
    y=yt;
    
    %now we populate the columns of zbar and varz:
    zbar(j)=mean(z);
    varz(j)=var(z);
    
    %We do a little plotting to see what is going on: comment this out
    %later to speed things up!
    
%     figure(1)
%     xp=w*[0:.01:1];
%     plot(w*[0,1],[0,0,],'k',xp,1+e*cos(2*pi*xp/w),'k',x,y,'ob')
%     axis('equal')
%     xlabel('x')
%     ylabel('y')
%     title(['x and y positions for t = ',num2str(t(j))])
%     
%     figure(2)
%     plot(x,z,'o')
%     xlabel('x')
%     ylabel('z')
%     title(['x and z positions for t = ',num2str(t(j))])
%     drawnow
    
end

%and we get the value of the Taylor dispersivity from the slope of the
%variance.  We fit the last 1/3 of the values:
ifit=[round(length(t)*2/3):length(t)];
tfit=t(ifit);
varzfit=varz(ifit);

xx=[ones(size(tfit)),tfit]\varzfit;

k=xx(2)/2; %The Taylor dispersivity

disp('The ratio of the Taylor dispersivity to that of a channel with no side-walls is:')
out=k*210  %The ratio to the channel theoretical result

figure(2) %We can comment this out if we wish.
plot(t,varz,t,[ones(size(t)),t]*xx,t,t/105)
legend('simulation','linear fit','e=0 theory','Location','NorthWest')
xlabel('t')
ylabel('variance in z')
title(['Variance for e = ',num2str(e),' and W = ',num2str(w)])
grid on
drawnow

    figure(3)
    xp=w*[0:.01:1];
    plot(w*[0,1],[0,0,],'k',xp,1+e*cos(2*pi*xp/w),'k',x,y,'ob')
    axis('equal')
    xlabel('x')
    ylabel('y')
    title(['x and y positions for t = ',num2str(t(j))])
    
    figure(4)
    plot(x,z,'o')
    xlabel('x')
    ylabel('z')
    title(['x and z positions for t = ',num2str(t(j))])
    drawnow