%========================================================================
%FMM demo to COMPUTE the interaction lists for a box (node) on a level
%========================================================================
%Usage: out=intera(node, level, io)    %e.g intera(3,2,1)    io = 0, 1
%   or  out=intera(node, level);       %e.g intera(3,2) % node=3
%   or  out=intera(node);              %e.g intera(3)   % for level=2
%
% Stuart C. Hawkins and Ke Chen (University of Liverpool, UK)
%
% SEE also index2.m
 
% Main internal notation ___________________________________________________
%  F0        : List of starting global index addresses of each level
%  ANCESTORS : List of positions of the node's ancestors at their level
% index ---- : Convert a local index (i,j) at level j to global index k
% vert ----- : Find the local index j for (i,j) from a global k
% horiz ---- : Find the local index i for (i,j) from a global k
% neighbours : Find the neighbour's index list for any node on a level
%            : (using near_neighbours)
% draw_intera: Plot near neighbours and interaction lists
% parent     : Find the parent's index id, given a child's index (i,j)
% children   : Find a parent's 4 children, given a children index list k
% Figure     : Automatic subplotting
% --------------------------------------------------------------------------

function out=INTERA_demo(node,levs,iout)
 if nargin<3, iout=0;       end  % plotting or not
 if nargin<2, levs=2;       end
 if nargin<1, help intera, return, end, hold off
 fprintf('-----intera.m------Finding near neighours for node=%d, ',node)
 fprintf('\n\t\ton a given finest level = %d\n',levs) 
  if node<1|node>4^levs, help intera, out=-1;
    fprintf('\tInput "node" out of range  /  1:%d\n',4^levs); return, end 
% get starting index of level l=levs' starting 
 F0=[];F0(levs)=index(1,1,levs); % position = 4^0+4^1+...+4^(l-1)
 out.node=node;

% update target --- node in "global" index (shift = f0-1)
 node=node-1+F0(levs);

% get neighbours of node across levels
        j=vert(node,levs);
        i=horiz(node,j,levs);  % local (i,j) for global node
    %-----------------% Find ancestors' id-------------------
        ANCESTORS=node; % the global index
      for lev=levs-1:-1:1
          F0(lev)=index(1,1,lev); % position shifted globally
        [ip jp lp]=parent(i,j,lev+1);
        ANCESTORS=[index(ip,jp,lp) ANCESTORS];
         i=ip; j=jp;
      end % lev
      for ll=1:levs
        Target=ANCESTORS(ll);
        loc=neighbours(Target,ll); %-generate out.near1 etc
        t=sprintf('out.near%d=loc-F0(ll)+1;',ll);eval(t)
      end % ll
for ll=levs-1:-1:1        %---------------------------------%
    t=sprintf('in_list=children(out.near%d+F0(ll)-1,ll,F0(ll));',ll);
        eval(t) %-----Find the interaction list (Step 1) ----- 
    t=sprintf('IL=setdiff(in_list,out.near%d);',ll+1);
        eval(t) %-----Find the interaction list (Step 2) -----
    t=sprintf('out.IL%d=IL;',ll+1);eval(t)
end %ll
  disp(' Main task of intera.m complete ...')
  out.lev=levs;
%---------------- Interaction list ---------------------------
if iout>0
    Figure(1)
    w=2; plot(0,0,'.');axis([0 w 0 w]); my_grid(1), hold on
       t=title(['Ancestor and coarest level ' num2str(1)]); 
            set(t,'fontsize',15); set(t,'color','r'),axis off
        target = ANCESTORS(1);  %coarest level
        draw_intera(out.near1,target,1,F0);
end

if iout>0 %------------------------------------------------%
for ll=levs-1:-1:1  
  Figure(ll+1)
  w=2^(ll+1); plot(0,0,'.');axis([0 w 0 w]); my_grid(1)
  if ll==levs-1
     t=title(['Box ' num2str(out.node) ' on Finest Level ' num2str(ll+1)]); 
          set(t,'fontsize',15); set(t,'color','r')
  else
     t=title(['Parent and Far-field Level ' num2str(ll+1)]); 
          set(t,'fontsize',15); set(t,'color','b')
  end, hold on
  t=sprintf('IL=out.IL%d;',ll+1);eval(t)
        jfine= vert(IL+F0(ll+1)-1,ll+1);
        ifine=horiz(IL+F0(ll+1)-1,jfine,ll+1);
    for k=1:length(ifine)
        t=text(ifine(k)-0.7,jfine(k)-0.7, num2str(IL(k)));
        set(t,'fontsize',30/(ll+1)), set(t,'color','r')
    end % k
        target = ANCESTORS(ll+1);  %fine level
    t=sprintf('draw_intera(out.near%d,target,ll+1,F0);',ll+1);
    eval(t)
