GuoProfileFunction.py

import numpy as np
import pandas as pd
from .calcUstar import calcUstar

def GuoProfileFunction(data, uStar):
    '''
    This function is to fit the measured data with Guo's theoretical formula.
    The following steps are accroding to Junke Guo (2017)'s paper.
    Input : data = output from velocityProfile function
    '''
    n = -0.5 # indication for turbulent mixing at eta = etaC
    # step 1
    depth_array = np.linspace(0.0000001, data[1]-0.0000001, 100)
    p = np.poly1d(np.polyfit(data[0].reset_index()['depth[cm]'], data[0].reset_index()['U'], 4)) # 5th order polynomial fitting
    uMax = p(depth_array).max()
    yMax = depth_array[np.argmax(p(depth_array))]

    # step 2
    etaMax = 2*yMax/data[1]
    lambda_ = np.sqrt(data[1]/yMax-1)

    # step 3
    numAlpha = 1 + lambda_**(0.33) + lambda_**(0.67)
    denAlpha = lambda_*(1+lambda_**(0.33))
    alpha = numAlpha/denAlpha

    # step 4
    etaC = 2/(1+lambda_**(n))

    kappa = 0.39
    term1 = np.log(etaMax/etaC)
    term2 = lambda_*np.log((2-etaMax)/(2-etaC))
    term3 = (1+lambda_)/2*np.log(1+alpha*(1-etaMax/etaC)**2)
    term4 = (1-lambda_**(1+n))*np.sqrt(alpha)*np.arctan(np.sqrt(alpha))*(1-etaMax/etaC)
    RHS = (term1+term2-term3-term4)/kappa

    uStar = uStar # obtained by method in Nezu's book
    LHS = uMax/uStar
    uC = (LHS - RHS)*uStar

    # step 5
    etaArray = 2*depth_array/data[1]
    term1 = np.log(etaArray/etaC)
    term2 = lambda_*np.log((2-etaArray)/(2-etaC))
    term3 = (1+lambda_)/2*np.log(1+alpha*(1-etaArray/etaC)**2)
    term4 = (1-lambda_**(1+n))*np.sqrt(alpha)*np.arctan(np.sqrt(alpha))*(1-etaArray/etaC)

    phiArray = (term1+term2-term3-term4)/kappa

    return (phiArray+uC/uStar)*uStar, depth_array