Из закона всемирного тяготения известно, что на поверхности тел шарообразной формы ускорение свободного падения постоянно по модулю и направлено к центру шара. Для тел неправильной формы это правило, очевидно, не выполняется. В этой статье я покажу способ расчёта и визуализации ускорения свободного падения для таких тел. Расчёт будем производить на JavaScript, визуализировать — на WebGL с использованием библиотеки three.js.

В итоге получим следующее (красным цветом отмечены области с большим ускорением свободного падения, синим — с малым):



Гифка


Ось вращения на гифке условная. Она не совпадает с осью вращения кометы.

Для расчетов гравитационного потенциала планет используется следующая формула:



Но, к сожалению, при малых r этот ряд сходится медленно, а при очень малых r вообще расходится. Поэтому для расчета гравитационного потенциала на поверхности небесных тел эта формула не годится. В общем случае, придется провести интегрирование по всему объёму тела.
Трёхмерное тело можно представить в виде наборов треугольных граней, с заданными координатами вершин. Далее будем предполагать, что плотность тела постоянна (для кометы это примерно соответствует истине). Гравитационный потенциал в заданной точке будет равен сумме потенциалов всех тетраэдров, с основанием в этих гранях и с вершиной в заданной точке (мы ищем не гравитационный потенциал, а ускорение свободного падения, которое есть градиент потенциала, но рассуждения остаются теми же).

Для начала нам нужно найти ускорение, которое производит масса тетраэдра на точку, находящуюся в его вершине. Для этого нужно произвести интегрирование по всему объёму тетраэдра. К сожалению данный интеграл не берётся в элементарных функциях, поэтому придётся пойти на хитрость.



Сила, действующая на точку в вершине тетраэдра приблизительно в три раза больше, чем сила, вызванная гравитацией шара, помещенного в середину основания тетраэдра (масса шара при этом равна массе тетраэдра). Т. е. F1?3*F2. Это равенство лучше выполняется при малом угле раствора тетраэдра. Для оценки величины отклонения равенства от строгого, я сгенерировал несколько тысяч рандомных тетраэдров и вычислил для них эту величину.



На графике по оси абсцисс показано отношение периметра треугольного основания к расстоянию от вершины до центра основания. По оси ординат — величина ошибки. Как видно, величину ошибки можно представить квадратичной функцией. За редким исключением все точки находятся ниже параболы (исключения нам погоду особо не испортят).

Величину гравитационного ускорения будем вычислять с заданной относительной погрешностью. Если ошибка вычисления превышает эту погрешность, то разбиваем наш тетраэдр на четыре части по основанию и вычисляем ускорение для этих частей по отдельности, затем суммируем.



Ускорение мы находим с точностью до некоторого постоянного, для данного тела, коэффициента, зависящего от массы (или плотности) тела, гравитационной постоянной и некоторых других параметров. Этот коэффициент мы учтём позже.

Пришло время кода
function force_pyramide(p1, p2, p3, rel) {
	// начальная точка в (0, 0, 0)
	// rel - относительная погрешность
	var volume = (1/6) * triple_product(p1, p2, p3); // объём тетраэдра
	if (volume == 0)
		return new vector3(0, 0, 0);
	if (!rel)
		rel = 0.01;
	var p0 = middleVector(p1, p2, p3); // вектор центра основания
	var len = p0.length();
	var per = perimeter(p1, p2, p3); // периметр основания
	var tan_per = per / len;
	var error = 0.015 * sqr(tan_per); // та самая эмпирическая квадратичная функция
	if (error > rel) {
		var p12 = middleVector(p1, p2);
		var p13 = middleVector(p1, p3);
		var p23 = middleVector(p2, p3);
		return sumVector(
			force_pyramide(p1, p12, p13, rel),
			force_pyramide(p2, p23, p12, rel),
			force_pyramide(p3, p13, p23, rel),
			force_pyramide(p12, p23, p13, rel)
		);
	}
	var ratio = 3 * volume * Math.pow(len, -3);
	return new vector3(p0.x, p0.y, p0.z).multiplyScalar(ratio);
}


Вычисление ускорения от трёхмерного тела производим суммированием по всем тетраэдрам, образованным гранями тела и заданной точкой (о чём я писал ранее).

