% ##########
% 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]
% Agar [% w/v]
% Benzoic acid [% w/v]
% Temperature [degC]
% This is a version of the phantom recipe generator found at:
% https: / / amri.ninds.nih.gov / phantomrecipe.html
%
% CREATED:
% 2025-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:
% Q Duan, JH Duyn, N Gudino, JA de Zwart, P van Gelderen, DK Sodickson, R Brown
% Characterization of a dielectric phantom for high-field magnetic resonance
% imaging applications. Med Phys 2014, 10:102303 (doi: 10.1118/1.4895823)
%
% MODIFICATION HISTORY:
% ##########
% ##### make this a script file #####
% see https://[Log in to view URL]
1;
function k_est = SANBA_Surcrose_Solution_Thermal_Conductivity(C,B)
% function k_est = SANBA_Surcrose_Solution_Thermal_Conductivity(C,B)
% Compute thermal conductivity for sucrose solution based on
% Proceedings of Tlze Soutlz African Sugar Technologists' Association -
% March 1966
% http://[Log in to view URL]
% C in C
% B = mass_sucrose/mass_all*100;
%
% version 0.1 Qi Duan 2016-03-07
% - creation
F = C*9/5+32;
k = (3.61e-4*F - 1.96e-3*B + 0.322) ; % in Btu(ft)/(sq ft)(hr)(F)
k_est = k*1.730734666; % in W/(m*K)
endfunction
function rou = SANBA_Surcrose_Solution_Density(DS)
p=[
-1.333170009997885e+05
7.745944445120688e+05
-1.953331667854245e+06
2.793537951018787e+06
-2.478692164622070e+06
1.395021682455995e+06
-4.808191248372806e+05
8.614918269335467e+04
-9.037724191534653e+02
-2.634624673009943e+03
5.885069691102377e+02
3.608707126742347e+02
1.000149435309792e+03
]';
rou=polyval(p,DS);
endfunction
function [M_NaCl,M_Sugar, M_Agar, M_BA, M_water] = SANBA_fitting(sig,epi_r, f, M_water, Agar_Percent, BA_Percent,X0)
if length(M_water(:)) == 1
M_water = M_water * ones(size(sig));
end
if length(Agar_Percent(:)) == 1
Agar_Percent = Agar_Percent * ones(size(sig));
end
if length(BA_Percent(:)) == 1
BA_Percent = BA_Percent * ones(size(sig));
end
M_Agar = M_water.*Agar_Percent;
M_BA = M_water.*BA_Percent;
M_NaCl = zeros(size(sig));
M_Sugar = M_NaCl;
for k = 1:length(sig(:))
[Sugar_ratio, Salt_ratio] = solve_SANBA_poly(sig(k),epi_r(k),f(k),X0(k,:));
M_NaCl(k) = Salt_ratio*M_water(k);
M_Sugar(k) = Sugar_ratio*M_water(k);
end
endfunction
function [Sugar_ratio, Salt_ratio] = solve_SANBA_poly(si,ei,fi,varargin)
X0 = [];
if nargin > 3
X0 = varargin{1};
end
if isempty(X0)
X0 = [0.65 0.0125];
end
%[tmp, fval] = fminsearch(@(x) err_fun(x,fi,si,ei), X0, [0 1e-4]); %this is for compatibility to octave. Seems working for version 0.3
%[tmp, fval] = fminsearch(@(x) err_fun(x,fi,si,ei), X0, optimset('TolX',1e-5));
%[tmp, fval] = fminsearch(@(x) err_fun(x,fi,si,ei), X0);
x = build_initial_simplex(X0);
[ tmp, n_feval ] = nelder_mead ( x,@(x) err_fun(x,fi,si,ei),0);
Sugar_ratio = tmp(1);
Salt_ratio =tmp(2);
endfunction
function errs = err_fun(x,fi,si,ei)
[s_fit,e_fit] = SANBA_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
function [sig,epi_r] = SANBA_get_property(Sugar_Water_Ratio, Salt_Water_Ratio, Freq)
%simple check
if (length(Sugar_Water_Ratio)~=length(Salt_Water_Ratio))||(length(Sugar_Water_Ratio)~=length(Freq))
error(' SANBA_get_property: Three inputs have to be the same size');
end
%S = load('SANBA_poly');
% data 20131107
S.f_scale = 1.000000000000000e-09;
S.sugar_scale = 1;
S.salt_scale = 100;
S.PM_e = [ 7 0 0
6 1 0
6 0 1
5 2 0
5 1 1
5 0 2
4 3 0
4 2 1
4 1 2
4 0 3
3 4 0
3 3 1
3 2 2
3 1 3
3 0 4
2 5 0
2 4 1
2 3 2
2 2 3
2 1 4
2 0 5
1 5 1
1 4 2
1 3 3
1 2 4
1 1 5
1 0 6
0 5 2
0 4 3
0 3 4
0 2 5
0 1 6
0 0 7
6 0 0
5 1 0
5 0 1
4 2 0
4 1 1
4 0 2
3 3 0
3 2 1
3 1 2
3 0 3
2 4 0
2 3 1
2 2 2
2 1 3
2 0 4
1 5 0
1 4 1
1 3 2
1 2 3
1 1 4
1 0 5
0 5 1
0 4 2
0 3 3
0 2 4
0 1 5
0 0 6
5 0 0
4 1 0
4 0 1
3 2 0
3 1 1
3 0 2
2 3 0
2 2 1
2 1 2
2 0 3
1 4 0
1 3 1
1 2 2
1 1 3
1 0 4
0 5 0
0 4 1
0 3 2
0 2 3
0 1 4
0 0 5
4 0 0
3 1 0
3 0 1
2 2 0
2 1 1
2 0 2
1 3 0
1 2 1
1 1 2
1 0 3
0 4 0
0 3 1
0 2 2
0 1 3
0 0 4
3 0 0
2 1 0
2 0 1
1 2 0
1 1 1
1 0 2
0 3 0
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.p_e =[
5.225569560526150e+00
2.455521740809368e-03
-4.009788611906068e-01
-5.373088464388277e-03
-6.499320303091425e-04
1.519227255273204e-01
-1.223174709494335e-03
-3.869093673445125e-04
-1.103257602639135e-02
2.009609707621798e-03
-1.159352150287421e-04
1.378476762903892e-04
-1.302665393315307e-03
-9.296005123066742e-03
-6.822483659803738e-02
-9.366869685483423e-07
7.142942536206241e-06
4.758884122555857e-05
-5.819640363357153e-04
-4.933838377608522e-03
-3.273200434820046e-02
5.634404479876636e-07
4.717829082704831e-06
3.757907888892905e-05
-9.892824107051136e-06
-1.071108822170348e-03
3.242326024303805e-03
-1.285096024289516e-08
1.630052860195072e-06
1.729453234736215e-05
1.735236717842798e-04
1.602794981466470e-03
-2.524883953702671e-02
-3.275962733877914e+01
9.090225583650569e-02
1.442774399597822e+00
9.358201368352827e-02
1.362352192117240e-01
-5.654962847445619e-01
1.259521305204740e-02
2.812217579386780e-03
1.903447241238020e-01
8.220949673523016e-01
3.948371505765507e-04
-1.367678302641368e-03
6.769223033649366e-03
1.112771623240412e-01
6.810743037878951e-01
-8.127962132916176e-07
-9.457347834421762e-05
-7.737808990459289e-04
-5.489425107651075e-04
2.706154550956394e-02
-1.315193694166136e-02
-4.628835366974855e-07
-1.729198205839712e-05
-3.439753872984422e-04
-3.134940878589152e-03
-2.847153409953584e-02
4.286410969697075e-01
7.876600583282652e+01
-1.892466378870984e+00
-2.161596437010127e+00
-4.859957102176131e-01
-1.118002420554065e+00
-3.299993054016922e+00
-2.792585418871490e-02
1.391218159560505e-02
-9.497267920162289e-01
-5.179904586476868e+00
-6.068206588939732e-05
6.950223704207972e-03
1.533577151250691e-02
-2.753789678052529e-01
-3.976676518511942e-01
1.145184553363301e-06
1.112353115929112e-04
2.402716579209772e-03
2.476972434269107e-02
2.026417394119680e-01
-2.932717916326574e+00
-8.514350943139550e+01
8.743339161839597e+00
1.003521500456503e+01
8.143907960904267e-01
3.198379162073299e+00
1.855045310452408e+01
4.996835096975421e-03
-1.541692916238931e-01
1.345728417727253e+00
4.007305115172384e+00
-1.710492918911572e-04
-8.910796609797059e-03
-1.072798986288166e-01
-7.473294091581615e-01
1.040871286641499e+01
2.477681868577220e+01
-1.501731727789784e+01
-2.862443459618763e+01
-1.811976336711438e-01
-2.642915440621548e+00
-1.379536889080945e+01
1.071854489055600e-02
2.766113823288153e-01
1.585283514667677e+00
-2.068256271812679e+01
2.920001419146715e+01
9.184832267645897e+00
9.921553492337058e+00
-2.523234891026405e-01
-2.126314959659203e+00
2.316882193368236e+01
-3.754960763673660e+01
-1.076784756998558e+00
-1.458828530327374e+01
8.318758984573770e+01];
S.PM_s=[
7 0 0
6 1 0
6 0 1
5 2 0
5 1 1
5 0 2
4 3 0
4 2 1
4 1 2
4 0 3
3 4 0
3 3 1
3 2 2
3 1 3
3 0 4
2 5 0
2 4 1
2 3 2
2 2 3
2 1 4
2 0 5
1 5 1
1 4 2
1 3 3
1 2 4
1 1 5
1 0 6
0 5 2
0 4 3
0 3 4
0 2 5
0 1 6
0 0 7
6 0 0
5 1 0
5 0 1
4 2 0
4 1 1
4 0 2
3 3 0
3 2 1
3 1 2
3 0 3
2 4 0
2 3 1
2 2 2
2 1 3
2 0 4
1 5 0
1 4 1
1 3 2
1 2 3
1 1 4
1 0 5
0 5 1
0 4 2
0 3 3
0 2 4
0 1 5
0 0 6
5 0 0
4 1 0
4 0 1
3 2 0
3 1 1
3 0 2
2 3 0
2 2 1
2 1 2
2 0 3
1 4 0
1 3 1
1 2 2
1 1 3
1 0 4
0 5 0
0 4 1
0 3 2
0 2 3
0 1 4
0 0 5
4 0 0
3 1 0
3 0 1
2 2 0
2 1 1
2 0 2
1 3 0
1 2 1
1 1 2
1 0 3
0 4 0
0 3 1
0 2 2
0 1 3
0 0 4
3 0 0
2 1 0
2 0 1
1 2 0
1 1 1
1 0 2
0 3 0
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.p_s=[
4.495221128078107e-01
1.245548825121805e-01
1.853032470845766e-01
-3.452142078488909e-04
-1.772777317812297e-03
2.570840914170591e-03
-2.440152363754440e-04
1.813624724911543e-05
8.891516857183723e-04
5.141462937972462e-03
-1.103726427004444e-05
-8.889203795166548e-06
6.780960743597113e-05
8.768218629287784e-04
-2.825822330414080e-04
6.292964056565373e-07
-6.587993463230426e-07
-1.963694427376431e-06
5.903564224294941e-05
5.986319292277040e-04
-1.657108229468983e-03
-5.563820816281884e-08
-1.196314967446678e-07
4.379755258110860e-07
2.168615919215417e-05
2.174720978400299e-04
-2.787199426905211e-03
-4.741504242777651e-09
4.626598503637192e-08
3.491023775597542e-08
4.020095183563364e-06
7.541809947779025e-05
-7.774437868091079e-04
-3.978070683849476e+00
-8.579246398054966e-01
-1.166041349488243e+00
1.016423784835303e-02
-5.640652953777102e-04
-8.106403505391500e-02
1.807999531816707e-03
5.360123008002256e-05
-1.428198578273746e-02
-3.611488316401906e-02
-2.091464323285701e-05
9.906228665176998e-05
-5.228544968804183e-04
-1.170107470763366e-02
1.497120433323589e-02
-1.800286386547424e-06
6.926820396694892e-06
1.250511617525732e-05
-3.749399949171893e-04
-5.168233248121258e-03
4.526236443539007e-02
6.547765644374506e-08
2.672539013391414e-07
-4.156298328838139e-06
-8.675991993111803e-05
-1.574828392144622e-03
1.745431864146911e-02
1.412645109895821e+01
2.323233076231719e+00
3.005049396608007e+00
-5.012048112153038e-02
5.190023467224773e-02
4.499961637582076e-01
-2.110806146955705e-03
-1.767425578487924e-03
7.841330882551424e-02
5.854522761913817e-02
1.319942578605043e-04
-3.724226798973041e-04
1.401253870140377e-03
4.650447209285513e-02
-2.706703838068150e-01
1.084552457528070e-06
-6.350981341925157e-06
-3.087485241603775e-06
7.364086633988089e-04
1.447368950812665e-02
-1.555777760426586e-01
-2.545196389276697e+01
-3.317036663828704e+00
-4.104031826500355e+00
5.711653607201386e-02
-1.796322900211750e-01
-8.872749035683030e-01
-2.396642217751786e-03
5.350155035682960e-03
-1.811910734553689e-01
6.596787453341868e-01
-9.314894097997961e-05
2.494425290890871e-04
-1.905283993142507e-03
-7.113616686739299e-02
7.001050531655403e-01
2.395449354684673e+01
3.325900081687352e+00
2.744940932958082e+00
3.602212169513015e-02
2.734781020137212e-01
-4.948634605295921e-01
2.766512112499329e-03
-1.913405623318942e-03
1.749610396558740e-01
-1.682898349832526e+00
-1.051637351998415e+01
-2.982742993099567e+00
6.341205258172584e-01
-5.381150940754970e-02
-2.107000396986206e-01
2.314669701756677e+00
1.162661949241928e+00
1.557412527577594e+00
-1.137205176161335e+00
2.795718634969978e-01
];
sig = zeros(size(Freq));
epi_r = sig;
x=[Sugar_Water_Ratio(:),Salt_Water_Ratio(:)];
sig(:) = ND_polyval(S.p_s,[x(:,1)*S.sugar_scale,x(:,2)*S.salt_scale, Freq(:)*S.f_scale],S.PM_s);
epi_r(:) = ND_polyval(S.p_e,[x(:,1)*S.sugar_scale,x(:,2)*S.salt_scale, Freq(:)*S.f_scale],S.PM_e);
endfunction
function [Y,EM,Porder] = ND_polyval(p,X,varargin)
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
function C=SANBA_Surcrose_Solution_Specific_Heat(T, M_Sugar,M_Water,varargin)
% function C=SANBA_Surcrose_Solution_Specific_Heat(T, M_Sugar,M_Water,M_other)
% Compute specific heat for sucrose solution based on
% Bubnik Z. et al. eds. (1995). Sugar Technologists Manual. Bartens, Berlin.
%
% T in C
%
% version 0.1 Qi Duan 2014-11-18
% - creation
% version 0.2 Qi Duan 2015-01-27
% - fix the bug on varargin indexing
if nargin > 3
M_other = varargin{1};
else
M_other = 0;
end
M_all = M_Sugar+M_Water+M_other;
M_dry = (M_Sugar+M_other);
DS = M_dry./M_all*100;
P = M_Sugar./M_dry*100;
C = 4.187-DS.*(0.0297-4.6e-5*P)+7.5e-5.*DS.*T;
endfunction
function XS0 = build_initial_simplex(x0)
%function to generate a simplex based on the initial guess
% x0 needs to be a N-D vector, output is N+1 by N
tmp = (x0(:))';
XS0(1,:) = tmp;
for k = 1:length(tmp)
XS0(k+1,:) = tmp;
if(tmp(k))==0
XS0(k+1,k) = 0.00025;
else
XS0(k+1,k) = XS0(k+1,k)+0.05;
end
end
endfunction
function [ x_opt, n_feval ] = nelder_mead ( x, function_handle, flag )
%*****************************************************************************80
%
%% NELDER_MEAD performs the Nelder-Mead optimization search.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 19 January 2009
%
% Author:
%
% Jeff Borggaard
%
% Reference:
%
% John Nelder, Roger Mead,
% A simplex method for function minimization,
% Computer Journal,
% Volume 7, Number 4, January 1965, pages 308-313.
%
% Parameters:
%
% Input, real X(M+1,M), contains a list of distinct points that serve as
% initial guesses for the solution. If the dimension of the space is M,
% then the matrix must contain exactly M+1 points. For instance,
% for a 2D space, you supply 3 points. Each row of the matrix contains
% one point; for a 2D space, this means that X would be a
% 3x2 matrix.
%
% Input, handle FUNCTION_HANDLE, a quoted expression for the function,
% or the name of an M-file that defines the function, preceded by an
% "@" sign;
%
% Input, logical FLAG, an optional argument; if present, and set to 1,
% it will cause the program to display a graphical image of the contours
% and solution procedure. Note that this option only makes sense for
% problems in 2D, that is, with N=2.
%
% Output, real X_OPT, the optimal value of X found by the algorithm.
%
%
% Define algorithm constants
%
rho = 1; % rho > 0
xi = 2; % xi > max(rho, 1)
gam = 0.5; % 0 < gam < 1
sig = 0.5; % 0 < sig < 1
%tolerance = 1.0E-06;
tolerance = 1.0E-10;
%max_feval = 250;
max_feval = 1000;
%
% Initialization
%
[ temp, n_dim ] = size ( x );
if ( temp ~= n_dim + 1 )
fprintf ( 1, '\n' );
fprintf ( 1, 'NELDER_MEAD - Fatal error!\n' );
error(' Number of points must be = number of design variables + 1\n');
end
if ( nargin == 2 )
flag = 0;
end
if ( flag )
xp = linspace(-5,5,101);
yp = xp;
for i=1:101
for j=1:101
fp(j,i) = feval(function_handle,[xp(i),yp(j)]);
end
end
figure ( 27 )
hold on
contour(xp,yp,fp,linspace(0,200,25))
if ( flag )
plot(x(1:2,1),x(1:2,2),'r')
plot(x(2:3,1),x(2:3,2),'r')
plot(x([1 3],1),x([1 3],2),'r')
pause
plot(x(1:2,1),x(1:2,2),'b')
plot(x(2:3,1),x(2:3,2),'b')
plot(x([1 3],1),x([1 3],2),'b')
end
end
index = 1 : n_dim + 1;
[f ] = evaluate ( x, function_handle );
n_feval = n_dim + 1;
[ f, index ] = sort ( f );
x = x(index,:);
%
% Begin the Nelder Mead iteration.
%
converged = 0;
diverged = 0;
while ( ~converged && ~diverged )
%
% Compute the midpoint of the simplex opposite the worst point.
%
x_bar = sum ( x(1:n_dim,:) ) / n_dim;
%
% Compute the reflection point.
%
x_r = ( 1 + rho ) * x_bar ...
- rho * x(n_dim+1,:);
f_r = feval(function_handle,x_r);
n_feval = n_feval + 1;
%
% Accept the point:
%
if ( f(1) <= f_r && f_r <= f(n_dim) )
x(n_dim+1,:) = x_r;
f(n_dim+1 ) = f_r;
if (flag)
title('reflection')
end
%
% Test for possible expansion.
%
elseif ( f_r < f(1) )
x_e = ( 1 + rho * xi ) * x_bar ...
- rho * xi * x(n_dim+1,:);
f_e = feval(function_handle,x_e);
n_feval = n_feval+1;
%
% Can we accept the expanded point?
%
if ( f_e < f_r )
x(n_dim+1,:) = x_e;
f(n_dim+1 ) = f_e;
if (flag), title('expansion'), end
else
x(n_dim+1,:) = x_r;
f(n_dim+1 ) = f_r;
if (flag), title('eventual reflection'), end
end
%
% Outside contraction.
%
elseif ( f(n_dim) <= f_r && f_r < f(n_dim+1) )
x_c = (1+rho*gam)*x_bar - rho*gam*x(n_dim+1,:);
f_c = feval(function_handle,x_c); n_feval = n_feval+1;
if (f_c <= f_r) % accept the contracted point
x(n_dim+1,:) = x_c;
f(n_dim+1 ) = f_c;
if (flag), title('outside contraction'), end
else
[x,f] = shrink(x,function_handle,sig); n_feval = n_feval+n_dim;
if (flag), title('shrink'), end
end
%
% F_R must be >= F(N_DIM+1).
% Try an inside contraction.
%
else
x_c = ( 1 - gam ) * x_bar ...
+ gam * x(n_dim+1,:);
f_c = feval(function_handle,x_c);
n_feval = n_feval+1;
%
% Can we accept the contracted point?
%
if (f_c < f(n_dim+1))
x(n_dim+1,:) = x_c;
f(n_dim+1 ) = f_c;
if (flag), title('inside contraction'), end
else
[x,f] = shrink(x,function_handle,sig); n_feval = n_feval+n_dim;
if (flag), title('shrink'), end
end
end
%
% Resort the points. Note that we are not implementing the usual
% Nelder-Mead tie-breaking rules (when f(1) = f(2) or f(n_dim) =
% f(n_dim+1)...
%
[ f, index ] = sort ( f );
x = x(index,:);
%
% Test for convergence
%
converged = f(n_dim+1)-f(1) < tolerance;
%
% Test for divergence
%
diverged = ( max_feval < n_feval );
if (flag)
plot(x(1:2,1),x(1:2,2),'r')
plot(x(2:3,1),x(2:3,2),'r')
plot(x([1 3],1),x([1 3],2),'r')
pause
plot(x(1:2,1),x(1:2,2),'b')
plot(x(2:3,1),x(2:3,2),'b')
plot(x([1 3],1),x([1 3],2),'b')
end
end
if ( 0 )
fprintf('The best point x^* was: %d %d\n',x(1,:));
fprintf('f(x^*) = %d\n',f(1));
end
x_opt = x(1,:);
if ( diverged )
fprintf ( 1, '\n' );
fprintf ( 1, 'NELDER_MEAD - Warning!\n' );
fprintf ( 1, ' The maximum number of function evaluations was exceeded\n')
fprintf ( 1, ' without convergence being achieved.\n' );
end
return
endfunction
function f = evaluate ( x, function_handle )
%*****************************************************************************80
%
%% EVALUATE handles the evaluation of the function at each point.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 19 January 2009
%
% Author:
%
% Jeff Borggaard
%
% Reference:
%
% John Nelder, Roger Mead,
% A simplex method for function minimization,
% Computer Journal,
% Volume 7, Number 4, January 1965, pages 308-313.
%
% Parameters:
%
% Input, real X(N_DIM+1,N_DIM), the points.
%
% Input, real FUNCTION_HANDLE ( X ), the handle of a MATLAB procedure
% to evaluate the function.
%
% Output, real F(1,NDIM+1), the value of the function at each point.
%
[ temp, n_dim ] = size ( x );
f = zeros ( 1, n_dim+1 );
for i = 1 : n_dim + 1
f(i) = feval(function_handle,x(i,:));
end
return
endfunction
function [ x, f ] = shrink ( x, function_handle, sig )
%*****************************************************************************80
%
%% SHRINK shrinks the simplex towards the best point.
%
% Discussion:
%
% In the worst case, we need to shrink the simplex along each edge towards
% the current "best" point. This is quite expensive, requiring n_dim new
% function evaluations.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 19 January 2009
%
% Author:
%
% Jeff Borggaard
%
% Reference:
%
% John Nelder, Roger Mead,
% A simplex method for function minimization,
% Computer Journal,
% Volume 7, Number 4, January 1965, pages 308-313.
%
% Parameters:
%
% Input, real X(N_DIM+1,N_DIM), the points.
%
% Input, real FUNCTION_HANDLE ( X ), the handle of a MATLAB procedure
% to evaluate the function.
%
% Input, real SIG, ?
%
% Output, real X(N_DIM+1,N_DIM), the points after shrinking was applied.
%
% Output, real F(1,NDIM+1), the value of the function at each point.
%
[temp,n_dim] = size(x);
x1 = x(1,:);
f(1) = feval(function_handle,x1);
for i=2:n_dim+1
x(i,:) = sig*x(i,:) + (1-sig)*x(1,:);
f(i) = feval(function_handle,x(i,:));
end
return
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," Agar [%% w/v]\n");
fprintf(1," Benzoic acid [%% w/v]\n");
fprintf(1," Temperature [degC]\n");
fprintf(1,"\n");
sig=input("Conductivity [S/m] : ")
epi_r=input("Permittivity : ")
f=input("Frequency [MHz] : ")
M_water=input("Volume [ml] : ")
Agar_Percent=input("Agar [% w/v] : ")
BA_Percent=input("Benzoic acit [% w/v]: ")
T=input("Temperature [degC] : ")
M_water0 = 2300;
y0 = (40*sig/0.77)/M_water0;
x0 = (80-epi_r)/0.014/M_water0;
% 2025-04-21: Input conversion from MHz and %
f=f*1.e6;
Agar_Percent=Agar_Percent/100.;
BA_Percent=BA_Percent/100.;
[M_NaCl,M_Sugar, M_Agar, M_BA, M_Water] = SANBA_fitting(sig,epi_r, f, M_water, Agar_Percent, BA_Percent, [x0,y0]);
Vol_est = zeros(size(M_BA));
Vol_est(:) = (M_Sugar(:)+M_Water(:))./SANBA_Surcrose_Solution_Density(M_Sugar(:)./(M_Sugar(:)+M_Water(:)))*1e3;%convert kg/m^3 to g/ml
[sig_exp,epi_exp] = SANBA_get_property(M_Sugar./M_Water, M_NaCl/M_Water, f);
rou_est = SANBA_Surcrose_Solution_Density(M_Sugar(:)./(M_Sugar(:)+M_Water(:)));
c_est = SANBA_Surcrose_Solution_Specific_Heat(T,M_Sugar(:),M_Water(:),(M_NaCl(:)+M_Agar(:)+M_BA(:)));
k_est = SANBA_Surcrose_Solution_Thermal_Conductivity(T,M_Sugar(:)./(M_Sugar(:)+M_Water(:)+M_NaCl(:)+M_Agar(:)+M_BA(:))*100);
err_flag = (sig_exp<0)|(epi_exp<0)|(epi_exp>80);
err_flag = err_flag | (M_NaCl<0)|(M_Sugar<0);
err_flag = err_flag | ((abs(epi_exp-epi_r)./epi_r)>0.1);
if sig == 0
err_flag = err_flag | abs(sig_exp)>1e-3;
else
err_flag = err_flag | ((abs(sig_exp-sig)./sig)>0.1);
end
err_flag = err_flag | (rou_est < 0) | (c_est< 0) | (k_est< 0) ;
if (err_flag)
fprintf(1,"\n");
fprintf("ERROR: Unable to reliably determine this recipe!\n");
fprintf(" Please adjust one or more of the parameters\n");
fprintf(1,"\n");
else
fprintf(1,"\n");
fprintf(1,"REQUIRED INGREDIENTS:\n")
fprintf(1,"NaCl: %f g\n",M_NaCl);
fprintf(1,"Sugar: %f g\n",M_Sugar);
fprintf(1,"Agar: %f g\n",M_Agar);
fprintf(1,"Benzoic acid: %f g\n",M_BA);
fprintf(1,"Water: %f ml\n",M_Water);
fprintf(1,"\n");
fprintf(1,"Est. volume: %f ml\n",Vol_est);
fprintf(1,"Est. specific heat: %f J/(g*K)\n",c_est);
fprintf(1,"Est. thermal conductivity: %f W/(m*K)\n",k_est);
fprintf(1,"\n");
end
To embed this project on your website, copy the following code and paste it into your website's HTML: