function zz = Examine_Intan_RHS_File_with_Octave__20230324_BMSEED()
close all; clc;
pkg load signal      %~~~ for Octave only  ~~~%

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
% Edited by BMSEED (update on 03/24/22)                                   %
%   - Slight modifications to make code more general and to plot/save     %
%       stimulation and analog input signals.                             %
%   - Added code to examine recordings from gen7a sMEAs.				  %
%                                                                         %
% Edited by BMSEED (update on 12/13/22)                                   %
%   - Added code to rearrange gen6 and glass MEAs so that all reference   %
%       channels are listed and plotted at the end of the list.	          %
%   - Added code to include vertical lines on the amplifier signals plot  %
%       to show what section of the electrodes' signals were used for     %
%       Vrms calculations.                                                %
%                                                                         %
% Edited by BMSEED (update on 12/09/22)                                   %
%   - Added code to save the impedance and Vrms values for all electrodes %
%       to a csv file and handle the short custom names ("E01", etc).	  %
%                                                                         %
% Edited by BMSEED (update on 02/22/22)                                   %
%   - Added code to identify the channel number/id on a more general way. %
%   - Added code (commented-out) to resize the amplifier data plots and   %
%       save these figures as png.                                        %
%                                                                         %
% Edited by BMSEED (update on 01/19/22)                                   %
%   - Updated code so that it is more general and can be used with any    %
%       type of sample (e.g. sMEA gen4e vs glass MEA)                     %
%   - Figure plotting the amplifier channels' data now displays the       %
%       impedance measurement prior to signal recording, as well as Vrms  %
%       measurement from the recording itself.                            %
%   - Exported output data will now exclude all stimulation channels, as  %
%       well as the analog input channel, that were not used (i.e.,       %
%       contain all zeros)                                                %
%   - Edited code so that it compatible in Matlab and Octave.             %
%                                                                         %
% Edited by BMSEED (update on 08/26/21)                                   %
%   - Added code to ask the user for the y-axis scale for the plots of    %
%       the gen4a-d sMEA channels' data.                                  %
%                                                                         %
% Edited by BMSEED (update on 06/01/21)                                   %
%   - Added code to store impedance magnitude values for all channels in  %
%       a variable called "Z_data".                                       %
%                                                                         %
% Edited by BMSEED (update on 02/08/19)                                   %
%   - Added code for situations when no analog input is present.  The     %
%       code creates a vector array of zeros with approriate length. This %
%       helps to eliminate formatting issues in the output text file.     %
%   - Minor fix to the output text file.  Eliminated the extra blank line %
%       that the previous version of the code was adding after each data  %
%       line.                                                             %
%                                                                         %
% Edited by BMSEED (last update on 9/12/17)                               %
%   - Output text files will now include stimulation signals for all 30   %
%       chanels and the analog input signal that triggers the stimulation.%
%   - Instead of a single output text file, multiple smaller text files   %
%       may be created to prevent excessively large files.  The user will %
%       define the time interval of each text file.                       %
%                                                                         %
% Function Description:                                                   %
%   This function will read an RHS data file and plot results.            %
%   Three figures will be created                                         %
%   for each trial:                                                       %
%     % Figure #1: raw biopotential (µV) signals                          %
%     % Figure #2: stimulation signals                                    %
%     % Figure #3: analog input signal (e.g. from the function generator) %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%




%% READING THE RHS DATA FILE USING THE M-FILE PROVIDED BY INTAN TECHNOLOG.
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% Loading variables from the signal (i.e. voltage) data file
[amplifier_ch, amplifier_data, fr_para, ~, ain_data, t, ...
    stim_data, file, ~] = read_Intan_RHS2000_file__BMSEED();
disp(['V File:  ' file]);
disp(' ');
Fs = fr_para.amplifier_sample_rate;     % sampling frequency


