Стабилизация ракеты в полёте

Intro, полную физику в игорах никто не реализует, пч она целиком кладётся на проц и обнуляет роль видюхи + для реальной физИки требуется длинную арифу крутить == такое чудо и супер-пупер комп не потянет тч в гамисах пиши физИку по фене шуеву == юзверям похЪ, им нужна красочность и динамика + умеренная требовательность к компу

 

UbIvItS, физика в играх использует float для относительно небольших уровней с квадратным радиусом 1-2 км, и double для больших уровней. Более высокая точность не требуется.
В общем вот такая эволюция кода.

Код (Lua):

1. local dir_angular = vector()
2. --------------------------------------------------------------------------------------------------------------------------
3. --v0.03
4. local dir_angular = vector()
5. local force_yaw_pitch = angular_yaw + vel_angular_pitch
6. if force_yaw_pitch~=0 then
7.     dir_angular:rotate_angle90(axis_yaw, axis_pitch, vel_angular_pitch/force_yaw_pitch*Deg90)
8. end
9. local force_angular = force_yaw_pitch + vel_angular_roll
10. if force_angular~=0 then
11.     dir_angular:rotate_angle90(dir_angular, axis_roll, vel_angular_roll/force_angular*Deg90)
12. end
13. ph_shell:set_angular_vel(dir_angular:mul(force_angular))
14. --------------------------------------------------------------------------------------------------------------------------
15. --v0.04
16. if angular_roll<0   then axis_roll:invert()  end
17. if angular_yaw<0    then axis_yaw:invert()   end
18. if angular_pitch<0  then axis_pitch:invert() end
19. if angular_yaw~=0 then
20.     dir_angular:rotate_angle90(axis_yaw, axis_pitch, math.atan(math.abs(angular_pitch/angular_yaw)))
21. else
22.     dir_angular:set(axis_pitch)
23. end
24. local force_angular =  math.sqrt(angular_yaw*angular_yaw + angular_pitch*angular_pitch)
25. if force_angular~=0 then
26.     dir_angular:rotate_angle90(dir_angular, axis_roll, math.atan(math.abs(angular_roll/force_angular)))
27. else
28.     dir_angular:set(axis_roll)
29. end
30. force_angular = math.sqrt(force_angular*force_angular + angular_roll*angular_roll)
31. ph_shell:set_angular_vel(dir_angular:mul(force_angular))
32. --------------------------------------------------------------------------------------------------------------------------
33. --v0.05
34. local force_angular = angular_roll
35. local force_yaw_pitch = math.sqrt(angular_yaw^2 + angular_pitch^2)
36. if force_yaw_pitch~=0 then
37.     local scale0, scale1 = angular_yaw/force_yaw_pitch, angular_pitch/force_yaw_pitch
38.     dir_angular.x = axis_yaw.x*scale0 + axis_pitch.x*scale1
39.     dir_angular.y = axis_yaw.y*scale0 + axis_pitch.y*scale1
40.     dir_angular.z = axis_yaw.z*scale0 + axis_pitch.z*scale1
41.     force_angular = math.sqrt(force_yaw_pitch^2 + angular_roll^2)
42.     scale0, scale1 = force_yaw_pitch/force_angular, angular_roll/force_angular
43.     dir_angular.x = dir_angular.x*scale0 + axis_roll.x*scale1
44.     dir_angular.y = dir_angular.y*scale0 + axis_roll.y*scale1
45.     dir_angular.z = dir_angular.z*scale0 + axis_roll.z*scale1
46. else
47.     dir_angular:set(axis_roll)
48. end
49. --------------------------------------------------------------------------------------------------------------------------
50. --v0.05a
51. local force_angular = math.sqrt(angular_yaw^2 + angular_pitch^2 + angular_roll^2)
52. if force_angular>0 then
53.     local scale0, scale1, scale2 = angular_yaw/force_angular, angular_pitch/force_angular, angular_roll/force_angular
54.     dir_angular.x = axis_yaw.x*scale0 + axis_pitch.x*scale1 + axis_roll.x*scale2
55.     dir_angular.y = axis_yaw.y*scale0 + axis_pitch.y*scale1 + axis_roll.y*scale2
56.     dir_angular.z = axis_yaw.z*scale0 + axis_pitch.z*scale1 + axis_roll.z*scale2
57. end
58. dir_angular:mul(force_angular)
59. --------------------------------------------------------------------------------------------------------------------------
60. --v0.05b
61. local force_angular = math.sqrt(angular_yaw^2 + angular_pitch^2 + angular_roll^2)
62. if force_angular>0 then
63.     local scale0, scale1, scale2 = angular_yaw/force_angular, angular_pitch/force_angular, angular_roll/force_angular
64.     dir_angular.x = (axis_yaw.x*scale0 + axis_pitch.x*scale1 + axis_roll.x*scale2)*force_angular
65.     dir_angular.y = (axis_yaw.y*scale0 + axis_pitch.y*scale1 + axis_roll.y*scale2)*force_angular
66.     dir_angular.z = (axis_yaw.z*scale0 + axis_pitch.z*scale1 + axis_roll.z*scale2)*force_angular
67. end
68. --------------------------------------------------------------------------------------------------------------------------
69. --v0.06
70. dir_angular.x = axis_yaw.x*angular_yaw + axis_pitch.x*angular_pitch + axis_roll.x*angular_roll
71. dir_angular.y = axis_yaw.y*angular_yaw + axis_pitch.y*angular_pitch + axis_roll.y*angular_roll
72. dir_angular.z = axis_yaw.z*angular_yaw + axis_pitch.z*angular_pitch + axis_roll.z*angular_roll
73. --------------------------------------------------------------------------------------------------------------------------

