From 94a169701908fedca92f128a9fda327aaa6c7cf3 Mon Sep 17 00:00:00 2001
From: xaver <xaver.roos@t-online.de>
Date: Sat, 14 Nov 2020 19:03:48 +0100
Subject: [PATCH] 14.11

---
 apteronotus_code/avg_test.py                  |  20 ++
 .../figure_apteronotus_gain_plot.py           |   9 +
 .../figure_apteronotus_gain_plot_5Hz.py       | 114 +++++++++
 ...igure_apteronotus_gaincurve_cutofff_tau.py |  10 +-
 .../figure_apteronotus_jar_filter_fit.py      |   2 +-
 .../figure_apteronotus_jar_plot.py            | 224 +++++++++++++++---
 .../figure_apteronotus_rms_gaincurve.py       |   6 +-
 apteronotus_code/sin_all.py                   | 182 +++++++-------
 apteronotus_code/sin_all_model.py             | 142 ++++++-----
 eigenmannia_code/eigenmannia_jar_stacked.py   |  22 +-
 eigenmannia_code/eigenmannia_jar_subplot.py   |   7 +-
 eigenmannia_code/figure_eigen_gain_plot.py    | 128 ++++++++++
 eigenmannia_code/figure_eigen_jar_plot.py     | 168 +++++++++++++
 13 files changed, 819 insertions(+), 215 deletions(-)
 create mode 100644 apteronotus_code/avg_test.py
 create mode 100644 apteronotus_code/figure_apteronotus_gain_plot_5Hz.py
 create mode 100644 eigenmannia_code/figure_eigen_gain_plot.py
 create mode 100644 eigenmannia_code/figure_eigen_jar_plot.py

diff --git a/apteronotus_code/avg_test.py b/apteronotus_code/avg_test.py
new file mode 100644
index 0000000..08727b2
--- /dev/null
+++ b/apteronotus_code/avg_test.py
@@ -0,0 +1,20 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pylab
+from IPython import embed
+from scipy.optimize import curve_fit
+from jar_functions import gain_curve_fit
+from jar_functions import avgNestedLists
+import matplotlib as mpl
+from matplotlib import cm
+
+ab = []
+a = [1, 1, None, 1]
+b = [2, 2, 2, 2]
+ab.append(a)
+ab.append(b)
+print(ab)
+print(np.mean(ab, axis = 0))
+#av = avgNestedLists(np.array(ab))
+#print(av)
+
diff --git a/apteronotus_code/figure_apteronotus_gain_plot.py b/apteronotus_code/figure_apteronotus_gain_plot.py
index 7dfe595..b77c1a3 100644
--- a/apteronotus_code/figure_apteronotus_gain_plot.py
+++ b/apteronotus_code/figure_apteronotus_gain_plot.py
@@ -30,6 +30,7 @@ gains = []
 taus = []
 f_cs = []
 predicts = []
+maxgains = []
 for ID in identifier:
     predict = []
 
@@ -47,10 +48,18 @@ for ID in identifier:
     f_cs.append(f_cutoff)
 
     # predict of gain
+    print('max gain:', np.max(gain))
+    maxgains.append(np.max(gain))
+
     for f in amf:
         G = np.max(gain) / np.sqrt(1 + (2 * ((np.pi * f * sinv[0]) ** 2)))
         predict.append(G)
     predicts.append(predict)
+print('absolute max gain:', np.max(maxgains))
+print('absolute min gain:', np.min(maxgains))
+print('absolute max f_c:', np.max(f_cs))
+print('absolute min f_c:', np.min(f_cs))
+
 
 sort = sorted(zip(f_cs, identifier))
 print(sort)
diff --git a/apteronotus_code/figure_apteronotus_gain_plot_5Hz.py b/apteronotus_code/figure_apteronotus_gain_plot_5Hz.py
new file mode 100644
index 0000000..c997137
--- /dev/null
+++ b/apteronotus_code/figure_apteronotus_gain_plot_5Hz.py
@@ -0,0 +1,114 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pylab
+from IPython import embed
+from scipy.optimize import curve_fit
+from matplotlib.mlab import specgram
+import os
+from jar_functions import gain_curve_fit
+
+plt.rcParams.update({'font.size': 12})
+
+identifier = [#'2018lepto1',
+              #'2018lepto4',
+              #'2018lepto5',
+              #'2018lepto76',
+              '2018lepto98',
+              #'2019lepto03',
+              #'2019lepto24',
+              #'2019lepto27',
+              #'2019lepto30',
+              #'2020lepto04',
+              #'2020lepto06',
+              '2020lepto16',
+              '2020lepto19',
+              '2020lepto20'
+              ]
+
+amfs = []
+gains = []
+taus = []
+f_cs = []
+predicts = []
+for ID in identifier:
+    predict = []
+
+    print(ID)
+    amf = np.load('5Hz_amf_%s.npy' %ID)
+    amfs.append(amf)
+    gain = np.load('5Hz_gain_%s.npy' %ID)
+    gains.append(gain)
+
+    sinv, sinc = curve_fit(gain_curve_fit, amf, gain, [2, 3])
+    #print('tau:', sinv[0])
+    taus.append(sinv[0])
+    f_cutoff = abs(1 / (2*np.pi*sinv[0]))
+    print('f_cutoff:', f_cutoff)
+    f_cs.append(f_cutoff)
+
+    # predict of gain
+    for f in amf:
+        G = np.max(gain) / np.sqrt(1 + (2 * ((np.pi * f * sinv[0]) ** 2)))
+        predict.append(G)
+    predicts.append(predict)
+
+sort = sorted(zip(f_cs, identifier))
+print(sort)
+# order of plotting: 2018lepto98, 2020lepto16, 2020lepto19, 2020lepto19, 2020lepto20
+# order of f_c: 2020lepto20, 2020lepto16, 2018lepto98, 2020lepto19
+
+fig = plt.figure(figsize=(8.27,11.69))
+ax0 = fig.add_subplot(221)
+fig.text(0.05, 0.5, 'gain [Hz/(mV/cm)]', ha='center', va='center', rotation='vertical')
+fig.text(0.5, 0.04, 'envelope frequency [Hz]', ha='center', va='center')
+
+ax0.set_xlim(0.0007, 1.5)
+ax0.set_ylim(0.001, 10)
+ax0.plot(amfs[0], gains[0],'o' , label = 'gain')
+ax0.plot(amfs[0], predicts[0], label = 'fit')
+ax0.axvline(x=f_cs[0], ymin=0, ymax=5, ls='-', alpha=0.5, label = 'cutoff frequency')
+ax0.set_xscale('log')
+ax0.set_yscale('log')
+ax0.axes.xaxis.set_ticklabels([])
+print('max[0]:', np.max(gain[0]))
+
+ax1 = fig.add_subplot(222)
+ax1.set_xlim(0.0007, 1.5)
+ax1.set_ylim(0.001, 10)
+ax1.plot(amfs[1], gains[1],'o' , label = 'gain')
+ax1.plot(amfs[1], predicts[1], label = 'fit')
+ax1.axvline(x=f_cs[1], ymin=0, ymax=5, ls='-', alpha=0.5, label = 'cutoff frequency')
+ax1.set_xscale('log')
+ax1.set_yscale('log')
+ax1.axes.yaxis.set_ticklabels([])
+ax1.axes.xaxis.set_ticklabels([])
+print('max[1]:', np.max(gain[1]))
+
+ax2 = fig.add_subplot(223)
+ax2.set_xlim(0.0007, 1.5)
+ax2.set_ylim(0.001, 10)
+ax2.plot(amfs[0], gains[0],'o' , label = 'gain')
+ax2.plot(amfs[0], predicts[0], label = 'fit')
+ax2.axvline(x=f_cs[0], ymin=0, ymax=5, ls='-', alpha=0.5, label = 'cutoff frequency')
+ax2.set_xscale('log')
+ax2.set_yscale('log')
+print('max[2]:', np.max(gain[2]))
+
+ax3 = fig.add_subplot(224)
+ax3.set_xlim(0.0007, 1.5)
+ax3.set_ylim(0.001, 10)
+ax3.plot(amfs[2], gains[2],'o' , label = 'gain')
+ax3.plot(amfs[2], predicts[2], label = 'fit')
+ax3.axvline(x=f_cs[2], ymin=0, ymax=5, ls='-', alpha=0.5, label = 'cutoff frequency')
+ax3.set_xscale('log')
+ax3.set_yscale('log')
+ax3.axes.yaxis.set_ticklabels([])
+print('max[3]:', np.max(gain[3]))
+
+#plt.legend(loc = 'lower left')
+
+plt.show()
+
+
+#np.save('f_c', f_c)
+#np.save('tau', tau)
\ No newline at end of file
diff --git a/apteronotus_code/figure_apteronotus_gaincurve_cutofff_tau.py b/apteronotus_code/figure_apteronotus_gaincurve_cutofff_tau.py
index b04b2dc..6a30387 100644
--- a/apteronotus_code/figure_apteronotus_gaincurve_cutofff_tau.py
+++ b/apteronotus_code/figure_apteronotus_gaincurve_cutofff_tau.py
@@ -7,7 +7,9 @@ from matplotlib.mlab import specgram
 import os
 from jar_functions import gain_curve_fit
 
