這個程式是之前的畢設，用Matlab寫的。學Python的時候用Python又寫了一遍。

估計沒別人看，所以就不詳細說明了。簡單來說就是淺水海域下，對點聲源產生噪聲的模擬，計算的是其修正係數delta=SL-RNL。具體可以去iso-17208看標準。需要解釋可以留言，對Matlab感興趣也可以留言。

Python程式碼如下：

import math
import pylab

print('如果不輸入值，則自動設為預設值')
todefault = input('是否直接使用所有預設值？（Y/N)')


def inputvalue(word, default='', de=todefault)
 
:
    if de == 'Y':
        return default
    else:
        input_value = input(word)
        if input_value == '':
            return default
        else:
            return input_value


def roughness(f, aoi):
    AoG = 3.1416 / 2 - aoi
    a = (0.006 + 0.0102 * VoW) ** (0-1)
    v31 = math.sin(AoG) * (1
 
 - (math.exp(-a * (AoG ** 2) / 4)) / (AoG * ((3.1416 * a) ** 0.5)))
    v32 = math.sin(AoG) / 2
    v3 = max(v31, v32)
    SL = -10 * math.log10(1 - v3) - 20 * math.log10(0.3 + 0.7 / (1 + 6 * (10 ** -11) * (VoW ** 4) * (f ** 2)))
    CoR = 10 ** (-SL / 20)
    return CoR


def pluslist(a,b):
    try:
        t = []
        for
 
 i in range(len(a)):
            t.append(a[i]+b[i])
    except:
        t = a
    return t


Rho = 1024
c = 1500
NoF = 100
NoG = 25
P = 1
try:
    ds = float(inputvalue('聲源深度（單位：米，預設4）：', '4'))
    H = float(inputvalue('水深（單位：米，預設60）：', '60'))
    dCPA = float(inputvalue('測量距離（單位：米，預設200）：', '200'))
    NoH = int(inputvalue('水聽器數量（預設3）：', '3'))
    HTop = float(inputvalue('頂部水聽器深度（單位：米，預設15）：', '15'))
    DeltaH = float(inputvalue('水聽器間距（單位：米，預設20）：', '20'))
    VoW = float(inputvalue('風速（單位：節，預設0）：', '0'))
    material = inputvalue('海底材質（預設細沙）：', '細沙')
    if material == '細沙':
        RhoSB = 1941
        cSB = 1749
        Att = 0.9
    else:
        print('未收錄該材料,請輸入引數：')
        RhoSB = float(input('密度(單位kg/m^3)：'))
        cSB = float(input('聲速(僅支援>1500m/s的材料）：'))
        Att = float(input('衰減(單位dB/波長）：'))
except:
    print('非法輸入，請重新執行')
    exit()
CGA = 3.1416/2 - math.asin(c / cSB)
f13 = [10**(19/20+x/10) for x in range(51)]
f = [10**(x/10/NoF+(1/NoF + 19)/20) for x in range(50*NoF)]
k = [2*3.1416*x/c for x in f]
dH = [-HTop-DeltaH*x for x in range(NoH)]
r = [((ds+dH[x])**2+dCPA**2)**0.5 for x in range(NoH)]


def reflection(aoi, c_water=c, r_water=Rho, c_seabed=cSB, r_seabed=RhoSB, attenuation=Att):
    alphaC = attenuation / 8.686 / c_seabed
    c_complex = 2 * 3.1416 / (2 * 3.1416 / cSB - 1j * alphaC)
    sinaoic = c_complex * math.sin(aoi) / c_water
    cosaoic = (1-sinaoic**2)**0.5
    z1 = c_water * r_water
    z2 = c_complex * r_seabed
    coefficient = (z2 * math.cos(aoi) - z1 * cosaoic) / (z2 * math.cos(aoi) + z1 * cosaoic)
    return coefficient


def pressure(noh, nog, fcount):
    pis = [-2*nog*H-ds, 2*nog*H+ds, 2*(nog+1)*H-ds, -2*(nog+1)*H+ds]
    rlocal = [((pis[x]-dH[noh])**2 + dCPA**2)**0.5 for x in range(4)]
    aoilocal = [math.asin(dCPA/x) for x in rlocal]
    prz = P*(((-reflection(aoilocal[0]) * roughness(f[fcount], aoilocal[0])) ** nog)
           * (2.71828**(0 - k[fcount] * rlocal[0] * 1j)) / rlocal[0] -
           ((-reflection(aoilocal[1]) * roughness(f[fcount], aoilocal[1])) ** nog) * roughness(f[fcount], aoilocal[1])
           * (2.71828**(0 - k[fcount] * rlocal[1] * 1j)) / rlocal[1] +
           ((-reflection(aoilocal[2]) * roughness(f[fcount], aoilocal[2])) ** (nog+1))
           * (2.71828**(0 - k[fcount] * rlocal[2] * 1j)) / rlocal[2] -
           ((-reflection(aoilocal[3]) * roughness(f[fcount], aoilocal[3])) ** (nog+1)) / roughness(f[fcount], aoilocal[3])
           * (2.71828**(0 - k[fcount] * rlocal[3] * 1j)) / rlocal[3])
    return prz


p4h = []
for i in range(NoH):
    sump = []
    for ii in range(len(f)):
        sump.append(((r[i] * abs(sum([pressure(i, x, ii) for x in range(NoG)])) / P)**2))
    p4h.append(sump)
ios = []
for i in range(NoH):
    poeh = []
    for ii in range(len(f13)-1):
        iii = ii * NoF
        poeh.append((sum(p4h[i][iii:iii + NoF])/NoF))
    ios = pluslist(poeh, ios)
x=[]
y=[]
for i in range(len(ios)):
    y.append(10*math.log10(ios[i]/NoH))
    x.append((f13[i]*f13[i+1])**0.5)
pylab.plot(x,y)
pylab.xlabel(u"frequency")
pylab.ylabel(u"correction factor")
pylab.title(u'delta-f curve')
pylab.semilogx()
file = open(r'./record.txt', 'w')
record=['x, y']
for i in range(len(x)):
    record.append(str(x[i])+' '+str(y[i]))
file.writelines(record)
pylab.show()