[MATLAB] Modeling gravity using polar coordinates Argh

[ellipsis]

I've been using MATLAB (ode45) to simulate the mechanics of a rocket under the forces of gravity, drag, and internal thrust.y I've recently refactored my simulation to include 2d space, orientation of the rocket, etc. (So I can try to make it orbit, finding optimal ascent profiles, etc)

Because I want to retain the simplicity of my 1d model, I've attempted to use polar coordinates for defining the rockets position, velocity, and acceleration. I've implemented it, but when I run the simulation and plot height vs theta... I get a disturbing picture.

Gravity does not affect the orbital velocity in any way. It always accelerates the height of the rocket downward. And yet for it to orbit, it needs to have a cyclic, non-zero height. Something is very wrong here, and I can't tell what it is.

The line defining gravitational acceleration is as follows. I define it as a vector, with the first component being orbital acceleration (how much theta's rate of change changes due to gravity), and the second component being height's acceleration. (height is y, here)

grav = [0; -S/(R+y)^2];

Since S (the standard gravitation parameter) is positive, the acceleration of height will always be negative. The height of the rocket must inevitably go to zero. My intuition of gravity no longer makes sense.

This is bothering me on a visceral level. I can't tell what I'm doing wrong. I now suspect the Earth will fall into Sol in the next few hours.

 

[jedishrfu]

It would help if you posted the whole MATLAB script as some folks have MATLAB and could run it to see what's happening.

 

[ellipsis]

All right, I will post the whole thing. It's large though. There's four files, p_step.m, simul.m, eventy.m, and designer.m

p_step.m is the function ode45 calls, and is basically a system of nonlinear differential equations.
simul.m is the main script which calls ode45 and plots the data.
eventy.m defines events (out_of_fuel, and ground_collision)
designer.m specifies rocket properties (thrust, efficiency, etc)

Spoiler

p_step.m:

function dx = p_step(~, Z)

  global T II IF;
  global TH;
  

  S = 3.5316E12;
  R = 6E5;
  D = 9.784758844E-4;
  H = 5E3;
  G = 9.82;
  
  x = Z(1);
  y = Z(2);
  TH = y*3.14/200000;
  
  vx = Z(3);
  vy = Z(4);
  th = Z(5);
  m = Z(6);
  
  atm = exp(-y/H);
  
  grav = [0; -S/(R+y)^2];
  drag = -D*atm*[vx; vy].^2;
  thrust = T/m*[sin(TH); cos(TH)];
  
  a = thrust + drag + grav;
  ax = a(1);
  ay = a(2);
  
  dm = T/(((IF-II)*atm-IF)*G);
  dth = 0;
  
  dx = [vx; vy; ax; ay; dth; dm];
  
end

simul.m:

  global T II IF ME MF;

  type = 8;

  T = engine_thrusts(type);
  II = engine_atm_isps(type);
  IF = engine_vac_isps(type);
  ME = engine_masses(type);
  
  G = 9.82;
  R = 1.5;
  
  MI = T/(G*R);
  MF = (MI+7*ME)/8;
  
  
  opt = odeset('Events',@eventy);
  
  [t,y,~,ey,ei] = ode45(@p_step,[0 1000],[0, 0, 0, 0, 0, MI], opt); 
  plot(y(:,1), y(:,2));
  
  if (numel(ei)~=0)
  switch ei(1)
  case 1
  max_dv = ey(3);
  case 2
  max_dv = -1;
  end
  else
  fprintf('Flight timed out\n');
  max_dv = 0;
  end
  
  fprintf('max-velocity:\t%.2f\n',max_dv);

eventy.m:

function [zero, terminal, negative] = eventy(~,x)

  global MF;
  
  zero = [x(6) - MF; sign(x(2))+1];
  terminal = [1; 1];
  negative = [0; 0];
  
end

designer.m:

global engine_names engine_thrusts engine_masses engine_atm_isps engine_vac_isps;

engine_names = ['LV-1 Liquid Fuel Engine'... % 1
  'Rockomax 48-7S'... % 2
  'LV-909 Liquid Fuel Engine'... % 3
  'LV-N Atomic Rocket Motor'... % 4
  'R.A.P.I.E.R. Engine'... % 5
  'Toroidal Aerospike Rocket'... % 6
  'LV-T45 Liquid Fuel Engine'... % 7
  'LV-T30 Liquid Fuel Engine'... % 8
  'Rockomax "Poodle" Liquid Engine'... % 9
  'Rockomax "Skipper" Liquid Engine'... % 10
  'Rockomax "Mainsail" Liquid Engine'... % 11
  'Kerbodyne KR-2L Advanced Engine'... % 12
  'S3 KS-25x4 Engine Cluster']; % 13

engine_thrusts = [4 30 50 60 175 175 200 215 220 650 1500 2500 3200];

engine_masses = [.03 .1 .5 2.25 1.2 1.5 1.5 1.25 2 3 6 6.5 9.75];

engine_atm_isps = [220 300 300 220 320 388 320 320 270 320 320 280 320];

engine_vac_isps = [290 350 390 800 360 390 370 370 390 370 360 380 360];

fuel_names = ['Kerbodyne S3-14400 Tank'...
  'Kerbodyne S3-7200 Tank'...
  'Rockomax Jumbo-64 Fuel Tank'...
  'Kerbodyne S3-3600 Tank'...
  'Rockomax X200-32 Fuel Tank'...
  'Rockomax X200-16 Fuel Tank'...
  'FL-T800 Fuel Tank'...
  'Rockomax X200-8 Fuel Tank'...
  'FL-T400 Fuel Tank'...
  'FL-T200 Fuel Tank'...
  'FL-T100 Fuel Tank'...
  'ROUND-8 Toroidal Fuel Tank'...
  'Oscar-B Fuel Tank'];

fuel_full_mass = [82 41 36 20.5 18 9 4.5 4.5 2.25 1.125 0.5625 0.136 0.078675];

fuel_empty_mass = [10 5 4 2.5 2 1 0.5 0.5 0.25 0.125 0.0625 0.025 0.015];

I don't think the problem is in the MATLAB code, though. I think it's a conceptual problem I'm having. If acceleration due to gravity is perpendicular (has no affect on) orbital velocity, why do orbiting objects get faster near the gravitational body?

 

[jedishrfu]

Thats large? Medium for programming is 10K lines and large is 100K lines...

 

[ellipsis]

Thats large? Medium for programming is 10K lines and large is 100K lines...

My program is plenty large! D:<

Can you help me understand orbits jedish, or are you satisfied with insulting the length and girth of my program? ;)