k = 8.34e-6 * 3600; A=105.6; clim = 32.383; V0=100;
c0=0.15; cs0 = 120;
cf = 11.913; csf = 0.1;

% Economical weights
wt0 = 3.6; wd0 = 10;

res = [];

% Minimum time
wt = 1; wd = 0;
cast = clim/exp(1);
lg = log(clim/cast); lg = 1;
dt1 = c0*V0/(k*A*clim)* (mfun('Li',clim/c0) - mfun('Li',clim/cast));
dt2 = c0*V0/(k*A*cast)*( log(cs0*cf/(cast*csf)) / lg );
Vw2 = k*A*lg * dt2;
Vw3 = c0*V0 * (1/cf - 1/cast);
time = dt1+dt2;
Volume = Vw2+Vw3;
res=[res; wt wd time Volume (wt0*time+wd0*Volume)*100/15];

% Almost minimum time
wt = 1; wd = 1e-1;
cast = clim * exp((wt - sqrt(wt^2+4*k*A*wt*wd))/(2*k*A*wd));
lg = log(clim/cast);

dt1 = c0*V0/(k*A*clim)* (mfun('Li',clim/c0) - mfun('Li',clim/cast));
dt2 = c0*V0/(k*A*cast)*( log(cs0*cf/(cast*csf)) / lg );
Vw2 = k*A*lg * dt2;
Vw3 = c0*V0 * (1/cf - 1/cast);
time = dt1+dt2;
Volume = Vw2+Vw3;
res=[res; wt wd time Volume (wt0*time+wd0*Volume)*100/15];

% Economical
wt = 3.6; wd = 10;
cast = clim * exp((wt - sqrt(wt^2+4*k*A*wt*wd))/(2*k*A*wd));
lg = log(clim/cast);

dt1 = c0*V0/(k*A*clim)* (mfun('Li',clim/c0) - mfun('Li',clim/cast));
dt2 = c0*V0/(k*A*cast)*( log(cs0*cf/(cast*csf)) / lg );
Vw2 = k*A*lg * dt2;
Vw3 = c0*V0 * (1/cf - 1/cast);
time = dt1+dt2;
Volume = Vw2+Vw3;
res=[res; wt wd time Volume (wt0*time+wd0*Volume)*100/15];

% Almost minimum diluant
wt = 0.01; wd = 1;
cast = clim * exp((wt - sqrt(wt^2+4*k*A*wt*wd))/(2*k*A*wd));
lg = log(clim/cast);

dt1 = c0*V0/(k*A*clim)* (mfun('Li',clim/c0) - mfun('Li',clim/cast));
dt2 = c0*V0/(k*A*cast)*( log(cs0*cf/(cast*csf)) / lg );
Vw2 = k*A*lg * dt2;
Vw3 = c0*V0 * (1/cf - 1/cast);
time = dt1+dt2;
Volume = Vw2+Vw3;
res=[res; wt wd time Volume (wt0*time+wd0*Volume)*100/15];

% Minimum diluant
wt = 0; wd = 1;
Vw2 = c0*V0/clim * log(cs0*cf/(clim*csf));
Vw3 = c0*V0 * (1/cf - 1/clim);
time = Inf;
Volume = Vw2+Vw3;
table=[res; wt wd time Volume (wt0*time+wd0*Volume)*100/15]
