toyfdtd1.m

This a translation of the ToyFDTD c code available from the web site http://www.borg.umn.edu/toyfdt

hy=zeros([nx,ny-1,nz]);
allocatedBytes = allocatedBytes+prod(size(hy));

% allocate the array of hz components:
hz=zeros([nx,ny,nz-1]);
allocatedBytes = allocatedBytes+prod(size(hz));
%MATLAB uses doubles as its maximum precision which is 8 bytes per element, (15 digits)
% note that MU and EPSILON (use 101 digits or an extended long)
allocatedBytes=allocatedBytes*8; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write some progress notes to standard output
% print out some identifying information 
disp('ToyFDTD1 version 1.03');
disp('Copyright (C) 1998,1999 Laurie E. Miller, Paul Hayes, ');
disp('Matthew O''Keefe');
disp(' ');
disp('ToyFDTD1 is free software published under the terms'); 
disp('of the GNU General Public License as published by the'); 
disp('Free Software Foundation.');  
% print out a bob command line,  
%     including the dimensions of the output files: 
disp(['bob -cmap chengGbry.cmap -s',num2str(nx+1),'x',num2str(ny+1),'x',num2str(nz),' *.bob\n']); 
disp(' '); 
% print out a viz command line: 
disp('viz ToyFDTD1c.viz'); 
disp(' '); 
% print out how much memory has been allocated:  
disp(['Dynamically allocated',num2str(allocatedBytes),' bytes']); 
disp(' '); 
% print out some simulation parameters: 
disp(['PLOT_MODULUS =', num2str(PLOT_MODULUS)]); 
disp('rectangular waveguide'); 
disp(['Stimulus =', num2str(FREQUENCY),' Hertz ']); 
disp('continuous plane wave'); 
disp(' '); 
disp('Meshing parameters:'); 
disp([num2str(nx),'x',num2str(ny),'x',num2str(nz),' cells']); 
disp(['dx=',num2str(dx),' dy=',num2str(dy),' dz=',num2str(dz),' meters']); 
disp([num2str(GUIDE_WIDTH),'x',num2str(GUIDE_HEIGHT ),'x',num2str(LENGTH_IN_WAVELENGTHS*lambda ),' meter^3 simulation region']);        
disp(' '); 
disp(['Time simulated will be ',num2str(totalSimulatedTime),' seconds,',num2str(MAXIMUM_ITERATION),' timesteps']);
disp(' '); 
disp('3D output scaling parameters:'); 
disp('Autoscaling every timestep'); 
disp(''); 
disp(' '); 
% print out some info on each timestep: 
disp('Following is the iteration number and current simulated time for each timestep/iteration of'); 
disp('the simulation.  For each timestep that 3D data is output to file, the maximum and minimum data'); 
disp('values are printed here with the maximum and minimum scaled values in parentheses.'); 
disp(' '); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% open and start writing the .viz file
%    this file will be handy to feed parameters to viz if desired
vizFilePointer = fopen('ToyFDTD1c.viz', 'w');
if(~vizFilePointer)
    error('Difficulty opening ToyFDTD1c.viz');