% Asking the user for input
prompt = {sprintf('%s\n%s\n\n\n%s\n%s\n%s\n%s\n%s\n%s', ...
                'Intan File: ', file, ...
                'Input the type of MEA (one of the following):', ...
                '    - sMEA_gen4', '    - sMEA_gen4e', ...
                '    - sMEA_gen6', '    - sMEA_gen7a', ...
                '    - glass_MEA'       ), ...
          sprintf('\n\n%s\n%s', 'For plotting the data:', ...
                            '    y-axis minimum value (uV):'), ...
          sprintf('%s',     '    y-axis maximum value (uV):'), ...
          sprintf('\n\n%s\n%s', 'For Vrms calculations:', ...
                            '    second marker to start calc (sec):'), ...
          sprintf('%s',     '    second marker to end calc (sec):')   };
dlgtitle = 'User Input';
dims = 1;
    %~~~  for Matlab --> dims = [1 50];  ~~~%
    %~~~  for Octave --> dims = 1;       ~~~%
definput = {'sMEA_gen4e','-100','100','1.0','1.5'};
answer = inputdlg(prompt,dlgtitle,dims,definput);
if isempty(answer) == 1
    answer = definput;
end


% Defining variables based on user input
MEA_type = answer{1};
y_limits = [str2double(answer{2})  str2double(answer{3})];
start_sec_mark = str2double(answer{4});
end_sec_mark = str2double(answer{5});


% Defining the number of channels based on the type of MEA being used
switch MEA_type
    case "sMEA_gen4"
        num_chs = 28+2;     % 28 recording channels, 2 reference channels
    case "sMEA_gen4e"
        num_chs = 28+4;     % 28 recording channels, 4 reference channels
    case "sMEA_gen6"
        num_chs = 56+4;     % 56 recording channels, 4 reference channels
    case "sMEA_gen7a"
        num_chs = 30+4;     % 30 recording channels, 4 reference channels
    case "glass_MEA"
        num_chs = 59+1;     % 59 recording channels, 1 reference channels
    otherwise
        msgbox('Invalid MEA type. Function execution will stop.', ...
               'Error', 'error');
        return
end


% Initializing variables prior to mapping/sorting data columns
id_num = zeros(num_chs,1);
data = zeros( length(t), num_chs);
stim = zeros( length(t), num_chs);
Z_data = zeros( num_chs, 1);
Vrms = zeros( num_chs, 1);

