% This program is an addition to ICCM_Dose_Calcs_4_Radio_Beta_Emitters.m.
% A brief introduction to this program is provided below:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%************************************************************************%
%                                                                        %
% Calculation of Absorbed Dose in Examination of Bystander Effects via   %
% --------------------------------------------------------------------   %
% Irradiation with Beta Emitting Radiopharmaceutical                     %
% --------------------------------------------------                     %
%                                                                        %
% This program is designed to calculate the absorbed dose to a set of    %
% donor cells used in the creation of Irradiated Cell Conditioned Medium %
% (ICCM)as per radiation bystander protocol utilized by Boyd et al in:   %
%                                                                        %
% Boyd M, Ross SC, et al. Radiation-Induced biological bystander effect  %
% elicited in vitro by targeted radiopharmaceuticals labeled with alpa-, %
% beta-, and auger electron-emitting radionuclides. J.Nucl.Med. Vol. 47. %
% pp. 1007-1015 (2006).                                                  %
%                                                                        %
% The application was designed with the purpose of performing this       %
% calculation of absorbed dose from the beta emitting                    %
% radiopharmaceutical I-131 MIBG.  This is accomplished through modelling%
% of in vitro conditions and procedures coupled with the appropriate     %
% utilization of the Vynckier-Wambersie (VW) point-source dose function  %
% (kernel) integrated over the appropriate geometry.                     %
%                                                                        %
% Having said this, this program is also believed to be suitable for     %
% other beta emitting radiopharmaceuticals potentially used for ICCM     %
% production with maximum energies ranging from 0.5 to 3.0 MeV and       %
% utilization of Boyd et al' s methodology outlined in the paper above.  %
%                                                                        %
% Created By: Michael Gow                                                %
%             Dept. of Medical Physics and Applied Radiation Sciences    %
%             McMaster University                                        %
%             Hamilton, ON CA                                            %
%             MMXI                                                       %
%                                                                        %
% Special thanks to Andrei Hanu (McMaster University) for his MATLAB     %
% insight in the early stages of this endeavour.                         %
%************************************************************************%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   *********************** BRIEF INTRODUCTION ***********************   %
% This specific M-file was designed in order to compare 2 methods
% used for calculating the average dose rate to the cells during ICCM
% production (aka 'Accumulation Phase' or 'Dose Rate in Slab').
% Effectively, the program computes the dose rate in the slab using both
% the "homogeneous-slab" model (aka Method #1) or the "multi-isoplane"
% model (aka Method #2).  For Method #1, analytical integration is
% performed via utilization of either the MuPAD (Matlab) or Maple symbolic
% engine.  The goal is to compare both methods to determine if their
% answers agree and which Method provides optimal performance for
% integration into the program ICCM_Dose_Calcs_4_Radio_Beta_emitters(.m or
% .exe).
%   ************************ END INTRODUCTION ************************   %

% Prior to the program starting ...

% Clears Memory
clear all

% Clears Command Window
clc

% Maintains 15 digits of accuracy unless otherwise flagged
format long

% ====================================================================== %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~Welcome Dialogue for User of Start-up~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% The following sets up and displays a welcome dialogue box for the user.
% The user can select to continue to proceed to perform the calculations
% or the user may exit the program.
welcome_text = {'Hello!'
                '                                              '
                'This program is designed to calculate the dose rate during Accumulation Phase for both "homogenous-slab" (Method #1) and "multi-isoplane" (Method #2) models.  This stage is after uptake of \beta emitting radiopharmaceutical used in the creation of Irradiated Cell Conditioned Medium (ICCM) as per radiation bystander protocol utilized by Boyd et al in:'
                '                                              '
                'Boyd M, Ross SC, et al. Radiation-Induced biological bystander effect elicited in vitro by targeted radiopharmaceuticals labeled with alpa-, beta-, and auger electron-emitting radionuclides. J.Nucl.Med. Vol. 47. pp. 1007-1015 (2006).'
                '                                               '
                'The following prompts will guide the user through data entry necessary to perform these calculations leveraging the Vynicker-Wambersie (VW) kernel for an infinite plane.  Full details of the modelling utilized in this application can be found in the associated report.'
                '                                               '
                'Created by: Michael Gow'
                '                   Dept. of Medical Physics and Applied Radiation Sciences'
                '                   McMaster University'
                '                   Hamilton, ON CA'
                '                   MMXI'
                '                                               '
                '                                               '
                '                                               '
                '                                               '
                '                                               '};

welcome_title = 'Welcome!';
str1 = 'Continue';
str2 = 'Exit';
options.Resize='on';
options.Interpreter='tex';
options.Default=str1;
welcome = questdlg(welcome_text, welcome_title, str1, str2, options);

    % Determine if the user hits the 'Cancel' button
    if strcmp(welcome,'Exit')
        disp('Bye Bye')
        return
    end

% ====================================================================== %

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~Identify Symbolic Engine for MATLAB install and optimize if possible~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% For the "homogenous-slab" model, we will need to utilize the Symbolic
% toolbox to perform the analytical integration required.  This section
% allow the user to choose the symbolic engine they wish to use.

% The program will determine the version of the Symbolic toolbox
% installed on the users version of MATLAB. The Maple Kernel, the only
% symbolic kernel utilized in the symbolic toolbox version less than 4.9,
% provides superior performance to the MuPAD symbolic kernel introduced in
% Symbolic Toolbox 4.9.  Although improvements have been made in recent
% versions (5.4 and 5.5 with Matlab R2010a and R2010b), if given the
% choice, the user should utilize the Maple symbolic engine when available
% (i.e. if Maple is installed on the host operating system).  The ability
% to choose between the MuPAD and Maple kernels was given in Symbolic
% toolbox version 4.9 to 5.3 (inclusive) through initiation of the
% 'symengine' command. The following section looks to determine if the
% version of Symbolic toolbox installed has this ability to change symbolic
% engines thus allowing the user to choose the Maple kernel as the symbolic
% engine if available.

v = ver('symbolic');
%Determine the Symbolic Toolbox Version installed in MATLAB

v_num = str2double(v.Version);
% Convert the verision string in 'v' generated above into a number

    if (v_num <= 5.3) && (v_num >= 4.9)
        % Change symbolic engine if applicable
        symengine
        % If available, Maple's Symbolic engine is
        % preferred as it has superior computation performance over
        % MATLAB's native MuPAD engine.
    end

% NOTE: If the user hits cancel, the calculations will be automatically
% performed with the default engine for that MATLAB session.
% ====================================================================== %

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~Input Variables to Determine Thickness of Slab of Cells~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% Now let's have the user start entering values:

slab_mess_text = {'The cellular distribution within the treated flasks is approximated by a symmetrical slab of cells of some surface area A and thickness 2\Delta [e.g. cylindrical slab].'
                 '                                                        '
                 'With this approximation, we can determine the thickness of the cells, 2\Delta, when adhered to the flask through knowledge of:'
                 '                                                        '
                 '1) Size (surface area) of the flask used for culturing of donor cells'
                 '2) Confluence of the cell population at the time of radiopharmaceutical treatment'
                 '3) Diameter of the cell line utilized in suspension (i.e. free float)'
                 '4) Number of cells plated in the donor flasks'
                 '5) Doubling time of cell line used'
                 '6) Time between cell plating and radiopharmaceutical treatment'
                 '                                                        '
                 '                                                        '
                 'Additionally, parameters involved in measurement of the radiopharmaceutical uptake which can be noted include:'
                 '                                                        '
                 '7) Size (surface area) of the culture well / flask for uptake measurements'
                 '8) Number of cells plated in the culture well / flask for uptake measurements'
                 '9) Time between plating and radiopharmaceutical application for uptake measurements'
                 '                                                        '
                 'Note:  When entering fractional numbers (i.e. less than 1) make sure to include a leading '' 0. '' to prevent an error message.'
                 '                                                        '
                 '                                                        '
                 '                                                      '};
slab_mess_title = 'Determine the thickness of the cells in the treatment flask:';
str1_slab = 'Continue';
str2_slab = 'Exit';
options.Resize='on';
options.Interpreter='tex';
options.Default=str1_slab;
cell_slab_note = questdlg(slab_mess_text, slab_mess_title, str1_slab, str2_slab, options);

    % Determine if the user hits the 'Cancel' button
    if strcmp (cell_slab_note,'Exit')
        disp('Bye Bye')
        return
    end

% Now, let's have the user input the necessary variables to calculate the
% cell thickness when adhered to the bottom of the treatment flask:

prompt_slab = {'1) s_f_l_a_s_k = Size (growth area) of the flask used in treatment (in cm^2):'
               '2a) %_c_o_n_f = Percentage of cell confluence in treatment flask (between 0 and 1):'
               '2b) Err_%_c_o_n_f = Error / Uncertainty in cell confluence in treatment flask (between 0 and %_c_o_n_f ):'
               '3) Name of cell line used:'
               '4a) Cell_D = Cellular diameter while in suspension (in \mum):'
               '4b) Err_D = Error in cellular diameter while in suspension (in \mum):'
               '5) #_c_e_l_l_s = Number of cells plated. NOTE: No commas, spaces, or scientific/exponential notation:'
               '6) DB_t_i_m_e = Doubling time of cell line used (in hours):'
               '7) T_2_p_l_a_t_e = Time between cell plating and ^1^3^1I treatment (in hours):'
               '8) s_w_e_l_l = Size (growth area) of the well / flask used in uptake measurements (in cm^2).  NOTE:  If a 6 - well plate was utilized, the typical surface area is 9.6 cm^2:'
               '9) #_u_p_t_a_k_e = Number of cells plated in well / flask in uptake measurements. NOTE: No commas, spaces, or scientific/exponential notation:'
               '10) T_2_u_p_t_a_k_e = Time between plating and application of radiopharmaceutical in uptake measurements (in hours):'};

title_slab = 'Input Parameters for determining cell thickness:';
num_lines_slab = [1 90]; % 1 row per prompt 100 characters wide
options.Resize='on';
options.Interpreter='tex';

% Define the standard parameters
def_slab = {'25.0','0.65','0.05','UVW/NAT','16.5','0.9','200000','18.0','24.0','9.6','50000','48'};

% 'answer_slab' is a four element vector containing the users input results
answer_slab = inputdlg(prompt_slab,title_slab,num_lines_slab,def_slab,options);

    % Determine if the user hits the 'Cancel' button
    if isempty(answer_slab) == 1
        disp('Bye Bye')
        return
    end

    % Checks to ensure that all input parameters have a value
    while isempty(answer_slab{1}) || isempty(answer_slab{2}) || isempty(answer_slab{3}) || isempty(answer_slab{4}) || isempty(answer_slab{5}) || isempty(answer_slab{6}) || isempty(answer_slab{7}) || isempty(answer_slab{8}) || isempty(answer_slab{9}) || isempty(answer_slab{10}) || isempty(answer_slab{11}) || isempty(answer_slab{12})
        empty_err_slab = errordlg({'You have missed entering information into one of the input parameters.'
                                   'All parameters require a value.  Please Try Again.'},'Error');
        uiwait(empty_err_slab)
        answer_slab = inputdlg(prompt_slab,title_slab,num_lines_slab,def_slab,options);
            if isempty(answer_slab) == 1
                disp('Bye Bye')
                return
            end
    end

% Now, we convert the strings entered by the user in to double precision
% floating point numbers.

% Size of the flask
S_flask = str2double(answer_slab{1});
% Percentage confluence
Per_conf = str2double(answer_slab{2});
% Error in Percentage Confluence
Err_per_conf = str2double(answer_slab{3});
% Cell diameter
Cell_diam = str2double(answer_slab{5});
% Error in cell diameter
Err_cell_diam = str2double(answer_slab{6});
% # of Cells Plated
Num_cells = str2double(answer_slab{7});
% Doubling time of cell line used
Db_time = str2double(answer_slab{8});
% Time between plating and treatment
Pl2treat_time = str2double(answer_slab{9});
% Size of well used for uptake experiment
S_well = str2double(answer_slab{10});
% # of cells plated for uptake measurement
Num_uptake = str2double(answer_slab{11});
% Time between plating and uptake experiments
Pl2treat_up_time = str2double(answer_slab{12});

% Also for convenience, lets assign a variable to the name of the cells
Cell_name = answer_slab{4};

    % Check to ensure that the user has entered valid numbers (e.g. not
    % negative).  ISNAN checks to ensure a number was entered as opposed to
    % a alpha or special character.
    while (S_flask <= 0) || isnan(S_flask) == 1 || (Per_conf <= 0) || (Per_conf > 1) || isnan(Per_conf) == 1 || (Err_per_conf < 0) || (isnan(Err_per_conf) == 1) || (Err_per_conf > Per_conf) || (Cell_diam <= 0) || (isnan(Cell_diam) == 1) || (Err_cell_diam < 0) || (isnan(Err_cell_diam) == 1) || (Err_cell_diam > Cell_diam) || (Num_cells <= 0) || (isnan(Num_cells) == 1) || (Db_time <= 0) || (isnan(Db_time) == 1) || (Pl2treat_time <= 0) || (isnan(Pl2treat_time) == 1) || (S_well <= 0) || (isnan(S_well) == 1) || (Num_uptake <= 0) || (isnan(Num_uptake) == 1) || (Pl2treat_up_time <= 0) || (isnan(Pl2treat_up_time) == 1)
        err_slab_neg = errordlg({'This error has resulted from one of the following:'
                                 '---------------------------------------------------------------------------'
                                 '                                        '
                                 'a) You have entered a negative value or a value of zero for one of the following parameters:'
                                 '-> size of flask, cell diameter, # of cells, doubling time, and/or time between plating and radiopharmaceutical application'
                                 '                                        '
                                 'b) The confluence entered is outside the acceptable range.'
                                 '                                        '
                                 'c) The error/uncertainty in the confluence is greater than the confluence value entered.'
                                 '                                        '
                                 'd) The error/uncertainty in the cell diameter is greater than the cell diameter value entered.'
                                 '                                        '
                                 'e) You have not entered valid numeric characters (with the exception of the cell line name).  Not entering a leading '' 0. '' for numbers less than 1 may cause this error.'
                                 '                                        '
                                 '                                        '
                                 'Please ensure all parameters are positive and have been entered correctly.'
                                 '                                        '
                                 '                                        '
                                 '                             '},'Error');
        uiwait(err_slab_neg)
        answer_slab = inputdlg(prompt_slab,title_slab,num_lines_slab,def_slab,options);

            % If the user enters 'Cancel'...
            if isempty(answer_slab) == 1
                disp('Bye Bye')
                return
            end

        % Convert to double floating point precision
        S_flask = str2double(answer_slab{1});
       	Per_conf = str2double(answer_slab{2});
        Err_per_conf = str2double(answer_slab{3});
        Cell_diam = str2double(answer_slab{5});
        Err_cell_diam = str2double(answer_slab{6});
        Num_cells = str2double(answer_slab{7});
        Db_time = str2double(answer_slab{8});
        Pl2treat_time = str2double(answer_slab{9});
        S_well = str2double(answer_slab{10});
        Num_uptake = str2double(answer_slab{11});
        Pl2treat_up_time = str2double(answer_slab{12});
    end

% Now, we calculate the cell thickness upon adhesion to the donor flasks
% as follows:

% *******
% First, calculate the volume of one cell in suspension.
% This volume is conserved upon adhesion to the flask:
% NOTE: We multiply by 1E-4 to convert micrometres to centimetres

Cell_vol = (4/3)*pi*(0.5*Cell_diam*0.0001)^3;

% Now, calculate the uncertainty in the cell volume in suspension.
% Remember, if we have Value +/- Error, the error (standard deviation) is
% calculated via: Error = Value*(Square Root of the Sum of the Squares) of
% all the relative errors (i.e. relative error = Error/Value).  So here, if
% we have a power of 3, we calculate the error as:
% (Value^3)*((3*(Err/Value)^2)^0.5)

Err_Cell_vol = ((4/3)*pi)*((0.5*Cell_diam*0.0001)^3*((3*(Err_cell_diam/Cell_diam)^2)^0.5));

% *******
% Second, determine surface area of the flask taken up by the cells.  For
% example, if a 25 cm^2 flask is 65% confluent, the total surface area
% taken up by the cells is 25*0.65 = 16.25 cm^2

Cell_surf_area = S_flask*Per_conf;

% Uncertainty in surface area

Err_Cell_surf_area = S_flask*Err_per_conf;

% *******
% Third, we need to determine how many cells make up this population
% covering surface area 'Cell_surf_area' above.  To do this, we need to
% take into account the doubling time of the cell line so that the total
% number cells at given time, t, between plating and treatment is as
% follows: C = Co*exp[(t*ln2)/DB] where DB is the doubling time of the cell
% line and Co is the number of cells originally plated.

Tot_num_cells = Num_cells*exp((Pl2treat_time*log(2))/Db_time);

% *******
% Fourth, so the total volume taken up by the cells at time of treatment
% is:

Tot_cell_vol = Cell_vol*Tot_num_cells;

% Uncertainty in total volume

Err_Tot_cell_vol = Err_Cell_vol*Tot_num_cells;

% *******
% Fifth, and final step is to calculate the thickness of the cells upon
% adhering to the flask.  We know the total volume of the cells (Step 4) as
% well as the surface area these cells covered in the flask (Step 2).
% Since Volume = Surface Area*Thickness, we have simply, Cell Slab
% Thickness = 2Delta = (Tot_cell_vol)/(Cell_surf_area):

Cell_slab_thick = Tot_cell_vol/Cell_surf_area;
Cell_slab_half_thick = 0.5*Cell_slab_thick;

% Uncertainty in slab thickness

Err_Cell_slab_thick = (Tot_cell_vol/Cell_surf_area)*((Err_Tot_cell_vol/Tot_cell_vol)^2 + (Err_Cell_surf_area/Cell_surf_area)^2)^0.5;
Err_Cell_slab_half_thick = 0.5*Err_Cell_slab_thick;
% *******

% For usage later, lets also calculate the maximum and minimum slab
% thickness
Cell_slab_thick_max = Cell_slab_thick + Err_Cell_slab_thick;
Cell_slab_thick_min = Cell_slab_thick - Err_Cell_slab_thick;
Cell_slab_half_thick_max = 0.5*Cell_slab_thick_max;
Cell_slab_half_thick_min = 0.5*Cell_slab_thick_min;

% *******
% Using the thickness calculated above, coupled with input variables
% provided for the uptake measurement experiments, we will determine the
% confluence of the wells / flask used for determining our
% radiopharmaceutical uptake percentages (to be entered by the user in the
% next section.  Remember, upon application of the radiopharmaceutical to
% the medium, we are assuming an instant, homogenous, static mix.  Because
% we assume it to be static, areas covered with radiopharmaceutical but no
% cells will not have the pharmaceutical taken up.  Thus, we can make a
% correlation between activity administered, well/flask confluence, and
% uptake percentage from our uptake experiment data which will allow us to
% more accurately determine the the amount of activity taken up by donor
% cells used for ICCM production. We will call this a confluence adjustment
% factor, Conf_adjust_uptake, obtained for an applied activity
% concentration, Aconc, as:
%
% [(% Uptake from exp.)/(Uptake Conf.)] =
%                                [(X actual uptake by donors)/(Donor Conf.]
%
% where (X actual uptake by donors) = Conf_adjust_uptake
%
% Typically, the wells in the uptake experiments are
% seeded and allowed to culture for long enough such that they are fully
% confluent at the time of uptake application and subsequent measurement
% (i.e. Uptake Conf. = 1).  In this case, the confluence adjustment is
% straight forward.  Simply multiple the uptake percentage by the donor
% cell confluence.  Intuitively, assuming a static, homogenous
% radiopharmaceutical distribution / coverage, this is obvious (i.e. cells
% will only uptake the activity that is covering them).  Otherwise the
% expression above will make the appropriate adjustment.
%
% Conf_adjust_uptake is calculated as follows:

% 1) Determine the number of cells after the time between seeding and start
% of uptake experiment (i.e. radiopharmaceutical application).

Tot_num_cells_uptake = Num_uptake*exp((Pl2treat_up_time*log(2))/Db_time);

% 2) Use the cell volume and donor cell thickness calculated above to
% determine the area covered in the well.

Conf_uptake = (Tot_num_cells_uptake*Cell_vol)/Cell_slab_thick;
Err_conf_uptake = Conf_uptake*(Err_Cell_slab_thick);

% We will identify Conf_adjust_uptake below after we define some additional
% variables....keep your eyes open ;)

% ====================================================================== %

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~Input variables to be used in Vynckier Wambersie Kernel~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% A note to the user about what is left to enter.
addn_param_text = {'Additional experimental parameters regarding the radiopharmaceutical and medium utilized required for the calculations are:'
                   '                                                     '
                   '1) The mean energy of \beta particles for the nuclide utilized'
                   '2) The maximum energy of \beta particles for the nuclide utilized'
                   '3) Density of the homogenous medium'
                   '4) The activity concentration applied'
                   '5) Percentage uptake of the radiopharmaceutical into the cell in vitro'
                   '6) Time for uptake phase'
                   '7) Time for bystander factor accumulation'
                   '                                                     '
                   'Note:  When entering fractional numbers (i.e. less than 1) make sure to include a leading '' 0. '' to prevent an error message.'
                   '                                                     '
                   '                                                    '};

addn_param_title = 'Additional Experimental Parameters';
str1_addn_param = 'Continue';
str2_addn_param = 'Exit';
options.Resize='on';
options.Interpreter='tex';
options.Default=str1_addn_param;
addn_param_note = questdlg(addn_param_text,addn_param_title,str1_addn_param,str2_addn_param,options);

    if strcmp (addn_param_note,'Exit')
        disp('Bye Bye')
        return
    end

% Now, lets have the user input the necessary variables
addn_param_input = {'<E_\beta> = Mean \beta energy per disintegration (in MeV):'
                    'E_\beta_m_a_x = Maximum \beta energy (in MeV). NOTE: Here, this value is limited to the range of 0.5 MeV to 3.0 MeV (exclusive).  The reason for this due to E_\beta_m_a_x''s use in evaluating additional variables in the VW kernel (see Vynckier and Wambersie, 1982, 1986):'
                    '\rho = Density of the homogeneous medium (in g/cm^3; 1 for standard medium):'
                    'A_c_o_n_c = Activity concentration of radiopharmaceutical applied (in MBq/cm^3). NOTE: With medium of unit density, then 1 MBq/ml = 1 MBq/cm^3 = 1 MBq/g:'
                    '%_u_p_t_a_k_e = Percentage uptake of radiopharmaceutical (between 0 and 1):'
                    'Err_%_u_p_t_a_k_e = Error / Uncertainty in uptake of radiopharmaceutical (between 0 and %_u_p_t_a_k_e):'};

addn_param_input_title = 'Additional Experimental Parameters';
addn_param_num_lines = [1 90];
options.Resize='on';
options.Interpreter='tex';

% Define the standard parameters
def_addn_param = {'0.19','0.61','1.00','1.00','0.314','0.037'};
addn_param_answer = inputdlg(addn_param_input,addn_param_input_title,addn_param_num_lines,def_addn_param,options);

    % Determine if the user hits the 'Cancel' button and if so displays
    % an exit prompt
    if isempty(addn_param_answer) == 1
        disp('Bye Bye')
        return
    end

    % Checks to ensure that all input parameters have a value
    while isempty(addn_param_answer{1}) || isempty(addn_param_answer{2}) || isempty(addn_param_answer{3}) || isempty(addn_param_answer{4}) || isempty(addn_param_answer{5}) || isempty(addn_param_answer{6})
        addn_param_empty_err = errordlg({'You have missed entering information into one of the input parameters.'
                                         'All parameters require a value (even 0).  Please Try Again.'},'Error');
        uiwait(addn_param_empty_err)
        addn_param_answer = inputdlg(addn_param_input,addn_param_input_title,addn_param_num_lines,def_addn_param,options);
            if isempty(addn_param_answer) == 1
                disp('Bye Bye')
                return
            end
    end

% Now, we convert the strings entered by the user in to double precision
% floating point numbers.

% Mean beta energy per disintegration (in Mev)
Mean_beta_eng = str2double(addn_param_answer{1});
% Maximum beta energy (in MeV)
Max_beta_eng = str2double(addn_param_answer{2});
% Density of medium (g/cm^3)
Dens_med = str2double(addn_param_answer{3});
% A_conc applied (in MBq/ml)
Act_conc = str2double(addn_param_answer{4});
% % of radiopharmaceutical uptake
Percent_up = str2double(addn_param_answer{5});
% Error / Uncertainty in % of radiopharmaceutical uptake
Err_percent_up = str2double(addn_param_answer{6});

    % Check to ensure that the user has not entered a negative number or
    % that the confluence level is not below the acceptable value.
    while (Mean_beta_eng <= 0) || (isnan(Mean_beta_eng) == 1) || (isnan(Max_beta_eng) == 1) || (Max_beta_eng <= 0.5) || (Max_beta_eng >= 3.0) || (Dens_med <= 0) || (isnan(Dens_med) == 1) || (Act_conc <= 0) || (isnan(Act_conc) == 1)|| (Percent_up <= 0) || (isnan(Percent_up) == 1) || (Percent_up > 1) || (Err_percent_up <= 0) || (isnan(Err_percent_up) == 1) || (Err_percent_up > Percent_up)
        err_slab_neg = errordlg({'This error has resulted from one of the following:'
                                 '---------------------------------------------------------------------------'
                                 '                                                  '
                                 'a) You have entered a negative value or a value of 0 for one of the following parameters:'
                                 '-> Average beta energy, density of the medium, percentage uptake, activity concentration, time allotted for uptake, and/or time allotted for bystander factor accumulation'
                                 '                                                  '
                                 'b) The maximum beta energy entered is outside the accepted range of 0.5 to 3.0 MeV (exclusive).'
                                 '                                                  '
                                 'e) The percentage of uptake is greater than 100% (i.e. greater than 1) or the error in uptake is larger than or equal to the entered uptake value.'
                                 '                                                  '
                                 'f) You have entered a non-numeric character.  Not entering a leading '' 0. '' for numbers less than 1 may cause this error.'
                                 '                                                  '
                                 '                                                  '
                                 'Please ensure all parameters are positive and have been entered correctly.'
                                 '                                                  '
                                 '                                                  '
                                 '                                                  '},'Error');

        uiwait(err_slab_neg)
        addn_param_answer = inputdlg(addn_param_input,addn_param_input_title,addn_param_num_lines,def_addn_param,options);

            % If the user enters 'Cancel'...
            if isempty(addn_param_answer) == 1
                disp('Bye Bye')
                return
            end

        % Convert to double floating point precision
        Mean_beta_eng = str2double(addn_param_answer{1});
        Max_beta_eng = str2double(addn_param_answer{2});
        Dens_med = str2double(addn_param_answer{3});
        Act_conc = str2double(addn_param_answer{4});
        Percent_up = str2double(addn_param_answer{5});
        Err_percent_up = str2double(addn_param_answer{6});

    end

% Now, we have Conf_adjust_uptake as follows:

if (Conf_uptake-Err_conf_uptake) >= S_well
    Conf_uptake = 1;
    %If the minimum area covered by the cells is greater than or equal to
    %the size of the well then we have 100% confluence at time of uptake
    %experiment and Conf_adjust_uptake is simply the percentage confluence:
    Conf_adjust_uptake = Per_conf;
    Err_conf_adjust_uptake = Err_per_conf;
else
    Conf_adjust_uptake = (Per_conf*Percent_up)/Conf_uptake;
    Err_conf_adjust_uptake = Conf_adjust_uptake*(sqrt((Err_per_conf/Per_conf)^2 + (Err_percent_up/Percent_up)^2));
end


% *******
% Now, let's calculate some of the necessary variables that form part of
% the VW kernel
% *******

% 1) Apparent absorption coefficient, v [cm^2/g -> density thickness].
% This relation is accurate provided the maximum beta energy of the
% radionuclide is between 0.5 and 3.5 MeV.

