Планетная физика шаткая

Question

Итак, на моем уроке физики мы работаем с законом Кеплера и другими подобными уравнениями. Я подумал, что пойму это лучше, если смогу сделать из этого пример. и у меня это работает ... Сорта. Предполагается, что движение планет или лун в данном случае является эллипсом, и они просто не находятся близко.

Я использую фактические массы, расстояния и скорости для планетоидов, поэтому я не уверен, что происходит не так

let G = 6.67408*Math.pow(10,-11)
let scl = 4000
let mars;
let phobos;
let deimos;

function setup() {
	createCanvas(600, 600)
	mars = new Planet(0, 0, 0, 0, 6.4171*pow(10,24))
	phobos = new Planet(-9234420, 0, 0, -2138 , 1.0659*pow(10,16))
	deimos = new Planet(23463200, 0, 0, 1351.3, 1.4762*pow(10,15))
	background(0)
}

function draw() {
	translate(width/2-mars.pos.x/scl, height/2-mars.pos.y/scl)
	background(0, 1)
	for(let i=0;i<10;i++){
		mars.show()
		phobos.show()
		deimos.show()
		mars.update()
		phobos.update()
		deimos.update()
		grav(phobos, mars)
		grav(mars, deimos)
		grav(phobos, deimos)
		
		if(deimos.pos.dist(mars.pos)/scl>(width/2)){
			scl = deimos.pos.dist(phobos.pos)/(0.99*width/2);
			background(0)
		}
		
	}
}

function force(p1, p2){
	d = p1.pos.dist(p2.pos)
	f = G*(p1.mass*p2.mass)/(d*d)
	return f
}

function grav(p1,p2){
	d1 = p2.pos.copy().sub(p1.pos).normalize()
	d2 = p1.pos.copy().sub(p2.pos).normalize()

	f = force(p1,p2)
	
	p1.acc.add(d1.mult(force/p1.mass))
	p2.acc.add(d2.mult(force/p2.mass))
}

function Planet(x, y, vx, vy, mass){
	this.mass = mass;
	this.pos = createVector(x, y)
	this.vel = createVector(vx, vy)
	this.acc = createVector(0, 0)
	
	this.show = function(){
		stroke(255)
		strokeWeight(width/100)
		point(this.pos.x/scl, this.pos.y/scl)
	}
	
	this.update = function(){
		this.pos = this.pos.add(this.vel);
		this.vel = this.vel.add(this.acc)
		this.acc = this.acc.mult(0)
	}
}

<script src="https://cdnjs.cloudflare.com/ajax/libs/p5.js/0.8.0/p5.js"></script>

как я сказал, можно подумать, что у них будут эллиптические орбиты, но это не так. так что я просто в растерянности.

Answer 1

Чтобы физика работала в этом коде, нам нужно внести пару исправлений.

Функция grav должна вычислять силу и использовать ее для обновления ускорения, однако есть ошибка, которая вызываетфункция силы, которая будет использоваться вместо вычисленной силы, а исходный код умножается на NaN вместо силы / массы, поскольку функция force не является числом.

function grav(p1,p2){
  d1 = p2.pos.copy().sub(p1.pos).normalize()
  d2 = p1.pos.copy().sub(p2.pos).normalize()

  f = force(p1,p2) // we need to store the calculated force
  p1.acc.add(d1.mult(f/(p1.mass))) // and use it here
  p2.acc.add(d2.mult(f/(p2.mass))) // and here
}

Код также должен вызываться update после каждого вызова функции grav:

grav(phobos, mars);
mars.update();
phobos.update();

grav(mars, deimos);
mars.update();
deimos.update();

grav(phobos, deimos);
phobos.update();
deimos.update();

let G = 6.67408*Math.pow(10,-11)
let scl = 4000
let mars;
let phobos;
let deimos;

function setup() {
  createCanvas(600, 600);
  mars = new Planet(0, 0, 0, 0, 6.4171*pow(10,24));
  phobos = new Planet(-9234420, 0, 0, -2138 , 1.0659*pow(10,16));
  deimos = new Planet(23463200, 0, 0, 1351.3, 1.4762*pow(10,15));
  background(0);
}

function draw() {
  translate(width/2-mars.pos.x/scl, height/2-mars.pos.y/scl);
  background(0, 1);
  for(let i=0;i<10;i++){
    mars.show();
    phobos.show();
    deimos.show();

    grav(phobos, mars);
    mars.update();
    phobos.update();

    grav(mars, deimos);
    mars.update();
    deimos.update();
    
    grav(phobos, deimos);
    phobos.update();
    deimos.update();
    
    if(deimos.pos.dist(mars.pos)/scl>(width/2)){
      scl = deimos.pos.dist(phobos.pos)/(0.99*width/2);
      background(0)
   }
  }
}

function force(p1, p2){
	d = p1.pos.dist(p2.pos)
	f = G*(p1.mass*p2.mass)/(d*d)
	return f
}

function grav(p1,p2){
  let d1 = p2.pos.copy().sub(p1.pos).normalize()
  let d2 = p1.pos.copy().sub(p2.pos).normalize()

  let f = force(p1,p2)
  p1.acc.add(d1.mult(f/(p1.mass)))
  p2.acc.add(d2.mult(f/(p2.mass)))
}

function Planet(x, y, vx, vy, mass){
  this.mass = mass;
  this.pos = createVector(x, y)
  this.vel = createVector(vx, vy)
  this.acc = createVector(0, 0)

  this.show = function(){
    stroke(255)
    strokeWeight(width/100)
    point(this.pos.x/scl, this.pos.y/scl)
  }

  this.update = function(){
    this.pos = this.pos.add(this.vel);
    this.vel = this.vel.add(this.acc)
    this.acc = this.acc.mult(0)
  }
}

    <script src="https://cdnjs.cloudflare.com/ajax/libs/p5.js/0.7.2/p5.min.js"></script>

С этими поправками мы получаем примерно эллиптические орбиты, однако после нескольких оборотов спутники отлетают от Марса.Это неудивительно, поскольку все начальные условия и значение G являются приблизительными.