% Test code to compare calculation of normals and surface areas
% by the scattering code and Matlab calculations, partly to 
% verify F117 formulation and partly for consistency check.

load duct; %sphere.dat; 
temp=duct ; %sphere;
x=temp(:,1);y=temp(:,2);z=temp(:,3);
n1=temp(:,4);n2=temp(:,5);n3=temp(:,6);areas=temp(:,7);

nx=120;nth=60;

x=reshape(x,nth,nx);
y=reshape(y,nth,nx);
z=reshape(z,nth,nx);
n1=reshape(n1,nth,nx);
n2=reshape(n2,nth,nx);
n3=reshape(n3,nth,nx);
areas=reshape(areas,nth,nx);



% surf(x,y,z,d3clip);shading interp;axis image;colorbar; 


XII=x; YII=y; ZII=z;

% Find surface normals via Matlab:
[NXX,NYY,NZZ]=surfnorm(XII,YII,ZII);   

% Now find the surface areas.  These can be obtained from a
% formula (see softsurfer.com, from search of google for "area of polygon"),
% under the assumption that the facets are PLANAR:
%   2*A = n.[(V2-V0)x(V3-V1)]  
%   ie the dot product of the normal with a cross-product
%   formed using opposite pairs of vertices.  I haven't checked this!

sizes=size(XII);
SIGMA=0*XII;
for i=1:sizes(2)-1,
    for j=1:sizes(1),
      j2=j+1;
        if j==sizes(1)
          j2=1;
        end
      V0=[XII(j,i) YII(j,i) ZII(j,i)];
      V1=[XII(j2,i) YII(j2,i) ZII(j2,i)];
      V2=[XII(j2,i+1) YII(j2,i+1) ZII(j2,i+1)];
      V3=[XII(j,i+1) YII(j,i+1) ZII(j,i+1)];
      NN=[NXX(j,i) NYY(j,i) NZZ(j,i)];
      SIGMA(j,i)=.5*dot( NN, cross( (V2-V0),(V3-V1) ) );
    end
end
     SIGMA=abs(SIGMA) ; % In case we've got the sign wrong


