Re: MATLAB 5.3.1 (1999) and latest Octave math discrepancy ..

[Will Flannery]

I have no problem showing the code.   First I'll describe what it does.

The program simulates the Juno space probe ... the probe is launched from a
low earth satellite orbit and goes into a 2 year orbit around the sun, after
2 years it returns to earth, trailing it by just a bit, and earth's
gravitational field pulls it, increasing it's velocity so that it breaks its
2 year orbit and is 'slingshot' toward Jupiter.  You can see a simulation
here ...

http://www.youtube.com/watch?v=sYp5p2oL51g

The program computes the gravitational forces of the earth and the sun on
the probe, resolves them into x and y components, an projects them over 2
second subintervals, it does the same thing for the earth's orbit around the
sun.

I did the program a long time ago using MATLAB, and just downloaded
Octave3.6.2gcc4.6.2 a few days ago.  I cut and pasted the program into
Octave and it ran no problem.

The thing looks good until the probe returns to earth, and there the earth's
gravity takes over, and the differences in the graphs, unnoticeable to this
point, get magnified greatly by earth's gravity and the trajectories diverge
significantly.

The code - for some reason way back when, when I was trying to get this
thing to work, I broke the code into two scripts, which I have preserved
just because I don't want to change anything.  I ran it with MATLAB and got
the result shown.

Note - the program can't store every state value for graphing, as it
requires too much storate, so it stores the state every 2000 iterations. 
This produces a lot of unessential code.

The main program ... 
_______________________  Orbit11  ____________________________

clear


G=6.7e-11;
mSun = 1.9e30;   
rEarthOrbit = 150e9;
rSun=695500000;
mEarth = 5.9742e24;   
rEarth = 6.378e6;

dt=2;
N=3*365*60*60*24 / dt

M=2000;
I=N/M+1;
Ii=1;

XE=zeros(I,1); YE=zeros(I,1); %VXE=zeros(I,1); VYE=zeros(I,1);
XP=zeros(I,1); YP=zeros(I,1); %VYP=zeros(I,1); VYP=zeros(I,1); 
RPE=zeros(I,1);VP=zeros(I,1);
rPEmin=200000;
rPE(1)=rPEmin;

t1=0; T(1)=0;

% set launch parameters for a earth orbit

xE1=rEarthOrbit; XE(1)=xE1;         % starting x coordinate
vxE1=0; vXE(1)=vxE1;            % starting x velocity
yE1=0;  yE(1)=yE1;          %2*rEarthOrbit*pi/(365*60*60*24)
vyE1=28500;  vYE(1)=vyE1;               % starting y velocity

% set launch parameters for probe 
xP1=rEarthOrbit + rEarth + 200000;  XP(1)=xP1;        
vxP1=0;   vXP(1)=vxP1;       
yP1=0;  yP(1)=yP1;
LV=5265;  % 5260 gives slingshot with collision  

vyP1=28500 + LV;  vYP(1)=vyP1; vP(1)=vyP1;   % earth's velocity + launch
velocity

Orbit11b

length(XP)
I2=30400;
I1=30100;


plot(XP,YP,'r')                         % plot the probe's orbit

hold on
plot(XE,YE)                         % plot the earth's orbit

__________ Orbit11b __________________

for i=1:N
   % project the earth's orbit

   t2 = t1+ dt;
   xE2=xE1+vxE1*dt;        % project x position
   yE2=yE1+vyE1*dt;        % project y position
   
   % resolve the earth sun gravity vector
   rE = sqrt(xE1^2 + yE1^2); % compute rE
   g = -G*mSun/rE^2;             % compute gravity magnitude
   ctheta = xE1/rE;             % compute cos theta
   stheta = yE1/rE;             % compute sin theta
   
   % now we're ready to project velocity
   vxE2= vxE1+ g*ctheta*dt;  % project x velocity
   vyE2 = vyE1+ g*stheta*dt;  % project y velocity
   
   % project the probe's orbit

   xP2=xP1+vxP1*dt;        % project x position
   yP2=yP1+vyP1*dt;        % project y position
   
   % resolve the probe sun gravity vector
   rP = sqrt(xP1^2 + yP1^2); % compute rP
   g = -G*mSun/rP^2;             % compute gravity magnitude
   ctheta = xP1/rP;             % compute cos theta
   stheta = yP1/rP;             % compute sin theta
   
   % resolve the probe earth gravity vector
   xPE = xP1 - xE1;
   yPE = yP1 - yE1;
   rPE = sqrt((xPE^2 + yPE)^2); % compute rPE
   gE = -G*mEarth/rPE^2;             % compute gravity magnitude
   comega = xPE/rPE;             % compute cos theta
   somega = yPE/rPE;             % compute sin theta

   
   % now we're ready to project velocity
   vxP2 = vxP1+ (g*ctheta + gE*comega)*dt;  % project x velocity
   vyP2 = vyP1+ (g*stheta + gE*somega)*dt;  % project y velocity
      
   if mod(i,M)== 0        
      Ii=Ii+1;
 
      Ii*M*dt/(60*60*24)
      XP(Ii)= xP2;
      %VXP(Ii)= vxP2;
      YP(Ii)= yP2;
      %VYP(Ii)= vyP2;
      
      XE(Ii)= xE2;
      %VXE(Ii)= vxE2;
      YE(Ii)= yE2;
      %VYE(Ii)= vyE2;
      
      VP(Ii)=sqrt(vxP2^2+vyP2^2);
      RPE(Ii)=rPEmin;
      rPEmin=rPE;
      end
   
   if(rPE<rPEmin)rPEmin=rPE;
   end
   
  

   t1=t2;
   xP1=xP2;
   vxP1=vxP2;
   yP1=yP2;
   vyP1=vyP2;
   xE1=xE2;
   vxE1=vxE2;
   yE1=yE2;
   vyE1=vyE2;
  _____________________________________ 


   







--
View this message in context: 
http://octave.1599824.n4.nabble.com/MATLAB-5-3-1-1999-and-latest-Octave-math-discrepancy-tp4651123p4651130.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.