calculate integration in matlab for environmental physics data series

Question

The following diagram is a schematic of a lake, and the equation illustrates how to calculate the effective heat flux of a lake.

where S is a vector of surface fluxes, q is short wave radiation, h is depth of the mixed layer, and z is the depth of the lake. For example:

q0 = 400+(1-400).*rand(100,1); % This is the short wave radiation
kd = 0.8; % extinction coefficient
h = 10; % depth of the surface mixed layer
for i = 1:length(q0); % loop for calculating short wave radiation at depth h
    qh(i) = q0(i).*exp(-kd*h); % here, qh is calculated according to the Lambert Beer law
end

given

dz = 0.5

and z varies from 0 (surface) to depth h in increments of dz i.e.

z = 0:dz:h

how would I calculate the last portion of this equation in matlab i.e. how to calculate q at depth z between the surface and h? which is expressed here as an integral?

Apologies if this should be on another one of the stack overflow forums but it seems more related to programming than pure physics or maths questions.

Answer 1

To integrate this correctly, you will need compute all values of q(z) in the range [0, h]. If q0 and qh are N-by-1 column vectors, this means that q should be an N-by-M matrix, where M is the number of sample points in the range [0, h].

First, lets define z properly:

z = linspace(0, h, 200); %// M=200, but it's an arbitrary number to your choosing

The computation of q can be reduced to:

q = q0 * exp(-kd * z);

and qh actually equals to the final column of q, i.e q(:, end).

The integral itself can be approximated to a sum and computed using sum:

dz = z(2) - z(1);
I = sum(q, 2) * dz;

P.S.

Since q(z) = e^{(-k_d ·z)}, it's simple enough for you to compute the integral analytically:

I = q0 * (1 - exp(-kd * h)) / kd;