terminal velocity

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

terminal velocity

Post by dafhi »

Finally. ( well, almost. identified a problem with my pre-loop logic )

[2020 Oct 13] - fast, accurate, framerate-independent

Code: Select all

/' -- fast terminal velocity - 2020 Oct 13 - by dafhi

  goal:  fps-independence and simple physics
  
  currently
  1. acceleration is not fps-independent
  2. necessary pre-calcs aren't straight-forward
 
    here's the process:
    
  dim sng dens  = .6 '' > 0 and less than 1
  
  dim sng v     = velocity_const / velo_sum(dens, phys_fps)
  dim sng a     = accel_const / accel_sum(dens, phys_fps)
  dim sng p     '' pos
  dim sng k     = dens ^ (1 / phys_fps)
  
  for i int = 0 to phys_fps*5 '' extend the curve
  
    '' horizontal and vertical scaling
    pset (i * 130 / phys_fps, h-1 - p * 99), col
    
    v += a
    v *= k
    p += v
  Next
  
'/

'' replaces Int()
#Define flo(x) (((x)*2.0-0.5)shr 1) '' http://www.freebasic.net/forum/viewtopic.php?p=118633

'' hacks
#undef int

#define int as integer
#define sng as single


function gs(k sng = .5, n sng = 2) sng

  /'
      -- geometric sum
      
    sum         =  k   + k^2 .. k^n
    sum - k*sum =  k   + k^2 .. k^n
                 - k^2 - k^3 .. k^(n+1)
    
    = k + k^2 .. k^n
        - k^2 .. k^n - k^(n+1)
  '/

  return k*(1-k^n) / (1-k)
  
end function


function sgs(k sng, fps sng) sng
  
  /' 
      -- sum of geometric sum
     
     inputs    particle_density, physics fps
     output
    position at time t
  
      - theory
      
    can summarize with a loop:
    [ vel+=accel, vel*=k, pos+=vel ]
   
    vel frame 1: k              'velocity = (velocity + accel) * k
    vel frame 2: k+kk           'notice k in front each time (acceleration)
    vel frame 3: k kk kkk       'geometric sum
    
    pos frame 1: k              'position = velocity sum = sum of geometric sum
    pos frame 2: k  k kk
    pos frame 3: k  k kk  k kk kkk
    
    the pattern becomes visible to the experienced human
    
    geometric sum        = k*(1-k^n) / (1-k)
    
    sum of geometric sum = k*(1-k^1)/(1-k) + .. k*(1-k^n)/(1-k)
    
    representation 'r'
    k/(1-k)
    
    series
    1-k^1 + .. 1-k^n
    = n - [k^1 + .. k^n]
    = n - geometric sum

    solution
    = r*( n - k*(1-k^n)/(1-k) )
    = r*( n - r*(1-k^n) )
  '/
  
  static sng r
  r=k/(1-k)
  
  return r*(fps - r*(1-k^fps))
  
end function

function velo_sum(k sng = .5, fps sng = 2) sng
  return gs( k^(1/fps), fps ) / k
end function

function accel_sum(k sng, fps sng) sng
  return sgs( k^(1/fps), fps ) / k
end function


const             w = 800, wm=w-1
const             h = 600, hm=h-1

sub fps_plot( phys_fps sng, col as ulong = -1 )
  
  dim sng dens  = .6
  
  dim sng v     = 2 / velo_sum(dens, phys_fps)
  dim sng a     = 0.1 / accel_sum(dens, phys_fps)
  dim sng p     '' pos
  
  dim sng k     = dens ^ (1 / phys_fps)
  
  for i int = 0 to phys_fps*5 '' extend the curve
  
    '' horizontal and vertical scaling
    pset (i * 130 / phys_fps, h-1 - p * 99), col
    
    v += a
    v *= k
    p += v
  Next
  
end sub


sub main
  
  screenres w,h,32
  
  const golden = (1+5^.5)/2
  dim sng phys_fps = .25
  for i int = 1 to 10
    fps_plot phys_fps, rgb(rnd*255,rnd*255,rnd*255)
    phys_fps *= golden
  next
  
end sub

main

sleep
old version (works perfectly):

[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 Oct 13, 2020 21:55, edited 37 times in total.
vdecampo
Posts: 2992
Joined: Aug 07, 2007 23:20
Location: Maryland, USA
Contact:

Re: Terminal Velocity

Post by vdecampo »

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: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: terminal velocity sim

Post by dodicat »

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: 1641
Joined: Jun 04, 2005 9:51

Re: terminal velocity sim

Post by dafhi »

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: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: terminal velocity sim

Post by dodicat »

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: 1641
Joined: Jun 04, 2005 9:51

Re: terminal velocity

Post by dafhi »

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