v = 14.5*(Max_beta_eng)^(-1.17);

% NOTE: By comparison, the original Loevinger kernel, calculated the
% apparent absorption coefficient, v, via: v = 18.6*(Emax - 0.036)^(-1.37)
% for the range 0.17 to 3 MeV *******

% 2) Dimensionless parameter c.  Remember we are only looking in the range
% of 0.5 to 3.0 MeV for the maximum beta energy

if (Max_beta_eng >= 0.5 && Max_beta_eng < 1.5)
    c = 1.5;
elseif (Max_beta_eng >= 1.5 && Max_beta_eng <= 3.0)
    c = 1;
end
% *******

% 3) Dimensionless f parameter. The parameter f/pv is the distance where
% the point kernel becomes zero.  This is the primary correction entered by
% VW into Loevinger's work. This relationship again holds true in the
% energy range between 0.5 and 3.5 MeV.

f = Dens_med*v*0.269*(Max_beta_eng)^1.31;
% *******

% 4) Alpha Coefficient. For explanation see Ref. 1-3 in "Brief
% Introduction".

alpha = ((3*c^2) - (c^2-1)*exp(1) + (3+f)*exp(1-f) - 4*exp(1-(f/2)))^-1;
% *******

% ====================================================================== %

fprintf(1,'______________________________________________________\n')
fprintf(1,'======================================================\n\n\n')
fprintf(1,'For an activity concentration of: %4.2f MBq/ml\n',Act_conc)
fprintf(1,'---------------------------------\n\n')

