hello everyone,
my friend and I have to hand out our HW today, and can't understand why after going over and over again on our codes – we can't get identical plots. the main difference is with x values. I attached the definition of the physical problem as a photo.
In addition, I am attaching the first code:
clear allclose allclct1=32.2;w0=6169;w11=2695;t2=138.2;w12=2082;w21=397;w22=100;T1=26950;T2=3973;x0=0;y0=2.0926*(10^7);xdot0=34.7;ydot0=196.9;g=32.15;GM=1.4077*(10^16); t_delta=1;counter=1;time=0;v_time=[];x=[];y=[];xdot=[];ydot=[];x2dot=[];y2dot=[];x(counter)=x0;y(counter)=y0;xdot(counter)=xdot0;ydot(counter)=ydot0;division1=T1/w0;x2dot(counter)=(-GM*(x(1)/((x(1)^2 + y(1)^2)^1.5)))+(g*division1)*(xdot(1)/((xdot(1)^2 + ydot(1)^2)^0.5));y2dot(counter)=(-GM*(y(1)/((x(1)^2 + y(1)^2)^1.5)))+(g*division1)*(ydot(1)/((xdot(1)^2 + ydot(1)^2)^0.5));r=sqrt(x0^2+y0^2);size=r;while (size>=r)time=time+t_delta;if time<t1 W=w0+(((w11-w0)/t1)*time); T=T1;elseif (time>=t1)&&(time<t2) W=w12+(((w21-w12)/(t2-t1))*(time-t1)); T=T2;else W=w22; T=0;endTvsW=T/W;counter=counter+1;x(counter)=x(counter-1)+xdot(counter-1)*t_delta;y(counter)=y(counter-1)+ydot(counter-1)*t_delta;xdot(counter)=xdot(counter-1)+x2dot(counter-1)*t_delta;ydot(counter)=ydot(counter-1)+y2dot(counter-1)*t_delta;x2dot(counter)=(-GM*(x(counter)/((x(counter)^2 + y(counter)^2)^1.5)))+(g*TvsW)*(xdot(counter)/((xdot(counter)^2 + ydot(counter)^2)^0.5));y2dot(counter)=(-GM*(y(counter)/((x(counter)^2 + y(counter)^2)^1.5)))+(g*TvsW)*(ydot(counter)/((xdot(counter)^2 + ydot(counter)^2)^0.5));v_time(counter)=time;size=sqrt(x(counter)^2+y(counter)^2);endplot(x,y,'.b')xlim([0 (10*10^6)])ylim([0 (3*10^7)])figureplot(v_time,x,'.r')hold onplot(v_time,y,'.g')figureplot(v_time,xdot,'.r')hold onplot(v_time,ydot,'.g')figure(1)hold onfor counter=0:0.001*pi:2*pi plot(r*cos(counter),r*sin(counter),'.m') hold onendhold off
and here is the second code:
close all;clear all;clc;%Basic data
GM = 1.4077 * 10^16;v0 = 200;gama = deg2rad(80);g = 32.17;dt=1;Time=100;%Starting conditions
x0 = 0;y0 = 2.0926 * 10^7;xdot0 = v0 * cos(gama);ydot0 = v0 * sin(gama);for n=1:Time T = 26950; w = 6169; t=0; s1(1,n) = 30.59 + rand * 2 * 1.61; s2(1,n) = 131.29 + rand * 2* 6.91; x_val(1,1) = x0; x_dot_val(1,1) = xdot0; %=x1dot
y_val(1,1) = y0; % =R world
y_dot_val(1,1) = ydot0; %=y1dot
while sqrt(x_val ^ 2 + y_val ^ 2) >= y0 v = ( x_dot_val ^ 2 + y_dot_val ^ 2 ) ^0.5; aT = ( g * T ) / w; x_2dot_val = ( -GM * x_val) / ((x_val ^ 2 + y_val ^ 2)^1.5) + (aT * x_dot_val) / v; y_2dot_val = ( -GM * y_val) / ((x_val ^ 2 + y_val ^ 2)^1.5) + (aT * y_dot_val) / v; x_val = x_val + dt * x_dot_val; x_dot_val = x_dot_val + dt * x_2dot_val; y_val = y_val + dt * y_dot_val; y_dot_val = y_dot_val + dt * y_2dot_val; if t<=s1(1,n) w = ((2695-6169)/(32.2-0))*dt + 6169; elseif t>s2(1,n) w = 100; T = 0; else w = ((397 - 2082) / (138.2 - 32.2)) * dt + 2082; T = 3973; end t=t+dt; x_vector(1,t) = x_val; y_vector(1,t) = y_val; x_dot_vector(1,t) = x_val; y_dot_vector(1,t) = y_val; x_2dot_vector(1,t) = x_val; y_2dot_vector(1,t) = y_val; end length_arc(1,n) = atan (x_val/y_val) * y0; %חישוב הקשת במעגל
figure(1); plot(x_vector,y_vector); hold on;end% g = viscircles(x,y,0,0,y0);
% plot(g);
th = 0:pi/50:2*pi;xunit = y0 * cos(th) ;yunit = y0 * sin(th);h = plot(xunit, yunit);hold offfigure(2);hist(length_arc,20);
Best Answer
The dt coded above is timestep, which is wrong if should be (t-s1) or (t-s0)