Agronomist: MQL recipes

Recipes

Growing degree days and frost risk for maize
Economic optimum nitrogen rate
July rainfall gamma law and irrigation need
Reference evapotranspiration with FAO-56 Penman-Monteith
Variety trial: ANOVA and least significant difference
Phosphorus response: Mitscherlich law and economic optimum
Field slope, aspect and USLE LS factor from a DEM
Seed lot germination: MGT, rate and T50
Pesticide dissipation: SFO and FOMC DT50 / DT90
NPK fertilizer blend with an LU decomposition

executed on MentDB

Growing degree days and frost risk for maize

From two weeks of daily max and min temperatures of a US station (degF), compute the daily means in degC with vector operations, the growing degree days above the 50 degF base, the day when 100 GDD are reached (interpolation of the cumulative curve) and the probability of frost tonight from a forecast minimum of 2.5 degC with a standard deviation of 1.8 degC.

Pipeline

01 · convert(Tmax + Tmin)/2 - 32, x 5/9vector add · vector scale · vector sub · matrix ones→02 · base50 degF in degCunit convert→03 · GDDCumulated max(0, Tmean - Tbase)math calc→04 · dateDay of 100 GDDnum interpolate→05 · frostP(T < 0) with erferf(x) · math eval

Script

admin

#1. Two weeks of daily max and min air temperature from a US station, in degF;
-> "[tx]" "[68, 72, 75, 79, 81, 77, 73, 70, 74, 80, 84, 86, 82, 78]";
-> "[tn]" "[46, 50, 52, 55, 58, 56, 51, 48, 50, 56, 60, 63, 61, 57]";
#2. Daily mean in degC: Tmean - 32 times 5/9, with vector operations;
-> "[tm]" (vector scale (vector sub (vector scale (vector add [tx] [tn]) 0.5) (vector scale (matrix row (matrix ones 1 14) 1) 32)) (bigdec div 5 9 12));
#3. Growing degree days above the base 50 degF for maize, cumulated;
-> "[tb]" (unit convert 50 "degF" "degC");
-> "[cum]" "[0"; -> "[g]" 0; -> "[days]" "[0";
for (-> "[i]" 0) (< [i] 14) (++ "[i]") {
  -> "[g]" (math calc (concat "[g]+max(0, " (json load select "rx_tm" [tm] (concat "/[" [i] "]")) "-[tb])"));
  -> "[cum]" (concat [cum] "," [g]); -> "[days]" (concat [days] "," (+ [i] 1));
};
-> "[cum]" (concat [cum] "]"); -> "[days]" (concat [days] "]");
#4. Day when 100 GDD are reached and frost risk tonight: Tmin forecast 2.5 degC, sd 1.8;
-> "[d100]" (num interpolate [cum] [days] 100 "pchip");
-> "[frost]" (math eval "0.5*(1+erf((0-2.5)/(1.8*sqrt(2))))");
concat "base " (math round [tb] 1) " degC, mean temperatures " (math round (stat mean [tm]) 2) " degC, GDD after 14 days = " (math round [g] 1) ", 100 GDD reached on day " (math round [d100] 1) ", frost probability tonight " (math round (math calc "100*[frost]") 1) " %"

Result

mentdb

base 10 degC, mean temperatures 18.77 degC, GDD after 14 days = 122.8, 100 GDD reached on day 11.9, frost probability tonight 8.2 %

Go further

Use the single-sine method: integrate the daily sine curve above the base with num integrate
Cap Tmax at 86 degF as in the modified GDD method
Draw the cumulative curve with plot interpolate

executed on MentDB

Economic optimum nitrogen rate

A nitrogen rate trial gives the grain yield for 7 rates. Fit the quadratic response, get the rate of maximum yield and the economic optimum rate (EONR) where the marginal yield equals the price ratio of nitrogen to grain, then subtract the soil nitrate measured around the plot, interpolated from a 3 x 3 grid of soil samples.

Pipeline

01 · fitY = a + b N + c N^2stat fit_model→02 · derivedY/dNnum poly_derive→03 · optimumRoots of dY/dN and dY/dN - rationum poly_roots · num poly_add→04 · yieldYield at the EONRnum poly_eval→05 · soilBilinear interpolation of soil nitratenum interpolate2

Script

admin

