How can I appear the streamline around the spheres?

11 次查看（过去 30 天）
Shreen El-Sapa
Shreen El-Sapa 2024-4-22
编辑： Voss 2024-10-29，20:12
Ran in:
clc
A=[ -0.2
0.1
-0.4
1.0
-2.1
3.9
-6.9
11.6
-18.7
29.2
-44.0
64.9
-93.5
132.6
-184.2
253.6
-343.1
462.3
-613.4
815.3
-1068.2
1413.3
-1842.6
2462.0
-3232.1
4517.3
-6127.8
10894.9
-17024.9
9869.9];
B=[ 279.2
34.9
-96.3
177.4
-248.7
283.6
-271.0
223.3
-160.5
102.6
-58.7
30.5
-14.4
6.3
-2.5
0.9
-0.3
0.1
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0];
C=[ -4.6
-0.3
0.3
-0.2
0.1
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0];
AA=[ -1.4
-0.1
0.4
-1.1
3.0
-7.5
17.7
-40.2
88.1
-188.2
390.9
-797.9
1594.7
-3150.4
6115.4
-11791.2
22393.6
-42435.9
79327.3
-148756.4
275335.2
-515215.3
950981.9
-1800517.1
3351792.8
-6647748.5
12804165.0
-32341874.0
71833152.0
-59216272.0];
BB=[ 15518.5
-406.6
1241.4
-2229.4
3366.8
-4306.2
4732.0
-4611.8
4004.3
-3152.0
2255.2
-1486.4
902.0
-510.4
268.4
-132.9
61.6
-27.2
11.3
-4.5
1.7
-0.6
0.2
-0.1
0.0
0.0
0.0
0.0
0.0
0.0];
CC=[ -15.6
0.3
-0.4
0.3
-0.2
0.1
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0];
%%%%%%%%%%%%%%%%%%%%
a = 1 ; %RADIUS
L=.22;
etta=0.01; %0.02;0.05;0.1d0
u2=1; delta=1.5; ettap=.01; alpha=2.d0; %C=a1/a2=0.1
xi=1./sqrt(etta);xi1=1./sqrt(ettap);
alpha1=sqrt((xi.^2+sqrt(xi.^4-4.*xi.^2.*alpha.^2))./sqrt(2));
alpha2=sqrt((xi.^2-sqrt(xi.^4-4.*xi.^2.*alpha.^2))./sqrt(2));
dd=5;
c =-a/L;
b =a/L;
m =a*40; % 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));
if ~isempty(I);
x(I,J) = 0; y(I,J) = 0;
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);
%for i1=1:length(x);
% for k1=1:length(x);
% if sqrt(x(i1,k1).^2+y(i1,k1).^2)>1./L;
% r(i1,k1)=0;r2(i1,k1)=0;
% end
% end
%end
warning off
psi1=0;
for i=2:7
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.^(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,psi1,25,'-k','LineWidth',1.1); %,psi2,'--k',psi2,':k'
%[DH1,h1]=contour(x,y,psi1);
%p1=contour(x,y,psi1,[0.3 0.3],'k','LineWidth',1.1); %,'ShowText','on'
%p2=contour(x,y,psi1,[0.4 0.4],'r','LineWidth',1.1);
%p3=contour(x,y,psi1,[0.5 0.5],'g','LineWidth',1.1);
%p4=contour(x,y,psi1,[0.6 0.6],'b','LineWidth',1.1);
%p5=contour(x,y,psi1,[0.7 0.7],'c','LineWidth',1.1);
%p6=contour(x,y,psi1,[0.8 0.8],'m','LineWidth',1.1);
%p7=contour(x,y,psi1,[0.9 0.9],'y','LineWidth',1.1);
p1=contour(x,y,psi1,[-0.001 0.001],'k','LineWidth',1.1); %,'ShowText','on'
p2=contour(x,y,psi1,[-0.005 .005],'r','LineWidth',1.1);
%p3=contour(x,y,psi1,[0.1 0.1],'g','LineWidth',1.1);
%p4=contour(x,y,psi1,[0.4 0.4],'b','LineWidth',1.1);
%p5=contour(x,y,psi1,[0.6 0.6],'c','LineWidth',1.1);
%p6=contour(x,y,psi1,[0.8 0.8],'m','LineWidth',1.1);
%%%%%%%%%%%%%%% $\frac{\textstyle a_1+a_2}{\textstyle h}=6.0,\;
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(x2,y2,'-k'),'LineWidth',1.1);
fill(x2,y2,'w')
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(x1,y1,'-k'),'LineWidth',1.1);
fill(x1,y1,'w')
%axis square;
axis('equal')
axis on
%xticklabels([])
%yticklabels([])
%legend('0.01','0.05','0.1','0.4','0.6','0.8','Location','northwest')
%title('$\frac{\beta_1}{a_1\mu}=\frac{a_1\beta_2}{\mu}=1.0,\;R_{H}=1.0,\;\frac{a_2}{a_1}=2.0$','Interpreter','latex','FontSize',12,'FontName','Times New Roman','FontWeight','Normal')
%title('$(a)\;\; R_{H}=1.0,\;\frac{\kappa}{\mu}=4.0$','Interpreter','latex','FontSize',12,'FontName','Times New Roman','FontWeight','Normal')
%%%%%%%%%%%%%%%%%%%%
view([90 90])

