описать, что за что отвечает код

import numpy as np
 
def f_y(x):
    return a*np.exp(-x*x) + b*np.exp(-x*x/2)
    
    
def int_y(y, bi):
    return np.sqrt(np.fabs(f_y(y) - bi))
 
 
def integr(f, yt, bi, eps=1e-3):
    
    n = int(np.ceil(yt/(eps**(0.25))/4)*4) 
    n = max(n, 40)
    
    while True:
        h = yt/n 
        y = np.linspace(0, yt, n+1) 
        fy = f(y, bi)
 
        Ih1_s1 = fy[1::2].sum()
        Ih1_s2 = fy[2:-1:2].sum()
        
        Ih1 = h/3 * (fy[0] + fy[40] + 4*(Ih1_s1) + 2*(Ih1_s2))  
 
        Ih2_s1 = fy[2:-1:4].sum()
        Ih2_s2 = fy[4:-4:4].sum()
        
        Ih2 = 2*h/3 * (fy[0] + fy[40] + 4*(Ih2_s1) + 2*(Ih2_s2))
        Ih_check = abs(Ih1-Ih2)/15 
        
        if Ih_check < eps:
            return Ih2
            
        n <<= 1
    
    
def fun(b, eps=1e-5):
    xa, xb = 0, 3
    xm = (xa+xb)/2 
    while abs(f_y(xm) - b) > eps:
            if f_y(xm) > b:
                xa = xm 
            else:
                xb = xm 
            xm = (xa+xb)/2
    return xm
 
 
def fun_V(bi, vmax=40, eps=.1, k=100):
    yt = fun(bi)
    IntF = integr(int_y, yt, bi)
    v = np.pi*(3/4+i)/IntF
    hb = (bmb-bma)/(10*m)
    sv = np.sign(v-V)
    
    while abs(v-V) > eps and k>0:
        bi -= hb*sv
        yt = fun(bi)
        IntF = integr(int_y, yt, bi)
        v = np.pi*(3/4+i)/IntF
        st = np.sign(v-V)
        if not (st+sv):
            sv = st 
            hb /= 10
        k -= 1
    #print(k)
    print(f'Vi = {v}')
    return yt, bi
    
    
a = .7
b = .3
bma = 0.05
bmb = 0.9
m = 5
bm = np.linspace(bma,bmb,m+1)
V = 20
ytm = [0]*(m+1)
for i in range(m+1):
    ytm[i], bm[i] = fun_V(bm[i])
    print(f'm = {i}  ytm[{i}] = {ytm[i]:.4f}  bm[{i}] = {bm[i]:.4f}')
    print()