%
%   Read and plot 4D field dump produced by "plot binary Ca fileName"
%          CalC version 7.4 * Victor Matveev * April 18, 2013
%
%-------------------------------------------------------------------------

[status result] = system('.\cwin674.exe simple.par') % Uncomment to run CalC from MATLAB
f  = fopen('Ca_simple.dat', 'rb'); % Read the dump file

G  = fread(f, 1, 'int');        % Read geometry (lowest 2 bits = # of dimensions)

N1 = fread(f, 1, 'int');        % Read number of x-nodes in binary plot
N2 = fread(f, 1, 'int');        % Read number of y-nodes in binary plot
N3 = fread(f, 1, 'int');        % Read number of z-nodes in binary plot

X1 = fread(f, N1, 'double');    % x-axis coordinates
X2 = fread(f, N2, 'double');    % y-axis coordinates
X3 = fread(f, N3, 'double');    % z-axis coordinates

[Y Z]  = meshgrid(X2, X3);

xNode  = 23;                    % Choose your x-slice
xCoord = X1(xNode);             % x-coordinate of the x-slice

N = N1 * N2 * N3;               % Number of nodes
state = 0;

while ~state
    t = fread(f, 1, 'double');  % Read Time
    if feof(f) 
        return;  
    end;
    A = fread(f, N, 'double');  % Read the data
    B = reshape(A, N1, N2, N3); % Reshape into 3D array
    if min(min(min(B))) == max(max(max(B)))
        fprintf(' Skipping time slice t=%g\n', t);
        continue;
    end;
    LogCa = log(0.01 + squeeze(B(xNode,:,:))'); % Take the x-slice, NOTE THE TRANSPOSITION!
    contourf(Y, Z, LogCa);      % plot the x-slice
    colormap hot;  shading flat;
    title([ '[Ca^{2+}] at x = ', num2str(xCoord), ' \mum; time = ', num2str(t), ' ms' ]);
    xlabel('Y (\mum)'); ylabel('Z (\mum)');
    colorbar;
    drawnow;
    pause(0.5);
    state = feof(f);
end

fclose(f);                   % close the file