-identifier = ['2018lepto98']
+plt.rcParams.update({'font.size': 12})
+
+identifier = ['2018lepto4']
 
 tau = []
 f_c = []
@@ -32,7 +34,7 @@ for ID in identifier:
         predict.append(G)
     print(np.max(gain))
 
-    fig = plt.figure()
+    fig = plt.figure(figsize=(8.27, 11.69/2))
     ax = fig.add_subplot()
     ax.plot(amf, gain,'o' , label = 'gain')
     ax.plot(amf, predict, label = 'fit')
@@ -40,9 +42,9 @@ for ID in identifier:
     ax.set_xscale('log')
     ax.set_yscale('log')
     ax.set_ylabel('gain [Hz/(mV/cm)]')
-    ax.set_xlabel('envelope_frequency [Hz]')
+    ax.set_xlabel('envelope frequency [Hz]')
     ax.set_title('gaincurve %s' %ID)
-    plt.legend(loc = 'lower left')
+    #plt.legend(loc = 'lower left')
     plt.show()
 
 
diff --git a/apteronotus_code/figure_apteronotus_jar_filter_fit.py b/apteronotus_code/figure_apteronotus_jar_filter_fit.py
index f41b8ea..b153fe3 100644
--- a/apteronotus_code/figure_apteronotus_jar_filter_fit.py
+++ b/apteronotus_code/figure_apteronotus_jar_filter_fit.py
@@ -82,7 +82,7 @@ for ident in identifier:
             # plt.show()
 
             # filter by running average
-            fig = plt.figure(figsize = (8,14))
+            fig = plt.figure(figsize = (8.27,11.69))
             fig.suptitle('JAR trace spectogram 2018lepto98:\n subtraction of mean and running average')
             ax = fig.add_subplot(211)
             ax.plot(time, jm, color = 'C0', label = '1)')
diff --git a/apteronotus_code/figure_apteronotus_jar_plot.py b/apteronotus_code/figure_apteronotus_jar_plot.py
index 4324f6e..af0e358 100644
--- a/apteronotus_code/figure_apteronotus_jar_plot.py
+++ b/apteronotus_code/figure_apteronotus_jar_plot.py
@@ -17,22 +17,37 @@ from jar_functions import gain_curve_fit
 
 def take_second(elem):      # function for taking the names out of files
     return elem[1]
-
-identifier = ['2018lepto1',
-              '2018lepto4',
-              '2018lepto5',
-              '2018lepto76',
+identifier = [#'2018lepto1',
+              #'2018lepto4',
+              #'2018lepto5',
+              #'2018lepto76',
+              '2018lepto98',
+              #'2019lepto03',
+              #'2019lepto24',
+              #'2019lepto27',
+              #'2019lepto30',
+              #'2020lepto04',
+              #'2020lepto06',
+              #'2020lepto16',
+              #'2020lepto19',
+              #'2020lepto20'
+              ]
+identifier_5Hz = [#'2018lepto1',
+              #'2018lepto4',
+              #'2018lepto5',
+              #'2018lepto76',
               '2018lepto98',
-              '2019lepto03',
-              '2019lepto24',
-              '2019lepto27',
-              '2019lepto30',
-              '2020lepto04',
-              '2020lepto06',
-              '2020lepto16',
-              '2020lepto19',
-              '2020lepto20'
+              #'2019lepto03',
+              #'2019lepto24',
+              #'2019lepto27',
+              #'2019lepto30',
+              #'2020lepto04',
+              #'2020lepto06',
+              #'2020lepto16',
+              #'2020lepto19',
+              #'2020lepto20'
               ]
+
 for ident in identifier:
 
     times = []
@@ -47,15 +62,15 @@ for ident in identifier:
 
     amf = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1]
 
-    data = sorted(np.load('%s files.npy' %ident), key = take_second)      # list with filenames in it
+    data = sorted(np.load('5Hz_%s files.npy' %ident), key = take_second)      # list with filenames in it
 
     for i, d in enumerate(data):
         dd = list(d)
         if dd[1] == '1' or dd[1] == '0.2' or dd[1] == '0.05' or dd[1] == '0.01' or dd[1] == '0.005' or dd[1] == '0.001':
-            jar = np.load('%s.npy' %dd)     # load data for every file name
+            jar = np.load('5Hz_%s.npy' %dd)     # load data for every file name
             jm = jar - np.mean(jar)         # low-pass filtering by subtracting mean
 
-            time = np.load('%s time.npy' %dd)       # time file
+            time = np.load('5Hz_%s time.npy' %dd)       # time file
             dt = time[1] - time[0]
 
             n = int(1/float(d[1])/dt)
@@ -73,7 +88,7 @@ for ident in identifier:
                 jms.append(jm)
                 times.append(time)
             else:
-                print('1:', dd)
+                #print('1:', dd)
                 amfreq1.append(dd[1])
                 jars1.append(jm - cutf)
                 jms1.append(jm)
@@ -81,51 +96,184 @@ for ident in identifier:
     if len(jars) != 6:
         continue
 
+    sample = 100000
+
     fig = plt.figure(figsize=(8.27,11.69))
     fig.suptitle('%s' %ident)
-    fig.text(0.06, 0.5, 'frequency [Hz]', ha='center', va='center', rotation='vertical')
+    fig.text(0.06, 0.5, 'fish frequency [Hz]', ha='center', va='center', rotation='vertical', color = 'C0')
+    fig.text(0.97, 0.5, 'stimulus amplitude [mV/cm]', ha='center', va='center', rotation='vertical', color = 'red')
     fig.text(0.5, 0.04, 'time [s]', ha='center', va='center')
 
     ax0 = fig.add_subplot(611)
-    ax0.plot(times[0], jms[0])
-    #ax0.plot(times[0], jars[0])
+    print('absolute frequency shift 0.001Hz:', np.max(jars[0]) - np.min(jars[0]))
+    ax0.plot(times[0], jars[0], zorder = 20)
+    #ax0.set_zorder(1)
     ax0.set_ylim(-12, 12)
