function out=rjfindlist(X,Y,Z,C)
%rjfindlists(X,Y,Z,C)
%
%A help file used to convert a matlab surface to a PovRay mesh2 object.
%It is called from rjmat2pov.
%
%The input is a surface specified by the matrixes X,Y and Z, and the colormatrix
%C (with pointers to the current colormap).
%The output is a set of the form {vertex,index,cindex,normal}:
%
%vertex: A matrix (list) of the surface vertices
%index:  A matrix (list) of the three indices for the vertices of every
%        triangle and a colorindex
%normal: A matrix (list) of the normals of the vertexes 

s1=size(X,1);
s2=size(X,2);

%Making a list of the vertex coordinates of the surface

n_vertices=s1*s2;
vertex=zeros(n_vertices,3);      %Initializing the vertex matrix
k=0;
for i=1:s1
    for j=1:s2
	k=k+1;
	vertex(k,:)=[X(i,j) Y(i,j) Z(i,j)]; %PovRay is supposedly lefthanded -
                                            %but doesn't seem so
                                            %practically - so I add no
                                            %minuses here ... 
     end
end

%Making a list with the vertex indices of the triangles of the surface
%and list with the color index of the triangle

n_triangles=(size(X,1)-1)*(size(X,2)-1)*2;
index=zeros(n_triangles,4);
%cindex=zeros(n_triangles);
k=0;
for i=1:s1-1
    for j=1:s2-1
	k=k+1;
	topleftindex=(i-1)*s2+j-1;   %N.B Povray has zero as starting index
	%index(k,:)=[topleftindex,topleftindex+s2,topleftindex+1,C(i,j)-1];k=k+1;
	%index(k,:)=[topleftindex+s2,topleftindex+s2+1, topleftindex+1,C(i,j)-1];
	index(k,:)=[topleftindex,topleftindex+s2,topleftindex+s2+1,C(i,j)-1];k=k+1;
	index(k,:)=[topleftindex,topleftindex+1, topleftindex+s2+1,C(i,j)-1];
    end
end


%Making a list of the normal vectors of the vertices. These are needed to use smooth
%triangles in PovRay -- corresponding to matlabs lighting phong)

normal=zeros(n_vertices,3);
k=0;
disp('Calculating normals')
%for i=1:s1
%    for j=1:s2
%	k=k+1;
%	normal(k,:)=rjfindnormal(i,j,X,Y,Z);
%    end
%end 
%
rjn=rjsurfacenormal(X,Y,Z);
for i=1:s1
    for j=1:s2
	k=k+1;
	temp=permute(rjn(i,j,:),[1 3 2]);
	tn=norm(temp);
	if tn~=0;
	    normal(k,:)=temp/tn;
	else
	    normal(k,:)=[1 0 0];   %Just a dummy for mesh2
	end
    end
end
disp('Done with normals')

out={vertex,index,normal};