end % ll
end % iout %-------------------------------------------------%

return %*************************************************

%========================================================
function [near]=neighbours(f,lev) 
% get coordinates of f (f shifted by f0-1)
  j=vert(f,lev); i=horiz(f,j,lev); 
% get (i,j) coords of neighbours of f
  [ii,jj,ll]=near_neighbours(i,j,lev); 
% get indexes of the neighbours of f from their (i,j) coords
  near=index(ii,jj,ll);

%========================================================
function [ii,jj,ll]=near_neighbours(i,j,l) 
% initialise arrays -- neighbour (i,j) list + lev list
  ii=[]; jj=[]; ll=[]; 
% loop through coordinates----Search in +/-1 Box
  for cj=max(1,j-1):min(2^l,j+1) 
      for ci=max(1,i-1):min(2^l,i+1) 
          if( cj~=j | ci~=i )
              ll=[ll;l]; ii=[ii;ci]; jj=[jj;cj];
          end 
      end 
  end

%===================================================
function l=level(f)  %linking f=partial sums to lev
  l=floor( log(3*(f-1)+1)/log(4) ); 
% Again the storage system = 4^0+4^1+4^2+...+4^LEV 

%===================================================
function j=vert(f,l)
  j=floor( (f-1-(4.^l-1)/3)./(2.^l)+1 );

%===================================================
function i=horiz(f,j,l)
  i=f-2.^l.*(j-1)-(4.^l-1)/3;

%===================================================
function f=index(i,j,l); % (i,j) to single f (shifted)
  f=i+2.^l.*(j-1)+(4.^l-1)/3;

%===================================================
function [ii,jj,ll]=parent(i,j,l)
ll=l-1; ii=ceil(i/2); jj=ceil(j/2);

%===================================================
function listo=children(k,lev,f0)
   k = k + 1 - f0;  % shifted back to [1, 4^lev]
  two=2^lev; twop=two*2;
   i = rem(k-1,two)+1; % horizontal index
  k1 = 4*k - 2*i; % vector
  listo = [ k1-1; k1+twop-1; k1+twop; k1]; 
%=====================================================
%                    loc  globe  l   shifts
function draw_intera(near,target,lev,F0);
          cc=[1.0,0.5,0.0]; % orange
        for k=1:length(near) % get coordinates of k
            j=vert(near(k)+F0(lev)-1,lev);
            i=horiz(near(k)+F0(lev)-1,j,lev); % draw patch
            patch([i-1 i-1 i i],[j j-1 j-1 j],cc); 
        end 
        %-------------% draw target------------------
          j=vert(target,lev); i=horiz(target,j,lev); 
        t=text(i-0.7,j-0.7,num2str(target-F0(lev)+1));
        set(t,'fontsize',30/lev), set(t,'color','b')
        t=plot(i-0.4,j-0.5,'r+'); set(t,'markersize',25/lev)
        axis off
%=====================================================
function Figure(in)
   m=4; n=1;  % for m x n subplotting
   total = m*n;
   fig = floor( (in-1)/total );     % sub = subplot no
   sub = in - total*fig; fig=fig+1; % fig = figure no
   figure(fig); 
   subplot(m,n,sub); hold on 
function my_grid(n)
  if nargin<1, n=1; end
  hold on
a=axis; x=a(1)/2;
    plot( [x x], a(3:4), 'k-')
  while x<a(2)
    x=x+n;
    plot( [x x], a(3:4), 'k-')
  end
y=a(3)/2;
    plot( a(1:2), [y y], 'k-')
  while y<a(4)
    y=y+n;
    plot( a(1:2), [y y], 'k-')
  end
hold off
