straight_undular_tanh_ramp_continuous.ipynb

Kernel: Python 3 (ipykernel)

import sys
import numpy as np
import matplotlib.pyplot as plt

from math import sqrt, cosh, tanh, log, floor, acosh, exp
from datetime import date

plt.style.use('thesis')

today_str = date.today().strftime("%b-%d-%Y")

Define some constant values

• $g$ - gravity (determines units)

• $h$ - still water depth

• $a$ - target wave amplitude

• $L$ - length of the push

• P997 - timestep between paddle moves in ms (typically 10 to 20)

g = 9810
h = 50
a = 5
stroke = 400
ramp_time = 0.1
P997 = 10

Define functions for paddle movement

starting from physical quantites, we can define some useful paramaters: $a$ the wave amplitude, $h$ the still water depth, $L$ the stroke distance, $dt$ the ramp up time

then we can calculate the paddle velocity

$$\frac{h_b}{h_0} = \frac{1}{2} \left( \sqrt{1 + 8 F_b^2} - 1 \right)$$$$F_b = \frac{1}{\sqrt{2}} \sqrt{\frac{h_b}{h_0}} \sqrt{\frac{h_b}{h_0}+1}$$$$\frac{u_b}{c_0} = F_b \frac{\sqrt{1+ 8 F_b^2}-3}{\sqrt{1+8F_b^2}-1}$$

and the total time for the movement

$$t_f = \frac{dt}{6} \left( \cosh^{-1} \left( \exp \left( \frac{6 L}{dt u_0} + \log(\cosh(3)) \right) \right) + 3 \right)$$$$\tau = \frac{6}{dt} t$$$$\tau_f = \frac{6}{dt} t_f$$

such that the position and velocity may be determined for the entire stroke by

$$u(t) = \frac{u_0}{2}\left(\tanh\left(\tau-3\right)-\tanh\left(\left(\tau-\tau_{f}\right)+3\right)\right)$$$$x(t) = \frac{u_0 dt}{12}\left(\log\left(\cosh\left(\tau-3\right)\right)-\log\left(\cosh\left(\tau-\tau_{f}+3\right)\right)+\log\left(\cosh\left(-\tau_{f}+3\right)\right)-\log\left(\cosh\left(-3\right)\right)\right)$$

def calc_tauf(v0,stroke,ramp_time):
    ''' returns tauf for the given paramaters '''
    from math import exp,log,cosh,acosh
    
    L = 6*stroke/(ramp_time*v0)
    x = L + log(cosh(3))
    # easy to overflow exp(x) our of floating point range
    # but acosh(exp(x)) can be evaulated with a taylor series
    # in that case take the taylor series at x to infinity
    try:
        tauf = 3 + acosh(exp(x))
    except OverflowError:
        tauf = 3 + log(2) + x # - exp(-2x)/4 - ...
    
    #return 6*(stroke/(ramp_time*v0) + 1)
    return tauf
    

def tophat_vel(t, v0, h, stroke, ramp_time):
    ''' returns the velocity of the paddle at time t '''
    from math import tanh
    
    tau  = 6*t/ramp_time
    tauf = calc_tauf(v0,stroke,ramp_time)
    
    return (v0/2)*(tanh(tau - 3) - tanh(tau - tauf + 3))


def tophat_pos(t, v0, h, stroke, ramp_time):
    ''' returns the position of the paddle at time t'''
    from math import cosh,log,acosh,exp
    
    tau  = 6*t/ramp_time
    tauf = calc_tauf(v0,stroke,ramp_time)
    
    # term1 = log(cosh(tau-3))
    # term2 = log(cosh(tau-tauf+3))
    
    try:
        term1 = log(cosh(tau-3))
    except OverflowError:
        term1 = abs(tau - 3) - log(2)
    
    try:
        term2 = log(cosh(tau-tauf+3))
    except OverflowError:
        term2 = abs(tau - tauf + 3) - log(2)
        
    try:
        term3 = log(cosh(-tauf+3))
    except OverflowError:
        term3 = abs(-tauf + 3) - log(2)
        
    try:
        val = term1 - term2 + term3 - log(cosh(-3))
    except:
        print(term1,term2,tauf)
        raise
    
    return (ramp_time*v0/12)*val

Determine the total time for the paddle movement

ah = (a+h)/h
Fb = np.sqrt(ah)*np.sqrt(ah+1)/np.sqrt(2)
Fb2 = Fb**2
v0 = np.sqrt(g*h)*Fb*(np.sqrt(1+8*Fb2)-3)/(np.sqrt(1+8*Fb2)-1)

print(Fb,v0/sqrt(g*h))

tauf = calc_tauf(v0,stroke,ramp_time)
tfinal = tauf*ramp_time/6

dt = P997/1000
nsteps = int(tfinal/dt)+1
time = np.linspace(0,nsteps*dt,nsteps)