% ====================================================================== %
tStart_1=tic; % START OF TIMER FOR METHOD 1

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~Calculate Dose Rate in Slab~% -> METHOD 1; homogenous slab
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% Briefly, we are dividing the cell into a bunch of infinitely thin plane
% sources located between x'=0 and 2delta.  For a position, x, the dose
% rate would be:
%
% Drate(x)=Aconc*int(K(x-x')) from x'=0 to 2delta and then take the average
% over all x.
%
% where K(x-x') is:
%
% K(x-x') =
% 0.288*Eavg*Ac*v*alpha*{c[1+ln(c/pv(x-x'))-exp(1-(pv(x-x'))/c)]+exp(1
% -pv(x-x'))-2*exp(1-pv(x-x')/c-f/2)+exp(1-f)}
%
% i.e. The dose rate kernel for an infinite, thin plane.
%
% As the dose rate is really a function of the magnitude of x-x', |x-x'|,
% we make the appropriate substitution of variables and work only the
% positive half of the kernel.  Effectively we get that:
%
% K(x-x') = K_low + K_high where K_low is K for x-x'<0 and K_high is K for
% x-x'>0.
%
% See the associated paper for full details.

syms x z real % IDENTIFY SYMBOLIC VARIABLES

K=(0.288*Mean_beta_eng*v*alpha)*(c+(c*log(c/(Dens_med*v*z)))-c*exp(1-((Dens_med*v*z)/c))+exp(1-(Dens_med*v*z))-2*exp(1-((Dens_med*v*z)/2)-(f/2))+exp(1-f));

    % Test to adjust K if c<= p*v*x.  We will set x=2*Delta to be
    % stringent.
    if c<=Dens_med*v*Cell_slab_thick
        K=(0.288*Mean_beta_eng*v*alpha)*(exp(1-(Dens_med*v*z))-2*exp(1-((Dens_med*v*z)/2)-(f/2))+exp(1-f));
    else
        if f<=Dens_med*v*Cell_slab_thick % Test for f to adjust K
            fprintf(1,'Dose rate is 0 Gy/h.')
        end
    end

% x-x' < 0
K_low=int(K,z,0,(Cell_slab_thick-x));
% x-x' > 0
K_high=int(K,z,0,x);

Drate_slab_avg_1_prime = ((Act_conc*Percent_up*Conf_adjust_uptake)/(S_flask*Per_conf*Cell_slab_thick))*(K_low+K_high);
Drate_slab_avg_1 = (1/Cell_slab_thick)*int(Drate_slab_avg_1_prime,x,0,Cell_slab_thick);

% Maximum
K_low=int(K,z,0,((Cell_slab_thick-Err_Cell_slab_thick)-x));
Drate_slab_avg_1_max_prime = ((Act_conc*(Percent_up+Err_percent_up)*(Conf_adjust_uptake+Err_conf_adjust_uptake))/(S_flask*(Per_conf-Err_per_conf)*(Cell_slab_thick-Err_Cell_slab_thick)))*(K_low+K_high);
Drate_slab_avg_1_max = (1/(Cell_slab_thick-Err_Cell_slab_thick))*int(Drate_slab_avg_1_max_prime,x,0,(Cell_slab_thick-Err_Cell_slab_thick));

% Minimum
K_low=int(K,z,0,((Cell_slab_thick+Err_Cell_slab_thick)-x));
Drate_slab_avg_1_min_prime = ((Act_conc*(Percent_up-Err_percent_up)*(Conf_adjust_uptake-Err_conf_adjust_uptake))/(S_flask*(Per_conf+Err_per_conf)*(Cell_slab_thick+Err_Cell_slab_thick)))*(K_low+K_high);
Drate_slab_avg_1_min = (1/(Cell_slab_thick+Err_Cell_slab_thick))*int(Drate_slab_avg_1_min_prime,x,0,(Cell_slab_thick+Err_Cell_slab_thick));

% Uncertainty
Drate_slab_avg_plus_err_1 = Drate_slab_avg_1_max - Drate_slab_avg_1;
Drate_slab_avg_minus_err_1 = Drate_slab_avg_1 - Drate_slab_avg_1_min;

tElapsed_1=toc(tStart_1); % END OF TIMER FOR METHOD 2

fprintf(1,'METHOD #1 ("homogenous slab"):\n\n')
fprintf(1,'The average dose rate to the slab of %s cells\n',Cell_name)
fprintf(1,'after uptake of the radiopharmaceutical is:\n** %d ', double(Drate_slab_avg_1))
fprintf(1,' +/- %d/%d Gy/h**\n\n', double(Drate_slab_avg_plus_err_1),double(Drate_slab_avg_minus_err_1))
fprintf(1,'Total time taken = %.5f sec\n\n', tElapsed_1)
fprintf(1,'~~~~~~~\n\n')

% ====================================================================== %
tStart_2=tic; % START OF TIMER FOR METHOD 2

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~Calculate Dose Rate in Slab~% -> METHOD 2; multi-isoplane
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% We will create two different activity distribution models (aka paradigms)
% and calculate the dose rate at 5 different locals within the slab of
% cells.  In general, the slab is located infinite in the y-, z- extents
% and has a thickness of 2delta along the x axis.  The bottom of the slab
% lies along x=0.

% 1)  Take the activity, Aconc, and distribute it equally in two infinite
% planes (i.e. 0.5*Aconc each) located at positions x=(0.5)delta and
% x=(1.5)delta.  Calculate the  dose rate (Drate) at the edges of the slab
% (x=0 and x=2delta) and in the centre (x=delta)
%
% 2)  Take the activity, Aconc, and distribute it equally in three infinite
% planes (i.e. (1/3)Aconc each) located at positions x=0,delta, and
% 2delta).  Calculate the Drate at x=(0.5)delta and x=(1.5)delta.
%
% The average dose rate to the slab will be derived from the values
% calculated at the various locations in paradigms 1) and 2).

