terminal velocity

Post your FreeBASIC tips and tricks here. Please don’t post your code without including an explanation.
dafhi
Posts: 1329
Joined: Jun 04, 2005 9:51

terminal velocity

Postby dafhi » Jan 29, 2012 2:43

[2017 Sep 17] - accurate, framerate-independent

Code: Select all

#ifndef pi
const   TwoPi = 8*atn(1)
const   Pi = 4*atn(1)
#EndIf

type float as single

type v3
  as float              x,y,z
  declare function      rndf as v3
  declare sub           rand
End Type
sub v3.rand
  y=2*(rnd-.5):  var r=sqr(1-y*y)
  z=rnd*twopi: x=r*cos(z): z=r*sin(z)
End Sub
function 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 function

operator -(r as v3) as v3: return type(-r.x, -r.y, -r.z): end operator
operator -(l as v3,r as v3) as v3: return type(l.x-r.x,l.y-r.y,l.z-r.z): end operator
operator +(l as v3,r as v3) as v3: return type(l.x+r.x, l.y+r.y, l.z+r.z): end operator
operator /(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 operator
operator *(l as v3,r as float) as v3: return type(l.x*r,l.y*r,l.z*r): end operator
operator *(l as float, r as v3) as v3: return type(l*r.x,l*r.y,l*r.z): end operator
operator *(l as v3,r as v3) as v3: return type(l.x*r.x,l.y*r.y,l.z*r.z): end operator


type 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, kj
End Type
constructor.tTermVel(k as single = 5, fps as single=67): ini k, fps
end constructor
sub tTermVel.ini(_k as single, _fps as single):  fps = _fps
  ifps = 1/fps:  k=_k
  new_kinetics:  j=0:  kj=1
End Sub
sub tTermVel.new_kinetics
  ik = 1/(1-k)
  iok = 1/(1-1/k)
End Sub
sub 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=ki
End Sub


const             w = 800, wm=w-1, midx = wm/2
const             h = 600, hm=h-1, midy = hm/2

screenres w,h,32

sub 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
  Next
end sub


var l = 1
var u = 7
for 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)
Next

sleep 2000


Here's the logic behind dpos()

Code: Select all

#ifndef pi
const   TwoPi = 8*atn(1)
const   Pi = 4*atn(1)
#EndIf

type float as single

type v3
  as float              x,y,z
  declare function      rndf as v3
  declare sub           rand
End Type
sub v3.rand
  y=2*(rnd-.5):  var r=sqr(1-y*y)
  z=rnd*twopi: x=r*cos(z): z=r*sin(z)
End Sub
function 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 function

operator -(r as v3) as v3: return type(-r.x, -r.y, -r.z): end operator
operator -(l as v3,r as v3) as v3: return type(l.x-r.x,l.y-r.y,l.z-r.z): end operator
operator +(l as v3,r as v3) as v3: return type(l.x+r.x, l.y+r.y, l.z+r.z): end operator
operator /(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 operator
operator *(l as v3,r as float) as v3: return type(l.x*r,l.y*r,l.z*r): end operator
operator *(l as float, r as v3) as v3: return type(l*r.x,l*r.y,l*r.z): end operator
operator *(l as v3,r as v3) as v3: return type(l.x*r.x,l.y*r.y,l.z*r.z): end operator


function 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 function

function 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 function


const             w = 800, wm=w-1, midx = wm/2
const             h = 600, hm=h-1, midy = hm/2

screenres w,h,32

sub 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
  Next
end sub


var l = 1
var u = 7
for 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)
Next

sleep


function  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

Postby vdecampo » Jan 29, 2012 3:40

FYI:

Terminal velocity without considering buoyancy is

vt = Sqrt(2mg/pAd)

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

Postby dodicat » Feb 06, 2012 12:40

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,yres
screen 20
screeninfo xres,yres
type v2
    as single x,y
end type
dim as v2 startpos,position,lastposition
startpos=type(100,0)
const k=1   
const g=9.81  'gravity
const m=5     'mass
const v=250   'initial velocity
const w=0     'windspeed
dim as single theta=60  'initial angle
theta=theta*(4*atn(1))/180 'degrees to radians

dim as single t

do
  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

Postby dafhi » Feb 07, 2012 1:09

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

Postby dodicat » Feb 07, 2012 21:13

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

Code: Select all

 
dim as integer xres,yres
screen 20
screeninfo xres,yres

const k=1   
const g=9.81  'gravity
const m=5     'mass
const v=200   'initial velocity
const w=0     'windspeed

dim as any pointer im=imagecreate(xres,yres,0)
type v2
    as single x,y
end type
#define r(f,l) Rnd * ((l) - (f)) + (f)
dim as integer n=100
dim as v2 startpos,position,lastposition
redim 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 angle
theta=theta*(4*atn(1))/180 'degrees to radians
dim as single t,y,zz
do
do
  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 n
b(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,1
loop

randomize
startpos=type(r(.4*xres,.6*xres))
theta=r(70,110)
theta=theta*(4*atn(1))/180 'degrees to radians
n=r(250,450)
y=0
setup(n)
lastposition=startpos
imagedestroy(im)
im=imagecreate(xres,yres,0)
t=0
loop until inkey=chr(27)
imagedestroy(im)


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

Re: terminal velocity

Postby dafhi » Sep 17, 2017 10:33

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

Return to “Tips and Tricks”

Who is online

Users browsing this forum: No registered users and 2 guests