How can I correct this code?
Mostrar comentarios más antiguos
% Data Initialization
A = [ -0.1880
0.0100
-0.0079
-0.0273
0.1587
-0.4069
0.8943
-1.7404
3.0826
-5.1443
8.1501
-12.4334
18.3324
-26.3551
36.9611
-50.9219
68.8245
-91.8411
120.6744
-157.1298
201.9740
-258.0539
326.0616
-410.6317
512.0566
-638.1943
788.2017
-976.1147
1198.2528
-1481.4827
1815.1080
-2255.8411
2774.5508
-3510.8018
4379.5371
-5837.7578
7571.4097
-12936.2158
19555.4707
-10974.3643];
B = [ 11665.0576
-4.7293
13.9157
-26.5743
37.1578
-42.1387
44.1123
-45.5787
44.7030
-39.4265
30.5014
-20.7728
12.5115
-6.7371
3.2495
-1.4121
0.5494
-0.1905
0.0572
-0.0141
0.0023
0.0001
-0.0003
0.0002
-0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000];
C = [ -22.4193
-0.0217
0.0677
-0.1272
0.1340
-0.0548
-0.0433
0.1017
-0.1076
0.0619
-0.0047
-0.0195
0.0135
-0.0030
-0.0010
0.0009
-0.0002
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000];
D = [ -0.0871
-0.0006
0.0002
0.0004
-0.0012
0.0016
-0.0018
0.0017
-0.0015
0.0013
-0.0010
0.0008
-0.0006
0.0004
-0.0003
0.0002
-0.0001
0.0001
-0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000];
E = [ -0.1880
0.0100
-0.0079
-0.0273
0.1587
-0.4069
0.8943
-1.7404
3.0826
-5.1443
8.1501
-12.4334
18.3324
-26.3551
36.9611
-50.9219
68.8245
-91.8411
120.6744
-157.1298
201.9740
-258.0539
326.0616
-410.6317
512.0566
-638.1943
788.2017
-976.1147
1198.2528
-1481.4827
1815.1080
-2255.8411
2774.5508
-3510.8018
4379.5371
-5837.7578
7571.4097
-12936.2158
19555.4707
-10974.3643];
F = [ 11665.0576
-4.7293
13.9157
-26.5743
37.1578
-42.1387
44.1123
-45.5787
44.7030
-39.4265
30.5014
-20.7728
12.5115
-6.7371
3.2495
-1.4121
0.5494
-0.1905
0.0572
-0.0141
0.0023
0.0001
-0.0003
0.0002
-0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000];
% Parameters
a = 1; % Radius
L = 0.1; % Length
dd = 6; % Separation distance
alpha = 10;
etta = 0.01;
ettap = 0.1; % Additional parameters
zalm = alpha * complex(1, 0) / sqrt(2.d0);
xi = 1.0 / etta;
xi1 = 1.0 / ettap;
zk1 = sqrt((xi + sqrt(xi^2 - 4.d0 * xi * zalm^2)) / 2.0);
zk2 = sqrt((xi - sqrt(xi^2 - 4.d0 * xi * zalm^2)) / 2.0);
c = -a / L;
b = a / L;
m = a * 50; % Number of intervals
% Create Mesh Grid
[x, y] = meshgrid(c + dd : (b - c) / m : b, c : (b - c) / m : b);
% Apply Boundary Conditions
[I, J] = find(sqrt(x.^2 + y.^2) < (a - 0.01));
if ~isempty(I)
x(I, J) = NaN; % Use NaN to avoid affecting the visualization
y(I, J) = NaN;
end
% Compute Values
r = sqrt(x.^2 + y.^2);
t = atan2(y, x);
r2 = sqrt(r.^2 + dd.^2 - 2 .* r .* dd .* cos(t));
zet = (r.^2 - r2.^2 - dd.^2) ./ (2 .* r2 .* dd);
% Calculate Psi1
psi1 = zeros(size(x));
for i = 2:5
Ai = A(i-1);
Bi = B(i-1);
Ci = C(i-1);
Di = D(i-1);
Ei = E(i-1);
Fi = F(i-1);
psi1 = psi1 + (Ai .* r.^(-i + 1) + r.^(1 / 2) .* besselk(i - 1 / 2, r .* zk1) .* Bi + ...
r.^(1 / 2) .* besselk(i - 1 / 2, r .* zk2) .* Ci) .* ...
gegenbauerC(i, -1 / 2, cos(t)) + ...
(Di .* r2.^(-i + 1) + r2.^(1 / 2) .* besselk(i - 1 / 2, r2 .* zk1) .* Ei + ...
r2.^(1 / 2) .* besselk(i - 1 / 2, r2 .* zk2) .* Fi) .* ...
gegenbauerC(i, -1 / 2, zet);
end
% Plot Contours with Different Colors
figure;
contourf(y, x, psi1, 20); % Create filled contours with automatic color mapping
colorbar; % Add colorbar to check values
axis equal;
title('$\delta=1.5,\frac{\beta\mu}{a_1}=1.0,\alpha=1.0,\ell=0.1$', ...
'Interpreter', 'latex', 'FontSize', 12, 'FontName', 'Times New Roman', 'FontWeight', 'Normal');
% Plot Spheres
hold on;
% Sphere 1 (Vertical)
t3 = linspace(0, 2 * pi, 1000);
rr2 = 2.5;
x2 = rr2 * cos(t3); % Keep x and y aligned for vertical spheres
y2 = rr2 * sin(t3);
plot(y2, x2, '-k', 'LineWidth', 1.1);
fill(y2, x2, [0.5 0.5 0.5]); % Fill with gray color
% Sphere 2 (Vertical)
t2 = linspace(0, 2 * pi, 1000);
h = dd;
rr = 2.5;
x1 = rr * cos(t2) + h; % Keep x and y aligned for vertical spheres
y1 = rr * sin(t2);
plot(y1, x1, '-k', 'LineWidth', 1.1);
fill(y1, x1, [0.5 0.5 0.5]); % Fill with gray color
% Axis Formatting
axis off;
7 comentarios
Torsten
el 19 de Mzo. de 2025
psi1 is complex-valued - you can't make a contour plot of complex-valued functions.
What you can plot is abs(psi1), real(psi1) and imag(psi1).
But you should first check whether psi1 is computed correctly in your code.
Shreen El-Sapa
el 20 de Mzo. de 2025
Editada: Cris LaPierre
el 20 de Mzo. de 2025
Fangjun Jiang
el 20 de Mzo. de 2025
Editada: Fangjun Jiang
el 20 de Mzo. de 2025
See above. This line was missing. You can format the text to be code and run it right away to check.
dd = 6; % Separation distance
Cris LaPierre
el 20 de Mzo. de 2025
Your updated code is missing the variable dd. Once I add that (copied from the code you originally posted), it doesn't error out, but there still is no contour plot. You code now results in psi1 being an array of NaNs.
Shreen El-Sapa
el 20 de Mzo. de 2025
Editada: Walter Roberson
el 20 de Mzo. de 2025
Shreen El-Sapa
el 20 de Mzo. de 2025
Editada: Walter Roberson
el 21 de Mzo. de 2025
Torsten
el 21 de Mzo. de 2025
Try
[DH1,h1]=contour(y,x,real(psi1));
colorbar
and you will see that the range of your values for psi1 goes up until 1e154.
So a variation of the contour colors is restricted to a very small region in the graphics.
Respuestas (0)
Categorías
Más información sobre Line Plots en Centro de ayuda y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!

