swDynamicHeight() fails with `integrate` error

[richardsc]

I can't tell exactly what about the profile is causing the error, but I'm getting this:

swDynamicHeight(ctd)
Error in integrate(integrand, 0, referencePressure) : 
  maximum number of subdivisions reached

It reads and plots fine:

Seems like there's a default subdivisions arg in integrate() (as well as a bunch of others), so I wonder if it might be a good idea to add a ... arg to swDynamicHeight()? I've emailed an .rda that contains the object. Note that it wasn't created from a .cnv but from doing as.ctd(), though I don't see why that should matter (and in fact, other profiles read the same way work fine, as this is part of a section).

[richardsc]

Hmmm. Adding a ... to swDynamicHeight() seems to fix it for that profile, but when I try and go through the rest of the ones I have, I get other weird integrate() errors, such as:

swDynamicHeight(sec@data$station[[7]], sub=200)
Error in integrate(integrand, 0, referencePressure, subdivisions = sub) (from 03.R#21) : 
  extremely bad integrand behaviour

Note that the sub arg is just my way of passing the subdivisions to integrate().

With a bit of noodling, I can get it to work by adjusting the number of subdivisions:

swDynamicHeight(sec@data$station[[7]], sub=150)
[1] 0.3016795

but I don't really know why.

I'll send you a file that has the entire section.

I wonder if it might make sense to have a try() around the integrate, and if it fails for some reason to spit a warning and then just do something more simple (like perhaps use integrateTrapezoid())

Also, this is not timely or high priority. I can code around it if I need to. I just report 'em when I see 'em :)

[dankelley]

I've done a couple of things in the develop branch (listed below) and it seems fine on the sample code on the 499.R test code.

1. The integration is now in scaled variables. That seemed to help several test cases. This is not uncommon: numerical procedures often are written to work best with variables that are O(1).
2. swDynamicHeight() now has arguments for subdivisions and relative tolerance. The default for the former is now 500, up a lot from the integrate() default of 100, but the default for the latter is the same as for integrate() because I imagine that is a value that is well thought out. (On a 64 bit machine it is 0.1 percent, which in fact is probably better than anyone would hope to get, a millimetre in height for the Gulf Stream signal say.)

@richardsc -- if this seems ok now, please close; otherwise, please let me know of further problems. Thanks!

[richardsc]

Great, thanks! Works for me now on the data I sent.