diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..04dc00b8d99217012156bafe86bf36669404d8f0
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,5 @@
+__pycache__/
+.ipynb_checkpoints/
+data/
+~tmp/
+*.pkl
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000000000000000000000000000000000000..014adadf8c150f53455d9e695fa285e11c6c7b9a
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,20 @@
+# Required python packages for the simulation
+
+# Install all packages with `pip install -r requirements.txt`.
+# To get a list of installed packages use `pip freeze`.
+
+
+# General tools
+matplotlib
+numpy
+pandas
+pint
+scipy
+
+
+# For SNLO helper
+pyautogui
+pyperclip
+
+
+# requires analysis.modules_JFK from the eg-labor repository.
diff --git a/snlo-helper/OPA-simulation.py b/snlo-helper/OPA-simulation.py
new file mode 100644
index 0000000000000000000000000000000000000000..4b1c315ec3311311c9a67513f6d46bb0b490d767
--- /dev/null
+++ b/snlo-helper/OPA-simulation.py
@@ -0,0 +1,797 @@
+# -*- coding: utf-8 -*-
+"""
+
+created on 25.05.2022 by Jan Frederic Kinder
+"""
+import time
+import pickle
+import pyautogui as gui
+import matplotlib.pyplot as plt
+import modules_JFK as m  # module from NLOQO/EG-Labor/DataAnalysis, check PYTHONPATH in Spyder
+import numpy as np
+from pyperclip import paste
+from socket import gethostname
+from collections import defaultdict
+from datetime import datetime  # get current time
+
+import snlohelper as snlo
+# imports all of snlohelper, risk of namespace-confusion...
+from snlohelper import *
+
+"""define some global filepath variables based upon the computer name:"""
+def get_results_path():
+    pc_name = gethostname()
+    if pc_name == 'WhiteShard':
+        return 'M:/'
+    elif pc_name == 'POPS':
+        return 'C:/Users/kinder.NLOQO/HESSENBOX-DA/OPA-Paper/analysis (python)/simulation/'
+    elif pc_name == 'Myres':
+        return 'D:/'
+    else:
+        print('data location not specified, use generic path: M:/')
+        return 'M:/'
+    
+def get_plots_path():
+    pc_name = gethostname()
+    if pc_name == 'WhiteShard':
+        return 'M:/'
+    elif pc_name == 'POPS':
+        return f'C:/Users/kinder.NLOQO/HESSENBOX-DA/OPA-Paper/analysis (python)/Ausgabe/'
+    elif pc_name == 'Myres':
+        return 'D:/'
+    else:
+        print('data location not specified, use generic path: M:/')
+        return 'M:/'
+results_path = get_results_path()#'C:/Users/kinder.NLOQO/HESSENBOX-DA/OPA-Paper/analysis (python)/simulation/'
+plots_path = get_plots_path()#f'C:/Users/kinder.NLOQO/HESSENBOX-DA/OPA-Paper/analysis (python)/Ausgabe/'
+
+
+#%% tests
+
+#%%% test for display of pulse durations:
+id_beam = snlo.import_snlo_file('ID_BEAM.dat')
+sig_beam = snlo.import_snlo_file('SIG_BEAM.dat')
+pmp_beam = snlo.import_snlo_file('PMP_BEAM.dat')
+
+# fit idler duration
+x = 1e9*id_beam.T[0]  # ns
+y = id_beam.T[1]/max(id_beam.T[1])
+idl = m.gaussian_fit(np.array([x, y]).T, False)[0]
+fwhm_ns_idl = abs(2 * np.sqrt(2 * np.log(2)) * idl[2])
+x_lists = [x.copy()]
+y_lists = [y.copy()]
+
+# fit signal duration
+x = 1e9*sig_beam.T[0]  # ns
+y = sig_beam.T[1]/max(sig_beam.T[1])
+sig = m.gaussian_fit(np.array([x, y]).T, False)[0]
+fwhm_ns_sig = abs(2 * np.sqrt(2 * np.log(2)) * sig[2])
+x_lists.append(x)
+y_lists.append(y)
+
+# fit pump duration
+x = 1e9*pmp_beam.T[0]  # ns
+y = pmp_beam.T[1]/max(pmp_beam.T[1])
+pmp = m.gaussian_fit(np.array([x, y]).T, False)[0]
+fwhm_ns_pmp = abs(2 * np.sqrt(2 * np.log(2)) * pmp[2])  
+x_lists.append(x)
+y_lists.append(y)
+
+fig, ax = plt.subplots(1, 3)
+fits = idl, sig, pmp
+fwhms = fwhm_ns_idl, fwhm_ns_sig, fwhm_ns_pmp
+titles = 'idler (Red1)', 'signal (Red2)', 'pump (Blue)'
+colors = 'red', 'green', 'blue'
+x_fit = np.arange(min(x), max(x), 0.1)
+
+for i in range(len(ax)):
+    m.plot_options(fig, ax[i], xlabel='time (ns)',  image_width=15, aspect_ratio=3)
+    ax[i].plot(x_lists[i], y_lists[i], '.', color=colors[i])
+    ax[i].set_title(titles[i], color=colors[i])
+    ax[i].plot(x_fit, m.gauss(x_fit, *fits[i]), '-', color=colors[i])
+    ax[i].text(0, 0, f'FWHM\n{fwhms[i]:.1f}ns', ha='center', color=colors[i])
+m.plot_options(fig, ax[0], image_width=15, aspect_ratio=3, xlabel='time (ns)', ylabel='normalized power')
+plt.tight_layout()
+plt.show()
+
+#%%% screenshot with offset (top left corner)
+off = 300, 300
+image = gui.screenshot(region=(off[0],off[1],500,550))
+plt.figure()
+plt.imshow(image)
+plt.show()
+
+#%%% import SNLO files
+# time, power, phase, Mx^2, My^2, Rad Curv x, Rad Curv y, X-tilt, w_x^2, w_y^2
+with open('C:/SNLO/ID_BEAM.dat', 'r') as f:
+    idler = f.read()
+with open('C:/SNLO/SIG_BEAM.dat', 'r') as f:
+    signal = f.read()
+with open('C:/SNLO/PMP_BEAM.dat', 'r') as f:
+    pump = f.read()
+print(idler)
+
+#%%% plot spectra
+with open('C:/SNLO/OPA2D_SP.dat', 'r') as f:
+    input = f.readlines()
+spectr = []
+for i in input:
+    element = i.split()
+    #pectr.append(float(j) for j in element)
+    spectr.append([float(element[0]), float(element[1]), float(element[2]), float(element[3])])
+spectr = np.array(spectr.copy())
+
+x = spectr.T[0]  # in MHz
+y_idl = spectr.T[1]
+y_sig = spectr.T[2]
+y_pmp = spectr.T[3]
+
+idl = m.gaussian_fit(np.array([x, y_idl]).T, False)[0]
+sig = m.gaussian_fit(np.array([x, y_sig]).T, False)[0]
+pmp = m.gaussian_fit(np.array([x, y_pmp]).T, False)[0]
+
+fwhm_idl = 2 * np.sqrt(2 * np.log(2)) * idl[2]
+fwhm_sig = 2 * np.sqrt(2 * np.log(2)) * sig[2]
+fwhm_pmp = 2 * np.sqrt(2 * np.log(2)) * pmp[2]
+print(f'idl: {fwhm_idl:.3f}MHz => {2*np.log(2)/np.pi/fwhm_idl*1e3:.3f}ns')
+print(f'sig: {fwhm_sig:.3f}MHz => {2*np.log(2)/np.pi/fwhm_sig*1e3:.3f}ns')
+print(f'pmp: {fwhm_pmp:.3f}MHz => {2*np.log(2)/np.pi/fwhm_pmp*1e3:.3f}ns')
+
+fig, ax = plt.subplots()
+m.plot_options(fig, ax, xlabel='detuning (MHz)', ylabel='intensity (arb. u.)')
+ax.plot(x,y_idl, 'r.', label='idler')
+ax.plot(x, y_sig, 'g.', label='signal')
+ax.plot(x, y_pmp, 'b.', label='pump')
+
+x_fit = np.arange(min(x), max(x), 1)
+ax.plot(x_fit, m.gauss(x_fit, *idl), 'r')
+ax.plot(x_fit, m.gauss(x_fit, *sig), 'g')
+ax.plot(x_fit, m.gauss(x_fit, *pmp), 'b')
+
+ax.legend(loc='upper left')
+plt.tight_layout()
+plt.show()
+
+
+#%%% import simulation files
+now = datetime.now()
+file_name = now.strftime("%Y-%m-%d") + '_' + 'OPA1_pump-dependence.pkl'
+result = pickle.load(open('simulation/' + file_name, 'rb'))
+
+
+#%% parameters for OPA1
+#%%% calc radius of curvature OPA1
+wavelength_nm = 1064.16
+ref_index = snlo.n_lnb(wavelength_nm, axis='o')
+waist_size_mm = 0.4897
+focus_pos_mm = 25  # focus center of crystal
+result = snlo.focus(wavelength_nm, ref_index, waist_size_mm, focus_pos_mm)
+print(f'radius of curvature: {result[2]:.3f}mm')
+
+
+#%%% parameters 2D mix LP for OPA1
+idler_nm = 3545
+signal_nm = 1520.6
+pump_nm = 1064.16  # check!
+seedpower_w = 0.96
+pump_j = 0.000000499
+pumpduration_ns = 8
+fwhm_seed_mm = 0.65/2*1.18  # 1/e^2 radius to FWHM
+fwhm_pump_mm = 0.83/2*1.18  # 1/e^2 radius to FWHM
+
+values = [
+    [idler_nm, signal_nm, pump_nm],
+    [snlo.n_lnb(idler_nm, axis='o'), snlo.n_lnb(signal_nm, axis='o'), snlo.n_lnb(pump_nm, axis='o')],
+    [0, 0, 0],
+    [0.01, 0.02, 0.01],  # coating specs input
+    [0.01, 0.02, 0.01],  # coating specs output
+    [0, 0, 0],  # crystal loss 1/mm
+    [seedpower_w, 0.0000000000000001, pump_j],
+    [0, 0, pumpduration_ns],
+    [fwhm_seed_mm, fwhm_pump_mm, fwhm_pump_mm],
+    [1, 1, 1],  # supergaussian coefficient
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],  # walkoff
+    [0, 0, 0],  # offset
+    [766.97, 4258, 23270],  # radius of curvature
+    [50, 32, 32],  # integration points
+    [50, 0, 0], # crystal length, grid, 0: auto: 2.4485x1.9588mm
+    [14.77],  # via QPM, LNB_M
+    [0],  # delta k
+    [1000],  # distance to image
+    [25]  # time steps
+]
+# todo: crystal reflectivities (data sheet)
+
+#%%% init 2D-mix-LP parameters
+snlo.open_function('2D-mix-LP')
+for index in range(len(values)):
+    if len(values[index])==1:
+        pos = list(snlo.dict_2dmixlp.values())[index][0]
+        snlo.set_value(pos, values[index][0])
+    else:
+        for element in range(len(values[index])):
+            snlo.set_value(list(snlo.dict_2dmixlp.values())[index][element], values[index][element])
+
+
+#%% OPA1 pump energy dependence
+# start with 1e-5
+energies = np.arange(1e-6,0.7e-3,1e-5)
+print(energies)
+print(len(energies))
+
+
+#%%% run OPA1 pump energy dependence
+results = []
+for energy in energies:
+    set_value(snlo.dict_2dmixlp['Energy/power (J or W)'][2], energy)
+    snlo.moveto(*dict_2dmixlp['Accept'])
+    #snlo.moveto(780, 800)
+    gui.click()
+    dict = run_2dmixlp()
+    m.set_key(dict, 'pump pulse energy (J)', energy)
+    dict2 = get_spectr()
+    results.append(merge_dict([dict, dict2]))
+    snlo.moveto(*dict_2dmixlp['Change Inputs'])
+        #(600, 250)
+    gui.click()
+
+#%% save OPA1 pump dependence
+now = datetime.now()
+file_name = now.strftime("%Y-%m-%d") + '_' + 'OPA1_pump-dependence.pkl'
+pickle.dump(results, open(results_path + file_name, 'wb'))
+
+
+#%%% plot results
+file_name = '2022-06-09_OPA1_pump-dependence.pkl'
+#results = snlo.merge_dict(pickle.load(open('simulation/' + file_name, 'rb')))
+results = snlo.merge_dict(pickle.load(open(results_path + file_name, 'rb')))
+
+x = np.array(results['pump pulse energy (J)']).T[0]*1e6
+y = np.array(results['Output pulse energy (mJ)']).T[0][0]*1e3
+y2 = np.array(results['duration_idl_ns']).T[0]
+y3 = np.array(results['duration_sig_ns']).T[0]
+y4 = np.array(results['duration_pmp_ns']).T[0]
+
+fig, [ax1, ax2] = plt.subplots(1,2)
+m.plot_options(fig, ax1, image_width=15, aspect_ratio=2, xlabel='pump pulse energy ($\mathrm{\mu J}$)', ylabel='idler pulse energy ($\mathrm{\mu J}$)')
+ax1.plot(x, y, 'r-', label='idler pulse energy')
+ax1.legend(loc='upper left')
+ax12 = ax1.twinx()
+m.plot_options(fig, ax12, image_width=15, aspect_ratio=2)
+ax12.plot(x, 1e-6*y/(1e-9*y2), 'b--', label='single pass gain')  # divided by 1W...
+ax12.set_ylabel('single pass gain', color="blue")
+ax12.legend(loc='lower right')
+
+m.plot_options(fig, ax2, ticks=['auto', 2], image_width=15, aspect_ratio=2, xlabel='pump pulse energy ($\mathrm{\mu J}$)',
+               ylabel='idler pulse duration (ns)')
+ax2.plot(x, y2, 'r-', label='idler')
+ax2.plot(x, y3, 'g--', label='sig')
+ax2.plot(x[y4<100], y4[y4<100], 'b-', label='pump')
+ax2.set_ylim([0, max(y4[y4<100]+0.5)])
+ax2.legend(loc='lower center', ncol=3, frameon=True, handlelength=1.5, mode="expand")
+plt.tight_layout()
+#m.save_plot('OPA1_pump_sim')
+plt.show()
+
+#%% conversion efficiency OPA1
+x = np.array(results['pump pulse energy (J)']).T[0]*1e6
+y = np.array(results['Output pulse energy (mJ)']).T*1e3
+eff = ((y[0]+y[1])/x)[0]
+
+fig, ax1 = plt.subplots()
+m.plot_options(fig, ax1, image_width=15, xlabel='pump pulse energy ($\mathrm{\mu J}$)', ylabel='pulse energy ($\mathrm{\mu J}$)')
+ax1.plot(x, y[0][0], 'r-', label='idler pulse energy')
+ax1.plot(x, y[1][0], 'g-', label='signal pulse energy')
+ax1.plot(x, y[2][0], 'b-', label='residual pump pulse energy')
+ax1.legend(loc='upper left')
+
+ax12 = ax1.twinx()
+m.plot_options(fig, ax12, image_width=15, aspect_ratio=2)
+ax12.plot(x[eff<1], eff[eff<1], 'k--', label='conversion efficiency')
+ax12.set_ylabel('conversion efficiency')
+
+ax12.legend(loc='lower right')
+
+plt.tight_layout()
+#m.save_plot('OPA1_pump_sim')
+plt.show()
+
+#%% OPA1 seed power dependence
+powers = np.arange(1e-6, 1.5, 20e-3)
+
+#%%% run seed dependence OPA1
+set_value(dict_2dmixlp['Energy/power (J or W)'][2], 0.5e-3)
+results = []
+for power in powers:
+    set_value(dict_2dmixlp['Energy/power (J or W)'][0], power)
+    gui.moveTo(*coords(780, 800))
+    gui.click()
+    dict1 = run_2dmixlp()
+    m.set_key(dict1, 'seed power (W)', power)
+    dict2 = get_spectr()
+    results.append(merge_dict([dict1, dict2]))
+    gui.moveTo(*coords(600, 250))
+    gui.click()
+# save OPA1 seed dependence
+now = datetime.now()
+file_name = now.strftime("%Y-%m-%d") + '_' + 'OPA1_seed-dependence.pkl'
+pickle.dump(results, open('simulation/' + file_name, 'wb'))
+
+#%%% import and plot OPA1 pump dependence
+file_name = '2022-06-09_OPA1_seed-dependence.pkl'
+results = merge_dict(pickle.load(open(results_path + file_name, 'rb')))
+
+x = np.array(results['seed power (W)']).T[0]
+y = np.array(results['Output pulse energy (mJ)']).T[0][0] * 1e3
+y2 = np.array(results['duration_idl_ns']).T[0]
+
+fig, [ax1, ax2] = plt.subplots(1, 2)
+m.plot_options(fig, ax1, xlabel='(cw) idler power (W)', ylabel='idler pulse energy ($\mathrm{\mu J}$)')
+ax1.plot(x, y, '-')
+
+m.plot_options(fig, ax2, xlabel='(cw) idler power (W)',
+               ylabel='idler pulse duration (ns)')
+ax2.plot(x, y2, '-')
+plt.tight_layout()
+# m.save_plot('OPA1_seed_sim')
+plt.show()
+
+#%% OPA1: Paper values
+pump_energy = 0.5e-3
+seed_power = 1
+set_value(dict_2dmixlp['Energy/power (J or W)'][2], pump_energy)
+set_value(dict_2dmixlp['Energy/power (J or W)'][0], seed_power)
+moveto(dict_2dmixlp['Accept'])
+gui.click()
+dict1 = run_2dmixlp()
+m.set_key(dict1, 'seed power (W)', seed_power)
+m.set_key(dict1, 'pump pulse energy (J)', pump_energy)
+dict2 = get_spectr()
+results = merge_dict([dict1, dict2])
+#%% gather OPA1 values
+duration_idl_ns = results['duration_idl_ns'][0]
+fwhm_idl_MHz = results['fwhm_idl_MHz'][0]
+energy_idl_mJ = results['Output pulse energy (mJ)'][0][0]
+y = np.array(results['Output pulse energy (mJ)'])*1e3
+efficiency = ((y[0, 0]+y[0, 1])/y[0, 2])
+gain = (1e-3*energy_idl_mJ)/(1e-9 * duration_idl_ns)/results['seed power (W)'][0]
+
+print(f'idler pulse energy: {energy_idl_mJ * 1e3:.3f}uJ')
+print(f'conversion efficiency: {efficiency*1e2:.3f}%')
+print(f'single pass gain: {gain:.1f}')
+print(f'idler duration:{duration_idl_ns:.3f}ns')
+print(f'idler bandwidth: {fwhm_idl_MHz:.3f}MHz')
+
+#%% parameters for OPA2
+#%%% calc radius of curvature OPA2 - pump
+idler_nm = 3545
+wavelength_nm = 1064.16
+ref_index = n(pm_theta(idler_nm, wavelength_nm), wavelength_nm, 25)
+waist_size_mm = 1.6/2*1.18 # 1/e^2 radius to FWHM
+focus_pos_mm = 15  # focus center of crystal
+result = focus(wavelength_nm, ref_index, waist_size_mm, focus_pos_mm)
+print(f'radius of curvature: {result[2]:.3f}mm')
+
+#%%% calc radius of curvature OPA2 - idler
+wavelength_nm = 3545
+ref_index = n_lnb(wavelength_nm, axis='o')
+waist_size_mm = 3.9/2*1.18  # 1/e^2 radius to FWHM
+focus_pos_mm = 15  # focus center of crystal
+result = focus(wavelength_nm, ref_index, waist_size_mm, focus_pos_mm)
+print(f'radius of curvature: {result[2]:.3f}mm')
+
+#%%% calc radius of curvature OPA2 - signal
+wavelength_nm = 1520.6
+ref_index = n_lnb(wavelength_nm, axis='o')
+waist_size_mm = 3.9/2*1.18  # 1/e^2 radius to FWHM
+focus_pos_mm = 15  # focus center of crystal
+result = focus(wavelength_nm, ref_index, waist_size_mm, focus_pos_mm)
+print(f'radius of curvature: {result[2]:.3f}mm')
+
+#%%% calc pulse duration of the idler seed
+file_name = '2022-05-23_OPA1_seed-dependence.pkl'
+results = merge_dict(pickle.load(open(results_path + file_name, 'rb')))
+i = 50
+print('pump pulse energy: 0.5mJ')
+power = results['seed power (W)'][i][0]
+print(f'idler seed power: {power}W')
+duration = results['duration_idl_ns'][i][0]
+print(f'idler duration: {duration:.3f}ns')
+energy = 1e3*results['Output pulse energy (mJ)'][i][0][0]
+print(f'idler pulse energy: {energy:.3f}uJ')
+
+#%%% parameters 2D-mix-LP for OPA2
+idler_nm = 3545
+signal_nm = 1520.6
+pump_nm = 1064.16
+seedenergy_j = 44.5e-6  # was 46.9uJ (seed: 1W, pump: 0.5mJ)
+pump_j = 23e-3
+pulseduration_idler_ns = 6.217  # abgelesen aus seed dep
+pumpduration_ns = 8
+fwhm_seed_mm = 3.9/2*1.18  # 1/e^2 radius to FWHM
+fwhm_pump_mm = 1.6/2*1.18  # 1/e^2 radius to FWHM
+
+values = [
+    [idler_nm, signal_nm, pump_nm],
+    [n_lnb(idler_nm, axis='o'), n_lnb(signal_nm, axis='o'), n(pm_theta(idler_nm, pump_nm), pump_nm, 25)],
+    [0, 0, 0],
+    [0.05, 0.005, 0.001],  # coating specs input, NC24, NC25
+    [0.05, 0.005, 0.001],  # coating specs output
+    [0, 0, 0],  # crystal loss 1/mm
+    [seedenergy_j, 1e-16, pump_j],
+    [pulseduration_idler_ns, 0, pumpduration_ns],
+    [fwhm_seed_mm, fwhm_seed_mm, fwhm_pump_mm],
+    [1, 1, 1],  # supergaussian coefficient
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 37.14],  # walkoff
+    [0, 0, 0],  # offset
+    [1631400, 9172700, 525250],  # todo: radius of curvature
+    [50, 32, 32],  # integration points
+    [30, 8, 8], # crystal length, grid in mm
+    [-4.04],  # via QPM, LNB_M
+    [0],  # delta k
+    [1000],  # distance to image
+    [25]  # time steps
+]
+# todo: import power time traces and calculate Gaussian FWHM for idler after OPA1 and OPA2-1
+
+#%%% init 2D-mix-LP parameters for OPA2
+open_function('2D-mix-LP')
+for index in range(len(values)):
+    if len(values[index])==1:
+        pos = list(dict_2dmixlp.values())[index][0]
+        set_value(pos, values[index][0])
+    else:
+        for element in range(len(values[index])):
+            set_value(list(dict_2dmixlp.values())[index][element], values[index][element])
+
+
+#%% OPA2 pump energy dependence
+# start with 1e-5
+energies = np.arange(1e-6, 25e-3, 250e-6)
+print(energies)
+print(len(energies))
+# set seed energy (was 20uJ)
+set_value(dict_2dmixlp['Energy/power (J or W)'][0], 4.45e-05)
+
+#%%% run OPA2 pump energy dependence
+results = []
+for energy in energies:
+    set_value(dict_2dmixlp['Energy/power (J or W)'][2], energy)
+    gui.moveTo(*coords(780, 800))
+    gui.click()
+    dict = run_2dmixlp()
+    m.set_key(dict, 'pump pulse energy (J)', energy)
+    dict2 = get_spectr()
+    results.append(merge_dict([dict, dict2]))
+    gui.moveTo(*coords(600, 250))
+    gui.click()
+
+#%% save OPA2 pump dependence
+now = datetime.now()
+file_name = now.strftime("%Y-%m-%d") + '_' + 'OPA2_pump-dependence.pkl'
+#pickle.dump(results, open('simulation/' + file_name, 'wb'))
+pickle.dump(results, open(results_path + file_name, 'wb'))
+
+
+
+#%%% plot results
+file_name = '2022-06-09_OPA2_pump-dependence.pkl'
+#results = merge_dict(pickle.load(open('simulation/' + file_name, 'rb')))
+results = merge_dict(pickle.load(open(results_path + file_name, 'rb')))
+
+
+x = np.array(results['pump pulse energy (J)']).T[0] * 1e3
+y = np.array(results['Output pulse energy (mJ)']).T[0][0]
+y2 = np.array(results['duration_idl_ns']).T[0]
+
+fig, [ax1, ax2] = plt.subplots(1, 2)
+m.plot_options(fig, ax1, xlabel='pump pulse energy (mJ)', ylabel='idler pulse energy (mJ)')
+ax1.plot(x,y, 'b-')
+
+m.plot_options(fig, ax2, xlabel='pump pulse energy (mJ)',
+               ylabel='idler pulse duration (ns)')
+ax2.plot(x[y2<100], y2[y2<100], '-')
+#ax2.set_ylim([0, 20])
+plt.tight_layout()
+#m.save_plot('OPA2_pump_sim')
+plt.show()
+
+#%% OPA2 seed energy dependence
+energies = np.arange(1e-7, 40e-6, 0.5e-6)
+print(energies)
+print(len(energies))
+
+#%%% run seed dependence OPA2
+set_value(dict_2dmixlp['Energy/power (J or W)'][2], 22.8e-3)
+results = []
+for energy in energies:
+    set_value(dict_2dmixlp['Energy/power (J or W)'][0], energy)
+    snlo.moveto(*dict_2dmixlp['Accept'])
+    #gui.moveTo(*coords(780, 800))
+    gui.click()
+    dict = run_2dmixlp()
+    m.set_key(dict, 'seed pulse energy (J)', energy)
+    dict2 = get_spectr()
+    results.append(merge_dict([dict, dict2]))
+    #gui.moveTo(*coords(600, 250))
+    snlo.moveto(*dict_2dmixlp['Change Inputs'])
+    gui.click()
+
+
+# save OPA2 seed dependence
+now = datetime.now()
+file_name = now.strftime("%Y-%m-%d") + '_' + 'OPA2_seed-dependence.pkl'
+pickle.dump(results, open(results_path + file_name, 'wb'))
+
+#%%% import OPA2 seed dependence
+file_name = '2022-06-09_OPA2_seed-dependence.pkl'
+results = merge_dict(pickle.load(open(results_path + file_name, 'rb')))
+
+x = np.array(results['seed pulse energy (J)']).T[0] * 1e6
+y = np.array(results['Output pulse energy (mJ)']).T[0][0]
+y2 = np.array(results['duration_idl_ns']).T[0]
+
+fig, [ax1, ax2] = plt.subplots(1, 2)
+m.plot_options(fig, ax1, xlabel='seed pulse energy ($\mathrm{\mu J}$)',
+               ylabel='idler pulse energy (mJ)')
+ax1.plot(x, y, 'b-')
+
+m.plot_options(fig, ax2, xlabel='seed pulse energy ($\mathrm{\mu J}$)',
+               ylabel='idler pulse duration (ns)')
+ax2.plot(x[y2<100], y2[y2<100], '-')
+plt.tight_layout()
+#m.save_plot('OPA2_seed_sim')
+plt.show()
+
+#%%% init 2D mix LP for second OPA2 crystal
+idler_nm = 3545
+signal_nm = 1520.6
+pump_nm = 1064.16
+seedenergy_j = 20e-6
+pump_j = 23e-3
+pumpduration_ns = 8
+fwhm_seed_mm = 3.9/2*1.18  # 1/e^2 radius to FWHM
+fwhm_pump_mm = 1.6/2*1.18  # 1/e^2 radius to FWHM
+
+values = [
+    [idler_nm, signal_nm, pump_nm],
+    [n_lnb(idler_nm, axis='o'), n_lnb(signal_nm, axis='o'), n(pm_theta(idler_nm, pump_nm), pump_nm, 25)],
+    [0, 0, 0],
+    [0.05, 0.005, 0.001],  # coating specs input, NC24, NC25
+    [0.05, 0.005, 0.001],  # coating specs output
+    [0, 0, 0],
+    [seedenergy_j, 1e-16, pump_j],
+    [pumpduration_ns, 0, pumpduration_ns],
+    [fwhm_seed_mm, fwhm_seed_mm, fwhm_pump_mm],
+    [1, 1, 1],  # supergaussian coefficient
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, -37.14],  # walkoff
+    [0, 0, 1.123],  # offset
+    [1631400, 9172700, 515780],  # todo: radius of curvature
+    [50, 32, 32],  # integration points
+    [30, 8, 8], # crystal length, grid
+    [-4.04],  # via QPM, LNB_M
+    [0],  # delta k
+    [1000],  # distance to image
+    [25]  # time steps
+]
+# variables: idler-energy, signal-energy, pump-energy, pulse-duratios (?)
+
+#%%% full input 2D-mix-LP - OPA2
+open_function('2D-mix-LP')
+for index in range(len(values)):
+    if len(values[index]) == 1:
+        pos = list(dict_2dmixlp.values())[index][0]
+        set_value(pos, values[index][0])
+    else:
+        for element in range(len(values[index])):
+            set_value(list(dict_2dmixlp.values())[index][element], values[index][element])
+
+#%% OPA2 pump energy dependence second crystal
+energies = np.arange(1e-6, 25e-3, 250e-6)
+print(energies)
+print(len(energies))
+# need pump dependence for first crystal
+file_name = '2022-06-15_OPA2_pump-dependence.pkl'
+result_xtal1 = pickle.load(open(results_path+ file_name, 'rb'))
+
+#print(np.array(result_xtal1[0]['Output pulse energy (mJ)']))
+
+#%%% run OPA2 pump energy dependence for the second crystal
+results = []
+for energy in energies:
+    index = list(energies).index(energy)
+    pump_energy = result_xtal1[index]['pump pulse energy (J)'][0]
+    input_energies = np.array(result_xtal1[index]['Output pulse energy (mJ)'])[0]
+    input_durations = result_xtal1[index]['duration_idl_ns'], result_xtal1[index]['duration_sig_ns'], result_xtal1[index][
+        'duration_pmp_ns']
+    
+    # ignore the rest, if the Gaussian fit failed after the first stage
+    if input_durations[0][0]>100 or input_durations[1][0]>100 or input_durations[2][0]>100:
+        continue
+    
+    set_value(dict_2dmixlp['Energy/power (J or W)'][0], 1e-3*input_energies[0])
+    set_value(dict_2dmixlp['Energy/power (J or W)'][1], 1e-3*input_energies[1])
+    set_value(dict_2dmixlp['Energy/power (J or W)'][2], 1e-3*input_energies[2])
+
+    set_value(dict_2dmixlp['Pulse duration (fwhm ns)'][0], input_durations[0][0])
+    set_value(dict_2dmixlp['Pulse duration (fwhm ns)'][1], input_durations[1][0])
+    set_value(dict_2dmixlp['Pulse duration (fwhm ns)'][2], input_durations[2][0])
+
+    snlo.moveto(*dict_2dmixlp['Accept'])
+    gui.click()
+    dict = run_2dmixlp()
+    m.set_key(dict, 'pump pulse energy (J)', energy)
+    m.set_key(dict, 'xtal1 pump pulse energy (J)', pump_energy)
+    dict2 = get_spectr()
+    results.append(merge_dict([dict, dict2]))
+    snlo.moveto(*dict_2dmixlp['Change Inputs'])
+    gui.click()
+
+# save OPA2 pump dependence
+now = datetime.now()
+file_name = now.strftime("%Y-%m-%d") + '_' + 'OPA2_pump-dependence_xtal2.pkl'
+pickle.dump(results, open(results_path + file_name, 'wb'))
+
+
+#%%% plot results OPA2 pump dep
+energies = np.arange(1e-6, 25e-3, 250e-6)
+
+file_name = '2022-06-09_OPA2_pump-dependence.pkl'
+#results = merge_dict(pickle.load(open('simulation/' + file_name, 'rb')))
+results = merge_dict(pickle.load(open(results_path + file_name, 'rb')))
+
+file_name = '2022-06-15_OPA2_pump-dependence_xtal2.pkl'
+#results2 = merge_dict(pickle.load(open('simulation/' + file_name, 'rb')))
+results2 = merge_dict(pickle.load(open(results_path + file_name, 'rb')))
+
+
+x1 = np.array(results['pump pulse energy (J)']).T[0]*1e3
+y1 = np.array(results['Output pulse energy (mJ)']).T[0][0]
+x2 = np.array(results2['pump pulse energy (J)']).T[0]*1e3
+y2 = np.array(results2['Output pulse energy (mJ)']).T[0][0]
+
+
+fig, [ax1, ax2] = plt.subplots(1, 2)
+m.plot_options(fig, ax1, xlabel='pump pulse energy (mJ)', ylabel='idler pulse energy (mJ)')
+ax1.plot(x1, y1, 'b--', label='only first crystal')
+ax1.plot(x2, y2, 'r-', label='both crystals')
+ax1.legend(loc='upper left')
+
+y1 = np.array(results['duration_idl_ns']).T[0]
+y2 = np.array(results2['duration_idl_ns']).T[0]
+
+limit = 10
+m.plot_options(fig, ax2, xlabel='pump pulse energy (mJ)', ylabel='idler pulse duration (ns)')
+ax2.plot(x1[y1<limit], y1[y1<limit], 'b--', label='only first crystal')
+ax2.plot(x2[y2<limit], y2[y2<limit], 'r-', label='both crystals')
+ax2.set_ylim([-1, 10])
+ax2.legend(loc='upper left')
+plt.tight_layout()
+#m.save_plot('OPA2_pump_sim_full')
+plt.show()
+
+#%% OPA2 seed energy dependence second crystal
+energies = np.arange(1e-7, 40e-6, 0.5e-6)
+print(energies)
+print(len(energies))
+# need pump dependence for first crystal
+file_name = '2022-06-15_OPA2_seed-dependence.pkl'
+#result_xtal1 = pickle.load(open('simulation/' + file_name, 'rb'))
+result_xtal1 = pickle.load(open(results_path + file_name, 'rb'))
+
+#print(np.array(result_xtal1[0]['Output pulse energy (mJ)']))
+
+#%%% run OPA2 seed energy dependence for the second crystal
+results = []
+for energy in energies:
+    index = list(energies).index(energy)
+    seed_energy = result_xtal1[index]['seed pulse energy (J)'][0]
+    input_energies = np.array(result_xtal1[index]['Output pulse energy (mJ)'])[0]
+    input_durations = result_xtal1[index]['duration_idl_ns'], result_xtal1[index]['duration_sig_ns'], \
+                      result_xtal1[index][
+                          'duration_pmp_ns']
+
+    # ignore the rest, if the Gaussian fit failed after the first stage
+    if input_durations[0][0]>100 or input_durations[1][0]>100 or input_durations[2][0]>100:
+        continue
+    
+    set_value(dict_2dmixlp['Energy/power (J or W)'][0], 1e-3 * input_energies[0])
+    set_value(dict_2dmixlp['Energy/power (J or W)'][1], 1e-3 * input_energies[1])
+    set_value(dict_2dmixlp['Energy/power (J or W)'][2], 1e-3 * input_energies[2])
+
+    set_value(dict_2dmixlp['Pulse duration (fwhm ns)'][0], input_durations[0][0])
+    set_value(dict_2dmixlp['Pulse duration (fwhm ns)'][1], input_durations[1][0])
+    set_value(dict_2dmixlp['Pulse duration (fwhm ns)'][2], input_durations[2][0])
+    snlo.moveto(*dict_2dmixlp['Accept'])
+    gui.click()
+    dict = run_2dmixlp()
+    m.set_key(dict, 'seed pulse energy (J)', energy)
+    m.set_key(dict, 'xtal1 seed pulse energy (J)', seed_energy)
+    dict2 = get_spectr()
+    results.append(merge_dict([dict, dict2]))
+    snlo.moveto(*dict_2dmixlp['Change Inputs'])
+    gui.click()
+
+# save OPA2 seed dependence
+now = datetime.now()
+file_name = now.strftime("%Y-%m-%d") + '_' + 'OPA2_seed-dependence_xtal2.pkl'
+#pickle.dump(results, open('simulation/' + file_name, 'wb'))
+pickle.dump(results, open(results_path + file_name, 'wb'))
+
+#%%% plot results
+file_name = '2022-06-09_OPA2_seed-dependence.pkl'
+results = merge_dict(pickle.load(open(results_path + file_name, 'rb')))
+
+file_name = '2022-06-15_OPA2_seed-dependence_xtal2.pkl'
+results2 = merge_dict(pickle.load(open(results_path + file_name, 'rb')))
+
+#x = np.array(result['seed pulse energy (J)']).T[0]*1e6
+#y1 = np.array(result['Output pulse energy (mJ)']).T[0][0]
+#y2 = np.array(result2['Output pulse energy (mJ)']).T[0][0]
+
+x1 = np.array(results['seed pulse energy (J)']).T[0]*1e6
+y1 = np.array(results['Output pulse energy (mJ)']).T[0][0]
+x2 = np.array(results2['seed pulse energy (J)']).T[0]*1e6
+y2 = np.array(results2['Output pulse energy (mJ)']).T[0][0]
+
+
+fig, [ax1, ax2] = plt.subplots(1, 2)
+m.plot_options(fig, ax1, xlabel='seed pulse energy ($\mathrm{\mu J}$)', ylabel='idler pulse energy (mJ)')
+ax1.plot(x1, y1, 'b--', label='only first crystal')
+ax1.plot(x2, y2, 'r-', label='both crystals')
+ax1.legend(loc='upper left')
+
+limit = 10
+m.plot_options(fig, ax2, xlabel='seed pulse energy ($\mathrm{\mu J}$)', ylabel='idler pulse duration (ns)')
+ax2.plot(x1[y1<limit], y1[y1<limit], 'b--', label='only first crystal')
+ax2.plot(x2[y2<limit], y2[y2<limit], 'r-', label='both crystals')
+
+ax2.legend(loc='upper left')
+
+plt.tight_layout()
+#m.save_plot('OPA2_seed_sim_full', path=plots_path)
+plt.show()
+
+#%% conversion efficiency OPA2
+file_name = '2022-06-15_OPA2_pump-dependence_xtal2.pkl'
+results2 = merge_dict(pickle.load(open(results_path + file_name, 'rb')))
+
+x = np.array(results2['pump pulse energy (J)']).T[0]*1e3
+y = np.array(results2['Output pulse energy (mJ)']).T
+eff = ((y[0]+y[1])/x)[0]
+
+fig, ax1 = plt.subplots()
+m.plot_options(fig, ax1, image_width=15, xlabel='pump pulse energy (mJ)', ylabel='pulse energy ($\mathrm{\mu J}$)')
+ax1.plot(x, y[0][0], 'r-', label='idler pulse energy')
+ax1.plot(x, y[1][0], 'g-', label='signal pulse energy')
+ax1.plot(x, y[2][0], 'b-', label='residual pump pulse energy')
+ax1.legend(loc='upper left')
+
+ax12 = ax1.twinx()
+m.plot_options(fig, ax12, image_width=15, aspect_ratio=2)
+ax12.plot(x[eff<1], eff[eff<1], 'k--', label='conversion efficiency')
+ax12.set_ylabel('conversion efficiency')
+
+ax12.legend(loc='lower right')
+
+plt.tight_layout()
+#m.save_plot('OPA1_pump_sim')
+plt.show()
\ No newline at end of file
diff --git a/snlo-helper/__init__.py b/snlo-helper/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/snlo-helper/lnbVsAgs.py b/snlo-helper/lnbVsAgs.py
new file mode 100644
index 0000000000000000000000000000000000000000..63eb2a902263dede0dc055af1714b1e162f5d998
--- /dev/null
+++ b/snlo-helper/lnbVsAgs.py
@@ -0,0 +1,104 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu May 25 13:01:56 2023
+
+@author: moneke
+"""
+
+# %% imports
+
+from time import perf_counter
+
+import numpy as np
+import matplotlib as plt
+
+from snlohelper import *
+
+# %% Setup
+
+
+sim = TwoDMixLP()
+sim.open()
+
+start = sim.get_configuration()
+
+
+# %%%
+
+# Parameters
+absorption = np.linspace(0, 41e-1, 30)
+
+# Loop
+results = []
+for a in absorption:
+    r = sim.configure_run_read({'Crystal loss (1/mm)': [a, None, None]})
+    results.append(r['Output pulse energy (mJ)'][0])
+    print(a, results[-1])
+
+
+# %% Test
+
+interval = 0#.5
+s0 = perf_counter()
+gui.click(*scale(133,293))
+gui.hotkey("ctrl", "home")
+gui.hotkey("ctrl", "shiftleft", "shiftright", "end")
+# gui.hotkey('ctrl', 'c')
+
+print(perf_counter()-s0)
+
+
+
+s0 = perf_counter()
+gui.rightClick(*scale(133,293))
+gui.press(6 * ['down'])
+gui.press('enter')
+# gui.hotkey('ctrl', 'c')
+
+print(perf_counter()-s0)
+
+
+
+# %%
+
+
+with gui.hold("shift", interval=interval):
+    with gui.hold("ctrl", interval=interval):
+        gui.press("end", interval=interval)
+# 
+
+
+#%%
+gui.hotkey("shift", "a")
+
+
+# %%
+
+results = []
+
+#%%
+
+sim.open()
+
+for data in (ags_o, ags_e):
+    results.append(sim.configure_run_read(data))
+
+
+
+# %%
+
+pe = np.linspace(0, 0.20, 40)
+for v in pe:
+    results.append(sim.configure_run_read({'Energy/power (J or W)': [None, None, v]})['Output pulse energy (mJ)'])
+
+
+# %%
+
+plt.figure()
+ax = plt.gca()
+ax.plot(pe, [r[0] for r in results], label="losless")
+ax.plot(pe, [0.00103, 0.00331, 0.00869, 0.0196, 0.0397, 0.074, 0.129, 0.213, 0.333, 0.494, 0.696, 0.934, 1.2, 1.47, 1.74, 2.0, 2.24, 2.46, 2.66, 2.84, 3.0, 3.14, 3.27, 3.4, 3.51, 3.62, 3.72, 3.82, 3.92, 4.03, 4.13, 4.24, 4.36, 4.49, 4.62, 4.77, 4.92, 5.09, 5.27, 5.46], label="absorption")
+ax.set_xlabel("Pumpenergie in J")
+ax.set_ylabel("MIR-Energie in mJ")
+ax.legend()
+
diff --git a/snlo-helper/plots/.gitkeep b/snlo-helper/plots/.gitkeep
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/snlo-helper/results/.gitkeep b/snlo-helper/results/.gitkeep
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/snlo-helper/snlohelper.py b/snlo-helper/snlohelper.py
new file mode 100644
index 0000000000000000000000000000000000000000..b67af3ef4c3f95cf79ca11ac63f88004893c5ba4
--- /dev/null
+++ b/snlo-helper/snlohelper.py
@@ -0,0 +1,503 @@
+# -*- coding: utf-8 -*-
+"""
+Helper for SNLO simulations ("autoclicker")
+
+This file clicks automatically on buttons of SNLO program.
+
+
+created on 25.05.2022 by Jan Frederic Kinder
+"""
+
+import logging
+import time
+import re
+
+import pyautogui as gui
+import matplotlib.pyplot as plt
+import numpy as np
+from pyperclip import paste
+
+import analysis.modules_JFK as m
+
+
+log = logging.getLogger(__name__)
+log.addHandler(logging.NullHandler())
+
+
+"""
+Setup of the screen and scaling.
+
+All positions in code are given on a Full HD (1920 * 1080) screen and dynamically adjusted to the
+screen resolution.
+"""
+
+
+def get_screenfactors(standard=[1920, 1080]):
+    """Get the scaling factor from Full HD to the current display resolution."""
+    width, height = gui.size()
+    return standard[0] / width, standard[1] / height
+
+
+def scale(x: float, y: float = None) -> tuple:
+    """Scale coordinates from the definition standard to the current screen."""
+    global factors
+    if y is None and isinstance(x, (list, tuple)):
+        x, y = x
+    return x / factors[0], y / factors[1]
+
+
+def standard_position() -> tuple:
+    """Get the mouse position in standard coordinates (x, y)."""
+    point = gui.position()
+    global factors
+    return point.x * factors[0], point.y * factors[1]
+
+
+"""
+Begin of dictionaries/lists etc.
+"""
+
+
+# coordinates of the functions (in FHD standard)
+_functions_coord = {
+    'Ref. Ind.': [66, 46],
+    'Qmix': [66, 66],
+    'Bmix': [66, 93],
+    'QPM': [66, 120],
+    'Opoangles': [66, 146],
+    'Ncpm': [66, 173],
+    'GVM': [66, 200],
+    'PW-mix-LP': [66, 233],
+    'PW-mix-SP': [66, 260],
+    'PW-mix-BB': [66, 286],
+    '2D-mix-LP': [66, 313],
+    '2D-mix-SP': [66, 340],
+    'PW-cav-LP': [66, 366],
+    'PW-OPO-SP': [66, 393],
+    'PW-OPO-BB': [66, 420],
+    '2D-cav-LP': [66, 446],
+    'Focus': [66, 473],
+    'Cavity': [66, 500],
+}
+
+
+# coordinates of the Focus-function (in FHD standard)
+off = 200, 200
+_dict_focus = {
+    'Wavelength (nm)': [366, off[1] + 20],
+    'Refractive Index': [366, off[1] + 40],
+    'Waist size (mm)': [366, off[1] + 60],
+    'Face to focus (mm)': [366, off[1] + 90],
+    'Dist. to focus (mm)': [366, 300],
+    'Rayleigh z in xtal (mm)': [366, 353],
+    'Beam size (mm)': [366, 380],
+    'Radius of curv. (mm)': [366, 393],
+    'Far field ang air (mrad)': [366, 413],
+    'fwhm': [219, 216],
+    '1e^2': [166, 216]
+}
+
+
+# coordinates of the 2DmixLP-function (in FHD standard)
+_configuration_2dmixlp = {
+    'Wavelengths (nm)': [[403, 186], [470, 186], [536, 186]],
+    'Indexes of refraction': [[403, 200], [470, 200], [536, 200]],
+    'Phases at input (rad)': [[403, 213], [470, 213], [536, 213]],
+    'Input face reflectivity (0-1)': [[403, 226], [470, 226], [536, 226]],
+    'Output face reflectivity (0-1)': [[403, 246], [470, 246], [536, 246]],
+    'Crystal loss (1/mm)': [[403, 260], [470, 260], [536, 260]],
+    'Energy/power (J or W)': [[416, 280], [483, 280], [550, 280]],
+    'Pulse duration (fwhm ns)': [[403, 293], [470, 293], [536, 293]],
+    'Beam diam. (fwhm mm)': [[403, 306], [470, 306], [536, 306]],
+    'Supergaussian coeff.': [[403, 326], [470, 326], [536, 326]],
+    'n2 red1 (sq cm/W)': [[403, 340], [470, 340], [536, 340]],
+    'n2 red2 (sq cm/W)': [[403, 360], [470, 360], [536, 360]],
+    'n2 blue (sq cm/W)': [[403, 373], [470, 373], [536, 373]],
+    'beta red1 (cm/W)': [[403, 393], [470, 393], [536, 393]],
+    'beta red2 (cm/W)': [[403, 406], [470, 406], [536, 406]],
+    'beta blue (cm/W)': [[403, 420], [470, 420], [536, 420]],
+    'Walkoff angles (mrad)': [[403, 440], [470, 440], [536, 440]],
+    'Offset in wo dir. (mm)': [[403, 453], [470, 453], [536, 453]],
+    'Rad. curv. (mm/air“)': [[403, 473], [470, 473], [536, 473]],
+    '# of integ/grid points': [[403, 486], [470, 486], [536, 486]],
+    'Crystal/grid sizes (mm)': [[403, 500], [470, 500], [536, 500]],
+    'Deff (pm/V)': [[403, 520]],
+    'Delta k (1/mm)': [[403, 533]],
+    'Dist. to image (mm)': [[403, 546]],
+    '# time steps': [[403, 566]],
+}
+
+
+_buttons_2dmixlp = {
+    'Accept': [506, 536],
+    'Run': [140, 220],
+    'Change Inputs': [373, 166],
+}
+
+
+# TODO: remove following three definitions from this file (already in modules_JFK)
+def n_lnb(lambda_nm, temp_degc=25, axis='o'):
+    """refractive index of 5% MgO doped congruent LNB (supplied by HC Photonics)
+    data and functions from Gayer et al. (Appl. Optics, 2008)"""
+    if axis == "e":
+        a1, a2, a3, a4, a5, a6 = 5.756, 0.0983, 0.2020, 189.32, 12.52, 1.32e-2
+        b1, b2, b3, b4 = 2.86e-6, 4.7e-8, 6.113e-8, 1.516e-4
+    elif axis == "o":
+        a1, a2, a3, a4, a5, a6 = 5.653, 0.1185, 0.2091, 89.61, 10.85, 1.97e-2
+        b1, b2, b3, b4 = 7.941e-7, 3.134e-8, -4.641e-9, -2.188e-6
+    f = (np.array(temp_degc) - 24.5) * (np.array(temp_degc) + 570.82)
+    lambda_um = np.array(lambda_nm) * 1e-3
+    # Sellmeier equation for lambda in micrometers!
+    return np.sqrt(
+        a1 + b1 * f + (a2 + b2 * f) / (lambda_um ** 2 - (a3 + b3 * f) ** 2)
+        + (a4 + b4 * f) / (lambda_um ** 2 - a5 ** 2) - a6 * lambda_um ** 2)
+
+
+def n_bbo(lambda_nm, axis='o'):
+    """refractive index of BBO
+    data and functions from Zhang et al. (Opt. Comm., 2000)
+    => no temperature dependence!"""
+    if axis == "o":
+        a1, a2, a3, a4, a5, a6 = 2.7359, 0.01878, 0.01822, 0.01471, 0.0006081, 0.00006740
+    elif axis == "e":
+        a1, a2, a3, a4, a5, a6 = 2.3753, 0.01224, 0.01667, 0.01627, 0.0005716, 0.00006305
+    lambda_um = np.array(lambda_nm) * 1e-3
+    # Sellmeier equation for lambda in micrometers!
+    return np.sqrt(
+        a1 + a2 / (lambda_um ** 2 - a3) - a4 * lambda_um ** 2 + a5 * lambda_um ** 4
+        - a6 * lambda_um ** 6)
+
+
+def n(theta, lambda_nm, temp_degc=25, crystal='LNB'):
+    """Return the refractive index under propagation with angle of theta with the crystal c-axis."""
+    if crystal == 'LNB':
+        return (1 / (np.sqrt(
+            np.cos(np.deg2rad(np.array(theta))) ** 2
+            / n_lnb(np.array(lambda_nm), np.array(temp_degc), 'o') ** 2
+            + np.sin(np.deg2rad(np.array(theta))) ** 2
+            / n_lnb(np.array(lambda_nm), np.array(temp_degc), 'e') ** 2)))
+    elif crystal == 'BBO':
+        return (1 / (np.sqrt(
+            np.cos(np.deg2rad(np.array(theta))) ** 2
+            / n_bbo(np.array(lambda_nm), 'o') ** 2
+            + np.sin(np.deg2rad(np.array(theta))) ** 2
+            / n_bbo(np.array(lambda_nm), 'e') ** 2)))
+
+    else:
+        print('no Sellmeier equation for ' + crystal + '!')
+
+
+def pm_theta(idler_nm, pump_nm, temp_degc=25, crystal='LNB'):
+    """ returns the angle, which provides critical phase matching in MgO:LNB"""
+    signal_nm = m.idler2signal(idler_nm, pump_nm)
+    if crystal == 'LNB':
+        n1 = n_lnb(idler_nm, temp_degc, 'o')
+        n2 = n_lnb(signal_nm, temp_degc, 'o')
+        n_pump_e = n_lnb(pump_nm, temp_degc, 'e')
+        n_pump_o = n_lnb(pump_nm, temp_degc, 'o')
+    elif crystal == 'BBO':
+        n1 = n_bbo(idler_nm, 'o')
+        n2 = n_bbo(signal_nm, 'o')
+        n_pump_e = n_bbo(pump_nm, 'e')
+        n_pump_o = n_bbo(pump_nm, 'o')
+    n_soll = (n1 / idler_nm + n2 / signal_nm) * pump_nm
+
+    return 180 / np.pi * np.arctan(
+        n_pump_e / n_pump_o
+        * np.sqrt((n_soll ** 2 - n_pump_o ** 2) / (n_pump_e ** 2 - n_soll ** 2)))
+
+
+def pm_theta_bbo(red1_nm, red2_nm):
+    """ returns the angle, which provides critical phase matching in MgO:LNB"""
+    blue_nm = 1 / (1 / red1_nm + 1 / red2_nm)
+    n1 = n_bbo(red1_nm, 'o')
+    n2 = n_bbo(red2_nm, 'e')
+    n_blue_e = n_bbo(blue_nm, 'e')
+    n_blue_o = n_bbo(blue_nm, 'o')
+    n_soll = (n1 / red1_nm + n2 / red2_nm) * blue_nm
+
+    return 180 / np.pi * np.arctan(
+        n_blue_e / n_blue_o
+        * np.sqrt((n_soll ** 2 - n_blue_o ** 2) / (n_blue_e ** 2 - n_soll ** 2)))
+
+
+# Helper functions
+def get_content(position: tuple) -> str:
+    """Get the content of the field at position via double click.
+
+    If there is a "-" in the text, the extraction fails!
+    """
+    gui.doubleClick(*scale(*position))
+    gui.hotkey('ctrl', 'c')
+    return paste()
+
+
+def get_value(position: tuple) -> float:
+    """move to position, retrieve value and return float"""
+    return float(get_content(position))
+
+
+def get_content_complete(position: tuple) -> str:
+    """Go to position and retrieve the content there, marking all."""
+    gui.click(*scale(*position))
+    gui.hotkey("ctrl", "home")
+    gui.hotkey("ctrl", "shiftleft", "shiftright", "end")  # both shift keys necessary if keylock on.
+    gui.hotkey('ctrl', 'c')
+    return paste()
+
+
+def get_value_complete(position: tuple) -> float:
+    """moves to position, retrieves value via context menu (slower) and returns float"""
+    return float(get_content_complete(position))
+
+
+def set_value(position: tuple, value) -> None:
+    """move to position, insert value as string"""
+    gui.doubleClick(*scale(*position))
+    gui.press('delete')
+    gui.doubleClick()
+    gui.write(str(value))
+
+
+# SNLO calls
+def open_function(key: str) -> None:
+    """opens function according to key"""
+    gui.click(*scale(*_functions_coord[key]))
+
+
+def focus(wavelength_nm: float, ref_index: float, fwhm_mm: float, focus_pos_mm: float) -> tuple:
+    """calls focus function to determine parameters for the radius of curvature.
+    inputs: wavelength_nm, ref_index, waist_size_mm (FWHM), focus_pos_mm
+    returns: zr, diameter (FWHM), radcurv, angle
+    """
+    open_function('Focus')
+    gui.leftClick(*scale(*_dict_focus['fwhm']))
+    set_value(_dict_focus['Wavelength (nm)'], wavelength_nm)
+    set_value(_dict_focus['Refractive Index'], ref_index)
+    set_value(_dict_focus['Waist size (mm)'], fwhm_mm)
+    set_value(_dict_focus['Face to focus (mm)'], focus_pos_mm)
+    # Setting 'Dist. to focus (mm)' equal to 'Face to focus (mm)' gives parameters in air at the
+    # input face
+    set_value(_dict_focus['Dist. to focus (mm)'], focus_pos_mm)
+    # readout
+    zr = get_value_complete(_dict_focus['Rayleigh z in xtal (mm)'])
+    diameter = get_value_complete(_dict_focus['Beam size (mm)'])
+    radcurv = get_value_complete(_dict_focus['Radius of curv. (mm)'])
+    angle = get_value_complete(_dict_focus['Far field ang air (mrad)'])
+    return zr, diameter, radcurv, angle
+
+
+class TwoDMixLP:
+    """The '2D-mix-LP' method."""
+
+    def open(self):
+        """Open the function."""
+        open_function("2D-mix-LP")
+
+    def accept(self):
+        """Click 'Accept'."""
+        gui.click(*scale(506, 536))
+
+    def run(self):
+        """Click 'Run'."""
+        gui.click(*scale(140, 220))
+
+    def change_inputs(self):
+        """Click 'Change Inputs'."""
+        gui.click(*scale(373, 166))
+
+    def configure(self, data: dict = None):
+        """Configure the values and leave the config window open."""
+        open_function("2D-mix-LP")
+        if data is None:
+            return
+        for key, value in data.items():
+            positions = _configuration_2dmixlp[key]
+            for i, val in enumerate(value):
+                if val is not None:
+                    set_value(positions[i], val)
+
+    def get_configuration(self) -> dict:
+        data = {}
+        for key, positions in _configuration_2dmixlp.items():
+            d = []
+            for pos in positions:
+                val = get_content_complete(pos)
+                try:
+                    d.append(float(val))
+                except ValueError:
+                    d.append(val)
+            data[key] = d
+        return data
+
+    def run_and_read(self, waiting_time: float = 1, max_tries: int = 10, interval: float = 0.5
+                     ) -> dict:
+        """Run an analysis and return the result."""
+        gui.click(*scale(*_buttons_2dmixlp['Run']))
+        time.sleep(waiting_time)
+
+        for _ in range(max_tries):
+            entries = get_content_complete((133, 293)).split("\r\n")
+            if len(entries) > 3:
+                break
+            time.sleep(interval)
+
+        # interpret results and save as dictionary:
+        return {
+            'Input peak irradiance (W/sq cm)': [float(i) for i in entries[0].split()[5:]],
+            'Input peak fluence (J/sq cm)': [float(i) for i in entries[1].split()[6:]],
+            'Input peak powers (W)': [float(i) for i in entries[2].split()[5:]],
+            'Output peak fluence (J/sq cm)': [float(i) for i in entries[3].split()[6:]],
+            'Output pulse energy (mJ)': [float(i) for i in entries[4].split()[5:]],
+            'So (W/sq cm)': float(entries[5].split()[4])
+        }
+
+    def configure_run_read(self, data: dict = None, **kwargs) -> dict:
+        """Configure and run an analysis and return the result."""
+        self.configure(data)
+        self.accept()
+        return self.run_and_read(**kwargs)
+
+
+def import_snlo_file(file_name, file_path='C:/SNLO/'):
+    """import a file generated by SNLO, specified by filename and return array (no header)"""
+    with open(file_path + file_name, 'r') as f:
+        input = f.readlines()
+    data = []
+    e_numbers = re.compile(r"-*[\d.]*?E-*[+]*\d*")
+    r""" Regex
+        -*        "-" or nothing
+        [\d.]*?   decimal number (0 or more characters)
+        -*        "-" or nothing
+        [+]*      "+" or nothing
+        \d*       decimal number (0 or more characters)
+    """
+    for i in input:
+        temp = []
+        # old implementation, suffers from strings without whitespace, e.g.:
+        #    '9.206897E-100-1.616264E-2'
+        #    element = i.split()
+        #    for item in element:
+        #        temp.append(float(item))
+        new_element = e_numbers.findall(i)
+        for item in new_element:
+            temp.append(float(item))
+        data.append(temp)
+    return np.array(data.copy())
+
+
+def get_spectr(call_function: bool = True, display_results: bool = False) -> dict:
+    """imports spectra and power vs. time data for all three fields of the mixing process.
+    call_function: for 2D-mix-LP you have to click on 'Spectra' to calculate the spectra.
+    display_results: gives plots for spectra and powers
+    returns dictionary with bandwidths and durations (given as FWHMs)"""
+    if call_function:
+        gui.click(*scale(400, 260))
+        time.sleep(0.5)
+
+    # detuning [MHz], Red1, Red2, Blue
+    spectr = import_snlo_file('OPA2D_SP.dat')
+    # time, power, phase, Mx^2, My^2, Rad Curv c, Rad Curv y, X-tilt, w_x^2, w_y^2
+    id_beam = import_snlo_file('ID_BEAM.dat')
+    sig_beam = import_snlo_file('SIG_BEAM.dat')
+    pmp_beam = import_snlo_file('PMP_BEAM.dat')
+
+    # fit idler duration
+    x = 1e9 * id_beam.T[0]  # ns
+    y = id_beam.T[1]
+    idl = m.gaussian_fit(np.array([x, y]).T, False)[0]
+    fwhm_ns_idl = abs(2 * np.sqrt(2 * np.log(2)) * idl[2])
+    x_lists = [x.copy()]
+    y_lists = [y.copy()]
+
+    # fit signal duration
+    x = 1e9 * sig_beam.T[0]  # ns
+    y = sig_beam.T[1]
+    sig = m.gaussian_fit(np.array([x, y]).T, False)[0]
+    fwhm_ns_sig = abs(2 * np.sqrt(2 * np.log(2)) * sig[2])
+    x_lists.append(x)
+    y_lists.append(y)
+
+    # fit pump duration
+    x = 1e9 * pmp_beam.T[0]  # ns
+    y = pmp_beam.T[1]
+    pmp = m.gaussian_fit(np.array([x, y]).T, False)[0]
+    fwhm_ns_pmp = abs(2 * np.sqrt(2 * np.log(2)) * pmp[2])
+    x_lists.append(x)
+    y_lists.append(y)
+    fits = idl, sig, pmp
+    if display_results:
+        # plot pulse durations
+        fig, ax = plt.subplots(1, 3)
+        fwhms = fwhm_ns_idl, fwhm_ns_sig, fwhm_ns_pmp
+        titles = 'idler (Red1)', 'signal (Red2)', 'pump (Blue)'
+        colors = 'red', 'green', 'blue'
+        for i in range(len(ax)):
+            m.plot_options(fig, ax[i], xlabel='time (ns)', image_width=15, aspect_ratio=3)
+            ax[i].plot(x_lists[i], y_lists[i], '.', color=colors[i])
+            ax[i].set_title(titles[i], color=colors[i])
+
+            x_fit = np.arange(min(x_lists[i]), max(x_lists[i]), 0.1)
+
+            ax[i].plot(x_fit, m.gauss(x_fit, *fits[i]), '-', color=colors[i])
+            ax[i].text(0, 0, f'FWHM\n{fwhms[i]:.1f}ns', ha='center', color=colors[i])
+            # ax[i].set_xlim([min(x_lists[i]), max(x_lists[i])])
+        m.plot_options(fig, ax[0], image_width=15, aspect_ratio=3, xlabel='time (ns)',
+                       ylabel='normalized power')
+        plt.tight_layout()
+        plt.show()
+
+    # fit spectra
+    x = spectr.T[0]  # in MHz
+    y_idl = spectr.T[1]
+    y_sig = spectr.T[2]
+    y_pmp = spectr.T[3]
+
+    idl = m.gaussian_fit(np.array([x, y_idl]).T, False)[0]
+    sig = m.gaussian_fit(np.array([x, y_sig]).T, False)[0]
+    pmp = m.gaussian_fit(np.array([x, y_pmp]).T, False)[0]
+
+    fwhm_spec_idl = abs(2 * np.sqrt(2 * np.log(2)) * idl[2])
+    fwhm_spec_sig = abs(2 * np.sqrt(2 * np.log(2)) * sig[2])
+    fwhm_spec_pmp = abs(2 * np.sqrt(2 * np.log(2)) * pmp[2])
+
+    y_lists = y_idl, y_sig, y_pmp
+
+    fits = idl, sig, pmp
+    if display_results:
+        # plot spectra durations
+        fig, ax = plt.subplots(1, 3)
+        fwhms = fwhm_spec_idl, fwhm_spec_sig, fwhm_spec_pmp
+        titles = 'idler (Red1)', 'signal (Red2)', 'pump (Blue)'
+        colors = 'red', 'green', 'blue'
+        for i in range(len(ax)):
+            m.plot_options(fig, ax[i], xlabel='detuning (MHz)', image_width=15, aspect_ratio=3)
+            ax[i].plot(x, y_lists[i], '.', color=colors[i])
+            ax[i].set_title(titles[i], color=colors[i])
+
+            x_fit = np.arange(min(x), max(x), 0.1)
+
+            ax[i].plot(x_fit, m.gauss(x_fit, *fits[i]), '-', color=colors[i])
+            ax[i].text(0, 0, f'FWHM\n{fwhms[i]:.1f}MHz', ha='center', color=colors[i])
+            ax[i].set_xlim([-2 * fwhms[i], 2 * fwhms[i]])
+        m.plot_options(fig, ax[0], image_width=15, aspect_ratio=3, xlabel='time (ns)',
+                       ylabel='normalized power')
+        plt.tight_layout()
+        plt.show()
+
+    return {
+        'fwhm_idl_MHz': fwhm_spec_idl,
+        'fwhm_sig_MHz': fwhm_spec_sig,
+        'fwhm_pmp_MHz': fwhm_spec_pmp,
+        'duration_idl_ns': fwhm_ns_idl,
+        'duration_sig_ns': fwhm_ns_sig,
+        'duration_pmp_ns': fwhm_ns_pmp
+    }
+
+
+if __name__ == "__main__":
+    if len(log.handlers) < 2:
+        log.addHandler(logging.StreamHandler())
+    log.setLevel(logging.INFO)
+    factors = get_screenfactors()
+    log.info(f"Setting screenfactors to {factors}.")
diff --git a/snlo-helper/upconversion.py b/snlo-helper/upconversion.py
new file mode 100644
index 0000000000000000000000000000000000000000..37030481e5d26016a02e597f6e5714d8cef6869e
--- /dev/null
+++ b/snlo-helper/upconversion.py
@@ -0,0 +1,355 @@
+# -*- coding: utf-8 -*-
+"""
+
+created on 25.05.2022 by Jan Frederic Kinder
+"""
+import time
+import pickle
+import pyautogui as gui
+import matplotlib.pyplot as plt
+import modules_JFK as m  # module from NLOQO/EG-Labor/DataAnalysis, check PYTHONPATH in Spyder
+import numpy as np
+import scipy
+from pyperclip import paste
+from socket import gethostname
+from collections import defaultdict
+from datetime import datetime  # get current time
+
+import snlohelper as snlo
+# imports all of snlohelper, risk of namespace-confusion...
+from snlohelper import *
+
+"""define some global filepath variables based upon the computer name:"""
+def get_results_path():
+    pc_name = gethostname()
+    if pc_name == 'WhiteShard':
+        return 'M:/'
+    elif pc_name == 'POPS':
+        return 'C:/Users/kinder.NLOQO/HESSENBOX-DA/OPA-Paper/analysis (python)/simulation/'
+    elif pc_name == 'Myres':
+        return 'D:/'
+    else:
+        print('data location not specified, use generic path: M:/')
+        return 'M:/'
+    
+def get_plots_path():
+    pc_name = gethostname()
+    if pc_name == 'WhiteShard':
+        return 'M:/'
+    elif pc_name == 'POPS':
+        return f'C:/Users/kinder.NLOQO/HESSENBOX-DA/OPA-Paper/analysis (python)/Ausgabe/'
+    elif pc_name == 'Myres':
+        return 'D:/'
+    else:
+        print('data location not specified, use generic path: M:/')
+        return 'M:/'
+results_path = get_results_path()#'C:/Users/kinder.NLOQO/HESSENBOX-DA/OPA-Paper/analysis (python)/simulation/'
+plots_path = get_plots_path()#f'C:/Users/kinder.NLOQO/HESSENBOX-DA/OPA-Paper/analysis (python)/Ausgabe/'
+
+
+#%% parameters for OPA1
+#%%% calc radius of curvature: Red1 (horizontal)
+wavelength_nm = 1064
+ref_index = snlo.n_bbo(wavelength_nm, axis='o')
+fwhm_mm = 1.11
+focus_pos_mm = 6  # focus center of crystal
+result = snlo.focus(wavelength_nm, ref_index, fwhm_mm, focus_pos_mm)
+print(f'Red1 ({wavelength_nm}nm): radius of curvature: {result[2]:.3f} mm')
+
+#%%% calc radius of curvature: Red1 (vertical)
+wavelength_nm = 1064
+ref_index = snlo.n_bbo(wavelength_nm, axis='o')
+fwhm_mm = 1.48
+focus_pos_mm = 6  # focus center of crystal
+result = snlo.focus(wavelength_nm, ref_index, fwhm_mm, focus_pos_mm)
+print(f'Red1 ({wavelength_nm}nm): radius of curvature: {result[2]:.3f} mm')
+
+#%%% calc radius of curvature: Red2 (horizontal)
+wavelength_nm = 1417
+ref_index = snlo.n_bbo(wavelength_nm, axis='o')
+fwhm_mm = 0.153
+focus_pos_mm = 6  # focus center of crystal
+result = snlo.focus(wavelength_nm, ref_index, fwhm_mm, focus_pos_mm)
+print(f'Red2 ({wavelength_nm}nm): radius of curvature: {result[2]:.3f} mm')
+
+#%%% calc radius of curvature: Red2 (vertical)
+wavelength_nm = 1417
+ref_index = snlo.n_bbo(wavelength_nm, axis='o')
+fwhm_mm = 0.151
+focus_pos_mm = 6  # focus center of crystal
+result = snlo.focus(wavelength_nm, ref_index, fwhm_mm, focus_pos_mm)
+print(f'Red2 ({wavelength_nm}nm): radius of curvature: {result[2]:.3f} mm')
+
+#%%% calc radius of curvature: Blue
+wavelength_nm = 608
+ref_index = snlo.n_bbo(wavelength_nm, axis='o')
+fwhm_mm = 0.15
+focus_pos_mm = 6  # focus center of crystal
+result = snlo.focus(wavelength_nm, ref_index, fwhm_mm, focus_pos_mm)
+print(f'Blue ({wavelength_nm}nm): radius of curvature: {result[2]:.3f} mm')
+
+#%%% parameters 2D mix LP for OPA1
+red1_nm = 1064
+red2_nm = 1417
+blue_nm = 1/(1/red1_nm + 1/red2_nm)
+theta = 28.2  # from QMIX, SFG pm calculation currently broken...
+seedpower_w = 70e-6
+pump_j = 9.16e-3
+pumpduration_ns = 8
+
+values = [
+    [red1_nm, red2_nm, blue_nm],
+    [snlo.n_bbo(red1_nm, axis='o'), snlo.n(theta, red2_nm, crystal='BBO'), snlo.n(theta, blue_nm, crystal='BBO')],
+    [0, 0, 0],  # phases at input
+    [0, 0, 0],  # coating specs input
+    [0, 0, 0],  # coating specs output
+    [0, 0, 0],  # crystal loss 1/mm
+    [pump_j, seedpower_w, 0],
+    [8, 0, 0],
+    [[1.11, 1.48], [0.153, 0.151], 0.15],
+    [1, 1, 1],  # supergaussian coefficient
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 62.92, 63.64],  # walkoff form QMIX
+    [0, 0, 0],  # offset
+    [[1898700, 6001000], [388.73, 368.99], 1959.4],  # radius of curvature
+    [30, 32, 32],  # integration points
+    [12, 5.55, 5.92], # crystal length, grid, 0: auto: 2.4485x1.9588mm
+    [1.55],  # via QPM, LNB_M
+    [0],  # delta k
+    [0],  # distance to image
+    [20]  # time steps
+]
+
+#%%% init 2D-mix-LP parameters
+snlo.open_function('2D-mix-LP')
+for index in range(len(values)):
+    if len(values[index])==1:
+        pos = list(snlo.dict_2dmixlp.values())[index][0]
+        snlo.set_value(pos, values[index][0])
+    else:
+        for element in range(len(values[index])):
+            # if there are different values for horizontal and vertical
+            if type(values[index][element])==list:
+                value = str(values[index][element][0])+' '+str(values[index][element][1])
+                gui.moveTo(*coords(*list(snlo.dict_2dmixlp.values())[index][element]))
+                gui.doubleClick()
+                gui.press('delete')
+                gui.doubleClick()
+                gui.write(value)
+            else:
+                snlo.set_value(list(snlo.dict_2dmixlp.values())[index][element], values[index][element])
+
+
+#%% OPA1 pump energy dependence
+# start with 1e-5
+energies = np.arange(1e-6,0.7e-3,1e-5)
+print(energies)
+print(len(energies))
+
+
+#%%% run OPA1 pump energy dependence
+results = []
+for energy in energies:
+    set_value(snlo.dict_2dmixlp['Energy/power (J or W)'][2], energy)
+    snlo.moveto(*dict_2dmixlp['Accept'])
+    #snlo.moveto(780, 800)
+    gui.click()
+    dict = run_2dmixlp()
+    m.set_key(dict, 'pump pulse energy (J)', energy)
+    dict2 = get_spectr()
+    results.append(merge_dict([dict, dict2]))
+    snlo.moveto(*dict_2dmixlp['Change Inputs'])
+        #(600, 250)
+    gui.click()
+
+#%% save OPA1 pump dependence
+now = datetime.now()
+file_name = now.strftime("%Y-%m-%d") + '_' + 'OPA1_pump-dependence.pkl'
+pickle.dump(results, open(results_path + file_name, 'wb'))
+
+
+#%%% plot results
+file_name = '2022-06-09_OPA1_pump-dependence.pkl'
+#results = snlo.merge_dict(pickle.load(open('simulation/' + file_name, 'rb')))
+results = snlo.merge_dict(pickle.load(open(results_path + file_name, 'rb')))
+
+x = np.array(results['pump pulse energy (J)']).T[0]*1e6
+y = np.array(results['Output pulse energy (mJ)']).T[0][0]*1e3
+y2 = np.array(results['duration_idl_ns']).T[0]
+y3 = np.array(results['duration_sig_ns']).T[0]
+y4 = np.array(results['duration_pmp_ns']).T[0]
+
+fig, [ax1, ax2] = plt.subplots(1,2)
+m.plot_options(fig, ax1, image_width=15, aspect_ratio=2, xlabel='pump pulse energy ($\mathrm{\mu J}$)', ylabel='idler pulse energy ($\mathrm{\mu J}$)')
+ax1.plot(x, y, 'r-', label='idler pulse energy')
+ax1.legend(loc='upper left')
+ax12 = ax1.twinx()
+m.plot_options(fig, ax12, image_width=15, aspect_ratio=2)
+ax12.plot(x, 1e-6*y/(1e-9*y2), 'b--', label='single pass gain')  # divided by 1W...
+ax12.set_ylabel('single pass gain', color="blue")
+ax12.legend(loc='lower right')
+
+m.plot_options(fig, ax2, ticks=['auto', 2], image_width=15, aspect_ratio=2, xlabel='pump pulse energy ($\mathrm{\mu J}$)',
+               ylabel='idler pulse duration (ns)')
+ax2.plot(x, y2, 'r-', label='idler')
+ax2.plot(x, y3, 'g--', label='sig')
+ax2.plot(x[y4<100], y4[y4<100], 'b-', label='pump')
+ax2.set_ylim([0, max(y4[y4<100]+0.5)])
+ax2.legend(loc='lower center', ncol=3, frameon=True, handlelength=1.5, mode="expand")
+plt.tight_layout()
+#m.save_plot('OPA1_pump_sim')
+plt.show()
+
+#%% conversion efficiency OPA1
+x = np.array(results['pump pulse energy (J)']).T[0]*1e6
+y = np.array(results['Output pulse energy (mJ)']).T*1e3
+eff = ((y[0]+y[1])/x)[0]
+
+fig, ax1 = plt.subplots()
+m.plot_options(fig, ax1, image_width=15, xlabel='pump pulse energy ($\mathrm{\mu J}$)', ylabel='pulse energy ($\mathrm{\mu J}$)')
+ax1.plot(x, y[0][0], 'r-', label='idler pulse energy')
+ax1.plot(x, y[1][0], 'g-', label='signal pulse energy')
+ax1.plot(x, y[2][0], 'b-', label='residual pump pulse energy')
+ax1.legend(loc='upper left')
+
+ax12 = ax1.twinx()
+m.plot_options(fig, ax12, image_width=15, aspect_ratio=2)
+ax12.plot(x[eff<1], eff[eff<1], 'k--', label='conversion efficiency')
+ax12.set_ylabel('conversion efficiency')
+
+ax12.legend(loc='lower right')
+
+plt.tight_layout()
+#m.save_plot('OPA1_pump_sim')
+plt.show()
+
+#%% seed power dependence
+powers1 = np.arange(10e-6, 100e-6, 5e-6)
+powers2 = np.array([1e-19, 1e-18, 1e-17, 1e-16, 1e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-10, 1e-9, 1e-8, 1e-7, 1e-6])
+powers = m.flatten([powers1, powers2])
+
+#%%% run seed dependence
+set_value(dict_2dmixlp['Energy/power (J or W)'][0], pump_j)
+results = []
+for power in powers:
+    set_value(dict_2dmixlp['Energy/power (J or W)'][1], power)
+    gui.moveTo(*coords(780, 800))
+    gui.click()
+    dict1 = run_2dmixlp()
+    m.set_key(dict1, 'seed power (W)', power)
+    dict2 = get_spectr()
+    results.append(merge_dict([dict1, dict2]))
+    gui.moveTo(*coords(600, 250))
+    gui.click()
+
+# save OPA1 seed dependence
+#%%
+now = datetime.now()
+file_name = now.strftime("%Y-%m-%d") + '_' + 'upconversion_seed-dependence.pkl'
+pickle.dump(results, open('results/'+file_name, 'wb'))
+
+#%%% import and plot seed dependence
+file_name = '2022-08-02_upconversion_seed-dependence.pkl'
+results = merge_dict(pickle.load(open('results/'+file_name, 'rb')))
+
+x = np.array(results['seed power (W)']).T[0] * 1e6
+y1 = np.array(results['Output pulse energy (mJ)']).T[1][0] * 1e6
+y2 = np.array(results['Output pulse energy (mJ)']).T[2][0] * 1e6
+
+fig, ax1 = plt.subplots()
+m.plot_options(fig, ax1, xlabel='(cw) seed power ($\mu$W)', ylabel='idler pulse energy ($\mathrm{\mu J}$)')
+ax1.plot(x, y1, '-', label='seed')
+ax1.plot(x, y2, '-', label='SFM')
+ax1.legend()
+#m.plot_options(fig, ax2, xlabel='(cw) idler power (W)',ylabel='idler pulse duration (ns)')
+#ax2.plot(x, y2, '-')
+plt.tight_layout()
+# m.save_plot('OPA1_seed_sim')
+plt.show()
+
+#%% import and plot seed dependence
+file_name = '2022-08-02_upconversion_seed-dependence.pkl'
+results = merge_dict(pickle.load(open('results/'+file_name, 'rb')))
+h = scipy.constants.h
+c = scipy.constants.c
+seed_nm = 1417
+pump_nm = 1064
+sfm_nm = 1/(1/seed_nm + 1/pump_nm)
+x = seed_nm*1e-3*np.array(results['Output pulse energy (mJ)']).T[1][0]/(h*c)
+y1 = seed_nm*1e-3*np.array(results['Output pulse energy (mJ)']).T[1][0]/(h*c)
+y2 = sfm_nm*1e-3*np.array(results['Output pulse energy (mJ)']).T[2][0]/(h*c)
+
+fig, ax1 = plt.subplots()
+m.plot_options(fig, ax1, xlabel='(cw) seed photons', ylabel='idler pulse energy ($\mathrm{\mu J}$)')
+
+ax1.plot(x, y1, '.', label='seed')
+ax1.plot(x, y2, '^', label='SFM')
+#ax1.set_xscale('log')
+#ax1.set_yscale('log')
+ax1.loglog()
+ax1.xaxis.set_major_formatter(FormatStrFormatter('%.0f'))
+#ax1.set_xlim(left=1e10)
+#ax1.set_ylim(bottom=1e10)
+ax1.legend()
+#m.plot_options(fig, ax2, xlabel='(cw) idler power (W)',ylabel='idler pulse duration (ns)')
+#ax2.plot(x, y2, '-')
+plt.tight_layout()
+plt.ticklabel_format(axis='both',  useOffset=False, style='plain')
+# m.save_plot('OPA1_seed_sim')
+plt.show()
+
+#%%
+#%% plot conversion efficiency
+fig, ax1 = plt.subplots()
+m.plot_options(fig, ax1, xlabel='(cw) seed photons', ylabel='conversion efficiency (%)')
+h = scipy.constants.h
+c = scipy.constants.c
+seed_nm = 1417
+pump_nm = 1064
+sfm_nm = 1/(1/seed_nm + 1/pump_nm)
+x = 1e-3*np.array(results['Output pulse energy (mJ)']).T[1][0]
+y1 = 100*sfm_nm*np.array(results['Output pulse energy (mJ)']).T[2][0]/(seed_nm*np.array(results['Output pulse energy (mJ)']).T[1][0])
+
+ax1.plot(x[y1<50], y1[y1<50], '-')
+#ax1.set_xscale('log')
+#ax1.set_yscale('log')
+#ax1.set_xlim(left=1e10)
+#ax1.set_ylim(bottom=1e10)
+#ax1.legend()
+#m.plot_options(fig, ax2, xlabel='(cw) idler power (W)',ylabel='idler pulse duration (ns)')
+#ax2.plot(x, y2, '-')
+plt.tight_layout()
+# m.save_plot('OPA1_seed_sim')
+plt.show()
+
+#%% get result
+pump_energy =  0.00916
+seed_power = 7e-05
+set_value(dict_2dmixlp['Energy/power (J or W)'][2], seed_power)
+set_value(dict_2dmixlp['Energy/power (J or W)'][0], pump_energy)
+moveto(dict_2dmixlp['Accept'])
+gui.click()
+dict1 = run_2dmixlp()
+m.set_key(dict1, 'seed power (W)', seed_power)
+m.set_key(dict1, 'pump pulse energy (J)', pump_energy)
+dict2 = get_spectr()
+results = merge_dict([dict1, dict2])
+#%% gather OPA1 values
+duration_idl_ns = results['duration_idl_ns'][0]
+fwhm_idl_MHz = results['fwhm_idl_MHz'][0]
+energy_sfm_mJ = results['Output pulse energy (mJ)'][0][2]
+y = np.array(results['Output pulse energy (mJ)'])*1e3
+input_energy_mJ = results['seed power (W)'][0]*8e-9*1e3
+efficiency = (results['Output pulse energy (mJ)'][0][2]/input_energy_mJ)
+#gain = (1e-3*energy_idl_mJ)/(1e-9 * duration_idl_ns)/results['seed power (W)'][0]
+
+print(f'idler pulse energy: {energy_idl_mJ * 1e3:.3f}uJ')
+print(f'conversion efficiency: {efficiency*1e2:.3f}%')
+print(f'single pass gain: {gain:.1f}')
+print(f'idler duration:{duration_idl_ns:.3f}ns')
+print(f'idler bandwidth: {fwhm_idl_MHz:.3f}MHz')