function f = simod( ss, mm, n, ref, sz, plf)

% Default values for room scale
% and reference distance
if nargin == 3
	sz = 1;
	ref = sqrt( sum( (ss-mm).^2));
	plf = 1;
end
if nargin == 4
	sz = 1;
	plf = 1;
end
if nargin == 5
	plf == 1;
end

% Make sure inputs are of
% the same dimensionality
if length( ss) ~= length( mm)
	error('Non-conformant inputs')
end

% Accept only odd size grids
% so that we have a center room
if ~mod( n, 2)
	n = n + 1;
	disp( 'Bumped order up by 1')
end

% Default room size (in meters)
xm = 1;
ym = 1;
zm = 1;

% Bring input to scale
ss = ss / sz;
mm = mm / sz;

% Check positions validity
for i = 1:length(ss)
  if ss(i) > 1
    error( 'Bad source coordinates');
  end
  if mm(i) > 1
    error( 'Bad mike coordinates');
  end  
end

% Prepare graphics
if n > 7
	plf = 0;
end
if plf
	clf
	hold on
	axis equal
end

% 2-D room case
if length( ss) == 2

	% Plot unfolded rooms and
	% get their origin points
	xo = zeros( 1, n);
	yo = zeros( 1, n);
	for i = 1:n
		xo(i) = i*xm-1;
		for j = 1:n
			if plf
				plot( [i*xm i*xm], [0 j*ym])
				plot( [0 i*xm], [j*ym j*ym])
			end
			yo(j) = j*ym-1;
			if plf
				plot( xo(i), yo(j), 'x')
			end
		end
	end

	% Place source and observer
	% In center room
	n2 = round( n/2);
	if plf
		plot( xo(n2)+ss(1), yo(n2)+ss(2), '*')
		plot( xo(n2)+mm(1), yo(n2)+mm(2), 'o')
	end

	% Place observers in virtual rooms
	off = 0;
	if ~mod( floor(n/2), 2)
		off = 1;
	end
	xi = zeros( 1, n);
	yi = zeros( 1, n);
	for i = 1:n
		if mod( i+off,2)
			xx = xm-ss(1);
		else
			xx = ss(1);
		end
		xi(i) = xo(i)+xx;
		for j = 1:n
			if mod( j+off,2)
				yy = ym-ss(2);
			else
				yy = ss(2);
			end
			yi(j) = yo(j)+yy;
			if plf
				plot( xi(i), yi(j), 'x')
			end
		end
	end	

	% Get distances of each virtual
	% source from the observer
	xs = xo(n2)+mm(1);
	ys = yo(n2)+mm(2);
	k = 1;
	d = zeros( 1, n*n);
	for i = 1:n
    for j = 1:n
			d(k) = sqrt( (xi(i)-xs)^2 + (yi(j)-ys)^2);
			k = k + 1;
		end
	end

elseif length( ss) == 3
	% 3-D room
	% Accept only odd size grids
	% so that we have a center
	if ~mod( n, 2)
		n = n + 1;
		disp( 'Bumped order up by 1')
	end

	% Plot unfolded rooms and
	% get their origin points
	xo = zeros( 1, n);
	yo = zeros( 1, n);
	zo = zeros( 1, n);
	for i = 1:n	
		xo(i) = i*xm-1;
		for j = 1:n
			yo(j) = j*ym-1;
			for k = 1:n
				if plf
					plot3( [0 i*xm],    [j*ym j*ym], [k*zm k*zm])
					plot3( [i*xm i*xm], [0 j*ym],    [k*zm k*zm])
					plot3( [i*xm i*xm], [j*ym j*ym], [0 k*zm])
				end
				zo(k) = k*zm-1;
%				plot3( xo(i), yo(j), zo(k), 'x')
			end
		end
	end
	if plf
		view( 166, 16)
	end

	% Place source and observer
	% In center room
	n2 = round( n/2);
	if plf
		plot3( xo(n2)+ss(1), yo(n2)+ss(2), zo(n2)+ss(3), 'g*')
		plot3( xo(n2)+mm(1), yo(n2)+mm(2), zo(n2)+mm(3), 'go')
	end

	% Place observers in virtual rooms
	off = 0;
	if ~mod( floor(n/2), 2)
		off = 1;
	end
	xi = zeros( 1, n);
	yi = zeros( 1, n);
	zi = zeros( 1, n);
	for i = 1:n
		if mod( i+off,2)
			xx = xm-ss(1);
		else
			xx = ss(1);
		end
		xi(i) = xo(i)+xx;
		for j = 1:n
			if mod( j+off,2)
				yy = ym-ss(2);
			else
				yy = ss(2);
			end
			yi(j) = yo(j)+yy;
			for k = 1:n
				if mod( k+off,2)
					zz = zm-ss(3);
				else
					zz = ss(3);
				end
				zi(k) = zo(k)+zz;
				if plf
					plot3( xi(i), yi(j), zi(k), 'rx')
				end
			end
		end
	end	

	% Get distances of each virtual
	% source from the observer
	xs = xo(n2)+mm(1);
	ys = yo(n2)+mm(2);
	zs = zo(n2)+mm(3);
	l = 1;
	d = ones( 1, n*n*n);
	for i = 1:n
		for j = 1:n
      for k = 1:n
				d(l) = sqrt( (xi(i)-xs)^2 + (yi(j)-ys)^2 + (zi(k)-zs)^2);
				l = l + 1;
			end
		end
	end
else
	if length( ss) == 1
		error( 'One dimensional rooms not supported')
	else
		error( sprintf('Have you seen many rooms with %d dimensions?', length(ss)))
	end
end

% Convert distance to time
% assume sampling rate of 44.1kHz and
% compute reflection tap addresses
t = round( (sz*d / 300) * 44100);

% Make the filter, assume the 
% intensity for each reflection
% is the inverse of the distance squared
% Flip sign to simulate bounce inversion
c = d / ref;
f(t) = square(linspace(0,pi*length(c),length(c))) ./ (c.^2);

% Release graphics
if plf
	hold off
end