function [h]=my_ISM(fs,mic,n,wallrefl,rm,src)

%   Generates a room impulse response(RIR) calculation by means of the
%   image method. This code is a modification of a prior existing code
%   which calculates an RIR with the assumption of all reflection
%   coefficients from all enclosing boundaries are equal. This modified 
%   code allows input of individual reflection coefficients for each
%   boundary to more accurately calculate the RIR and implements a phase
%   inversion for each sound reflection off each boundary which results in
%   more accurate representation of true RIRs recorded in real world
%   acoustic environments.
%
%   Input Arguments:
%
%   fs       = sample rate.
%   mic      = row vector giving the x,y,z coordinates of the microphone position.  
%   n        = This function will calculate (2*n+1)^3 virtual sources 
%   wallrefl = row vector [x1,x2,y1,y2,z1,z2] reflection coefficients
%            for the walls, for this function only positive values of r (0<r<1)
%            are assumed.
%   rm       = row vector giving the dimensions of the room.  
%   src      = row vector giving the x,y,z coordinates of the sound source.
%
%   Example:
%
%   fs=44100;
%   mic=[4.5 3 1.3];
%   n=24;
%   wallrefl=[0.75 0.65 0.8 0.2 0.45 0.7];
%   rm=[5 3.5 2.7];
%   src=[0.5 0.5 1.3];
%
%   All distances are in meters.
%
%   The output is scaled such that the largest value of the 
%   absolute value of the output vector is equal to one.
%
%   This code is a modification of an existing code: Version 3.4.1
%   Copyright © 2003 Stephen G. McGovern. 


nn=-n:1:n;                        % The index for the sequence: a row
%                                   vector from -'n'th virtual source to
%                                   the 'n'th virtual source where 0 is the
%                                   real source.

rms=nn+0.5-0.5*(-1).^nn;          % Part of equation for calculating 
%                                   coordinate position (x, y or z)

srcs=(-1).^(nn);                  % Part of equation for calculating 
%                                   coordinate position (x, y or z)

xi=srcs*src(1)+rms*rm(1)-mic(1);  % Distance between 'i'th virtual source
%                                   and mic position in the x axis. 

yj=srcs*src(2)+rms*rm(2)-mic(2);  % Distance between 'j'th virtual source
%                                   and mic position in the y axis. 

zk=srcs*src(3)+rms*rm(3)-mic(3);  % Distance between 'k'th virtual source
%                                   and mic position in the z axis. 

[i,j,k]=meshgrid(xi,yj,zk);       % converts vectors to 3D matrices for
%                                   calculation of 3D virtual sources

d=sqrt(i.^2+j.^2+k.^2);           % straight line distance calculation
%                                   to each virtual source using
%                                   Pythagoras' theorem

time=round(fs*d/343)+1;           % time domain IR function of each virtual
%                                   source

[e,f,g]=meshgrid(nn, nn, nn);     % converts vectors to 3D matrices for
%                                   calculation of 3D virtual sources

wallrefl = -abs(wallrefl);        % implementation of phase inversion as
%                                   per Lehmann & Johansson's ISM structure

%  Combined reflection coefficient of all reflections of virtual sources in x axis
rx=wallrefl(1).^(abs(0.5.*e-0.25+0.25.*(-1).^e)).*wallrefl(2).^(abs(0.5.*e+0.25-0.25.*(-1).^e));

%  Combined reflection coefficient of all reflections of virtual sources in y axis
ry=wallrefl(3).^(abs(0.5.*f-0.25+0.25.*(-1).^f)).*wallrefl(4).^(abs(0.5.*f+0.25-0.25.*(-1).^f));

%  Combined reflection coefficient of all reflections of virtual sources in z axis
rz=wallrefl(5).^(abs(0.5.*g-0.25+0.25.*(-1).^g)).*wallrefl(6).^(abs(0.5.*g+0.25-0.25.*(-1).^g));

c=rx.*ry.*rz;                     % Total reflection coefficient for every
                                  % virtual source

ee=c./d;                          % Total magnitude of each echo

h=full(sparse(time(:),1,ee(:)));  % Impulse Response: Total magnitude multiplied
%                                   by unit impulse response function.

h=h/max(abs(h));                  % Scale output