function [p,t] = volumemesh3(Size,Depth,CellSize,NCell,NDepth)
%   VOLUMESH3 Creates a symmetric volume rectangular mesh of 
%   intercepting strips 
%
%   Arguments:
%   Size of the square mesh     (m)             Size
%   Depth of the mesh           (m)             Depth
%   Size of the rectangular cell(m)             CellSize
%   Number of discretization rectangles 
%   per cell length                             NCell
%   Number of discretization rectangles 
%   per cell depth                              NDepth
%
%   Example: volumemesh3(1,0.05,0.2,3,3)
%
%   Copyright 2002 AEMM. Revision 2002/05/06 Appendix C


%   Identify strip parameters
    W   =Depth;
    L   =Size;
    Ny  =Size/CellSize*NCell;
    Nx  =NDepth;
    N   =Size/CellSize;
    
%   Clone strips along the y axis 
[p,t] = strip(W,L,Nx,Ny);
p(2,:)=p(2,:)-Size/2+CellSize; %first strip
tbase   =t;
pbase   =p;
for k=1:N-2
    pbase1(1,:)=pbase(1,:);
    pbase1(2,:)=pbase(2,:)+k*CellSize;
    pbase1(3,:)=pbase(3,:);
    p=[p pbase1];    
    t=[t tbase+k*length(pbase)]; 
end
%   Clone strips along the x axis 
[pbase,tbase] = strip(W,L,Nx,Ny);
pbase(2,:)  =pbase(1,:);
pbase(1,:)  =-Size/2+CellSize; %first strip
for k=1:N-1
    pbase1(1,:)=pbase(1,:)+(k-1)*CellSize;
    pbase1(2,:)=pbase(2,:);
    pbase1(3,:)=pbase(3,:);
    p=[p pbase1];    
    t=[t tbase+(N-2+k)*length(pbase)]; 
end

%   Introduce indexing (used in rwg1)
t(4,:)=1; t(4,1:length(t)/2)=0; 

epsilon=5e-6;   %Important  
%   Eliminate double points
P(1:3,1)=p(1:3,1);
for L=2:length(p)
   [M,N]=size(P);
   Point=p(1:3,L);
   Index=find(sum(abs(P-repmat(Point,1,N)))<epsilon);
   if (isempty(Index))
       P=[P Point];
   end  
end

%   Reenumerate triangle vertexes
[M,N]=size(P);
for L=1:length(t)
    n=t(:,L);
    for i=1:3
        Index=find(sum(abs(P-repmat(p(1:3,n(i)),1,N)))<epsilon);
        if(~isempty(Index))
            t(i,L)=Index;
        end
    end
end
clear p; p=P;

save mesh p t;
viewer mesh