% ##########
% NOTES:
%    Before clicking Run, please make sure to enter values for the following parameters
%    into the 'Program input' field on the top-right:
%        Conductivity [S/m]
%        Permittivity
%        Frequency [MHz]
%        Volume [ml]
%    This is a version of the phantom recipe generator found at:
%        https: / / amri.ninds.nih.gov / phantomrecipe.html
%
% CREATED:
%    2025-03-27..04-21, JAdZ
%
% LICENSE:
%    Copyright (C) <= 2025 NIH, NYU and others
%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, see <http://[Log in to view URL]
%
% REFERENCE:
%    C Ianniello, JA de Zwart, Q Duan, CM Deniz, L Alon, JS Lee, R Lattanzi, R Brown
%    Synthesized tissue-equivalent dielectric phantoms using salt and polyvinylpyrrolidone
%    solutions. Magn Reson Med 2018, 80:413 (doi: 10.1002/mrm.27005)
%
% MODIFICATION HISTORY:
% ##########




% ##### make this a script file #####
% see https://[Log in to view URL]
1;



% ##### PVP_fitting() #####
function [M_NaCl,M_PVP, M_water] = PVP_fitting(sig,epi_r, f, M_water)

if length(M_water(:)) == 1
    M_water = M_water * ones(size(sig));
end

M_NaCl = zeros(size(sig));
M_PVP = M_NaCl;

for k = 1:length(sig(:))
    [PVP_ratio, Salt_ratio] = solve_PVP_poly(sig(k),epi_r(k),f(k));
    M_NaCl(k) = Salt_ratio*M_water(k);
    M_PVP(k) = PVP_ratio*M_water(k);
end

endfunction



% ##### solve_PVP_poly() #####
function [PVP_ratio, Salt_ratio] = solve_PVP_poly(si,ei,fi)
tmp = fminsearch(@(x) err_fun(x,fi,si,ei), [0.65 0.0125]);
PVP_ratio = tmp(1);
Salt_ratio =tmp(2);
end

function errs = err_fun(x,fi,si,ei)
 [s_fit,e_fit] = PVP_get_property(x(1), x(2), fi);
x1 = [ei,si];
x2 = [e_fit,s_fit];

diffs = (x1-x2)./x1;
errs = sqrt(sum(diffs.^2));

endfunction



% ##### ND_polyval() #####
function [Y,EM,Porder] = ND_polyval(p,X,varargin)

% 2016/08/31, Jacco A. de Zwart
%
% Copied from Qi Duan's SANBA_get_recipe_togo.m, since it is also needed in
% PVP_get_recipe.m

if(~isvector(p))
    error('p must be a vector');
end
p=p(:);
Porder=[];
PM=[];

Nterm = length(p);
[Nsample,Ndim]=size(X);

switch nargin
    case 2
        
        %first, we need to figure out the order
        Porder = ND_polyFindPorder(Nterm,Ndim);
        
        %calculate the power matrix
        PM = flipud(perms_sum([0 Porder],Ndim));
    case 3
        PM = varargin{1};
        if (length(p)~=size(PM,1))
            error('p and PM are not compatible');
        end
        if (size(X,2)~=size(PM,2))
            error('X and PM are not compatible');
        end
       
    otherwise
        error('Incorrect number of input variables');
        
end

buff = ones(Nterm,Ndim,Nsample);

%alternatively, presumable faster since usually Nsample>Nterm>Ndim
for l = 1:Ndim
    dmatrix = X(:,l);
    for k = 1:Nterm
        pwr = dmatrix.^PM(k,l);
        buff(k,l,:) = pwr;
    end
end
%Y = sum(squeeze(prod(buff,2)).*repmat(p(:),[1,Nsample]),1);

%compute encoding matrix
EM = squeeze(prod(buff,2)); %Nterm x Nsample
EM = EM'; %Nsample x Nterm
Y = EM*p;
Y = Y(:);

endfunction


% ##### PVP_get_property() #####
function [sig,epi_r] = PVP_get_property(PVP_Water_Ratio, Salt_Water_Ratio, Freq)
if (length(PVP_Water_Ratio)~=length(Salt_Water_Ratio))||(length(PVP_Water_Ratio)~=length(Freq))
    error(' PVP_get_property: Three inputs have to be the same size');
end

% Cannot load S from anywhere, hard-code locally instead
%S = load('PVP_poly');
S.PM_e=[2 2 3; 2 1 4; 2 0 5; 1 2 4; 1 1 5; 1 0 6; 0 2 5; 0 1 6; 0 0 7; 2 2 2; ...
    2 1 3; 2 0 4; 1 2 3; 1 1 4; 1 0 5; 0 2 4; 0 1 5; 0 0 6; 2 2 1; 2 1 2; ...
    2 0 3; 1 2 2; 1 1 3; 1 0 4; 0 2 3; 0 1 4; 0 0 5; 2 2 0; 2 1 1; 2 0 2; ...
    1 2 1; 1 1 2; 1 0 3; 0 2 2; 0 1 3; 0 0 4; 2 1 0; 2 0 1; 1 2 0; 1 1 1; ...
    1 0 2; 0 2 1; 0 1 2; 0 0 3; 2 0 0; 1 1 0; 1 0 1; 0 2 0; 0 1 1; 0 0 2; ...
    1 0 0; 0 1 0; 0 0 1; 0 0 0];