% The expression for dose rate from an infinite disk at distance x from
% source is given by:
%
% Drate(x,0,inf)=0.288*Eavg*Ac*v*alpha*{c[1+ln(c/pvx)-exp(1-(pvx)/c)]+exp(1
% -pvx)-2*exp(1-pvx/c-f/2)+exp(1-f)}
%
% with [  ] defined as 0 for pvx >= c and Drate(x,0,inf) defined as 0 for
% pvx >= f [Note: For definitions of v,c,alpha,and f, see below.]
%
% References:
% [1] Vynckier S and Wambersie A. Dosimetry of beta sources in
% radiotherapy: Absorbed dose distributions around plane sources. Rad.
% Prot. Dos. Vol 14(2), pp. 169-173 (1986).
%
% [2] Appendix C: Calculation of beta-ray dose distributions by integration
% of the beta-ray point-source dose function. Oxford University Press. J.
% ICRU. Vol. 4(2), pp. 155-163 (2004).

% Note that we have used scalar .* and ./ here so as to ensure a more
% object-oriented design. Remember that the activity in the infinite plane
% is only the appropriate fraction of the percentage of uptake.

Drate_x=@(x)((Act_conc.*Percent_up.*Conf_adjust_uptake)/(S_flask*Per_conf)).*0.288.*Mean_beta_eng.*v.*alpha.*(c+(c.*log(c./(Dens_med.*v.*x)))-(c.*exp(1-((Dens_med.*v.*x)./c)))+exp(1-(Dens_med.*v.*x))-2.*exp(1-((Dens_med.*v.*x)./2)-(f./2))+exp(1-f));