Скрытый текст
function force_object(p0, obj, rel) {
	var result = new vector3(0, 0, 0);
	for (var i = 0; i < obj.faces.length; i++) {
		p1 = subVectors(obj.vertices[obj.faces[i][0] - 1], p0);
		p2 = subVectors(obj.vertices[obj.faces[i][1] - 1], p0);
		p3 = subVectors(obj.vertices[obj.faces[i][2] - 1], p0);
		var f = force_pyramide(p1, p2, p3, rel);
		result = addVectors(result, f);
	}
	return result;
}


Вычисляем гравитационное и центробежное ускорения на комете (вращение кометы вызывает ускорение порядка 25% от гравитационного, поэтому пренебрегать им нельзя). Относительную погрешность я выставил в 0,01. Казалось бы это много, но по факту при вычислении погрешность получается раза в три меньше и при визуализации разница абсолютно не заметна (т. к. минимальная разница в цвете пикселей составляет 1/256?0,004). А если выставить погрешность меньше, то время расчёта увеличивается.

При погрешности, равной 0,01, расчёт выполняется 1-2 секунды, поэтому оформляем его через setInterval, во избежание подвисания браузера.

Код
var rel = 0.01;
// относительная погрешность
var scale = 1000;
// модель задана в километрах, переводим в метры
var grav_ratio = 6.672e-11 * 1.0e+13 / (sqr(scale) * volume);
// масса кометы 1е+13 кг
var omega2 = sqr(2 * Math.PI / (12.4 * 3600));
// период вращения кометы 12,4 часов

function computeGrav() {
	info.innerHTML = 'Расчёт: ' + (100 * item / object3d.vertices.length).toFixed() + '%';
	for (var i = 0; i < 50; i++) {
		var p0 = object3d.vertices[item];
		grav_force[item] = force_object(p0, object3d, rel).multiplyScalar(grav_ratio);
		// гравитационное ускорение
		circular_force[item] = new vector3(omega2 * p0.x * scale, omega2 * p0.y * scale, 0);
		// центробежное ускорение, ось вращения совпадает с осью z
		abs_grav_force[item] = grav_force[item].length();
		abs_circular_force[item] = circular_force[item].length();
		item++;
		if (item >= object3d.vertices.length) {
			console.timeEnd('gravity calculate');
			clearInterval(timerId);
			accelSelect();
			init();
			animate();
			break;
		}
	}
}

var item = 0;
console.time('gravity calculate');
var timerId = setInterval(computeGrav, 1);


Теперь стоит рассказать о формате работы с 3D телом. Тут всё просто. Это объект, состоящий из трёх массивов: вершины, нормали и треугольные грани.

var obj = {
	vertices: [],
	normals: [],
	faces: []
};

Массив вершин хранит объекты-векторы координат вершин. Массив нормалей — векторы нормалей вершин. Массив граней — массивы из трёх чисел, хранящие номера вершин.

Модель кометы Чурюмова-Герасименко я скачал с сайта ESA и упростил её в 3Ds Max. В исходной модели было несколько десятков тысяч вершин и граней, но для нашей задачи это слишком много, т. к. вычислительная сложность алгоритма зависит от произведения количества вершин на количество граней. Модель сохранена в формате obj.

В библиотеке three.js есть функция OBJLoader для загрузки этого формата, но при загрузке остается только информация о вершинах, а информация о гранях теряется, что не подходит для нашей задачи. Поэтому я её несколько модифицировал.

Код
function objLoad(url) {
	var result;
	var xhr = new XMLHttpRequest();
	xhr.open('GET', url, false);
	xhr.onreadystatechange = function(){
		if (xhr.readyState != 4) return;
		if (xhr.status != 200) {
			result = xhr.status + ': ' + xhr.statusText + ': ' + xhr.responseText;
		}
		else {
			result = xhr;
		}
	}
	xhr.send(null);
	return result;
}

