% The data below are total production and ef ratios generated from the Laws model at the 
% indicated temperatures.  The units of production are mg N per cubic meter per day.
% This Matlab program calculates ef ratios as a function of total production and temperature.


% In order to interpolate these data with Matlab, I need the production values to be 
% equally spaced and the same for each temperature.  In order to do that, I first 
% did a one-dimensional interpolation (using the program interp1) of the ef ratios at 
% each temperature so that they corresponded to total production values equally spaced 
% between 0.07 and 94.  The two-dimensional interpolation was then done using the program
% interp2.  Both interp1 and interp2 are Matlab routines.  Note that the interpolation 
% uses temperature and the log of total production as the two independent variables.

% In the example given in the last line, the interpolated ef value is 0.32 at a 
% temperature of 19 degrees and total production of 10 mg N per cubic meter per day.


% Data at 0 degrees C

p0=[94.6 47.3 23.7 12 6.12 3.21 1.75 .965 .586 .295 .143 .0672];
f0=[.6765 .6765 .6765 .6674 .6536 .6237 .5718 .5183 .4268 .1061 .1092 .1163];
f2l0=[1512 1511 1510 1509 1508 1507 1419 1372 1302 977 942 907]/10000;
n0=14*75/1000./(1-f2l0);

% Data at 2.5 degrees C

p025=[95.1 47.6 23.8 12.2 6.22 3.25 1.84 1.42 1.19 .592 .29 .139 .0638];
f025=[.6729 .6729 .6729 .6571 .6432 .6158 .5447 .5282 .1046 .1056 .108 .1124 .1225];
f2l025=[1629 1628 1627 1581 1570 1558 1465 1433 1177 1037 1009 960 937]/10000;
n025=14*75/1000./(1-f2l025);

% Data at 5 degrees C

p5=[94.9 47.4 23.8 12.3 6.33 3.44 1.9 1.19 .587 .283 .132 .0578];
f5=[.6745 .6745 .6714 .651 .6322 .5822 .526 .1051 .1066 .1103 .1181 .1352];
f2l5=[1686 1685 1684 1683 1665 1609 1572 1221 1123 1042 993 944]/10000;
n5=14*75/1000./(1-f2l5);


% Data at 7.5 degrees C

p75=[94.5 47.9 24.2 12.6 6.5 3.57 2.80 2.38 1.18 .577 .275 .124 .0534];
f75=[.6774 .6679 .6611 .636 .6157 .5598 .5361 .1051 .1058 .1083 .1137 .1256 .1464];
f2l75=[1885 1884 1883 1842 1786 1702 1647 1195 1174 1153 1088 1014 940]/10000;
n75=14*75/1000./(1-f2l75);

% Data at 10 degrees C

p10=[96 48 24.8 12.9 6.93 3.77 3.01 2.37 1.16 .559 .2626 .112 .0494];
f10=[.6669 .6669 .645 .6224 .5769 .53 .497 .1053 .1074 .1119 .119 .1395 .1581];
f2l10=[2257 2256 2255 2012 1872 1744 1721 1395 1228 1130 1129 1060 893]/10000;
n10=14*75/1000./(1-f2l10);


% Data at 12.5 degrees C

p125=[118 58.8 29.2 14.6 7.36 4.13 2.83 2.35 1.15 .544 .244 .104 .0472];
f125=[.5439 .5439 .5472 .5472 .5433 .4846 .1061 .1064 .1091 .1148 .1281 .1502 .1657];
f2l125=[2163 2162 2159 2158 2088 1902 1350 1344 1260 1233 1177 1035 1034]/10000;
n125=14*75/1000./(1-f2l125);

% Data at 15 degrees C

p15=[145 71.3 35.7 18.2 9.15 4.68 3.28 2.31 1.12 .513 .224 .0985 .0462];
f15=[.4405 .4485 .4485 .4400 .4373 .428 .1067 .108 .1115 .1219 .1393 .1586 .1692];
f2l15=[2223 2222 2221 2120 2065 2023 1477 1388 1340 1258 1177 991 670]/10000;
n15=14*75/1000./(1-f2l15);