-    #plt.text(-0.1, 1.05, "A)", fontweight=550, transform=ax0.transAxes)
+
+    lower0 = 0
+    upper0 = 2000
+    x0 = np.linspace(lower0, upper0, sample)
+    y0 = (sin_response(np.linspace(lower0, upper0, sample), 0.001, -np.pi/2, .35) + 0.5)
+    ax0_0 = ax0.twinx()
+    ax0_0.set_ylim(-0.2, 1.2)
+    ax0_0.plot(x0, y0, color = 'red', zorder = 1, alpha = 0.5)
+    #ax0_0.set_zorder(2)
 
     ax1 = fig.add_subplot(612)
-    ax1.plot(times[1], jms[1])
-    #ax1.plot(times[1], jars[1])
+    print('absolute frequency shift 0.005 Hz:', np.max(jars[1]) - np.min(jars[1]))
+    ax1.plot(times[1], jars[1])
     ax1.set_ylim(-12, 12)
-    #plt.text(-0.1, 1.05, "B)", fontweight=550, transform=ax1.transAxes)
+
+    lower1 = 0
+    upper1 = 400
+    x1 = np.linspace(lower1, upper1, sample)
+    y1 = (sin_response(np.linspace(lower1, upper1, sample), 0.005, -np.pi / 2, .35) + 0.5)
+    ax1_0 = ax1.twinx()
+    ax1_0.set_ylim(-0.2, 1.2)
+    ax1_0.plot(x1, y1, color='red', alpha = 0.5)
 
     ax2 = fig.add_subplot(613)
-    ax2.plot(times[2], jms[2])
-    #ax2.plot(times[2], jars[2])
+    print('absolute frequency shift 0.01 Hz:', np.max(jars[2]) - np.min(jars[2]))
+    ax2.plot(times[2], jars[2])
     ax2.set_ylim(-12, 12)
-    #plt.text(-0.1, 1.05, "C)", fontweight=550, transform=ax2.transAxes)
+
+    lower2 = 0
+    upper2 = 400
+    x2 = np.linspace(lower2, upper2, sample)
+    y2 = (sin_response(np.linspace(lower2, upper2, sample), 0.01, np.pi / 2, -0.35) + 0.5)
+    ax2_0 = ax2.twinx()
+    ax2_0.set_ylim(-0.2, 1.2)
+    ax2_0.plot(x2, y2, color='red', alpha = 0.5)
 
     ax3 = fig.add_subplot(614)
-    ax3.plot(times[3], jms[3])
-    #ax3.plot(times[3], jars[3])
+    print('absolute frequency shift 0.02 Hz:', np.max(jars[3]) - np.min(jars[3]))
+    ax3.plot(times[3], jars[3])
     ax3.set_ylim(-12, 12)
-    #plt.text(-0.1, 1.05, "D)", fontweight=550, transform=ax3.transAxes)
+
+    lower3 = 0
+    upper3 = 200
+    x3 = np.linspace(lower3, upper3, sample)
+    y3 = (sin_response(np.linspace(lower3, upper3, sample), 0.05, np.pi / 2, -0.35) + 0.5)
+    ax3_0 = ax3.twinx()
+    ax3_0.set_ylim(-0.2, 1.2)
+    ax3_0.plot(x3, y3, color='red', alpha = 0.5)
 
     ax4 = fig.add_subplot(615)
-    ax4.plot(times[4], jms[4])
-    #ax4.plot(times[4], jars[4])
+    print('absolute frequency shift 0.5 Hz:', np.max(jars[4]) - np.min(jars[4]))
+    ax4.plot(times[4], jars[4])
     ax4.set_ylim(-12, 12)
-   # plt.text(-0.1, 1.05, "E)", fontweight=550, transform=ax4.transAxes)
+
+    lower4 = 0
+    upper4 = 200
+    x4 = np.linspace(lower4, upper4, sample)
+    y4 = (sin_response(np.linspace(lower4, upper4, sample), 0.2, np.pi / 2, -0.35) + 0.5)
+    ax4_0 = ax4.twinx()
+    ax4_0.set_ylim(-0.2, 1.2)
+    ax4_0.plot(x4, y4, color='red', alpha = 0.5)
 
     ax5 = fig.add_subplot(616)
-    ax5.plot(times[5], jms[5])
-    #ax5.plot(times[5], jars[5])
+    print('absolute frequency shift 1 Hz:', np.max(jars[5]) - np.min(jars[5]))
+    ax5.plot(times[5], jars[5])
     ax5.set_ylim(-12, 12)
-    #plt.text(-0.1, 1.05, "F)", fontweight=550, transform=ax5.transAxes)
 
+    lower5 = 0
+    upper5 = 200
+    x5 = np.linspace(lower5, upper5, sample)
+    y5 = (sin_response(np.linspace(lower5, upper5, sample), 1, np.pi / 2, -0.35) + 0.5)
+    ax5_0 = ax5.twinx()
+    ax5_0.plot(x5, y5, color='red', lw = 0.5, alpha = 0.5)
+    ax5_0.set_ylim(-0.2, 1.2)
+    '''
+    
+for ident in identifier_5Hz:
+
+    times = []
+    jars = []
+    jms = []
+    amfreq = []
+
+    times1 = []
+    jars1 = []
+    jms1 = []
+    amfreq1 = []
+
+    amf = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1]
+    print('5Hz')
+    data = sorted(np.load('5Hz_%s files.npy' % ident), key=take_second)  # list with filenames in it
+
+    for i, d in enumerate(data):
+        dd = list(d)
+        if dd[1] == '1' or dd[1] == '0.2' or dd[1] == '0.05' or dd[1] == '0.01' or dd[1] == '0.005' or dd[1] == '0.001':
+            jar = np.load('5Hz_%s.npy' % dd)  # load data for every file name
+            jm = jar - np.mean(jar)  # low-pass filtering by subtracting mean
+
+            time = np.load('5Hz_%s time.npy' % dd)  # time file
+            dt = time[1] - time[0]
+
+            n = int(1 / float(d[1]) / dt)
+            cutf = mean_noise_cut(jm, n=n)
+            cutt = time
+            if dd[1] == '0.001':
+                amfreq1.append(dd[1])
+                jars1.append(jm - cutf)
+                jms1.append(jm)
+                times1.append(time)
+            if dd[1] not in amfreq:
+                #print(dd)
+                amfreq.append(dd[1])
+                jars.append(jm - cutf)
+                jms.append(jm)
+                times.append(time)
+            else:
+                #print('1:', dd)
+                amfreq1.append(dd[1])
+                jars1.append(jm - cutf)
+                jms1.append(jm)
+                times1.append(time)
+    if len(jars) != 6:
+        continue
+
+    ax6 = fig.add_subplot(622)
+    ax6.plot(times1[0], jars1[0])
+    print('5Hz_absolute frequency shift 0.001Hz:', np.max(jms1[0]) - np.min(jms1[0]))
+    ax6.set_ylim(-12, 12)
+    ax6.axes.set_yticklabels([])
+    plt.text(-0.05, 1.15, "B)", fontweight=550, transform=ax6.transAxes)
+
+    ax7 = fig.add_subplot(624)
+    ax7.plot(times1[1], jars1[1])
+    print('5Hz_absolute frequency shift 0.005Hz:', np.max(jms1[1]) - np.min(jms1[1]))
+    ax7.set_ylim(-12, 12)
+    ax7.axes.set_yticklabels([])
+
+    ax8 = fig.add_subplot(626)
+    ax8.plot(times1[2], jars1[2])
+    print('5Hz_absolute frequency shift 0.05Hz:', np.max(jms1[2]) - np.min(jms1[2]))
+    ax8.set_ylim(-12, 12)
+    ax8.axes.set_yticklabels([])
+
+    ax9 = fig.add_subplot(6,2,8)
+    ax9.plot(times1[3], jars1[3])
+    print('5Hz_absolute frequency shift 0.02Hz:', np.max(jms1[3]) - np.min(jms1[3]))
+    ax9.set_ylim(-12, 12)
+    ax9.axes.set_yticklabels([])
+
+    ax10 = fig.add_subplot(6,2,10)
+    ax10.plot(times1[4], jars1[4])
+    print('5Hz_absolute frequency shift 0.5Hz:', np.max(jms1[4]) - np.min(jms1[4]))
+    ax10.set_ylim(-12, 12)
+    ax10.axes.set_yticklabels([])
+
+    ax11 = fig.add_subplot(6,2,12)
+    ax11.plot(times1[5], jars1[5])
+    print('5Hz_absolute frequency shift 1Hz:', np.max(jms1[5]) - np.min(jms1[5]))
+    ax11.set_ylim(-12, 12)
+    ax11.axes.set_yticklabels([])
+    '''
     plt.subplots_adjust(left=0.125,
                         bottom=0.1,
                         right=0.9,
                         top=0.9,
-                        wspace=0.2,
+                        wspace=0.1,
                         hspace=0.35)
     plt.show()
