Index exceeds matrix dimensions

Hi, my tutor is sick of my questions, please help!

    for j=2:31
        dz(j)=abs(depth(j)-depth(j-1))
    end

    for i=1:360
    for j = 1:180
    for k=1:31
        intFe68 = Fe68(i,j,k)  *dy(360) * dx(j) * dz(j)
    end
    end
    end

    nansum(intFe68)

Fe68 has dimensions 360x180x31. It is iron data from a model. I can contour a hydrothermal plume by plotting Fe68 against depth and longitude, below I am trying to get one value for the size of the plume in m^3 but matlab doesn't like the second loop due to the 'index exceeding the matrix dimensions'. Im also not sure if I need the first loop, or what to put in the second loop if I do, dz(j), dz(:), dz(k)..

Any help would be very appreciated!

Alex

[Chad Greene]

Hi Alex,

What exactly do you mean by "integrate up to get one value"? Perhaps describing the context and what you're trying to do will help. This looks like maybe a month's worth of ocean iron content data? What exactly are you trying to calculate?

[Alexandra Ridgway]

Hi Chad, it's iron data from a model. I can plot a hydrothermal plume by plotting Fe68 against longitude however I want to quantify the size of the plume. (Im not sure if that would involve integrating...) Does that help at all? Thank you for trying to help me!

Answer 1

[Walter Roberson]

In your loop

    for j=2:31
        dz(j)=abs(depth(j)-depth(j-1))
    end

you initialize dz up to dz(31)

In

    for i=1:360
    for j = 1:180
    for k=1:31
        intFe68 = Fe68(i,j,k)  *dy(360) * dx(j) * dz(j)
    end
    end
    end

you access dz(j) with j from 1 to 180, so you need up to dz(180) to have been defined.

But notice you have k = 1 : 31 so perhaps you should be accessing dz(k) ?

[Alexandra Ridgway]

Hi Walter, My depth values increase in different increments, hence the reasoning for the first loop to work out the difference between each value. If I put dz(k) will it not just use the numbers 1-31 rather than the difference in the values between each depth?

Thank you :)

Alexandra

[Walter Roberson]

If you only have 31 different depths, then why are you trying to access 180 of them? dz(j) when j = 1:180

Perhaps your code should be

          intFe68(i,j,k) = Fe68(i,j,k)  *dy(360) * dx(j) * dz(k)

but I have to wonder then whether dy(360) should be dy(i) ?

          intFe68(i,j,k) = Fe68(i,j,k) * dy(i) * dx(j) * dz(k)

possibly?

Once you figure out the correct multiplications then we could proceed to optimize the code.

... though with your references to getting one value and your use of "intFe68", I wonder whether you should be using trapz() instead of what you are doing now.

[Alexandra Ridgway]

What does the trapz() function do? And right, ideally, as it makes sense in my head I would like to do this,

for l=2:31
    dz(l)=abs(depth(l)-depth(l-1))
end

for i=1:360
for j = 1:180
for k=1:31
    intFe68(i,j,k) = Fe68(i,j,k)  *dy(i) * dx(j) * dz(:)
end
end
end

(then maybe trapz(intFe68) ?)

Thank you again!

(I also have written another question which may be simpler to solve first)

[Walter Roberson]

Are you trying to do numeric integration of the entire Fe68 array under the assumption that the distance between first coordinates is dy, the distance between second coordinates is dx, and the distance between third coordinates is dz ? If you are, then you should be using trapz() which is the numeric integration routine.

dy = gradient(y_coordinates);
dx = gradient(x_coordinates);
dz = gradient(z_coordinates);

numeric_integral = trapz( dy, trapz(dx, trapz(dz, Fe68, 3), 2), 1);

[Alexandra Ridgway]

Hi Walter, x is longitude, y is latitude and z is depth. With that knowledge does it make sense that that would be what I'm trying to do?! Im so confused sorry! Thank you for your help. (and would this code be entirely instead of what I've done already?)

[Walter Roberson]

If Fe68(i,j,k) is the concentration of iron at i, j, k, then you could estimate the total iron by using the gradient and trapz that I showed above, without any explicit looping required.

[Alexandra Ridgway]

Hi, yes I understand, however numeric_integral is returning NaN values only :/ Thank you for your help.

[Walter Roberson]

NaN would be returned if any element in Fe68 is NaN, or if any element in the the gradients was NaN, and could be potentially be returned if there are infinities in the data.