Two variable function computation

Hello,

I am trying to reproduce the Matlab plot I found in a research paper. The function in P(t,z) (it is time and space dependent). P(t,z) have two parameters (alpha and beta) which are frequency dependent.

I have two loops for computing P(t,z). My code is below, however, I can't find a similar solution which what is mentioned. This is my code

clc
clear all
close all

f=linspace(0,20*10^6,100);
t=1./f;
x=linspace(0,0.15*10^-2,100);

alpha=(0.45*10^-11)*f.^2;
beta=(0.15*10^-4)*f;

figure(1)
subplot(2,1,1)
plot(f/10^6,alpha/100)
ylabel('alpha')
subplot(2,1,2)
plot(f/10^6,beta/100)
ylabel('beta')
p = zeros(length(x),length(t));

for it= 1:length(t)
    for ix=1:length(x)
        p(ix,it)= exp(-alpha(it).*x(ix)).*exp(-i.*1.*(f(it).*t(it)-beta(it).*t(it)));
    end
end

figure(2)
plot3(x,t,p)

Any help will be appriciated.

Thanks

[Walter Roberson]

Note that your p is complex, so you should abs() before plotting the magnitude.

You are changing alpha and omega and beta over time, which is not part of the formula. That p formula is for some particular alpha and beta and omega.

Your definition of t is non-linear, but the plot on the right appears to use linear time.

In the right hand plot, the Distance does not appear to have a scaling factor, but you have used 10^-2 scaling factor in your definition of x.

[friet]

Hello,

Thanks for your comment. I have made some modification based on your kind comment. I have selected omega as 10MHz and select the corresponding values of alpha and beta. The factor 10^-2 is change from cm to m. All the units are in for x are in m and time in sec.

clc
clear all
close all

f=linspace(0,20*10^6,100);
t=linspace(0,1*10^-6,100);

x=linspace(0,0.15*10^-2,100);

alpha=(0.45*10^-11)*f.^2;
beta=(0.15*10^-4)*f;

omega=10*10^6;
a=5.5*100;
b=1.8*100;

figure(1)
subplot(2,1,1)
plot(f/10^6,alpha/100)
ylabel('alpha')
subplot(2,1,2)
plot(f/10^6,beta/100)
ylabel('beta')
p = zeros(length(x),length(t));

for it= 1:length(t)
    for ix=1:length(x)
          p(ix,it)= exp(-a.*x(ix)).*exp(-i.*1.*(omega.*t(it) - b.*t(it)));
    end
end

tt= repmat(t,100,1);
xx= repmat(x,100,1);
figure(2)
waterfall(xx,tt,abs(p))

But still, can't find the same solution. By the way, how can I change the origins of the 3d plot, the distance and time to be start at 0 and the magnitude origin separately as it was plotted in the original plot.

Thanks

Answer 1

[Honglei Chen]

You may want to check out waterfall

https://www.mathworks.com/help/matlab/ref/waterfall.html

HTH

Answer 2

[Walter Roberson]

Your expression just does not act like that, at least not over those parameter ranges.

See attached that explores a number of possibilities. (Note the vectorized implementations of the equations.) I put up several isosurface plots to explore the way the surface grows along the 4 different parameters. The last couple of plots explore the possibility that your x is too small or too large. The last one does suggest that possibly the x range should be much smaller; for example if you were to use x/200 then it would have room for only one peak.

On the isosurface plots notice that I plot the real component of p rather than the absolute magnitude. The absolute magnitude are completely boring; plotting the real component at least allows you to track through a couple of phase cycles.