% We need the appropriate 'if' statements for the conditions on an
% infinite thin plan VW kernel.  If c>=p*v*x, the quantity in the square
% brackets of Eq. [9] becomes 0 and the DPK changes to:

Drate_x_c_cond=@(x)((Act_conc.*Percent_up.*Conf_adjust_uptake)/(S_flask*Per_conf)).*0.288.*Mean_beta_eng.*v.*alpha.*(exp(1-(Dens_med.*v.*x))-2.*exp(1-((Dens_med.*v.*x)./2)-(f./2))+exp(1-f));

% Checks on this condition will be done on each evaluation below.

% Paradigm 1):

% Bottom of the slab:
Drate_bottom = (0.5*Drate_x((0.5)*Cell_slab_half_thick)) + (0.5*Drate_x((1.5)*Cell_slab_half_thick));
    if (c<=(Dens_med*v*((0.5)*Cell_slab_half_thick)))
        % the smallest condition
        Drate_bottom = (0.5*Drate_x_c_cond((0.5)*Cell_slab_half_thick)) + (0.5*Drate_x_c_cond((1.5)*Cell_slab_half_thick));
    else
        if (f<=(Dens_med*v*((0.5)*Cell_slab_half_thick)))
            % the smallest condition
            fprintf(1,'The dose rate is 0 Gy/h. No need to go any further...\n')
        end
    end