function objParse(url) {
	var txt = objLoad(url).responseText;
	var lines = txt.split('\n');
	var result;
	
	// v float float float
	var vertex_pattern = /v( +[\d|\.|\+|\-|e|E]+)( +[\d|\.|\+|\-|e|E]+)( +[\d|\.|\+|\-|e|E]+)/;
	// f vertex vertex vertex ...
	var face_pattern1 = /f( +-?\d+)( +-?\d+)( +-?\d+)( +-?\d+)?/;
	// f vertex/uv vertex/uv vertex/uv ...
	var face_pattern2 = /f( +(-?\d+)\/(-?\d+))( +(-?\d+)\/(-?\d+))( +(-?\d+)\/(-?\d+))( +(-?\d+)\/(-?\d+))?/;
	// f vertex/uv/normal vertex/uv/normal vertex/uv/normal ...
	var face_pattern3 = /f( +(-?\d+)\/(-?\d+)\/(-?\d+))( +(-?\d+)\/(-?\d+)\/(-?\d+))( +(-?\d+)\/(-?\d+)\/(-?\d+))( +(-?\d+)\/(-?\d+)\/(-?\d+))?/;
	// f vertex//normal vertex//normal vertex//normal ...
	var face_pattern4 = /f( +(-?\d+)\/\/(-?\d+))( +(-?\d+)\/\/(-?\d+))( +(-?\d+)\/\/(-?\d+))( +(-?\d+)\/\/(-?\d+))?/;

	var obj = {
		vertices: [],
		normals: [],
		faces: []
	};

	for (var i = 0; i < lines.length; i++) {
		var line = lines[i].trim();
		if ((result = vertex_pattern.exec(line)) !== null) {
			obj.vertices.push(new vector3(
				parseFloat(result[1]),
				parseFloat(result[2]),
				parseFloat(result[3])
			));
		}else if ((result = face_pattern1.exec(line)) !== null) {
			obj.faces.push([
				parseInt(result[1]),
				parseInt(result[2]),
				parseInt(result[3])
			]);
			if (result[4])
				obj.faces.push([
					parseInt(result[1]),
					parseInt(result[3]),
					parseInt(result[4])
				]);
		}else if ((result = face_pattern2.exec(line)) !== null) {
			obj.faces.push([
				parseInt(result[2]),
				parseInt(result[5]),
				parseInt(result[8])
			]);
			if (result[11])
				obj.faces.push([
					parseInt(result[2]),
					parseInt(result[8]),
					parseInt(result[11])
				]);
		}else if ((result = face_pattern3.exec(line)) !== null) {
			obj.faces.push([
				parseInt(result[2]),
				parseInt(result[6]),
				parseInt(result[10])
			]);
			if (result[14])
				obj.faces.push([
					parseInt(result[2]),
					parseInt(result[10]),
					parseInt(result[14])
				]);
		}else if ((result = face_pattern4.exec(line)) !== null) {
			obj.faces.push([
				parseInt(result[2]),
				parseInt(result[5]),
				parseInt(result[8])
			]);
			if (result[11])
				obj.faces.push([
					parseInt(result[2]),
					parseInt(result[8]),
					parseInt(result[11])
				]);
		}
	}
	obj.normals = computeNormalizeNormals(obj);
	return obj;
}


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

Сцена создаётся просто.

scene = new THREE.Scene();

С камерой тоже никаких проблем.

var fieldWidth = 500; // размеры холста в пикселях
var fieldHeight = 500;
camera = new THREE.PerspectiveCamera(50, fieldWidth / fieldHeight, 0.01, 10000);
scene.add(camera);
cameraZ = 3 * boundingSphereRadius;
// boundingSphereRadius - радиус сферы, описанной вокруг нашего тела
camera.position.x = 0;
camera.position.y = 0;
camera.position.z = cameraZ;
camera.lookAt(new THREE.Vector3(0, 0, 0)); // камера смотрит в начало координат

Для создания камеры мы использовали функцию THREE.PerspectiveCamera(fov, aspect, near, far), где:

fov — высота поля зрения камеры в градусах;
aspect — отношение горизонтального угла зрения камеры к вертикальному;
near — расстояние до ближнего плана (то, что находится ближе, не будет рендериться);
far — расстояние до дальнего плана.

Ставим свет.

var ambientLight = new THREE.AmbientLight(0xffffff);
scene.add(ambientLight);

Создаём визуализатор.

renderer = new THREE.WebGLRenderer({antialias: true});
renderer.setClearColor(0xffffff); // цвет фона
renderer.setPixelRatio(window.devicePixelRatio);
renderer.setSize(fieldWidth, fieldHeight); // задаём размер холста
container.appendChild(renderer.domElement); // привязываемся к холсту

