Underwater optical wireless communication (Matlab simulation)

I'm having trouble on getting the result as shown in the picture. I'm using the equation and the table to do the simulation. I wonder where is my mistake in my coding, I wasn't able to get the same result as the picture. My code is shown in coding.txt.

(Edit by John) I've added the short code here:

bo=0.0429;
s=-0.103e-3;
h=11.87;
zmax=115.4;
cch_max=0.708;
std=h./sqrt((2*pi)*(cch_max-bo-s.*zmax));
z=0:1:250;
ccz=(bo+s.*z)+(h./(std.*sqrt(2*pi)))*exp((-(z-zmax).^2)./(2*std.^2));
plot(ccz,z);

[John D'Errico]

I was trying to read your question, then look back at your short code fragment. Too hard to go back and forth, and too easy just to insert the code into your question. So I did that for you, which I imagine will make it easier for someone to answer your question.

Now that I can more easily look at your code, I see that you use std as a variable. That may be a dangerous thing because std is such a useful function. Generally try to avoid the use of variable names that step on top of (overload) existing function names.

Answer 1

[John D'Errico]

Looking at your code, it seems essentially correct. It is a fairly standard normal distribution, with mode at zmax, standard deviation of sigma, and a y offset of (bo+s.*z). You could write this more simply just using the normpdf function if you have the stats toolbox.

However, I see one significant inconsistency.

The figure you have shown (a) has z going from -250 to 0. zmax is listed in the table as a POSITIVE value, thus 115.4. But zmax is the point where the normal PDF will beat its maximum value. And in the figure, that value is clearly negative. As well, you use z going from 0 to 250. So it is not clear if they simply dropped the sign on zmax by accident or on purpose, but the sign of zmax is important. If you change the sign of zmax, it might change the computation of sigma.

So you should probably resolve the question of the sign of z and zmax.

[Salad]

>> bo=0.0429;
>> s=-0.103e-3;
>> h=11.87;
>> zmax=-115.4;
>> cch_max=0.708;
>> std_1=h./sqrt((2*pi)*(cch_max-bo-s.*zmax));
>> z=0:-1:-250;
>> ccz=(bo+s.*z)+(h./(std_1.*sqrt(2*pi)))*exp((-(z-zmax).^2)./(2*(std_1.^2)));
>> plot(ccz,z);

I changed the sign of the z and the zmax to negative, and it gave me

So I wonder what went wrong, because of the given equation itself, the parameter in the table or my coding?