заряды

ball.go

// Структура для представления шара в двумерном мире
type Ball struct {
    // Поля для хранения координат центра и радиуса шара
    X float64 // координата x центра шара
    Y float64 // координата y центра шара
    R float64 // радиус шара
    Speed float64 // скорость шара в метрах в секунду
    Density float64 // плотность, масса на единицу площади
    Charge float64 // электрический заряд, может быть как положительным так и отрицательным
}


// Метод для структуры Ball, который обновляет его скорость и положение
func (b *Ball) Update(balls []Ball, dt float64) {
	// Константа Кулона в СИ
	const k = 8.987742438e9
	// Вектор силы, действующей на мяч
	var fx, fy float64
	// Цикл по всем другим мячам
	for _, other := range balls {
		// Пропускаем себя
		if b == &other {
			continue
		}
		// Расстояние между мячами
		dx := other.X - b.X
		dy := other.Y - b.Y
		r := math.Sqrt(dx*dx + dy*dy)
		// Единичный вектор направления
		ux := dx / r
		uy := dy / r
		// Кулоновская сила взаимодействия
		f := k * b.Charge * other.Charge / (r * r)
		// Проекции силы на оси
		fx += f * ux
		fy += f * uy
	}
	// Масса мяча
	m := b.Density * math.Pi * b.R * b.R
	// Ускорение мяча
	ax := fx / m
	ay := fy / m
	// Обновление скорости и положения мяча с помощью неявного метода Рунге-Кутты 4-го порядка
	// Инициализация векторов k и l
	k := make([][2]float64, 4)
	l := make([][2]float64, 4)
	// Начальное приближение для k и l
	k[0][0] = b.SpeedX
	k[0][1] = b.SpeedY
	l[0][0] = ax
	l[0][1] = ay
	// Функция для решения системы нелинейных уравнений методом Ньютона
	solve := func(k, l [][2]float64, i int) {
		// Погрешность
		eps := 1e-6
		// Максимальное число итераций
		maxIter := 100
		// Счетчик итераций
		iter := 0
		// Функция для вычисления невязки
		residual := func(k, l [][2]float64, i int) float64 {
			// Вектор невязки
			var r [2]float64
			// Цикл по всем другим мячам
			for _, other := range balls {
				// Пропускаем себя
				if b == &other {
					continue
				}
				// Расстояние между мячами на следующем шаге
				dx := other.X + other.SpeedX*dt - (b.X + k[i][0]*dt)
				dy := other.Y + other.SpeedY*dt - (b.Y + k[i][1]*dt)
				r := math.Sqrt(dx*dx + dy*dy)
				// Единичный вектор направления
				ux := dx / r
				uy := dy / r
				// Кулоновская сила взаимодействия
				f := k * b.Charge * other.Charge / (r * r)
				// Проекции силы на оси
				fx += f * ux
				fy += f * uy
			}
			// Ускорение мяча на следующем шаге
			ax := fx / m
			ay := fy / m
			// Вычисление невязки
			r[0] = k[i][0] - b.SpeedX - (l[i-1][0] + 2*l[i][0])*dt/6
			r[1] = k[i][1] - b.SpeedY - (l[i-1][1] + 2*l[i][1])*dt/6
			r[2] = l[i][0] - ax - (k[i-1][0] + 2*k[i][0])*dt/6
			r[3] = l[i][1] - ay - (k[i-1][1] + 2*k[i][1])*dt/6
			// Возвращаем норму невязки
			return math.Sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2] + r[3]*r[3])
		}
		// Функция для вычисления Якобиана
		jacobian := func(k, l [][2]float64, i int) [4][4]float64 {
			// Матрица Якобиана
			var J [4][4]float64
			// Шаг дифференцирования
			h := 1e-6
			// Цикл по строкам матрицы
			for j := 0; j < 4; j++ {
				// Копируем векторы k и l
				kc := make([][2]float64, len(k))
				lc := make([][2]float64, len(l))
				copy(kc, k)
				copy(lc, l)
				// Прибавляем шаг к j-му элементу
				if j < 2 {
					kc[i][j] += h
				} else {
					lc[i][j-2] += h
				}
				// Вычисляем невязку с прибавленным шагом
				rp := residual(kc, lc, i)
				// Вычитаем шаг из j-го элемента
				if j < 2 {
					kc[i][j] -= 2 * h
				} else {
					lc[i][j-2] -= 2 * h
				}
				// Вычисляем невязку с вычтенным шагом
				rm := residual(kc, lc, i)
				// Вычисляем производную по j-му элементу
				dr := (rp - rm) / (2 * h)
				// Записываем производную в матрицу Якобиана
				J[j] = dr
			}
			// Возвращаем матрицу Якобиана
			return J
		}
		// Функция для решения системы линейных уравнений методом Гаусса
		// Функция для решения системы линейных уравнений методом Гаусса
		solveLinear := func(A [4][4]float64, b [4]float64) [4]float64 {
			// Вектор решения
			var x [4]float64
			// Прямой ход
			for i := 0; i < 4; i++ {
				// Ищем максимальный элемент в i-м столбце
				max := math.Abs(A[i][i])
				maxRow := i
				for j := i + 1; j < 4; j++ {
					if math.Abs(A[j][i]) > max {
						max = math.Abs(A[j][i])
						maxRow = j
					}
				}
				// Меняем местами i-ю и maxRow-ю строки
				for j := i; j < 4; j++ {
					A[i][j], A[maxRow][j] = A[maxRow][j], A[i][j]
				}
				b[i], b[maxRow] = b[maxRow], b[i]
				// Обнуляем элементы под главной диагональю
				for j := i + 1; j < 4; j++ {
					c := A[j][i] / A[i][i]
					for k := i; k < 4; k++ {
						A[j][k] -= c * A[i][k]
					}
					b[j] -= c * b[i]
				}
			}
			// Обратный ход
			for i := 3; i >= 0; i-- {
				// Вычисляем i-е неизвестное
				x[i] = b[i] / A[i][i]
				// Подставляем в предыдущие уравнения
				for j := i - 1; j >= 0; j-- {
					b[j] -= A[j][i] * x[i]
				}
			}
			// Возвращаем вектор решения
			return x
		}
		// Цикл по итерациям метода Ньютона
		for {
			// Вычисляем невязку
			r := residual(k, l, i)
			// Проверяем условие остановки
			if r < eps || iter > maxIter {
				break
			}
			// Вычисляем Якобиан
			J := jacobian(k, l, i)
			// Решаем систему линейных уравнений J * delta = -r
			delta := solveLinear(J, [-r[0], -r[1], -r[2], -r[3]])
			// Обновляем приближение для k и l
			k[i][0] += delta[0]
			k[i][1] += delta[1]
			l[i][0] += delta[2]
			l[i][1] += delta[3]
			// Увеличиваем счетчик итераций
			iter++
		}
		// Записываем найденные значения k и l
		k[i] = [k[i][0], k[i][1]]
		l[i] = [l[i][0], l[i][1]]
	}
	// Вызываем функцию solve для i = 1, 2, 3
	solve(k, l, 1)
	solve(k, l, 2)
	solve(k, l, 3)
	// Обновляем скорость и положение мяча с помощью метода Рунге-Кутты 4-го порядка
	b.SpeedX += (k[0][0] + 2*k[1][0] + 2*k[2][0] + k[3][0]) * dt / 6
	b.SpeedY += (k[0][1] + 2*k[1][1] + 2*k[2][1] + k[3][1]) * dt / 6
	b.X += (l[0][0] + 2*l[1][0] + 2*l[2][0] + l[3][0]) * dt / 6
	b.Y += (l[0][1] + 2*l[1][1] + 2*l[2][1] + l[3][1]) * dt / 6
}