gr1.m

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

% Grid

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);

for mon=1:mbox,
   h = waitbar(0,['First grid, Gridding box number:',num2str(mon)]);
   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 (ciib(x(mon,:),y(mon,:),z(mon,:),x0,y0,z0)==1), 
               gridd(p,1)=gridd(p,1)+ep1(mon); 
               gridd(p,2)=gridd(p,2)+ep2(mon);
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Second grid, Gridding box number:',num2str(mon)]);

   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 (ciib(x(mon,:),y(mon,:),z(mon,:),x0,y0,z0)==1), 
	       a4(p)=a4(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Third grid, Gridding box number:',num2str(mon)]);
   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 (ciib(x(mon,:),y(mon,:),z(mon,:),x0,y0,z0)==1), 
	       a4(p)=a4(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Fourth grid, Gridding box number:',num2str(mon)]);
   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 (ciib(x(mon,:),y(mon,:),z(mon,:),x0,y0,z0)==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,['Thith grid, Gridding box number:',num2str(mon)]);

   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 (ciib(x(mon,:),y(mon,:),z(mon,:),x0,y0,z0)==1), 
	       a5(p)=a5(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Sixth grid, Gridding box number:',num2str(mon)]);
   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 (ciib(x(mon,:),y(mon,:),z(mon,:),x0,y0,z0)==1), 
	       a5(p)=a5(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Seventh grid, Gridding box number:',num2str(mon)]);
   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 (ciib(x(mon,:),y(mon,:),z(mon,:),x0,y0,z0)==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 box number:',num2str(mon)]);
   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 (ciib(x(mon,:),y(mon,:),z(mon,:),x0,y0,z0)==1), 
	       a6(p)=a6(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Ninth grid, Gridding box number:',num2str(mon)]);
   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 (ciib(x(mon,:),y(mon,:),z(mon,:),x0,y0,z0)==1), 
	       a6(p)=a6(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Tenth grid, Gridding box number:',num2str(mon)]);
   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 (ciib(x(mon,:),y(mon,:),z(mon,:),x0,y0,z0)==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);
end;



%%%%%%%%% Elipsoid %%%%%%%%%%

for mon=1:melipsoid,
   p=0;
   for k=1: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;
	    p=i+nx*((j-1)+ny*(k-1));
            if (ciie(xec(mon),yec(mon),zec(mon),ae(mon),be(mon),ce(mon),tes(mon),tee(mon),pes(mon),pee(mon),x0,y0,z0)==1), 
               gridd(p,1)=gridd(p,1)+epel1(mon); 
               gridd(p,2)=gridd(p,2)+epel2(mon);
            end;
         end;
      end;
   end;
   p=0;
   a4=zeros(nx*ny*nz,1);
   for k=1: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=i+nx*((j-1)+ny*(k-1));
            p=p+1;
            if (ciie(xec(mon),yec(mon),zec(mon),ae(mon),be(mon),ce(mon),tes(mon),tee(mon),pes(mon),pee(mon),x0,y0,z0)==1), 
	       a4(p)=a4(p)+1;
            end;
         end;
      end;
   end;
   p=0;
   for k=1: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=i+nx*((j-1)+ny*(k-1));
            p=p+1;
            if (ciie(xec(mon),yec(mon),zec(mon),ae(mon),be(mon),ce(mon),tes(mon),tee(mon),pes(mon),pee(mon),x0,y0,z0)==1), 
	       a4(p)=a4(p)+1;
            end;
         end;
      end;
   end;
   p=0;
   for k=1:nz,
      z0=dz*(k-0.5);
      for j=1:ny,
         y0=dy*(j-0.5);
         for i=1:nx-1,
            x0=dx*(i);
%	    p=i+nx*((j-1)+ny*(k-1));
            p=p+1;
	    p1=i+nx*((j-1)+ny*(k-1));
            if (ciie(xec(mon),yec(mon),zec(mon),ae(mon),be(mon),ce(mon),tes(mon),tee(mon),pes(mon),pee(mon),x0,y0,z0)==1), 
	       a4(p)=a4(p)+1;
            end;
	    if (a4(p)==3), 
	       if (max(gridd(p1,:))>1.01),  
		  a1(p)=1; 
               else, 
		  a1(p)=0; 
 	       end; 
            end; 
         end;
      end;
   end;
   p=0;
   a5=zeros(nx*ny*nz,1);
   for k=1: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=i+nx*((j-1)+(k-1)*ny);
             p=p+1;
            if (ciie(xec(mon),yec(mon),zec(mon),ae(mon),be(mon),ce(mon),tes(mon),tee(mon),pes(mon),pee(mon),x0,y0,z0)==1), 
	       a5(p)=a5(p)+1;
            end;
         end;
      end;
   end;
   p=0;
   for k=1: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=i+nx*((j-1)+(k-1)*ny);
            p=p+1;
            if (ciie(xec(mon),yec(mon),zec(mon),ae(mon),be(mon),ce(mon),tes(mon),tee(mon),pes(mon),pee(mon),x0,y0,z0)==1), 
	       a5(p)=a5(p)+1;
            end;
         end;
      end;
   end;
   p=0;
   for k=1:nz,
      z0=dz*(k-0.5);
      for j=1:ny-1,
         y0=dy*(j);
         for i=1:nx,
            x0=dx*(i-0.5);
%            p=i+nx*((j-1)+(k-1)*ny);
            p=p+1;
	    p1=i+nx*((j-1)+ny*(k-1));
            if (ciie(xec(mon),yec(mon),zec(mon),ae(mon),be(mon),ce(mon),tes(mon),tee(mon),pes(mon),pee(mon),x0,y0,z0)==1), 
	       a5(p)=a5(p)+1;
            end;
	    if (a5(p)==3), 
	       if (max(gridd(p1,:))>1.01),  
		  a2(p)=1; 
               else, 
		  a2(p)=0; 
 	       end; 
            end; 
         end;
      end;
   end;
   p=0;
   a6=zeros(nx*ny*nz,1);
   for k=1: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=i+nx*((j-1)+(k-1)*ny);
            p=p+1;
            if (ciie(xec(mon),yec(mon),zec(mon),ae(mon),be(mon),ce(mon),tes(mon),tee(mon),pes(mon),pee(mon),x0,y0,z0)==1), 
	       a6(p)=a6(p)+1;
            end;
         end;
      end;
   end;
   p=0;
   for k=1: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=i+nx*((j-1)+(k-1)*ny);
            p=p+1;
            if (ciie(xec(mon),yec(mon),zec(mon),ae(mon),be(mon),ce(mon),tes(mon),tee(mon),pes(mon),pee(mon),x0,y0,z0)==1), 
	       a6(p)=a6(p)+1;
            end;
         end;
      end;
   end;
   p=0;
   for k=1:nz-1,
      z0=dz*(k);
      for j=1:ny,
         y0=dy*(j-0.5);
         for i=1:nx,
            x0=dx*(i-0.5);
%            p=i+nx*((j-1)+(k-1)*ny);
            p=p+1;
	    p1=i+nx*((j-1)+ny*(k-1));
            if (ciie(xec(mon),yec(mon),zec(mon),ae(mon),be(mon),ce(mon),tes(mon),tee(mon),pes(mon),pee(mon),x0,y0,z0)==1), 
	       a6(p)=a6(p)+1;
            end;
	    if (a6(p)==3), 
	       if (max(gridd(p1,:))>1.01),  
		  a3(p)=1; 
               else, 
		  a3(p)=0; 
 	       end; 
            end; 
         end;
      end;
   end;
end;


%%%%%%%%%%%%%% Trancated Cone  %%%%%%%%%%%%