%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%************************************************************************%
%                                                                        %
% 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 program will allow the user to input variables required
% to carry out the calculations.  These calculations will be based on the
% Vynckier-Wambersie point source dose function (kernel) integrated for the
% appropriate geometries.  In essence, we will use the dose plane kernel
% derived from VW's function.  The following articles are essential
% reference material:
%
% [1] Vynckier S and Wambersie A. Dosimetry of beta sources in radiotherapy
% I. The beta point source dose function. Phys. Med. Bio. Vol 27(11), pp.
% 1339-1347 (1982).
%
% [2] 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).
%
% [3] 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).
%
% And of course the primary article outlining the I-131 MIBG treatment
% methodology can be found in:
%
% [4] Boyd M et al. Radiation-Induced Biologic Bystander Effect Elicited In
% Vitro by Targeted Radiopharmaceuticals Labeled with alpha-, beta-, and
% Auger Electron - Emitting Radionuclides. J. Nucl. Med. Vol. 47, pp.
% 1007-1015 (2006).
%
% Lastly, reference material indicating a linear uptake of MIBG
% radiopharmaceutical across various neuroblastoma cell lines (see uptake
% section below) is:
%
% [5] Armour A, Mairs RJ, Gaze MN, and Wheldon TE. Modification of
% meta-iodobenzylguanidine uptake in neuroblastoma cells by elevated
% temperature. Br. J. Cancer. Vol. 70, pp. 445-8 (1994).
%
% The Vynckier Wambersie Point Source Dose Function (Dose Point Kernel) is
% as follows:
%
% J(x)=
% (B/(pvx)^2)*{c[1-((pvx)/c)*exp(1-(pvx)/c)]+(pvx)*exp(1-pvx)
% -(pvx)*exp(1-((pvx)/2)-f/2)}
%
% with [  ] defined as 0 for pvx >= c
% and J(x) defined as 0 for pvx >= f
%
% J(x) = is the absorbed dose at a distance x (in cm) from the point source
% in unit of Gy/(MBq*h); v = apparent absorption coefficient in cm^2/g; p =
% density of the medium in g/cm^3; c = dimensionless parameter, which gives
% the value of the first term inside the curly brackets {} at x=0 and the
% value of pvx at which this term becomes and remains zero; B =
% normalization constant evaluated by the requirement that the total energy
% absorbed by a very large sphere (infinite) must just be the energy
% emitted. Carrying out the integration of J(x) over all x yields:
%
% B = 0.046*p^2*v^3*Eavg*alpha
%
% where Eavg is the average beta energy of the radionuclide and alpha =
% [3c^2-(c^2-1)e+(3+f)exp(1-f)-4exp(1-(f/2))]^-1
%
% The parameter f/pv is the distance where the point kernel becomes zero.
% This is the supplementary term that was added the Loevinger function in
% order to make the beta point kernel equal to zero at f/pv.  The distance
% f/pv is close but always smaller than the CSDA Rmax for the radionuclide,
% since only a very small fraction of the electrons possess the maximum
% energy and additionally the probability that these electrons reach Rmax
% becomes negligible.
%
% ---------------------------------------------------------------------- %
%
% The irradiation to the donor cells used in ICCM production occurs in two
% phases  (see figure in associated report):
%
% 1) Dose Rate during Uptake Phase (Gy/hr):  Upon application of the
% radiopharmaceutical, we assume an instant, homogenous, static mixture
% into the medium.  According to Ref. (4), at the time of treatment, 1 ml
% of medium was applied followed by the radiopharmaceutical with the
% appropriate activity.  We will make the assumption that the now
% radioactive medium is a source of infinite extent but finite thickness.
% In fact, the thickness is equal to 0.04 cm.
%
% (Note: This is easily obtained by noting that flask itself is 25cm^2, and
% 1 ml of liquid covering the inside of the flask, we have the volume of
% fluid = surface area x thickness.  With 1ml of fluid equivalent to 1
% cm^3, we get the fluid thickness to be 1cm^3/25cm^2 = 0.04cm.  Hence this
% is the thickness of the source).
%
% There are a couple items to mention at this juncture
%
% 1) We can not approximate the source as infinite in along all extents.
% That is, we need to take thickness into account.  The slab can only be
% considered as infinitely thick if it's thickness is greater the range of
% the most energetic beta particle (Ref. 2).  For I-131, Emax is 0.61 MeV
% and a CSDA range of approximately 0.2 cm (Note: CSDA values for electrons
% / beta particles can be found at the National Institute of Standards and
% Technology website:
% http://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html\). So, we will
% need to take this into account in the uptake phase.
%
% 2)  We can approximate the extent (i.e. length and width) of the source
% as "infinite" provided the radius of the source is >= 0.5*Rmax where Rmax
% is the CSDA range of the most energetic beta particle (Ref. 3).  From
% above, we have that the CSDA range for I-131 is ~ 0.227cm.  Under these
% conditions, the VW plane kernel can be modelled as infinite in extent
% down to a cylindrical slab of cells of approximate surface area ~ 0.04
% cm^2.  With the confluence in Ref [4] being 65% (+/- 5%), this covers an
% area of 16.25 cm^2. Thus, rather than model each donor flask
% independently with an appropriate amount of "infinite VW pixels" of
% surface area 0.04 cm^2 (for 16.25 cm^2, this results in over 400
% "infinite VW pixels") and summing the resultant dose in each pixel for
% the total dose, we perform a simple transformation, assuming all cells
% make up a single infinite slab. This single slab will provide the same
% dose as addition of the appropriate number of "infinite VW pixels" with
% identical performance but far less complication.
%
% Now the dose rate at the surface of an infinite thick plane source is
% equal to the dose rate of a semi - infinite medium or half of the dose
% rate of an infinite medium. This is due to symmetry reasons as the dose
% rate in the inside of an infinite medium, with a homogeneously
% distributed beta source is equal to the energy dissipated per unit mass
% and time (see Ref. [1]-[3]). The dose rate is a function of the distance
% to the source, x, height of the source, h, and diameter of the source, d,
% we have the dose rate at the surface in contact with the infinite medium
% (i.e. Drate(x=0,h=inf,d=inf) is:
%
% Drate(0,inf,inf) [Gy/s] = 0.5*Am[Bq/g]*Eavg[J]                        (1)
% where Am is the activity per unit mass and Eavg is the average energy per
% beta disintegration.  In more convinent units, we can represent this as:
%
% Drate(0,inf,inf) [Gy/h] = 0.288*Am[MBq/g]*Eavg(MeV)                   (2)
%
% Now, since we will be assuming that the density of the medium is
% equivalent to that of water, this is also equivalent to:
%
% Drate(0,inf,inf) [Gy/h] = 0.288*Ac[MBq/ml]*Eavg(MeV)                  (3)
% where Ac is the activity concentration of the radiopharmaceutical
% applied.  {Remember 1 ml water = 1 cm^3 water, and density of water is 1
% g/cm^3 and accordingly 1 MBq/ml = 1MBq/cm^3 = 1MBq/g}
%
% Now, the dose rate a distance x from an infinite medium is calculated as
% follows (see Ref. [2]):
%
% Drate(x,inf,inf)=Drate(0,inf,inf)*alpha*{c^2[3-exp(1-pvx)-((pvx)/c)*(2+ln
% (c/(pvx))]+exp(1-pvx)-4*exp(1-((pvx)/c)-f/2)+(3+f-pvx)*exp(1-f)}      (4)
%
% with [ ] defined as 0 for pvx >= c and Drate(x,inf,inf) defined as 0 for
% pvx >= f
%
% So in order to take into account the does rate @ x of a source of
% infinite extent but finite thickness, we have by symmetry (Ref. [2]):
%
% Drate(x,h,inf) = Drate (x,inf,inf) - Drate(x+h,inf,inf)
%
% So, we calculate the dose rate at the onset of radiopharmaceutical
% treatment at two points
%
% a) in contact with the edge of the source and,
% b) at the opposite side of the slab a distance equal to the cell slab
% thickness away.
%
% So for a), we would get:
%
% Drate(0,0.04,inf) = Drate(0,inf,inf) - Drate(0.04,inf,inf)            (5)
%
% and for b), we would get:
%
% Drate(Cell_slab_thick,0.04,inf) = Drate(Cell_slab_thick,inf,inf) -
%                                   Drate(Cell_slab_thick+0.04,inf,inf) (6)
%
% and finally:
%
% Drate_uptake = (Drate(0,0.4,inf)+Drate(Cell_slab_thick,0.04,inf))/2   (7)
%
% Alternatively, you can take the average integral as follows:
%
% Drate_uptake = (1/Cell_slab_thick)*int(Drate(x,h,inf),dx)             (8)
%                 between 0 and Cell_slab_thick
%
% The integration is more computationally expensive with little
% gain in precision within uncertainties due to the thickness of the cells
% upon adhesion to the flask.  Additionally, because the thickness of the
% slab is small, there is effectively a less than 2% difference in dose
% rate between the top and bottom of the slab and therefore a less than 1%
% difference between the arithmetic mean and any point in the slab.  Thus,
% this program utilizes Eq. (7).
%
% ~~~~~~~~~~~~~~~~~~~
%
% 2) Dose Rate from Internal Source (Gy/hr) (Accumulation Phase):  After
% the Uptake phase, fresh medium replaces the treated medium and the
% radiopharmaceutical source is now inside the slab of cells.  To estimate
% the dose rate inside the cells, there are few different approaches. Here,
% we will break up the activity into infinite thin disks, strategically
% placed within the slab.
%
% From Ref. 2 and 3, we have for an infinite plane source:
%
% 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)}                                    (9)
%
% 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.]
%
% We will calculate the dose rate at various locations in the slab as
% through two scenarios (Note: the slab is located at x=0 to 2delta):
%
% 1)  Break the activity into 2 infinite thin plane sources of equal
% activity (0.5 Aconc) located at x=(0.5)delta and x=(1.5)delta.  Calculate
% Drate(0,0,inf), Drate(delta,o,inf), and Drate(2delta,0,inf).
%
% 2)  Break the activity into 3 infinite thin plane sources of equal
% activity ((1/3) Aconc) located at x=0,delta, and 2delta.  Calculate
% Drate ((0.5)delta,0,inf) and Drate((1.5)delta,0,inf).
%
% We can will then calculate the average dose rate of the above values in
% order to obtain the average dose rate top the slab.
%
% The another approach is (briefly):
%
% Divide 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.
% This method is computationally very expensive, requiring symbolic
% integration through either the Maple symbolic engine (best if possible)
% or Matlab's MuPAD symbolic engine.
%
% This method yields results in agreement for average dose rates to
% the slab with a decrease in computational performance and/or
% increases in complexity with no increase in precision or accuracy to the
% method employed here.
%
% ~~~~~~~~~~~~~~~~~~~
%
% 3) Total Absorbed Dose for ICCM production (Gy):
%
% Since we have an approximate linear uptake of radiopharmaceutical over
% the given uptake phase, and we collect bystander factors in the fresh
% medium for a specified period after uptake, we simply need to calculate:
%
% a) The instant dose rate upon application of the radiopharmaceutical to
% the medium (Step 1).
% b) The dose rate from the internal source alone (Step 2).
%
% From here, we can use the trapezoid rule of integration for the specified
% periods of time for uptake and bystander factor collection in order to
% determine the total dose to the donor cells providing the irradiated cell
% conditioned medium (ICCM).  Please refer to the associated report for
% complete details.
%
%   ************************ 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 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 being used for the calculation of absorbed dose from the beta emitting radiopharmaceutical ^1^3^1I meta-iodobenzylgaunidine (MIBG).  This is accomplished through modelling of the 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 that accumulate at extra-nuclear sites (like MIBG) with maximum energies ranging from 0.5 to 3.0 MeV and utilization of Boyd et al''s methodology outlined in the paper above.',
                '                                               '
                'The following sections will guide the user through data entry necessary to perform these calculations.  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

