I want to know if my code is well written based on a set of criteria.

Question

Daniel 2024-3-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2095466-i-want-to-know-if-my-code-is-well-written-based-on-a-set-of-criteria

评论： Daniel 2024-3-18

To estimate DIP_s, it is necessary to approximate the expected DIP concentration within the current time step used in the integration(i.e DIP_interm . If the time step is one day, this concentration (DIP_interm) is the sum of the present DIP and the result of Equation 1 with DIP_fert= 0. The difference between DIP_interm and DIP_O is an estimate of DIP_s.

dDIPdt = DIP_Fert + DIP_resp_phy - DIP_grow_phy - (0.1 * DIP) ...............................................................EQN(1)

DIP_req = DIP_grow_phy-(DIP_resp_phy + DIP_msc);

DIP_fert can then be calculated by applying the following rules:

if DIP_interm < DIP_O

DIP_Fert = DIP_req + (DIP_O - DIP_interm);

elseif 0 <= DIP_s <= DIP_req

DIP_Fert = DIP_req - DIP_s;

elseif DIP_s > DIP_req

DIP_Fert = 0;

end

Finally, Equation 1 is re-evaluated by inserting the value of DIP_fert obtained based on the above rules so as to obtain the rate change of DIP. The sum of DIP_fert calculated at each time step over the entire integration period is a rough estimate of the overall fertilizer N requirements for that period

----------------------------------------------------------------------------------------------------------------------------------------------------------------------------

% Define simulation period

startdate = datetime(2023, 01, 01, 'Format', 'uuuu-MM-dd');

finishdate = datetime(2023, 06, 30, 'Format', 'uuuu-MM-dd');

tspans = (1:1:days(finishdate - startdate) + 1);

%Dissolved inorganic phosphorous (DIP)

T_w = 20;

T_min = 18;

T_max = 27;

T_opt = 33;

if T_w < T_opt

Tau = exp(-4.6 * ((T_opt - T_w) / (T_opt - T_min)).^4);

else

Tau = exp(-4.6 * ((T_w - T_opt) / (T_max - T_opt)).^4);

end

V_n = 1;

GPP_lambda = 7.0; % assumed as & in g m^-3 d^-1

GPP = GPP_lambda * Tau * V_n;

DIP_IN = 0.2;

DIP_Fert = 0;

Redfield_ratio_p = 0.025;

DIP_grow_phy = GPP*Redfield_ratio_p;

DIP_resp_phy = 0.5 *DIP_grow_phy;

% Initialize DIP

DIP = zeros(length(tspans), 1); % Initialize DIP vector

% Initialize arrays to store DIP_Fert

DIP_Fert_values = zeros(length(tspans), 1);

for t = 1:length(tspans)-1

dDIPdt = @(t, DIP) calculate_dDIPdt(t, DIP, DIP_Fert, DIP_resp_phy, DIP_grow_phy);

[t_integrated, DIP_integrated] = ode45(@(t, DIP) dDIPdt(t, DIP), tspans(t:t+1), DIP_IN);

DIP(t+1) = DIP_integrated(end);

% Recalculate DIP_Fert based on conditions

DIP_O = 0.1;

DIP_msc = 0.1*DIP(t);

DIP_req = DIP_grow_phy-(DIP_resp_phy + DIP_msc);

DIP_interm = DIP(t) + (DIP_resp_phy + DIP_msc - DIP_grow_phy);

DIP_s =max(0,(DIP_interm - DIP_O));

if DIP_interm < DIP_O

DIP_Fert = DIP_req + (DIP_O - DIP_interm);

elseif all((0 <= DIP_s) & (DIP_s <= DIP_req))

DIP_Fert = DIP_req - DIP_s;

elseif DIP_s > DIP_req

DIP_Fert = 0;

end

% Check if DIP_Fert is negative and set it to zero

if DIP_Fert < 0

DIP_Fert = 0;

end

% Update DIP_Fert array

DIP_Fert_values(t+1) = DIP_Fert;

% Update initial condition for the next iteration

DIP_IN = DIP(t+1);

end

% Output DIP_Fert values at different time steps

disp('DIP_Fert values:');

DIP_Fert values:

disp(DIP_Fert_values);

0 0.1131 0 0 0 0 0 0 0 0 0 0.0163 0.0330 0.0294 0.0072 0 0.0020 0.0200 0.0341 0.0262 0.0030 0 0.0066 0.0242 0.0326 0.0200 0 0 0.0133 0.0303 0.0304 0.0112 0 0 0.0179 0.0344 0.0289 0.0050 0 0.0032 0.0211 0.0337 0.0246 0.0020 0 0.0087 0.0261 0.0319 0.0173 0 0 0.0147 0.0315 0.0299 0.0093 0 0.0008 0.0189 0.0345 0.0278 0.0040 0 0.0047 0.0224 0.0332 0.0226 0.0008 0 0.0111 0.0282 0.0311 0.0141 0 0 0.0164 0.0330 0.0294 0.0071 0 0.0020 0.0201 0.0340 0.0261 0.0030 0 0.0067 0.0243 0.0325 0.0199 0 0 0.0134 0.0303 0.0304 0.0110 0 0 0.0180 0.0345 0.0289 0.0050 0 0.0033 0.0212 0.0336 0.0245 0.0020 0 0.0087 0.0261 0.0319 0.0172 0 0 0.0148 0.0316 0.0299 0.0092 0 0.0008 0.0189 0.0344 0.0277 0.0040 0 0.0047 0.0225 0.0332 0.0225 0.0007 0 0.0112 0.0283 0.0311 0.0140 0 0 0.0165 0.0331 0.0294 0.0070 0 0.0021 0.0201 0.0340 0.0260 0.0029 0 0.0068 0.0244 0.0325 0.0197 0 0 0.0135 0.0304 0.0303 0.0110 0 0 0.0180 0.0345 0.0289 0.0049 0 0.0033 0.0212 0.0336 0.0245 0.0019 0 0.0088 0.0262 0.0319 0.0171 0 0 0.0148 0.0316 0.0299 0.0092

% Plot DIP

plot(DIP);

xlabel('Time');

ylabel('Concentration');

legend('DIP',Location='northwest');

function dDIPdt = calculate_dDIPdt(~, DIP, DIP_Fert, DIP_resp_phy, DIP_grow_phy)

dDIPdt = DIP_Fert + DIP_resp_phy - DIP_grow_phy +- (0.1 * DIP);

end