plot stream over two spheres

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

[Adam Danz]

Shreen El-Sapa's answer moved here

I can not get the streamlines

Answer 1

[Cris LaPierre]

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.

[Shreen El-Sapa]

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]

Ok. Is this what you wanted? If not, you're going to have to describe better what it is you expect.

[Shreen El-Sapa]

But this program plots streamlines directly without determining what is each value of streamlines over the contour? can you help me or anyone?

thanks

[Cris LaPierre]

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]

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