function [M,T] = A366311(s)
% A366311: given a a side length "s", gives 
%  M: The 2*s-1 x 2*s-1 Square Spiral representation of A366311
%  T: Table of the first (2*s-1)^2 terms of A366311
s=s+3
mx=(2*(s-3)-1)^2
r=s;
c=r;
M=zeros(2*r-1,2*c-1); %initialize matrix
M(s,s)=1; %initialize central point
s=s-3;
L=[ceil((1:(s-1)*4)/2),2*s-2]; %length in each direction
D=mod(1:s*4-3,4); %direction
P=zeros((2*s-1)^2-1,1);
q=1; %counter
T=[(1:(2*s-1)^2)' zeros((2*s-1)^2,1)];
T(1,2)=1;
for b=1:length(L)
    for f=1:L(b)
        P(q)=D(b);
        if P(q)==1
            c=c+1;
        elseif P(q)==2
            r=r-1;
        elseif P(q)==3
            c=c-1;
        else
            r=r+1;
        end
        S=[]; % surroundings 
        S=[M(r+2,c) M(r-2,c) M(r+2,c+2) M(r-2,c+2) M(r+2,c-2) M(r-2,c-2) M(r,c+2) M(r,c-2)];
        if M(r+1,c)>0 && M(r+1,c)==M(r+2,c) 
            S = [S M(r+3,c)];
        end
        if M(r-1,c)>0 && M(r-1,c)==M(r-2,c) 
            S = [S M(r-3,c)];
        end
        if M(r+1,c+1)>0 && M(r+1,c+1)==M(r+2,c+2) 
            S = [S M(r+3,c+3)];
        end
        if M(r-1,c+1)>0 && M(r-1,c+1)==M(r-2,c+2) 
            S = [S M(r-3,c+3)];
        end
        if M(r+1,c-1)>0 && M(r+1,c-1)==M(r+2,c-2) 
            S = [S M(r+3,c-3)];
        end
        if M(r-1,c-1)>0 && M(r-1,c-1)==M(r-2,c-2) 
            S = [S M(r-3,c-3)];
        end
        if M(r,c+1)>0 && M(r,c+1)==M(r,c+2) 
            S = [S M(r,c+3)];
        end
        if M(r,c-1)>0 && M(r,c-1)==M(r,c-2) 
            S = [S M(r,c-3)];
        end
        S=nonzeros(S);
        a=1;
        while any(S(:)==a)
            a=a+1;
        end
        M(r,c)=a;
        q=q+1;
        T(q,2)=a;
    end
end
M((1:3),:)=[];
M((end-2:end),:)=[];
M(:,(1:3))=[];
M(:,(end-2:end))=[];
end