{NOTE this file only contains the essential procedures with the new formulas presented in:
 Roodbergen, K.J., Sharp, G.P., and Vis, I.F.A. (2008), Designing the layout structure of 
 manual order-picking areas in warehouses. IIE Transactions, forthcoming.
 This file is therefore not a ready-to-run program. You must create a shell first that can
 call the procedures below. Most variable names are chosen consistent with the notation in 
 the article.}


function bico(n,k:integer):extended;

var
  b:extended;
  i,temp:integer;

begin
  if (k=0) or (n=k)
  then bico:=1
  else begin
    if k>n-k
    then begin
      b:=1;
      temp:=k;
      for i:=n downto n-k+1 do
      begin
        b:=b*i/temp;
        temp:=temp-1;
      end;
      bico:=b;
    End
    else begin
      b:=1;
      temp:=n-k;
      for i:=n downto k+1 do
      begin
        b:=b*i/temp;
        temp:=temp-1;
      end;
      bico:=b;
    End;
  end;
end;

function B(u,m:integer;p:extended):extended;
begin
  b:=bico(m,u)*power(p,u)*power(1-p,m-u);
end;


procedure Sshape_new_k_blocks(S:extended;n,k:integer;m:integer;wc,wa:extended;var formlength:extended);

var
  e:array[1..11] of extended;
  a:array[1..maxaisles,1..maxblocks] of extended;
  Eij3,Eij5,Eij6,Eij7,Eij8,Eij9,temp,sumL,p:extended;
  i,j,u,g,gg,h,l:integer;

begin
  y:=S/n;

{===============================EQUATION=1===============================}
  E[1]:=n*k*(1-power((n*k-1)/(n*k),m))*(y/k+wc);

{===============================EQUATION=2===============================}
  if (n mod 2 =1) and (m/(n*k)>1)
  then e[2]:=2*y*m/(m+n*k)-y
  else e[2]:=0;

{===============================EQUATION=3===============================}
  e[3]:=0;
  for i:=1 to n do for j:=1 to k do
  begin
    A[i,j]:=(power(i/n,m)-power((i-1)/n,m))*(power(j/k,m)-power((j-1)/k,m));
    Eij3:=(j-1)*power((i*j-1)/(i*j),m)+(j-1)*power((i*j-i+1)/(i*j),m);
    E[3]:=e[3]+A[i,j]*Eij3;
  end;
  E[3]:=(y/k+wc)*E[3];

{===============================EQUATION=4===============================}
  E[4]:=0;
  for i:=1 to n do
    E[4]:=E[4]+wa*(n-i)*(power(i/n,m)-power((i-1)/n,m));

{==============================EQUATION 5================================}
  E[5]:=0;
  for i:=2 to n do for j:=1 to k do
  begin
    Eij5:=0;
    for g:=1 to i-1 do
    begin
      temp:=0;
      for u:=1 to m do
        temp:=temp+B(u,m,(i-1)/(i*j))*(power((i-g)/(i-1),u)-power((i-g-1)/(i-1),u));
      Eij5:=Eij5+(i-g)*temp;
    end;
    E[5]:=E[5]+A[i,j]*Eij5;
  end;
  E[5]:=wa*E[5];

{==============================EQUATION 6 and 7===========================}
  E[6]:=0;
  E[7]:=0;
  for i:=3 to n do for j:=1 to k do
  begin
    Eij6:=0;
    for u:=1 to m-1 do
    begin
      suml:=0;
      for l:=1 to i-2 do
        SumL:=SumL+(i-1-l)*(power((l+1)/(i-1),u)-2*power(l/(i-1),u)+power((l-1)/(i-1),u))*l;
      SumL:=B(u,m,(i-1)/(i*j))*SumL;
      Eij6:=Eij6+SumL;
    end;
    Eij7:=(1/3)*((i-2)-Eij6)*(1-power((i*j-i+1)/(i*j),m));;
    E[6]:=E[6]+A[i,j]*(j-1)*Eij6;
    E[7]:=E[7]+A[i,j]*(j-1)*Eij7;
  end;
  E[6]:=wa*E[6];
  E[7]:=wa*E[7];


{=============================EQUATION 8 and 9===========================}
  E[8]:=0;
  E[9]:=0;
  for i:=1 to n do for j:=1 to k do
  begin
    p:=0;
    for h:=2 to j do if ((h mod 2 = 0) and (h<=m)) then
    begin
      temp:=0;
      for g:= 1 to h-1 do
      begin
        if g mod 2 = 0 then gg:=-1 else gg:=1;
        temp:= temp + gg*bico(h,g)*power((h-g)/h,m)
      end;
      p := p + bico(j,h)*power(h/j,m)*(1-temp);
    end;
    Eij8:=0;
    Eij9:=0;
    for g:=1 to i-1 do for u:=1 to m-1 do
    begin
      Eij8 := Eij8 + B(u,m,(i-1)/(i*j))*(n-g)*(power(g/(i-1),u)-power((g-1)/(i-1),u));
      Eij9 := Eij9 + B(u,m,(i-1)/(i*j))*(n-g)*(power((i-g)/(i-1),u)-power((i-g-1)/(i-1),u));
    end;
    E[8]:=E[8]+p*A[i,j]*Eij8;
    E[9]:=E[9]+(1-p)*A[i,j]*Eij9;
  end;
  E[8]:=wa*E[8];
  E[9]:=wa*E[9];

{==============================EQUATION 10===============================}
  E[10]:=0;
  for i:=1 to n do for j:=1 to k do
    E[10]:=E[10] + A[i,j]*(n-i/2)*power((i*j-i+1)/(i*j),m);
  E[10]:=wa*E[10];

{==============================EQUATION 11===============================}

  E[11]:=0;
  for i:=1 to n do for g:=1 to n do
    E[11]:=E[11]+A[i,1]*(g-1)*(power(g/n,m)-power((g-1)/n,m));
  E[11]:=wa*E[11];

{============================================================================}

  if k>1 then formlength:=E[1]+E[2]+E[3]+E[4]+E[5]+E[6]+E[7]+E[8]+E[9]+E[10]+E[11]
         else formlength:=E[1]+E[2]+2*(E[4]+E[5]);
  if m=1 then formlength:=y+k*wc+(n-1)*wa;
end;