gr1.m

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

for mon=1:mtcone,
   h = waitbar(0,['First grid, Gridding Cone 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 (ciitc(xtcs(mon),ytcs(mon),ztcs(mon),xtce(mon),ytce(mon),ztce(mon),rtcs(mon),rtce(mon),ptcs(mon),ptce(mon),x0,y0,z0)==1), 
               gridd(p,1)=gridd(p,1)+eptc1(mon);
               gridd(p,2)=gridd(p,2)+eptc2(mon);
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Second grid, Gridding Cone 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 (ciitc(xtcs(mon),ytcs(mon),ztcs(mon),xtce(mon),ytce(mon),ztce(mon),rtcs(mon),rtce(mon),ptcs(mon),ptce(mon),x0,y0,z0)==1), 
	       a4(p)=a4(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Third grid, Gridding Cone 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 (ciitc(xtcs(mon),ytcs(mon),ztcs(mon),xtce(mon),ytce(mon),ztce(mon),rtcs(mon),rtce(mon),ptcs(mon),ptce(mon),x0,y0,z0)==1), 
	       a4(p)=a4(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Fourth grid, Gridding Cone 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 (ciitc(xtcs(mon),ytcs(mon),ztcs(mon),xtce(mon),ytce(mon),ztce(mon),rtcs(mon),rtce(mon),ptcs(mon),ptce(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; 
		  a4(p)=1;
               else
		  a1(p)=0; 
 	       end; 
            end; 
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Fith grid, Gridding Cone 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 (ciitc(xtcs(mon),ytcs(mon),ztcs(mon),xtce(mon),ytce(mon),ztce(mon),rtcs(mon),rtce(mon),ptcs(mon),ptce(mon),x0,y0,z0)==1), 
	       a5(p)=a5(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Sixth grid, Gridding Cone 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 (ciitc(xtcs(mon),ytcs(mon),ztcs(mon),xtce(mon),ytce(mon),ztce(mon),rtcs(mon),rtce(mon),ptcs(mon),ptce(mon),x0,y0,z0)==1), 
	       a5(p)=a5(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Seventh grid, Gridding Cone 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 (ciitc(xtcs(mon),ytcs(mon),ztcs(mon),xtce(mon),ytce(mon),ztce(mon),rtcs(mon),rtce(mon),ptcs(mon),ptce(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 Cone 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 (ciitc(xtcs(mon),ytcs(mon),ztcs(mon),xtce(mon),ytce(mon),ztce(mon),rtcs(mon),rtce(mon),ptcs(mon),ptce(mon),x0,y0,z0)==1), 
	       a6(p)=a6(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Ninth grid, Gridding Cone 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 (ciitc(xtcs(mon),ytcs(mon),ztcs(mon),xtce(mon),ytce(mon),ztce(mon),rtcs(mon),rtce(mon),ptcs(mon),ptce(mon),x0,y0,z0)==1), 
	       a6(p)=a6(p)+1;
            end;
         end;
      end;
   end;
   close(h);
   h = waitbar(0,['Tenth grid, Gridding Cone 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 (ciitc(xtcs(mon),ytcs(mon),ztcs(mon),xtce(mon),ytce(mon),ztce(mon),rtcs(mon),rtce(mon),ptcs(mon),ptce(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;


%%%%%%%%%%%%%%%%%%%%%  Trancated Triangular %%%%%%%%%


for mon=1:mttrai,
   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;
            if (ciitt(xtts(mon,:),ytts(mon,:),ztts(mon,:),xtte(mon,:),ytte(mon,:),ztte(mon,:),x0,y0,z0)==1), 
               gridd(p,1)=gridd(p,1)+eptt1(mon);
               gridd(p,2)=gridd(p,2)+eptt2(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=p+1;
            if (ciitt(xtts(mon,:),ytts(mon,:),ztts(mon,:),xtte(mon,:),ytte(mon,:),ztte(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=p+1;
            if (ciitt(xtts(mon,:),ytts(mon,:),ztts(mon,:),xtte(mon,:),ytte(mon,:),ztte(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=p+1;
            if (ciitt(xtts(mon,:),ytts(mon,:),ztts(mon,:),xtte(mon,:),ytte(mon,:),ztte(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;
   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=p+1;
            if (ciitt(xtts(mon,:),ytts(mon,:),ztts(mon,:),xtte(mon,:),ytte(mon,:),ztte(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=p+1;
            if (ciitt(xtts(mon,:),ytts(mon,:),ztts(mon,:),xtte(mon,:),ytte(mon,:),ztte(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=p+1;
            if (ciitt(xtts(mon,:),ytts(mon,:),ztts(mon,:),xtte(mon,:),ytte(mon,:),ztte(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;
   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=p+1;
            if (ciitt(xtts(mon,:),ytts(mon,:),ztts(mon,:),xtte(mon,:),ytte(mon,:),ztte(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=p+1;
            if (ciitt(xtts(mon,:),ytts(mon,:),ztts(mon,:),xtte(mon,:),ytte(mon,:),ztte(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=p+1;
            if (ciitt(xtts(mon,:),ytts(mon,:),ztts(mon,:),xtte(mon,:),ytte(mon,:),ztte(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;
end;


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',[abs(sqrt(gridd(p,1)^2+gridd(p,2)^2))/abs(sqrt(m(1)^2+m(2)^2)) 0.5 1 ]);
            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;
rotate3d on;
plot3d(XX,30,20)
colormap('gray');