% The following is for an exit 'while' loop that will be prompted to the
% user upon completion of the calculations at the end of the program
exit_note='Yes & Continue';
while strcmp(exit_note,'Yes & Continue')

% Now, let's provide the user with the option to save the output of their
% calculations to a Microsoft Excel file.

save_text = {'Do you have Microsoft Excel 97 or greater installed?'
            '                                               '
            'If so, you can save the calculation output to an Excel file in one of the following formats:'
            '                                               '
            'a) .xls -> Excel 97-2003'
            'b) .xlsx, .xlsm, .xlsb -> Excel 2007, 2010'
            '                                               '
            'Note: If you have Excel 2003 but wish to save to .xlsx format, please ensure that you have the Excel 2007 Compatibility Pack installed on your local machine.'
            '                                               '
            'Select ''Yes'' to enter the file details or ''No'' to proceed without saving to a file (results will simply print to the command window)'
            '                                               '
            '                                               '};
save_title = 'Save to Excel?';
str1 = 'Yes';
str2 = 'No';
str3 = 'Exit';
options.Resize='on';
options.Interpreter='tex';
options.Default=str1;
save = questdlg(save_text, save_title, str1, str2, str3, options);

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

    if strcmp (save,'Yes')
        prompt_choices = {'Select an Excel file format:'};
        format_choices = {'.xls','.xlsx','.xlsm','.xlsb'};
        size = [180 60];
        excel_format = listdlg('PromptString',prompt_choices,...
                               'SelectionMode','single',...
                               'ListString',format_choices,...
                               'Name','Excel Format',...
                               'ListSize',size);
        if excel_format == 1 || excel_format == 2 || excel_format == 3 || excel_format == 4
            filename_text = {'Please enter the name of the file.  Do not enter the extension (e.g only ''Results'' NOT ''Results.xls'').  This file will be saved to your current active directory (i.e. the directory where you are currently running the .exe or .m file).'};
            title_filename = 'Filename for Calculation Results';
            num_lines_filename = 1;
            options.Resize='on';
            options.Interpreter='tex';

            % Define the standard parameters
            def_filename = {'results'};
            answer_filename = inputdlg(filename_text,title_filename,num_lines_filename,def_filename,options);

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

                % We are going to make sure the user does not enter an
                % invalid character for Windows file systems
                check=0;
                wrong_char1 = strfind(answer_filename{1},'.');
                    if isempty(wrong_char1)==0
                        check=1;
                    end
                wrong_char2 = strfind(answer_filename{1},''\'');
                    if isempty(wrong_char2)==0
                        check=1;
                    end
                wrong_char3 = strfind(answer_filename{1},'/');
                    if isempty(wrong_char3)==0
                        check=1;
                    end
                wrong_char4 = strfind(answer_filename{1},':');
                    if isempty(wrong_char4)==0
                        check=1;
                    end
                wrong_char5 = strfind(answer_filename{1},'*');
                    if isempty(wrong_char5)==0
                        check=1;
                    end
                wrong_char6 = strfind(answer_filename{1},'?');
                    if isempty(wrong_char6)==0
                        check=1;
                    end
                wrong_char7 = strfind(answer_filename{1},'"');
                    if isempty(wrong_char7)==0
                        check=1;
                    end
                wrong_char8 = strfind(answer_filename{1},'<');
                    if isempty(wrong_char8)==0
                        check=1;
                    end
                wrong_char9 = strfind(answer_filename{1},'>');
                    if isempty(wrong_char9)==0
                        check=1;
                    end
                wrong_char10 = strfind(answer_filename{1},'|');
                    if isempty(wrong_char10)==0
                        check=1;
                    end
                wrong_char11 = strfind(answer_filename{1},'%');
                    if isempty(wrong_char11)==0
                        check=1;
                    end

                while isempty(answer_filename{1}) || check == 1
                    err_filename_dlg = errordlg({'You need to enter a valid file name!'
                                                 'It can not be empty or contain a one of the following characters:'
                                                 '                                     '
                                                 '\  /  :  *  ?  "  <  >  |  .  %'},'Error');
                    uiwait(err_filename_dlg);
                    check = 0;
                    answer_filename = inputdlg(filename_text,title_filename,num_lines_filename,def_filename,options);
                        if isempty(answer_filename) == 1
                            disp('Bye Bye')
                            return
                        end
                end
        end
    end

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

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

