function [x y D Gn nu Ek Ep Ewp nt C]=ppack(Lx,Ly,gam,Dn,ispack,plotit,ic,fx,fy)
% polypack pack particles in a lees-edwards cell
% Usage: [x y]=ppack(Lx,Ly,gam,Dn,ispack,plotit);
%
% Pack disks of diameter list Dn in a periodic lees-edwards cell
% defined by Lx x Ly at angle gam.

% revision history:
% 07/22/2010 Mark D. Shattuck <mds> md.m
%            MD Demo for HandsOn 2010
%            000 mds Initial conditions and visualization
% 07/24/2010 001 mds Add Euler
% 07/24/2010 002 mds Add Interaction detection and Force Law
% 07/24/2010 003 mds Add Velocity Verlet
% 07/24/2010 004 mds Add Nplotskip
% 08/22/2010 006 mds Add saving
% 10/11/2010 007 mds Add add fast force
% 11/18/2010 000 mds pack.m Packing finder redo of fastpack
% 04/30/2012 000 mds polypack.m pack in polygon
% 11/27/2012 000 mds added 'fix' ic and plotit
% 02/16/2013 001 mds add new stop and hunt
% 02/11/2019 000 mds convert to periodic BC

%% Experimental Parmeters

N=length(Dn);    % number particles
Ds=min(Dn);      % smallest particle
%Dl=max(Dn);      % largest particle
Gn=Dn/Ds;        % particle ratios
K=100;   % spring constant for harmonic force law
B=1;     % damping

Ptol=1e-6;    % potential energy per particle<Ptol in final state
Pthresh=1e-6; % potential per particle<Pthresh before growth
Pstep=1e-5;   % potential per particle for add in growth

if(~exist('ispack','var')|| isempty(ispack))
  ispack=true;
end

if(~exist('ic','var') || isempty(ic))
  ic='Rand';
end

Atot=Lx*Ly;

%% Display Parameters
if(~exist('plotit','var') || isempty(plotit))
  plotit=true;
end
Nplotskip=100;  % number of timesteps to skip before plotting

%% Simulation Parmeters
dt=.5e-1;        % time step
Nt=20000;       % max time per pack