Всё вместе
var container;
var camera, scene, renderer;
var axis;
var mesh;
var boundingSphereRadius, cameraZ;
var lines = [];
var angleGeneral = 0;

function init() {
	container = document.getElementById('container');
	scene = new THREE.Scene();
	
	var fieldWidth = 500;
	var fieldHeight = 500;
	camera = new THREE.PerspectiveCamera(50, fieldWidth / fieldHeight, 0.01, 10000);
	scene.add(camera);
	var ambientLight = new THREE.AmbientLight(0xffffff);
	scene.add(ambientLight);

	loadModel();
	
	cameraZ = 3 * boundingSphereRadius;
	camera.position.x = 0;
	camera.position.y = 0;
	camera.position.z = cameraZ;
	camera.lookAt(new THREE.Vector3(0, 0, 0));
	
	axis = new THREE.Vector3(0.6, 0.8, 0); // 

	renderer = new THREE.WebGLRenderer({antialias: true});
	renderer.setClearColor(0xffffff);
	renderer.setPixelRatio(window.devicePixelRatio);
	renderer.setSize(fieldWidth, fieldHeight);
	container.appendChild(renderer.domElement);
}


Для работы с моделями в three.js удобно использовать класс BufferGeometry.

var geometry = new THREE.BufferGeometry();

Задаём материал.

var material = new THREE.MeshLambertMaterial( {
	side: THREE.FrontSide, // отображаться будет только передняя сторона грани
	vertexColors: THREE.VertexColors // грани окрашиваются градиентом по цветам их вершин
});

Загружаем координаты вершин.

var vertices = new Float32Array(9 * object3d.faces.length);
for (var i = 0; i < object3d.faces.length; i++) {
	for (var j = 0; j < 3; j++) {
		vertices[9*i + 3*j    ] = object3d.vertices[object3d.faces[i][j] - 1].x;
		vertices[9*i + 3*j + 1] = object3d.vertices[object3d.faces[i][j] - 1].y;
		vertices[9*i + 3*j + 2] = object3d.vertices[object3d.faces[i][j] - 1].z;
	}
}
geometry.addAttribute('position', new THREE.BufferAttribute(vertices, 3));

Вычисляем цвет вершин.

var colors = addColor();
geometry.addAttribute('color', new THREE.BufferAttribute(colors, 3));

function addColor() {
	var colorMin = [0.0, 0.6, 1.0]; // цвет точек с минимальным ускорением
	var colorMax = [1.0, 0.0, 0.0]; // то же с максимальным
	var colors = new Float32Array(9 * object3d.faces.length);
	for (var i = 0; i < object3d.faces.length; i++) {
		for (var j = 0; j < 3; j++) {
			var intensity = (abs_force[object3d.faces[i][j] - 1] - forceMin) / (forceMax - forceMin);
			colors[9*i + 3*j    ] = colorMin[0] + (colorMax[0] - colorMin[0]) * intensity;
			colors[9*i + 3*j + 1] = colorMin[1] + (colorMax[1] - colorMin[1]) * intensity;
			colors[9*i + 3*j + 2] = colorMin[2] + (colorMax[2] - colorMin[2]) * intensity;
		}
	}
	return colors;
}

Также нам понадобится отобразить направления ускорения свободного падения в каждой вершине.

Код
var lines;

function addLines() {
	lines = new Array();
	
	for (var i = 0; i < object3d.vertices.length; i++) {
		var color;
		var ratio;
		if (angle_force[i] < 90) {
			color = 0x000000; // Вектор ускорения направлен в сторону тела
			ratio = -0.2 * boundingSphereRadius / forceMax;
		}
		else {
			color = 0xdddddd; // Вектор ускорения направлен наружу
			ratio = 0.2 * boundingSphereRadius / forceMax;
		}
		
		var material = new THREE.LineBasicMaterial({
			color: color
		});
		
		var geometry = new THREE.Geometry();
		var point1 = object3d.vertices[i];
		var point2 = new THREE.Vector3();
		point2.copy(force[i]);
		point2.addVectors(point1, point2.multiplyScalar(ratio));
		geometry.vertices.push(
			new THREE.Vector3(point1.x, point1.y, point1.z),
			new THREE.Vector3(point2.x, point2.y, point2.z)
		);
	
		var line = new THREE.Line(geometry, material);
		rotateAroundWorldAxis(line, axis, angleGeneral);
		if (hair.checked) 
			scene.add(line);
		lines.push(line);
	}
}


