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

Тема в разделе "WASM.A&O", создана пользователем Intro, 29 сен 2019.

  1. UbIvItS

    UbIvItS Well-Known Member

    Публикаций:
    0
    Регистрация:
    5 янв 2007
    Сообщения:
    6.201
    Intro, полную физику в игорах никто не реализует, пч она целиком кладётся на проц и обнуляет роль видюхи + для реальной физИки требуется длинную арифу крутить == такое чудо и супер-пупер комп не потянет :laugh1::laugh2::laugh3: тч в гамисах пиши физИку по фене шуеву == юзверям похЪ, им нужна красочность и динамика + умеренная требовательность к компу:grin:
     
  2. Intro

    Intro Active Member

    Публикаций:
    0
    Регистрация:
    29 авг 2009
    Сообщения:
    579
    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 и остальные векторы известны.
    Ну смоделировать работу трёхосных акселерометров и гироскопа. У нас на скорость может что-то внешние повлиять, коллизия например.
     
  3. UbIvItS

    UbIvItS Well-Known Member

    Публикаций:
    0
    Регистрация:
    5 янв 2007
    Сообщения:
    6.201
    Intro, дабл и флоут - такие лапочки :) https://docs.python.org/3/tutorial/floatingpoint.html
    если делать мало-мальски приличный контроль аппроксимации - ошибка останется на высоком уровне, а вот скорости заметно упадут, поэтому в научно-валидных расчётах пользуют именно длинную арифу. :)
     
  4. Intro

    Intro Active Member

    Публикаций:
    0
    Регистрация:
    29 авг 2009
    Сообщения:
    579
    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. --------------------------------------------------------------------------------------------------------------------------
    Как видно правила прямого треугольника. Так же сделал оптимизацию, и получилось весьма хорошо и быстро.
     
  5. UbIvItS

    UbIvItS Well-Known Member

    Публикаций:
    0
    Регистрация:
    5 янв 2007
    Сообщения:
    6.201
    да, при чём там питоха??? == это Вопрос арифы флоутов: 580000 - .00000000001 == 580000. 128 бит - тожЪ слёзы == меньше 256 для приличных расчётов брать не стоит.. в большинстве случаев гораздо дешевле и сильно лучше себя стенды показывают :)
     
  6. alex_dz

    alex_dz Active Member

    Публикаций:
    0
    Регистрация:
    26 июл 2006
    Сообщения:
    389
    при чем тут музикальная арфа?
     
  7. UbIvItS

    UbIvItS Well-Known Member

    Публикаций:
    0
    Регистрация:
    5 янв 2007
    Сообщения:
    6.201
    арфа??? арифа == арифметика :)
     
  8. alex_dz

    alex_dz Active Member

    Публикаций:
    0
    Регистрация:
    26 июл 2006
    Сообщения:
    389
    пердец сокращения
    рехнуцца моно
     
  9. UbIvItS

    UbIvItS Well-Known Member

    Публикаций:
    0
    Регистрация:
    5 янв 2007
    Сообщения:
    6.201
    зато коротко: алго(с) == алгоритм.. хотя, конечно, сокращения тожЪ делают свои проблемы :)
     
  10. Intro

    Intro Active Member

    Публикаций:
    0
    Регистрация:
    29 авг 2009
    Сообщения:
    579
    Нафлудили тут понимаш, ладно.
    С сопротивлением при угле атаки тела вращения типа ракета разобрались. А вот какое будет сопротивлением у тела более сложной формы, но так же симметричное типа самолёт, вертолёт, квадрокоптер? Очевидно от угла атаки 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. В общем, тут предполагается что тело у нас близко к параллелепипеду, и мы приблизительно вычисляем площадь миделя, которая приблизительно равна текущему коэффициенту сопротивления воздуха. Хотя для более точных расчётах можно взять табличный данные по всем этим углам атаки и скорости. Но нам и так сойдёт. К тому же у нас просто нет этих таблиц, это надо обдувать либо реальной трубе, либо в виртуальной.
    В общем, надо ещё раз оптимизировать, там как видно можно убрать тригонометрию.