#1. Nitrogen rate trial: N in kg/ha, grain yield in t/ha;
-> "[N]" "[0, 40, 80, 120, 160, 200, 240]";
-> "[Y]" "[5.1, 6.6, 7.8, 8.6, 9.0, 9.1, 8.9]";
#2. Quadratic response Y = a + b N + c N^2;
json load "rx_nf" (stat fit_model "a+b*x+c*x^2" "x" [N] [Y] "{\"a\": 5, \"b\": 0.04, \"c\": -0.0001}");
-> "[p]" (concat "[" (json select "rx_nf" /parameters/a) ", " (json select "rx_nf" /parameters/b) ", " (json select "rx_nf" /parameters/c) "]");
#3. Economic optimum: dY/dN = price of N / price of grain = 1.10 / 200 t per kg;
-> "[d]" (num poly_derive [p]);
-> "[nmax]" (json load select "rx_r1" (num poly_roots [d]) "/[0]");
-> "[eonr]" (json load select "rx_r2" (num poly_roots (num poly_add [d] "[-0.0055]")) "/[0]");
#4. Soil nitrate test of the field at the plot position, from a 3 x 3 grid of samples;
-> "[soil]" (num interpolate2 "[0, 50, 100]" "[0, 40, 80]" "[[28, 34, 31], [36, 42, 38], [30, 35, 33]]" 62 25 "linear");
concat "max yield N = " (math round [nmax] 0) " kg/ha, EONR = " (math round [eonr] 0) " kg/ha, yield there " (math round (num poly_eval [p] [eonr]) 2) " t/ha, soil nitrate " (math round [soil] 1) " kg/ha, fertilizer to apply " (math round (- [eonr] [soil]) 0) " kg N/ha, r2 = " (math round (json select "rx_nf" /r2) 4)

Result

mentdb

max yield N = 191 kg/ha, EONR = 166 kg/ha, yield there 9.07 t/ha, soil nitrate 38.2 kg/ha, fertilizer to apply 128 kg N/ha, r2 = 0.9997

Go further

Fit a quadratic-plateau model with stat fit_model and min
Draw the trial and the curve with plot fit
Test the sensitivity of the EONR to the grain price with a for loop

executed on MentDB

July rainfall gamma law and irrigation need

Twenty years of July rainfall are fitted with a gamma law by maximum likelihood (ln k - digamma(k) = ln mean - mean ln x, standard error with trigamma). The script gives the probability that the rain does not cover the crop need (ETc minus soil reserve), the 80 % dependable rainfall, and the irrigation volume and pumping time for 12 ha.

Pipeline

01 · meanExact meanbigdec mean→02 · MLEShape k from digamma, error from trigammanum solve · digamma(x) · polygamma(n, x)→03 · probabilityP(rain < need) = P(k, need/theta)gamma_inc(a, x)→04 · dependable20 % quantile of the gamma lawstat dist_inv→05 · volumemm x ha in m3unit calc

Script

admin

#1. July rainfall over 20 years at the farm, in mm;
-> "[r]" "[62, 38, 95, 41, 120, 57, 33, 78, 49, 66, 104, 29, 71, 54, 88, 45, 36, 83, 60, 52]";
-> "[n]" 20; -> "[m]" (bigdec mean [r] 10);
#2. Gamma law by maximum likelihood: ln k - digamma k = ln mean - mean of ln x;
-> "[ln]" "["; -> "[sep]" "";
json load "rx_rain" [r];
json parse_array "rx_rain" "/" "[x]" { -> "[ln]" (concat [ln] [sep] (math log [x])); -> "[sep]" ","; };
-> "[s]" (math calc (concat "ln([m])-" (stat mean (concat [ln] "]"))));
-> "[k]" (num solve (concat "ln(x)-digamma(x)-" [s]) "x" 0.5 50);
-> "[th]" (math calc "[m]/[k]");
-> "[sek]" (math calc "sqrt([k]/([n]*([k]*polygamma(1, [k])-1)))");
#3. Crop need: ETc 4.6 mm/day over 31 days minus 55 mm of soil reserve;
-> "[need]" (math calc "4.6*31-55");
-> "[pdef]" (math calc "gamma_inc([k], [need]/[th])");
#4. Dependable rain at 80 %: 20 % quantile, irrigation volume for 12 ha;
-> "[p20]" (stat dist_inv "gamma" (concat "[" [k] ", " [th] "]") 0.2);
-> "[vol]" (unit calc (concat (math calc "[need]-[p20]") " mm * 12 ha") "m^3");
concat "gamma law k = " (math round [k] 2) " +/- " (math round [sek] 2) ", theta = " (math round [th] 1) " mm; P(July rain < " (math round [need] 1) " mm) = " (math round (math calc "100*[pdef]") 1) " %; 80 % dependable rain " (math round [p20] 1) " mm, irrigation " (math round [vol] 0) " m3 = " (math round (/ [vol] 45) 1) " h of pumping at 45 m3/h"

