plot stream over two spheres

Question

0 votes

A  =[       -4.7107
    0.0012
    -0.0056
    0.0132
    -0.0253
    0.0435
    -0.0689
    0.1031
    -0.1473
    0.2040
    -0.2737
    0.3607
    -0.4647
    0.5927
    -0.7425
    0.9265
    -1.1387
    1.4014
    -1.7014
    2.0810
    -2.5114
    3.0805
    -3.7224
    4.6475
    -5.6872
    7.5039
    -9.5388
    16.4146
    -25.4535
    14.3236];
B=[      -3.3794
    0.0005
    -0.0009
    0.0006
    -0.0003
    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];
C=[        6.8417
    -0.0007
    0.0007
    -0.0003
    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];
D=[      -4.7100
    -0.0012
    0.0058
    -0.0132
    0.0253
    -0.0435
    0.0689
    -0.1032
    0.1473
    -0.2040
    0.2737
    -0.3608
    0.4648
    -0.5928
    0.7426
    -0.9266
    1.1388
    -1.4016
    1.7016
    -2.0813
    2.5118
    -3.0810
    3.7230
    -4.6482
    5.6881
    -7.5050
    9.5403
    -16.4171
    25.4574
    -14.3258];
E=[       -3.3789
    -0.0005
    0.0009
    -0.0006
    0.0003
    -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];
F=[       6.8407
    0.0007
    -0.0008
    0.0003
    -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];
a = 1 ;   %RADIUS
L=.1;
dd=4;
kappa=1;gam=0.3;arh=1;  %a2=1;u2=1;beta1=beta2=1
al=kappa.*(2+kappa)./(gam.*(1+kappa));
alpha1=real(((al.^2+arh.^2)./2+((al.^2+arh.^2).^2-(2.*kappa.*arh.^2./gam).^(1./2))./2).^(1./2));
alpha2=real(((al.^2+arh.^2)./2-((al.^2+arh.^2).^2-(2.*kappa.*arh.^2./gam).^(1./2))./2).^(1./2));
c =-a/L;
b =a/L;
m =a*100; % 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-.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);
%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
qr1=0;
for i=2:7
    Ai=A(i-1);Bi=B(i-1);Ci=C(i-1);Di=D(i-1);Ei=E(i-1);Fi=F(i-1);
    qr1=qr1-(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))-(Di.*r2.^(-i-1)+r2.^(-3./2).*besselk(i-1./2,r2.*alpha1).*Ei+r2.^(-3./2).*besselk(i-1./2,r2.*alpha2).*Fi).*legendreP(i-1,zet);
end
hold on
[DH1,h1]=contour(x,y,qr1,3,'-k');
%axis square;
title('$(a)$ $\ell=0.1,\;\alpha=1.0$','Interpreter','latex','FontSize',10,'FontName','Times New Roman','FontWeight','Normal')
%%%%%%%%%%%%%%%%  $\frac{\textstyle a_1+a_2}{\textstyle h}=6.0,\;
hold on
t3 = linspace(0,2*pi,1000);
h2=0;
k2=0;
rr2=1;
x2 = rr2*cos(t3)+h2;                   
y2 = rr2*sin(t3)+k2;
set(plot(x2,y2,'-k'),'LineWidth',1.1);
fill(x2,y2,'w')
%axis square;
hold on
t2 = linspace(0,2*pi,1000);
h=dd;
k=0;
rr=2;
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 off

1 Comment
Show -1 older comments Hide -1 older comments

Adam Danz on 20 Nov 2021

Shreen El-Sapa's answer moved here

I can not get the streamlines

Sign in to comment.

Sign in to answer this question.

Follow Question

Answer 1

Cris LaPierre on 20 Nov 2021

1 vote

qr1 is all NaNs. Assuming the countours are supposed to be your streamlines, you should check your equation. I'm not sure your for loop is doing what you intended. At the least, there is an issue with your calculation.

