function [Smw Fmw2] = dens2mstp(dens_ratio)
% DENS2MSTP Convert physical density to mass stopping power and Fmw2
%
% Usage:
%     [Smw Fmw2] = dens2mstp ( density_matrix )
% INPUT:
%     density_matrix = density matrix (g/cm^3) with any dimensionality
% OUTPUT:
%     Smw = stopping power ratio
%     Fmw2 = squared F value ratio to water
%
% Given an array of densities relative to water (assumed to be from a CT
% scan), the mass stopping power ratio (relative to water) array is
% calculated. The conversion is done by assuming that every voxel in the CT
% consists of either air, adipose tissue, muscle, or bone. A lookup table
% of mass stopping power ratios between each material and water is then
% constructed by computing the average mass stopping power ratio over some
% proton energy region, usually between 50 MeV and 250 MeV. The mass
% stopping powers are looked-up in tables obtained from the NIST website.
%
% This code was intended to be used as a way of converting CT
% data to mass stopping power data for proton dose calculations. The energy
% of the proton beam is assumed to be between 50 MeV and 250 MeV, a regime
% in which the approximation that stopping power ratios are energy
% independent holds.
%
% Note that the value returned is the stopping power ratio. The stopping power 
% has units of MeV/cm. This is an important distinction!
%
% Copyleft: Ryan Flynn, Xiaohu Mo
% Version: 1.2

% Set up the boundary points for the density-to-material lookup table.
air_low = 0;
air_high = 0.05; fat_low = air_high;
fat_high = 0.99; wat_low = fat_high;
wat_high = 1.01; mus_low = wat_high;
mus_high = 1.25; bon_low = mus_high;
bon_high = 20.0;

% check the input argument
if ~isempty(find(dens_ratio < air_low, 1))
    error('Cannot have negative densities.  Redefine input argument.');
end

if ~isempty(find(dens_ratio > bon_high, 1))
    error('Extremely high densities are present in the input array, data may not be scaled properly');
end

% Product of Fmw and rho for tabluated materials
fmw_air_to_water = 990.8 * 1.205E-3;
fmw_fat_to_water = 0.972 * 0.92;
fmw_mus_to_water = 0.972 * 1.04;
fmw_bon_to_water = 0.656 * 1.85;
fmw_wat_to_water = 1.0;

% Stopping power ranges that will be used to calculate the average stopping
% power ratios between each material and water.
Tlo = 1;      % lower energy boundary in MeV
Thi = 250;    % upper energy boundary in MeV

% load the stopping powers
load('MSTPliquid_water.mat');
water = MSTP;
load('MSTPair.mat');
air = MSTP;
load('MSTPadipose.mat');
fat = MSTP;
load('MSTPskeletal_muscle.mat');
mus = MSTP;
load('MSTPcompact_bone.mat');
bon = MSTP;

% calculate the spacing between energy bins
dT = [water.T(1) diff(water.T)];

% find the indices of integration that will be used
indices_low = find(water.T <= Tlo);
ind_low = indices_low(end);

indices_high = find(water.T >= Thi);
ind_high = indices_high(1);

% Calculate the average stopping power ratio for the material and water
% between the energy range specified above.
smw_air_to_water = sum(air.tot(ind_low:ind_high)./water.tot(ind_low:ind_high).*dT(ind_low:ind_high))/(Thi-Tlo);
smw_fat_to_water = sum(fat.tot(ind_low:ind_high)./water.tot(ind_low:ind_high).*dT(ind_low:ind_high))/(Thi-Tlo);
smw_wat_to_water = 1.0;
smw_mus_to_water = sum(mus.tot(ind_low:ind_high)./water.tot(ind_low:ind_high).*dT(ind_low:ind_high))/(Thi-Tlo);
smw_bon_to_water = sum(bon.tot(ind_low:ind_high)./water.tot(ind_low:ind_high).*dT(ind_low:ind_high))/(Thi-Tlo);

% material mask
region_air = find (dens_ratio >= air_low & dens_ratio < air_high);
region_fat = find (dens_ratio >= fat_low & dens_ratio < fat_high);
region_wat = find (dens_ratio >= wat_low & dens_ratio < wat_high);
region_mus = find (dens_ratio >= mus_low & dens_ratio < mus_high);
region_bon = find (dens_ratio >= bon_low);

