Asked 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);

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)

Sign in to comment.

Answer by Eliraz Nahum
on 30 Sep 2018

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

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## madhan ravi (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/421469-problem-figuring-out-2-different-solutions#comment_616040

## Eliraz Nahum (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/421469-problem-figuring-out-2-different-solutions#comment_616043

Sign in to comment.