دوشنبه, ۱ اسفند ۱۴۰۱، ۰۴:۰۰ ب.ظ
رشد مستطیل در یک مسیر سینوسی روی یک صفحه فضایی
شما می توانید با تغییر معادله صفحه زیر در سط 4،نمودار سه بعدی رشد یافته روی مسیر را روی آن صفحه مشاهده کنین.
f(p,q)=2*p +2*q;
مثلا می توانین معادله صفحه در سط 4 را به فرم 2+ f(p,q)=2*p +0.5*q; تغییر دهین و نتیجه را مشاهده کنین.
حتی می تونین، به جای صفحه معادله یک رویه را وارد کنین، مثلا سطر 4 را با معادله رویه f(p,q)=2*p +.5*q+2*p*q; جایگزین کنین. و نتیجه شکل زیر را مشاهده کنین.
کد اجرایی:
clc;
clear;
close all;
syms p q X Y Z yd(X) yl(X) f(p,q) ft(p,q) T xp(T) yp(T) zp(T)
f(p,q)=2*p +2*q;
[p1,q1]=meshgrid( 0:.1:1 ,-1:.1:1 );
surf(p1,q1,double(f(p1,q1)),FaceColor='g');alpha(.2);hold on;
t=linspace(0,1,100);
x=t;
y=sin(2*pi*x);yprim=2*pi*cos(2*pi*x);
plot3(x,y,f(x,y),'-r');hold on
n=length(t);
t1=linspace(0,1,200);
r=.08;
x1=r*sin(2*pi*t1);
y1=r*cos(2*pi*t1);
xM=[];
yM=[];
zM=[];
for s=1:length(t) -1
%plot equLine for slop and its 90 slop lines
%%
% yprims =double((y(s+1)- y(s))/(x(s+1)- x(s)));
yd(X)=yprim(s)*(X-x(s))+y(s);
yl(X)=(-1./yprim(s))*(X-x(s))+y(s);
tr=[x(s),x(s)+.00005 ];
d1=diff( double(tr));
d2=diff( double(yd(tr)));
d3=diff( double(f(tr,yd(tr))));
rd=[d1,d2,d3];
xp(T)=x(s)+rd(1)*T;
yp(T)=y(s)+rd(2)*T;
zp(T)=double(f(x(s),y(s)))+rd(3)*T;
tr1=double(solve((xp-x(s))^2+(yp-y(s))^2-r^2,T));
% plot3(xp(tr1),yp(tr1) ,zp(tr1),'b-' );hold on
%%
d1=diff( double(tr));
d2=diff( double(yl(tr)));
d3=diff( double(f(tr,yl(tr))));
r2=[d1,d2,d3];
% r2=[1,double(yl(x(s)+1))-double(yl(x(s))),double(f(x(s)+1,double(yl(x(s)+1))))-double(f(x(s),double(yl(x(s)))))];
xl(T)=x(s)+r2(1)*T;
yl(T)=y(s)+r2(2)*T;
zl(T)=double(f(x(s),y(s)))+r2(3)*T;
tr2=(double(solve((xl-x(s))^2+(yl-y(s))^2-r^2,T)));
% plot3(double(xl(tr2)),double(yl(tr2)) ,double(zl(tr2)),'c-' );hold on
df=(double((yprim(s))))/abs(double((yprim(s))));
if df==-1
tr2=tr2(end:df:1);
end
%%
% plot3(x1+x(s),y1+y(s) ,double(f(x1+x(s),y1+y(s))),'b-' );hold on
%%
r3=double([subs(gradient(f ),[p,q],[x(s),y(s)])',-1]);
xr=double(xl(tr2));
yr=double(yl(tr2));
zr=double(zl(tr2));
xr0(T)=xr(1)+r3(1)*T;
yr0(T)=yr(1)+r3(2)*T;
zr0(T)=zr(1)+r3(3)*T;
xr1(T)=xr(2)+r3(1)*T;
yr1(T)=yr(2)+r3(2)*T;
zr1(T)=zr(2)+r3(3)*T;
tr3=[-.02,.02];
% plot3(double([xr0(tr3),xr1([tr3(2),tr3(1)]),xr0(tr3(1))]),double([yr0(tr3),yr1([tr3(2),tr3(1)]),yr0(tr3(1))]),double([zr0(tr3),zr1([tr3(2),tr3(1)]),zr0(tr3(1))]));hold on
% plot3(xr0(tr3),yr0(tr3),zr0(tr3),'*')
xM=[xM;double([xr0(tr3),xr1([tr3(2),tr3(1)]),xr0(tr3(1))])];
yM=[yM;double([yr0(tr3),yr1([tr3(2),tr3(1)]),yr0(tr3(1))])];
zM=[zM;double([zr0(tr3),zr1([tr3(2),tr3(1)]),zr0(tr3(1))])];
% peripendecular plane
% peripendecular passpoints
rpX=diff(double(xl(tr2))');
rpY=diff(double(yl(tr2))') ;
rpZ=diff(double(zl(tr2)'));
B=[rpX,rpY,rpZ] ;
C=double([subs(gradient(f ),[p,q],[x(s),y(s)])',-1]);
C=cross(B,C);
refPlane=C(1)*(X-x(s ))+C(2)*(Y-y(s))+C(3)*(Z-f(x(s ),y(s))) ;
verPlaneCoff='XYZ';
[~,ind]=find(C);
ind=ind(end);
in=find(~ismember([1,2,3],ind));
ft(p, q)=subs(refPlane-C(ind)*sym(verPlaneCoff(ind)),[sym(verPlaneCoff(in(1))),sym(verPlaneCoff(in(2)))],[p,q]);
ft=-(ft/C(ind));
plot3(double([xr0(tr3),xr1([tr3(2),tr3(1)]),xr0(tr3(1))]),double([yr0(tr3),yr1([tr3(2),tr3(1)]),yr0(tr3(1))]),double([zr0(tr3),zr1([tr3(2),tr3(1)]),zr0(tr3(1))]),'r');hold on
% surf(double([xr0(tr3);xr1(tr3)]),double([yr0(tr3);yr1(tr3)]),double([zr0(tr3);zr1(tr3)]),FaceColor='g');hold on
% surf(double([xr0(tr3);xr1(tr3)]),double([yr0(tr3);yr1(tr3)]),double(ft(double([xr0(tr3);xr1(tr3)]),double([yr0(tr3);yr1(tr3)]))),FaceColor='b');hold on
plot3(double([xr0(tr3);xr1(tr3)])',double([yr0(tr3);yr1(tr3)])',double(ft(double([xr0(tr3);xr1(tr3)]),double([yr0(tr3);yr1(tr3)])))','b');hold on
plot3(double([xr0(tr3);xr1(tr3)]),double([yr0(tr3);yr1(tr3)]),double(ft(double([xr0(tr3);xr1(tr3)]),double([yr0(tr3);yr1(tr3)]))),'b');hold on
view(3)
xlabel('X');
ylabel('Y');
zlabel('Z');
% axis equal
pause(.0001)
end
axis equal
plot3(xM ,yM ,zM )
h=surf(xM ,yM ,zM ,'FaceColor','g');alpha(h,.1)
