problem figuring out 2 different solutions

1 view (last 30 days)
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 all
close all
clc
t1=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;
end
TvsW=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);
end
plot(x,y,'.b')
xlim([0 (10*10^6)])
ylim([0 (3*10^7)])
figure
plot(v_time,x,'.r')
hold on
plot(v_time,y,'.g')
figure
plot(v_time,xdot,'.r')
hold on
plot(v_time,ydot,'.g')
figure(1)
hold on
for counter=0:0.001*pi:2*pi
plot(r*cos(counter),r*sin(counter),'.m')
hold on
end
hold 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 off
figure(2);
hist(length_arc,20);
  2 Comments
madhan ravi
madhan ravi on 30 Sep 2018
"I attached the definition of the physical problem as a photo" you didn't
Eliraz Nahum
Eliraz Nahum on 30 Sep 2018
thanks for informing me. I attached now.

Sign in to comment.

Accepted Answer

Bruno Luong
Bruno Luong on 30 Sep 2018
Edited: Bruno Luong on 30 Sep 2018
                 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

The dt coded above is timestep, which is wrong if should be (t-s1) or (t-s0)

More Answers (1)

Eliraz Nahum
Eliraz Nahum on 30 Sep 2018
thank you very much!!! you have no idea how you helped us.

Products


Release

R2018a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!