请先登录，再回答此问题。

回答（1 个）

Gautam
Gautam 2024-10-23，5:06
Hello Shreen,
You can plot the streamlines by using a series of “contour” functions the way you've been attempting to.
hold on
p11=contour(y,x,psi1,[1 1],'k','LineWidth',1.1);
p21=contour(y,x,psi1,[3 3],'r','LineWidth',1.1);
p31=contour(y,x,psi1,[4 4],'g','LineWidth',1.1);
p41=contour(y,x,psi1,[5 5],'b','LineWidth',1.1);
p51=contour(y,x,psi1,[10 10],'c','LineWidth',1.1);
p61=contour(y,x,psi1,[50 50],'m','LineWidth',1.1);
p71=contour(y,x,psi1,[100 100],'y','LineWidth',1.1);
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,'w')
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,'w')
axis('equal')
axis off
hold off
  1 个评论
shreen elsapa
shreen elsapa 2024-10-29，20:05
编辑：Voss 2024-10-29，20:12
Ran in:
@Gautam
Can you tell me about this code ?
clc
A =[ -4.56964
0.11715
-0.49082
1.14190
-2.15372
3.63497
-5.67271
8.40030
-11.89370
16.34167
-21.79783
28.55245
-36.60938
46.44017
-57.94571
71.92785
-88.07967
107.83102
-130.44882
158.64197
-190.75777
232.40297
-279.74454
346.19617
-421.83228
549.04956
-694.61816
1157.37988
-1701.72839
928.16284];
B=[ 3.88442
-0.03321
0.02748
-0.00907
0.00189
-0.00029
0.00003
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000];
AA=[ -4.49021
-0.11744
0.56476
-1.16020
2.19876
-3.70012
5.77588
-8.55282
12.10988
-16.63875
22.19414
-29.07158
37.27501
-47.28453
58.99926
-73.23563
89.68111
-109.79158
132.82062
-161.52635
194.22607
-236.62846
284.83081
-352.49063
429.50192
-559.03229
707.24756
-1178.42310
1732.66882
-945.03845];
BB=[ 3.82074
0.03326
-0.03161
0.00921
-0.00193
0.00029
-0.00004
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000];
% Parameters
a = 1 ; %RADIUS
L=.1;
alm1=.1; alpha=1; %u2=1;a2=1;beta1=1;beta2=1;delta=3; RH=alpha magnetic; alm1=alpha frequency
alpha1=((alpha^2-alm1^2 * complex(1, 0))^(0.5));
dd = 6;
c =-a/L;
b =a/L;
m =a*200; % NUMBER OF INTERVALS
%[x,y]=meshgrid((c+dd:(b-c)/m:b),(c:(b-c)/m:b)');
[x,y]=meshgrid((c+dd:(b-c)/m:b),(0:(b-c)/m:b)');
[I, J]=find(sqrt(x.^2+y.^2)<(a-0.1));
if ~isempty(I)
x(I,J) = 0; y(I,J) = 0;
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:7
Ai=A(i-1);Bi=B(i-1);AAi=AA(i-1);BBi=BB(i-1);
psi1=psi1+(Ai.*r.^(-i+1)+r.^(1./2).*besselk(i-1./2,r.*alpha1).*Bi ).*gegenbauerC(i,-1./2, cos(t))+(AAi.*r2.^(-i+1)+r2.^(1./2).*besselk(i-1./2,r2.*alpha1).*BBi).*gegenbauerC(i,-1./2,zet);
end
hold on
[DH1,h1]=contour(x,y,psi1,25,'-k','LineWidth',1.1); %,psi2,'--k',psi2,':k'
%[DH1,h1]=contour(x,y,psi1);
%p1=contour(x,y,psi1,[0.3 0.3],'k','LineWidth',1.1); %,'ShowText','on'
%p2=contour(x,y,psi1,[0.4 0.4],'r','LineWidth',1.1);
%p3=contour(x,y,psi1,[0.5 0.5],'g','LineWidth',1.1);
%p4=contour(x,y,psi1,[0.6 0.6],'b','LineWidth',1.1);
%p5=contour(x,y,psi1,[0.7 0.7],'c','LineWidth',1.1);
%p6=contour(x,y,psi1,[0.8 0.8],'m','LineWidth',1.1);
%p7=contour(x,y,psi1,[0.9 0.9],'y','LineWidth',1.1);
%%%%%%%%%%%%%%%
hold on
t3 = linspace(0,pi,1000);
h2=0;
k2=0;
rr2=2;
x2 = rr2*cos(t3)+h2;
y2 = rr2*sin(t3)+k2;
set(plot(x2,y2,'-k'),'LineWidth',1.1);
fill(x2,y2,'w')
hold on
t2 = linspace(0,pi,1000);
h=dd;
k=0;
rr=1;
x1 = rr*cos(t2)+h;
y1 = rr*sin(t2)+k;
set(plot(x1,y1,'-k'),'LineWidth',1.1);
fill(x1,y1,'w')
%axis square;
axis('equal')
box on
axis on
xticklabels([])
yticklabels([])
legend('0.01','0.05','0.1','0.4','0.6','0.8','Location','northwest')
Warning: Ignoring extra legend entries.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Vector Fields 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by