diff --git a/misc/binary_reader_example.m b/misc/binary_reader_example.m
new file mode 100644
index 000000000..1d8e78938
--- /dev/null
+++ b/misc/binary_reader_example.m
@@ -0,0 +1,48 @@
+% Example of using binary_reader_wrapper
+% Need to change "format" to 2 when post-processing the example cases
+
+clear; clc; close all;
+
+mfcPath = '..';
+
+%% 1D
+
+binDir = fullfile(mfcPath, 'examples', '1D_acoustic_dipole', 'binary');
+ti = 0;
+tf = 250;
+tDelta = 1;
+
+pres = binary_reader_wrapper(binDir, ti, tf, tDelta, 1);
+
+figure;
+contourf(pres)
+
+%% 2D
+
+binDir  = fullfile(mfcPath, 'examples', '2D_acoustic_support5', 'binary');
+ti = 0;
+tf = 200;
+tDelta = 10;
+
+pres = binary_reader_wrapper(binDir, ti, tf, tDelta, 2);
+
+figure;
+for i = 1:9
+    subplot(3,3,i)
+    contourf(squeeze(pres(:, :, 2*i))');
+end
+
+%% 3D
+
+binDir  = fullfile(mfcPath, 'examples', '3D_acoustic_support7', 'binary');
+ti = 0;
+tf = 100;
+tDelta = 5;
+
+pres = binary_reader_wrapper(binDir, ti, tf, tDelta, 3);
+
+figure;
+for i = 1:9
+    subplot(3,3,i)
+    contourf(squeeze(pres(:, :, 25, 2*i))');
+end
\ No newline at end of file
diff --git a/misc/binary_reader_wrapper.m b/misc/binary_reader_wrapper.m
new file mode 100644
index 000000000..72023300e
--- /dev/null
+++ b/misc/binary_reader_wrapper.m
@@ -0,0 +1,200 @@
+% Wrapper for reading MFC binary output file
+% Works for 1D/2D/3D with multiple processors
+
+function pres = binary_reader_wrapper(binDir, ti, tf, t_delta, dim)
+    % Add more output variables (like rho or xCoords) as desired
+
+    % Total time steps
+    tArr = ti : t_delta : tf;
+    tArrLen = length(tArr);
+
+    if (dim ~= 1)
+        [iProcList, m, n, p, xIdxs, yIdxs, zIdxs, xCoords, yCoords, zCoords] = getProcIdx(binDir, tArr, dim);
+    else
+        % For 1D, the folder 'root' contains all the data
+        if (isempty(binDir))
+            error(strcat("ERROR: invalid binDir (not a binary/root folder): ", binDir))
+        end
+        filename = fullfile(binDir, 'root', [num2str(tArr(1)), '.dat']);
+        dat = f_binary_reader(filename, 'n', 'real*8', 50);
+        m = dat.m;
+        iProcList = 1;
+        xIdxs{1} = 1:m+1;
+    end
+
+    % Loop through files for each time step
+    for tIdx = 1:tArrLen
+        if (mod(tIdx, 10) == 0 || tIdx == 1)
+            disp(['Reading time step ', num2str(tIdx), ' of ', num2str(tArrLen)]);
+        end
+        
+        for iProc = 1 : length(iProcList)
+            if (dim ~= 1)
+                filename = fullfile(binDir, ['p', num2str(iProcList(iProc)-1)], [num2str(tArr(tIdx)), '.dat']);
+            else
+                filename = fullfile(binDir, 'root', [num2str(tArr(tIdx)), '.dat']);
+            end
+            dat = f_binary_reader(filename, 'n', 'real*8', 50);
+
+            xIdx = xIdxs{iProc};
+            nx = m(iProc)+1;
+            if (dim >= 2)
+                yIdx = yIdxs{iProc};
+                ny = n(iProc)+1;
+            end
+            if (dim == 3)
+                zIdx = zIdxs{iProc};
+                nz = p(iProc)+1;
+            end
+            
+            % Add more variables (like rho) as desired
+            if (dim == 1)
+                pres(xIdx, tIdx) = dat.pres;
+            elseif (dim == 2)
+                pres(xIdx, yIdx, tIdx) = reshape(dat.pres, nx, ny);
+            elseif (dim == 3)
+                pres(xIdx, yIdx, zIdx, tIdx) = reshape(dat.pres, nx, ny, nz);
+            end
+        end
+    end
+end
+
+%% Helper Functions
+
+function [iProcList, m, n, p, xIdxs, yIdxs, zIdxs, xCoords, yCoords, zCoords] = getProcIdx(binDir, tArr, dim)
+% Set up global index mapping for multiple processors (only for 2D/3D, so dim == 2 or 3)
+
+% Returns:
+%   iProcList - list of folder indices corresponding to valid processors
+%   m, n, p   - number of cells in each proc (lists)
+%   xIdxs, yIdxs, zIdxs       - global index arrays for each proc
+%   xCoords, yCoords, zCoords - global cell center coordinate arrays
+
+    % List folders according to processor number used for the simulation
+    p_folders = dir( fullfile(binDir, 'p*') ); 
+    nProcFolders = length(p_folders);
+    if (nProcFolders == 0)
+        disp("ERROR: No p_* folders found in binDir:")
+        disp(binDir)
+        error('No processor folders found.');
+    end
+
+    % First get the coord range for each proc (using the first time step)
+    iProcList = [];
+    validProc = 0;
+    for procNum = 1:nProcFolders
+        filename = fullfile(binDir, ['p', num2str(procNum-1)], [num2str(tArr(1)), '.dat']);
+        dat = f_binary_reader(filename, 'n', 'real*8', 50);
+        if ((dat.m == 0) || (dim >= 2 && dat.n == 0) || (dim == 3 && dat.p == 0))
+            continue
+        end
+        validProc = validProc + 1;
+        m(validProc) = dat.m;
+        n(validProc) = dat.n;
+        xCoord{validProc} = dat.x_cb;
+        yCoord{validProc} = dat.y_cb;
+        if (dim == 3)
+            p(validProc) = dat.p;
+            zCoord{validProc} = dat.z_cb;
+        end
+        iProcList(end+1) = procNum;
+    end
+    nProc = validProc;
+    assert(nProc == length(iProcList));
+
+    % Assign empty arrays if not used
+    if dim < 3
+        p = [];
+        zCoord = {};
+    end
+
+    % For each proc we build a list of processors that come before it 
+        % in a given direction by comparing the coordinates
+    xProcOrder = cell(1, nProc);
+    yProcOrder = cell(1, nProc);
+    zProcOrder = {};
+    if (dim == 3)
+        zProcOrder = cell(1, nProc);
+    end
+
+    for iProc = 1 : nProc
+        if (dim == 2)
+            xProcOrder{iProc} = getProcOrder(xCoord, iProc, yCoord);
+            yProcOrder{iProc} = getProcOrder(yCoord, iProc, xCoord);
+        elseif (dim == 3)
+            xProcOrder{iProc} = getProcOrder(xCoord, iProc, yCoord, zCoord);
+            yProcOrder{iProc} = getProcOrder(yCoord, iProc, xCoord, zCoord);
+            zProcOrder{iProc} = getProcOrder(zCoord, iProc, xCoord, yCoord);
+        end
+    end
+
+    % Using the order, we compute the global indices using the proc orders
+    for iProc = 1 : nProc
+        xLocal{iProc} = 1:(m(iProc)+1);
+        yLocal{iProc} = 1:(n(iProc)+1);
+        if (dim == 3)
+            zLocal{iProc} = 1:(p(iProc)+1);
+        end
+    end
+
+    xIdxs = getGlobalIdx(xLocal, xProcOrder);
+    yIdxs = getGlobalIdx(yLocal, yProcOrder);
+    zIdxs = {};
+    if (dim == 3)
+        zIdxs = getGlobalIdx(zLocal, zProcOrder);
+    end
+
+    % Get global cell center locations
+    for iProc = 1 : nProc
+        xCoords( xIdxs{iProc} ) = xCoord{iProc}(1:end-1) + diff(xCoord{iProc})/2;
+        yCoords( yIdxs{iProc} ) = yCoord{iProc}(1:end-1) + diff(yCoord{iProc})/2;
+        if (dim == 3)
+            zCoords( zIdxs{iProc} ) = zCoord{iProc}(1:end-1) + diff(zCoord{iProc})/2;
+        end
+    end
+    % Assign empty outputs if not used
+    if dim < 3
+        zCoords = [];
+    end
+end
+
+function order = getProcOrder(coordCell, iProc, varargin)
+% For the processor, compute a list of processors that come it based on their coordinates
+
+    order = [];
+    currentCoord = coordCell{iProc}(1);
+    nProc = numel(coordCell);
+    for j = 1:nProc
+        if j == iProc
+            continue;
+        end
+        if coordCell{j}(1) < currentCoord
+            valid = true;
+            for k = 1:length(varargin)
+                if varargin{k}{iProc}(1) ~= varargin{k}{j}(1)
+                    valid = false;
+                    break;
+                end
+            end
+            if valid
+                order(end+1) = j;
+            end
+        end
+    end
+end
+
+function globalIdxs = getGlobalIdx(localIdxs, procOrder)
+% Compute the global indices for each processor
+
+    nProc = length(localIdxs);
+    globalIdxs = cell(1, nProc);
+    for i = 1:nProc
+        offset = 0;
+        if ~isempty(procOrder{i})
+            for j = procOrder{i}
+                offset = offset + length(localIdxs{j});
+            end
+        end
+        globalIdxs{i} = localIdxs{i} + offset;
+    end
+end
\ No newline at end of file