% analyze2.m
% Uses
%	lookskin1.m
%		getSpectrum.m
%		skin1.Master.Scope
%		std1.Master.Scope
%	loadTakataniGraham.m
%	loadwaterHaleQuerry.m
%	fitspectrum2.m
%	getRr.m

clear
global cnt
cnt = 0;

disp('analyze skin1')


lookskin1  % loads spectrum --> nm, skin, std, R
nm0 = nm;

%%%%%%
% Load LIBRARY of spectra
%%%%
loadTakataniGraham  % hemoglobin (oxy, deoxy)
%
loadwaterHaleQuerry  % water
nmHQ = waterHaleQuerry(:,1);
muaHQ = waterHaleQuerry(:,2)
clear waterHaleQuerry
%%%%%%%
drawnow

r = 0.20; % cm = separation of fibers

i = 1; while nm(i)<450; i = i+1; end; i450 = i;
i = 1; while nm(i)<800; i = i+1; end; i800 = i;
u = (i450:i800);
i=1; while nm(u(i))<500, i=i+1;end; u500 = i;while nm(u(i))<600, i=i+1;end; u600 = i;
v = (u500:u600); % for use as nm(u(v)) in emphasizing fit to the 500-600 nm range, to specify oxy/deoxy hemoglobin


%%%%%%
% fit data
%%%%
M = Mskin./Mstd;  % repeating what was done in lookskin1.m

W = 0.65; % water content

B      = 0.001; % mean blood volume fraction in skin
S      = 0.655; % tissue oxygen saturation (mixture of arterial and venous blood)
a      = 1.0;   % scale the skin scattering (Rayleigh + Mie)
Mel    = 0.05;  % aveage volume fraction of melanosomes in 60-um-thick epidermis
const  = 0.50;  % const = Rstd*Gstd/Gskin
nmOff  = 0;     % correct for any offset error in spectrometer wavelength registration
start  = [B S a Mel const nmOff];

% It sometimes helps to run the fitting several times, and perturb the
% blood content a little to force a readjustment.
start = fminsearch('fitspectrum2', start, [], nm(u), M(u), nmTG, muaoxyTG, muadeoxyTG, nmHQ, muaHQ, u, v)
start = fminsearch('fitspectrum2', start, [], nm(u), M(u), nmTG, muaoxyTG, muadeoxyTG, nmHQ, muaHQ, u, v)
start(1) = start(1)/2;
result = fminsearch('fitspectrum2', start, [], nm(u), M(u), nmTG, muaoxyTG, muadeoxyTG, nmHQ, muaHQ, u, v)

B     = result(1);
S     = result(2);
a     = result(3);
Mel   = result(4);
const = result(5);
nmOff = result(6);

%%%%%%%%
% plot the data with fits
%%%

%%%%%
% model
%%
nm = nm0 + nmOff;
muaoxy = interp1(nmTG, muaoxyTG, nm);
muadeoxy = interp1(nmTG, muadeoxyTG, nm);
muawater = interp1(nmHQ, muaHQ, nm);
muamel = 6.6e11*nm.^-3.33;

W = 0.65; % water content
%b = 1.0;  % scattering coefficient, musp = a*nm.^-b

Mie = 4.59e3*nm.^-0.913;
Ray = 1.74e12*nm.^-4;
musp = a*(Mie + Ray);

r = 0.2;
n = 1.4;

R = const*M;  % R = (Mskin/Mstd)(Rstd*Gstd/Gskin)

% epidermal melanin
Lepi = 2*0.0060;
Tepi = exp(-Mel*muamel*Lepi);

% +scatter only
mua = 1e-6;
pRS = const*getRr(mua, musp, r, n); 

% +scatter +Water
mua = W*muawater;
pRw  = const*getRr(mua, musp, r, n); 

% +scatter +Water +Mel
pRwM = const*Tepi.*getRr(mua, musp, r, n); 

% +scatter +Water +Mel +Blood
mua  = B*(S*muaoxy + (1-S)*muadeoxy) + W*muawater;
pRwMB = const*Tepi.*getRr(mua, musp, r, n); 
%%
% end model
%%%%%