1.074709263010234 0.09770084209183953

print(tfinal)
print(tfinal*np.sqrt(g/h))

5.945817376981531
83.28390257317153

Generate wave form

ppmac = np.empty(time.shape)
vpmac = np.empty(time.shape)

for idx,t in enumerate(time):
    ppmac[idx] = tophat_pos(t, v0, h, stroke, ramp_time)
    vpmac[idx] = tophat_vel(t, v0, h, stroke, ramp_time)
    

# over ride small values at beginning and end
ppmac[0] = 0.
vpmac[0] = 0.
ppmac[-1] = stroke
vpmac[-1] = 0.

Check wave form

fig, axes = plt.subplots(nrows=1,ncols=2)

ax = axes[0]

ax.plot(time,ppmac/10,'k')
ax.set(xlabel='time (s)',ylabel='position (cm)')

ax = axes[1]

ax.plot(time,vpmac/10,'k')
ax.set(xlabel='time (s)',ylabel='velocity (cm/s)')

plt.tight_layout()

plt.savefig(f'undular_tanh_paddle_motions_h={h}_a={a}_L={stroke}_tau={ramp_time:.2f}.png')

c0 = np.sqrt(g*h)
t0 = h/c0

fig, axes = plt.subplots(nrows=1,ncols=2)

ax = axes[0]

ax.plot(time/t0,ppmac/h,'k')
ax.set(xlabel=r'$t/t_0$',ylabel=r'$\xi/h_0$')

ax = axes[1]

ax.plot(time/t0,vpmac/c0,'k')
ax.set(xlabel=r'$t/t_0$',ylabel=r'$(\partial \xi/ \partial t)/c_0$')

plt.tight_layout()

plt.savefig(f'undular_tanh_paddle_motions_h={h}_a={a}_L={stroke}_tau={ramp_time:.2f}_nondim.png')

xlabels = ('time (s)','time (s)','time (s)','time (s)','paddle position (mm)','paddle position (mm)')
ylabels = ('paddle position (mm)','paddle position (mm)','paddle velocity (mm/s)','paddle velocity (mm/s)','paddle velocity (mm/s)','paddle velocity (mm/s)')
titles  = ('paddle position vs time','paddle position vs time','paddle velocity vs time','paddle velocity vs time','paddle velocity vs paddle position','paddle velocity vs paddle position')
abscissa= (time ,time ,time ,time ,ppmac,ppmac)
ordinate= (ppmac,ppmac,vpmac,vpmac,vpmac,vpmac)
ranges  = ((0,30),(-30,-1),(0,30),(-30,-1),(0,30),(-30,-1))

for xlabel, ylabel, title, absc, ordi, rng in zip(xlabels,ylabels,titles,abscissa,ordinate,ranges):
    fig, ax = plt.subplots(figsize=(8,6))
    ax.plot(absc[rng[0]:rng[1]],ordi[rng[0]:rng[1]],'.')
    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)
    ax.set_title(title)
    plt.show()
    plt.close(fig)

xlabels = ('time (s)','time (s)','paddle position (mm)')
ylabels = ('paddle position (mm)','paddle velocity (mm/s)','paddle velocity (mm/s)')
titles  = ('paddle position vs time','paddle velocity vs time','paddle velocity vs paddle position')
abscissa= (time,time,ppmac)
ordinate= (ppmac,vpmac,vpmac)

for xlabel, ylabel, title, absc, ordi in zip(xlabels,ylabels,titles,abscissa,ordinate):
    fig, ax = plt.subplots(figsize=(8,6))
    ax.plot(absc[rng[0]:rng[1]],ordi[rng[0]:rng[1]],'.')
    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)
    ax.set_title(title)
    plt.show()
    plt.close(fig)

Check some paddle motion properties

print('the max paddle stroke is:',ppmac.max(),'mm')
print('the max paddle velocity is:',vpmac.max(),'mm/s')

print('the initial paddle velocity is:',ppmac[0],'mm')
print('the final paddle velocity is:',ppmac[-1],'mm')

print('the initial paddle velocity is:',vpmac[0],'mm/s')
print('the final paddle velocity is:',vpmac[-1],'mm/s')

the max paddle stroke is: 400.0 mm
the max paddle velocity is: 68.42547372540042 mm/s
the initial paddle velocity is: 0.0 mm
the final paddle velocity is: 400.0 mm
the initial paddle velocity is: 0.0 mm/s
the final paddle velocity is: 0.0 mm/s

Finally, write out the motion control programs

the motion control program for paddles 1 through 8

print(f"undular_h{h}_a{a}_A")
print(f"undular_h{h}_a{a}_B")

undular_h50_a5_A
undular_h50_a5_B

original_stdout = sys.stdout