Result

mentdb

gamma law k = 6.89 +/- 2.13, theta = 9.1 mm; P(July rain < 87.6 mm) = 85 %; 80 % dependable rain 42.5 mm, irrigation 5414 m3 = 120.3 h of pumping at 45 m3/h

Go further

Check the fit with stat ks_test against the gamma law
Draw the rainfall with plot histogram
Repeat for each month and sum the yearly irrigation need

executed on MentDB

Reference evapotranspiration with FAO-56 Penman-Monteith

From one day of weather data, compute the wind at 2 m, the psychrometric constant at the site altitude, the vapour pressures and the slope of the saturation curve (numerical derivative), the reference ET0 and the crop water need of maize at mid-season.

Pipeline

01 · windu2 from 10 m, log profileunit convert→02 · vapoures, ea from e0(T)math eval→03 · slopeDelta = de0/dTnum derivative→04 · ET0Penman-Monteith, then ETc = Kc.ET0math calc

Script

admin

#1. Weather of one July day in Reunion lowlands: Tmax 30, Tmin 18 degC, RH 65 %, wind 9 km/h at 10 m, Rn 15 MJ/m2/d, altitude 100 m;
-> "[Tx]" 30; -> "[Tn]" 18; -> "[RH]" 65; -> "[Rn]" 15; -> "[z]" 100;
#2. Wind at 2 m from the log profile, psychrometric constant;
-> "[u2]" (math calc (concat (unit convert 9 "km/h" "m/s") "*4.87/ln(67.8*10-5.42)"));
-> "[g]" (math calc "0.000665*101.3*((293-0.0065*[z])/293)^5.26");
#3. Saturation vapour pressure e0 of T = 0.6108.exp of 17.27T over T+237.3, slope Delta, actual vapour pressure;
-> "[e0]" "0.6108*exp(17.27*T/(T+237.3))";
-> "[es]" (/ (+ (math eval [e0] (concat "{\"T\": " [Tx] "}")) (math eval [e0] (concat "{\"T\": " [Tn] "}"))) 2);
-> "[Tm]" (/ (+ [Tx] [Tn]) 2);
-> "[D]" (num derivative [e0] "T" [Tm]);
-> "[ea]" (math calc "[RH]/100*[es]");
#4. FAO-56 Penman-Monteith ET0 with G = 0, then maize at mid-season Kc = 1.20;
-> "[ET0]" (math calc "(0.408*[D]*[Rn]+[g]*900/([Tm]+273)*[u2]*([es]-[ea]))/([D]+[g]*(1+0.34*[u2]))");
concat "u2 = " (math round [u2] 2) " m/s, gamma = " (math round [g] 5) " kPa/degC, es = " (math round [es] 3) ", ea = " (math round [ea] 3) " kPa, Delta = " (math round [D] 4) " kPa/degC; ET0 = " (math round [ET0] 2) " mm/d, ETc maize = " (math round (* 1.2 [ET0]) 2) " mm/d, 7-day net irrigation = " (math round (* 7 (* 1.2 [ET0])) 1) " mm"

Result

mentdb

u2 = 1.87 m/s, gamma = 0.06658 kPa/degC, es = 3.154, ea = 2.05 kPa, Delta = 0.1791 kPa/degC; ET0 = 5.25 mm/d, ETc maize = 6.3 mm/d, 7-day net irrigation = 44.1 mm

Go further

Estimate Rn from the sunshine hours and the extraterrestrial radiation (FAO-56 equations 21 to 40)
Compare with the Hargreaves equation when only temperatures are available

executed on MentDB

Variety trial: ANOVA and least significant difference

Four wheat varieties are compared on four plots each. The one-way ANOVA tests the variety effect, the Fisher LSD at 5 % says which means differ, the CV checks the quality of the trial, and Kruskal-Wallis gives a non-parametric check.

Pipeline

