How to run a function through multiple arrays
4 views (last 30 days)
Show older comments
I have a function that requires input from arrays representing data from around the globe (see below). The function is designed to solve for one value at a time, but I need to solve for all the points in my global data set (i.e. I have arrays for temp, s, alk, dic that correspond to specific locations around the globe). What is the best way to loop this function through all the arrays?
function[pco2,pH,co2,hco3,co3] = co3eq(temp,s,z,alk,dic)
% Function to calculate pCO2 (microatm), pH, H2CO3* (micromol/kg), bicarbonate ion (micromol/kg), and
% carbonate ion (micromol/kg) concentrations from Alkalinity (micromol/kg) and
% DIC (micromol/kg). This function requires temperature (degrees C),
% salinity (ppt), and depth (m), and uses the equations in Zeebe and
% Wolf-Gladrow (2000).
t = temp + 273.15;
Pr = z/10;
alk = alk*.000001;
dic = dic*.000001;
R = 83.131;
% Calculate total borate from chlorinity
tbor = .000416 * s / 35.0;
% Calculate Henry's Law coeffieicent, K0 (Weiss, 1974)
U1 = -60.2409 + 93.4517 * (100/t) + 23.3585*log(t/100); % note that *log* in matlab is the natural log, ln
U2 = s * (.023517 - .023656 * (t/100) + .0047036 * (t/100)^2);
KH = exp(U1+U2);
% Calculate KB from temp and salinity (Dickson, 1990)
KB = exp((-8966.9 - 2890.53 * s^0.5 - 77.942 * s + 1.728 * s^1.5 - 0.0996 * s^2)/t ...
+ 148.0248 + 137.1942 * s^0.5 + 1.62142 * s - (24.4344 + 25.085 * s^0.5 + ...
0.2474 * s) * log(t) + 0.053105 * s^0.5 * t);
% Calculate K1 and K2 (Luecker et al., 2000)
K1 = 10^(-(3633.86/t - 61.2172 + 9.67770 * log(t) - 0.011555 * s + 0.0001152 * s^2));
K2 = 10^(-(471.78/t + 25.92990 - 3.16967 * log(t) - 0.01781*s + 0.0001122 * s^2));
% Pressure variation of K1, K2, and KB (Millero, 1995)
dvB = -29.48 + 0.1622 * temp - .002608 * (temp)^2;
dv1 = -25.50 + 0.1271 * temp;
dv2 = -15.82 - 0.0219 * temp;
dkB = -.00284;
dk1 = -.00308 + 0.0000877 * temp;
dk2 = .00113 - .0001475 * temp;
KB = (exp(-(dvB / (R * t)) * Pr + (0.5 * dkB / (R * t)) * Pr^2)) * KB;
K1 = (exp(-(dv1 / (R * t)) * Pr + (0.5 * dk1 / (R * t)) * Pr^2)) * K1;
K2 = (exp(-(dv2 / (R * t)) * Pr + (0.5 * dk2 / (R * t)) * Pr^2)) * K2;
% temperature dependence of KW (DoE, 1994)
KW1 = 148.96502 - 13847.26 / t - 23.65218 * log(t);
KW2 = (118.67 / t - 5.977 + 1.0495 * log(t)) * s^.5 - 0.01615 * s;
KW = exp(KW1 + KW2);
% solve for H ion (Zeebe and Wolf-Gladrow, 2000)
a1 = 1;
a2 = (alk + KB + K1);
a3 = (alk * KB - KB * tbor - KW + alk * K1 + K1 * KB + K1 * K2 - dic * K1);
a4 = (-KW * KB + alk * KB * K1 - KB * tbor * K1 - KW * K1 + alk * K1 * K2 ...
+ KB * K1 * K2 - dic * KB * K1 - 2 * dic * K1 * K2);
a5 = (-KW * KB * K1 + alk * KB * K1 * K2 - KW * K1 * K2 - KB * tbor * K1 * K2 ...
- 2 * dic * KB * K1 * K2);
a6 = -KB * KW * K1 * K2;
p = [a1 a2 a3 a4 a5 a6];
r = roots(p);
h = max(real(r));
% calculate bicarbonate, carbonate, and aqueous CO2 using DIC, Alk, and H+
format short g;
hco3 = dic/ (1 + h/K1 + K2/h) * 1e6;
co3 = dic / (1 + h/K2 + h * h / (K1 * K2)) * 1e6;
co2 = dic / (1 + K1/h + K1 * K2 / (h * h)) * 1e6;
pco2 = co2 / KH ;
pH = -log10(h);
% calculate B(OH)4 and OH
BOH4 = KB * tbor / (h + KB);
OH = KW / h;
% recalculate DIC and Alk to check calculations
Ct = (hco3 + co3 + co2) * 1e6;
At = (hco3 + 2*co3 + BOH4 + OH - h) * 1e6;
0 Comments
Answers (1)
Walter Roberson
on 18 Apr 2016
Everything before the roots() can be vectorized, slightly rewritten to be able to work on arrays of values (such as might be produced with ndgrid). However, to find the roots of an array of numbers, you are going to need to use arrayfun() with 'Uniform', 0. Once you have that, cellfun to find the max of the real parts of the roots, and you are back to a grid of values, and can proceed in vectorized manner.
2 Comments
Walter Roberson
on 18 Apr 2016
[r] = arrayfun(@roots,p,'UniformOutput',false);
However, that is not really what you want. You want
r = arrayfun(@A1,A2,A3,A4,A5,A6) roots([A1,A2,A3,A4,A5,A6]), a1, a2, a3, a4, a5, a6, 'Uniform', 0);
Your arrays are probably 2D, so when you used
p = [a1, a2, a3, a4, a5, a6]
then you got out a wide 2D array, and then your arrayfun was iterating over each element of that 2D array. If you had used
p = cat(3, a1, a2, a3, a4, a5, a6);
then you could at least worked a stack at a time, roots(squeeze(p(J,K,:)) but that is not that convenient. The work-arounds to be able to nicely do one "element" at a time such as by way of cellfun are more costly then the above code.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!