% fill the mass stopping power ratio
Smw = zeros ( size ( dens_ratio ) );
Smw ( region_air ) = smw_air_to_water * dens_ratio ( region_air );
Smw ( region_fat ) = smw_fat_to_water * dens_ratio ( region_fat );
Smw ( region_wat ) = smw_wat_to_water * dens_ratio ( region_wat );
Smw ( region_mus ) = smw_mus_to_water * dens_ratio ( region_mus );
Smw ( region_bon ) = smw_bon_to_water * dens_ratio ( region_bon );

% fill the Fmw2
% divide by density to appropriately scale Fmw
Fmw2 = zeros ( size ( dens_ratio ) );
Fmw2 ( region_air ) = fmw_air_to_water ./ dens_ratio ( region_air );
Fmw2 ( region_fat ) = fmw_fat_to_water ./ dens_ratio ( region_fat );
Fmw2 ( region_wat ) = fmw_wat_to_water ./ dens_ratio ( region_wat );
Fmw2 ( region_mus ) = fmw_mus_to_water ./ dens_ratio ( region_mus );
Fmw2 ( region_bon ) = fmw_bon_to_water ./ dens_ratio ( region_bon );
Fmw2 = Fmw2.^2;


% % plot the stopping power ratio as a function of energy for bone and water
% f1 = figure; set(gcf,'color',[1 1 1]);
% p = plot(water.T(ind_low:ind_high),bone.tot(ind_low:ind_high)./water.tot(ind_low:ind_high));
% set(p,'LineWidth',2);
% set(gca,'YLim',[0.9 0.95]);
% set(gca,'FontSize',14);
% title('Mass stopping power ratio between bone and water','FontSize',14);
% xlabel('Energy (MeV)','FontSize',14);
% ylabel('Ratio','FontSize',14);
% print -dbitmap MSTP_ratio_with_energy.png;
% 
% % plot the density ratio to mass stopping power ratio lookup table:
% xmax = 2.5;
% N = 50;
% dx = xmax/N;
% x = [0:N-1]*dx;
% y = zeros(size(x));
% y(x >= air_low & x <= air_high) = air_to_water;
% y(x > adipose_low & x <= adipose_high) = adipose_to_water;
% y(x > muscle_low & x <= muscle_high) = muscle_to_water;
% y(x > bone_low) = bone_to_water;
% 
% z = y.*x;
% 
% f2 = figure; set(gcf,'color',[1 1 1]);
% plot(x(x >= air_low & x <= air_high),y(x >= air_low & x <= air_high),'.b',...
%     x(x > adipose_low & x <= adipose_high),y(x > adipose_low & x <= adipose_high),'+b',...
%     x(x > muscle_low & x <= muscle_high),y(x > muscle_low & x <= muscle_high),'xb',...
%     x(x > bone_low),y(x > bone_low),'*b');
% l = legend('air','adipose','muscle','bone');
% set(l,'FontSize',14);
% set(gca,'FontSize',14);
% xlabel('Density ratio','FontSize',14);
% ylabel('Mass stopping power ratio','FontSize',14);
% title('Change in mass stopping power ratio with density','FontSize',14);
% print -dbitmap MSTP_ratio_with_density.png;
% 
% f3 = figure; set(gcf,'color',[1 1 1]);
% plot(x(x >= air_low & x <= air_high),z(x >= air_low & x <= air_high),'.b',...
%     x(x > adipose_low & x <= adipose_high),z(x > adipose_low & x <= adipose_high),'+b',...
%     x(x > muscle_low & x <= muscle_high),z(x > muscle_low & x <= muscle_high),'xb',...
%     x(x > bone_low),z(x > bone_low),'*b');
% l = legend('air','adipose','muscle','bone');
% set(l,'FontSize',14,'position',[0.6258 0.2433 0.2250 0.2460]);
% set(gca,'FontSize',14);
% xlabel('Density ratio','FontSize',14);
% ylabel('Stopping power ratio','FontSize',14);
% title('Change in stopping power ratio with density','FontSize',14);
% print -dbitmap STP_ratio_with_density.png;