I want to plot the half only in the attached figure. How can I do that?
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%subplot(1,2,1)
A =[ -16.
0.
0.
0.
0.
1.
-2.
7.
-20.
54.
-145.
384.
-989.
2525.
-6332.
15811.
-38887.
95760.
-232600.
569166.
-1374036.
3371349.
-8148474.
20348640.
-49802872.
131423232.
-334120288.
1149872256.
-3565926912.
4013070592.];
B=[ -350738.
26.
-101.
194.
-291.
368.
-403.
395.
-347.
280.
-207.
142.
-91.
54.
-31.
16.
-8.
4.
-2.
1.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.];
C=[ 35.
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=[ -4.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
-1.
1.
-1.
1.
-1.
1.
-2.
2.
-3.
3.
-4.
5.
-9.
14.
-8.];
BB=[ -95.
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.];
CC=[ 5.
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=.1;
akm=4;gamma=0.3;arh=1; %beta1=beta2=1,a1=1,a2=2,arh=1,delta=0.1,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=5;
c =-a/L;
b =a/L;
m =a*60; % 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.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,15,'-k'); %,psi2,'--k',psi2,':k'
p11=contour(x,y,psi1,[1 1],'k','LineWidth',1.1);
p21=contour(x,y,psi1,[3 3],'r','LineWidth',1.1);
p31=contour(x,y,psi1,[4 4],'g','LineWidth',1.1);
p41=contour(x,y,psi1,[5 5],'b','LineWidth',1.1);
p51=contour(x,y,psi1,[10 10],'c','LineWidth',1.1);
p61=contour(x,y,psi1,[50 50],'m','LineWidth',1.1);
p71=contour(x,y,psi1,[100 100],'y','LineWidth',1.1);
%axis square;
%beta1=beta2=1,a1=1,a2=2,arh=1,delta=0.1,u2=1
%akm=4;
%title('$\frac{\beta_1}{a_1\mu}=\frac{a_1\beta_2}{\mu}=1.0,\;R_{H}=1.0,\;\frac{a_2}{a_1}=0.5$','Interpreter','latex','FontSize',12,'FontName','Times New Roman','FontWeight','Normal')
%title('$(b)$','Interpreter','latex','FontSize',12,'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=2;
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=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 off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 commentaires
Réponses (1)
Gautam
le 23 Oct 2024
Hello @Shreen El-Sapa
You can make the following changes to the contour plot to plot only the top half
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);
And change the variables "t2" and "t3" to:
t2 = linspace(0,pi,1000);
t3 = linspace(0,2*pi,1000);

0 commentaires
Voir également
Catégories
En savoir plus sur Line Plots dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!