% Now let's have the user start entering values but first, a message to the
% user saying what we are doing...

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].  Provided the treatment flask is largely confluent, then the majority of cells are in close contact and such an approximation holds true.'
                 '                                                       '
                 '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: If entering a fractional number (i.e. less than 1), make sure to enter a leading '' 0. '' where applicable.'
                 '                                                       '
                 '                                                       '
                 '                                                       '};
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 variable to calculate the
% cell thickness when adhered to the bottom of the treatment flask

prompt_slab = {'1) s_f_l_a_s_k = Size (surface 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 (surface 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 100]; % 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 enter a leading '' 0. '' for fractional numbers 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;
% *******

% Now, lets print the result!
fprintf(1,'======================================================\n')
fprintf(1,'The thickness of the slab of %s cells is:\n--> %d', Cell_name,Cell_slab_thick)
fprintf(1,' +/- %.0d cm\n', Err_Cell_slab_thick)
fprintf(1,'\nTherefore, the half-thickness of the slab of %s\n',Cell_name)
fprintf(1,'cells is:\n--> %d',Cell_slab_half_thick)
fprintf(1,' +/- %.0d cm\n======================================================\n\n', Err_Cell_slab_half_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*% Uptake from exp.
%
% 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, homogeneous
% 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 = {'In order to perform the necessary calculations, we need to account for the following additional experimental parameters regarding the radiopharmaceutical and medium utilized:'
                   '                                                       '
                   '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 concentrations applied'
                   '5) Percentage uptake of the radiopharmaceutical into the cell in vitro'
                   '6) Time for uptake phase'
                   '7) Time for bystander factor accumulation'
                   '                                                       '
                   'NOTE: If entering a fractional number (i.e. less than 1), make sure to enter a leading '' 0. '' where applicable.'
                   '                                                       '
                   '                                                       '};

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):'
                    'Smallest A_c_o_n_c = Smallest activity concentration of radiopharmaceutical applied (in MBq/ml; With medium of unit density, then 1 MBq/ml = 1 MBq/cm^3 = 1 MBq/g). NOTE: If the experiment involved only one A_c_o_n_c, enter it here:'
                    'Largest A_c_o_n_c = Largest activity concentration of radiopharmaceutical applied (in MBq/ml; With medium of unit density, then 1 MBq/ml = 1 MBq/cm^3 = 1 MBq/g). NOTE: If the experiment involved only one A_c_o_n_c, enter 0 here:'
                    'A_c_o_n_c Increments = Enter the incremental steps utilized for applying the range of A_c_o_n_c''s across the experiment (ex. Experiments with 1,2, & 3 Mq/ml A_c_o_n_c''s have A_c_o_n_c Increments = 1).  If your increments were not equal, then each A_c_o_n_c value should be entered one at a time into "Smallest A_c_o_n_c" field and "A_c_o_n_c Increments" should remain 0.  Lastly, if you enter a value for Largest A_c_o_n_c, then the appropriate step value must be entered:'
                    '%_u_p_t_a_k_e = Percentage uptake of radiopharmaceutical (between 0-1):'
                    'Err_%_u_p_t_a_k_e = Error / Uncertainty in uptake of radiopharmaceutical (between 0 and %_u_p_t_a_k_e):'
                    'T_u_p = Time allotted for radiopharmaceutical uptake (in hours):'
                    'T_a_c_c_u_m = Time allotted for bystander factors to accumulate after application of fresh medium (in hours):'};

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

% Define the standard parameters
def_addn_param = {'0.19','0.61','1.00','1.00','0.00','0','0.314','0.037','2','1'};

% "answer" is the 10 element column vector that contains the user input
% results.
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}) || isempty(addn_param_answer{7}) || isempty(addn_param_answer{8}) || isempty(addn_param_answer{9}) || isempty(addn_param_answer{10})
        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});
% Smallest A_conc or only A_conc applied (in MBq/ml)
Small_act_conc = str2double(addn_param_answer{4});
% Largest A_conc(in MBq/ml)
Large_act_conc = str2double(addn_param_answer{5});
% A_conc increments
Increm_act_conc = str2double(addn_param_answer{6});
% % of radiopharmaceutical uptake
Percent_up = str2double(addn_param_answer{7});
% Error / Uncertainty in % of radiopharmaceutical uptake
Err_percent_up = str2double(addn_param_answer{8});
% Time allocated for radiopharmaceutical uptake
Time_up = str2double(addn_param_answer{9});
% Time allocated for bystander factor accumulation
Time_accum = str2double(addn_param_answer{10});

    % 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) || (Small_act_conc < 0) || (isnan(Small_act_conc) == 1)|| (Large_act_conc < 0) || (isnan(Large_act_conc) == 1) || (Small_act_conc > Large_act_conc && Large_act_conc ~= 0) || (Increm_act_conc < 0) || (isnan(Increm_act_conc) ==1) || (Increm_act_conc > Large_act_conc) || (Increm_act_conc == 0 && Large_act_conc ~= 0) || (Percent_up <= 0) || (isnan(Percent_up) == 1) || (Percent_up > 1) || (Err_percent_up <= 0) || (isnan(Err_percent_up) == 1) || (Err_percent_up > Percent_up) || (Time_up <= 0) || (isnan(Time_up) == 1) || (Time_accum <= 0) || (isnan(Time_accum) == 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 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).'
                                 '                                                  '
                                 'c) The "Smallest" activity concentration entered is bigger than the "Largest" activity concentration entered.'
                                 '                                                  '
                                 'd) The incremental steps entered for application of the activity concentrations is bigger than the "Largest" activity concentration of a "Large" activity concentration was entered without an "Increment" value.'
                                 '                                                  '
                                 '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 enter a leading '' 0. '' for fractional numbers 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});
        Small_act_conc = str2double(addn_param_answer{4});
        Large_act_conc = str2double(addn_param_answer{5});
        Increm_act_conc = str2double(addn_param_answer{6});
        Percent_up = str2double(addn_param_answer{7});
        Err_percent_up = str2double(addn_param_answer{8});
        Time_up = str2double(addn_param_answer{9});
        Time_accum = str2double(addn_param_answer{10});
    end

% Now, we have Conf_adjust_uptake as follows:

if (Conf_uptake-Err_conf_uptake) >= S_well
    Conf_uptake = 1;
    Err_conf_uptake = 0;
    %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:
end

Conf_adjust_uptake = Per_conf/Conf_uptake;
Err_conf_adjust_uptake = (Per_conf/Conf_uptake)*((Err_per_conf/Per_conf)^2 + (Err_conf_uptake/Conf_uptake)^2)^0.5;


% *******
% 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;
% *******

