## terminal velocity

dafhi
Posts: 1329
Joined: Jun 04, 2005 9:51

### terminal velocity

[2017 Sep 17] - accurate, framerate-independent

Code: Select all

`#ifndef piconst   TwoPi = 8*atn(1)const   Pi = 4*atn(1)#EndIftype float as singletype v3  as float              x,y,z  declare function      rndf as v3  declare sub           randEnd Typesub v3.rand  y=2*(rnd-.5):  var r=sqr(1-y*y)  z=rnd*twopi: x=r*cos(z): z=r*sin(z)End Subfunction v3.rndf as v3  var y=2*(rnd-.5):  var r=sqr(1-y*y)  var z=rnd*twopi:  return type(r*cos(z), y, r*sin(z))End functionoperator -(r as v3) as v3: return type(-r.x, -r.y, -r.z): end operatoroperator -(l as v3,r as v3) as v3: return type(l.x-r.x,l.y-r.y,l.z-r.z): end operatoroperator +(l as v3,r as v3) as v3: return type(l.x+r.x, l.y+r.y, l.z+r.z): end operatoroperator /(l as v3,r as float) as v3: dim as float s = 1/r: return type(l.x*s,l.y*s,l.z*s): end operatoroperator *(l as v3,r as float) as v3: return type(l.x*r,l.y*r,l.z*r): end operatoroperator *(l as float, r as v3) as v3: return type(l*r.x,l*r.y,l*r.z): end operatoroperator *(l as v3,r as v3) as v3: return type(l.x*r.x,l.y*r.y,l.z*r.z): end operatortype tTermVel  as single     v, a  as single     k=.5, fps=67  declare sub   ini(k as single = .5, fps as single=67)  declare sub   dpos(i as single)  declare sub   new_kinetics  as single     ifps=1/67  declare constructor(k as single = 5, fps as single=67) private:  as single     ik, iok, i, j, kjEnd Typeconstructor.tTermVel(k as single = 5, fps as single=67): ini k, fpsend constructorsub tTermVel.ini(_k as single, _fps as single):  fps = _fps  ifps = 1/fps:  k=_k  new_kinetics:  j=0:  kj=1End Subsub tTermVel.new_kinetics  ik = 1/(1-k)  iok = 1/(1-1/k)End Subsub tTermVel.dpos(i as single)  var ki = k^i  v = k*(kj-ki)  a = ((i-j)*k + v*iok)*ik  v *= ik  j=i: kj=kiEnd Subconst             w = 800, wm=w-1, midx = wm/2const             h = 600, hm=h-1, midy = hm/2screenres w,h,32sub insta_plot(fps as single, col as ulong=-1, frames as ubyte=30)  dim as v3       p = type(100, 200, 0)  dim as v3       v = type(20,20) * 3  dim as v3       a = type(0,-1,0) * 3  var             k = .85  dim as tTermVel t=type(k, fps)  var             ifps = 1/fps, age = rnd*ifps  for n as long = 1 to fps*frames    t.dpos age    p += v*t.v + a*t.a    pset ( p.x, hm - p.y ), col    age += ifps  Nextend subvar l = 1var u = 7for fps as single = u to l step -1  var col0 = 70  var col = col0 + (255-col0)*(u+1-fps)/(u+1-l)  insta_plot fps, rgb(col,col,col)Nextsleep 2000`

Here's the logic behind dpos()

Code: Select all