Вот так. Сейчас надо найти формулу, где надо восстановить angular_yaw, angular_pitch, angular_roll, если dir_angular и остальные векторы известны.
Ну смоделировать работу трёхосных акселерометров и гироскопа. У нас на скорость может что-то внешние повлиять, коллизия например.

 

Intro, дабл и флоут - такие лапочки https://docs.python.org/3/tutorial/floatingpoint.html
если делать мало-мальски приличный контроль аппроксимации - ошибка останется на высоком уровне, а вот скорости заметно упадут, поэтому в научно-валидных расчётах пользуют именно длинную арифу.

 

UbIvItS, питухоний? Это бред для питухонцев у которых проблемы с математикой, даже с элементарной. У любого моделирования есть погрешность, для научных часто double не хватает, используют четвертную. А для игр используется даже половинная в два байта, правда с этим типом работает В/К, для некоторых рассчётов.

Насчёт восстановить angular_yaw, angular_pitch, angular_roll. Вот:

Код (Lua):

1. --------------------------------------------------------------------------------------------------------------------------
2. --v0
3. mag = dir_angular:magnitude()
4. dir_angular:normalize()
5. angle = math.acos(axis_yaw:dotproduct(dir_angular))
6. angular_yaw = math.cos(angle)*mag
7. angle = math.acos(axis_pitch:dotproduct(dir_angular))
8. angular_pitch = math.cos(angle)*mag
9. angle = math.acos(axis_roll:dotproduct(dir_angular))
10. angular_roll = math.cos(angle)*mag
11. --------------------------------------------------------------------------------------------------------------------------
12. --v1
13. mag = dir_angular:magnitude()
14. dir_angular:normalize()
15. angular_yaw = axis_yaw:dotproduct(dir_angular)*mag
16. angular_pitch = axis_pitch:dotproduct(dir_angular)*mag
17. angular_roll = axis_roll:dotproduct(dir_angular)*mag
18. --------------------------------------------------------------------------------------------------------------------------
19. --v2
20. angular_yaw = axis_yaw:dotproduct(dir_angular)
21. angular_pitch = axis_pitch:dotproduct(dir_angular)
22. angular_roll = axis_roll:dotproduct(dir_angular)
23. --------------------------------------------------------------------------------------------------------------------------

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

 