% Getting the mapping to go from element number to sMEA #
n = length(amplifier_ch);
cnt = 1;
for k = 1:n

    % Ignoring ground channels
    temp = amplifier_ch(k).custom_channel_name;
    if length(temp) >= 4    % to take into account short custom channel names
      if strcmp(temp(1:4), 'Plex') == 1 | strfind(temp, 'GND') > 0
          continue
      end
	%else
	%  if strcmp(temp(1:2), 'IB') == 1
	%	  continue
	%  end
    end


    % Looking for the first numerical character in channel name string
    char_type_array = isstrprop(temp,'digit');
    num_char_id = find(char_type_array == 1, 2);%1, 'first'); % 2);

    if char_type_array(1) == 1
        num_char_id = num_char_id(2);     % if the first character is a number
    else
        num_char_id = num_char_id(1);
    end

    if temp(num_char_id - 1) == 'r' || temp(num_char_id - 1) == 'R'
        % if the channel is a reference electrode
        id = temp( num_char_id-1:num_char_id   );
    elseif temp(num_char_id - 1) == 'C' || temp(num_char_id - 1) == 'P'
        % if the channel is labeled "C" or "P" in a gen7a sMEA
        id = temp( num_char_id-1:num_char_id+1 );
    else
        % if the channel is a recording electrode (i.e., not a reference)
        id = temp(   num_char_id:num_char_id+1 );
    end


    % Converting string array to numerical array
    if      strcmp(id,'r1') == 1 || strcmp(id,'R1') == 1    % Ref. ch. #1
        switch MEA_type
            case "sMEA_gen4",   id_num(cnt,1) = 29;         % sMEA gen4
            case "sMEA_gen4e",  id_num(cnt,1) = 29;         % sMEA gen4e
            case "sMEA_gen6",   id_num(cnt,1) =  7;         % sMEA gen6
            case "sMEA_gen7a",  id_num(cnt,1) = 31;         % sMEA gen7a
            otherwise,          id_num(cnt,1) = 15;         % glass MEA
        end

    elseif  strcmp(id,'r2') == 1 || strcmp(id,'R2') == 1    % Ref. ch. #2
        switch MEA_type
            case "sMEA_gen4",   id_num(cnt,1) = 30;         % sMEA gen4
            case "sMEA_gen4e",  id_num(cnt,1) = 30;         % sMEA gen4e
            case "sMEA_gen7a",  id_num(cnt,1) = 32;         % sMEA gen7a
            otherwise,          id_num(cnt,1) = 25-1;       % sMEA gen6
        end

    elseif  strcmp(id,'r3') == 1 || strcmp(id,'R3') == 1    % Ref. ch. #3
        switch MEA_type
            case "sMEA_gen4e",  id_num(cnt,1) = 31;         % sMEA gen4e
            case "sMEA_gen7a",  id_num(cnt,1) = 33;         % sMEA gen7a
            otherwise,          id_num(cnt,1) = 40-2;       % sMEA gen6
        end

    elseif  strcmp(id,'r4') == 1 || strcmp(id,'R4') == 1    % Ref. ch. #4
        switch MEA_type
            case "sMEA_gen4e",  id_num(cnt,1) = 32;         % sMEA gen4e
            case "sMEA_gen7a",  id_num(cnt,1) = 34;         % sMEA gen7a
            otherwise,          id_num(cnt,1) = 57-3;       % sMEA gen6
        end

    else          % if not a reference channel
        id_num(cnt,1) = str2double(id);

        % additional renumbering for gen6 sMEAs
        if strcmp(MEA_type,'sMEA_gen6') == 1
            if str2double(id) >= 17 && str2double(id) <= 31
                  id_num(cnt,1) = str2double(id)-1;
            elseif  str2double(id) >= 33 && str2double(id) <= 47
                  id_num(cnt,1) = str2double(id)-2;
            elseif  str2double(id) >= 49 && str2double(id) <= 63
                  id_num(cnt,1) = str2double(id)-3;
            end
        end
        
        % additional renumbering for gen7a sMEAs
        if strcmp(MEA_type,'sMEA_gen7a') == 1
            if id(1) == 'C'
                  id_num(cnt,1) = str2double(id(2:3));
            elseif  id(1) == 'P'
                  id_num(cnt,1) = 10 + str2double(id(2:3));
            end
        end

    end


    % Re-organizing the channels data to match the sMEA ID
    data(:,id_num(cnt,1)) = amplifier_data( k, :)';
    stim(:,id_num(cnt,1)) = stim_data( k, :)';
    a_in = ain_data';
    Z_data(id_num(cnt,1),1) = amplifier_ch(k).electrode_impedance_magnitude;


    % Updating the arrays with the channel labels
    if strcmp(id(1),'r') == 1 || strcmp(id(1),'R') == 1
        id_label{id_num(cnt,1),1} = ['REF ' id(2)];
        id_header{id_num(cnt,1),1} = ['REF_' id(2)];
    else
        id_label{id_num(cnt,1),1} = ['E' id];
        id_header{id_num(cnt,1),1} = ['E' id];     % ['E_' id]
    end

    % Updating specific labels on gen6 sMEAs and glass MEAs
    if      strcmp(MEA_type,'sMEA_gen6') == 1 && strcmp(id,'07') == 1
                    id_label{id_num(cnt,1),1}  = 'REF 1';
                    id_header{id_num(cnt,1),1} = 'REF_1';
    elseif  strcmp(MEA_type,'sMEA_gen6') == 1 && strcmp(id,'25') == 1
                    id_label{id_num(cnt,1),1}  = 'REF 2';
                    id_header{id_num(cnt,1),1} = 'REF_2';
    elseif  strcmp(MEA_type,'sMEA_gen6') == 1 && strcmp(id,'40') == 1
                    id_label{id_num(cnt,1),1}  = 'REF 3';
                    id_header{id_num(cnt,1),1} = 'REF_3';
    elseif  strcmp(MEA_type,'sMEA_gen6') == 1 && strcmp(id,'57') == 1
                    id_label{id_num(cnt,1),1}  = 'REF 4';
                    id_header{id_num(cnt,1),1} = 'REF_4';
    elseif  strcmp(MEA_type,'glass_MEA') == 1 && strcmp(id,'15') == 1
                    id_label{id_num(cnt,1),1}  = 'REF 1';
                    id_header{id_num(cnt,1),1} = 'REF_1';
    end


    % Updating the counter
    cnt = cnt + 1;