% Data at 20 degrees C

p20=[225 115 55.8 27.4 13.9 6.61 4.48 2.19 .987 .439 .198 .0918 .0466];
f20=[.2849 .2778 .2869 .2915 .2870 .3025 .1116 .1141 .1266 .1424 .1577 .1702 .1678];
f2l20=[2262 2261 2260 2259 2159 2158 1547 1546 1456 1337 1128 760 486]/10000;
n20=14*75/1000./(1-f2l20);

% Data at 25 degrees C

p25=[335 166 82.6 40.9 19.5 9.32 5.12 4.01 1.6 .927 .394 .186 .0925 .0467];
f25=[.191 .1923 .1936 .1957 .2051 .2146 .127 .1247 .125 .1349 .1586 .1677 .1689 .1672];
f2l25=[2663 2635 2634 2580 2572 2449 1822 1821 1742 1626 1351 851 565 377]/10000;
n25=14*75/1000./(1-f2l25);

% Data at 28 degrees C

p28=[418 207 103 50.4 24.4 11.4 6.47 3.7 1.89 .833 .379 .187 .093 .0464];
f28=[.1532 .1544 .1558 .1588 .1639 .1754 .1546 .1353 .1336 .15 .165 .1674 .168 .1685];
f2l28=[3400 3372 3223 3058 2872 2674 2274 2023 1819 1581 1088 721 488 343]/10000;
n28=14*75/1000./(1-f2l28);

% Data at 30 degrees C

p30=[531 262 128 61.3 28.3 19.4 15 11.9 9.74 7.04 3.94 1.79 .79 .371 .189 .0942 .0468];
f30=[.1205 .1221 .1253 .1306 .1414 .1550 .1667 .168 .154 .142 .127 .1397 .1582 .1683 .1657 .1659 .1669];
f2l30=[3540 3490 3395 3256 3116 3115 3081 2942 2698 2419 2072 1816 1442 986 640 447 314]/10000;
n30=14*75/1000./(1-f2l30);

% Interpolation for Production-T-ef plot

Pmax=log(94);
Pmin=log(0.07);
D=(Pmax-Pmin)/100;
P=Pmin:D:Pmax;

ef0=interp1(log(p0),f0,P);
ef025=interp1(log(p025),f025,P);
ef5=interp1(log(p5),f5,P);
ef75=interp1(log(p75),f75,P);
ef10=interp1(log(p10),f10,P);
ef125=interp1(log(p125),f125,P);
ef15=interp1(log(p15),f15,P);
ef20=interp1(log(p20),f20,P);
ef25=interp1(log(p25),f25,P);
ef28=interp1(log(p28),f28,P);
ef30=interp1(log(p30),f30,P);

% Interpolation for DIN-T-ef plot

nmax=14*75/1000./(1-0.03);
nmin=14*75/1000./(1-0.36);
D=(nmax-nmin)/100;
N=nmin:D:nmax;

EF0=interp1(n0,f0,N);
EF025=interp1(n025,f025,N);
EF5=interp1(n5,f5,N);
EF75=interp1(n75,f75,N);
EF10=interp1(n10,f10,N);
EF125=interp1(n125,f125,N);
EF15=interp1(n15,f15,N);
EF20=interp1(n20,f20,N);
EF25=interp1(n25,f25,N);
EF28=interp1(n28,f28,N);
EF30=interp1(n30,f30,N);

clf;

ef=[ef0' ef025' ef5' ef75' ef10' ef125' ef15' ef20' ef25' ef28' ef30'];
EF=[EF0' EF025' EF5' EF75' EF10' EF125' EF15' EF20' EF25' EF28' EF30'];
T=[0 2.5 5 7.5 10 12.5 15 20 25 28 30];
LP=length(ef0);

% Here is a sample interpolation

EFI=interp2(T,P,ef,19,log(10))