MATLAB Answers

0

problem figuring out 2 different solutions

Asked by Eliraz Nahum on 30 Sep 2018
Latest activity Answered by Eliraz Nahum on 30 Sep 2018
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

"I attached the definition of the physical problem as a photo" you didn't
thanks for informing me. I attached now.

Sign in to comment.

Products


Release

R2018a

2 Answers

Answer by Bruno Luong
on 30 Sep 2018
Edited by Bruno Luong
on 30 Sep 2018
 Accepted Answer

                 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)

  0 Comments

Sign in to comment.


Answer by Eliraz Nahum on 30 Sep 2018

thank you very much!!! you have no idea how you helped us.

  0 Comments

Sign in to comment.