% analyze.m
% updated Apr 18, 2008 by SLJ (so runs on current MATLAB)
% Uses
%	lookskin1a.m
%		getSpectrum.m
%		skin1.Master.Scope
%		std1.Master.Scope
%	loadTakataniGraham.m
%	loadwaterHaleQuerry.m
%	fitspectrum.m
%	getRr.m

clear
global cnt
cnt = 0;

disp('analyze skin1')

lookskin1  % loads spectrum --> nm, skin, std, R

%%%%%%
% Load LIBRARY of spectra
%%%%
loadTakataniGraham  % hemoglobin (oxy, deoxy)
muaoxy = interp1(nmTG, muaoxyTG, nm);
muadeoxy = interp1(nmTG, muadeoxyTG, nm);
clear nmTH muaoxyTG muadeoxyTG
%
loadwaterHaleQuerry  % water
muawater = interp1(waterHaleQuerry(:,1), waterHaleQuerry(:,2), nm);
clear waterHaleQuerry
%%%%%%%
drawnow

r = 0.20; % cm = separation of fibers

i = 1; while nm(i)<460; i = i+1; end; i400 = i;
i = 1; while nm(i)<800; i = i+1; end; i900 = i;
u = (i400:i900);

%%%%%%
% fit data
%%%%
W = 0.65; % water content

B     = 0.001;
S     = 0.655;
a     = 1.0;
P     = 0.85;
Pwidth = 50;
Pnm = 460;
GG = 0.50;
start = [B S a P Pwidth Pnm GG];

result = fminsearch('fitBS', start, [], nm(u), R(u), muaoxy(u), muadeoxy(u), muawater(u));

B     = result(1)
S     = result(2)
a     = result(3)
P     = result(4);
Pwidth = result(5);
Pnm   = result(6);
GG    = result(7);

%%%%%%%%
% plot the data with fits
% P = porphyrins, i.e. bilirubin and possible PPIX
%%%
muaP = exp(-(abs(nm - Pnm)/Pwidth).^2);

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

r = 0.2;
n = 1.37;
mua   = B*(S*muaoxy + (1-S)*muadeoxy) + W*muawater + P*muaP;
pRwBb = GG*getRr(mua, musp, r, n); % +Blood +bili

mua  = B*(S*muaoxy + (1-S)*muadeoxy) + W*muawater;
pRwB = GG*getRr(mua, musp, r, n); % +Blood -bili

mua = W*muawater;
pRw = GG*getRr(mua, musp, r, n); % -Blood -bili

figure(2);clf
set(figure(2),'position',[645    40   577   757],'color','w')
sz = 18;
subplot(2,1,1); hold off
ymax = 1;
plot(nm(u), R(u), 'm-','linewidth',2)
hold on
plot(nm, pRwBb, 'b-','linewidth',2)
plot(nm, pRwB, 'r-','linewidth',2)
plot(nm, pRw, 'k-','linewidth',2)
set(gca,'fontsize',sz,'linewidth',2)
xlabel('wavelength [nm]')
ylabel('T = Tstd GG M/Mstd')
axis([400 900 0 1])
text(450, ymax*0.9, sprintf('B = %0.5f', B),'fontsize',sz)
text(450, ymax*0.8, sprintf('S = %0.3f', S),'fontsize',sz)
text(450, ymax*0.7, sprintf('a = %0.3f', a),'fontsize',sz)
text(450, ymax*0.6, sprintf('P = %0.3f', P),'fontsize',sz)
text(450, ymax*0.5, sprintf('widthP = %0.1f', Pwidth),'fontsize',14)
text(450, ymax*0.4, sprintf('nmP = %0.1f', Pnm),'fontsize',14)
text(450, ymax*0.3, sprintf('GG = %0.2f', GG),'fontsize',sz)

subplot(2,1,2); hold off
semilogy(nm, W*muawater, 'k-','linewidth',2)
hold on
plot(nm, B*S*muaoxy, 'r-','linewidth',2)
plot(nm, B*(1-S)*muadeoxy, 'b-','linewidth',2)
plot(nm, P*muaP, 'm-','linewidth',2)
text(470,P*1.1,'P','fontsize',sz,'color','m')
text(580,0.3,'oxy','fontsize',sz,'color','r')
text(500,0.04,'deoxy','fontsize',sz,'color','b')
text(750,0.025,'water','fontsize',sz,'color','k')
set(gca,'fontsize',sz,'linewidth',2)
xlabel('wavelength [nm]')
ylabel('OD = -log_e(M)')
axis([400 900 1e-3 1])

figure(3);clf
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,160,sprintf('\mu_s'' = %0.2f(Mie + Ray)', a),'fontsize',sz,'color','r')
text(600,130,'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('musp [cm^-^1]')
axis([400 900 0 200])
%
subplot(2,1,2); hold off
nm2 = 100*2.^(1:50)';
Mie2 = 2e5*nm2.^-1.5;
Ray2 = 2e12*nm2.^-4;
musp2 = a*(Mie2 + Ray2);
loglog(nm2, musp2,'r-','linewidth',2,'color','r')
hold on
plot(nm2, a*Mie2,'b-','linewidth',2,'color','b')
plot(nm2, a*Ray2,'k-','linewidth',2,'color','k')
text(600,400,sprintf('\mu_s'' = %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('musp [cm^-^1]')
axis([100 10000 1 1000])