jamesorr/co2sys-matlab

CO2SYS software for MATLAB (or octave) to compute variables of ocean CO2 system

co2sys-matlab's Issues

Deprecate mention + Add link to https://github.com/jonathansharp/CO2-System-Extd?

It would be good in my opinion if a link to CO2-System-Extd was provided here (somewhere in the ReadMe), especially if CO2SYS-MATLAB is not maintained anymore, in which case a note at the top of the ReadMe would also be welcome.

(I missed that there was another more up-to-date version out there myself.)

Broken ReadMe reference link

The link to http://cdiac.ornl.gov/oceans/co2rprt.html in the ReadMe is broken.

[Suggestion] More robust root finding for pH evaluation

Not sure this is useful to anyone but me, but I found out that sometimes the Newton solver for evaluating pH hangs when the other tracer values are outside of some acceptable range. I think this might be due to the current backtracking implementation where the Newton step is halved until the pH step deltapH is less than 1:

CO2SYS-MATLAB/src/CO2SYS.m

Lines 1525 to 1530 in b62ad9e

	deltapH = Residual./Slope; % this is Newton's method
	% to keep the jump from being too big;
	while any(abs(deltapH) > 1)
	FF=abs(deltapH)>1; deltapH(FF)=deltapH(FF)./2;
	end
	pHx = pHx + deltapH; % Is on the same scale as K1 and K2 were calculated...

However, I think that even an overshoot in deltapH of less than 1 can still be too large and result in some never-ending cycle (illustrated below).

I think this is because the criteria for reducing the step size should be on the next residuals. Thus, maybe this backtracking line search could be replaced with one for which the criteria would be on the next residual instead?

Another option would be to write a function that returns the residual for a given pH and then use an existing Newton-solver that already has that criteria implemented (e.g., C. T. Kelley's nsold.m, which uses a quadratic approximation of the residual along the Newton-step direction for a rather efficient backtracking IMHO).

Separate the Julia version?

I could not submit an issue on @mvdh7's CO2System.jl repository likely because it is a fork of this repository. Would it make sense to split the Julia version completely and cut the umbilical cord tying it to this MATLAB version?

FWIW, for the Julia version, I wanted to suggest converting the main function, which takes scalars or arrays as input, into a function that takes only scalars as inputs. I think a scalar function, called e.g., co2sys(args...), would have many advantages over the current version (in Julia). Note that a scalar function would be easy to broadcast to arrays with Julia's dot syntax, e.g., for a simple example:

julia> f(x,y) = x + y # define the function on scalars
f (generic function with 1 method)

julia> f.([1, 2, 3], 4) # easily broadcast to those inputs that are arrays
3-element Vector{Int64}:
 5
 6
 7

Another potential gain is that the code would be simpler (e.g., less lines of code because no need for all the array treatment, the allocations, and so on). Furthermore, it might be computationally more efficient because, e.g., no need to reallocate the boolean indices array F dozens of times. (Although efficiency would need to be benchmarked against the current vectorized version I guess.)

Don't hesitate to close if these are stupid ideas! 😅

Broken van Heuven et al. (2011) link in ReadMe

https://doi.org/10.3334/CDIAC/otg.CO2SYS_MATLAB_v1.1 is not working.