5 Comments
Show 3 older comments Hide 3 older comments

Shreen El-Sapa on 20 Nov 2021

Edited: Cris LaPierre on 20 Nov 2021

Open in MATLAB Online

edited code to make it more readable

A=[-4.7107 0.0012 -0.0056 0.0132 -0.0253 0.0435 -0.0689 0.1031 -0.1473 0.2040 -0.2737 0.3607 -0.4647 0.5927 -0.7425 0.9265 -1.1387 1.4014 -1.7014 2.0810 -2.5114 3.0805 -3.7224 4.6475 -5.6872 7.5039 -9.5388 16.4146 -25.4535 14.3236];

B=[-3.3794 0.0005 -0.0009 0.0006 -0.0003 0.0001 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=[6.8417 -0.0007 0.0007 -0.0003 0.0001 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=[-4.7100 -0.0012 0.0058 -0.0132 0.0253 -0.0435 0.0689 -0.1032 0.1473 -0.2040 0.2737 -0.3608 0.4648 -0.5928 0.7426 -0.9266 1.1388 -1.4016 1.7016 -2.0813 2.5118 -3.0810 3.7230 -4.6482 5.6881 -7.5050 9.5403 -16.4171 25.4574 -14.3258];

BB=[-3.3789 -0.0005 0.0009 -0.0006 0.0003 -0.0001 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];

CC=[6.8407 0.0007 -0.0008 0.0003 -0.0001 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=.13;

akm=1;gamma=0.3;arh=0.1;

alphaa=sqrt(((2+akm).*(akm./gamma)).^2+arh.^2);

betaa=(2.*akm.*arh.^2./gamma).^(0.25);

alpha1=sqrt((alphaa.^2+sqrt(alphaa.^4-4.*betaa.^4))./2);

alpha2=sqrt((alphaa.^2-sqrt(alphaa.^4-4.*betaa.^4))./2);

dd=5;

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-.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:3;

Ai=A(i-1);

Bi=B(i-1);

Ci=C(i-1);

AAi=AA(i-1);

BBi=BB(i-1);

CCi=C(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);

end

[DH1,h1]=contour(x,y,psi1,20,'-k');

title('$(a)$ $\ell=0.1,\;\alpha=1.0$','Interpreter','latex','FontSize',10,'FontName','Times New Roman','FontWeight','Normal')

hold on

t3 = linspace(0,2*pi,1000);

h2=0;

k2=0;

rr2=1;

x2 = rr2*cos(t3)+h2;

y2 = rr2*sin(t3)+k2;

set(plot(x2,y2,'-k'),'LineWidth',1.1);

fill(x2,y2,'w')

t2 = linspace(0,2*pi,1000);

h=dd;

k=0;

rr=2;

x1 = rr*cos(t2)+h;

y1 = rr*sin(t2)+k;

set(plot(x1,y1,'-k'),'LineWidth',1.1);

fill(x1,y1,'w')

axis off

hold off

Cris LaPierre on 20 Nov 2021

Personally, I use the streamlines function to create streamlines, not contour. However, even with streamlines, you will need to calculate and input the vector field components u and v.

Shreen El-Sapa on 20 Nov 2021

Capture.GIF

I used the equations below. Can I used it to plot the streamlines? as you say?

Sign in to comment.

plot stream over two spheres

1 Comment
Show -1 older comments Hide -1 older comments

Accepted Answer

5 Comments
Show 3 older comments Hide 3 older comments

More Answers (0)

Categories

Tags

Community Treasure Hunt

plot stream over two spheres

1 Comment Show -1 older comments Hide -1 older comments

Accepted Answer

5 Comments Show 3 older comments Hide 3 older comments

More Answers (0)

Categories

Tags

See Also

Community Treasure Hunt

1 Comment
Show -1 older comments Hide -1 older comments

5 Comments
Show 3 older comments Hide 3 older comments