% 5) Thickness of radiopharmaceutical medium, h (in cm).  Based the
% protocol utilized by Boyd et al, the radiopharmaceutical was added to 1
% ml of medium during the uptake phase.  Based on the size of the flask
% used, S_flask, we have A*t = Vol (of 1 ml) = 1 cm^3.  Therefore, we have
% 1 cm^3/25cm^2 =:

h = 1/S_flask;
% ====================================================================== %

% Let's start the calculations.  If the user enters a range of activity
% concentrations, we have the following:

if Large_act_conc ~= 0

% Start of 'for' loop to cycle through the activity concentrations entered
% by the user and perform the calculations...

i = 0; % This a counter for us to use in the 'for' loop as a matrix index

% For performance, we pre-allocate memory space for vectors we will be
% creating throughout the calculations:
Activity = zeros(1,(Large_act_conc/Increm_act_conc));
Drate_uptake_avg = zeros(1,(Large_act_conc/Increm_act_conc));
Drate_uptake_x_2delta_plus_err = zeros(1,(Large_act_conc/Increm_act_conc));
Drate_uptake_x_2delta_minus_err = zeros(1,(Large_act_conc/Increm_act_conc));
Drate_slab_avg = zeros(1,(Large_act_conc/Increm_act_conc));
Drate_slab_avg_max = zeros(1,(Large_act_conc/Increm_act_conc));
Drate_slab_avg_min = zeros(1,(Large_act_conc/Increm_act_conc));
Drate_slab_avg_plus_err = zeros(1,(Large_act_conc/Increm_act_conc));
Drate_slab_avg_minus_err = zeros(1,(Large_act_conc/Increm_act_conc));
Dose = zeros(1,(Large_act_conc/Increm_act_conc));
Dose_err_plus = zeros(1,(Large_act_conc/Increm_act_conc));
Dose_err_minus = zeros(1,(Large_act_conc/Increm_act_conc));

for A = Small_act_conc:Increm_act_conc:Large_act_conc
i = i + 1; % start matrix index counter and add 1 for each run through

fprintf(1,'For an activity concentration of: %4.2f MBq/ml\n',A)
fprintf(1,'~~~~~~~\n\n')

Activity(i) = A;

% ====================================================================== %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~Calculate Dose Rate for Uptake~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% We will now calculate the dose rate during the uptake phase as noted in
% the section the "Brief Introduction".  The 'if' statements seen below
% take into consideration the appropriate conditions from Eq. (4):

% First, if the f <= p*v*x ...
if f <= (Dens_med*v*Cell_slab_thick_min)
% Note we have used the minimum cell thickness calculated to be the
% most stringent

    fprintf(1,'The dose rate is 0 Gy/h. No need to go any further...\n')