S.PM_s=[2 2 3; 2 1 4; 2 0 5; 1 2 4; 1 1 5; 1 0 6; 0 2 5; 0 1 6; 0 0 7; 2 2 2; ...
	2 1 3; 2 0 4; 1 2 3; 1 1 4; 1 0 5; 0 2 4; 0 1 5; 0 0 6; 2 2 1; 2 1 2; ...
	2 0 3; 1 2 2; 1 1 3; 1 0 4; 0 2 3; 0 1 4; 0 0 5; 2 2 0; 2 1 1; 2 0 2; ...
	1 2 1; 1 1 2; 1 0 3; 0 2 2; 0 1 3; 0 0 4; 2 1 0; 2 0 1; 1 2 0; 1 1 1; ...
	1 0 2; 0 2 1; 0 1 2; 0 0 3; 2 0 0; 1 1 0; 1 0 1; 0 2 0; 0 1 1; 0 0 2; ...
	1 0 0; 0 1 0; 0 0 1; 0 0 0];
S.PVP_scale=1.;
S.f_scale=1.e-9;
S.p_e=transpose([ 0.011805335 ; -0.14464843 ; 0.0069357918 ; 0.013468636 ; -0.0059084284 ; 0.010146790 ; ...
    0.00078412937 ; 0.0040864598 ; -0.015136267 ; -0.095121764 ; 1.5438998 ; -0.035230146 ; ...
    -0.16863082 ; 0.22784519 ; -0.17019150 ; -0.025323176 ; -0.069588966 ; 0.25411302 ; 0.28599051 ; ...
    -5.5439930 ; -0.14294256 ; 0.72280651 ; -1.9631792 ; 1.0517561 ; 0.22518252 ; 0.48922931 ; ...
    -1.7269313 ; -0.082226131 ; 8.3682848 ; 0.85457320 ; -1.3976375 ; 6.3753179 ; -2.9454405 ; ...
    -0.81507715 ; -1.8610658 ; 6.1196438 ; -4.7693457 ; 0.34965290 ; 0.51793472 ; -8.2129545 ; ...
    4.5307713 ; 1.3414729 ; 4.2276048 ; -12.137641 ; 15.144375 ; 5.5387159 ; -10.556753 ; -0.53304948 ; ...
    -5.8848859 ; 13.082945 ; -50.563847 ; 1.5125963 ; -7.1994831 ; 79.442701 ]);
S.p_s=transpose([ 0.00050292956 ; 0.0044025354 ; 0.00081702467 ; -0.00035627814 ; -0.00040857246 ; -0.00058451886 ; ...
    8.6110845e-06 ; -1.5125634e-05 ; -0.00041641555 ; -0.0042955366 ; -0.033546825 ; -0.012998685 ; ...
    0.0017412356 ; -0.00044430190 ; 0.0054871013 ; 0.00031705808 ; 0.00059590591 ; 0.0079475010 ; ...
    0.00098822067 ; 0.027044634 ; 0.077991047 ; 0.0052494236 ; 0.024306980 ; -0.0031549297 ; ...
    -0.0025991520 ; -0.0060862995 ; -0.058010820 ; -0.093563581 ; -0.11462494 ; -0.12628383 ; ...
    0.0086047419 ; -0.026771431 ; -0.13466063 ; 0.0019014825 ; 0.027826888 ; 0.20578417 ; 1.7840338 ; ...
    -0.0081733085 ; 0.15066118 ; 0.064693028 ; 0.29797058 ; -0.012923270 ; -0.083016310 ; -0.37726294 ; ...
    -0.38131759 ; -3.2915352 ; 0.22245810 ; -0.055920294 ; 0.14949451 ; 0.59238066 ; 0.31376929 ; ...
    1.5957381 ; -0.15722930 ; 0.050235394 ]);
S.salt_scale=100;

sig = zeros(size(Freq));
epi_r = sig;

x=[PVP_Water_Ratio(:),Salt_Water_Ratio(:)];

sig(:) = ND_polyval(S.p_s,[x(:,1)*S.PVP_scale,x(:,2)*S.salt_scale, Freq(:)*S.f_scale],S.PM_s);
epi_r(:) = ND_polyval(S.p_e,[x(:,1)*S.PVP_scale,x(:,2)*S.salt_scale, Freq(:)*S.f_scale],S.PM_e);

endfunction


% ##### MAIN #####
fprintf(1,"Make sure to enter values for the following in the 'Program input' box\n");
fprintf(1,"directly above this text (newline separated):\n");
fprintf(1,"    Conductivity [S/m]\n");
fprintf(1,"    Permittivity\n");
fprintf(1,"    Frequency [MHz]\n");
fprintf(1,"    Volume [ml]\n");
fprintf(1,"\n");
cond=input("Conductivity [S/m]: ")
perm=input("Permittivity      : ")
freq=input("Frequency [MHz]   : ")
volu=input("Volume [ml]       : ")
%[nacl,pvp,water] = PVP_fitting(0.6,50.,297400000.,1000.)
[nacl,pvp,water] = PVP_fitting(cond,perm,freq*1.e6,volu);

fprintf(1,"\n");
fprintf(1,"REQUIRED INGREDIENTS:\n")
fprintf(1,"NaCl: %f g\n",nacl);
fprintf(1,"PVP: %f g\n",pvp);
fprintf(1,"Water: %f g\n",water);
fprintf(1,"\n");