% This code tries to fit as much of a quadratic surface as possible to 
% a collection of points.

clear;

format rat; % sets rational format on

for N=4:11,

  for i=2:N+1,
    t = (i-2)*2*pi/N;
    xi(i,1) = 0.5*cos(t);
    xi(i,2) = 0.5*sin(t);
  end;
%  xi = xi;


  %p=[0 1 0 1 2 0 1 2 2;0 0 1 1 0 2 2 1 2];

  %p=[0 1 1 2 1 0 2 2 0;0 1 2 1 0 1 2 0 2];

  % powers of xis
  p=[
    0 1 0 1 2 0 2 1 2 3 0 3 1 3 2 3
    0 0 1 1 0 2 1 2 2 0 3 1 3 2 3 3
    ];

  for j=1:N+1,
    xi1 = xi(j,1);
    xi2 = xi(j,2);
    for k=1:16,
      A(j,k) = xi1^p(1,k)*xi2^p(2,k);
    end;

    if j ~= 2
      b(j,1) = 0;
    else
      b(j,1) = 1;
    end;
  end;
  % calculate coefficients
  c = A\b;

% visualise
%n=2;
%xi1 = 0.0:0.05:1;
%xi2 = xi1;
%[xxi1,xxi2]=meshgrid(xi1,xi2);

%[n1,n2]=size(xxi1);
%phi = zeros(n1,n2);

%for j=1:min(N+1,9),
%for j=1:9,
%  phi = phi + c(j).*xxi1.^p(1,j).*xxi2.^p(2,j);
%end;

%theta = (n-2)*2*pi/N;
%psi1 =  (cos(theta).*(xxi1-0.5)+sin(theta).*(xxi2-0.5));
%psi2 = (-sin(theta).*(xxi1-0.5)+cos(theta).*(xxi2-0.5));

%for j=1:16,
%  phi = phi + c(j).*psi1.^p(1,j).*psi2.^p(2,j);
%end;


%clf;
%grid on;
%hold on;
%ptphi = zeros(N+1,1);
%ptphi(n) = 1;
%plot3(xi(:,1)+0.5,xi(:,2)+0.5,ptphi,'ok','linewidth',4);
%view(3);
%surf(xi1,xi2,phi);
%hold off;

  fprintf('Number of perimeter grid points: %d\n',N);
  for k=1:16,
      fprintf('%5.2fxi1^(%d)xi2^(%d) +',c(k),p(1,k),p(2,k));
  end;
  fprintf('\n');

  display(c);
%  for k=1:16,
%      display(c(k));
%      fprintf('%20.10f\n',c(k));
%  end;

end; % end of N loop

% apex basis function
%for j=1:N+1,
%  xi1 = xi(j,1)+0.5;
%  xi2 = xi(j,2)+0.5;
%  for k=1:9,
%    A(j,k) = xi1^p(1,k)*xi2^p(2,k);
%  end;

%  if j ~= 1
%    b(j,1) = 0;
%  else
%    b(j,1) = 1;
%  end;
%end;
%cApex = A\b;