01 · ANOVAF, p, residual mean squarestat anova→02 · LSDt(0.975, df).sqrt(2.MSE/r)stat dist_inv→03 · qualityMeans and CVstat mean→04 · checkRank-based teststat kruskal

Script

admin

#1. Wheat variety trial: 4 varieties, 4 plots each, grain yield in t/ha;
-> "[g]" "[[7.2, 6.8, 7.5, 7.0], [8.1, 7.9, 8.4, 8.0], [6.5, 6.9, 6.2, 6.6], [7.8, 7.4, 8.0, 7.7]]";
json load "rx_tu1_av" (stat anova [g]);
#2. Fisher least significant difference at 5 % from the residual mean square;
-> "[ms]" (json select "rx_tu1_av" "/ms_within"); -> "[df]" (json select "rx_tu1_av" "/df_within");
-> "[t]" (stat dist_inv "student" (concat "[" [df] "]") 0.975);
-> "[lsd]" (math calc "[t]*sqrt(2*[ms]/4)");
#3. Variety means, coefficient of variation of the trial, non-parametric check;
-> "[m]" (concat (stat mean "[7.2, 6.8, 7.5, 7.0]") ", " (stat mean "[8.1, 7.9, 8.4, 8.0]") ", " (stat mean "[6.5, 6.9, 6.2, 6.6]") ", " (stat mean "[7.8, 7.4, 8.0, 7.7]"));
-> "[cv]" (math calc (concat "100*sqrt([ms])/" (stat mean (concat "[" [m] "]"))));
concat "F = " (math round (json select "rx_tu1_av" "/f") 2) ", p = " (math round_sig (json select "rx_tu1_av" "/p_value") 3) ", eta2 = " (math round (json select "rx_tu1_av" "/eta2") 3) "; means " [m] " t/ha; LSD 5 % = " (math round [lsd] 3) " t/ha, CV = " (math round [cv] 1) " %; Kruskal-Wallis p = " (math round_sig (json load select "rx_tu1_kw" (stat kruskal [g]) "/p_value") 3)

Result

mentdb

F = 26.34, p = 0.0000145, eta2 = 0.868; means 7.125, 8.1, 6.55, 7.725 t/ha; LSD 5 % = 0.409 t/ha, CV = 3.6 %; Kruskal-Wallis p = 0.00477

Go further

Add the block effect of a randomized complete block design (two-way ANOVA)
Correct the multiple comparisons with Tukey HSD instead of the LSD

executed on MentDB

Phosphorus response: Mitscherlich law and economic optimum

A wheat phosphorus trial is fitted with the Mitscherlich law y = A(1 - exp(-c(x + b))), where b is the phosphorus supplied by the soil. The script gives the relative yield of the unfertilized plot, the economic optimum rate for the price ratio, and plots the trial with the results in the title.

Pipeline

01 · fitMitscherlich, 3 parametersstat fit_model→02 · parametersA, c, bjson load select→03 · RY0-expm1(-c.b)math expm1→04 · optimumln(A.c/ratio)/c - bmath log→05 · plotTrial scatterpa xy_scatter

Script

admin

#1. Phosphorus trial on wheat: P rate in kg P2O5/ha and grain yield in t/ha;
-> "[x]" "[0, 30, 60, 90, 120, 150]";
-> "[y]" "[4.1, 5.6, 6.5, 7.0, 7.2, 7.3]";
#2. Mitscherlich law y = A.1 - exp -c.x + b, b = P supplied by the soil;
-> "[fit]" (stat fit_model "A*(1-exp(-c*(x+b)))" "x" [x] [y] "{\"A\": 7.5, \"c\": 0.02, \"b\": 30}");
-> "[a]" (json load select "rx_p4d_mit" [fit] "/parameters/A");
-> "[c]" (json select "rx_p4d_mit" "/parameters/c");
-> "[b]" (json select "rx_p4d_mit" "/parameters/b");
#3. Relative yield of the unfertilized plot -expm1 -c.b, and economic optimum when dy/dx = price of P / price of grain;
-> "[ry0]" (- 0 (math expm1 (* (- 0 [c]) [b])));
-> "[ratio]" (/ 1.1 220);
-> "[xopt]" (- (/ (math log (/ (* [a] [c]) [ratio])) [c]) [b]);
#4. Scatter of the trial, the fitted law in the title;
pa xy_scatter "[[0, 4.1], [30, 5.6], [60, 6.5], [90, 7.0], [120, 7.2], [150, 7.3]]" (concat "Yield t/ha vs P2O5 kg/ha: A = " (math round [a] 2) ", c = " (math round [c] 4) ", soil b = " (math round [b] 1) ", RY0 = " (math round (* 100 [ry0]) 1) " %, optimum = " (math round [xopt] 0) " kg/ha")

