gr2.m

This software gave volume MoM solution by the CG-FFT method

% Grid

[FILENAME, PATHNAME] = uigetfile('*', 'Enter File Which contains the points on the surface');
if (FILENAME==0), break ; end;

s=['load ',PATHNAME,FILENAME,' -ascii'];
eval(s);
s=['pat=',FILENAME,';'];
eval(s);
pat(:,1)=pat(:,1)-min(pat(:,1));
pat(:,2)=pat(:,2)-min(pat(:,2));
pat(:,3)=pat(:,3)-min(pat(:,3));
%pat=pat(1:100:length(pat),:);


% pat should be a matrix with 3 coloms for x, y and z on the arbitrary surface

disp('start')
figure(2);
clf;
set (gcf,'Position',[0.5 0.1 0.35 0.4],'Name','Gridded Data View','NumberTitle','off','Units','normalized');
set(gcf,'position',[0.5 0.1 0.35 0.4],'Units','normalized');

plot3(nx*dx,0,0);
hold;
plot3(0,ny*dy,0);
plot3(0,0,nz*dz);

gridd=zeros(nx*ny*nz,2);
a1=zeros(nx*ny*nz,1);
a2=zeros(nx*ny*nz,1);
a3=zeros(nx*ny*nz,1);
a4=zeros(nx*ny*nz,1);
a5=zeros(nx*ny*nz,1);
a6=zeros(nx*ny*nz,1);

