How can I plot the complete two circles vertical not horizontal ?
Show older comments
clc
A =[ -1.
0.
0.
0.
0.
0.
0.
0.
0.
-1.
1.
-1.
1.
-1.
2.
-2.
3.
-3.
4.
-5.
5.
-7.
8.
-10.
12.
-16.
20.
-34.
53.
-30.];
B=[ 3262.
131.
-375.
563.
-639.
602.
-486.
345.
-218.
124.
-64.
31.
-13.
5.
-2.
1.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.];
C=[ 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=[ 8.
0.
1.
-1.
3.
-8.
21.
-52.
126.
-307.
738.
-1771.
4215.
-10047.
23743.
-56327.
132493.
-313806.
736630.
-1749066.
4111518.
-9852368.
23316548.
-57140296.
137506208.
-357384160.
896199040.
-3046175232.
9340706816.
-10404635648.];
BB=[ -76625208.
858156.
-3341452.
4741591.
-7006134.
8310705.
-9026788.
8857093.
-7988619.
6701862.
-5230164.
3847242.
-2655485.
1743048.
-1080089.
641116.
-360810.
195865.
-101116.
50743.
-24261.
11394.
-5098.
2281.
-969.
430.
-179.
97.
-47.
8.];
CC=[ 29.
0.
1.
-1.
1.
-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.];
a = 1 ; %RADIUS
L=.1;
akm=2;gamma=0.3;arh=10; %beta1=beta2=1,a1=1,a2=2,arh=10,delta=0.5,u2=-1
alphaa=sqrt(((2+akm).*akm./(gamma.*(2+akm))).^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=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);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,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.01 0.01],'k','LineWidth',1.1); %,'ShowText','on'
p2=contour(x,y,psi1,[0.05 .05],'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);
%clabel(DH1,h1,'FontSize',10,'Color','red')
%%%%%%%%%%%%%%% $\frac{\textstyle a_1+a_2}{\textstyle h}=6.0,\;
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
%set(gca,'XTick',[], 'YTick', [])
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')
%%%%%%%%%%%%%%%%%%%%
2 Comments
Shreen El-Sapa
on 22 Apr 2024
Mathieu NOE
on 23 Apr 2024
hello again
you have to change the range of t2 and t3 from 0 / pi to 0 / 2*pi (btw they are identical so you could simply use one variable and not two) , like :
t3 = linspace(0,pi,1000); => t3 = linspace(0,2*pi,1000);
Accepted Answer
More Answers (0)
Categories
Find more on Mathematics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!

