Energy-Pathways-Group / GLOceanKit

Tools for physical oceanography in Matlab and Objective-C
21 stars 11 forks source link

Problems with InternalModes and constant stratification #16

Open NeilTandon opened 6 years ago

NeilTandon commented 6 years ago

I haven't been able to get InternalModes to work with constant stratification. The first issue I encountered was because the first argument fed to InternalModesConstantStratification is [N0, min(rho(zIn))], but rho is typically a function handle, so min(rho) returns an error. I attempted to fix this by replacing min(rho) with min(rho(zIn)) . So that error is gone, but now I get the following output:

Initialization detected that you are using constant stratification. The modes will now be computed using the analytical form. If you would like to override this behavior, specify the method parameter.
Using the analytical form for constant stratification N0=0.006777884
Error using .*
Matrix dimensions must agree.

Error in InternalModesConstantStratification/BaroclinicModesWithEigenvalue (line 146)
            G = A .*  sin(k_z .* (self.z + self.Lz));

Error in InternalModesConstantStratification/ModesAtFrequency (line 79)
            [F,G] = self.BaroclinicModesWithEigenvalue(k_z,h);

Error in InternalModes/ModesAtFrequency (line 419)
            [F,G,h,k] = self.internalModes.ModesAtFrequency( omega );

Error in batch85b_modal_decomp_ECCO (line 163)
    [Ftmp,G,h,k] = im.ModesAtFrequency(omega);
NeilTandon commented 6 years ago

Also, if I force my code to use the wkbAdaptiveSpectral method, everything works fine.

JeffreyEarly commented 6 years ago

I'm confused, what exactly are you feeding into InternalModesConstantStratification?

I'm trying,

z = linspace(-5000, 0, 128)';
im = InternalModesConstantStratification([5.2e-3 1025], [-5000 0], z, 35);
[F,G,h,k] = im.ModesAtFrequency(im.f0);
JeffreyEarly commented 6 years ago

I did mess with that interface recently, so it wouldn't surprise that something is broken---it's also possible I fixed this recently!

NeilTandon commented 6 years ago

Sorry, I should have provided more detail. Here is a simple code that produces the error I mentioned above:

N0 = 5.2e-3;
g = 9.81;
rho0 = 1025;
rho = @(z) -(N0*N0*rho0/g)*z + rho0;
L = -5000;
zOut = linspace(L,0,200)';
latitude = 33;
im = InternalModes(rho,[L 0],zOut,latitude);
[F,G,h,k] = im.ModesAtFrequency(5*im.f0);

The error is produced by the last command, im.ModesAtFrequency. When using constant stratification, is that the right command for obtaining F and G?

JeffreyEarly commented 6 years ago

Thank you! That let me repeat the error. I just updated the code with a fix---trying pulling the changes and see if that fixes it for you.

NeilTandon commented 6 years ago

Re-downloaded the code, but still getting the same error. :( BTW, I'm using MATLAB R2013b, in case that matters.

JeffreyEarly commented 6 years ago

Ah ha! Matlab 2013b is the problem---I'm using the new (as of 2016b) array expansion operators. Shoot---any chance you can easily upgrade?

NeilTandon commented 6 years ago

Ah, ok. Unfortunately, it would not be easy to upgrade, since I'm running MATLAB on a computer that I don't administer. Can you recommend a workaround? My current workaround is to just force the code to use the wkbAdaptiveSpectral method, which works fine in my MATLAB version.

JeffreyEarly commented 6 years ago

Yes, there's definitely no problem overriding it to use wkbAdaptiveSpectral. The only real compromise is speed... and that you have to watch out for garbage modes that might occur, but it's generally quite good, especially in smooth stratification profiles.