h = waitbar(0,['First grid, Gridding Arbitrary Surface Object:']);
   p=0;
   for k=1:nz,
      waitbar(k/nz);
      z0=dz*(k-0.5);
      for j=1:ny,
         y0=dy*(j-0.5);
         for i=1:nx,
            x0=dx*(i-0.5);
            p=p+1;
            if (ciias(pat,x0,y0,z0,dx,dy,dz)==1), 
               gridd(p,1)=gridd(p,1)+1.;  % Epsilon is for perfect conductor
               gridd(p,2)=gridd(p,2)+1000.; % Epsilon is for perfect conductor 
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Second grid, Gridding Arbitrary Surface Object:']);

   p=0;
   a4=zeros(nx*ny*nz,1);
   for k=1:nz,
      waitbar(k/nz);
      z0=dz*(k-0.5);
      for j=1:ny,
         y0=dy*(j-0.5);
         for i=1:nx-1,
            x0=dx*(i+0.51);
            p=p+1;
            if (ciias(pat,x0,y0,z0,dx,dy,dz)==1), 
      	       a4(p)=a4(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Third grid, Gridding Arbitrary Surface Object:']);
   p=0;
   for k=1:nz,
      waitbar(k/nz);
      z0=dz*(k-0.5);
      for j=1:ny,
         y0=dy*(j-0.5);
         for i=1:nx-1,
            x0=dx*(i-0.51);
            p=p+1;
            if (ciias(pat,x0,y0,z0,dx,dy,dz)==1), 
	             a4(p)=a4(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Fourth grid, Gridding Arbitrary Surface Object:']);
   p=0;
   for k=1:nz,
      waitbar(k/nz);
      z0=dz*(k-0.5);
      for j=1:ny,
         y0=dy*(j-0.5);
         for i=1:nx-1,
            x0=dx*(i);
            p=p+1;
            if (ciias(pat,x0,y0,z0,dx,dy,dz)==1), 
      	       a4(p)=a4(p)+1;
            end;
	         if (a4(p)==3), 
	             p1=i+nx*((j-1)+ny*(k-1));
	             if (max(gridd(p1,:))>1.01),  
            		  a1(p)=1; 
                else, 
		              a1(p)=0; 
      	       end; 
           end; 
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Fith grid, Gridding Arbitrary Surface Object:']);

   p=0;
   a5=zeros(nx*ny*nz,1);
   for k=1:nz,
      waitbar(k/nz);
      z0=dz*(k-0.5);
      for j=1:ny-1,
         y0=dy*(j+0.51);
         for i=1:nx,
            x0=dx*(i-0.5);
            p=p+1;
            if (ciias(pat,x0,y0,z0,dx,dy,dz)==1), 
      	       a5(p)=a5(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Sixth grid, Gridding Arbitrary Surface Object:']);
   p=0;
   for k=1:nz,
      waitbar(k/nz);
      z0=dz*(k-0.5);
      for j=1:ny-1,
         y0=dy*(j-0.51);
         for i=1:nx,
            x0=dx*(i-0.5);
            p=p+1;
            if (ciias(pat,x0,y0,z0,dx,dy,dz)==1), 
               a5(p)=a5(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Seventh grid, Gridding Arbitrary Surface Object:']);
   p=0;
   for k=1:nz,
      waitbar(k/nz);
      z0=dz*(k-0.5);
      for j=1:ny-1,
         y0=dy*(j);
         for i=1:nx,
            x0=dx*(i-0.5);
            p=p+1;
            if (ciias(pat,x0,y0,z0,dx,dy,dz)==1), 
      	       a5(p)=a5(p)+1;
            end;
     	      if (a5(p)==3), 
               p1=i+nx*((j-1)+ny*(k-1));
	            if (max(gridd(p1,:))>1.01),  
		            a2(p)=1; 
               else, 
       		      a2(p)=0; 
 	            end; 
            end; 
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Eighth grid, Gridding Arbitrary Surface Object:']);
   p=0;
   a6=zeros(nx*ny*nz,1);
   for k=1:nz-1,
      waitbar(k/(nz-1));
      z0=dz*(k+0.51);
      for j=1:ny,
         y0=dy*(j-0.5);
         for i=1:nx,
            x0=dx*(i-0.5);
            p=p+1;
            if (ciias(pat,x0,y0,z0,dx,dy,dz)==1), 
 	            a6(p)=a6(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Ninth grid, Gridding Arbitrary Surface Object:']);
   p=0;
   for k=1:nz-1,
      waitbar(k/(nz-1));
      z0=dz*(k-0.51);
      for j=1:ny,
         y0=dy*(j-0.5);
         for i=1:nx,
            x0=dx*(i-0.5);
            p=p+1;
            if (ciias(pat,x0,y0,z0,dx,dy,dz)==1), 
      	       a6(p)=a6(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Tenth grid, Gridding Arbitrary Surface Object:']);
   p=0;
   for k=1:nz-1,
      waitbar(k/(nz-1));
      z0=dz*(k);
      for j=1:ny,
         y0=dy*(j-0.5);
         for i=1:nx,
            x0=dx*(i-0.5);
            p=p+1;
            if (ciias(pat,x0,y0,z0,dx,dy,dz)==1), 
      	       a6(p)=a6(p)+1;
            end;
	         if (a6(p)==3), 
      	       p1=i+nx*((j-1)+ny*(k-1));
	             if (max(gridd(p1,:))>1.01),  
            		 a3(p)=1; 
                else, 
		             a3(p)=0; 
 	             end; 
            end; 
         end;
      end;
   end;
   close(h);


m=max(gridd);
p=0;

XX=zeros(nx,ny,nz);
for k=1:nz,
   z0=dz*(k-1);
   for j=1:ny,
      y0=dy*(j-1);
      for i=1:nx,
         x0=dx*(i-1);
         p=p+1;
         if (max(gridd(p,:)>1.01)),
            plot3(x0,y0,z0,'Color',[0 1 abs(sqrt(gridd(p,1)^2+gridd(p,2)^2))/abs(sqrt(m(1)^2+m(2)^2))]);
            XX(i,j,k)=1;
%	    a1(p)=1;
%	    a2(p)=1;
%	    a3(p)=1;
         end;
	 gridd(p,2)=gridd(p,2)+0.0001;
      end;
   end;
end;
grid on;

%figure
%plot3d(XX,30,20)
%colormap('gray');