figure(2);clf
set(figure(2),'position',[645    40   577   757],'color','w')
sz = 18;
subplot(2,1,1); hold off
plot(nm(u), M(u), 'ko','linewidth',1)
hold on
plot(nm, pRS, 'k--','linewidth',2)
plot(nm, pRw, 'b-','linewidth',2)
plot(nm, pRwM, 'k-','linewidth',2)
plot(nm, pRwMB, 'r-','linewidth',2)
%
set(gca,'fontsize',sz,'linewidth',2)
xlabel('wavelength [nm]')
ylabel('M = M_s_k_i_n / M_s_t_d')
title('M = const T_e_p_i getRr(\mu_a, \mu_s'', r, n)')
axis([400 900 0 1])
x = 410; ymax = 1; dy = .07;
text(x, ymax - dy, sprintf('B = %0.4f', B),'fontsize',sz)
text(x, ymax - 2*dy, sprintf('S = %0.3f', S),'fontsize',sz)
text(x, ymax - 3*dy, sprintf('a = %0.3f', a),'fontsize',sz)
text(x, ymax - 4*dy, sprintf('const = %0.2f', const),'fontsize',sz)
text(x, ymax - 5*dy, sprintf('Mel = %0.3f', Mel),'fontsize',sz)
text(x, ymax - 6*dy, sprintf('nmOff = %0.3f',nmOff),'fontsize',sz)

subplot(2,1,2); hold off
semilogy(nm, W*muawater, 'b--','linewidth',2)
hold on
plot(nm, B*S*muaoxy, 'r-','linewidth',2)
plot(nm, B*(1-S)*muadeoxy, 'b-','linewidth',2)
plot(nm, Mel*muamel, 'k-','linewidth',2)
plot(nm, B*(S*muaoxy+(1-S)*muadeoxy)+W*muawater,'m-','linewidth',2)
plot(nm, Tepi,'k-','linewidth',2)
text(630, 1.4e-2,'total','fontsize',sz,'color','m')
text(480,0.07,'oxy','fontsize',sz,'color','r')
text(480,0.015,'deoxy','fontsize',sz,'color','b')
text(750,0.01,'water','fontsize',sz,'color','b')
text(650,0.5,'T.epidermal.melanin','fontsize',sz,'color','k')
set(gca,'fontsize',sz,'linewidth',2)
xlabel('wavelength [nm]')
ylabel('\mu_a [cm^-^1]')
axis([400 900 1e-3 1])

figure(3);clf
set(figure(3),'position',[300    10   577   757],'color','w')
subplot(2,1,1); hold off
plot(nm, musp,'r-','linewidth',2)
hold on
plot(nm, a*Mie,'b-','linewidth',2)
plot(nm, a*Ray,'k-','linewidth',2)
text(600,140,['\mu_s''' sprintf(' = %0.2f(Mie + Ray)', a)],'fontsize',sz,'color','r')
text(600,120,'Mie = 4.x10^5 nm^-^1^.^5','fontsize',sz,'color','b')
text(600,100,'Ray = 2x10^1^2 nm^-^4','fontsize',sz,'color','k')
set(gca,'fontsize',sz)
xlabel('wavelength [nm]')
ylabel('\mu_s'' [cm^-^1]')
axis([400 900 0 150])
%
subplot(2,1,2); hold off
%
nm2 = (100:1e4)';
Mie2 = 4.59e3*nm2.^-0.913;
Ray2 = 1.74e12*nm2.^-4;
musp2 = a*(Mie2 + Ray2);
%
loglog(nm, musp,'r-','linewidth',4,'color','r')
hold on
plot(nm, a*Mie,'b-','linewidth',4,'color','b')
plot(nm, a*Ray,'k-','linewidth',4,'color','k')
%
loglog(nm2, musp2,'r--','linewidth',2,'color','r')
plot(nm2, a*Mie2,'b--','linewidth',2,'color','b')
plot(nm2, a*Ray2,'k--','linewidth',2,'color','k')
%
text(600,400,['\mu_s''' sprintf(' = %0.2f(Mie + Ray)', a)],'fontsize',sz,'color','r')
text(600,200,'Mie = 4.x10^5 nm^-^1^.^5','fontsize',sz,'color','b')
text(600,100,'Ray = 2x10^1^2 nm^-^4','fontsize',sz,'color','k')
set(gca,'fontsize',sz)
xlabel('wavelength [nm]')
ylabel('\mu_s'' [cm^-^1]')
axis([100 1e4 1 1e4])