It seems that Octave crashes when the nr of measurements is >500 (roughly).
I've been playing a bit with the software and reading the documentation. Both are very impressive.
Currently the number of modes are predetermined. If my memory doesn't fail me, the fact that spherical harmonics are orthogonal functions means that every mode can be determined in isolation. Adding another node in de fit should not alter the already calculated coefficients. Just like with adding (co)sines with higher frequencies to fit e.g. a square wave and unlike polynomial fitting. Does that also hold when there is noise on the measurements and you're doing a least squares fit?
I was wondering if we could keep adding modes and see the residue go down until and stop when the solution starts diverging. If that works, one could do that per frequency as for low frequencies it is unlikely that you need a lot of modes.