Result

mentdb

The chart is shown in the Plot tab of the editor.

Go further

Compare with the quadratic-plateau model of the nitrogen recipes (stat fit_model).
Recompute the optimum for several grain prices with a for loop.

executed on MentDB

Field slope, aspect and USLE LS factor from a DEM

From a 3x3 window of a 10 m elevation model, the script computes the Horn gradients, the slope angle, the aspect (direction of the steepest descent) and the USLE topographic factor LS of a 60 m slope, to plan contour farming against erosion.

Pipeline

01 · windowElevationsmatrix get→02 · slopeatan(hypot(dz/dx, dz/dy))math hypot→03 · aspectatan2 of the descentmath atan2→04 · angleSlope anglemath atan→05 · LSWischmeier-Smith LSmath sin

Script

admin

#1. Elevations in m of a 3x3 window of a 10 m DEM, rows north to south: a b c / d e f / g h i;
-> "[z]" "[[212.4, 211.8, 211.0], [211.6, 210.9, 210.1], [210.7, 210.0, 209.3]]";
-> "[cs]" 10;
#2. Horn gradients: dz/dx eastward and dz/dy northward;
-> "[dx]" (/ (- (+ (matrix get [z] 1 3) (+ (* 2 (matrix get [z] 2 3)) (matrix get [z] 3 3))) (+ (matrix get [z] 1 1) (+ (* 2 (matrix get [z] 2 1)) (matrix get [z] 3 1)))) (* 8 [cs]));
-> "[dy]" (/ (- (+ (matrix get [z] 1 1) (+ (* 2 (matrix get [z] 1 2)) (matrix get [z] 1 3))) (+ (matrix get [z] 3 1) (+ (* 2 (matrix get [z] 3 2)) (matrix get [z] 3 3)))) (* 8 [cs]));
#3. Slope angle atan of the gradient norm, aspect = direction of the steepest descent, clockwise from north;
-> "[sl]" (math atan (math hypot [dx] [dy]));
-> "[asp]" (math calc (concat "mod(" (math calc (concat "(" (math atan2 (- 0 [dx]) (- 0 [dy])) ")*180/pi")) "+360, 360)"));
#4. USLE topographic factor LS for a 60 m slope length, exponent m = 0.5 for a slope of 5 % or more;
-> "[s]" (math sin [sl]);
-> "[ls]" (* (math power (/ 60 22.13) 0.5) (+ (* 65.41 (* [s] [s])) (+ (* 4.56 [s]) 0.065)));
concat "slope = " (math round (* 100 (math tan [sl])) 2) " % = " (math round (math calc "[sl]*180/pi") 2) " deg, aspect = " (math round [asp] 1) " deg, LS = " (math round [ls] 3)

Result

mentdb

slope = 11.36 % = 6.48 deg, aspect = 140.4 deg, LS = 2.328

Go further

Multiply by R, K, C and P to get the USLE soil loss A in t/ha/year.
Run the window over a whole grid with two for loops and matrix set.

executed on MentDB

Seed lot germination: MGT, rate and T50

In a germination test of 50 maize seeds, the day of emergence of each germinated seed gives the germination percentage, the daily frequency table, the mean germination time, the mean germination rate (inverse of the harmonic mean time) and the T50 class.

Pipeline

01 · tableDaily and cumulative countsstat frequency→02 · peakPeak daystat mode→03 · rateHarmonic mean timestat harm_mean→04 · T50First class >= 50 %json select→05 · spreadVariance of the timesstat pvariance

Script

admin

