brews / co2syspy

A Python interpretation of CO2SYS
GNU General Public License v3.0
3 stars 0 forks source link

Broadcasting and array input #7

Open brews opened 6 years ago

brews commented 6 years ago

Originally posted by @aholdsworth in #4 :

The code only seems to work for scalar inputs. I would have liked to use it for 3d fields, but it's not yet working for even the following simple test case. This is example 6 found here (https://github.com/jamesorr/CO2SYS-MATLAB/blob/master/notebooks/CO2SYS-Matlab_tutorial.ipynb) for the matlab version:


par1type = 1; # The first parameter supplied is of type "1", which is "alkalinity"
par1 = 2300; # value of the first parameter
par2type = 2; # The first parameter supplied is of type "1", which is "DIC"
par2 =np.arange(2000,2200,5) # value of the second parameter, which is a long vector of different DIC's!
sal = 35; # Salinity of the sample
tempin = 10; # Temperature at input conditions
presin = 0; # Pressure at input conditions
tempout = 0; # Temperature at output conditions - doesn't matter in this example
presout = 0; # Pressure at output conditions - doesn't matter in this example
sil = 50; # Concentration of silicate in the sample (in umol/kg)
po4 = 2; # Concentration of phosphate in the sample (in umol/kg)
pHscale = 1; # pH scale at which the input pH is reported ("1" means "Total Scale") - doesn't matter in this example
k1k2c = 4; # Choice of dissociation constants K1 and K2 ("4" means "Mehrbach refit")
kso4c = 1; # Choice of KSO4 ("1" means "Dickson")
print('par2', par2)

A =co2sys.CO2SYS(par1,par2,par1type,par2type,sal,tempin,tempout,presin,presout,sil,po4,pHscale,k1k2c,kso4c)

We've only written co2syspy for scalar input and broadcasting support would be nice.

In the meantime, as a quick fix, the simplest option would be to use something like a for-loop for the different input.

aholdsworth commented 6 years ago

Thanks for getting back to me on this. Unfortunately, a for loop is prohibitively inefficient in my case as the inputs are very large. I'll let you know if I modify your code to support broadcasting. I haven't decided how to proceed just yet.