С визуализацией закончили. Смотрим, что получилось.

Демо 1.

Как видим, ускорение свободного падения на комете от 40 до 80 тысяч раз меньше, чем на Земле. Максимальное отклонение вектора ускорения от нормали составляет порядка 60-70 градусов. Мелкие и крупные камни на этих участках вероятно не задерживаются и потихоньку скатываются в области, где угол не такой большой.

Можно поиграться с телами другой формы. Демо 2.

Видим, что для куба максимальное ускорение (в центре грани) составляет 1,002 ускорения для шара такой же массы, но на самом деле эта величина чуть меньше единицы (сыграло роль то, что мы считали с относительной погрешностью 0,01). Для куба, в отличие от тетраэдра, существуют точные формулы расчёта ускорений и точное значение для центра грани составляет (для куба со стороной, равной 1):



Для шара того же объёма:



Их отношение составляет 0.999376 и лишь слегка не дотягивает до единицы.

В заключение вопрос. Существуют ли тела, у которых хотя бы в одной точке отношение абсолютной величины гравитационного ускорения к ускорению на поверхности шара той же массы и объёма больше единицы? Если да, то для какого тела это отношение максимально? Во сколько раз это отношение больше единицы, в разы или на доли процентов?

Комментарии (38)


  1. p53
    23.11.2015 11:01

    2,6 %


    1. Milliard
      23.11.2015 12:07

      А форма тела какая?


      1. p53
        23.11.2015 17:58
        +1

        Вот такая.


        <zanuda>
        Если уж вспоминать баяны классику, то хотя бы с верной формулировкой. Нужно требовать постоянную плотность. Если попросить, чтобы лишь массы и объёмы совпадали, то планета, большая часть массы которой сосредоточена в точке на поверхности, сможет иметь сколь угодно большое g, и задача не будет иметь смысла.
        </zanuda>


        1. Milliard
          23.11.2015 18:07

          В теле поста писал про постоянную плотность, прошу прощения, что не продублировал в вопросе.
          Про книжку не знал.


  1. Aclz
    23.11.2015 13:19

    Гравитационный потенциал в заданной точке будет равен сумме потенциалов всех тетраэдров, с основанием в этих гранях и с вершиной в заданной точке

    А разве для этого не должно выполняться условие, что тетраэдр при этом обязательно должен полностью лежать в объеме тела? Для вогнутых тел, типа этой кометы, такое может в некоторых точках не выполняться.


    1. Milliard
      23.11.2015 13:34

      Для тех граней, которые обращены к заданной точке передней стороной, потенциал будет положительным, а которые задней — отрицательным.
      В качестве аналогии — подсчет площади плоской фигуры, заданной конечным числом вершин. Строим из любой точки векторы ко всем вершинам и ищем площади образованных треугольников через векторное произведение. Очевидно, что для ребер, вершины которых располагаются по часовой стрелке относительно исходной точки, знак векторного произведения будет один, которые против часовой — противоположный. В сумме будем иметь площадь. Аналогично вычисляется объём тела через смешанное произведение векторов.


      1. Aclz
        23.11.2015 15:45

        Точно. Весьма элегантно.


  1. khdavid
    23.11.2015 14:17
    +1

    Я не очень понимаю, почему нельзя было без сферических функций и размышлений про тетраеэдры?
    Просто берешь дробишь тело на очень маленькие кубики массой m. В каждой точке потенциал складывается от потенциалов этих кубов. От каджого куба потенциал примерно пропорционален m/r. В пределе, когда размер кубов стремится к нулю, вы получите правильно значение потенциала.


    1. Milliard
      23.11.2015 15:06

      Можно и на кубики. Но во-первых, задача разбиения тела на кубики не так тривиальна, как кажется на первый взгляд. Во-вторых, ближние кубики будут иметь намного большее влияние на общий потенциал, чем дальние. Тогда придётся либо всё тело разбивать на очень маленькие кубики (что значительно увеличивает время расчёта), либо разбивать только ближние кубики (но это разбиение придётся проводить для каждой точки, что опять же неблагоприятно влияет на время расчёта). к тому же кубики, примыкающие к границам тела будут не просто кубиками, а усечёнными кубиками и учёт их вклада в потенциал потребует много вычислительных ресурсов.


      1. khdavid
        23.11.2015 15:17

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

        2) Вклад ближних кубиков ( или тетраедров ) в потенциал пренебрежетельно мал по сравнению со вкладом дальних, это можно понять двумя способами
        a) Интеграл от функции 1/r по трехмерному пространству сходится в нуле, из этого все следует.
        b) Или же так. Берем маленькую сферу шириной dr и радиуса r. Потенциал от нее в центре этой сферы будет пропорционалень 1/r * ( r^2 * dr). Видим, чем ближе эта сфера к центру, тем меньше от нее вклад.


        1. p53
          25.11.2015 17:01

          Когда вы будете брать градиент потенциала, то внезапно ближайшие кубики начнут вносить вклад, и ещё огого какой.


          1. khdavid
            26.11.2015 11:35

            Еще раз, потенциал от кубика растет как 1/r, при уменьшении радиуса, но зато количество кубиков уменьшается квадратично: как r^2. Поэтому вклад ближних кубиков будет пропорционален r, то есть пренебрежительно мал по сравнению с вкладом дальних кубиков.
            Про градиент я не понял. Никаких проблем с ним не должно быть.


            1. Aclz
              26.11.2015 11:51

              Про градиент я не понял. Никаких проблем с ним не должно быть.
              Сам потенциал не интересен, интересна сила тяжести (которую можно получить через градиент потенциала). А он уже, в свою очередь, от расстояния зависит согласно принципу тех же обратных квадратов, что и рост количества кубиков. Т.е. вроде бы вклад ближних кубиков равен вкладу дальних. Но при этом дальние кубики находятся не в одной точке, а разбросаны влево-вправо, поэтому латеральная составляющая создаваемой ими силы взаимно гасит друг друга. Поэтому вклад дальних кубиков будет даже меньше, чем ближних.


            1. Milliard
              26.11.2015 12:55

              Накидал код. Кубик со стороной 1 разбиваем на 1000 кубиков (10 по каждой стороне) и считаем силу притяжения из вершины кубика. Это первый вариант. Во втором варианте самый ближний кубик разбиваем еще на 10 кубиков и считаем, что в этом случае получается. Как видим, относительная погрешность в первом случае 2,0%, во втором 0,17%.


              1. p53
                26.11.2015 14:17

                Скорее всего имелось ввиду что-то типа этого, интегрирование потенциала, а затем взятие градиента. Но при таком способе вычислиельная трудность растёт как куб разбиения, а ошибка уменьшается линейно. Так себе способ.


                1. khdavid
                  26.11.2015 14:47

                  Чтобы расчеты были точнее можно отдельно один раз очень точно посчить графитационное поле генерируемое одним кубиком.
                  А от кубической зависимости времени расчетов от мелкости разбиения имхо никуда не уйти.


                  1. p53
                    26.11.2015 16:39

                    Как же не уйти, если в статье интеграл берётся по площади = квадратичная зависимость.


                    1. khdavid
                      26.11.2015 17:04

                      В статье интегрируется по объему


                      1. p53
                        26.11.2015 18:31

                        Вы невнимательно читали статью. Фактически поверхность тела разбивается на малые площадки (треугольники), а само тело представляется как алгебраическая сумма конусов с малыми телесными углами (вытянутых тетраэдров). Тогда от каждого такого конуса вклад в силу тяжести составит:
                        image

                        Конечно, мы нормальную проекцию площадки вычисляем через объём конуса image, но это не отменяет того, что интегрируем мы по площади (хотя совсем корректно было бы сказать, что мы интегрируем по телесному углу, так как близкие треугольники нам приходится дополнительно разбивать, чтобы image).

                        Кстати, это сразу нам даёт идею для альтернативного алгоритма методом Монте-Карло: надо из точки пустить n лучей, и из расстояний от точки до пересечения лучей с телом можно будет посчитать итоговую силу.


                        1. Aclz
                          26.11.2015 20:31

                          Кстати, это сразу нам даёт идею для альтернативного алгоритма методом Монте-Карло
                          Вы предлагаете вместо прокладывания лучей на вершины бросать рандомные лучи? А в чем преимущество? Искать расстояние до объекта в простом случае придется перебором всех граней и проверкой пересечения луча с каждой, как итог, сложность решения будет расти весьма быстро с увеличением объектов (грубо: количество вершин х количество лучей х количество граней). Либо могу лишь предположить для облегчения поиска пересечения перевод координат вершин граней относительно каждой рассматриваемой точки в полярные координаты и сортировку по ним, что будет быстрее в плане поиска пересечения, но объем работ также будет б0льшим, чем однопроходный алгоритм (сложность: количество вершин х количество граней + тесселяция в зависимости от объекта)…


                          1. p53
                            27.11.2015 07:34

                            А в чем преимущество?

                            Я не утверждал, что для многогранника Монте-Карло будет чем-то лучше. Я написал «альтернативный алгоритм» =)


    1. Aclz
      23.11.2015 15:59

      Видимо, дело в том, что а) форма объекта уже дана в виде множества координат вершин треугольников б) метод расчета потетраэдрно работает.


  1. Karroplan
    23.11.2015 14:34

    IMHO, рассуждение "Трёхмерное тело можно представить в виде наборов треугольных граней, с заданными координатами вершин. Далее будем предполагать, что плотность тела постоянна (для кометы это примерно соответствует истине). Гравитационный потенциал в заданной точке будет равен сумме потенциалов всех тетраэдров, с основанием в этих гранях и с вершиной в заданной точке (мы ищем не гравитационный потенциал, а ускорение свободного падения, которое есть градиент потенциала, но рассуждения остаются теми же)" верно только для выпуклых тел.

    В случае с невыпуклым телом (например, комета Чурюмова-Герасименко) нельзя так просто брать и перебором обходить все грани, т.к. появятся некорректные тетраэдры.

    Соотвественно, расчет неверен. Особенно вот это место: "Максимальное отклонение вектора ускорения от нормали составляет порядка 60-70 градусов". В обычном случае, угол между нормалью к поверхности и ускорением свободного падения 180 (потому как нормаль нуружу, а ускорение внутрь), и мы видим на комете есть места где похоже этот угол может доходить до 0.


    1. Milliard
      23.11.2015 15:15

      По поводу невыпуклых тел ответил выше. Любое невыпуклое тело можно разбить на конечное число выпуклых. Если для каждого из выпуклых расчёт верен, то почему он должен быть неверен для суммы выпуклых?
      Под отклонением вектора ускорения от нормали, я предполагал отклонение от нормали, направленной вниз, т.е.отклонение 0 градусов означает, что вектор ускорения направлен вертикально вниз. Наверное, стоило написать «отклонение от вертикали», чтобы было понятнее.


      1. Karroplan
        23.11.2015 15:45
        -1

        Любое невыпуклое тело можно разбить на конечное число выпуклых. Если для каждого из выпуклых расчёт верен, то почему он должен быть неверен для суммы выпуклых — это тоже неверно из-за того, что вы не захотели нормально посчитать интеграл по объему.

        Давайте на вашей-же аналогии про площадь (которая точно так же неверна). Вот вам картинка:

        Выбираем в качестве референсной точки для расчета плоащади — точку А. Считаем площади треугольников:
        ABC — ok
        ACD — ok
        ADE — опа, хрень какая-то, такого треугольника в природе не существует
        AFE — опять фигня какая-то
        AFA == 0
        AAA == 0
        AGG == 0
        AGH — еще один несуществующий треугольник
        и т.д.
        и соооовсем неважно какой там знак будет у векторного произведения. Эти площади просто некорректно считать.
        Также и с вашим разбиением.


        1. Aclz
          23.11.2015 16:02
          +4

          ADE — опа, хрень какая-то, такого треугольника в природе не существует
          AFE — опять фигня какая-то

          На примере площади ADEFA. Она равна сумме площадей: +(-) ADA + ADE — AFE +(-) AFA, что после сокращения первого и последнего членов (равных нулю, поэтому у них идет неопределенность по знаку), даст ADE — AFE. Визуально видна корректность полученного результата.


          1. Karroplan
            23.11.2015 16:46

            ок, уговорили!


  1. Aclz
    23.11.2015 17:26
    +1

    В заключение вопрос. Существуют ли тела, у которых хотя бы в одной точке отношение абсолютной величины гравитационного ускорения к ускорению на поверхности шара той же массы и объёма больше единицы?

    Навскидку: предположим, что наблюдатель массы m1 стоит на «северном полюсе» шара радиуса R. Тогда сила тяжести, создаваемая бесконечно малой единицей объема массой dm2 на южном полюсе будет равна F = G (m1 * dm2) / 4R^2. Если обозначить числитель за A, то получим F = A / 4R^2.

    Если посчитать вклад в результирующую силу тяжести, создаваемый такой же единицей объема шара, находящейся на экваторе, то ее тангенциальная составляющая (остающаяся после уравновешивания ее поля полем такой же, но находящейся на противоположной стороне шара частицы) будет равна A / 2v2R^2

    Если посчитать аналогичную составляющую единицы объема, лежащей на 30й параллели северной широты, то значение силы тяжести уже будет равно A / 2R^2 (в 2 раза больше, чем на южном полюсе).

    Соответственно, чем выше поднимаемся по широтам, тем вклад материала поверхности шара в результирующий вектор будет больше (ну по крайней мере в рамках данных точек).
    Таким образом, «отщипнув» единицу материи с южного полюса шара и «прилепив» его на ту же 30ю параллель (сверху нашлёпкой), можно утверждать, что результирующая сила тяжести нового тела будет больше, чем сила тяжести, создаваемая исходным шаром.

    Можно придаться вычислениям, а можно посчитать экстремум функции по сечению сферы и запустить итерационный процесс переброски частичек с наименьшим «вкладом» в окресности частичек с наибольшим и посмотреть, какая форма объекта получится.


    1. Milliard
      23.11.2015 17:52

      Вполне логично. И какой формы будет тело с максимальной силой тяжести на полюсе?


      1. Aclz
        23.11.2015 18:11

        О, пока вычерчивал, уже ответили выше цитатой из какой-то книжки. На бумажке получил примерно ту же лепешку.


  1. Aclz
    25.11.2015 08:50

    А вы ради любопытства не сверялись, Филы садили в зону максимальной силы тяжести по вашим данным?


    1. Milliard
      25.11.2015 11:50

      Филы сел на «макушку головы» кометы. Ускорение там не максимальное, но отклонение от вертикали зато небольшое.


      1. Aclz
        25.11.2015 12:00

        Хм. Мест с отклонением по вертикали-то и так предостаточно, зато гравитация в выбранной точке чуть ли не минимальная из всех возможных вариантов.


        1. Milliard
          25.11.2015 12:23

          Для Фил гравитация не очень важна. Если бы гарпуны закрепись как планировалось, то небольшие вариации силы тяжести не повлияли бы на его работу.


  1. Aclz
    25.11.2015 09:01

    И еще вопрос:

    Для оценки величины отклонения равенства от строгого, я сгенерировал несколько тысяч рандомных тетраэдров и вычислил для них эту величину.
    Как было вычислено точное значение силы тяжести на вершине тетраэдра? Взятием интеграла чисто математически или каким-то машинным обсчетом? (например, разделением треугольника на большое количество маленьких с приравниванием каждого к точечной единице массы)


    1. Milliard
      25.11.2015 11:59

      Точное значение считал разбиением тетраэдра на множество узких тетраэдров с вершиной в заданной точке. Гравитационную силу от них нельзя приравнять к точечной массе, но для них существует точная формула (чем уже, тем точнее). Для некоторых тетраэдров интеграл берётся (например, если в основании равносторонний треугольник), по ним тоже сравнивал.
      И, для контроля, сравнивал вычисленные значения для куба с точным (для куба интеграл тоже существует в элементарных функциях). До шестого знака после запятой значения сходились. При большей точности вычисления становились долгими, но, полагаю, и шести знаков хватит с лихвой.


      1. Aclz
        25.11.2015 12:09

        Для некоторых тетраэдров интеграл берётся (например, если в основании равносторонний треугольник)

        Даже если равносторонний треугольник не перпендикулярен направлению на точку? Если да, то можете поделиться ссылкой? Попытки решения задачи для тетраэдра в общем виде как для потенциала, так и для силы, у меня при интегрировании всякий раз превращалась в такой ацкий трехэтажный матан из логарифмов и арксекансов, что я пока отставил эту затею.


        1. Milliard
          25.11.2015 12:20

          Нет, только если перпендикулярен. Для общего случая у меня тоже нет решения.