\ No newline at end of file
diff --git a/apteronotus_code/figure_apteronotus_rms_gaincurve.py b/apteronotus_code/figure_apteronotus_rms_gaincurve.py
index bfa4d67..4487e69 100644
--- a/apteronotus_code/figure_apteronotus_rms_gaincurve.py
+++ b/apteronotus_code/figure_apteronotus_rms_gaincurve.py
@@ -19,7 +19,7 @@ plt.rcParams.update({'font.size': 12})
 def take_second(elem):      # function for taking the names out of files
     return elem[1]
 
-identifier = ['2020lepto19']
+identifier = ['2018lepto5']
 for ident in identifier:
 
     predict = []
@@ -124,10 +124,10 @@ for ident in identifier:
     # condition needed to be fulfilled: RMS < threshold or RMS < mean(RMS)
     idx_arr = (rootmeansquare_arr < threshold_arr) | (rootmeansquare_arr < np.mean(rootmeansquare_arr))
 
-    fig = plt.figure(figsize = (8,14))
+    fig = plt.figure(figsize = (8.27, 11.69))
     fig.suptitle('gaincurve and RMS %s' %ident)
     ax0 = fig.add_subplot(2, 1, 1)
-    ax0.plot(amfreq_arr[idx_arr], mgain_arr[idx_arr], 'o')
+    ax0.plot(amfreq_arr, mgain_arr, 'o')
     ax0.set_yscale('log')
     ax0.set_xscale('log')
     ax0.set_ylabel('gain [Hz/(mV/cm)]')
diff --git a/apteronotus_code/sin_all.py b/apteronotus_code/sin_all.py
index be66ef9..404b04a 100644
--- a/apteronotus_code/sin_all.py
+++ b/apteronotus_code/sin_all.py
@@ -7,125 +7,143 @@ from jar_functions import gain_curve_fit
 from jar_functions import avgNestedLists
 import matplotlib as mpl
 from matplotlib import cm
+import math
 
-#plt.rcParams.update({'font.size': 18})
+#plt.rcParams.update({'font.size': 16})
 
-identifier_uniform = ['2018lepto1',
+identifier = [#'2018lepto1',
               #'2018lepto4',
-              '2018lepto5',
-              '2018lepto76',
+              #'2018lepto5',
+              #'2018lepto76',
               '2018lepto98',
               #'2019lepto03',
-              '2019lepto24',
+              #'2019lepto24',
               #'2019lepto27',
               #'2019lepto30',
               #'2020lepto04',
-              '2020lepto06',
-              #'2020lepto16',
-              #'2020lepto19',
-              #'2020lepto20'
-              ]
-identifier = ['2018lepto1',
-              '2018lepto4',
-              '2018lepto5',
-              '2018lepto76',
-              '2018lepto98',
-              '2019lepto03',
-              '2019lepto24',
-              '2019lepto27',
-              '2019lepto30',
-              '2020lepto04',
-              '2020lepto06',
+              #'2020lepto06',
               '2020lepto16',
               '2020lepto19',
               '2020lepto20'
               ]
 
-custom_amf = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1]
-#colors = ['dimgray', 'dimgrey', 'gray', 'grey', 'darkgray', 'darkgrey', 'silver', 'lightgray', 'lightgrey', 'gainsboro', 'whitesmoke']
-colorss = ['g', 'b', 'r', 'y', 'c', 'm', 'k']
-all = []
-new_all = []
-for ident in identifier:
-    data = np.load('gain_%s.npy' %ident)
-    all.append(data)
-for ident in identifier_uniform:
-    data = np.load('gain_%s.npy' % ident)
-    new_all.append(data)
-
-av = avgNestedLists(all)
-new_av = avgNestedLists(new_all)
-
 fig = plt.figure(figsize=(8.27, 11.69/2))
 ax = fig.add_subplot(111)
-ax.plot(custom_amf, av, 'o', label = 'normal')
-#ax.plot(amf, new_av, 'o', label = 'uniform')
-
-sinv, sinc = curve_fit(gain_curve_fit, custom_amf, av, [2, 3])
-predict = []
-for f in custom_amf:
-    G = np.max(av) / np.sqrt(1 + (2 * ((np.pi * f * sinv[0]) ** 2)))
-    predict.append(G)
 
+custom_amf = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1]
 tau = []
+IDs = []
 f_c = []
 fit = []
 fit_amf = []
+
+all_gains = []
 for ID in identifier:
+
     print(ID)
-    amf = np.load('amf_%s.npy' %ID)
-    gain = np.load('gain_%s.npy' %ID)
+    IDs.append(ID)
+    gain_10 = np.zeros(10)
+    amf = np.load('5Hz_amf_%s.npy' % ID)
+    gain = np.load('5Hz_gain_%s.npy' % ID)
+    b = 0
+    for aa, a in enumerate(custom_amf):
+        if a in amf:
+            gain_10[aa] = gain[b]
+            b += 1
+        else:
+            gain_10[aa] = None
+    print(gain_10)
+    #print(amf)
+    all_gains.append(gain_10)
+
 
     sinv, sinc = curve_fit(gain_curve_fit, amf, gain)
-    #print('tau:', sinv[0])
+    print('tau:', sinv[0])
     tau.append(sinv[0])
     f_cutoff = abs(1 / (2*np.pi*sinv[0]))
     print('f_cutoff:', f_cutoff)
     f_c.append(f_cutoff)
     fit.append(gain_curve_fit(amf, *sinv))
     fit_amf.append(amf)