% Top of the slab (by symmetry):
Drate_top = Drate_bottom;

% Centre of the slab:
Drate_centre = 2*(0.5*Drate_x((0.5)*Cell_slab_half_thick));
    if (c<=(Dens_med*v*((0.5)*Cell_slab_half_thick)))
        Drate_centre = 2*(0.5*Drate_x_c_cond((0.5)*Cell_slab_half_thick));
    else
        if (f<=(Dens_med*v*((1/3)*Cell_slab_half_thick)))
            fprintf(1,'The dose rate is 0 Gy/h. No need to go any further...\n')
        end
    end

% Couple notes to remember is that delta = Cell_slab_half_thick and we are
% multiplying by 0.5 as each infinite plane of activity contains half the
% overall activity taken up by the cells.

% Paradigm 2):

% Where x=(0.5)delta
Drate_delta_one_half = 2*((1/3)*Drate_x((0.5)*Cell_slab_half_thick))+((1/3)*Drate_x((1.5)*Cell_slab_half_thick));
    if (c<=(Dens_med*v*((0.5)*Cell_slab_half_thick)))
        % the smallest condition
        Drate_one_third = 2*((1/3)*Drate_x_c_cond((0.5)*Cell_slab_half_thick))+((1/3)*Drate_x_c_cond((1.5)*Cell_slab_half_thick));
    else
        if (f<=(Dens_med*v*((0.5)*Cell_slab_half_thick)))
            % the smallest condition
            fprintf(1,'The dose rate is 0 Gy/h. No need to go any further...\n')
        end
    end