да, при чём там питоха??? == это Вопрос арифы флоутов: 580000 - .00000000001 == 580000. 128 бит - тожЪ слёзы == меньше 256 для приличных расчётов брать не стоит.. в большинстве случаев гораздо дешевле и сильно лучше себя стенды показывают

 

при чем тут музикальная арфа?

 

арфа??? арифа == арифметика

 

пердец сокращения
рехнуцца моно

 

зато коротко: алго(с) == алгоритм.. хотя, конечно, сокращения тожЪ делают свои проблемы

 

Нафлудили тут понимаш, ладно.
С сопротивлением при угле атаки тела вращения типа ракета разобрались. А вот какое будет сопротивлением у тела более сложной формы, но так же симметричное типа самолёт, вертолёт, квадрокоптер? Очевидно от угла атаки dir_vel:dotproduct(axis_roll), а так же от крена.

Код (Lua):

1. --------------------------------------------------------------------------------------------------------------------------
2. --v0
3. --сила сопротивления воздуха...
4. local lvel, dir_vel = vector(), vector()
5. ph_shell:get_linear_vel(lvel)
6. local vel, k_air = lvel:magnitude(), 1
7. local atack_pitch, atack_roll = 0, 0
8. if vel>0.005 then
9.     --зависит от угла атаки тангажа и крена, зависимость синусоидальная.
10.     dir_vel:set(lvel):normalize()
11.     atack_pitch = math.acos(clamp(dir_vel:dotproduct(axis_roll), -1.0, 1.0))                --тангаж
12.     if atack_pitch>0.001 then
13.         local pitch90 = vector():slerp(axis_roll, dir_vel, Deg90/atack_pitch)
14.         atack_roll = math.acos(clamp(pitch90:dotproduct(axis_pitch), -1.0, 1.0)) - Deg90    --крен
15.     end
16.     local angle_sin = math.sin(atack_pitch)
17.     local k_max = k_front + angle_sin*(k_top-k_front)
18.     local k_min = k_front + angle_sin*(k_side-k_front)
19.     k_air = k_min + math.cos(atack_roll)*(k_max-k_min)
20. end
21. local k = -vel*k_air*self.mass100*fDelta
22. --------------------------------------------------------------------------------------------------------------------------
23. --v1
24. --сила сопротивления воздуха...
25. local lvel, dir_vel = vector(), vector()
26. ph_shell:get_linear_vel(lvel)
27. local vel, k_air = lvel:magnitude(), 1
28. local atack_pitch, atack_roll = 0, 0
29. if vel>0.005 then
30.     --зависит от угла атаки тангажа и крена, зависимость синусоидальная.
31.     dir_vel:set(lvel):normalize()
32.     atack_pitch = math.acos(clamp(dir_vel:dotproduct(axis_roll), -1.0, 1.0))                --тангаж
33.     if atack_pitch>0.001 then
34.         local pitch90 = vector():crossproduct(dir_vel, axis_roll):normalize()
35.         atack_roll = math.acos(clamp(pitch90:dotproduct(axis_pitch), -1.0, 1.0))            --крен
36.     end
37.     local angle_sin = math.sin(atack_pitch)
38.     local k_max = k_front + angle_sin*(k_top-k_front)
39.     local k_min = k_front + angle_sin*(k_side-k_front)
40.     k_air = k_min + math.cos(math.abs(atack_roll))*(k_max-k_min)
41. end
42. local k = -vel*k_air*self.mass100*fDelta

Я тут конечно выбрал сложный путь, и в v0 нагородил много кода, потом вспомнил crossproduct, которая сильно проще slerp. В общем, тут предполагается что тело у нас близко к параллелепипеду, и мы приблизительно вычисляем площадь миделя, которая приблизительно равна текущему коэффициенту сопротивления воздуха. Хотя для более точных расчётах можно взять табличный данные по всем этим углам атаки и скорости. Но нам и так сойдёт. К тому же у нас просто нет этих таблиц, это надо обдувать либо реальной трубе, либо в виртуальной.
В общем, надо ещё раз оптимизировать, там как видно можно убрать тригонометрию.