{NOTE this file only contains the essential procedures with the new formulas presented in:
 Roodbergen, K.J., and Vis, I.F.A. (2006), A model for warehouse layout, IIE Transactions 38(10),799-811.
 This file is therefore not a ready-to-run program. You must create a shell first that can call the procedures below.
 All variable names are chosen consistent with the notation in the article (except for LG[i] which is denoted as Di 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 SumC(n,m : integer; y,wc : extended) : extended;

Var
  SumCtemp,X : extended;
  g,i,ii     : integer;

Begin
  SumCtemp:=0;
  for g:=1 to n do if (g mod 2 <> 0) and (g<=m)
  Then begin
    X:=0;
    For i:=1 To g-1 Do
    Begin
      if i mod 2 = 0 then ii:=-1 else ii:=1;
      X:=X+ii*bico(g,g-i)*power((g-i)/g,m);
    End;
    X:=1-X;
    SumCtemp:=SumCtemp+bico(n,g)*power(g/n,m)*X*(2*y*((m/g)/((m/g)+1))-y);
  End;
  SumC:=SumCtemp;
End;


{---------------------------------------------------------------------------}

Procedure Sshape_formula(S       : extended;
                         n       : integer;
                         d       : extended;
                         m       : integer;
                         y,wc,wa : extended;
                         var T   : extended);

var
  i : integer;
  EA,EDy,EDx1,EDx2 : extended;

Begin
  if n=1 then EA:=1
         else EA:=n*(1-power((n-1)/n,m));
  EDy:=(y+wc)*EA+sumC(n,m,y,wc);

  EDx1:=0;
  for i:=1 to n do
    EDx1:=EDx1+(abs(i-d)+abs(i-n+d-1))*(power(i/n,m)-power((i-1)/n,m));
  EDx2:=n-1;
  for i:=1 to n-1 do
    EDx2:=EDx2-2*power(i/n,m);
  T:=EDy+wa*(EDx1+EDx2);
End;

{---------------------------------------------------------------------------}

Procedure Lgap_formula(S       : extended;
                       n       : integer;
                       d       : extended;
                       m       : integer;
                       y,wc,wa : extended;
                       var T   : extended);

var
  i : integer;
  EA2,temp,EDy,EDx1,EDx2 : extended;
  E,LG : array[0..50] of real;

Begin
{note it is better to move the next 51 lines elsewhere so you only have to call this once}
  LG[0]:=  0;
  LG[1]:=0.499771;  E[1]:=1.000000;
  LG[2]:=0.777786;  E[2]:=1.332917;
  LG[3]:=0.958318;  E[3]:=1.499837;
  LG[4]:=1.086759;  E[4]:=1.600193;
  LG[5]:=1.183417;  E[5]:=1.666608;
  LG[6]:=1.259226;  E[6]:=1.714524;
  LG[7]:=1.320537;  E[7]:=1.749988;
  LG[8]:=1.371426;  E[8]:=1.777744;
  LG[9]:=1.414194;  E[9]:=1.800016;
  LG[10]:=1.450892; E[10]:=1.818164;
  LG[11]:=1.482815; E[11]:=1.833412;
  LG[12]:=1.510678; E[12]:=1.846085;
  LG[13]:=1.535439; E[13]:=1.857141;
  LG[14]:=1.557522; E[14]:=1.866621;
  LG[15]:=1.577467; E[15]:=1.875213;
  LG[16]:=1.595420; E[16]:=1.882601;
  LG[17]:=1.611719; E[17]:=1.888940;
  LG[18]:=1.626579; E[18]:=1.894618;
  LG[19]:=1.640204; E[19]:=1.900006;
  LG[20]:=1.652786; E[20]:=1.904615;
  LG[21]:=1.664514; E[21]:=1.909254;
  LG[22]:=1.675317; E[22]:=1.913237;
  LG[23]:=1.685335; E[23]:=1.916682;
  LG[24]:=1.694703; E[24]:=1.919979;
  LG[25]:=1.703493; E[25]:=1.923019;
  LG[26]:=1.711752; E[26]:=1.925929;
  LG[27]:=1.719490; E[27]:=1.928514;
  LG[28]:=1.726739; E[28]:=1.930860;
  LG[29]:=1.733616; E[29]:=1.933349;
  LG[30]:=1.740216; E[30]:=1.935555;
  LG[31]:=1.746336; E[31]:=1.937472;
  LG[32]:=1.752165; E[32]:=1.939296;
  LG[33]:=1.757786; E[33]:=1.941219;
  LG[34]:=1.763056; E[34]:=1.942816;
  LG[35]:=1.768070; E[35]:=1.944392;
  LG[36]:=1.772884; E[36]:=1.945919;
  LG[37]:=1.777472; E[37]:=1.947348;
  LG[38]:=1.781877; E[38]:=1.948592;
  LG[39]:=1.786086; E[39]:=1.949936;
  LG[40]:=1.790084; E[40]:=1.951239;
  LG[41]:=1.793939; E[41]:=1.952350;
  LG[42]:=1.797688; E[42]:=1.953575;
  LG[43]:=1.801250; E[43]:=1.954658;
  LG[44]:=1.804646; E[44]:=1.955654;
  LG[45]:=1.807956; E[45]:=1.956453;
  LG[46]:=1.811147; E[46]:=1.957352;
  LG[47]:=1.814230; E[47]:=1.958241;
  LG[48]:=1.817174; E[48]:=1.959106;
  LG[49]:=1.820044; E[49]:=1.959975;
  LG[50]:=1.822789; E[50]:=1.960857;

  if ((n=1) or (m=1))
  then EA2:=0
  else EA2:=(n*(1-power((n-1)/n,m))-power(1/n,m-1))/(1-power(1/n,m-1));

  temp:=0;
  for i:= 1 to m do temp:=temp+ bico(m,i)*power(1/n,i)*power((n-1)/n,m-i)*(y*LG[i]+wc*E[i]);
  EDy:=power(1/n,m-1)*(2*y*m/(m+1)+wc)
       + (1-power(1/n,m-1))*2*(y+wc)
       + (1-power(1/n,m-1))/(1-power((n-1)/n,m))*(EA2-2)*temp;

  EDx1:=0;
  for i:=1 to n do
    EDx1:=EDx1+(abs(i-d)+abs(i-n+d-1))*(power(i/n,m)-power((i-1)/n,m));
  EDx2:=n-1;
  for i:=1 to n-1 do
    EDx2:=EDx2-2*power(i/n,m);
  T:=EDy+wa*(EDx1+EDx2);
End;