with open("tryA.pmc", "w") as text_file:
    sys.stdout = text_file
    print("; WAVE TYPE - undular bore with tanh shape ramp up velocity")
    print(";              using continuous integration of tanh")
    print(f"; GRAVITATION CONSTANT - {g}")
    print(f"; WATER DEPTH - {h}")
    print(f"; TARGET AMPLITUDE - {a}")
    print(f"; TARGET STROKE LENGTH - {stroke}")
    print(f"; RAMP UP TIME - {ramp_time}")
    print("; FILE A")
    print(f"; GENERATION DATE - {today_str}")
    print(";")
    print("CLOSE")
    print("DELETE GATHER")
    print("DELETE TRACE")
    print("; M7024 = 1 turns off the TTL from output 1")
    print("M7024=1")
    print("; assign motors and scaling to coordinate system")
    print("&1 A")
    print("; in this case, all motors have the same motion, assign all to X")
    print("; to prevent motor motion, simply comment out the assigment",)
    print("; 1 inch is 50,000 counts therefore 1mm is 1968.50393700787 counts")
    print("; the units of X,Y,Z,U,V,W,A,B depend on the scaling defined here")
    print(f"#1->{50000/25.4:.6f}X")
    print(f"#2->{50000/25.4:.6f}X")
    print(f"#3->{50000/25.4:.6f}X")
    print(f"#4->{50000/25.4:.6f}X")
    print(f"#5->{50000/25.4:.6f}X")
    print(f"#6->{50000/25.4:.6f}X")
    print(f"#7->{50000/25.4:.6f}X")
    print(f"#8->{50000/25.4:.6f}X")
    print("; redefine motion program 2 ")
    print("OPEN PROG 2 CLEAR")
    print(f"P997={P997}")
    print("HOMEZ1..8")
    print("ABS")
    print("LINEAR")
    print("TS0")
    print("TA1000")
    print("TM5000")
    print("DWELL5000")
    print("PVT(P997)")
    print("; M7024=0 will turn on the TTL signal from output 1")
    print("M7024=0")
    print("; define motion like X(position):(velocity)")
    for pos, vel in zip(ppmac,vpmac):
        print(f"X({pos:.6f}):({vel:.6f})")
        
    print("DWELL15000 ; dwell time in ms")
    print("; M7024 = 1 turns off the TTL from output 1")
    print("M7024=1")
    print("LINEAR")
    print("TA100")
    print("TS0")
    print("TM1000")
    print("; Comment below to prevent paddles moving back to zero after dwelling")
    print("X0")
    print("CLOSE")
    print("CLOSE")
    sys.stdout = original_stdout

the motion control program for paddles 9 through 16

with open("tryB.pmc", "w") as text_file:
    sys.stdout = text_file
    print("; WAVE TYPE - undular bore with tanh shape ramp up velocity")
    print(";              using continuous integration of tanh")
    print(f"; GRAVITATION CONSTANT - {g}")
    print(f"; WATER DEPTH - {h}")
    print(f"; TARGET AMPLITUDE - {a}")
    print(f"; TARGET STROKE LENGTH - {stroke}")
    print(f"; RAMP UP TIME - {ramp_time}")
    print("; FILE B")
    print(f"; GENERATION DATE - {today_str}")
    print(";")
    print("CLOSE")
    print("DELETE GATHER")
    print("DELETE TRACE")
    print("; assign motors and scaling to coordinate system")
    print("&2 A")
    print("; in this case, all motors have the same motion, assign all to X")
    print("; to prevent motor motion, simply comment out the assigment")
    print("; 1 inch is 50,000 counts therefore 1mm is 1968.50393700787 counts")
    print("; the units of X,Y,Z,U,V,W,A,B depend on the scaling defined here")
    print(f"#9->{50000/25.4:.6f}X")
    print(f"#10->{50000/25.4:.6f}X")
    print(f"#11->{50000/25.4:.6f}X")
    print(f"#12->{50000/25.4:.6f}X")
    print(f"#13->{50000/25.4:.6f}X")
    print(f"#14->{50000/25.4:.6f}X")
    print(f"#15->{50000/25.4:.6f}X")
    print(f"#16->{50000/25.4:.6f}X")
    print("; redefine motion program 3 ")
    print("OPEN PROG 3 CLEAR")
    print(f"P997={P997}")
    print("HOMEZ9..16")
    print("ABS")
    print("LINEAR")
    print("TS0")
    print("TA1000")
    print("TM5000")
    print("DWELL5000")
    print("PVT(P997)")
    print("; define motion like X(position):(velocity)")
    for pos, vel in zip(ppmac,vpmac):
        print(f"X({pos:.6f}):({vel:.6f})")
        
    print("DWELL15000 ; dwell time in ms")
    print("LINEAR")
    print("TA100")
    print("TS0")
    print("TM1000")
    print("; Comment below to prevent paddles moving back to zero after dwelling")
    print("X0")
    print("CLOSE")
    print("CLOSE")
    sys.stdout = original_stdout