% Where x=(1.5)delta, by symmetry
Drate_delta_one_and_one_half = Drate_delta_one_half;

% Total dose rate in the slab is:

Drate_slab_avg_2 = (Drate_bottom+Drate_top+Drate_centre+Drate_delta_one_half+Drate_delta_one_and_one_half)/5;

    % Now, we need to determine the +/- uncertainties.  Maximum Drate (i.e.
    % +) occurs when Aconc is highest, slab thickness is smallest, and
    % confluence is greatest.  Conversely, minimum Drate occurs when Aconc
    % is lowest and slab thickness is largest and confluence is lowest.
    % Therefore...

    % Maximum Drate:
    Drate_x_max=@(x)((Act_conc.*(Percent_up+Err_percent_up).*(Conf_adjust_uptake+Err_conf_adjust_uptake))/(S_flask*(Per_conf-Err_per_conf))).*0.288.*Mean_beta_eng.*v.*alpha.*(c+(c.*log(c./(Dens_med.*v.*x)))-(c.*exp(1-((Dens_med.*v.*x)./c)))+exp(1-(Dens_med.*v.*x))-2.*exp(1-((Dens_med.*v.*x)./2)-(f./2))+exp(1-f));
    %With c condition...
    Drate_x_max_c_cond=@(x)((Act_conc.*(Percent_up+Err_percent_up).*(Conf_adjust_uptake+Err_conf_adjust_uptake))/(S_flask*(Per_conf-Err_per_conf))).*0.288.*Mean_beta_eng.*v.*alpha.*(exp(1-(Dens_med.*v.*x))-2.*exp(1-((Dens_med.*v.*x)./2)-(f./2))+exp(1-f));

    % Paradigm 1):

    % Bottom of the slab:
    Drate_bottom_max = (0.5*Drate_x_max((0.5)*(Cell_slab_half_thick_min))) + (0.5*Drate_x_max((1.5)*(Cell_slab_half_thick_min)));
        if (c<=(Dens_med*v*((0.5)*Cell_slab_half_thick_min)))
            Drate_bottom_max = (0.5*Drate_x_max_c_cond((0.5)*(Cell_slab_half_thick_min))) + (0.5*Drate_x_max_c_cond((1.5)*(Cell_slab_half_thick_min)));
        else
            if (f<=(Dens_med*v*((0.5)*Cell_slab_half_thick_min)))
                fprintf(1,'The dose rate is 0 Gy/h. No need to go any further...\n')
            end
        end

    % Top of the slab (by symmetry):
    Drate_top_max = Drate_bottom_max;

    % Centre of the slab:
    Drate_centre_max = 2*(0.5*Drate_x_max((0.5)*(Cell_slab_half_thick_min)));
        if (c<=(Dens_med*v*((0.5)*Cell_slab_half_thick_min)))
            Drate_centre_max = 2*(0.5*Drate_x_max_c_cond((0.5)*(Cell_slab_half_thick_min)));
        else
            if (f<=(Dens_med*v*((0.5)*Cell_slab_half_thick_min)))
                frintf(1,'The dose rate is 0 Gy/h. No need to go any further...\n')
            end
        end

    % Paradigm 2):

    % Where x=(0.5)delta
    Drate_delta_one_half_max = 2*((1/3)*Drate_x_max((0.5)*(Cell_slab_half_thick_min)))+((1/3)*Drate_x_max((1.5)*(Cell_slab_half_thick_min)));
        if (c<=(Dens_med*v*((1/3)*Cell_slab_half_thick)))
            Drate_delta_one_half_max = 2*((1/3)*Drate_x_max_c_cond((0.5)*(Cell_slab_half_thick_min)))+((1/3)*Drate_x_max_c_cond((1.5)*(Cell_slab_half_thick_min)));
        else
            if (f<=(Dens_med*v*((1/3)*Cell_slab_half_thick)))
                fprintf(1,'The dose rate is 0 Gy/h. No need to go any further...\n')
            end
        end

    % Where x=(1.5)delta, by symmetry
    Drate_delta_one_and_one_half_max = Drate_delta_one_half_max;

    % Total dose rate in the slab is:

    Drate_slab_avg_max = (Drate_bottom_max+Drate_top_max+Drate_centre_max+Drate_delta_one_half_max+Drate_delta_one_and_one_half_max)/5;

    % Minimum Drate:

    Drate_x_min=@(x)((Act_conc.*(Percent_up-Err_percent_up).*(Conf_adjust_uptake-Err_conf_adjust_uptake))/(S_flask*(Per_conf+Err_per_conf))).*0.288.*Mean_beta_eng.*v.*alpha.*(c+(c.*log(c./(Dens_med.*v.*x)))-(c.*exp(1-((Dens_med.*v.*x)./c)))+exp(1-(Dens_med.*v.*x))-2.*exp(1-((Dens_med.*v.*x)./2)-(f./2))+exp(1-f));
    % With c condition...
    Drate_x_min_c_cond=@(x)((Act_conc.*(Percent_up-Err_percent_up).*(Conf_adjust_uptake-Err_conf_adjust_uptake))/(S_flask*(Per_conf+Err_per_conf))).*0.288.*Mean_beta_eng.*v.*alpha.*(exp(1-(Dens_med.*v.*x))-2.*exp(1-((Dens_med.*v.*x)./2)-(f./2))+exp(1-f));

    % Paradigm 1):

    % Bottom of the slab:
    Drate_bottom_min = (0.5*Drate_x_min((0.5)*(Cell_slab_half_thick_max))) + (0.5*Drate_x_min((1.5)*(Cell_slab_half_thick_max)));
        if (c<=(Dens_med*v*((0.5)*Cell_slab_half_thick_max)))
            Drate_bottom_min = (0.5*Drate_x_min_c_cond((0.5)*(Cell_slab_half_thick_max))) + (0.5*Drate_x_min_c_cond((1.5)*(Cell_slab_half_thick_max)));
        else
            if (f<=(Dens_med*v*((0.5)*Cell_slab_half_thick_max)))
                fprintf(1,'The dose rate is 0 Gy/h. No need to go any further...\n')
            end
        end

    % Top of the slab (by symmetry):
    Drate_top_min = Drate_bottom_min;

    % Centre of the slab:
    Drate_centre_min = 2*(0.5*Drate_x_min((0.5)*(Cell_slab_half_thick_max)));
        if (c<=(Dens_med*v*((0.5)*Cell_slab_half_thick_max)))
            Drate_centre_min = 2*(0.5*Drate_x_min_c_cond((0.5)*(Cell_slab_half_thick_max)));
        else
            if (f<=(Dens_med*v*((0.5)*Cell_slab_half_thick_max)))
                fprintf(1,'The dose rate is 0 Gy/h. No need to go any further...\n')
            end
        end

    % Paradigm 2):

    % Where x=(0.5)delta
    Drate_delta_one_half_min = 2*((1/3)*Drate_x_min((0.5)*(Cell_slab_half_thick_max)))+((1/3)*Drate_x_min((1.5)*(Cell_slab_half_thick_max)));
        if (c<=(Dens_med*v*((0.5)*Cell_slab_half_thick_max)))
            Drate_one_third_min = 2*((1/3)*Drate_x_min_c_cond((0.5)*(Cell_slab_half_thick_max)))+((1/3)*Drate_x_min_c_cond((1.5)*(Cell_slab_half_thick_max)));
        else
            if (f<=(Dens_med*v*((0.5)*Cell_slab_half_thick_max)))
                fprintf(1,'The dose rate is 0 Gy/h. No need to go any further...\n')
            end
        end

    % Where x=(4/3)delta, by symmetry
    Drate_delta_one_and_one_half_min = Drate_delta_one_half_min;

    % Total dose rate in the slab is:

    Drate_slab_avg_min = (Drate_bottom_min+Drate_top_min+Drate_centre_min+Drate_delta_one_half_min+Drate_delta_one_and_one_half_min)/5;

% Therefore the +/- errors are:

Drate_slab_avg_plus_err_2 = Drate_slab_avg_max - Drate_slab_avg_2;
Drate_slab_avg_minus_err_2 = Drate_slab_avg_2 - Drate_slab_avg_min;

tElasped_2=toc(tStart_2); % END OF TIMER FOR METHOD 2

fprintf(1,'METHOD #2 ("multi-isoplane"):\n\n')
fprintf(1,'The average dose rate to the slab of %s cells\n',Cell_name)
fprintf(1,'after uptake of the radiopharmaceutical is:\n** %d ', Drate_slab_avg_2)
fprintf(1,' +/- %d/%d Gy/h**\n\n', Drate_slab_avg_plus_err_2,Drate_slab_avg_minus_err_2)
fprintf(1,'Total time taken = %.5f sec\n\n', tElasped_2)
fprintf(1,'~~~~~~~\n\n')
fprintf(1,'______________________________________________________\n')
fprintf(1,'======================================================\n\n\n')
% ====================================================================== %




