Не могу понять в чем бага....blitz

PooH

прога про три планеты, их гравитационное взаимодействие
на производительность пох
пишется на блице3д
для незнающих:
planet(i)\a[k] - это поле а[k] списка planet(i)
Короче, вроде все правильно, только прога при g=50 не меняет картинку (ускорения = 0 а при g>60 все работает слишком быстро - т.е. планеты разлетаются за миллисекунды
 

  Graphics3D 640,480,32,2

SetBuffer BackBuffer







Const n=3

Global g#=60

Global dt#=1



;начальные координаты, скорости и массы

Dim pos(n-1,n-1)

pos(0,0)=50

pos(0,1)=0

pos(0,2)=0

pos(1,0)=0

pos(1,1)=50

pos(1,2)=0

pos(2,0)=0

pos(2,1)=0

pos(2,2)=50



Dim v(n-1,n-1)

v(0,0)=0.01

v(0,1)=0.00

v(0,2)=0.02



v(1,0)=0.

v(1,1)=0.

v(1,2)=0.



v(2,0)=0.

v(2,1)=0.

v(2,2)=0.



Dim mass#(n-1)

mass(0)=1

mass(1)=1

mass(2)=1

;---------------------



;создание камеры и освещения

pivot=CreatePivot

cam=CreateCamera(pivot)

MoveEntity cam,0,0,-15000

light=CreateLight

;-----------------------------





;создаем планеты

Type planets

Field x[n-1]

Field v#[n-1]

Field a#[n-1]

Field mass#

Field mesh

End Type

 

Dim planet.planets(n-1)



For i=0 To n-1

planet(i)=New planets

planet(i)\x[0]=pos(i,0)

planet(i)\x[1]=pos(i,1)

planet(i)\x[2]=pos(i,2)

planet(i)\v[0]=v(i,0)

planet(i)\v[1]=v(i,1)

planet(i)\v[2]=v(i,2)

planet(i)\mass=mass(i)

planet(i)\mesh=CreateSphere

Next

;-------------------------------



;расчитываем центр масс для камеры

M#=planet(0)\mass+planet(1)\mass+planet(2)\mass

xp=(pos(0,0)*mass(0)+pos(1,0)*mass(1)+pos(2,0)*mass(2/M

yp=(pos(0,1)*mass(0)+pos(1,1)*mass(1)+pos(2,1)*mass(2/M

zp=(pos(0,2)*mass(0)+pos(1,2)*mass(1)+pos(2,2)*mass(2/M

PositionEntity pivot,xp,yp,zp

;-----------------------





Stop

;основной цикл

Repeat

oldtime=MilliSecs

CalculateWorld_runge









;TurnEntity pivot,0,1,0

movemesh

UpdateWorld

RenderWorld

Flip

Delay(100)

dt=(MilliSecs-oldtime)



Until KeyHit(1)



End

;-------------------------------







Function CalculateWorld_runge



planet(0)\a[0] = (F(0,0,1)+F(0,0,2/(planet(0)\mass)



;planet(0)\a[0] = (50.9)/(planet(0)\mass)



planet(0)\a[1] = (F(1,0,1)+F(1,0,2/(planet(0)\mass)

planet(0)\a[2] = (F(2,0,1)+F(2,0,2/(planet(0)\mass)



DebugLog "planet(0)\a[0]=" + Str(planet(0)\a[0])

;DebugLog "planet(0)\a[1]=" + Str(planet(0)\a[1])

;DebugLog "planet(0)\a[2]=" + Str(planet(0)\a[2])



planet(1)\a[0]=(F(0,1,0)+F(0,1,2/planet(1)\mass

planet(1)\a[1]=(F(1,1,0)+F(1,1,2/planet(1)\mass

planet(1)\a[2]=(F(2,1,0)+F(2,1,2/planet(1)\mass

;DebugLog "planet(1)\a[0]=" + Str(planet(1)\a[0])

;DebugLog "planet(1)\a[1]=" + Str(planet(1)\a[1])

;DebugLog "planet(1)\a[2]=" + Str(planet(1)\a[2])



planet(2)\a[0]=(F(0,2,1)+F(0,2,1/planet(2)\mass

planet(2)\a[1]=(F(1,2,1)+F(1,2,1/planet(2)\mass

planet(2)\a[2]=(F(2,2,1)+F(2,2,1/planet(2)\mass



For i=0 To n-1

planet(i)\v[0]=planet(i)\v[0] + planet(i)\a[0]*dt

planet(i)\v[1]=planet(i)\v[1] + planet(i)\a[1]*dt

planet(i)\v[2]=planet(i)\v[2] + planet(i)\a[2]*dt





;DebugLog "planet("+Str(i)+")\v[1]=" + Str(planet(i)\v[1])

;DebugLog "planet("+Str(i)+")\v[2]=" + Str(planet(i)\v[2])





planet(i)\x[0]=planet(i)\x[0]+planet(i)\v[0]*dt

planet(i)\x[1]=planet(i)\x[1]+planet(i)\v[1]*dt

planet(i)\x[2]=planet(i)\x[2]+planet(i)\v[2]*dt

Next

DebugLog "dt="+Str(dt)

DebugLog "planet("+Str(1)+")\v[0]=" + Str(planet(1)\v[0])



End Function



;функция расчета силы для i-той проекции планет n1 и n2

Function F(i,n1,n2)

;DebugLog "f___" +Str(G* planet(n1)\mass * planet(n2)\mass/L(n1,n2

F# = G * planet(n1)\mass * planet(n2)\mass / L(n1,n2)

;DebugLog "g="+Str(g) +"   "+"planet("+Str(n1)+") mass=" + Str(planet(n1)\mass)

;DebugLog "g="+Str(g) +"   "+"planet("+Str(n2)+") mass=" + Str(planet(n2)\mass)

;DebugLog "F"+Str(i)+"_"+Str(n1)+"_"+Str(n2)+"=" + Str(F) 

f#=F*(planet(n1)\x[i]-planet(n2)\x[i])/L(n1,n2)



DebugLog "dx="+Str(planet(n1)\x[i]-planet(n2)\x[i])+"	"+"f=" + Str(F)

Return -f

End Function





Function L(n1,n2) 

l#planet(n1)\x[0]-planet(n2)\x[0])^2+(planet(n1)\x[1]-planet(n2)\x[1])^2+(planet(n1)\x[2]-planet(n2)\x[2])^2)^0.5

Return l

End Function 

;-------------------------------





Function movemesh

For i=0 To n-1

PositionEntity planet(i)\mesh,planet(i)\x[0],planet(i)\x[1],planet(i)\x[2]

Next

End Function 

PooH

оч надо найти причину до завтрашнего вечера
просьба бэсик говно - не писать

PooH

короче, вроде понятно, что бага в округлении, но я никак не могу догнать, как ее обойти...

мож, кто какие костыли придумает?

pinstripe

Потести с двумя планетами с простыми начальными данными, чтобы первые несколько кадров ты мог сам вручную просчитать.
Кроме того прикинь на сколько у тебя смещаются планеты за один кадр. Если это расстояние больше 1/10 от расстояния между планетами то точность алгоритма будет плохой.
Ну и последнее, если две планеты сильно сблизятся, то из-за дискретности алгоритма они приобретут слишком большое ускорение на один кадр и очень быстро разлетятся. Как решение можно уменьшать dt.

pinstripe

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

PooH

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

pinstripe

Можно придать им начальные скорости так, чтобы они стабильно вращались не сближаясь. Для двух это легко. Для трех тоже можно, по крайней мере у Солнца Земли и Луны получается. Еще можно уменьшать dt с уменьшением расстояния между планетами, но если они слипаются, а не просто пролетают рядом, то это не поможет.

PooH

угу я понял
для двух тел все работает замечательно - тела вертятся вокруг центра масс
ща допишу перемещение камеры с центром масс и заценю для трех

PooH

чет я седня уже не могу писать....или центр масс гуляет
завтра продолжу...
большая просьба писавшим на этом....посоветовать что-нибудь

agaaaa

для начала я бы сделал все константы (в т.ч. и нулевые) флоатами - т.е. добавил бы .0 после каждой целочисленной константы

danilov

В таких задачах есть смысл работу сил на каждом дискретном шаге честно посчитать.
В случае пересечения 2х планет (проход через 0 на отрезке [-a, b] считать, что для всякого x работа при прохождении отрезка [-x, x] равна 0.
Конечно, будет погрешность за счёт поочерёдного перемещения. Тут уж не знаю, что делать.