end



% Rearranging data so that all reference electrodes are listed towards the end
[~, i] = sort(id_label);      % sorting the labels alphabetically
id_label = id_label(i);       % rearranging the id labels
id_header = id_header(i);     % rearranging the id headers for output file

data = data(:,i);             % rearranging the columns of the data 2d array
stim = stim(:,i);             % rearranging the columns of the stim 2d array
Z_data = Z_data(i,1);         % rearranging the impedance data



% V_rms calculations
for k = 1:num_chs
    %~~~  For Octave, simply use the rms() function  ~~~%
    Vrms(k,1) = rms( data(find(t == start_sec_mark):find(t == end_sec_mark), k) );

    %~~~  For Matlab, calcuate Vrms "manually"  ~~~%
    %{
    sum_sq = 0;
    N = (end_sec_mark-start_sec_mark)*Fs + 1;
    for kk = find(t == start_sec_mark):find(t == end_sec_mark)
        sum_sq = sum_sq + (data(kk,k))^2;
    end
    Vrms(k,1) = sqrt(sum_sq/N);
    %}
end

% defining the "time" variable
time = t';


% Cleaning the workspace
clear  fr_para  temp id  amplifier_ch  amplifier_data  t_amplifier
clear  ain_data  stim_data  k  cnt  t  N  sum_sq




%% EXPORTING THE Z AND VRMS MEASUREMENTS
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% preparing the data to export
Z_data_kOhms = Z_data/1000;                 % impedance in kohms
output_Z_Vrms_data{1,1} = 'Electrode';      % headers
output_Z_Vrms_data{1,2} = 'Z (kohms)';
output_Z_Vrms_data{1,3} = 'Vrms (uV)';
output_Z_Vrms_data = [output_Z_Vrms_data; ...
                      id_label num2cell(Z_data_kOhms) num2cell(Vrms)];

% defining the output file's name
output_Z_Vrms_file = [file(1:end-4) '_Z_and_Vrms.csv'];

%~~~  For Octave, load the "io" package and the "cell2csv" function  ~~~%
pkg load io
cell2csv(output_Z_Vrms_file, output_Z_Vrms_data);

%~~~  For Matlab, use other functions (... needs to be tested ...)  ~~~%
% T = cell2table(output_Z_Vrms_data);
% writetable(T, output_Z_Vrms_file);




%% PLOTTING THE AMPLIFIER DATA, STIMULUS DATA, AND ANALOG INPUT
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% To start measuring the time elapsed to plot figures
tStartFig = tic;

% Preparing for the figures
step = 1;
    %~~~  for Matlab --> step = 1;  ~~~%
    %~~~  for Octave --> step = 10;  ~~~%
switch MEA_type
    case "sMEA_gen4",   fig_num = 1;
                        fig_rows = 5;
                        fig_plots = 30;
    case "sMEA_gen4e",  fig_num = 1;
                        fig_rows = 6;
                        fig_plots = 32;
    case "sMEA_gen7a",  fig_num = 1;
                        fig_rows = 6;
                        fig_plots = 34;
    otherwise,          fig_num = 2;      % sMEA gen6, glass MEA
                        fig_rows = 5;
                        fig_plots = 30;
end


% Defining the start and end time markers used for Vrms calculations
start_sec_indx = find(time == start_sec_mark);    % for Vrms start marker
end_sec_indx = find(time == end_sec_mark);        % for Vrms end marker