`#ifndef piconst   TwoPi = 8*atn(1)const   Pi = 4*atn(1)#EndIftype float as singletype v3  as float              x,y,z  declare function      rndf as v3  declare sub           randEnd Typesub v3.rand  y=2*(rnd-.5):  var r=sqr(1-y*y)  z=rnd*twopi: x=r*cos(z): z=r*sin(z)End Subfunction v3.rndf as v3  var y=2*(rnd-.5):  var r=sqr(1-y*y)  var z=rnd*twopi:  return type(r*cos(z), y, r*sin(z))End functionoperator -(r as v3) as v3: return type(-r.x, -r.y, -r.z): end operatoroperator -(l as v3,r as v3) as v3: return type(l.x-r.x,l.y-r.y,l.z-r.z): end operatoroperator +(l as v3,r as v3) as v3: return type(l.x+r.x, l.y+r.y, l.z+r.z): end operatoroperator /(l as v3,r as float) as v3: dim as float s = 1/r: return type(l.x*s,l.y*s,l.z*s): end operatoroperator *(l as v3,r as float) as v3: return type(l.x*r,l.y*r,l.z*r): end operatoroperator *(l as float, r as v3) as v3: return type(l*r.x,l*r.y,l*r.z): end operatoroperator *(l as v3,r as v3) as v3: return type(l.x*r.x,l.y*r.y,l.z*r.z): end operatorfunction gsum(k as single=.5, n as single=2) as double  /' --- geometric sum ---    1. s    =  k+k^2 .. k^n    2. s-ks =  k+k^2 .. k^n            -    k^2 .. k^n - k^(n+1)               ____________________            =  k  -  k^(n+1)        '/  return (k-k^(n+1)) / (1-k)end functionfunction vgsum(k as single, i as single) as single  /' --- velo grav sum ---    Loop   [ vel+=accel, vel*=k, pos+=vel ]       vel frame 1:  k    pos frame 1: k    vel frame 2:  kk+k    pos frame 2: kk+k+k     - solving for pos -    p      = 2k1 + 1k2    pk     = 2k2 + 1k3    p(1-k) = 2k1 + 1k2 - 2k2 - 1k3    f      = 2k1    g      =      - 1k2 - 1k3    g/k    = - 1k - 1k2    g(1-1/k) = - k3 + k    g        = (k-k3) / (1 - 1/k)    f+g      = 2k + (k-k3) / (1-1/k)    p        =(2k + (k-k3) / (1-1/k)) / (1-k)  '/  return ( i*k + (k-k^(i+1)) / (1-1/k) ) / (1-k)end functionconst             w = 800, wm=w-1, midx = wm/2const             h = 600, hm=h-1, midy = hm/2screenres w,h,32sub insta_plot(fps as single, col as ulong=-1, frames as ubyte=30)  dim as v3       p = type(100, 200, 0)  dim as v3       v = type(20,20) * 3  dim as v3       a = type(0,-1,0) * 3  var             k = .85  var             ifps = 1/fps  var             age = rnd*ifps  var             ik = 1/(1-k)  var             iok = 1/(1-1/k)  var             j=csng(0)  for n as long = 1 to fps*frames    p += v*(gsum(k,age)-gsum(k,j)) + _  'pos delta is based upon absolute pos: current - previous         a*(vgsum(k,age)-vgsum(k,j))                                            'logic for absolute pos is shown in function tvel(),                                        'at the end of this source code    pset ( p.x, hm - p.y ), col        j = age    age += ifps  Nextend subvar l = 1var u = 7for fps as single = u to l step -1  var col0 = 70  var col = col0 + (255-col0)*(u+1-fps)/(u+1-l)  insta_plot fps, rgb(col,col,col)Nextsleepfunction  tvel(a as single, k as single, t as single, v0 as single=0) as single  /' --- terminal velocity - solved using "recurrence relations"     'pos(time+1) = pos(time) + (vel(time) + accel)*k        fn+1 = fn+(vn+a)k  <--- shorthand        fn+2 = fn+1 + (vn+1 + a)k         = fn+(vn+a)k + (vn+1 + a)k    vn+1 = (vn+a)k    fn+2 = fn+(vn+a)k + ( (vn+a)k + a)k         = fn + (vn+a)k + (vn+a)kk + ak    fn+3 = fn+1 + (vn+1+a)k + (vn+1+a)kk + ak         = fn+(vn+a)k + ((vn+a)k+a)k + ((vn+a)k+a)kk + ak         = fn + vnk + ak + vnkk + akk + ak + vnkkk + akkk + akk + ak         = fn + vn*(k+kk+kkk) + a*(3k + 2kk + 1kkk)                  = fn + vn*geom sum   + a*poly  '/    'return v0*gsum(k,j) + a*vgsum(k,j)  dim as single kpow = k-k^(t+1)  return ( v0*kpow + a*( t*k + kpow/(1-1/k) ) )/(1-k)  End Function`
Last edited by dafhi on Sep 23, 2017 3:35, edited 27 times in total.
vdecampo
Posts: 2982
Joined: Aug 07, 2007 23:20
Location: Maryland, USA
Contact:

### Re: Terminal Velocity

FYI:

Terminal velocity without considering buoyancy is

where
Vt = terminal velocity,
m = mass of the falling object,
g = acceleration due to gravity,
d = drag coefficient,
ρ = density of the fluid through which the object is falling, and
A = projected area of the object.

-Vince
dodicat
Posts: 6234
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: terminal velocity sim

Hi dafhi