#1. Germination test of a maize seed lot: day of radicle emergence of the 46 germinated seeds out of 50;
-> "[d]" "[3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 9, 9, 10]";
#2. Daily counts and cumulative germination;
json load "rx_p4d_germ" (stat frequency [d]);
#3. Germination %, mean germination time, mean germination rate = harmonic mean of the times;
-> "[gp]" (* 100 (/ 46 50));
-> "[mgt]" (stat mean [d]);
-> "[hm]" (stat harm_mean [d]);
#4. Day when half of the germinated seeds are out: first class with a cumulative frequency >= 0.5;
-> "[t50]" 0;
for (-> "[i]" 0) (< [i] (json count "rx_p4d_germ" "/")) (++ "[i]") {
	if (and (== [t50] 0) (>= (json select "rx_p4d_germ" (concat "/[" [i] "]/cumulative_frequency")) 0.5)) {-> "[t50]" (json select "rx_p4d_germ" (concat "/[" [i] "]/value"));};
};
concat "germination = " [gp] " %, peak day = " (stat mode [d]) ", MGT = " (math round [mgt] 2) " d, mean rate = 1/" (math round [hm] 2) " = " (math round (/ 1 [hm]) 3) " per day, T50 class = day " [t50] ", spread = " (math round (stat pvariance [d]) 2) " d2"

Result

mentdb

germination = 92 %, peak day = [5], MGT = 5.61 d, mean rate = 1/5.17 = 0.194 per day, T50 class = day 5, spread = 2.63 d2

Go further

Compare two seed treatments with stat mann_whitney on the germination days.
Fit a cumulative logistic curve with stat fit logistic to get a continuous T50.

executed on MentDB

Pesticide dissipation: SFO and FOMC DT50 / DT90

Herbicide residues in soil are fitted with the single first-order model and with the FOMC (Gustafson-Holden) model of the FOCUS kinetics guidance. The script gives both DT50 and the FOMC DT90, and tells if the decline is bi-phasic.

Pipeline

01 · SFOExponential fitstat fit→02 · FOMCC0/(1 + t/beta)^alphastat fit_model→03 · parametersalpha and betajson load select→04 · endpointsbeta(2^(1/alpha) - 1)math power

Script

admin

#1. Soil dissipation of a herbicide: residue in % of the applied dose;
-> "[t]" "[0, 3, 7, 14, 28, 56, 90]";
-> "[r]" "[100, 74, 57, 42, 29, 18, 12]";
#2. Single first-order model SFO and FOMC Gustafson-Holden model C0 / 1 + t/beta ^ alpha;
-> "[k]" (json load select "rx_p4d_sfo" (stat fit [t] [r] "exp") "/b");
-> "[fomc]" (stat fit_model "100/(1+x/be)^al" "x" [t] [r] "{\"al\": 1, \"be\": 10}");
-> "[al]" (json load select "rx_p4d_fomc" [fomc] "/parameters/al");
-> "[be]" (json select "rx_p4d_fomc" "/parameters/be");
#3. FOCUS endpoints: DT50 = beta.2^1/alpha - 1, DT90 = beta.10^1/alpha - 1;
-> "[dt50]" (* [be] (- (math power 2 (/ 1 [al])) 1));
-> "[dt90]" (* [be] (- (math power 10 (/ 1 [al])) 1));
concat "SFO DT50 = " (math round (/ (math log 2) (math abs [k])) 1) " d, FOMC alpha = " (math round [al] 3) ", beta = " (math round [be] 2) ", DT50 = " (math round [dt50] 1) " d, DT90 = " (math round [dt90] 1) " d, r2 FOMC = " (math round (json select "rx_p4d_fomc" "/r2") 4) ", FOMC needed: " (> (/ [dt90] [dt50]) 3.32)

Result

mentdb

SFO DT50 = 31.7 d, FOMC alpha = 0.744, beta = 6.23, DT50 = 9.6 d, DT90 = 131.2 d, r2 FOMC = 0.9996, FOMC needed: 1

Go further

Add the chi2 error % of FOCUS for each model (stat fit_model chi2).
Use the DT90 to compute the plateau concentration after repeated yearly applications.

executed on MentDB

NPK fertilizer blend with an LU decomposition

A 120-60-90 N-P2O5-K2O recommendation must be met with urea, DAP and potash, with 50 kg of ammonium nitrate already in the hopper. The script solves the nutrient balance with an LU decomposition and shows the nutrients brought by each product with a diagonal matrix.

Pipeline

01 · needRemaining needsvector sub→02 · decomposeP.A = L.Umatrix lu→03 · solveMasses of the productsmatrix solve→04 · supplycomposition.diag(masses)matrix diag→05 · totalTotal productstat sum

Script

admin

json unload "rx_p4d_fomc"

Result

mentdb

Go further

Add a cost vector and minimise it with an extra product: least squares or matrix pinv.
Convert P2O5 and K2O into elemental P and K with the factors 0.436 and 0.830.