% Plotting the amplfier data
for j = 1:fig_num
  %figure(j+0);
  figure(j+0, 'position',[1,62,1918,940]);

  for k = 1:fig_plots
      subplot(fig_rows,6,k); hold on;
		%plot(time(start_sec_indx:step:end_sec_indx), data(start_sec_indx:step:end_sec_indx,30*(j-1)+k), ...
	  plot(time(1:step:end), data(1:step:end,30*(j-1)+k), ...
           'Color', 'b', 'LineWidth', 1);
      plot(time(start_sec_indx)*ones(1,2), [y_limits(1) y_limits(2)], 'r:');
      plot(  time(end_sec_indx)*ones(1,2), [y_limits(1) y_limits(2)], 'r:');
      title([id_label{30*(j-1)+k} ': ' ...
             num2str(round(Z_data(30*(j-1)+k)/1000)) 'k\Omega, ' ...
             num2str(round(Vrms(30*(j-1)+k))) '\muV'])
      ylim([y_limits(1) y_limits(2)]);
      %xlim([start_sec_mark end_sec_mark]);
	  %xlim([33 45]);
  end

  h = figure(j+0);
  if fig_num == 1,    file2 = file(1:end-4);
  else,               file2 = [file(1:end-4) '__part' num2str(j+0)];  
  end
  
  saveas(h, file2, 'png');
end
% close all



%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
% Code will stop here         %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                             %
%if 1 > 0, return end
%                             %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%



% Plotting the stimulus data
for j = 1:fig_num
  %figure(j+2);
  figure(j+2, 'position',[1,62,1918,940]);
  
  for k = 1:fig_plots
      subplot(fig_rows,6,k); hold on;
      if ~any(stim(:,30*(j-1)+k))       % if all elements are 0
          plot([time(1) time(end)], [0 0], 'Color', 'b', 'LineWidth', 1);
      else
          plot(time(1:step:end), stim(1:step:end,30*(j-1)+k), ...
               'Color', 'b', 'LineWidth', 1);
      end
      title(id_label{30*(j-1)+k});
      %xlim([0 2]); ylim([-50 50]);
  end
  
  h = figure(j+2);
  if fig_num == 1,    file3 = [file(1:end-4) '_stim'];
  else,               file3 = [file(1:end-4) '_stim__part' num2str(j+0)];
  end
  
  %saveas(h, file3, 'png');
end
%close all


% Plotting the analog input data, if any
if isempty(a_in) == 0
    figure(5);     % Plot of the analog input (in V), if any
    plot(time, a_in, 'Color', 'b', 'LineWidth', 1);
    title('Analog Input 1');
    ylabel('signal (V)');
    xlabel('time (s)');

    h = figure(5);
    file4 = [file(1:end-4) '_analog_input_1'];
    %saveas(h, file4, 'png');
    %close all
    
else
    a_in = zeros(length(time),1);
end



% To stop measuring the time elapsed to plot figures
tEndFig = toc(tStartFig);
fprintf('\n%s%s%s\n', 'Time elapsed to plot figures: ', ...
                      num2str(tEndFig), ' sec');



%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
% Code will stop here         %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                             %
%if 1 > 0, return end
%                             %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%




%% PREPARING TO EXPORT DATA TO A TEXT FILE
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% Compiling the time and amplifier data into a single var
exp_data = [time  data];

% Creating the header for amplifier channel data columns
data_header_txt = "";     % stores headers from all iterations
data_formatting = "";     % stores formatting info from all iterations
for k = 1:num_chs
    data_header_txt = [data_header_txt '\t' id_header{k} ' '];
    data_formatting = [data_formatting '\t% .4f '];
end


% Removing stimulation channels with all zeros (i.e., unused)
stim_data_cols_w_all_0s = all(stim == 0);
exp_stim = stim(:, ~stim_data_cols_w_all_0s);    % trimmed stim data

% Creating the header for the stim data columns
stim_header_txt = "";     % stores headers from all iterations
for k = 1:length(stim_data_cols_w_all_0s)
    if stim_data_cols_w_all_0s(k) == 0    % if stim data is not all 0's
        stim_header_txt = [stim_header_txt '\tSt_' id_header{k} ' '];
        data_formatting = [data_formatting '\t% .4f '];
    end
end


% Updating the compiled data, header, and formatting variables, depending
%   if the analog input channel is empty or not
if any(a_in)            % if there is at least one non-zero element
    exp_data = [exp_data  exp_stim  a_in];
    a_in_header_txt = '\tAI_01 ';
    data_formatting = [data_formatting '\t% .4f '];
else
    exp_data = [exp_data  exp_stim];
    a_in_header_txt = "";
    % no need to updata the "data_formatting" variable
end


% Output file variables
outfile_header = ['Time '  data_header_txt  stim_header_txt  ...
                  a_in_header_txt  '\n'];
outfile_formatting = ['%f '  data_formatting  '\n'];


%~~~  For Matlab only, concatenate the string array elements  ~~~%
%{
temp = "";
for j = 1:length(outfile_header)
    temp = strcat(temp, outfile_header(j));
end
outfile_header = temp;

temp = "";
for j = 1:length(outfile_formatting)
    temp = strcat(temp, outfile_formatting(j));
end
outfile_formatting = temp;
%}




%% EXPORTING DATA TO A TEXT FILE
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% Getting user input on the filename
clear answers
prompt = {'Enter base filename (exclude ".txt"):', ...
      'Time interval, in seconds, per output file (whole number only):'};
dlg_title = 'Output File Name';
dims = 1;
    %~~~  for Matlab --> dims = [1 65];  ~~~%
    %~~~  for Octave --> dims = 1;       ~~~%
defaultans = {[file(1:end-4) '_data'], '5'};
options.Resize = 'on';
answers = inputdlg(prompt, dlg_title, dims, defaultans);
    %~~~  for Matlab --> answers = inputdlg(prompt, dlg_title, dims, defaultans, options);  ~~~%
    %~~~  for Octave --> answers = inputdlg(prompt, dlg_title, dims, defaultans);  ~~~%

% Error checking the user input
if isempty(answers) == 1
    msgbox('No output filename specified. Code execution will stop.', ...
               'Error', 'error');
    return
end
file = answers{1};                % converting to string variable
intv = str2double(answers{2});    % time interval (in sec) for separation


% Parameters to create the multiple data files
rmnder = mod(exp_data(:,1)-exp_data(1,1), intv);    % remainder
indx = find(rmnder == 0);                           % indices to break file
% time_str = "";        %~~~  for Matlab only  ~~~%
for k = 1:length(indx)
    time_str(k) = {num2str( time(indx(k)) )};
end


% To start measuring the time elapsed to export data
tStartExport = tic;


% Creating the output text files
disp('Writing output files:  ');
for k = 1:length(indx)

    % Extracting the subset of the data that will be exported
    if k < length(indx)
        temp_data = exp_data( indx(k):step:indx(k+1), :);
            %~~~  for storing the full data --> (     indx(k):indx(k+1), :)  ~~~%
            %~~~  for Octave troubleshooting --> (indx(k):step:indx(k+1), :)  ~~~%
    else
        temp_data = exp_data(       indx(k):step:end, :);
            %~~~  for storing the full data --> (      indx(k):end, :)  ~~~%
            %~~~  for Octave troubleshooting --> ( indx(k):step:end, :)  ~~~%
    end


    % Specifying the output file name
    if k < length(indx)
        temp_file = [file '_' num2str(k) '_-_' ...
                     time_str{k} 'to' time_str{k+1} 'sec'];
    else
        temp_file = [file '_' num2str(k) '_-_' ...
                     time_str{k} 'sectoEND'];
    end
    disp(['    ' temp_file '.txt']);


    % Writing the output data file
    fid = fopen( [temp_file '.txt'], 'wt' );
    fprintf(fid, outfile_header);
    fprintf(fid, outfile_formatting, temp_data');
    fclose(fid);
end


% To stop measuring the time elapsed to export data
tEndExport = toc(tStartExport);




%% FUNCTION OUTPUT
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%fprintf('\n%s%s%s\n', 'Time elapsed to plot figures: ', ...
%                      num2str(tEndFig), ' sec');
fprintf(  '%s%s%s\n', 'Time elapsed to export data: ',  ...
                      num2str(tEndExport), ' sec');
fprintf(  '\n%s\n\n', 'Code execution completed.');



end