else % Second, condition on c to alter value of the square brackets

    % First, as stated in the "Brief Introduction" above, we need to
    % calculate this in two parts, a and b (see Eqs. 4,5,and 6).  Step one
    % is calculate Drate(0,0.04,inf) = Drate(0,inf,inf) -
    % Drate(0.04,inf,inf) where Drate(0,inf,inf) [Gy/h] =
    % 0.288*Ac[MBq/ml]*Eavg(MeV):

    Drate_zero_inf_inf = 0.288*A*Per_conf*Mean_beta_eng;
    % Remember, we take into account the confluence adjustment as activity
    % not covering cells does not contribute to the absorbed dose.  This is
    % simply the percentage confluence during the Uptake phase.  During the
    % Accumulation phase, this is the Conf_adjust_uptake factor calculated
    % previously (which is simply Per_conf if Uptake experiments were
    % confluent.

    Drate_h_inf_inf = (Drate_zero_inf_inf)*alpha*(c^2*(3-exp(1-(Dens_med*v*h))-((Dens_med*v*h)/c)*(2+log(c/(Dens_med*v*h))))+exp(1-(Dens_med*v*h))-4*exp(1-((Dens_med*v*h)/2)-(f/2))+(3+f-(Dens_med*v*h))*exp(1-f));

        if c <= (Dens_med*v*h)
            Drate_zero_inf_inf = 0.288*A*Per_conf*Mean_beta_eng;
            Drate_h_inf_inf = (Drate_zero_inf_inf)*alpha*(exp(1-(Dens_med*v*h))-4*exp(1-((Dens_med*v*h)/2)-(f/2))+(3+f-(Dens_med*v*h))*exp(1-f));
        end

    Drate_zero_h_inf = Drate_zero_inf_inf - Drate_h_inf_inf;

    % Step 2 is to calculate
    % Drate(Cell_slab_thick,0.04,inf) = Drate(Cell_slab_thick,inf,inf) -
    %                                   Drate(Cell_slab_thick+0.04,inf,inf)

    Drate_2delta_inf_inf = (Drate_zero_inf_inf)*alpha*(c^2*(3-exp(1-(Dens_med*v*Cell_slab_thick))-((Dens_med*v*Cell_slab_thick)/c)*(2+log(c/(Dens_med*v*Cell_slab_thick))))+exp(1-(Dens_med*v*Cell_slab_thick))-4*exp(1-((Dens_med*v*Cell_slab_thick)/2)-(f/2))+(3+f-(Dens_med*v*Cell_slab_thick))*exp(1-f));
    Drate_2delta_h_inf_inf = (Drate_zero_inf_inf)*alpha*(c^2*(3-exp(1-(Dens_med*v*(Cell_slab_thick+h)))-((Dens_med*v*(Cell_slab_thick+h))/c)*(2+log(c/(Dens_med*v*(Cell_slab_thick+h)))))+exp(1-(Dens_med*v*(Cell_slab_thick+h)))-4*exp(1-((Dens_med*v*(Cell_slab_thick+h))/2)-(f/2))+(3+f-(Dens_med*v*(Cell_slab_thick+h)))*exp(1-f));

        if c <= (Dens_med*v*Cell_slab_thick)
            Drate_2delta_inf_inf = (Drate_zero_inf_inf)*alpha*(exp(1-(Dens_med*v*Cell_slab_thick))-4*exp(1-((Dens_med*v*Cell_slab_thick)/2)-(f/2))+(3+f-(Dens_med*v*Cell_slab_thick))*exp(1-f));
            Drate_2delta_h_inf_inf = (Drate_zero_inf_inf)*alpha*(exp(1-(Dens_med*v*(Cell_slab_thick+h)))-4*exp(1-((Dens_med*v*(Cell_slab_thick+h))/2)-(f/2))+(3+f-(Dens_med*v*(Cell_slab_thick+h)))*exp(1-f));
        end

    Drate_2delta_h_inf = Drate_2delta_inf_inf - Drate_2delta_h_inf_inf;

    % Since the thickness of the cells is so small (typically ~2%), we
    % expect that the difference between of dose rate to the cells x = o
    % and x = 2delta to be nearly equal.  As such, for performance, rather
    % than integrate as per Eq. (8) above, we will simply calculate the
    % average Drate by adding Drate_zero_h_inf and Drate_2delta_h_inf and
    % dividing by 2 [Eq. (7)].

    Drate_uptake_avg(i) = (Drate_zero_h_inf + Drate_2delta_h_inf)/2;

    % Now we need to take into account the variability / uncertainty in the
    % slab thickness.  This obviously, only effects the value calculated at
    % 2delta.  So we will set up a couple maximum and minimum variables.
    % This will enable to have the full breadth of possible dose rate
    % calculations by appropriately optimizing the values from the
    % experimental parameters for each scenario.

    % Maximum:

    % Maximize activity for uptake
    Drate_zero_inf_inf_max = 0.288*A*(Per_conf+Err_per_conf)*Mean_beta_eng;
    % Maximize dose rate by slab being of minimal thickness
    Drate_uptake_x_2delta_max_1=(Drate_zero_inf_inf_max)*alpha*(c^2*(3-exp(1-(Dens_med*v*Cell_slab_thick_min))-((Dens_med*v*Cell_slab_thick_min)/c)*(2+log(c/(Dens_med*v*Cell_slab_thick_min))))+exp(1-(Dens_med*v*Cell_slab_thick_min))-4*exp(1-((Dens_med*v*Cell_slab_thick_min)/2)-(f/2))+(3+f-(Dens_med*v*Cell_slab_thick_min))*exp(1-f));
    Drate_uptake_x_2delta_max_2=(Drate_zero_inf_inf_max)*alpha*(c^2*(3-exp(1-(Dens_med*v*(Cell_slab_thick_min+h)))-((Dens_med*v*(Cell_slab_thick_min+h))/c)*(2+log(c/(Dens_med*v*(Cell_slab_thick_min+h)))))+exp(1-(Dens_med*v*(Cell_slab_thick_min+h)))-4*exp(1-((Dens_med*v*(Cell_slab_thick_min+h))/2)-(f/2))+(3+f-(Dens_med*v*(Cell_slab_thick_min+h)))*exp(1-f));
    % Drate_uptake_x_2delta_max_2 is the correction for thickness, h

        if c <= (Dens_med*v*Cell_slab_thick_min)
            Drate_uptake_x_2delta_max_1=(Drate_zero_inf_inf_max)*alpha*(exp(1-(Dens_med*v*Cell_slab_thick_min))-4*exp(1-((Dens_med*v*Cell_slab_thick_min)/2)-(f/2))+(3+f-(Dens_med*v*Cell_slab_thick_min))*exp(1-f));
            Drate_uptake_x_2delta_max_2=(Drate_zero_inf_inf_max)*alpha*(exp(1-(Dens_med*v*(Cell_slab_thick_min+h)))-4*exp(1-((Dens_med*v*(Cell_slab_thick_min+h))/2)-(f/2))+(3+f-(Dens_med*v*(Cell_slab_thick_min+h)))*exp(1-f));
        end


    Drate_uptake_x_2delta_max = Drate_uptake_x_2delta_max_1 - Drate_uptake_x_2delta_max_2;

    % Minimum:

    % Minimize activity for uptake
    Drate_zero_inf_inf_min = 0.288*A*(Per_conf-Err_per_conf)*Mean_beta_eng;
    % Minimize dose rate by slab being of maximal thickness
    Drate_uptake_x_2delta_min_1=(Drate_zero_inf_inf_min)*alpha*(c^2*(3-exp(1-(Dens_med*v*Cell_slab_thick_max))-((Dens_med*v*Cell_slab_thick_max)/c)*(2+log(c/(Dens_med*v*Cell_slab_thick_max))))+exp(1-(Dens_med*v*Cell_slab_thick_max))-4*exp(1-((Dens_med*v*Cell_slab_thick_max)/2)-(f/2))+(3+f-(Dens_med*v*Cell_slab_thick_max))*exp(1-f));
    Drate_uptake_x_2delta_min_2=(Drate_zero_inf_inf_min)*alpha*(c^2*(3-exp(1-(Dens_med*v*(Cell_slab_thick_max+h)))-((Dens_med*v*(Cell_slab_thick_max+h))/c)*(2+log(c/(Dens_med*v*(Cell_slab_thick_max+h)))))+exp(1-(Dens_med*v*(Cell_slab_thick_max+h)))-4*exp(1-((Dens_med*v*(Cell_slab_thick_max+h))/2)-(f/2))+(3+f-(Dens_med*v*(Cell_slab_thick_max+h)))*exp(1-f));
    % Drate_uptake_x_2delta_min_2 is the correction for thickness, h

        if c <= (Dens_med*v*Cell_slab_thick_max)
            Drate_uptake_x_2delta_min_1=(Drate_zero_inf_inf_min)*alpha*(exp(1-(Dens_med*v*Cell_slab_thick_max))-4*exp(1-((Dens_med*v*Cell_slab_thick_max)/2)-(f/2))+(3+f-(Dens_med*v*Cell_slab_thick_max))*exp(1-f));
            Drate_uptake_x_2delta_min_2=(Drate_zero_inf_inf_min)*alpha*(exp(1-(Dens_med*v*(Cell_slab_thick_max+h)))-4*exp(1-((Dens_med*v*(Cell_slab_thick_max+h))/2)-(f/2))+(3+f-(Dens_med*v*(Cell_slab_thick_max+h)))*exp(1-f));
        end


    Drate_uptake_x_2delta_min = Drate_uptake_x_2delta_min_1 - Drate_uptake_x_2delta_min_2;

    % Therefore, 2Delta Plus error
    Drate_uptake_x_2delta_plus_err(i) = Drate_uptake_x_2delta_max - Drate_2delta_h_inf;

    % and 2Delta Minus error
    Drate_uptake_x_2delta_minus_err(i) = Drate_2delta_h_inf - Drate_uptake_x_2delta_min;

    % Now, lets see if the plus and minus error are in fact equal to one
    % significant digit.  MATLAB does not have a trivial way to accomplish
    % this, so the following is the work around:

    % First, define a string for the plus error output:
    str_plus=sprintf('%.0d',Drate_uptake_x_2delta_plus_err(i));
    % Second, define a string for the plus error output:
    str_minus=sprintf('%.0d',Drate_uptake_x_2delta_minus_err(i));
    % The '.0d' rounds the values to the nearest whole integer in e -
    % notation. For example, 2.13e-005 becomes 2e-005.

    % Convert strings to double precision numbers:
    str_plus_num = str2double(str_plus);
    str_minus_num = str2double(str_minus);

    % Now, let's create a variable that will compare the first value in
    % str_plus and str_minus and see if they are the same:

    str_error_test=strncmp(str_plus,str_minus,1);
    % If this test is true, we only need to print one of the uncertainty
    % values (either the plus or the minus) as they are equal to one
    % significant digit.  If this is false, we will print out both values.

    % Now, lets print the result!

    fprintf(1,'The average dose rate to the %s cells at \n',Cell_name)
    fprintf(1,'the initiation of the uptake phase is:\n** %d', Drate_uptake_avg(i))

        if (str_error_test == 1)
            fprintf(1,' +/- %.0d', str_plus_num)
            fprintf(1,' Gy/h ** [1]\n\n')
        elseif (str_error_test == 0)
            fprintf(1,' +/- (%.0d/%.0d)',str_plus_num,str_minus_num)
            fprintf(1,' Gy/h ** [1]\n\n')
        end

    fprintf(1,'~~~~~~~\n\n')
end

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

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~Calculate Dose Rate in Slab~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% For our methodology employed here, see "Brief Introduction" above.
% 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*Asurf 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).

% If all activity were isolated to a single plane, the surface activity of
% the plane source would be:
%
% Asurf = (A*%Uptake*Conf_adjust)/(Size of flask*%Confluence)
%
% Therefore, for paradigm 1, Atot = 0.5Asurf + 0.5Asurf
%        and for paradigm 2, Atot = (1/3)Asurf + (1/3)Asurf + (1/3)Asurf

% The expression for dose rate from an infinite disk at distance x from
% source is given in Eq. (9) and provided in Ref. [2] and [3].  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)((A.*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));