-    ax.axvline(x=f_cutoff, ymin=0, ymax=5, color='C0', ls='-', alpha=0.5, label='cutoff frequency')
-
+    ax.axvline(x=f_cutoff, ymin=0, ymax=5, color='C0', ls='-', alpha=0.5)
+
+f_c_ID = zip(ID, f_c)
+
+mean = []
+
+g0 = []
+g1 = []
+g2 = []
+g3 = []
+g4 = []
+g5 = []
+g6 = []
+g7 = []
+g8 = []
+g9 = []
+for g in all_gains:
+    if math.isnan(g[0]) is False:
+        g0.append(g[0])
+    if math.isnan(g[1]) is False:
+        g1.append(g[1])
+    if math.isnan(g[2]) is False:
+        g2.append(g[2])
+    if math.isnan(g[3]) is False:
+        g3.append(g[3])
+    if math.isnan(g[4]) is False:
+        g4.append(g[4])
+    if math.isnan(g[5]) is False:
+        g5.append(g[5])
+    if math.isnan(g[6]) is False:
+        g6.append(g[6])
+    if math.isnan(g[7]) is False:
+        g7.append(g[7])
+    if math.isnan(g[8]) is False:
+        g8.append(g[8])
+    if math.isnan(g[9]) is False:
+        g9.append(g[9])
+print(g0)
+print(np.mean(g0))
+print(g1)
+print(np.mean(g1))
+print(g2)
+print(np.mean(g2))
+print(g3)
+print(np.mean(g3))
+print(g4)
+print(np.mean(g4))
+print(g5)
+print(np.mean(g5))
+print(g6)
+print(np.mean(g6))
+print(g7)
+print(np.mean(g7))
+print(g8)
+print(np.mean(g8))
+print(g9)
+print(np.mean(g9))
+
+mean.append(np.mean(g0))
+mean.append(np.mean(g1))
+mean.append(np.mean(g2))
+mean.append(np.mean(g3))
+mean.append(np.mean(g4))
+mean.append(np.mean(g5))
+mean.append(np.mean(g6))
+mean.append(np.mean(g7))
+mean.append(np.mean(g8))
+mean.append(np.mean(g9))
+
+print('maximum of mean:', np.max(mean))
+
+ax.plot(custom_amf, mean, 'o')
 
 # uniformed: 2018lepto1, 2018lepto5,  2018lepto76,  2018lepto98, 2020lepto06, 2019lepto24, 2020lepto06
-
-tau_uniform = []
-f_c_uniform = []
-fit_uniform = []
-fit_amf_uniform = []
-for ID in identifier_uniform:
-    print(ID)
-    amf = np.load('amf_%s.npy' %ID)
-    gain = np.load('gain_%s.npy' %ID)
-
-    sinv, sinc = curve_fit(gain_curve_fit, amf, gain)
-    print('tau:', sinv[0])
-    tau_uniform.append(sinv[0])
-    f_cutoff = abs(1 / (2*np.pi*sinv[0]))
-    print('f_cutoff:', f_cutoff)
-    print('alpha:', sinv[1])
-    f_c_uniform.append(f_cutoff)
-    fit_uniform.append(gain_curve_fit(amf, *sinv))
-    fit_amf_uniform.append(amf)
-
-
-colors = plt.cm.flag(np.linspace(0,1,len(fit_uniform)))
-
-#for ff ,f in enumerate(fit_uniform):
-#    ax.plot(fit_amf_uniform[ff], fit_uniform[ff],color =  colorss[ff])
-#    ax.axvline(x=f_c_uniform[ff], ymin=0, ymax=5, ls = '-', alpha = 0.5, color= colorss[ff])#colors_uniform[ff])
-
 ax.set_xscale('log')
 ax.set_yscale('log')
 ax.set_ylabel('gain [Hz/(mV/cm)]')
 ax.set_xlabel('envelope frequency [Hz]')
 ax.set_xlim(0.0007, 1.5)
 ax.set_ylim(0.001, 10)
-ax.plot(custom_amf, predict)
-#ax.plot(f_c, np.full(len(identifier), 0.0015), 'o', color = 'C0', alpha = 0.5, label = 'normal cutoff frequencies')
-#ax.plot(f_c_uniform, np.full(len(identifier_uniform), 0.002), 'o', alpha = 0.5, c = 'C0', label = 'uniform cutoff frequencies')
-#ax.legend(loc = 'center left')
-
-
 plt.show()
 embed()
diff --git a/apteronotus_code/sin_all_model.py b/apteronotus_code/sin_all_model.py
index e86bbe7..ed3885d 100644
--- a/apteronotus_code/sin_all_model.py
+++ b/apteronotus_code/sin_all_model.py
@@ -5,28 +5,24 @@ from IPython import embed
 from scipy.optimize import curve_fit
 from jar_functions import gain_curve_fit
 from jar_functions import avgNestedLists
+from matplotlib import gridspec
 
-
-identifier = ['2018lepto1',
-              '2018lepto4',
-              '2018lepto5',
-              '2018lepto76',
-              '2018lepto98',
-              '2019lepto03',
-              '2019lepto24',
-              '2019lepto27',
-              '2019lepto30',
-              '2020lepto04',
-              '2020lepto06',
-              '2020lepto16',
-              '2020lepto19',
-              '2020lepto20'
-              ]
+#plt.rcParams.update({'font.size': 16})
 
 amf = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1]
 
-custom_f = np.logspace(-2, -1, 10)
-custom_alpha = np.logspace(1.5, 1, 10)
+low_lim = 0.005
+high_lim = 6
+
+# subplot 1
+fig = plt.figure(figsize=(8.27,11.69))
+gs = gridspec.GridSpec(2, 2)
+ax0 = fig.add_subplot(gs[0,:])
+fig.text(0.06, 0.5, 'gain [Hz/(mV/cm)]', ha='center', va='center', rotation='vertical')
+fig.text(0.5, 0.04, 'envelope frequency [Hz]', ha='center', va='center')
+
+custom_f = np.logspace(-2, -1, 4)
+custom_alpha = np.logspace(0.6, 0.1, 4)
 c_gain = []
 custom_tau = abs(1 / (2 * np.pi * custom_f))
 for t, a in zip(custom_tau, custom_alpha):
@@ -35,70 +31,68 @@ for t, a in zip(custom_tau, custom_alpha):
         custom_g = gain_curve_fit(am, t, a)
         custom_gain.append(custom_g)
     c_gain.append(custom_gain)
-fig = plt.figure()
-ax = fig.add_subplot(111)
-ax.set_xscale('log')
-ax.set_yscale('log')
-for cc, c in enumerate(c_gain): 
-    ax.plot(amf, c)
-    ax.axvline(x=custom_f[cc], ymin=0, ymax=5, alpha=0.8)  # colors_uniform[ff])
+col = ['blue', 'orange', 'green', 'purple']
+for cc, c in enumerate(c_gain):
+    ax0.plot(amf, c, c = col[cc])
+    ax0.axvline(x=custom_f[cc], c = col[cc], ymin=0, ymax=5, alpha=0.5)
 
-plt.show()
+mean = avgNestedLists(c_gain)
+
+ax0.set_xscale('log')
+ax0.set_yscale('log')
+ax0.set_ylim(low_lim, high_lim)
+ax0.plot(amf, mean, lw = 3, c = 'r')
+
+# subplot 2
+ax1 = fig.add_subplot(gs[1,0])
+
+custom_f = np.logspace(-2, -1, 10)
+custom_alpha = np.logspace(0.6, 0.1, 10)
+c_gain = []
+custom_tau = abs(1 / (2 * np.pi * custom_f))
+for t, a in zip(custom_tau, custom_alpha):
+    custom_gain = []
+    for am in amf:
+        custom_g = gain_curve_fit(am, t, a)
+        custom_gain.append(custom_g)
+    c_gain.append(custom_gain)
+col = ['blue', 'orange', 'green']
+for cc, c in enumerate(c_gain):
+    ax1.plot(amf, c, c = 'C0')
+    ax1.axvline(x=custom_f[cc], c = 'C0', ymin=0, ymax=5, alpha=0.5)  # colors_uniform[ff])
 
 mean = avgNestedLists(c_gain)
-fig = plt.figure()
-ax = fig.add_subplot(111)
-ax.set_xscale('log')
-ax.set_yscale('log')
-ax.plot(amf, mean)
-plt.show()
 
-all = []
+ax1.set_xscale('log')
+ax1.set_yscale('log')
+ax1.set_ylim(low_lim, high_lim)
+ax1.plot(amf, mean, lw = 3, c = 'r')
 
-for ident in identifier:
-    data = np.load('gain_%s.npy' %ident)
-    all.append(data)
+# subplot 3
+ax2 = fig.add_subplot(gs[1,1])
 
