import numpy as np
from scipy.signal import butter,freqs
from math import pi
f=np.geomspace(8,200,1000);w=2*pi*f;rho=1.18;c=343;N=4;Sd=.0504;Re=3.9;Le=.0035;Fs=22;Qes=.48;Qms=4.08;Vas=.0804;Vb=.35;Fb=20;S=.03
Cms=Vas/(rho*c*c*Sd*Sd);Mms=1/((2*pi*Fs)**2*Cms);Rms=2*pi*Fs*Mms/Qms;Bl=np.sqrt(2*pi*Fs*Mms*Re/Qes)
Cab=Vb/(rho*c*c);Map=1/((2*pi*Fb)**2*Cab);Rp=2*pi*Fb*Map/20;Rl=7/(2*pi*Fb*Cab)
Zp=Rp+1j*w*Map;Za=1/(1j*w*Cab+1/Zp+1/Rl);Ze=Re+1j*w*Le;Zm=Rms+1j*w*Mms+1/(1j*w*Cms)
b,a=butter(4,2*pi*17,btype='high',analog=True);hp=freqs(b,a,w)[1]
def sim(V=44.72136,HP=True,closed=False):
 za=1/(1j*w*Cab) if closed else Za
 v=Bl*(V/2)*(hp if HP else 1)/(Ze*(Zm+N*Sd*Sd*za)+Bl*Bl)
 up=np.zeros_like(v) if closed else N*Sd*v*za/Zp
 U=N*Sd*v-up
 spl=20*np.log10(np.abs(1j*w*rho*U/(2*pi))/2e-5)
 return spl,np.sqrt(2)*abs(v)/w*1000,np.sqrt(2)*abs(up)/S
if __name__=='__main__':
 for hz in [15,20,25,30,40,60,80]:
  k=np.argmin(abs(f-hz)); print(hz,[(round(z[k],2)) for z in sim()])
 print('xmax',max(sim()[1]),'port',max(sim()[2]))