%% Calculated Parameters
Gnm=bsxfun(@plus,Gn,Gn')/2;  % min separation matrix ratio Dnm=Ds*Gnm
size(Gnm)

%% Initial Conditions
switch ic
  case 'Rand'
    x=Lx*rand(1,N);
    y=Ly*rand(1,N);
  case 'fix'
    x=fx;
    y=fy;
end

% find initial minimum distance D

% find initial D
if(N>1)
  dx=repmat(x,N,1);
  dx=dx-dx';
  dx=dx-round(dx/Lx)*Lx;  % Periodic x
  dy=repmat(y,N,1);
  dy=dy-dy';
  dy=dy-round(dy/Ly)*Ly;  % Periodic x
  dnm=sqrt(dx.^2+dy.^2); % Distance between all particles
  dnm(1:N+1:end)=inf;    % set diagonal to D
  D=min(dnm(:)./Gnm(:));
else
  D=min(Lx,Ly);
end

Dn=D*Gn;        % Diameter list

vx=zeros(1,N);
vy=zeros(1,N);

ax_old=0*x;
ay_old=0*y;


%% Save Variables
Ek=zeros(Nt,1);    % Kinetic Energy
Ep=zeros(Nt,1);    % particle-particle potential
Ewp=zeros(Nt,1);   % wall-particle potential

%% Setup Plotting
clf;
Nh=N+2*ceil((Lx+Ly)*sqrt(N/Lx/Ly))+4;
h=zeros(1,Nh);
cc=jet(N);
for np=1:N
  h(np)=rectangle('Position',[x(np)-.5*D y(np)-.5*D D D],'Curvature',[1 1],'edgecolor','k','facecolor',cc(np,:));
end
for np=N+1:Nh
  h(np)=rectangle('Position',[2*Lx 2*Ly 1 1],'Curvature',[1 1],'edgecolor','k');
end
hold on;
plot([0 gam*Lx],[0 Ly],'k:');
axis('equal');
axis([0 Lx 0 Ly]);
%axis([-3 3 -3 3]);

%% Main Loop
Dh=-1;
Dl=-1;
nu=pi/(4*Atot)*sum(Dn.^2);
for nt=1:Nt
  % plot particles
  if(plotit && rem(nt-1,Nplotskip)==0)
    fmy=floor(y/Ly);
    xp=mod(x-fmy*gam,Lx);
    yp=mod(y,Ly);
    for np=1:N
      set(h(np),'Position',[xp(np)-.5*Dn(np) yp(np)-.5*Dn(np) Dn(np) Dn(np)]);
    end
    pmy=round(yp/Ly);
    dm=yp-pmy*Ly;
    ii=find(abs(dm)<Dn/2);
    Nc=N;
    Ny=length(ii);
    for ni=1:Ny
      np=ii(ni);
      set(h(Nc+ni),'Position',[xp(np)-.5*Dn(np)+sign(dm(np))*gam dm(np)+(1-pmy(np))*Ly-.5*Dn(np) Dn(np) Dn(np)],'facecolor',cc(np,:));
    end
    Nc=Nc+Ny;
    pmx=round(xp/Lx);
    dm=xp-pmx*Lx;
    ii=find(abs(dm)<Dn/2);
    Nx=length(ii);
    for ni=1:Nx
      np=ii(ni);
      set(h(Nc+ni),'Position',[dm(np)+(1-pmx(np))*Lx-.5*Dn(np) yp(np)-.5*Dn(np) Dn(np) Dn(np)],'facecolor',cc(np,:));
    end
    Nc=Nc+Nx;
    for np=Nc+1:Nh
      set(h(np),'Position',[-2 -2+np-Nc-1 1 1]);
    end
    drawnow;
    disp(nu);
  end
  
  Ek(nt)=sum((vx.^2+vy.^2))/2;
  
  x=x+vx*dt+ax_old.*dt.^2/2;  % first step in Verlet integration
  y=y+vy*dt+ay_old.*dt.^2/2;
  
  % Interaction detector and Force Law
  Fx=zeros(1,N);
  Fy=zeros(1,N);
  
  C=0;
  Gmin=9e9;
  for nn=1:N
    for mm=nn+1:N
      dy=y(mm)-y(nn);
      im=round(dy/Ly);
      dy=dy-im*Ly;  % Periodic x
      Dnm=(Dn(nn)+Dn(mm))/2;
      if(abs(dy)<Dnm)
        dx=x(mm)-x(nn);
        dx=dx-round(dx/Lx-im*gam)*Lx-im*gam*Lx;
        dnm=sqrt(dx.^2+dy.^2);
        Gmin=min(Gmin,dnm/Dnm);
        if(dnm<Dnm)
          C=C+1;
          F=-K*(Dnm/dnm-1);
          Ep(nt)=Ep(nt)+K/2*(Dnm-dnm).^2;
          Fx(nn)=Fx(nn)+F.*dx;  % particle-particle Force Law
          Fx(mm)=Fx(mm)-F.*dx;
          Fy(nn)=Fy(nn)+F.*dy;  % particle-particle Force Law
          Fy(mm)=Fy(mm)-F.*dy;
        end
      end
    end
  end
  
  Fx=Fx-B*vx;
  Fy=Fy-B*vy;
  
  if(Ep(nt)+Ewp(nt)<Pthresh*N && Dh<0)
    if(~ispack)
      if(D==Ds)
        break;
      end
    end
    D=D+sqrt(Pstep/K);
    if(~ispack)
      D=min(D,Ds);
    end
    Dn=D*Gn;
  elseif(Ep(nt)+Ewp(nt)<Ptol*N/100 && Dh>0)
    Dl=D;
    D=(Dl+Dh)/2;
    Dn=D*Gn;
  elseif(Ek(nt)<1e-20*N && C>N && Ep(nt)+Ewp(nt)>Ptol*N)
    if(Dh<0)
      Dh=D;
    end
    if(Dl>0)
      Dh=D;
      D=(Dl+Dh)/2;
    else
      D=D-sqrt(Pstep/K);
    end
    Dn=D*Gn;
  elseif(Ek(nt)<1e-25*N && C>N)
    break;
  end
  
  nu=pi/(4*Atot)*sum(Dn.^2);
  
  ax=Fx;
  ay=Fy;
  
  vx=vx+(ax_old+ax).*dt/2;  % second step in Verlet integration
  vy=vy+(ay_old+ay).*dt/2;
  
  ax_old=ax;
  ay_old=ay;
end

nu=pi/(4*Atot)*sum(Dn.^2);


%save(sprintf('%s%02d',basename,n,'.mat'),'x','y','D','Lx','Ly','nu','B','N','K','G','Ns');
fprintf(1,' %f\n',nu);