-av = avgNestedLists(all)
-embed()
+custom_f = np.logspace(-2.75, -0.25, 10)
+custom_alpha = np.logspace(0.6, 0.1, 10)
+c_gain = []
+custom_tau = abs(1 / (2 * np.pi * custom_f))
+for t, a in zip(custom_tau, custom_alpha):
+    custom_gain = []
+    for am in amf:
+        custom_g = gain_curve_fit(am, t, a)
+        custom_gain.append(custom_g)
+    c_gain.append(custom_gain)
+col = ['blue', 'orange', 'green']
+for cc, c in enumerate(c_gain):
+    ax2.plot(amf, c, c = 'C0')
+    ax2.axvline(x=custom_f[cc], c = 'C0', ymin=0, ymax=5, alpha=0.5)  # colors_uniform[ff])
 
-fig = plt.figure(figsize=(8.27,11.69/2))
-ax = fig.add_subplot(111)
-ax.plot(amf, av, 'o', c = 'C0', label = 'gain')
-
-#plt.show()
-
-tau = []
-f_c = []
-fit = []
-fit_amf = []
-for ID in identifier:
-    print(ID)
-    amf = np.load('amf_%s.npy' %ID)
-    gain = np.load('gain_%s.npy' %ID)
-
-    sinv, sinc = curve_fit(gain_curve_fit, amf, gain, [2, 3])
-    #print('tau:', sinv[0])
-    tau.append(sinv[0])
-    f_cutoff = abs(1 / (2*np.pi*sinv[0]))
-    print('f_cutoff:', f_cutoff)
-    f_c.append(f_cutoff)
-    fit.append(gain_curve_fit(amf, *sinv))
-    fit_amf.append(amf)
-col = plt.cm.magma(np.linspace(0,0.8,len(fit)))
-for ff ,f in enumerate(fit):
-    ax.plot(fit_amf[ff], fit[ff], c = col[ff])
-    ax.axvline(x=f_c[ff], ymin=0, ymax=5, alpha=0.8, c = col[ff])  # colors_uniform[ff])
-
-ax.set_xscale('log')
-ax.set_yscale('log')
-ax.set_title('gain average all fish')
-ax.set_ylabel('gain [Hz/(mV/cm)]')
-ax.set_xlabel('envelope_frequency [Hz]')
-ax.set_ylim(0.0008, )
-#ax.plot(f_c, np.full(len(identifier), 0.0015), 'o', alpha = 0.5, c = 'darkorange', label = 'cutoff frequencies')
-ax.legend(loc = 'center left')
+mean = avgNestedLists(c_gain)
 
+ax2.set_xscale('log')
+ax2.set_yscale('log')
+ax2.set_ylim(low_lim, high_lim)
+ax2.set_yticklabels([])
+ax2.plot(amf, mean, lw = 3, c = 'r')
 plt.show()
 
 embed()
diff --git a/eigenmannia_code/eigenmannia_jar_stacked.py b/eigenmannia_code/eigenmannia_jar_stacked.py
index e281e80..f0bf0bf 100644
--- a/eigenmannia_code/eigenmannia_jar_stacked.py
+++ b/eigenmannia_code/eigenmannia_jar_stacked.py
@@ -12,12 +12,12 @@ from jar_functions import get_time_zeros
 from jar_functions import import_data_eigen
 from scipy.signal import savgol_filter
 
-plt.rcParams.update({'font.size': 18})
+plt.rcParams.update({'font.size': 12})
 
 base_path = 'D:\\jar_project\\JAR\\eigenmannia\\deltaf'
 
 #2015eigen8 no nix files
