'===============================================================
' wave-rider buoy, water pump. Numerical model. Written in FreeBasic
' https://www.freebasic.net/  ( compatible with MS VisualBasic )
'===============================================================
' constants needed for simulation
Const As Double TwoPi = 8 * Atn( 1 )    ' radians per wave cycle
Const As Double g = 9.8     ' acceleration due to gravity
Const As Double density = 1020.4081632653 ' kg/m3 of water = 1025 for sea water
' when density is 1e4/g = 1020.4081632653 the numbers come out neatly

'---------------------------------------------------------------
' model parameters
Dim As Double wave_amplit = 1.5 ' metre, amplitude = ( peak to peak ) / 2
Dim As Double wave_period = 10  ' seconds

' assume buoy is made tall enough to always apply a +/- buoyant force
'   some extra height is needed to accelerate the inertia of the system and buoy
Dim As Double buoy_area = 1     ' section in sq metre
Dim As Double buoy_mass = 1000  ' kg total, will be approx half buoyancy
Dim As Double buoy_lag  = 0.5   ' metre, height before it moves

' set force needed to move buoy to buoy_lag height difference
Dim As Double restraint = buoy_lag * buoy_area * density * g   ' newton

' simulation time
Dim As Double step_time = 0.001 ' seconds
Dim As Double end_time  = 10    ' seconds

'---------------------------------------------------------------
' first the theoretical energy balance
Dim As Double dh = 2 * ( wave_amplit - buoy_lag )
dh = dh * 2             ' differential, rise then fall, full wave rectifier
If dh < 0 Then dh = 0   ' insufficient wave height
dh = dh * restraint / wave_period
Print Using " Power =######.## watt"; dh

'---------------------------------------------------------------
' setup the screen for graphics
Dim As Integer sx, sy, sz
Screeninfo sx, sy, sz   ' get maximum dimensions available
sx /= 2 ' half size =
sy /= 2 '   quarter pixels
sz = 4  ' 16 colours
Screenres sx, sy, sz    ' make graphics window
' scale, mapping
Window ( -0.05, -wave_amplit * 1.02 )-( end_time + 0.05, wave_amplit * 1.02 )
' axes and grid
For y As Double = -wave_amplit To wave_amplit + 0.1 Step wave_amplit
    Line( 0, y )-( end_time, y ), 1 ' horizontal grid lines
Next y
For x As Double = 0 To end_time Step 0.25 * wave_period
    Line( x, -wave_amplit )-( x, +wave_amplit ), 1 ' vertical grid lines
Next x

'---------------------------------------------------------------
' state variables
Dim As Double run_time, theta   ' time and wave phase
Dim As Double wave_height       ' relative to sea level
Dim As Double buoy_height = wave_amplit - buoy_lag   ' start relative to buoy neutral waterline
Dim As Double diff_height, buoy_force   ' buoy displacement and force
Dim As Double temp, joules      ' work done = accumulate energy extracted
Dim As Double last

'---------------------------------------------------------------
' simulation in time
For run_time = 0 To end_time Step step_time
    theta = TwoPi * Frac( run_time / wave_period )  ' wave phase    
    wave_height = wave_amplit * Cos( theta )' new value for wave surface
    diff_height = wave_height - buoy_height
    If  Abs( diff_height ) > buoy_lag Then  ' move buoy to extract energy
        temp = buoy_height    ' last value
        If wave_height > buoy_height Then   ' movement of buoy
            buoy_height = wave_height - buoy_lag
        Else
            buoy_height = wave_height + buoy_lag
        End If
        temp = Abs( buoy_height - temp )    ' how far the buoy moves
        joules = joules + restraint * temp
        ' else, no change in buoy_height since restrained
    End If
    Pset( run_time, wave_height ), 11   ' plot (x,y) points for this time
    Pset( run_time, buoy_height ), 13
    Pset( run_time, wave_height - buoy_height ), 12
    Line( run_time - step_time, last )-( run_time, temp / step_time ), 14
    Last = temp / step_time
Next

'---------------------------------------------------------------
' place coloured legend on screen
Draw String( 0.3, -wave_amplit * 0.2 ), " Wave ", 11
Draw String( 0.8, -wave_amplit * 0.2 ), " Buoy ", 13
Draw String( 0.3, -wave_amplit * 0.3 ), " Difference ", 12
Draw String( 0.3, -wave_amplit * 0.4 ), "   Energy ", 14

Bsave "wave_power.bmp", 0   ' block save graphics screen to a file.bmp

Print
Locate 52
Color 7
Print Using " ######.## watt"; joules / wave_period

'===============================================================
Sleep
Print " Done "
'===============================================================