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

for mon=1:mbox,
   p=0;
   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 (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;
end;

for mon=1:melipsoid,
   p=0;
   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 (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;
end;

for mon=1:mtcone,
   p=0;
   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 (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;
end;

for mon=1:mttrai,
   p=0;
   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 (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;
end;


m=max(gridd);
p=0;

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))]);
	    a1(p)=1;
	    a2(p)=1;
	    a3(p)=1;
         end;
	 gridd(p,2)=gridd(p,2)+0.0001;
      end;
   end;
end;
grid on;