byte TP = 0b10101010; // Every other port receives the inverted signal

void setup() {

DDRC = 0b11111111; // Set all analog ports to be outputs

// Initialize Timer1

noInterrupts(); // Disable interrupts

TCCR1A = 0;

TCCR1B = 0;

TCNT1 = 0;

OCR1A = 200; // Set compare register (16MHz / 200 = 80kHz square wave -> 40kHz full wave)

TCCR1B |= (1 << WGM12); // CTC mode

TCCR1B |= (1 << CS10); // Set prescaler to 1 ==> no prescaling

TIMSK1 |= (1 << OCIE1A); // Enable compare timer interrupt

interrupts(); // Enable interrupts

}

ISR(TIMER1_COMPA_vect) {

PORTC = TP; // Send the value of TP to the outputs

TP = ~TP; // Invert TP for the next run

}

void loop() {

// Nothing left to do here :)

}

#Analysis script for numerical solution of differential equation

import numpy as np

import scipy as sp

from ipywidgets import interact

from scipy.integrate import odeint

import matplotlib.pyplot as plt

%matplotlib inline

# Import data

file_dir = './'

filename= '/home/nano/ultra.csv'

data = np.genfromtxt(filename,delimiter=',')

time_data = data[:,0]

z_data = data[:,1]

time_data_clean = time_data[np.isfinite(z_data)]

z_data_clean = z_data[np.isfinite(z_data)]

# plot raw data

plt.plot(time_data_clean,z_data_clean,'.')

plt.xlabel('time (sec)',fontsize=15)

plt.ylabel('fluid level (cm)',fontsize=15)

plt.show()

def pend(y, t, b, a, m, g, f):

    x, omega = y

    dydt = [omega, b/m*omega + a/m*x + g + f/m]

    return dydt

#h(mm), b(空氣阻力係數), a(超音波震盪係數), m(ug), g(重力加速度), f(空氣浮力)

def plot_osc(h = 13.3, b = -0.2, a = -11.5, m = 130, g = 9.8, f = -1.4, T = 0.7):

    N_data = np.arange(0,T,0.001)

    plt.title('Ultrasonic oscillations\nh=%2.2f, b=%2.2f, a=%2.2f, m=%2.2f, g=%2.2f, f=%2.2f\nfilename = %s'%(h,b,a,m,g,f,filename),fontsize=15)

    h=h/100

    b=b/100000

    a=a/10000

    m=m/1000000000

    g=g*100

    f=f/10000

    x0 = [h, 0.0]

    sol = sp.integrate.odeint(pend, x0, N_data, args = (b, a, m, g, f))

    plt.plot(time_data_clean,z_data_clean,'.')

    plt.xlabel('time (sec)',fontsize=15)

    plt.ylabel('fluid level (cm)',fontsize=15)

    plt.plot(N_data, sol[:, 0], 'b', label='y(t)')

    plt.legend(loc='best')

    plt.xlabel('t')

    plt.grid()

    plt.show()

#plot_osc()

interact(plot_osc, h = (0.0, 50.0,0.01), b = (-0.5, 0.00,0.01), a = (-15,0,0.01), m = (0.00,150.00,0.01), g = (9.8,9.8,0.01), f = (-2,0.0,0.01), T = (0.0, 1.5, 0.1));