% Again, 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)(A.*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(i) = (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 % Uptake is highest, slab thickness is smallest, and
    % confluence is smallest and vice verse for the minimum dose rate.

    % Maximum Drate:
    Drate_x_max=@(x)((A.*(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)((A.*(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(i) = (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)((A.*(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)(A.*(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=(1.5)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(i) = (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(i) = Drate_slab_avg_max(i) - Drate_slab_avg(i);
Drate_slab_avg_minus_err(i) = Drate_slab_avg(i) - Drate_slab_avg_min(i);

    % Now, lets see if the plus and minus error are in fact equal to one
    % significant digit.  MATLAB does not have a trivial way to accomplish
    % this, so the following is the work around:

    % First, define a string for the plus error output:
    str_plus=sprintf('%.0d',Drate_slab_avg_plus_err(i));
    % Second, define a string for the plus error output:
    str_minus=sprintf('%.0d',Drate_slab_avg_minus_err(i));
    % The '.0d' rounds the values to the nearest whole integer in e -
    % notation. For example, 2.13e-005 becomes 2e-005.

    % Convert strings to double precision numbers:
    str_plus_num = str2double(str_plus);
    str_minus_num = str2double(str_minus);

    % Now, let's create a variable that will compare the first value in
    % str_plus and str_minus and see if they are the same:

    str_error_test=strncmp(str_plus,str_minus,1);
    % If this test is true, we only need to print on of the uncertainty
    % values (either the plus or the minus) as they equal to the number of
    % significant digits allotted for uncertainty values.  If this is
    % false, we will print out both values.

% Now, lets print the result!

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(i))

    if (str_error_test == 1)
        fprintf(1,' +/- %.0d', str_plus_num)
        fprintf(1,' Gy/h ** [2]\n\n')
    elseif (str_error_test == 0)
        fprintf(1,' +/- (%.0d/%.0d)',str_plus_num,str_minus_num)
        fprintf(1,' Gy/h ** [2]\n\n')
    end

fprintf(1,'~~~~~~~\n\n')

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

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~Calculate the Total Dose to the Cells~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% Because there is believed to be a linear uptake of the pharmaceutical
% over the 2 hour time frame coupled with negligible egress of the activity
% following uptake (<1%), we can use a simple Trapezoid integration rule of
% the dose rate over the 3 hour time window to determine the total dose to
% the cells.  Essentially, we are finding the "area under the curve" of a
% Dose rate vs. time plot to get the total dose.

Dose(i) = ((Drate_uptake_avg(i)+Drate_slab_avg(i))/2)*(2)+(Drate_slab_avg(i))*(1);

    % Now calculate the uncertainties.

    Dose_err_plus(i) = (sqrt(((Drate_uptake_x_2delta_plus_err(i))^2+(Drate_slab_avg_plus_err(i))^2))/2)*(2)+Drate_slab_avg_plus_err(i)*(1);
    Dose_err_minus(i) = (sqrt(((Drate_uptake_x_2delta_minus_err(i))^2+(Drate_slab_avg_minus_err(i))^2))/2)*(2)+Drate_slab_avg_minus_err(i)*(1);

    % Now, lets see if the plus and minus error are in fact equal to one
    % significant digit.  MATLAB does not have a trivial way to accomplish
    % this, so the following is the work around:

    % First, define a string for the plus error output:
    str_plus=sprintf('%.0d',Dose_err_plus(i));
    % Second, define a string for the plus error output:
    str_minus=sprintf('%.0d',Dose_err_minus(i));
    % See MATLAB Help 'sprintf' for function parameter explanation. The
    % '.0d' rounds the values to the nearest whole integer in e - notation.
    % For example, 2.13e-005 becomes 2e-005.

    % Convert strings to double precision numbers:
    str_plus_num = str2double(str_plus);
    str_minus_num = str2double(str_minus);

    % Now, let's create a variable that will compare the first value in
    % str_plus and str_minus and see if they are the same:

    str_error_test=strncmp(str_plus,str_minus,1);
    % If this test is true, we only need to print on of the uncertainty
    % values (either the plus or the minus) as they equal to the number of
    % significant digits allotted for uncertainty values.  If this is
    % false, we will print out both values.

% Now, lets print the result!

fprintf(1,'Therefore, the dose to the %s cells after\n',Cell_name)
fprintf(1,'treatment with %.1d MBq/ml of radiopharmaceutical with\n',A)
fprintf(1,'%.0d hour(s)for uptake and %.0d hour(s) for accumulation is:\n\n',Time_up,Time_accum)

    if (str_error_test == 1)
        fprintf(1,'***********************************\n')
        fprintf(1,'*** %d',Dose(i))
        fprintf(1,' +/- %.0d', str_plus_num)
        fprintf(1,' Gy ***\n')
        fprintf(1,'***********************************\n')
    elseif (str_error_test == 0)
        fprintf(1,'********************************************\n')
        fprintf(1,'*** %d',Dose(i))
        fprintf(1,' +/- (%.0d/%.0d)',str_plus_num,str_minus_num)
        fprintf(1,' Gy ***\n')
        fprintf(1,'********************************************\n')
    end

fprintf(1,'______________________________________________________\n')
fprintf(1,'======================================================\n')
end

else if Large_act_conc == 0 % The user only enters a single value
i=1;
Activity(i) = Small_act_conc;
A = Activity(i);
fprintf(1,'For an activity concentration of: %4.2f MBq/ml\n',A)
fprintf(1,'~~~~~~~\n\n')

% ====================================================================== %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~Calculate Dose Rate for Uptake~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% First, if the f <= p*v*x ...
if f <= (Dens_med*v*Cell_slab_thick_min)
    fprintf(1,'The dose rate is 0 Gy/h. No need to go any further...\n')
else % Second, condition on c to alter value of the square brackets

    Drate_zero_inf_inf = 0.288*A*Per_conf*Mean_beta_eng;
    % Remember, we take into account the confluence adjustment as activity
    % not covering cells does not contribute to the absorbed dose.  This is
    % simply the percentage confluence during the Uptake phase.  During the
    % Accumulation phase, this is the Conf_adjust_uptake factor calculated
    % previously (which is simply Per_conf if Uptake experiments were
    % confluent.

    Drate_h_inf_inf = (Drate_zero_inf_inf)*alpha*(c^2*(3-exp(1-(Dens_med*v*h))-((Dens_med*v*h)/c)*(2+log(c/(Dens_med*v*h))))+exp(1-(Dens_med*v*h))-4*exp(1-((Dens_med*v*h)/2)-(f/2))+(3+f-(Dens_med*v*h))*exp(1-f));

        if c <= (Dens_med*v*h)
            Drate_zero_inf_inf = 0.288*A*Per_conf*Mean_beta_eng;
            Drate_h_inf_inf = (Drate_zero_inf_inf)*alpha*(exp(1-(Dens_med*v*h))-4*exp(1-((Dens_med*v*h)/2)-(f/2))+(3+f-(Dens_med*v*h))*exp(1-f));
        end

    Drate_zero_h_inf = Drate_zero_inf_inf - Drate_h_inf_inf;

    % Step 2 is to calculate
    % Drate(Cell_slab_thick,0.04,inf) = Drate(Cell_slab_thick,inf,inf) -
    %                                   Drate(Cell_slab_thick+0.04,inf,inf)

    Drate_2delta_inf_inf = (Drate_zero_inf_inf)*alpha*(c^2*(3-exp(1-(Dens_med*v*Cell_slab_thick))-((Dens_med*v*Cell_slab_thick)/c)*(2+log(c/(Dens_med*v*Cell_slab_thick))))+exp(1-(Dens_med*v*Cell_slab_thick))-4*exp(1-((Dens_med*v*Cell_slab_thick)/2)-(f/2))+(3+f-(Dens_med*v*Cell_slab_thick))*exp(1-f));
    Drate_2delta_h_inf_inf = (Drate_zero_inf_inf)*alpha*(c^2*(3-exp(1-(Dens_med*v*(Cell_slab_thick+h)))-((Dens_med*v*(Cell_slab_thick+h))/c)*(2+log(c/(Dens_med*v*(Cell_slab_thick+h)))))+exp(1-(Dens_med*v*(Cell_slab_thick+h)))-4*exp(1-((Dens_med*v*(Cell_slab_thick+h))/2)-(f/2))+(3+f-(Dens_med*v*(Cell_slab_thick+h)))*exp(1-f));

        if c <= (Dens_med*v*Cell_slab_thick)
            Drate_2delta_inf_inf = (Drate_zero_inf_inf)*alpha*(exp(1-(Dens_med*v*Cell_slab_thick))-4*exp(1-((Dens_med*v*Cell_slab_thick)/2)-(f/2))+(3+f-(Dens_med*v*Cell_slab_thick))*exp(1-f));
            Drate_2delta_h_inf_inf = (Drate_zero_inf_inf)*alpha*(exp(1-(Dens_med*v*(Cell_slab_thick+h)))-4*exp(1-((Dens_med*v*(Cell_slab_thick+h))/2)-(f/2))+(3+f-(Dens_med*v*(Cell_slab_thick+h)))*exp(1-f));
        end

    Drate_2delta_h_inf = Drate_2delta_inf_inf - Drate_2delta_h_inf_inf;

    Drate_uptake_avg(i) = (Drate_zero_h_inf + Drate_2delta_h_inf)/2;

    % Now we need to take into account the variability / uncertainty in the
    % slab thickness.  This obviously, only effects the value calculated at
    % 2delta.  So we will set up a couple maximum and minimum variables.
    % This will enable to have the full breadth of possible dose rate
    % calculations by appropriately optimizing the values from the
    % experimental parameters for each scenario.

    % Maximize activity for uptake
    Drate_zero_inf_inf_max = 0.288*A*(Per_conf+Err_per_conf)*Mean_beta_eng;

    % Maximize dose rate by slab being of minimal thickness
    Drate_uptake_x_2delta_max_1=(Drate_zero_inf_inf_max)*alpha*(c^2*(3-exp(1-(Dens_med*v*Cell_slab_thick_min))-((Dens_med*v*Cell_slab_thick_min)/c)*(2+log(c/(Dens_med*v*Cell_slab_thick_min))))+exp(1-(Dens_med*v*Cell_slab_thick_min))-4*exp(1-((Dens_med*v*Cell_slab_thick_min)/2)-(f/2))+(3+f-(Dens_med*v*Cell_slab_thick_min))*exp(1-f));
    Drate_uptake_x_2delta_max_2=(Drate_zero_inf_inf_max)*alpha*(c^2*(3-exp(1-(Dens_med*v*(Cell_slab_thick_min+h)))-((Dens_med*v*(Cell_slab_thick_min+h))/c)*(2+log(c/(Dens_med*v*(Cell_slab_thick_min+h)))))+exp(1-(Dens_med*v*(Cell_slab_thick_min+h)))-4*exp(1-((Dens_med*v*(Cell_slab_thick_min+h))/2)-(f/2))+(3+f-(Dens_med*v*(Cell_slab_thick_min+h)))*exp(1-f));
    % Drate_uptake_x_2delta_max_2 is the correction for thickness, h

        if c <= (Dens_med*v*Cell_slab_thick_min)
            Drate_uptake_x_2delta_max_1=(Drate_zero_inf_inf_max)*alpha*(exp(1-(Dens_med*v*Cell_slab_thick_min))-4*exp(1-((Dens_med*v*Cell_slab_thick_min)/2)-(f/2))+(3+f-(Dens_med*v*Cell_slab_thick_min))*exp(1-f));
            Drate_uptake_x_2delta_max_2=(Drate_zero_inf_inf_max)*alpha*(exp(1-(Dens_med*v*(Cell_slab_thick_min+h)))-4*exp(1-((Dens_med*v*(Cell_slab_thick_min+h))/2)-(f/2))+(3+f-(Dens_med*v*(Cell_slab_thick_min+h)))*exp(1-f));
        end


    Drate_uptake_x_2delta_max = Drate_uptake_x_2delta_max_1 - Drate_uptake_x_2delta_max_2;

    % Minimize activity for uptake
    Drate_zero_inf_inf_min = 0.288*A*(Per_conf-Err_per_conf)*Mean_beta_eng;

    % Minimize dose rate by slab being of maximal thickness
    Drate_uptake_x_2delta_min_1=(Drate_zero_inf_inf_min)*alpha*(c^2*(3-exp(1-(Dens_med*v*Cell_slab_thick_max))-((Dens_med*v*Cell_slab_thick_max)/c)*(2+log(c/(Dens_med*v*Cell_slab_thick_max))))+exp(1-(Dens_med*v*Cell_slab_thick_max))-4*exp(1-((Dens_med*v*Cell_slab_thick_max)/2)-(f/2))+(3+f-(Dens_med*v*Cell_slab_thick_max))*exp(1-f));
    Drate_uptake_x_2delta_min_2=(Drate_zero_inf_inf_min)*alpha*(c^2*(3-exp(1-(Dens_med*v*(Cell_slab_thick_max+h)))-((Dens_med*v*(Cell_slab_thick_max+h))/c)*(2+log(c/(Dens_med*v*(Cell_slab_thick_max+h)))))+exp(1-(Dens_med*v*(Cell_slab_thick_max+h)))-4*exp(1-((Dens_med*v*(Cell_slab_thick_max+h))/2)-(f/2))+(3+f-(Dens_med*v*(Cell_slab_thick_max+h)))*exp(1-f));
    % Drate_uptake_x_2delta_min_2 is the correction for thickness, h

        if c <= (Dens_med*v*Cell_slab_thick_max)
            Drate_uptake_x_2delta_min_1=(Drate_zero_inf_inf_min)*alpha*(exp(1-(Dens_med*v*Cell_slab_thick_max))-4*exp(1-((Dens_med*v*Cell_slab_thick_max)/2)-(f/2))+(3+f-(Dens_med*v*Cell_slab_thick_max))*exp(1-f));
            Drate_uptake_x_2delta_min_2=(Drate_zero_inf_inf_min)*alpha*(exp(1-(Dens_med*v*(Cell_slab_thick_max+h)))-4*exp(1-((Dens_med*v*(Cell_slab_thick_max+h))/2)-(f/2))+(3+f-(Dens_med*v*(Cell_slab_thick_max+h)))*exp(1-f));
        end


    Drate_uptake_x_2delta_min = Drate_uptake_x_2delta_min_1 - Drate_uptake_x_2delta_min_2;

    % Therefore, 2Delta Plus error
    Drate_uptake_x_2delta_plus_err(i) = Drate_uptake_x_2delta_max - Drate_2delta_h_inf;

    % and 2Delta Minus error
    Drate_uptake_x_2delta_minus_err(i) = Drate_2delta_h_inf - Drate_uptake_x_2delta_min;

    % Now, lets see if the plus and minus error are in fact equal to one
    % significant digit.  MATLAB does not have a trivial way to accomplish
    % this, so the following is the work around:

    % First, define a string for the plus error output:
    str_plus=sprintf('%.0d',Drate_uptake_x_2delta_plus_err(i));
    % Second, define a string for the plus error output:
    str_minus=sprintf('%.0d',Drate_uptake_x_2delta_minus_err(i));
    % The '.0d' rounds the values to the nearest whole integer in e -
    % notation. For example, 2.13e-005 becomes 2e-005.

    % Convert strings to double precision numbers:
    str_plus_num = str2double(str_plus);
    str_minus_num = str2double(str_minus);

    % Now, let's create a variable that will compare the first value in
    % str_plus and str_minus and see if they are the same:

    str_error_test=strncmp(str_plus,str_minus,1);
    % If this test is true, we only need to print one of the uncertainty
    % values (either the plus or the minus) as they are equal to one
    % significant digit.  If this is false, we will print out both values.

    % Now, lets print the result!

    fprintf(1,'The average dose rate to the %s cells at \n',Cell_name)
    fprintf(1,'the initiation of the uptake phase is:\n** %d', Drate_uptake_avg(i))

        if (str_error_test == 1)
            fprintf(1,' +/- %.0d', str_plus_num)
            fprintf(1,' Gy/h ** [1]\n\n')
        elseif (str_error_test == 0)
            fprintf(1,' +/- (%.0d/%.0d)',str_plus_num,str_minus_num)
            fprintf(1,' Gy/h ** [1]\n\n')
        end

    fprintf(1,'~~~~~~~\n\n')
end

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

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~Calculate Dose Rate in Slab~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

Drate_x=@(x)((A.*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));

% Again, 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)((A.*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

% 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(i) = (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 and slab thickness is smallest.
    % Conversely, minimum Drate occurs when Aconc is lowest and slab
    % thickness is largest. Therefore...

    % Maximum Drate:
    Drate_x_max=@(x)((A.*(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)((A.*(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(i) = (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)((A.*(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)((A.*(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(i) = (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(i) = Drate_slab_avg_max(i) - Drate_slab_avg(i);
Drate_slab_avg_minus_err(i) = Drate_slab_avg(i) - Drate_slab_avg_min(i);

    % First, define a string for the plus error output:
    str_plus=sprintf('%.0d',Drate_slab_avg_plus_err(i));
    % Second, define a string for the plus error output:
    str_minus=sprintf('%.0d',Drate_slab_avg_minus_err(i));
    % The '.0d' rounds the values to the nearest whole integer in e -
    % notation. For example, 2.13e-005 becomes 2e-005.

    % Convert strings to double precision numbers:
    str_plus_num = str2double(str_plus);
    str_minus_num = str2double(str_minus);

    % Now, let's create a variable that will compare the first value in
    % str_plus and str_minus and see if they are the same:

    str_error_test=strncmp(str_plus,str_minus,1);
    % If this test is true, we only need to print on of the uncertainty
    % values (either the plus or the minus) as they equal to the number of
    % significant digits allotted for uncertainty values.  If this is
    % false, we will print out both values.

% Now, lets print the result!

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(i))

    if (str_error_test == 1)
        fprintf(1,' +/- %.0d', str_plus_num)
        fprintf(1,' Gy/h ** [2]\n\n')
    elseif (str_error_test == 0)
        fprintf(1,' +/- (%.0d/%.0d)',str_plus_num,str_minus_num)
        fprintf(1,' Gy/h ** [2]\n\n')
    end

fprintf(1,'~~~~~~~\n\n')

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

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~Calculate the Total Dose to the Cells~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

Dose(i) = ((Drate_uptake_avg(i)+Drate_slab_avg(i))/2)*(2)+(Drate_slab_avg(i))*(1);

    % Now calculate the uncertainties.

    Dose_err_plus(i) = (sqrt(((Drate_uptake_x_2delta_plus_err(i))^2+(Drate_slab_avg_plus_err(i))^2))/2)*(2)+Drate_slab_avg_plus_err(i)*(1);
    Dose_err_minus(i) = (sqrt(((Drate_uptake_x_2delta_minus_err(i))^2+(Drate_slab_avg_minus_err(i))^2))/2)*(2)+Drate_slab_avg_minus_err(i)*(1);

    % First, define a string for the plus error output:
    str_plus=sprintf('%.0d',Dose_err_plus(i));
    % Second, define a string for the plus error output:
    str_minus=sprintf('%.0d',Dose_err_minus(i));
    % See MATLAB Help 'sprintf' for function parameter explanation. The
    % '.0d' rounds the values to the nearest whole integer in e - notation.
    % For example, 2.13e-005 becomes 2e-005.

    % Convert strings to double precision numbers:
    str_plus_num = str2double(str_plus);
    str_minus_num = str2double(str_minus);

    % Now, let's create a variable that will compare the first value in
    % str_plus and str_minus and see if they are the same:

    str_error_test=strncmp(str_plus,str_minus,1);

% Now, lets print the result!

fprintf(1,'Therefore, the dose to the %s cells after\n',Cell_name)
fprintf(1,'treatment with %.1d MBq/ml of radiopharmaceutical with\n',A)
fprintf(1,'%.0d hour(s) for uptake and %.0d hour(s) for accumulation is:\n\n',Time_up,Time_accum)

    if (str_error_test == 1)
        fprintf(1,'***********************************\n')
        fprintf(1,'*** %d',Dose(i))
        fprintf(1,' +/- %.0d', str_plus_num)
        fprintf(1,' Gy ***\n')
        fprintf(1,'***********************************\n')
    elseif (str_error_test == 0)
        fprintf(1,'********************************************\n')
        fprintf(1,'*** %d',Dose(i))
        fprintf(1,' +/- (%.0d/%.0d)',str_plus_num,str_minus_num)
        fprintf(1,' Gy ***\n')
        fprintf(1,'********************************************\n')
    end

fprintf(1,'______________________________________________________\n')
fprintf(1,'======================================================\n')

    end % End of 'if' statement if the user only enter a single value
end % End of 'if' statement if the user enters a range of values

% Now, if the user entered 'Yes' to whether they wanted to save the results
% to a Microsoft Excel file, the following will create the file.

if strcmp(save,'Yes') == 1
    filename = char(strcat(answer_filename,format_choices(excel_format)));
    headings = {'Activity'
                'Uptake Dose Rate (Gy/h)'
                'Plus Error Uptake Dose Rate'
                'Minus Error Uptake Dose Rate'
                'Accum. Dose Rate (Gy/h)'
                'Plus Error Accum. Dose Rate'
                'Minus Error Accum. Dose Rate'
                'Absorbed Dose (Gy)'
                'Plus Error Absorbed Dose'
                'Minus Error Absorbed Dose'}';
   results_num = [Activity;
                   Drate_uptake_avg;
                   Drate_uptake_x_2delta_plus_err;
                   Drate_uptake_x_2delta_minus_err;
                   Drate_slab_avg;
                   Drate_slab_avg_plus_err;
                   Drate_slab_avg_minus_err;
                   Dose;
                   Dose_err_plus;
                   Dose_err_minus]';
    results = num2cell(results_num);
    output = [headings; results];
    xlswrite(filename, output, 'Sheet1', 'A1'); % Write the results
%     If you want a title on the Excel output, uncomment the following but
%     it will be a 1 - 1.5 sec performance hit. Remember to change the cell
%     output for the results 'output' above from 'A1' to something else.
%     title_string = sprintf('Cell line = %s',Cell_name);
%     title = {title_string}; Converts previous string to a single element
%     xlswrite(filename, title, 'Sheet1', 'A1');
end

% The following is a completion dialogue box that allows the user to either
% redo the calculations or exit the application.

exit_text = {'Calculations are complete. Would you like to do another series?'
             ' '
             ''' Yes '' -  Restarts program'
             ''' No ''  -  Exits program immediately'
             ' '};
exit_title = 'Again?';
str1_exit = 'Yes & Continue';
str2_exit = 'No & Exit';
options.Resize='on';
options.Interpreter='tex';
options.Default=str1_exit;
exit_note = questdlg(exit_text, exit_title, str1_exit, str2_exit, options);
end % for 'while' loop initiating the exit test function

% *************************** END OF PROGRAM *************************** %

