diff --git a/data/simulations/lif-mu011-D0010-chi2.npz b/data/simulations/lif-mu011-D0010-chi2.npz
index 2642cf9..67b74f2 100644
Binary files a/data/simulations/lif-mu011-D0010-chi2.npz and b/data/simulations/lif-mu011-D0010-chi2.npz differ
diff --git a/lifsuscept.py b/lifsuscept.py
index b9b1362..409bfe1 100644
--- a/lifsuscept.py
+++ b/lifsuscept.py
@@ -55,18 +55,25 @@ def susceptibility2(omega1, omega2, chi1_1, chi1_2, r0, mu, D, tau_ref, vR, vT):
     return a1 * b1/c + a2*a3*b2/c - a2*a3*b3/c
 
 
-def susceptibilities(frange, mu, D, tau_ref, vR, vT):
-    print('compute LIF susceptibilites:')
+def susceptibilities(frange1, frange2, mu, D, tau_ref, vR, vT):
+    print(f'compute LIF susceptibilites for mu={mu:4.1f} and D={D:g}:')
+    print(f'    mean firing rate ...')
     r0 = firingrate(mu, D, tau_ref, vR, vT)
-    chi1_data = np.zeros(len(frange), dtype=complex)
-    chi2_data = np.zeros((len(frange), len(frange)), dtype=complex)
-    for f2 in range(len(frange)):
-        print(f'    step {f2 + 1:4d} of {len(frange):4d}')
-        omega2 = 2.0*np.pi*frange[f2]
+    # chi1:
+    print(f'    chi1 ...')
+    chi1_data = np.zeros(len(frange1), dtype=complex)
+    for f2 in range(len(frange1)):
+        omega2 = 2.0*np.pi*frange1[f2]
         chi1_2 = susceptibility1(omega2, r0, mu, D, tau_ref, vR, vT)
         chi1_data[f2] = chi1_2
-        for f1 in range(len(frange)):
-            omega1 = 2.0*np.pi*frange[f1]
+    # chi2:
+    chi2_data = np.zeros((len(frange2), len(frange2)), dtype=complex)
+    for f2 in range(len(frange2)):
+        print(f'    chi2 step {f2 + 1:4d} of {len(frange2):4d}')
+        omega2 = 2.0*np.pi*frange2[f2]
+        chi1_2 = susceptibility1(omega2, r0, mu, D, tau_ref, vR, vT)
+        for f1 in range(len(frange2)):
+            omega1 = 2.0*np.pi*frange2[f1]
             chi1_1 = susceptibility1(omega1, r0, mu, D, tau_ref, vR, vT)
             chi2 = susceptibility2(omega1, omega2, chi1_1, chi1_2, r0, mu, D, tau_ref, vR, vT)
             chi2_data[f2, f1] = chi2
@@ -76,19 +83,26 @@ def susceptibilities(frange, mu, D, tau_ref, vR, vT):
 def load_lifdata(mu, D, vT=1, vR=0, tau_ref=0):
     file_path = sims_path / f'lif-mu{10*mu:03.0f}-D{10000*D:04.0f}-chi2.npz'
     if not file_path.exists():
-        freqs = np.linspace(0.01, 1.0, 200)
-        r0, chi1, chi2 = susceptibilities(freqs, mu, D, tau_ref, vR, vT)
-        np.savez(file_path, mu=mu, D=D, vT=vT, vR=vR, tau_mem=1, tau_ref=tau_ref,
-                 r0=r0, freqs=freqs, chi1=chi1, chi2=chi2)
+        freqs1 = np.linspace(0.01, 1.0, 500)
+        freqs2 = np.linspace(0.01, 1.0, 200)
+        r0, chi1, chi2 = susceptibilities(freqs1, freqs2, mu, D,
+                                          tau_ref, vR, vT)
+        np.savez(file_path, mu=mu, D=D, vT=vT, vR=vR,
+                 tau_mem=1, tau_ref=tau_ref, r0=r0,
+                 freqs1=freqs1, chi1=chi1, freqs2=freqs2, chi2=chi2)
     data = np.load(file_path)
     r0 = float(data['r0'])
-    freqs = data['freqs']
+    freqs1 = data['freqs1']
     chi1 = data['chi1']
+    freqs2 = data['freqs2']
     chi2 = data['chi2']
-    return r0, freqs, chi1, chi2
+    print(f'LIF with mu={mu:4.1f} and D={D:g}')
+    return r0, freqs1, chi1, freqs2, chi2
 
 
 def plot_gain(ax, s, r0, freqs, chi1):
+    ax.axvline(r0, **s.lsGrid)
+    ax.axvline(2*r0, **s.lsGrid)
     ax.plot(freqs, np.abs(chi1), **s.lsM1)
     ax.set_xlabel('$f$')
     ax.set_ylabel('$|\\chi_1(f)|$', labelpad=6)
@@ -96,6 +110,8 @@ def plot_gain(ax, s, r0, freqs, chi1):
     ax.set_ylim(0, 14)
     ax.set_xticks_delta(0.2)
     ax.set_yticks_delta(3)
+    ax.text(r0, 14.2, '$r$', ha='center')
+    ax.text(2*r0, 14.2, '$2r$', ha='center')
 
 
 def plot_chi2(ax, s, r0, freqs, chi2):
@@ -130,15 +146,16 @@ if __name__ == '__main__':
     mu = 1.1
     D = 0.001
 
-    r0, freqs, chi1, chi2 = load_lifdata(mu, D)
+    r0, freqs1, chi1, freqs2, chi2 = load_lifdata(mu, D)
 
     s = plot_style()
     plt.rcParams['axes.labelpad'] = 2
-    fig, (axg, axc) = plt.subplots(1, 2, cmsize=(s.plot_width, 0.38*s.plot_width))
-    fig.subplots_adjust(leftm=8, rightm=8.5, topm=1, bottomm=3.5, wspace=0.4)
+    fig, (axg, axc) = plt.subplots(1, 2,
+                                   cmsize=(s.plot_width, 0.38*s.plot_width))
+    fig.subplots_adjust(leftm=8, rightm=8.5, topm=1.5, bottomm=3.5, wspace=0.4)
     fig.set_align(autox=False)
-    plot_gain(axg, s, r0, freqs, chi1)
-    plot_chi2(axc, s, r0, freqs, chi2)
+    plot_gain(axg, s, r0, freqs1, chi1)
+    plot_chi2(axc, s, r0, freqs2, chi2)
     fig.tag()
     fig.savefig()
     print()