end
fprintf(vizFilePointer, '#Viz V1.0\n');
fprintf(vizFilePointer, 'time:%g %g\n', currentSimulatedTime, dt);
fprintf(vizFilePointer, 'color: chengGbry.cmap\n');
fprintf(vizFilePointer, '\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% main loop:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic; % start the timer
for  iteration =[1:MAXIMUM_ITERATION]
    % Output section: 
    % time in simulated seconds that the simulation has progressed:
    currentSimulatedTime = dt*iteration;  
    % print to standard output the iteration number and current simulated time:
    fprintf(stdout,'#%d, sim time= %g sec\n',iteration,currentSimulatedTime);
    % 3D data output every PLOT_MODULUS timesteps:
    %     The first time through the main loop all the data written to 
    %     file will be zeros.  If anything is nonzero, there's a bug.  :>
    if ( mod(iteration,PLOT_MODULUS) == 0)
        % bob output section
        % create the filename for this iteration, 
        %     which includes the iteration number:
        filename=sprintf('c_%06d.bob', iteration);
        % open a new data file for this iteration:
        openFilePointer = fopen(filename, 'wb');
        if(~openFilePointer)
            % if the file fails to open, print an error message to standard output:
            error(sprintf('Difficulty opening c_%06d.bob', iteration));
        end
        % find the max and min values to be output this timestep:  
        mmin = min(min(min(ez)));
        mmax = max(max(max(ez)));
        % update the tracking variables for minimum and maximum 
        %      values for the entire simulation:
        simulationMin = min(simulationMin,mmin);
        simulationMax = max(simulationMax,mmax);
        %  set nnorm to be max or min, whichever is greater in magnitude:
        %  if everything is zero, give nnorm a tiny value 
        %      to avoid division by zero:
        nnorm=norm([mmax,mmin,eps],inf); % infinity norm
        scalingValue = 127.0/nnorm;
        % write to standard output the minimum and maximum values 
        %     from this iteration and the minimum and maximum values
        %     that will be written to the bob file this iteration:
        fprintf(stdout,'\t%g(%d) < ez BoB <%g(%d)\n',mmin,ceil(127.0 + scalingValue*mmin),mmax,ceil(127.0 + scalingValue*mmax));
        % scale each ez value in the mesh to the range of integers 
        %     from zero through 254 and write them to the output 
        %     file for this iteration:
        fprintf(openFilePointer,'%c',(127 + fix(scalingValue*ez)));
        % close the output file for this iteration:
        fclose(openFilePointer);
        % write the dimensions and name of the output file for this 
        %     iteration to the viz file:
        fprintf(vizFilePointer, '%dx%dx%d%s\n', nx+1, ny+1, nz, filename);
        % write identification of the corners of the mesh and the max
        %     and min values for this iteration to the viz file:
        fprintf(vizFilePointer, 'bbox: 0.0 0.0 0.0 %g %g %g %g %g\n',dx*nx, dy*ny, dz*nz, mmin, mmax);
        if (exist('interactive'))
            if (exist('h'))
                close(h); % must close the figure otherwise we can run out of memory
            end
            if(interactive==1)
                sliceomatic(ez); 
                beep;
                disp(['Press enter to continue']);
                pause;
            else
                surfc(ez(:,:,floor(nz/2)));
                drawnow;
            end
            h=gcf;
        end
    end %output section
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
    % Compute the stimulus: a plane wave emanates from the x=0 face:
    %     The length of the guide lies in the x-direction, the width of the 
    %     guide lies in the y-direction, and the height of the guide lies
    %     in the z-direction.  So the guide is sourced by all the ez 
    %     components on the stimulus face.  
    stimulus = sin(omega*currentSimulatedTime);
    for ii=[1:2]
        ez(ii,:,:) = stimulus*ones([ny+1,nz]);
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Update the interior of the mesh:
    %    all vector components except those on the faces of the mesh.  
    %
    % Update all the H field vector components within the mesh:
    %     Since all H vectors are internal, all H values are updated here.
    %     Note that the normal H vectors on the faces of the mesh are not 
    %     computed here, and in fact were never allocated -- the normal 
    %     H components on the faces of the mesh are never used to update 
    %     any other value, so they are left out of the memory allocation  
    %     entirely.   
    kindx=[1:nz];
    % Update the hx values:
    for ii=[1:(nx-1)]
        for jj=[1:ny]
                hx(ii,jj,kindx) = hx(ii,jj,kindx)+(dtmudz*(ey(ii+1,jj,kindx+1) - ey(ii+1,jj,kindx))-dtmudy*(ez(ii+1,jj+1,kindx) - ez(ii+1,jj,kindx)));
        end
    end    
    % Update the hy values:
    for ii=[1:nx]
        for jj=[1:(ny-1)]
                hy(ii,jj,kindx) = hy(ii,jj,kindx) +(dtmudx*(ez(ii+1,jj+1,kindx) - ez(ii,jj+1,kindx)) -dtmudz*(ex(ii,jj+1,kindx+1) - ex(ii,jj+1,kindx)));
        end
    end
    % Update the hz values:
    kindx=[1:(nz-1)];
    for ii=[1:nx]          
        for jj=[1:ny]
                hz(ii,jj,kindx)=  hz(ii,jj,kindx)+ (dtmudy*(ex(ii,jj+1,kindx+1) - ex(ii,jj,kindx+1))-dtmudx*(ey(ii+1,jj,kindx+1) - ey(ii,jj,kindx+1)));
        end
    end
    % Update the E field vector components.  
    % The values on the faces of the mesh are not updated here; they 
    %      are handled by the boundary condition computation 
    %      (and stimulus computation).  
    % Update the ex values:
    kindx=[2:nz];
    for ii=[1:nx]
        for jj=[2:ny]
                ex(ii,jj,kindx) = ex(ii,jj,kindx)+ (dtepsdy*(hz(ii,jj,kindx-1) - hz(ii,jj-1,kindx-1))-dtepsdz*(hy(ii,jj-1,kindx) - hy(ii,jj-1,kindx-1)));
        end
    end      
    % Update the ey values:
    for ii=[2:nx]
        for jj=[1:ny]
                ey(ii,jj,kindx)= ey(ii,jj,kindx)+(dtepsdz*(hx(ii-1,jj,kindx) - hx(ii-1,jj,kindx-1))-dtepsdx*(hz(ii,jj,kindx-1) - hz(ii-1,jj,kindx-1)));
        end
    end
    % Update the ez values:
    kindx=[1:nz];
    for ii=[2:nx]
        for jj=[2:ny]
                ez(ii,jj,kindx) =  ez(ii,jj,kindx)+(dtepsdx*(hy(ii,jj-1,kindx) - hy(ii-1,jj-1,kindx))-dtepsdy*(hx(ii-1,jj,kindx) - hx(ii-1,jj-1,kindx)));
        end
    end  
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
    % Compute the boundary conditions:
    
    % OK, so I'm yanking your chain on this one.  The PEC condition is 
    % enforced by setting the tangential E field components on all the
    % faces of the mesh to zero every timestep (except the stimulus 
    % face).  But the lazy/efficient way out is to initialize those 
    % vectors to zero and never compute them again, which is exactly 
    % what happens in this code.  
    
end% end mainloop
xtime=toc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/ 
% Output section:  
% The output routine is repeated one last time to write out  
% the last data computed.    
% time in simulated seconds that the simulation has progressed:
currentSimulatedTime = dt*iteration;  
% print to standard output the iteration number and current simulated time:
fprintf(stdout,'#%d  sim time=%g  sec, execute time= %g sec\n', iteration, currentSimulatedTime, xtime);
% 3D data output for the last timestep: 
% create the filename for this iteration, 
%     which includes the iteration number:
filename=sprintf('c_%06d.bob', iteration);
% open a new data file for this iteration:
openFilePointer = fopen(filename, 'wb');
if(~openFilePointer)
    % if the file fails to open, print an error message to standard output:
    error(sprintf('Difficulty opening c_%06d.bob', iteration));
end
% find the max and min values to be output this timestep:  
mmin = min(min(min(ez)));
mmax = max(max(max(ez)));
% update the tracking variables for minimum and maximum 
%      values for the entire simulation:
simulationMin = min(simulationMin,mmin);
simulationMax = max(simulationMax,mmax);
%  set nnorm to be max or min, whichever is greater in magnitude:
%  if everything is zero, give nnorm a tiny value 
%      to avoid division by zero:
nnorm=norm([mmax,mmin,eps],inf); % infinity norm
scalingValue = 127.0/nnorm;
% write to standard output the minimum and maximum values 
%     from this iteration and the minimum and maximum values
%     that will be written to the bob file this iteration:
fprintf(stdout,'\t%g(%d) < ez BoB <%g(%d)',mmin,ceil(127.0 + scalingValue*mmin),mmax,ceil(127.0 + scalingValue*mmax));
% scale each ez value in the mesh to the range of integers 
%     from zero through 254 and write them to the output 
%     file for this iteration:
fprintf(openFilePointer,'%c',(127 + fix(scalingValue*ez)));
% close the output file for this iteration:
fclose(openFilePointer);
% write the dimensions and name of the output file for this 
%     iteration to the viz file:
fprintf(vizFilePointer, '%dx%dx%d%s\n', nx+1, ny+1, nz, filename);
% write identification of the corners of the mesh and the max
%     and min values for this iteration to the viz file:
fprintf(vizFilePointer, 'bbox: 0.0 0.0 0.0 %g %g %g %g %g\n',dx*nx, dy*ny, dz*nz, mmin, mmax);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
% close the viz file for this simulation:
fclose(vizFilePointer);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
% write some progress notes to standard output:
fprintf(stdout,'\n'); 
fprintf(stdout,'\n'); 
% print out how much memory has been allocated:  
fprintf(stdout,'Dynamically allocated%d bytes\n', allocatedBytes); 
fprintf(stdout,'\n'); 
% print out some simulation parameters: 
fprintf(stdout,'PLOT_MODULUS =%d\n', PLOT_MODULUS); 
fprintf(stdout,'rectangular waveguide\n'); 
fprintf(stdout,'Stimulus =%g Hertz ', FREQUENCY); 
fprintf(stdout,'continuous plane wave\n'); 
fprintf(stdout,'\n'); 
fprintf(stdout,'Meshing parameters:\n'); 
fprintf(stdout,'%dx%dx%d cells\n', nx, ny, nz); 
fprintf(stdout,'dx=%g, dy=%g, dz=%g meters\n', dx, dy, dz); 
fprintf(stdout,'%g x %g x %g meter^3 simulation region\n',  GUIDE_WIDTH, GUIDE_HEIGHT, LENGTH_IN_WAVELENGTHS*lambda);           
fprintf(stdout,'\n'); 
fprintf(stdout,'Time simulated was %g seconds,%d timesteps\n',    totalSimulatedTime, MAXIMUM_ITERATION); 
fprintf(stdout,'\n'); 
fprintf(stdout,'3D output scaling parameters:\n'); 
fprintf(stdout,'Autoscaling every timestep\n'); 
fprintf(stdout,'\n'); 
% print out simulationMin and simulationMax: 
fprintf(stdout,'Minimum output value was %g \n', simulationMin); 
fprintf(stdout,'Maximum output value was %g \n', simulationMax); 
fprintf(stdout,'\n'); 
% print out a bob command line, including the dimensions  
%      of the output files: 
fprintf(stdout,'bob -cmap chengGbry.cmap -s%dx%dx%d *.bob\n',   nx+1, ny+1, nz); 
fprintf(stdout,'\n'); 
% print out a viz command line: 
fprintf(stdout,'viz ToyFDTD1c.viz\n'); 
fprintf(stdout,'\n'); 
if (exist('interactive'))
    if (exist('h'))
        close(h); % must close the figure otherwise we can run out of memory
    end
    sliceomatic(ez); 
    beep;
    pause;
    h=gcf;          
end