Доброго времени суток всем!
Дорогие друзья, возможно, моя просьба покажется очень странной, но все же решусь задать, потому что бьюсь над ней уже не одну неделю :)
Нужно смоделировать биполярный акустический импульс, который возникает в воде от каскада заряженных частиц. Сделать это в Матлабе, и смоделировать его распространение. Есть статьи, где описаны формулы и параметризация, но результат у меня получается какой-то явно неправильный.
Пожалуйста, помогите!
Вот примерно так выглядит формула:
Задача цилиндрически симметрична, поэтому у меня 2 переменные: z и x (радиальная).
Собственно, код
Код:
x1= 0:0.1:10;
E = 10^10; % total deposited energy
a = 0.125;
b = 0.002;
A = (2*E*(a*b)^0.5/3.14)^(0.5);
al = 1.2*0.001;
Cp = 3.8*1000;
cs = 150000;
t = 0:(1/fs1):T; % array of time moments - for each moment i calcutate the z coordinate
for k = 1:3;
e1 = zeros(size(x1));
e2 = zeros(size(x1));
for i=1:length(x1)
z1(i)=z-(cs^2*t(k)^2 - (x-x1(i))^2)^(0.5);
z2(i)=z+(cs^2*t(k)^2 - (x-x1(i))^2)^(0.5);
e1(i) = A*exp(-a*(x1(i))^2)*exp(-b*(z1(i))^2);
e2(i) = A*exp(-a*(x1(i))^2)*exp(-b*(z2(i))^2);
end
figure(k);
plot(x1,z1);
P(k) = sum(e1)*0.02- sum(e2)*0.02;
end
То есть сначала я нахожу такие координаты z, для которых
Распределение энергии у меня получается, оно гауссово с заданной амплитудой (интегральным значением)
А вот что касается самого давления от времени.. Получается построить его зависимость от каждой из координат (вполне правдоподобно), когда просто выкидываю дельта-функцию.
Как и можно ли работать с этой самой дельта-функцией в матлабе, не знаю..
Один из неработающих вариантов кода:
Код:
for k=1:1:101
for i=1:1:101
R = ((x-x1).^2+(z-z1(i)).^2).^(0.5);
d = dirac(t(k)-R./c);
Ev(i) = sum(exp(-b*z1(i)^2).*exp(-a*x1.^2)*0.1.*d./R);
I(i) = sum(exp(-b*z1(i)^2).*exp(-a*x1.^2)*0.1./R);
end
% t(i) = R./c;
P(k) = al*A^2*sum(I*4)/Cp/4/3.14;
Ev(k) = al*A^2*sum(Ev*4)/Cp/4/3.14;
% t = 0:0.01:t;
end
end
Посоветуете что-нибудь? Хотя бы идею...
