clear all

% Plot of technetium radioactivity from urine of a MIBI scan patient 
% downloaded 4/17/2020 from http://www.vk2zay.net/article/239.
%
% Executes properly on MATLAB version 9.4.0.813654 (R2018a).
% Executes mostly properly on OCTAVE version 4.2.1.
% Joshua C. Vaughan 4/17/2020

basepath = 'F:\Josh\CHEM452\build_files\MATLAB_scripts\data_files';
nm_in = [basepath '\MIBI_scan_radioactive_technetium_urine_data.txt'];

A = dlmread(nm_in)/1955;
% A shows radioactivity which is directly the rate of change but based on
% first-order kinetics is also proportional to the amount of radioactive
% species (technetium) present. A is divided here by the approximate
% maximum.

n = length(A);                      % Number of samples along time axis.
% Web page didn't provide time vector that goes along with the 
% radioactivity vector. It seems there are ~395 samples/hr based on the 
% known ~6 hour half-life of Technetium-99m.

t = [linspace(0,n/395,n)]';          % Time vector 1. Units are hours

% Fit data with an exponential function.
s = fitoptions('Method','NonlinearLeastSquares',...
    'Lower',[0.5 0],...                                 % parameter lower bound [A_min B_min]
    'Upper',[2 inf],...                                 % parameter upper bound [A_max B_max]
    'StartPoint',[1 log(2)/6.35]);                      % estimate of fit parameters
f = fittype('A*exp(-x*B)','options',s);                 % form of the fit (exponential, here)
FitStruct = fit(t,A,f);                               % do the actual fit

% Evaluate fitted function for plotting
A_fit = feval(FitStruct,t); % fitted function
% disp(['half life = ' num2str(log(2)/FitStruct.B) ' h'])


% create plot

ax_lw = 2;
pl_lw = 2;
ax_fs = 24;
la_fs = 24;
ti_fs = 24;
fontname = 'Arial';

figure(1); clf
plot(t,A,'k','linewidth',pl_lw); hold on
plot(t,A_fit,'r','linewidth',pl_lw); hold off
xlabel('time (h)','fontsize',la_fs)
ylabel('^{99m}Tc concentration (arb.)','fontsize',la_fs)
set(gca,'xtick',[0:10:40])
set(gca,'ytick',[0:0.2:1])
set(gcf,'color',[1 1 1])
set(gca,'linewidth',ax_lw)
set(gca,'fontsize',ax_fs)
set(gca,'tickdir','out')
set(gca,'fontname',fontname)
xlim([0 40])
ylim([0 1])
h = legend('radioactivity','exponential fit');
set(h,'box','off')
set(h,'fontsize',20)
set(h,'position',[0.5272    0.2757    0.3600    0.1330]);
title(['^{99m}Tc \rightarrow ^{99}Tc + \gamma'],'fontsize',ti_fs)
set(gca,'LooseInset',get(gca,'TightInset'));

% Load Sestamibi chemical structure image
img_nm_in = [basepath '\sestamibi_structure.png'];
S = imread(img_nm_in);

axes('Position',[.45 .45 .4 .4]);
imshow(S);
% Kinda pixelated but just good enough to use.
