How can I correct this code?

5 views (last 30 days)
Shreen El-Sapa
Shreen El-Sapa on 19 Mar 2025
Commented: Torsten on 21 Mar 2025
Ran in:
% 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
Error using contourf (line 55)
Input arguments must be real. Use the function REAL to get the real part of the inputs.
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 Comments
Show 5 older comments
Shreen El-Sapa
Shreen El-Sapa on 20 Mar 2025
Edited: Walter Roberson on 21 Mar 2025
Ran in:
clc
A =[ 0.4383
-0.4080
1.4130
-3.7497
8.7052
-18.3631
35.6095
-64.5347
110.4802
-180.8033
284.0856
-432.3704
638.2991
-920.6697
1296.9564
-1795.9568
2440.4854
-3274.5806
4326.2554
-5663.6147
7318.2793
-9397.8857
11932.9268
-15098.8994
18913.5957
-23675.0488
29361.4824
-36506.2852
44985.0234
-55821.1367
68631.0156
-85581.2891
105597.8906
-134030.5938
167689.9219
-224158.2500
291518.9688
-499382.0625
756808.3125
-425743.2188];
B=[ 784.3210
-41.3373
140.2780
-272.2328
368.2668
-379.0945
310.1994
-209.5842
118.3485
-54.6368
18.1556
-1.5761
-3.6361
3.8433
-2.6156
1.4482
-0.6952
0.2997
-0.1176
0.0427
-0.0144
0.0045
-0.0013
0.0004
-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];
C=[ 3.9761
-0.1226
0.3397
-0.5729
0.5937
-0.3249
0.0574
0.0449
-0.0379
0.0126
-0.0014
-0.0005
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
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000];
AA=[ -0.1536
0.0249
-0.0437
0.0587
-0.0688
0.0729
-0.0710
0.0645
-0.0553
0.0453
-0.0356
0.0271
-0.0200
0.0144
-0.0101
0.0070
-0.0048
0.0032
-0.0021
0.0014
-0.0009
0.0006
-0.0004
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];
BB=[ 6.3555
-0.2784
-0.1038
1.4075
-2.6351
2.8560
-2.3294
1.5033
-0.8018
0.3637
-0.1427
0.0493
-0.0152
0.0042
-0.0011
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
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000];
CC=[ -1.2158
0.1330
-0.2657
0.1976
-0.0445
-0.0175
0.0124
-0.0028
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
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000E+00
0.0000E+00
0.0000E+00];
a = 1 ; %RADIUS
L=.11;
etta=0.01; ettap=0.01; alpha=4;delta=1.5;
lambda= complex(0,alpha)/sqrt(2);
xi=sqrt(1./etta); xi1=sqrt(1./ettap);
alpha1=sqrt((xi+sqrt(xi^2-4.*xi.*lambda.^2))./2);
alpha2=sqrt((xi-sqrt(xi^2-4.*xi.*lambda.^2))./2);
dd=6; %h
c =-a/L;
b =a/L;
m =a*50; % NUMBER OF INTERVALS
[x,y]=meshgrid((c+dd:(b-c)/m:b),(c:(b-c)/m:b)');
[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
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);
warning on
psi1=0;
for i=2:25
Ai=A(i-1);Bi=B(i-1);Ci=C(i-1);AAi=AA(i-1);BBi=BB(i-1);CCi=CC(i-1);
%psi1=-psi1-(Ai.*r.^(-i-1)+r.^(-3./2).*besselk(i-1./2,r.*alpha1).*Bi+r.^(-3./2).*besselk(i-1./2,r.*alpha2).*Ci).*legendreP(i-1,cos(t))-(AAi.*r2.^(-i-1)+r2.^(-3./2).*besselk(i-1./2,r2.*alpha1).*BBi+r2.^(1./2).*besselk(i-1./2,r2.*alpha2).*CCi).*legendreP(i-1,zet);
psi1=psi1+(Ai.*r.^(-i+1)+r.^(1./2).*besselk(i-1./2,r.*alpha1).*Bi+r.^(1./2).*besselk(i-1./2,r.*alpha2).*Ci).*gegenbauerC(i,-1./2, cos(t))+(AAi.*r2.^(-i+1)+r2.^(1./2).*besselk(i-1./2,r2.*alpha1).*BBi+r2.^(1./2).*besselk(i-1./2,r2.*alpha2).*CCi).*gegenbauerC(i,-1./2,zet);
end
hold on
%[DH1,h1]=contour(x,y,real(psi1),'--k','LineWidth',1.1,'ShowText','on'); %,psi2,'--k',psi2,':k'
[DH1,h1]=contour(y,x,real(psi1),'-k');
%p1=contour(y,x,real(psi1),[0.001 0.001],'-c','LineWidth',1.3,'ShowText','on'); %,'ShowText','on'
%p2=contour(y,x,real(psi1),[0.005 0.005],'-r','LineWidth',1.3);
%p3=contour(y,x,real(psi1),[0.05 0.05],'-b','LineWidth',1.3);
%p4=contour(y,x,real(psi1),[.1 .1],'-k','LineWidth',1.3);
%p5=contour(y,x,real(psi1),[.2 .2],'--c','LineWidth',1.3);
%p6=contour(y,x,real(psi1),[.4 .4],'--r','LineWidth',1.3);
%p7=contour(y,x,real(psi1),[.6 .6],'--b','LineWidth',1.3);
%p8=contour(y,x,real(psi1),[.8 .8],'--k','LineWidth',1.3);
%p9=contour(y,x,real(psi1),[1 1],'-.c','LineWidth',1.3);
%p10=contour(y,x,real(psi1),[2 2],'-.r','LineWidth',1.3);
%p11=contour(y,x,real(psi1),[3 3],'-.b','LineWidth',1.3);
%p12=contour(y,x,real(psi1),[4 4],'-.k','LineWidth',1.3);
%axis equal; % Maintain aspect ratio
%axis tight; % Adjust axis limits to fit the data
%%%%%%%%%%%%%%%
hold on
t3 = linspace(0,2*pi,1000);
h2=0;
k2=0;
rr2=2;
x2 = rr2*cos(t3)+h2;
y2 = rr2*sin(t3)+k2;
set(plot(y2,x2,'-k'),'LineWidth',1.1);
fill(y2,x2,[0.7 0.7 0.7])
hold on
t2 = linspace(0,2*pi,1000);
h=dd;
k=0;
rr=1;
x1 = rr*cos(t2)+h;
y1 = rr*sin(t2)+k;
set(plot(y1,x1,'-k'),'LineWidth',1.1);
fill(y1,x1,[0.7 0.7 0.7])
axis off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Torsten
Torsten on 21 Mar 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.

Sign in to comment.

Sign in to answer this question.

Answers (0)

Categories

Find more on Line Plots in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!