%externalroutine route(%realarray(1)%name ax,bx,%integer n,%real awid, %c
  bwid,twid,tsep,ay,%realname by,%integer v)
%externalrealfnspec sin(%real x)
%externalrealfnspec cos(%real x)
%routinespec genpolys(%integer i,%realarray(1)%name ax,bx, %c
  %realarray(3)%name apx,apy,bpx,bpy,%integerarray(2)%name apv,bpv, %c
  %integerarray(1)%name apn,bpn,%real atw,btw)
%routinespec getsep(%integer i,%realarray(3)%name apx,apy,bpx,bpy,%c
  %integerarray(2)%name apv,bpv,%integerarray(1)%name apn,bpn)
%realarray px,py(0:n-1,1:v),apx,apy,bpx,bpy(1:n,0:n-1,1:v+1)
%integerarray apn,bpn,sign(1:n),apv,bpv(1:n,0:n-1)
%real atw,btw,maxsep,d,r,s,t
%constreal pi=3.14159265
%integer i,j,k,l,m
  atw=(awid-twid)/2
  btw=(bwid-twid)/2
  d=twid+tsep
  ! generate polygon vertex positions
  r=cos(pi/(4*v))
  %for i=1,1,v %cycle
    s=d*cos((2*i-1)*pi/(4*v))/r
    t=d*sin((2*i-1)*pi/(4*v))/r
    %for j=0,1,n-1 %cycle
      px(j,i)=j*s
      py(j,i)=j*t
    %repeat
  %repeat
  ! determine ax:bx
  %for i=1,1,n %cycle
    %if ax(i)=bx(i) %then sign(i)=0 %and %continue
    %if ax(i)<bx(i) %then sign(i)=1 %else sign(i)=-1
  %repeat
  ! generate all boundary polygons
  %for i=1,1,n %cycle    ;! for each pair of ports
    %if sign(i)=0 %then %continue
    %if sign(i)=1 %then %c
      genpolys(i,ax,bx,apx,apy,bpx,bpy,apv,bpv,apn,bpn,atw,btw) %else %c
      genpolys(i,bx,ax,bpx,bpy,apx,apy,bpv,apv,bpn,apn,btw,atw)
  %repeat
  ! get minimum separation possible
  maxsep=0
  %for i=1,1,n %cycle    ;! for each pair of ports
    %if sign(i)=0 %then %continue
    %if sign(i)=1 %then getsep(i,apx,apy,bpx,bpy,apv,bpv,apn,bpn) %c
      %else getsep(i,bpx,bpy,apx,apy,bpv,apv,bpn,apn)
  %repeat
  by=ay+maxsep
  %return
!----------------------------------------------------------------------
%routine genpolys(%integer i,%realarray(1)%name ax,bx, %c
  %realarray(3)%name apx,apy,bpx,bpy,%integerarray(2)%name apv,bpv, %c
  %integerarray(1)%name apn,bpn,%real atw,btw)
%integer j,k,l,m
  m=0
  %for j=i,1,n %cycle
    %for k=1,1,v %cycle
      apx(i,m,k)=ax(j)-atw-px(m,k)
      apy(i,m,k)=py(m,k)
      %if apx(i,m,k)>bx(i) %then %start
        %if k=1 %then apn(i)=m-1 %and ->dob
        apv(i,m)=k
        ->na
      %finish
    %repeat
    apx(i,m,v+1)=ax(j)-atw-px(m,v)
    apy(i,m,v+1)=py(m,v)
    apv(i,m)=v+1
na: m=m+1
  %repeat
  apn(i)=n-i
dob:m=0
  %for j=i,-1,1 %cycle
    %for k=1,1,v %cycle
      bpx(i,m,k)=bx(j)+btw+px(m,k)
      bpy(i,m,k)=py(m,k)
      %if bpx(i,m,k)<ax(i) %then %start
        %if k=1 %then bpn(i)=m-1 %and %return
        bpv(i,m)=k
        ->nb
      %finish
    %repeat
    bpx(i,m,v+1)=bx(j)+btw+px(m,v)
    bpy(i,m,v+1)=py(m,v)
    bpv(i,m)=v+1
nb: m=m+1
  %repeat
  bpn(i)=i-1
%end
!----------------------------------------------------------------------
%routine getsep(%integer i,%realarray(3)%name apx,apy,bpx,bpy, %c
  %integerarray(2)%name apv,bpv,%integerarray(1)%name apn,bpn)
%real sep,t1,t2,t3
%integer ap,av,bp,bv
  %for ap=0,1,apn(i) %cycle
    %for av=1,1,apv(i,ap)-1 %cycle
      t1=apx(i,ap,av)
      %for bp=0,1,bpn(i) %cycle
        %if bp#0 %and bpx(i,bp,v+1)<=t1<=bpx(i,bp,1) %then %start
          %for bv=2,1,v+1 %cycle
            t2=bpx(i,bp,bv)
            %if t2<=t1 %then %exit
          %repeat
          t3=(t1-t2)/(bpx(i,bp,bv-1)-t2)
          sep=apy(i,ap,av)+bpy(i,bp,bv)*(1-t3)+bpy(i,bp,bv-1)*t3
        %finish %else sep=apy(i,ap,av)
        %if sep>maxsep %then maxsep=sep
      %repeat
    %repeat
  %repeat
  %for bp=0,1,bpn(i) %cycle
    %for bv=1,1,bpv(i,bp)-1 %cycle
      t1=bpx(i,bp,bv)
      %for ap=0,1,apn(i) %cycle
        %if ap#0 %and apx(i,ap,1)<=t1<=apx(i,ap,v+1) %then %start
          %for av=2,1,v+1 %cycle
            t2=apx(i,ap,av)
            %if t2>=t1 %then %exit
          %repeat
          t3=(t2-t1)/(t2-apx(i,ap,av-1))
          sep=bpy(i,bp,bv)+apy(i,ap,av)*(1-t3)+apy(i,ap,av-1)*t3
        %finish %else sep=bpy(i,bp,bv)
        %if sep>maxsep %then maxsep=sep
      %repeat
    %repeat
  %repeat
%end
%end
%endoffile