Did a thing years ago kinda similar.
Here I didn't use any viscosity formulas (Stoke's law or the like), but just assumed that the drag on an object through a medium acted directly opposite to it's velocity.
i.e. the starter equation was:
m*(dv/dt)=-(mg+kv)
m=mass,g=gravity, k is some constant, dv and dt are differentials.
Then I added a wind factor, so here's a body being fired upwards under the effect of gravity and air resistance applied as per above.

Code: Select all

`dim as integer xres,yresscreen 20screeninfo xres,yrestype v2    as single x,yend typedim as v2 startpos,position,lastpositionstartpos=type(100,0)const k=1   const g=9.81  'gravityconst m=5     'massconst v=250   'initial velocityconst w=0     'windspeeddim as single theta=60  'initial angletheta=theta*(4*atn(1))/180 'degrees to radiansdim as single tdo  t=t+.01      position.x=startpos.x+(m/k)*(1-exp(-(k/m)*t))*(V*cos(theta)-w)+w*t    position.y=startpos.y+(m/k)*(1-exp(-(k/m)*t))*(V*sin(theta)+g*m/k)-(g*m/k)*t    screenlock    locate 2,2    print "Speed "; 100*sqr((position.x-lastposition.x)^2+(position.y-lastposition.y)^2)    circle (position.x,yres-position.y),5,,,,,f    lastposition=position    screenunlock    sleep 1,1    loop until inkey=chr(27)  `
dafhi
Posts: 1329
Joined: Jun 04, 2005 9:51

### Re: terminal velocity sim

Interesting that k seems to represent the density of the surrounding medium, and that you've used e.
It gives a satisfying result.

I hope to be working with sequence sums again soon because I love exploring this stuff
dodicat
Posts: 6234
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: terminal velocity sim

The question is Dafhi:
Is this stuff Rocket Science or a bad Hair Day?

Code: Select all

` dim as integer xres,yresscreen 20screeninfo xres,yresconst k=1   const g=9.81  'gravityconst m=5     'massconst v=200   'initial velocityconst w=0     'windspeeddim as any pointer im=imagecreate(xres,yres,0)type v2    as single x,yend type#define r(f,l) Rnd * ((l) - (f)) + (f)dim as integer n=100dim as v2 startpos,position,lastpositionredim shared as v2 b()redim shared as single ang(),vel(),rad()#macro setup(n)redim b(1 to n)redim rad(1 to n)redim vel(1 to n)redim ang(1 to n)for z as integer=1 to n    rad(z)=(8-3)*(z-1)/(n-1)+3    ang(z)=r(0,360)    ang(z)=ang(z)*(4*atn(1))/180    vel(z)=r(10,50)  next z   #endmacro    setup(n)startpos=type(100,0)dim as single theta=60  'initial angletheta=theta*(4*atn(1))/180 'degrees to radiansdim as single t,y,zzdodo  t=t+.03      position.x=startpos.x+(m/k)*(1-exp(-(k/m)*t))*(V*cos(theta)-w)+w*t    position.y=startpos.y+(m/k)*(1-exp(-(k/m)*t))*(V*sin(theta)+g*m/k)-(g*m/k)*t    screenlock    cls    put(0,0),im    pset im,(position.x,yres-position.y),2    circle (position.x,yres-position.y),5,4,,,,f        if position.y<lastposition.y then        screenunlock        startpos=type(position.x,position.y)        t=0        do            t=t+.03            screenlock            cls            put(0,0),im            for z as integer=1 to nb(z)=type(startpos.x+(m/k)*(1-exp(-(k/m)*t))*(vel(z)*cos(ang(z))-w)+w*t,_       startpos.y+(m/k)*(1-exp(-(k/m)*t))*(vel(z)*sin(ang(z))+g*m/k)-(g*m/k)*t)       circle(b(z).x,yres-b(z).y),rad(z),z,,,,f       pset im,(b(z).x,yres-b(z).y),z       if y<b(z).y then y=b(z).y:zz=z   next z   screenunlock   sleep 1,1   if inkey=chr(27) then exit do,do,do   if b(zz).y<0 then exit do,do   loop         end if    lastposition=position    screenunlock    sleep 1,1looprandomizestartpos=type(r(.4*xres,.6*xres))theta=r(70,110)theta=theta*(4*atn(1))/180 'degrees to radiansn=r(250,450)y=0setup(n)lastposition=startposimagedestroy(im)im=imagecreate(xres,yres,0)t=0loop until inkey=chr(27)imagedestroy(im) `
dafhi
Posts: 1329
Joined: Jun 04, 2005 9:51

### Re: terminal velocity

updated. new format will let me do all kinds of neat tricks