-identifier = [#'2013eigen13',
+identifier = ['2013eigen13',
               '2015eigen16','2015eigen17', '2015eigen19', '2020eigen22','2020eigen32']
 
 
@@ -32,8 +32,8 @@ for ID in identifier:
         delta_f, duration = parse_stimuli_dat(stimuli_dat)
         dur = int(duration[0][0:2])
         print(delta_f)
-        if delta_f != [4.0]:
-            continue
+        #if delta_f != [-4.0]:
+        #   continue
         data, pre_data, dt = import_data_eigen(datapath)
 
         #hstack concatenate: 'glue' pre_data and data
@@ -74,24 +74,24 @@ for ID in identifier:
             if i > 15 and i < 55:
                 j.append(jar4[idx])
 
-        r = np.median(j) - np.median(b)
+        r = (np.median(j) - np.median(b)) / 4       # divided by 4 cause of data at 4th harmonic, therefore response 4 times higher
         print('response:', r)
         deltaf.append(delta_f[0])
         response.append(r)
 
-        plt.figure(figsize = (14,8))
+        plt.figure(figsize = (8.27,11.69/2))
         plt.imshow(spec4, cmap='jet', origin='lower', extent=(times[0], times[-1], lim0, lim1), aspect='auto', vmin=-80, vmax=-10)
-        plt.plot(cut_time_jar, jar4, 'k', label = 'peak detection trace', lw = 2)
-        plt.hlines(y=lim0 + 5, xmin=10, xmax=70, lw=4, color='yellow', label='stimulus duration')
+        plt.plot(cut_time_jar, jar4, color = 'k', label = 'peak detection trace', lw = 2)
         plt.hlines(y=lim0 + 5, xmin=0, xmax=10, lw=4, color='red', label='pause')
+        plt.hlines(y=lim0 + 5, xmin=10, xmax=70, lw=4, color='gold', label='stimulus duration')
         plt.title('spectogram %s, deltaf: %sHz' %tuple(ID_delta_f))
         plt.xlim(times[0],times[-1])
         #embed()
-        #plt.xticks((times[0], 10, 20, 30, 40, 50, 60, times[-1]), [0, 10, 20, 30 ,40, 50, 60, 70])
+        plt.xticks((times[0], 10, 20, 30, 40, 50, 60, times[-1]), [0, 10, 20, 30 ,40, 50, 60, 70])
         plt.xlabel('time [s]')
         plt.ylabel('frequency [Hz]')
         plt.legend(loc = 'best')
-        plt.show()
+        #plt.show()
         delta_f_ID = [str(delta_f[0]).split('.')[0], ID]
 
         plt.close()
@@ -99,4 +99,4 @@ for ID in identifier:
 
     res_df = sorted(zip(deltaf,response))
 
-    #np.save('res_df_%s_new' %ID, res_df)
+    np.save('res_df_%s_new' %ID, res_df)
diff --git a/eigenmannia_code/eigenmannia_jar_subplot.py b/eigenmannia_code/eigenmannia_jar_subplot.py
index 4742f0e..6ee11ea 100644
--- a/eigenmannia_code/eigenmannia_jar_subplot.py
+++ b/eigenmannia_code/eigenmannia_jar_subplot.py
@@ -12,7 +12,7 @@ from jar_functions import get_time_zeros
 from jar_functions import import_data_eigen
 from scipy.signal import savgol_filter
 
-#plt.rcParams.update({'font.size': 18})
+plt.rcParams.update({'font.size': 10})
 
 base_path = 'D:\\jar_project\\JAR\\eigenmannia\\deltaf'
 
@@ -107,6 +107,7 @@ ax0.set_xlim(times[0],times[-1])
 ax0.set_ylabel('frequency [Hz]')
 ax0.axes.xaxis.set_ticklabels([])
 ax0.set_title('∆F -2 Hz')
+plt.xticks((1.7, 10, 20, 30, 40, 50, 60, times[-1]))
 
 ax1 = fig.add_subplot(222)
 ax1.imshow(specs[2], cmap='jet', origin='lower', extent=(times[0], times[-1], sub_lim0[2], sub_lim1[2]), aspect='auto', vmin=-80, vmax=-10)
@@ -116,6 +117,7 @@ ax1.axes.xaxis.set_ticklabels([])
 #ax1.axes.yaxis.set_ticklabels([])
 ax1.set_title('∆F 2 Hz')
 ax1.get_shared_y_axes().join(ax0, ax1)
+plt.xticks((1.7, 10, 20, 30, 40, 50, 60, times[-1]))
 
 ax2 = fig.add_subplot(223)
 ax2.imshow(specs[1], cmap='jet', origin='lower', extent=(times[0], times[-1], sub_lim0[1], sub_lim1[1]), aspect='auto', vmin=-80, vmax=-10)
@@ -124,6 +126,7 @@ ax2.set_xlim(times[0],times[-1])
 ax2.set_ylabel('frequency [Hz]')
 ax2.set_xlabel('time [s]')
 ax2.set_title('∆F -10 Hz')
+plt.xticks((1.7, 10, 20, 30, 40, 50, 60, times[-1]), [0, 10, 20, 30 ,40, 50, 60, 70])
 
 ax3 = fig.add_subplot(224)
 ax3.imshow(specs[3], cmap='jet', origin='lower', extent=(times[0], times[-1], sub_lim0[3], sub_lim1[3]), aspect='auto', vmin=-80, vmax=-10)
@@ -132,7 +135,7 @@ ax3.set_xlim(times[0],times[-1])
 ax3.set_xlabel('time [s]')
 #ax3.axes.yaxis.set_ticklabels([])
 ax3.set_title('∆F 10 Hz')
-
+plt.xticks((1.7, 10, 20, 30, 40, 50, 60, times[-1]), [0, 10, 20, 30 ,40, 50, 60, 70])
 plt.subplots(sharex = True, sharey = True)
 plt.show()
 
diff --git a/eigenmannia_code/figure_eigen_gain_plot.py b/eigenmannia_code/figure_eigen_gain_plot.py
new file mode 100644
index 0000000..cc93b3b
--- /dev/null
+++ b/eigenmannia_code/figure_eigen_gain_plot.py
@@ -0,0 +1,128 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pylab
+from IPython import embed
+from scipy.optimize import curve_fit
+from matplotlib.mlab import specgram
+import os
+from jar_functions import gain_curve_fit
+
+plt.rcParams.update({'font.size': 10})
+
+identifier = ['2015eigen8',
+              '2015eigen15',
+              '2015eigen16',
+              '2015eigen17',
+              '2015eigen19'
+              ]
+
+amfs = []
+gains = []
+maxgains = []
+mingains = []
+taus = []
+f_cs = []
+predicts = []
+for ID in identifier:
+    predict = []
+
+    print(ID)
+    amf = np.load('eigen_amf_%s.npy' % ID)
+    amfs.append(amf)
+    gain = np.load('eigen_gain_%s.npy' % ID)
+    gains.append(gain)
+    print(np.max(gain))
+    sinv, sinc = curve_fit(gain_curve_fit, amf, gain, [2, 3])
+    # print('tau:', sinv[0])
+    taus.append(sinv[0])
+    f_cutoff = abs(1 / (2 * np.pi * sinv[0]))
+    print('f_cutoff:', f_cutoff)
+    f_cs.append(f_cutoff)
+
+    print('min gain:', np.min(gain))
+    print('max gain:', np.max(gain))
+    maxgains.append(np.max(gain))
+    mingains.append(np.min(gain))
+    # predict of gain
+    for f in amf:
+        G = np.max(gain) / np.sqrt(1 + (2 * ((np.pi * f * sinv[0]) ** 2)))
+        predict.append(G)
+    predicts.append(predict)
+
+print('max of absolute max gain:', np.max(maxgains))
+print('min of absolute max gain:', np.min(maxgains))
+print('max of absolute min gain:', np.max(mingains))
+print('min of absolute min gain:', np.min(mingains))
+print('absolute max f_c:', np.max(f_cs))
+print('absolute min f_c:', np.min(f_cs))
+
+sort = sorted(zip(f_cs, identifier))
+print(sort)
+# order of plotting: 2018lepto1, 2018lepto5, 2018lepto76, 2018lepto98, 2019lepto24, 2020lepto06
+# order of f_c: 2019lepto24, 2020lepto06, 2018lepto98, 2018lepto76, 2018lepto1, 2018lepto5
+
+fig = plt.figure(figsize=(8.27, 11.69))
+# ax0 = plt.subplot2grid(shape=(3,4), loc=(0,0), colspan = 2)
+# ax1 = plt.subplot2grid((3,4), (0,2), colspan = 2)
+# ax2 = plt.subplot2grid((3,4), (1,0), colspan = 2)
+# ax3 = plt.subplot2grid((3,4), (1,2), colspan = 2)
+# ax4 = plt.subplot2grid((3,4), (2,0), colspan = 2)
+
+ax0 = fig.add_subplot(321)
+fig.text(0.05, 0.5, 'gain [Hz/(mV/cm)]', ha='center', va='center', rotation='vertical')
+fig.text(0.5, 0.04, 'envelope frequency [Hz]', ha='center', va='center')
+
+ax0.set_xlim(0.0007, 1.5)
+ax0.set_ylim(0.001, 10)
+ax0.plot(amfs[1], gains[1], 'o', label='gain')
+ax0.plot(amfs[1], predicts[1], label='fit')
+ax0.axvline(x=f_cs[1], ymin=0, ymax=5, ls='-', alpha=0.5, label='cutoff frequency')
+ax0.set_xscale('log')
+ax0.set_yscale('log')
+ax0.axes.xaxis.set_ticklabels([])
+
+ax1 = fig.add_subplot(322)
+ax1.set_xlim(0.0007, 1.5)
+ax1.set_ylim(0.001, 10)
+ax1.plot(amfs[0], gains[0], 'o', label='gain')
+ax1.plot(amfs[0], predicts[0], label='fit')
+ax1.axvline(x=f_cs[0], ymin=0, ymax=5, ls='-', alpha=0.5, label='cutoff frequency')
+ax1.set_xscale('log')
+ax1.set_yscale('log')
+ax1.axes.yaxis.set_ticklabels([])
+ax1.axes.xaxis.set_ticklabels([])
+
+ax2 = fig.add_subplot(323)
+ax2.set_xlim(0.0007, 1.5)
+ax2.set_ylim(0.001, 10)
+ax2.plot(amfs[4], gains[4], 'o', label='gain')
+ax2.plot(amfs[4], predicts[4], label='fit')
+ax2.axvline(x=f_cs[4], ymin=0, ymax=5, ls='-', alpha=0.5, label='cutoff frequency')
+ax2.set_xscale('log')
+ax2.set_yscale('log')
+ax2.axes.xaxis.set_ticklabels([])
+
+ax3 = fig.add_subplot(324)
+ax3.set_xlim(0.0007, 1.5)
+ax3.set_ylim(0.001, 10)
+ax3.plot(amfs[2], gains[2], 'o', label='gain')
+ax3.plot(amfs[2], predicts[2], label='fit')
+ax3.axvline(x=f_cs[2], ymin=0, ymax=5, ls='-', alpha=0.5, label='cutoff frequency')
+ax3.set_xscale('log')
+ax3.set_yscale('log')
+ax3.axes.yaxis.set_ticklabels([])
+
+ax4 = fig.add_subplot(325)
+ax4.set_xlim(0.0007, 1.5)
+ax4.set_ylim(0.001, 10)
+ax4.plot(amfs[3], gains[3], 'o', label='gain')
+ax4.plot(amfs[3], predicts[3], label='fit')
+ax4.axvline(x=f_cs[3], ymin=0, ymax=5, ls='-', alpha=0.5, label='cutoff frequency')
+ax4.set_xscale('log')
+ax4.set_yscale('log')
+
+# plt.legend(loc = 'lower left')
+plt.show()
+
+# np.save('f_c', f_c)
+# np.save('tau', tau)
diff --git a/eigenmannia_code/figure_eigen_jar_plot.py b/eigenmannia_code/figure_eigen_jar_plot.py
new file mode 100644
index 0000000..68b72e3
--- /dev/null
+++ b/eigenmannia_code/figure_eigen_jar_plot.py
@@ -0,0 +1,168 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pylab
+from IPython import embed
+from scipy.optimize import curve_fit
+from scipy.optimize import curve_fit
+from matplotlib.mlab import specgram
+import os
+
+from jar_functions import import_data
+from jar_functions import import_amfreq
+from jar_functions import sin_response
+from jar_functions import mean_noise_cut
+from jar_functions import gain_curve_fit
+
+#plt.rcParams.update({'font.size': 10})
+
+def take_second(elem):      # function for taking the names out of files
+    return elem[1]
+
+identifier = ['2015eigen8',
+              '2015eigen15',
+              '2015eigen16',
+              '2015eigen17',
+              '2015eigen19'
+              ]
+for ident in identifier:
+
+    times = []
+    jars = []
+    jms = []
+    amfreq = []
+
+    times1 = []
+    jars1 = []
+    jms1 = []
+    amfreq1 = []
+
+    amf = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1]
+
+    data = sorted(np.load('eigen_%s files.npy' %ident), key = take_second)      # list with filenames in it
+
+    for i, d in enumerate(data):
+        dd = list(d)
+        if dd[1] == '1' or dd[1] == '0.2' or dd[1] == '0.05' or dd[1] == '0.01' or dd[1] == '0.005' or dd[1] == '0.001':
+            jar = np.load('eigen_%s.npy' %dd)     # load data for every file name
+            jm = jar - np.mean(jar)         # low-pass filtering by subtracting mean
+
+            time = np.load('eigen_%s time.npy' %dd)       # time file
+            dt = time[1] - time[0]
+
+            n = int(1/float(d[1])/dt)
+            cutf = mean_noise_cut(jm, n = n)
+            cutt = time
+            if dd[1] == '0.001':
+                amfreq1.append(dd[1])
+                jars1.append(jm - cutf)
+                jms1.append(jm)
+                times1.append(time)
+            if dd[1] not in amfreq:
+                print(dd)
+                amfreq.append(dd[1])
+                jars.append(jm - cutf)
+                jms.append(jm)
+                times.append(time)
+            else:
+                #print('1:', dd)
+                amfreq1.append(dd[1])
+                jars1.append(jm - cutf)
+                jms1.append(jm)
+                times1.append(time)
+    if len(jars) != 6:
+        continue
+
+    ssample = 100000
+
+    fig = plt.figure(figsize=(8.27, 11.69))
+    fig.suptitle('%s' % ident)
+    fig.text(0.06, 0.5, 'fish frequency [Hz]', ha='center', va='center', rotation='vertical', color='C0')
+    fig.text(0.97, 0.5, 'stimulus amplitude [mV/cm]', ha='center', va='center', rotation='vertical', color='red')
+    fig.text(0.5, 0.04, 'time [s]', ha='center', va='center')
+
+    ax0 = fig.add_subplot(611)
+    print('absolute frequency shift 0.001Hz:', np.max(jars[0]) - np.min(jars[0]))
+    ax0.plot(times[0], jars[0], zorder=20)
+    # ax0.set_zorder(1)
+    ax0.set_ylim(-12, 12)
+
+    lower0 = 0
+    upper0 = 2000
+    x0 = np.linspace(lower0, upper0, sample)
+    y0 = (sin_response(np.linspace(lower0, upper0, sample), 0.001, np.pi / 2, .35) + 0.5)
+    ax0_0 = ax0.twinx()
+    ax0_0.set_ylim(-0.2, 1.2)
+    ax0_0.plot(x0, y0, color='red', zorder=1, alpha=0.5)
+    # ax0_0.set_zorder(2)
+
+    ax1 = fig.add_subplot(612)
+    print('absolute frequency shift 0.005 Hz:', np.max(jars[1]) - np.min(jars[1]))
+    ax1.plot(times[1], jars[1])
+    ax1.set_ylim(-12, 12)
+
+    lower1 = 0
+    upper1 = 400
+    x1 = np.linspace(lower1, upper1, sample)
+    y1 = (sin_response(np.linspace(lower1, upper1, sample), 0.005, np.pi / 2, .35) + 0.5)
+    ax1_0 = ax1.twinx()
+    ax1_0.set_ylim(-0.2, 1.2)
+    ax1_0.plot(x1, y1, color='red', alpha=0.5)
+
+    ax2 = fig.add_subplot(613)
+    print('absolute frequency shift 0.01 Hz:', np.max(jars[2]) - np.min(jars[2]))
+    ax2.plot(times[2], jars[2])
+    ax2.set_ylim(-12, 12)
+
+    lower2 = 0
+    upper2 = 400
+    x2 = np.linspace(lower2, upper2, sample)
+    y2 = (sin_response(np.linspace(lower2, upper2, sample), 0.01, np.pi / 2, 0.35) + 0.5)
+    ax2_0 = ax2.twinx()
+    ax2_0.set_ylim(-0.2, 1.2)
+    ax2_0.plot(x2, y2, color='red', alpha=0.5)
+
+    ax3 = fig.add_subplot(614)
+    print('absolute frequency shift 0.02 Hz:', np.max(jars[3]) - np.min(jars[3]))
+    ax3.plot(times[3], jars[3])
+    ax3.set_ylim(-12, 12)
+
+    lower3 = 0
+    upper3 = 200
+    x3 = np.linspace(lower3, upper3, sample)
+    y3 = (sin_response(np.linspace(lower3, upper3, sample), 0.05, np.pi / 2, 0.35) + 0.5)
+    ax3_0 = ax3.twinx()
+    ax3_0.set_ylim(-0.2, 1.2)
+    ax3_0.plot(x3, y3, color='red', alpha=0.5)
+
+    ax4 = fig.add_subplot(615)
+    print('absolute frequency shift 0.5 Hz:', np.max(jars[4]) - np.min(jars[4]))
+    ax4.plot(times[4], jars[4])
+    ax4.set_ylim(-12, 12)
+
+    lower4 = 0
+    upper4 = 200
+    x4 = np.linspace(lower4, upper4, sample)
+    y4 = (sin_response(np.linspace(lower4, upper4, sample), 0.2, np.pi / 2, 0.35) + 0.5)
+    ax4_0 = ax4.twinx()
+    ax4_0.set_ylim(-0.2, 1.2)
+    ax4_0.plot(x4, y4, color='red', alpha=0.5)
+
+    ax5 = fig.add_subplot(616)
+    print('absolute frequency shift 1 Hz:', np.max(jars[5]) - np.min(jars[5]))
+    ax5.plot(times[5], jars[5])
+    ax5.set_ylim(-12, 12)
+
+    lower5 = 0
+    upper5 = 200
+    x5 = np.linspace(lower5, upper5, sample)
+    y5 = (sin_response(np.linspace(lower5, upper5, sample), 1, np.pi / 2, 0.35) + 0.5)
+    ax5_0 = ax5.twinx()
+    ax5_0.plot(x5, y5, color='red', lw=0.5, alpha=0.5)
+    ax5_0.set_ylim(-0.2, 1.2)
+    plt.subplots_adjust(left=0.125,
+                        bottom=0.1,
+                        right=0.9,
+                        top=0.9,
+                        wspace=0.1,
+                        hspace=0.35)
+    plt.show()
\ No newline at end of file