Hi Jort,
I understand that you are having a problem with updating the values in your code through a loop.
It seems that the main logic for the centrifugation process is contained in two files: “odepaper.m”, which defines the ordinary differential equation (ODE) system, and “centrifugation.m”, which sets up the parameters and solves the ODE using ode45.
I noticed that the loop you are referring to in “odepaper.m” attempts to modify the concentration c(i) by giving back excess to the previous cells. However, this modification is done after the calculation of the derivatives “dcdt” and does not affect the output of the function. Consequently, the ode45 solver, which calls “odepaper”, does not see the changes made in the loop.
To address this, you may need to modify the function to take into account the excess concentration during the calculation of the derivatives dcdt.
The key change is to move the loop that handles excess concentrations before the calculation of the derivatives. This will ensure that the excess concentrations are properly redistributed among the cells before the derivatives are computed.
Here’s the modified “odepaper.m” file,
function [dcdt] = odepaper(t, c, Nr, dr, vgs, cmax, cini, n)
dcdt = zeros(Nr,1);
r = zeros(Nr,1);
v = zeros(Nr,1);
excess = 0;
%% Mass replacement
for i = Nr:-1:2 % Force system to give back excess to previous cells
c(i) = c(i) + excess;
if c(i) > cmax
excess = c(i) - cmax;
c(i) = cmax;
c(i-1) = c(i-1) + excess;
else
excess = 0;
end
end
%% Loop
for i = 2:Nr
r(i) = r(i-1) + dr;
v(i) = (1 - c(i))^n * vgs;
vin = (1 - cini)^n * vgs;
rin = 0.5 * dr;
if i == 1 % Initial conditions
dcdt(1) = (1/dr) * ((rin * cini * vin) / r(i) - c(i) * v(i));
elseif i == Nr % Wall conditions
dcdt(Nr) = (1/dr) * ((r(i-1) * c(i-1) * v(i-1)) / r(i)) + excess;
else % Everything in between
dcdt(i) = (1/dr) * ((r(i-1) * c(i-1) * v(i-1)) / r(i) - c(i) * v(i)) + excess;
end
end
end
I moved the section handling excess concentrations to the beginning of the function, so that it modifies the concentrations "c" before the derivatives "dcdt" are computed.
