'''
Author: Ugo Lo Cicero
Email: ugo.locicero@inaf.it

v. 03.2 31/01/2023
- estratta la parte di calcolo dalla funzione transmittanceM, trasferita nella funzione _transmittanceCalc
- ottimizzate per numba le funzioni Forouhi19, _updateM e _transmittanceCalc (ridotto il tempo medio di esecuzione di "residuals()" da 0.38 s a 0.15 s per il file TF110-2466-PI415_UV-MIR.dat).
- scambiato l'ordine delle dimensioni della matrice M per ottimizzazione numba: la prima dimensione adesso è la lunghezza d'onda
- ottimizzato update_plots per evitare di calcolare più volte l'indice di rifrazione
- aggiunto il parametro riPI a residual e transmittanceM per poter passare lindice di rifrazione del PI se già calcolato

v. 03.1 17/01/2023
- implementati i vincoli W < 4*np.sqrt(2)*(Ep - Eg) e Eh>Eg per i termini E>Eg del poliimide.
- aggiunti bottoni Fix all e Unfix all: IL COMPORTAMENTO NON SEMBRA AFFIDADBILE!!
- introdotta variabile ε per valori prossimi allo zero
- compilazione con libreria numba per la funzione _updateM

v. 03 04/01/2023
- implementato salvataggio e caricamento dei parametri in .json (mantenuta la modalità pickle).
	Il .json è leggibile. Il salvataggio in .json è dei parametri correnti (params), mentre in .pkl è dei parametri di output del fit (oggetto out). Il .pkl contiene anche altre informazioni dell'oggetto out.
- implementate le nuove equazioni:
	B = -np.sqrt(4*Ep**2-8*Ep*X+4*X**2+W**2)+4*Ep-2*X
	C = Ep**2-2*Ep*X+B*X
	A = P*(2*Ep-B)
	La X rappresenta una soglia di energia, è =Eg nel caso degli amorfi E>Eg, è =0 per amorfi con E<Eg, oppure una soglia di interbanda per i metalli
- implementazione matriciale dei calcoli
- implementato il controllo Wlim < 4*np.sqrt(2)*(Ep - Eg) con eventuale stampa in console
- aggiunto bottone per aggiungere termini per amorfi E<Eg
- quando si caricano parametri da file la GUI si aggiorna mostrando i nuovi parametri, anche se il numero di parametri è cambiato
- cliccando su un'etichetta viene mostrata la linea di Ep corrispondente, cliccando nuovamente si cancella; la linea è fissa, non si sposta se varia il parametro; funziona solo con la visualizzazione in Energia

v. 02 03/01/2023
- modifica di get_parameters_from_fit per bloccare i segnali e velocizzare il caricamento
- aggiunto bottone di commutazione tra visualizzazione in λ o in E
- la generazione dei plot è racchiusa in una funzione piuttosto che essere globale
- emissione suono alla fine del fit (Edo)
- cambio formato visualizzazione parametri GUI :3e (o :2e) (Edo)

v. 01.1 20/12/2021
- aggiunto bottone Toggle Fix per commutare tutti i Fix
- aggiunti Al e Al2O3 da dati tabulati combinati

v. 01 17/12/2021

TODO
- implementare i vincoli (Ep>X, W<Wlim(Ep-X))
- utilizzare il metodo delle etichette per potere gestire più materiali e la altre tipologie di termini
	es. PI_L_E: contributo per il PI per E<Eg, PI_H_E: contributo per il PI per E>Eg,  Al_inter_E: contributo per l'Al interbanda, Al_intra_E: contributo per l'Al intrabanda, Al_inel_E: contributo per l'Al inelastico. Usare funzione split('_') per separare le parti dell'etichetta.
- aggiunzione termine: aggiungere due menu a tendina, una per selezionare il materiale e una per selezionare la tipologia di termine aggiunto. Queste informazioni vengono passate a add_fit_term() che dovrà generare il nuovo termine di conseguenza, anche IMPLEMENTANDO I VINCOLI

- dividere il bottone save in due bottoni: "Save init." (salva params), "Save out" (salva out.params in json o out in pkl)
- verificare che il salvataggio in .json e .pkl non tronchi i valori
- modificare visualizzazione linee verticali Ep spostando il disegno in update_plots e gestendo la visualizzazione per lunghezza d'onda
x visualizzazione linee ed etichette per identificare i contributi nei plot (disattivabili - ad esempio si può commutare cliccando sull'etichetta della Ep)
- gestione dinamica dei parametri per potere aggiungere e rimuovere parametri direttamente da GUI
x gestione dei parametri per consentire il caricamento di file parametri salvati anche in caso di cambiamento del numero di contributi
x modifica di get_parameters_from_fit per bloccare i segnali e velocizzare il caricamento
x switch da GUI tra visualizzazione in energia e in lunghezza d'onda
x posizionamento e dimensionamento automatico finestre all'avvio

'''


import numpy as np
import lmfit as lm
import math, copy, time, pickle, sys, traceback, winsound
import copy
import matplotlib as mpl
from matplotlib import pyplot as plt
from matplotlib.widgets import Slider, Button, TextBox
from PyQt5 import QtCore, QtWidgets
from PyQt5.QtCore import Qt
from PyQt5.QtWidgets import (QApplication, QMainWindow, QDesktopWidget, QCheckBox, QGridLayout, QGroupBox, QFrame, QMenu, QPushButton, QRadioButton, QVBoxLayout, QHBoxLayout, QWidget, QSlider, QLabel, QLineEdit, QScrollArea, QFileDialog)
from numba import jit, njit	# versione numba 0.56.4
#from numba.typed import Dict


# PI 415
#filename = r"G:\Drive condivisi\ATHENA\EXPERIMENTAL_STUFF\_ DATI PER CAMPIONE\Polyimide 415 nm\Dati combinati\TF110-2466_PI415nm_UV-MIR.dat"
#filename = r"G:\Drive condivisi\ATHENA\EXPERIMENTAL_STUFF\UV-VIS-NIR TRASMISSION\Lambda1050\20211217\tf110_2466_PI415.asc"
filename = r"G:\Drive condivisi\ATHENA\Software\UV-Vis-IR-fit\misure\TF110-2466-PI415_UV-MIR.dat"

# PI 150
#filename = r"G:\Drive condivisi\ATHENA\EXPERIMENTAL_STUFF\_ DATI PER CAMPIONE\Polyimide 150 nm\Dati combinati\TF110-2482_PI_150nm_UV_MIR.dat"

# PI 45
#filename = r"G:\Drive condivisi\ATHENA\EXPERIMENTAL_STUFF\_ DATI PER CAMPIONE\Polyimide 45 nm\TF110-2476\Dati combinati\TF110-2476_PI45nm_UV-MIR.dat"
#filename = r"G:\Drive condivisi\ATHENA\EXPERIMENTAL_STUFF\UV-VIS-NIR TRASMISSION\Lambda1050\20211217\TF110_2476_PI45.asc"

# PI 45/ Al 30
#filename = "G:\Drive condivisi\ATHENA\EXPERIMENTAL_STUFF\_ DATI PER CAMPIONE\Polyimide 45 nm   Aluminum 30 nm\TF110-2462\Dati combinati\TF110-2462_PI45nm_Al30nm_UV-MIR.dat"
#filename = r"misure\tran_UV-FIR_45-30_TF110-2462_in_house.dat"

# PI 150 / Al 30
#filename = "G:\Drive condivisi\ATHENA\EXPERIMENTAL_STUFF\_ DATI PER CAMPIONE\Polyimide 150 nm   Aluminum 30 nm\TF110-2460\Dati combinati\TF110-2460_PI150nm_Al30nm_UV-MIR.dat"

xtype = 'nm'	#'cm^-1'
ytype = 'T'		#'T' 'A'
skiprows = 0
dpi = 100
wav_min = 190			#minima lunghezza d'onda del fit
wav_max = 50000			#massima lunghezza d'onda del fit
plot_ref_index = False
plot_transmittance = True
plot_residuals = True
res_view_mode = 'percent' 		# Può essere 'diff' o 'percent'
viewE = True					# True to view the plots in energy instead of wavelength
xscale = 'log'				# Scala lineare o logaritmica per gli assi x dei grafici
yscale = 'linear'
exclude_from_fit = [] #[(5893,7213)]	# regioni escluse dal fit (in nm)
show_excluded_areas = True
EW = 1239.8		# conversione eV-nm
ε = 1e-12		# valore piccolo per evitare divisioni per zero

# FIT PARAMETERS
params = lm.Parameters()
excl_pars = ['p', 'q', 'CB0', 'CB1', 'CB2', 'DEh', 'DWh']			# parametri esclusi dal GUI
# Forouhi19
# Picchi di assorbimento in eV
#ppEl = []
#ppEl = [0.425, 4.6, ]
#ppEl = [0.053, 0.0679, 0.0899, 0.103, 0.135, 0.153, 0.17, 0.23, 0.425, 4.6, ]
#ppEl = [0.053, 0.0679, 0.0899, 0.103, 0.135, 0.153, 0.17, 0.21, 0.23, 0.425, 4.6, ]
#ppW = [EW/e for e in ppE]
#ppEh = [4]#, 5, 6.1]

# Picchi di assorbimento E: energia picco (eV), P: massimo k, W: FWHM (eV)
# Per ogni parametro è indicata una tupla con (valore, minimo, massimo)
# Picchi con E<Eg:
ppEl =[
#{'E':(0.05,0.001,6.5)	, 'P':(1,0,100)	, 'W':(0.01,0,1)	, 'lbl':'00'},
# {'E':(0.05301,0.04,0.08)	, 'P':(0.0735,0,0.21)	, 'W':(0.006509,ε,0.03)	, 'lbl':'00'},
# {'E':(0.06664,0.05,0.1)	, 'P':(0.09499,0,0.27)	, 'W':(0.01265,ε,0.03)	, 'lbl':'01'},
# {'E':(0.08852,0.07,0.13)	, 'P':(0.09846,0,0.27)	, 'W':(0.01235,ε,0.03)	, 'lbl':'02'},
# {'E':(0.1044,0.08,0.15)	, 'P':(0.09854,0,0.30)	, 'W':(0.01316,ε,0.03)	, 'lbl':'03'},
# {'E':(0.1353,0.1,0.2)	, 'P':(0.111,0,0.38)	, 'W':(0.01231,ε,0.03)	, 'lbl':'04'},
# {'E':(0.1517,0.11,0.23)	, 'P':(0.05758,0,0.19)	, 'W':(0.008374,ε,0.03)	, 'lbl':'05'},
# {'E':(0.168,0.13,0.26)	, 'P':(0.3797,0,1.4)	, 'W':(0.008688,ε,0.03)	, 'lbl':'06'},
# {'E':(0.2092,0.16,0.32)	, 'P':(0.1011,0,0.30)	, 'W':(0.03291,ε,0.06)	, 'lbl':'07'},
# {'E':(0.3256,0.17,0.50)	, 'P':(0.03555,0,0.12)	, 'W':(0.1161,ε,0.25)	, 'lbl':'08'},
# {'E':(0.4473,0.32,2.00)	, 'P':(0.03447,0,0.12)	, 'W':(0.143,ε,1.0)	, 'lbl':'09'},

{'E':(3.521,2.5,6.9)	, 'P':(11.25,0,50)	, 'W':(2.695e-2,ε,0.1)	, 'lbl':'gap'},

{'E':(0.05123704,0.04,0.08)	, 'P':(0.07607364,0,1)	, 'W':(0.00957728,ε,5e-2)	, 'lbl':'00'},
{'E':(0.06088548,0.04,0.08)	, 'P':(0.03264536,0,1)	, 'W':(0.01299649,ε,5e-2)	, 'lbl':'01'},
{'E':(0.0625,0.04,0.08)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'02'},
{'E':(0.0655,0.04,0.08)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'03'},
{'E':(0.0681,0.04,0.08)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'04'},
{'E':(0.0705,0.04,0.08)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'05'},
{'E':(0.0728,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'06'},
{'E':(0.0760,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'07'},
{'E':(0.0788,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'08'},
{'E':(0.0833,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'09'},
{'E':(0.0862,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'10'},
{'E':(0.0913,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'11'},
{'E':(0.0932,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,5e-2)	, 'lbl':'12'},
{'E':(0.0946,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'13'},
{'E':(0.0986,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'14'},
{'E':(0.1030,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'15'},
{'E':(0.1070,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'16'},
{'E':(0.1100,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'17'},
{'E':(0.1160,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,5e-2)	, 'lbl':'18'},
{'E':(0.1210,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,5e-2)	, 'lbl':'19'},
{'E':(0.1265,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'20'},
{'E':(0.1337,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'21'},
{'E':(0.1360,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'22'},
{'E':(0.1390,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'23'},

{'E':(0.1456,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'24'},
{'E':(0.1511,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'25'},
{'E':(0.1568,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'26'},
{'E':(0.1679,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'27'},
{'E':(0.1760,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'28'},
{'E':(0.1786,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,5e-2)	, 'lbl':'29'},
{'E':(0.1823,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'30'},
{'E':(0.1877,0.04,0.20)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'31'},
{'E':(0.2005,0.04,0.40)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'32'},
{'E':(0.2124,0.04,0.40)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'33'},
{'E':(0.2197,0.04,0.40)	, 'P':(1e-1,0,1)	, 'W':(1e-3,ε,1e-2)	, 'lbl':'34'},

]
# Picchi con E>Eg
ppEh =[
#{'E':(3.08,3,4.5)		, 'P':(0.4664,0,1)		, 'W':(0.11,ε,1)	, 'lbl':'00'},
#{'E':(4.275,3,4.5)		, 'P':(0.4664,0,1)		, 'W':(3.301,ε,10)	, 'lbl':'00'},
##{'E':(5.121,4.5,6)		, 'P':(0.3463,0,3)		, 'W':(32.56,ε,60)	, 'lbl':'01'},
#{'E':(5.121,4.5,6)		, 'P':(0.3463,0,3)		, 'W':(11.6,ε,60)	, 'lbl':'01'},
#{'E':(6.366,6.1,6.7)	, 'P':(0.2077,0,0.4)	, 'W':(1.243,ε,10)	, 'lbl':'02'},
#{'E':(3.445,3,4.5)		, 'P':(0.03357,0,0.17)		, 'W':(0.2474,ε,2.1)	, 'lbl':'03'}
]

# ottimizzazione per accelerare il disegno
mpl.rcParams['path.simplify'] = True
mpl.rcParams['path.simplify_threshold'] = 1
mpl.rcParams["legend.loc"] = 'upper right'		# fissa la posizione della legenda per non doverla ricalcolare sempre

p = len(ppEl)
q = len(ppEh)

params.add('n_inf'			, value=1.689	, min=1,	max=2	, vary=False) #1.64278914
params.add('Eg'				, value=3.06	, min=0,	max=4	, vary=False)
params.add('thick_PI'		, value=432.8	, min=390,	max=450	, vary=False)
params.add('thick_Al'		, value=0		, min=0,	max=40	, vary=False)
params.add('thick_Al2O3_1'	, value=0		, min=0,	max=7	, vary=False)
params.add('thick_Al2O3_2'	, value=0		, min=0,	max=7	, vary=False)

for j in range(0,p):
	# params.add('El'+str(j), value=ppEl[j]['E'][0], min=ppEl[j]['E'][1], max=ppEl[j]['E'][2], vary=False)
	# params.add('Pl'+str(j), value=ppEl[j]['P'][0], min=ppEl[j]['P'][1], max=ppEl[j]['P'][2], vary=False)
	# params.add('Wl'+str(j), value=ppEl[j]['W'][0], min=ppEl[j]['W'][1], max=ppEl[j]['W'][2], vary=False)
	#params.add('El'+ppEl[j]['lbl'], value=ppEl[j]['E'][0], min=ppEl[j]['E'][1], max=ppEl[j]['E'][2], vary=False)
	params.add('El'+ppEl[j]['lbl'], value=ppEl[j]['E'][0], min=ppEl[j]['E'][0]*(1-0.005), max=ppEl[j]['E'][0]*(1+0.005), vary=False)
	params.add('Pl'+ppEl[j]['lbl'], value=ppEl[j]['P'][0], min=ppEl[j]['P'][1], max=ppEl[j]['P'][2], vary=False)
	params.add('Wl'+ppEl[j]['lbl'], value=ppEl[j]['W'][0], min=ppEl[j]['W'][1], max=ppEl[j]['W'][2], vary=False)

for j in range(0,q):
	params.add('DEh'+str(j), value=ppEh[j]['E'][0]-params['Eg'].value, min=ε)	#vincolo Eh>Eg => Eh=DEh+Eg, DEh>0 
	params.add('Eh'+str(j), value=ppEh[j]['E'][0], min=ppEh[j]['E'][1], max=ppEh[j]['E'][2], vary=False, expr='DEh'+str(j)+'+Eg')
	params.add('Ph'+str(j), value=ppEh[j]['P'][0], min=ppEh[j]['P'][1], max=ppEh[j]['P'][2], vary=False)
	params.add('DWh'+str(j), 5.65685*(ppEh[j]['E'][0] - params['Eg'].value) - ppEh[j]['W'][0], min=ε)	#vincolo W < 4*np.sqrt(2)*(Ep - Eg)
	params.add('Wh'+str(j), value=ppEh[j]['W'][0], min=ppEh[j]['W'][1], max=ppEh[j]['W'][2], vary=False, expr='5.65685*(Eh'+str(j)+' - Eg) - DWh'+str(j))


# params.add('A0'		, value=0.059		, min=0	, max=4		, vary=True)
# params.add('B0'		, value=8.391		, min=0	, max=20	, vary=True)
# params.add('C0'		, value=17.811		, min=0	, max=25	, vary=True)
# params.add('CB0'	, expr='C0-B0**2/4'	, min=0				, vary=False)
# params.add('A1'		, value=1.457		, min=0	, max=4		, vary=True)
# params.add('B1'		, value=7.763		, min=0	, max=20	, vary=True)
# params.add('C1'		, value=19.814 		, min=0	, max=50	, vary=True)



# E<Eg: 
# Emax = sqrt(Cj)
# kmax = A/(2*sqrt(C)-B)

# E>Eg: 
# Emax = Eg+sqrt(Eg**2-Bj*Eg+Cj) (per la stima del valore iniz. si può fare Emax=sqrt(Cj))

# E2FWHM: (-B+4*sqrt(C)+sqrt((B-4*sqrt(C))**2-4*C))/2
# E1FWHM: (-B+4*sqrt(C)-sqrt((B-4*sqrt(C))**2-4*C))/2
# FWHM(eV): sqrt((B-4*sqrt(C))**2-4*C)
# B = 4*sqrt(C)-sqrt(FWHM**2+4*C)
# W1FWHM: 1239.8/(-B+4*sqrt(C)+sqrt((B-4*sqrt(C))**2-4*C))*2
# W2FWHM: 1239.8/(-B+4*sqrt(C)-sqrt((B-4*sqrt(C))**2-4*C))*2
# FWHM (nm): 2*1239.8*(1/(-B+4*sqrt(C)-sqrt((B-4*sqrt(C))**2-4*C))-1/(-B+4*sqrt(C)+sqrt((B-4*sqrt(C))**2-4*C)))


'''
PI
Nella zona a k=0 l'indice di rifrazione è legato al minimo della trasmittanza:
Tmin = (2*n/(1+n**2))**2
n = 1/np.sqrt(Tmin)+np.sqrt((1-Tmin)/T)

La distanza fa due massimi è legata allo spessore d:
Δυ = 1/(2*n*d)
d = 1/(2*n*Δυ)

Per trovare i picchi (in nm):
#wav_sel = (wav>503) & (wav<555)
wav_sel = (wav>730) & (wav<800)
#wav_sel = (wav>1400) & (wav<1600)
wav_peak = wav[wav_sel]
peak = data[wav_sel]
wav_peak[np.where(peak==np.max(peak))]

Campione nominale 415:
picchi massimi: λ=526; 760; 1488.0952381
Tmin = 0.754 => n = 1.44
Lunghezza d'onda massimi: 760 nm (1.316e-3 nm^-1)  1485 nm (6.734e-4 nm^-1)
Δυ = 6.426e-4
Assumendo d=415 nm => n=1.875
Usando invece l'n calcolato prima viene: d=540.3

I picchi legati all'interferenza hanno una Δλ (e Δν) non costante e si riducono andando verso Eg. Ciò denota la presenza di un termine di assorbimento, valido anche per E<Eg, ma con picco su energia maggiore. Si ottiene un andamento che approssima (da migliorare) usando: n_inf=1.614, PI_thickness=437.7 nm, e un termine con E=3.847 eV, P=16, W=3.576e-02 eV
'''

def Forouhi19w(wave, parvals, sel=''):
	values_list = [float(i) for i in parvals.values()]
	return Forouhi19(wav,(list(parvals.keys()),values_list))

@njit
def Forouhi19(wave, parlist, sel=''):
	'''
	wave: array lunghezze d'onda
	parlist: tupla (k,v), con k lista dei nomi dei parametri e v lista dei corrispondenti valori
	(uso una tupla di liste e costruisco il dizionario dentro perché "numba" non supporta dizionari come parametri di funzione)
	sel: stringa di selezione del materiale (es. Al_ o PI_)
	
	Eg: energy gap
	n_inf: n per Energia infinita. Se ci sono abbastanza dati tende a zero
	Aj, Bj, Cj: parametri per i vari contributi
	Posizione dei picchi:
	E<Eg: Emax = sqrt(Cj)
	E>Eg: Emax = Eg+sqrt(Eg**2-Bj*Eg+Cj) (per la stima del valore iniz. si può fare Emax=sqrt(Cj))

	picchi di assorbimento in nm per E<Eg:
	pp = np.array([440, 1200, 5800, 7500, 9170, 12e3, 14e3, 18e3, 23e3])
	'''
	
	parvals = dict()
	for key, val in zip(parlist[0],parlist[1]):
		parvals[key] = val
	#parvals = pars.valuesdict()
	#p = parvals['p']
	#q = parvals['q']
	#p = 0	# n. termini E<Eg (Fourouhi2019), viene calcolato in base alle etichette dei parametri
	#q = 0	# n. termini E>Eg (Fourouhi2019), viene calcolato in base alle etichette dei parametri
	n_inf = parvals['n_inf']
	Eg = parvals['Eg']
	
	# Classifica e conta i termini in base alla stringa:
	# sel è la selezione per materiale dei prametri da utilizzare, li prende tutti se non specificata
	# El, Pl, Wl -> low energy (contributo di tipo 'p')
	# Eh, Ph, Wh -> high energy (contributo di tipo 'q')
	Ep_arr = []
	X_arr = []		# vettore contenente le Eg. Viene popolato in base alle etichette
	P_arr = []
	W_arr = []
	keys = []
	for key,val in parvals.items():
		if sel in key[:len(sel)]:	# la stringa di selezione del materiale deve essere all'inizio
			key = key[len(sel):]
			if ('El' in key[:2]): # or ('Pl' in key) or ('Wl' in key):
				#p += 1
				Ep_arr.append(val)
				X_arr.append(0.0)	# i contributi 'l' hanno Eg=0
			elif ('Eh' in key[:2]): # or ('Ph' in key) or ('Wh' in key):
				#q += 1
				Ep_arr.append(val)
				print(type(Eg))
				X_arr.append(Eg)	# i contributi 'h' hanno Eg=Eg
			elif ('P' in key[:1]):
				P_arr.append(val)
			elif ('W' in key[:1]):
				W_arr.append(val)
				keys.append(key)
	
	# Converte le liste in array numpy per fare i calcoli
	Ep = np.array(Ep_arr)
	X = np.array(X_arr)
	P = np.array(P_arr)
	W = np.array(W_arr)
	
	#n_contribs = p+q
	#A, B, C = np.empty(n_contribs), np.empty(n_contribs),np.empty(n_contribs)
	#for j in range(n_contribs):
	
	# Effettuo un calcolo matriciale. B, C, A sono vettori che contengono i parametri per tutti i contributi
	B = -np.sqrt(4*Ep**2-8*Ep*X+4*X**2+W**2)+4*Ep-2*X
	C = Ep**2-2*Ep*X+B*X
	A = P*(2*Ep-B)
	'''
	for j in range(p):
		# A[j] = (parvals['A'+str(j)])
		# B[j] = (parvals['B'+str(j)])
		# C[j] = (parvals['C'+str(j)])
		Ej = (parvals['El'+str(j)])
		Wj = (parvals['Wl'+str(j)])
		Pj = (parvals['Pl'+str(j)])
		C[j] = Ej**2
		B[j] = 4*math.sqrt(C[j])-math.sqrt(Wj**2+4*C[j])
		A[j] = Pj*(2*math.sqrt(C[j])-B[j])
	for i in range(0,q):
		j = p+i
		# A[j] = (parvals['A'+str(j)])
		# B[j] = (parvals['B'+str(j)])
		# C[j] = (parvals['C'+str(j)])
		Ej = (parvals['Eh'+str(i)])
		Wj = (parvals['Wh'+str(i)])
		Pj = (parvals['Ph'+str(i)])
		#B[j] = 2*Ej
		#C[j] = Wj*B[j]**2/4
		#A[j] = Pj
		C[j] = Ej**2
		
		#C = B**2/4 + W/10 +1e-5
		B[j] = 2*math.sqrt(C[j]-Wj/10-1e-5)
		#B[j] = 4*math.sqrt(C[j])-math.sqrt(Wj**2+4*C[j])
		#B[j] = math.sqrt(4*C[j]/Wj)
		Emax = Eg + math.sqrt(Eg**2-B[j]*Eg+C[j])
		A[j] = Pj * (Emax**2-B[j]*Emax+C[j])/math.sqrt(Eg**2-B[j]*Eg+C[j])
		# Nel caso Eg=0 diventa Emax=sqrt(C); A = P*(C-B*sqrt(C)+C)/sqrt(C) = P*(sqrt(C)-B+sqrt(C)) = P*(2*sqrt(C)-B), identica al caso per E<Eg
	'''
	#Q = (C-B**2/4)**0.5
	Q1 = (C-B**2/4)
	bad = np.where(Q1<=0)[0]
	if bad.size>0:
		#print(bad)
		for ind in bad:
			print(f"Bad parameter values for: {keys[int(ind)]}")
			### il seguente codice mostra il limite di W e il valore attuale e fa un beep
			### non è supportato da numba, deve essere commentato se si attiva numba
			# Wbad = W[ind]
			# Xbad = X[ind]
			# Epbad = Ep[ind]
			# Wlim = 4*np.sqrt(2)*(Epbad - Xbad)
			# print(f"W limit: {Wlim}, W = {Wbad}")
			# winsound.MessageBeep()
	Q = Q1**0.5
	D = A/Q * (X-B/2)
	D1 = -A*B/(2*Q)
	F = A/Q * (C- X*B/2)
	F1 = A*C/Q
	
	E = np.expand_dims(EW/np.asarray(wave), axis=1)
	
	
	denom = E**2 - B*E + C
	k_arr = np.where(E>X, A * (E-X)/denom, 0)
	n_arr = np.where(E>X, (D * E + F)/ denom, 0)

	k = np.sum(k_arr,axis=1)
	n = np.sum(n_arr,axis=1)
	#print(f"k.shape: {k.shape}")

	'''
	E = EW/np.asarray(wave)
	k = np.zeros_like(E)
	n = np.zeros_like(E)
	# E < Eg
	#low = E<=Eg
	low = E>0
	for j in range(p):
		denom = E[low]**2 - B[j]*E[low] + C[j]
		k[low] += (A[j] * E[low])/ denom
		n[low] += (D1[j] * E[low] + F1[j])/ denom
		
	# E > Eg
	high = E>Eg
	for j in range(p,p+q):
		denom = E[high]**2 - B[j]*E[high] + C[j]
		k[high] += (A[j] * (E[high]-Eg))/ denom
		n[high] += (D[j] * E[high] + F[j])/ denom
	'''
	k[k<0] = 0		# imposta a zero i k negativi
	n += n_inf
	#print(f"N: {n - 1j*k}")
	return(n - 1j*k)


def refInd_Al(wave):
	w,n,k = np.genfromtxt('ref_index Al.csv',delimiter=';').T
	n_w = np.interp(wave, w, n)
	k_w = np.interp(wave, w, k)
	return(n_w - 1j*k_w)


def refInd_Al2O3(wave):
	w,n,k = np.genfromtxt('ref_index Al2O3.csv',delimiter=';').T
	n_w = np.interp(wave, w, n)
	k_w = np.interp(wave, w, k)
	return(n_w - 1j*k_w)


@njit
def _updateM(M_multi=np.empty((0,0,0),dtype='complex128'),wavelengths=None,refindex=1,thickness=0):
	n_wave = wavelengths.size
	delta = ((2*np.pi)/wavelengths)*refindex*thickness
	M_layer	= np.zeros((n_wave,2,2),dtype=np.complex128)
	delta_imag_max = 700
	delta_ovfl = np.abs(delta.imag)>delta_imag_max
	delta[delta_ovfl] = delta[delta_ovfl].real + delta_imag_max
	M_layer[:,0,0] = np.cos(delta)
	M_layer[:,0,1] = 1j*np.sin(delta)/refindex
	M_layer[:,1,0] = 1j*np.sin(delta)*refindex
	M_layer[:,1,1] = np.cos(delta)

	if M_multi.shape == (0,0,0):
		M_multi = M_layer
	for ii in range(0,n_wave):
		M_multi[ii,:,:] = np.dot(M_multi[ii,:,:],M_layer[ii,:,:])
		#M_multi[:,:,ii] = np.dot(np.dot(M_Al1[:,:,ii],M_Poly[:,:,ii]),M_Al2[:,:,ii])
	return(M_multi)


@njit
def _transmittanceCalc(M_multi,inc_index,exit_index):
	inc_n = np.real(inc_index)
	exit_n, exit_k = np.real(exit_index), -np.imag(exit_index)
	X = M_multi[:,0,0]
	Y = - 1j*inc_n*exit_n*M_multi[:,0,1]
	W = exit_n*M_multi[:,1,1]
	V = - 1j*M_multi[:,1,0] - exit_k*M_multi[:,1,1]
	
	t_ampl = 2*inc_index/ ((X+W) + 1j*(Y+V))
	r_ampl = ((X-W) + 1j*(Y-V))/((X+W) + 1j*(Y+V))
	
	# Calcolo dell'intensita` della trasmissivita' e della riflettivita'
	transmission = (exit_n/inc_n)* (np.real(t_ampl)**2 + np.imag(t_ampl)**2)
	reflection	 = (np.real(r_ampl)**2 + np.imag(r_ampl)**2)
	return(transmission,reflection)


def transmittanceM(wav,inc_index=1,exit_index=1,parlist=([],[]),riPI=None):
	#nc = polyModel2(wav,params)
	M_multi = np.empty((0,0,0),dtype='complex128')
	
	#parvals = dict(params.valuesdict())
	parvals = dict()
	for key, val in zip(parlist[0],parlist[1]):
		parvals[key] = val
	
	# Al2O3 layer #1
	nc = refInd_Al2O3(wav)
	thickness = parvals['thick_Al2O3_1']
	M_multi=_updateM(M_multi,wav,nc,thickness)
	# Al layer
	nc = refInd_Al(wav)
	thickness = parvals['thick_Al']
	M_multi=_updateM(M_multi,wav,nc,thickness)
	# Al2O3 layer #2
	nc = refInd_Al2O3(wav)
	thickness = parvals['thick_Al2O3_2']
	M_multi=_updateM(M_multi,wav,nc,thickness)
	# PI layer
	if riPI is None:
		riPI = Forouhi19(wav,parlist)
	thickness = parvals['thick_PI']
	M_multi=_updateM(M_multi,wav,riPI,thickness)
	
	# Calcolo della trasmissivita' e della riflettivita'
	transmission,reflection = _transmittanceCalc(M_multi,inc_index,exit_index)
	return(transmission,reflection)


def residual(pars, x, data=None, res_mode='rel', eps=None, model=None, riPI=None):
	'''
	eps: vettore incertezza dati
	'''
	print('.',end='',flush=True)
	parvals = pars.valuesdict()
	#n = parvals['n']
	#k = parvals['k']
	#th = parvals['thick']
	#model = transmittance_layer(x,th,n,k,)
	#model = transmittance_layer(x,th,poly_params=pars)
	values_list = [float(i) for i in parvals.values()]
	if model == None:
		model = transmittanceM(x,parlist=(list(parvals.keys()),values_list),riPI=riPI)[0]
	
	#model = transmittanceM(x,params=pars)[0]
	if data is None:
		return model
	if res_mode == 'diff':
		res = model - data
	elif res_mode == 'log_diff':
		res = np.log(model)-np.log(data)
	elif res_mode == 'rel':
		res = (model - data)/data
	if eps is None:
		return res
	
	return res / eps



def fit():
	global out
	print('Starting fit')
	fit_mask = (wav>=wav_min) & (wav<=wav_max)
	for exc in exclude_from_fit:
		fit_mask = fit_mask & ((wav<=exc[0]) | (wav>=exc[1]))	# esclude le aree in exclude_from_fit
	#out = lm.minimize(residual, params, args=(wav[fit_mask], data[fit_mask]),nan_policy='propagate')
	out = lm.minimize(residual, params, args=(wav[fit_mask], data[fit_mask]))
	print(lm.fit_report(out))
	out.params.pretty_print()
	winsound.MessageBeep()
	#winsound.Beep(440, 800)		# emette un suono appena finito il fit
	mask = plots['mask']
	trans_fit = residual(out.params,wav[mask])
	plots['lfit'].set_ydata(trans_fit)
	if res_view_mode == 'diff':
			res = trans_fit-data[mask]
			plots['ax_res'].set_ylim((-1,1))
	elif res_view_mode == 'percent':
		res = (trans_fit-data[mask])/data[mask]*100
		plots['ax_res'].set_ylim((-100,100))
	plots['lres'].set_ydata(res)
	plots['ftran'].canvas.draw()
	plots['ftran'].canvas.flush_events()
	print('Fit completed')


def plot_excluded_areas(ax):
	# Plot ranges excluded from fit
	if show_excluded_areas == False: return
	lims = ax.get_xlim()
	if not getattr(ax,'aspan_list',None) is None:
		for asp in ax.aspan_list:
			asp.remove()
	else:
		for ex in exclude_from_fit:
			if viewE:
				ax.axvspan(EW/ex[1], EW/ex[0], alpha=0.2, color='gray')
			else:
				ax.axvspan(ex[0], ex[1], alpha=0.2, color='gray')
		ax.xlim_init = ax.get_xlim()
	ax.aspan_list = list()
	rangecol = 'red'
	rangealpha = 0.05
	if viewE:
		#wav_max_set = 1e-3 if wav_max==0 else wav_max
		sp1 = max(EW/wav_max,lims[0])
		sp2 = min(EW/wav_min,lims[1])
		ax.aspan_list.append(ax.axvspan(lims[0], sp1, alpha=rangealpha, color=rangecol))
		ax.aspan_list.append(ax.axvspan(sp2, lims[1], alpha=rangealpha, color=rangecol))
	else:
		sp1 = max(ax.xlim_init[0],wav_min)
		sp2 = min(ax.xlim_init[1],wav_max)
		ax.aspan_list.append(ax.axvspan(ax.xlim_init[0], sp1, alpha=rangealpha, color=rangecol))
		ax.aspan_list.append(ax.axvspan(sp2, ax.xlim_init[1], alpha=rangealpha, color=rangecol))
	ax.set_xlim(lims)
	ax.get_figure().canvas.draw()
	#ax.get_figure().canvas.flush_events()



def load_data_file():
	global wav, data
	wav,data = np.genfromtxt(filename,skip_header=skiprows).T
	if xtype == 'cm^-1': 
		wav = 10e6/wav
		wav = np.fliplr(wav)
		data = np.fliplr(data)
	if ytype == 'A': data = 10**(-data)	


def do_plots():	# filename, xtype, ytype
	global plots
	plots = {}
	plots['vlines'] = {}	# crea un dizionario in cui inserire le linee verticali
	# PLOT transmittance
	mask = wav>0	#wav>410
	plots['mask'] = mask
	subpl = 2 if plot_residuals else 1
	plots['ftran'], ax_tran_arr = plt.subplots(subpl,1,dpi=dpi,squeeze=0,figsize=(12,9.4))
	#place figure upper left
	mngr = plots['ftran'].canvas.manager
	#mngr = plt.get_current_fig_manager()
	geom = mngr.window.geometry()
	x,y,dx,dy = geom.getRect()
	mngr.window.setGeometry(0,40,dx, dy)

	x_arr = EW/wav[mask] if viewE else wav[mask]
	ax_tran = ax_tran_arr[0,0]
	plots['ax_tran'] = ax_tran
	ldata = ax_tran.plot(x_arr, data[mask], label='data')[0]
	plots['lfit'] = ax_tran.plot(x_arr, np.zeros_like(x_arr), label='fit')[0]
	trans_model = residual(params,wav[mask])
	plots['lmodel'] = ax_tran.plot(x_arr, trans_model,':', label='model')[0]
	# plt.plot(wav, transmittance_layer(wav,1.4,0,415), label='fixed transmittance model')
	# plt.plot(wav, transmittance_layer(wav,poly_params=params), label='transmittance model')
	ax_tran.set_xscale(xscale)
	ax_tran.set_yscale(yscale)
	plot_excluded_areas(ax_tran)
	x_label = "Energy [eV]" if viewE else "Wavelenght [nm]"
	ax_tran.set_xlabel(x_label)
	ax_tran.set_ylabel("Transmittance")
	ax_tran.legend()
	ax_tran.grid(True)
	ax_tran.set_ylim((1e-8,1.1))
	if plot_residuals:
		ax_res = ax_tran_arr[1,0]
		plots['ax_res'] = ax_res
		if res_view_mode == 'diff':
			res = trans_model-data[mask]
			ax_res.set_ylabel('Model-data')
		elif res_view_mode == 'percent':
			res = (trans_model-data[mask])/data[mask]*100
			ax_res.set_ylabel('Residuals %')
		plots['lres'] = ax_res.plot(x_arr, res,'.', label='residuals')[0]
		ax_res.set_xscale(xscale)
		ax_res.set_yscale(yscale)
		plot_excluded_areas(ax_res)
		ax_res.grid(True)
	plots['ftran'].tight_layout()

	# PLOT refractive index
	if plot_ref_index:
		plots['fri'], [axn, axk] = plt.subplots(2,1,dpi=dpi)
		plots['axn'] = axn
		plots['axk'] = axk
		axn.set_title("n")
		axk.set_title("k")
		axn.grid(True)
		axk.grid(True)
		x_arr = EW/wav if viewE else wav
		#axn.plot(wav, np.real(Forouhi19(wav,out.params)), label='n')
		#axk.plot(wav, -np.imag(Forouhi19(wav,out.params)), label='k')
		ri_PI = Forouhi19w(wav,params.valuesdict())
		plots['ln'] = axn.plot(x_arr, np.real(ri_PI), label='n PI')[0]
		plots['lk'] = axk.plot(x_arr, -np.imag(ri_PI), label='k PI')[0]
		ri_Al = refInd_Al(wav)
		ln = axn.plot(x_arr, np.real(ri_Al), label='n Al')[0]
		lk = axk.plot(x_arr, -np.imag(ri_Al), label='k Al')[0]
		ri_Al2O3 = refInd_Al2O3(wav)
		ln = axn.plot(x_arr, np.real(ri_Al2O3), label='n Al2O3')[0]
		lk = axk.plot(x_arr, -np.imag(ri_Al2O3), label='k Al2O3')[0]
		axn.legend()
		axk.legend()
		
		axn.set_xscale(xscale)
		axk.set_xscale(xscale)
		axn.set_yscale(yscale)
		axk.set_yscale(yscale)
		x_label = "Energy [eV]" if viewE else "Wavelenght [nm]"
		plt.xlabel(x_label)
	plt.show(block=False)

def toggle_view():
	''' Switch from Energy plots to Wavelength plots '''
	global viewE
	viewE = not viewE
	if plot_transmittance:
		plt.close(plots['ftran'])
	if plot_ref_index:
		plt.close(plots['fri'])
	do_plots()
	try:
		plots['lfit'].set_ydata(residual(out.params,wav[plots['mask']]))
	except(NameError):
		pass
	# Update GUI fit limit values
	wind.wave_range_dict['txt_wav_min'].setText(str(round(EW/wav_min,3) if viewE else wav_min))
	wind.wave_range_dict['txt_wav_max'].setText(str(round(EW/wav_max,3) if viewE else wav_max))

def update_plots():
	mask = plots['mask']
	riPI = Forouhi19w(wav,params.valuesdict())
	trans_model = residual(params,wav[mask],riPI=riPI)
	plots['ax_tran'].set_xscale(xscale)
	plots['ax_tran'].set_yscale(yscale)
	plots['lmodel'].set_ydata(trans_model)
	#plots['ax_tran'].relim()
	#plots['ax_tran'].autoscale_view(True,True,True)
	if plot_ref_index:
		plots['lk'].set_ydata(-np.imag(riPI))
		plots['ln'].set_ydata(np.real(riPI))
		#plots['fri'].canvas.blit(axn.bbox)
		#plots['fri'].canvas.blit(axk.bbox)
		plots['fri'].canvas.draw()
		plots['fri'].canvas.flush_events()
		
	if plot_residuals:
		plots['ax_res'].autoscale_view(False,False,False)
		if res_view_mode == 'diff':
			res = trans_model-data[mask]
			plots['ax_res'].set_ylim((-1,1))
		elif res_view_mode == 'percent':
			res = (trans_model-data[mask])/data[mask]*100
			plots['ax_res'].set_ylim((-100,100))
		plots['lres'].set_ydata(res)
		#plots['ax_res'].relim()
		plots['ax_res'].set_xlim(plots['ax_tran'].set_xlim())
		plots['ax_res'].set_xscale(xscale)
		#plots['ax_res'].set_yscale(yscale)
	#plots['fri'].canvas.draw_idle()
	#ftran.canvas.blit(plots['ax_tran'].bbox)
	#plots['ftran'].canvas.draw()
	plots['ftran'].canvas.draw_idle()
	plots['ftran'].canvas.flush_events()


def update_fit_range(valmin, valmax):
	global wav_min, wav_max
	wav_min = EW/valmin if viewE else valmin
	wav_max = EW/valmax if viewE else valmax
	plot_excluded_areas(plots['ax_tran'])
	plot_excluded_areas(plots['ax_res'])


def save_parameters():
	#np.save('fit_results.npy', out, allow_pickle=True)
	filename = gui_SaveFile()
	if '.json' in filename:
		with open(filename, 'w') as f:
			params.dump(f)
			print("Saved current parameters to lmfit dump file")
	elif '.pkl' in filename:
		with open(filename, 'wb') as f:
			pickle.dump(out, f)
			print("Saved output parameters to pickle dump file")


def load_parameters():
	global out, params
	filename = gui_OpenFile()
	print(f"Loading fit results from {filename}")
	#out = np.load('fit_results.npy',allow_pickle=True)
	if '.json' in filename:
		with open(filename, 'r') as f:
			params.load(f)
		paramsToLoad = params
		print("Parameters loaded")
	elif '.pkl' in filename:
		with open(filename, 'rb') as f:
			out = pickle.load(f)
		print("Fit results loaded")
		paramsToLoad = out.params
	get_parameters_from_fit(paramsToLoad=paramsToLoad)
	paramsToLoad.pretty_print()
	plots['lfit'].set_ydata(residual(paramsToLoad,wav[plots['mask']]))

	# Aggiorna la GUI (può essere cambiato il numero di parametri)
	wind.close()
	wind.__init__()
	wind.show()
	update_plots()


####### GUI #################################

def gui_OpenFile(dir=None):
	"""Select a file via a dialog and return the file name."""
	if dir is None: dir ='./'
	#fname = QFileDialog.getOpenFileName(None, "Select fit output file...", 
	#			dir, filter="Pickle data (*.pkl);; All files (*)")
	fname = QFileDialog.getOpenFileName(None, "Select fit output file...", 
				dir, filter="JSON data (*.json);;Pickle data (*.pkl);; All files (*)")
	return fname[0]


def gui_SaveFile(dir=None):
	"""Select a file via a dialog and return the file name."""
	if dir is None: dir ='./'
	#fname = QFileDialog.getSaveFileName(None, "Save output file...", 
	#			dir, filter="Pickle data (*.pkl)")
	fname = QFileDialog.getSaveFileName(None, "Save output file...", 
				dir, filter="JSON data (*.json);;Pickle data (*.pkl)")
	return fname[0]


def update_from_gui(key,src=None):
	'''
	Set a parameter taking it from the GUI
	'''
	slider = wind.slide_dict[key]
	txt_min = wind.txt_min_dict[key]
	txt_max = wind.txt_max_dict[key]
	txt_val = wind.txt_val_dict[key]
	ckb_fix = wind.ckb_fix_dict[key]
	if src == 'txt':
		txt_val.blockSignals(True)
		val = float(txt_val.text())
		txt_val.blockSignals(False)
	else:
		val = slider.value()
	txt_val.setText(f"{val:.3e}")
	val_min, val_max = float(txt_min.text()), float(txt_max.text())
	slider.setMinimum(val_min)
	slider.setMaximum(val_max)
	slider.setValue(val)
	params[key].set(value=val, min=val_min, max=val_max, vary= not ckb_fix.isChecked())
	# Aggiorna i parametri che implementano i vincoli:
	if key[:2] == 'Eh' or key[:2] == 'Wh':
		DEh = val-params['Eg'].value
		params['DEh'+key[2:]].set(value=DEh)
		DWh = 5.65685*(params['Eh'+key[2:]].value - params['Eg'].value) - params['Wh'+key[2:]].value
		params['DWh'+key[2:]].set(value=DWh)
		Wmax = 4*np.sqrt(2)*(params['Eh'+key[2:]].value - params['Eg'].value) # W < 4*np.sqrt(2)*(Ep - Eg)
		keyWh = 'Wh'+key[2:]
		params[keyWh].max = Wmax	# aggiorna il massimo di Wh in modo che rispetti il limite
		wind.slide_dict[keyWh].setMaximum(Wmax)
		wind.txt_max_dict[keyWh].setText(f"{Wmax:.2e}")
	elif key == 'Eg':	# se stiamo modificando Eg dobbiamo aggiornare tutti i parametri che ne sono vincolati
		Eg = val
		for key1 in params.keys():
			if key1[:2] == 'Eh' or key1[:2] == 'Wh':
				DEh = val-params['Eg'].value
				params['DEh'+key1[2:]].set(value=DEh)
				DWh = 5.65685*(params['Eh'+key1[2:]].value - params['Eg'].value) - params['Wh'+key1[2:]].value
				params['DWh'+key1[2:]].set(value=DWh)
				## Il codice sottostante serve per aggiornare i parametri Eh e Wh nella GUI in modo che rispettino i vincoli. Tuttavia ciò potrabbe portare a variazioni incontrollate di questi parametri al variare di Eg, meglio avere soltanto un avviso.
				# if key1[:2] == 'Eh'
					# if params[key1].min   <= val:	params[key1].min   = Eg+ε
					# if params[key1].value <= val:	params[key1].value = Eg+ε
					# if params[key1].max   <= val:	params[key1].max   = Eg+ε
				# elif key1[:2] == 'Wh':
					# Wmax = 4*np.sqrt(2)*(params['Eh'+key1[2:]].value - Eg) # W < 4*np.sqrt(2)*(Ep - Eg)
					# if params[key1].min   >= Wmax:	params[key1].min   = Wmax-ε
					# if params[key1].value >= Wmax:	params[key1].value = Wmax-ε
					# if params[key1].max   >= Wmax:	params[key1].max   = Wmax-ε
		#get_parameters_from_fit(paramsToLoad=params)
	update_plots()



def toggleFix():
	for key,par in params.items():
		if not key in excl_pars:
			par.vary = not par.vary
			ckb_fix = wind.ckb_fix_dict[key]
			ckb_fix.toggle()


def fix(setTo = True):
	for key,par in params.items():
		if not key in excl_pars:
			par.vary = setTo
			ckb_fix = wind.ckb_fix_dict[key]
			ckb_fix.setChecked(setTo)


def get_parameters_from_fit(paramsToLoad=None):
	'''
	Get parameters from fit and load them in the GUI
	'''
	global params
	if not paramsToLoad:
		try:
			#print(out.params)
			#params.__dict__.update(out.params.__dict__)
			print(params['Pl02'])
			print(out.params['Pl02'])
			params = copy.deepcopy(out.params)	# fit results are copied in the live parameters
			print(params['Pl02'])
			print('Fit parameters loaded')
			#pars = out.params					# read only from the fit results to populate GUI
		except(NameError):
			print("Fit output not available")
			return
	else:
		#params.__dict__.update(paramsToLoad.__dict__)
		params = copy.deepcopy(paramsToLoad)
		print('Input parameters loaded')
	for key,val in wind.parvals.items():
		slider = wind.slide_dict[key]
		txt_min = wind.txt_min_dict[key]
		txt_max = wind.txt_max_dict[key]
		txt_val = wind.txt_val_dict[key]
		ckb_fix = wind.ckb_fix_dict[key]
		# block signals to speed up loading
		txt_val.blockSignals(True)
		slider.blockSignals(True)
		txt_min.blockSignals(True)
		txt_max.blockSignals(True)
		txt_val.setText(f"{params[key].value:.3e}")
		#print(f"{key} set to: {params[key].value:.3e}") 
		slider.setValue(params[key].value)
		slider.setMinimum(params[key].min)
		slider.setMaximum(params[key].max)
		txt_min.setText(f"{params[key].min:.2e}")
		txt_max.setText(f"{params[key].max:.2e}")
		ckb_fix.setChecked(not params[key].vary)
		# re-enable signals
		txt_val.blockSignals(False)
		slider.blockSignals(False)
		txt_min.blockSignals(False)
		txt_max.blockSignals(False)
	update_plots()
	#print('get_parameters_from_fit done')


def load_parameters_from_file():
	''' Load parameters from a file, both for fitting and in the GUI '''
	load_parameters()
	get_parameters_from_fit()


def plot_vertical_lines(key):
	global plots
	if viewE == False: return
	if not key[0] in ['E', 'P', 'W']: return
	linekey = 'E'+key[1:]
	Eval = float(wind.txt_val_dict[linekey].text())
	lbl = wind.lbl_dict[linekey]
	if plots['vlines'].get(linekey,None) is None:	# la specifica linea non esiste, viene creata
		lbl.setStyleSheet("background-color: lightcoral; border: 1px solid black;")
		plots['vlines'][linekey] = plots['ax_tran'].vlines(Eval,0,1,'r')
	else:	# la linea esiste, va eliminata
		lbl.setStyleSheet("background-color: transparent; border: 0px;")
		plots['vlines'][linekey].remove()
		del plots['vlines'][linekey]
	plots['ftran'].canvas.draw()
	plots['ftran'].canvas.flush_events()


def add_fit_term():
	'''
	Sequenza:
	- finestra di input per indicare il tipo di termine, a che materiale appartiene, i parametri relativi (i.e. E, P, W, X). Opzionalmente click del mouse sul grafico per selezionare l'energia.
	- aggiunge i parametri alla lista opportuna (es. ppEh o ppEl)
	- aggiunge i parametri all'oggetto params:
		params.add('El'+str(j), value=ppEl[j]['E'][0], min=ppEl[j]['E'][1], max=ppEl[j]['E'][2], vary=False)
	- ricostruisce la finestra di controllo:
		wind.close()
		wind.__init__()
		wind.show()
	In alternativa alla finestra di input si può fare un menu a tendina in cui si seleziona tipo termine e materiale, il parametro viene aggiunto con valori di default e l'utente li modifica direttamente dalla finestra di controllo.
	
	
	params.add('El'+ppEl[-1]['lbl'], value=ppEl[-1]['E'][0], min=ppEl[-1]['E'][1], max=ppEl[-1]['E'][2], vary=False)
	params.add('Pl'+ppEl[-1]['lbl'], value=ppEl[-1]['P'][0], min=ppEl[-1]['P'][1], max=ppEl[-1]['P'][2], vary=False)
	params.add('Wl'+ppEl[-1]['lbl'], value=ppEl[-1]['W'][0], min=ppEl[-1]['W'][1], max=ppEl[-1]['W'][2], vary=False)
	'''
	ppEl.append({'E':(0.05301,1e-3,6.5), 'P':(0.1,0,1), 'W':(1e-3,0,1e-2), 'lbl':str(len(ppEl))})
	# Termine amorfo E<Eg
	params.add('El'+str(len(ppEl)-1), value=0.3, min=1e-3, max=3.5, vary=True)
	params.add('Pl'+str(len(ppEl)-1), value=1, min=0, max=1, vary=True)
	params.add('Wl'+str(len(ppEl)-1), value=0.02, min=0, max=1e-2, vary=True)
	
	wind.close()
	wind.__init__()
	wind.show()
	update_plots()

class DoubleSlider(QSlider):
	def __init__(self, *args, **kwargs):
		super().__init__(*args, **kwargs)
		self.decimals = 5
		self._max_int = 10 ** self.decimals
		super().setMinimum(0)
		super().setMaximum(self._max_int)
		self._min_value = 0.0
		self._max_value = 1.0
	@property
	def _value_range(self):
		return self._max_value - self._min_value
	def value(self):
		return float(super().value()) / self._max_int * self._value_range + self._min_value
	def setValue(self, value):
		super().setValue(int((value - self._min_value) / self._value_range * self._max_int))
	def setMinimum(self, value):
		if value > self._max_value:
			raise ValueError("Minimum limit cannot be higher than maximum")
		self._min_value = value
		#self.setValue(self.value())
	def setMaximum(self, value):
		if value < self._min_value:
			raise ValueError("Minimum limit cannot be higher than maximum")
		self._max_value = value
		#self.setValue(self.value())
	def minimum(self):
		return self._min_value
	def maximum(self):
		return self._max_value



#class Window(QWidget):
class Window(QMainWindow):
#class Window(QScrollArea):
	def __init__(self, parent=None):
		super(Window, self).__init__(parent)

		# Crea un dizionario parametri escludendone alcuni su cui non si vuole controllo
		parvals_all = params.valuesdict()
		self.parvals = dict()
		for key,val in parvals_all.items():
			key_no_digits = ''.join([i for i in key if not i.isdigit()])
			if not key_no_digits in excl_pars:
				self.parvals[key] = val
		npar = len(self.parvals)

		self.setWindowTitle("Parameters control")
		#self.resize(475, 850)
		self.resize(720, 980)
		ag = QDesktopWidget().availableGeometry()
		sg = QDesktopWidget().screenGeometry()
		x = ag.width() - self.geometry().width()
		y = 0
		self.move(x, y)
		
		# Crea dizionari per rendere accessibili vari elementi del GUI in seguito
		self.lbl_dict = {}
		self.slide_dict = {}
		self.txt_min_dict = {}
		self.txt_max_dict = {}
		self.txt_val_dict = {}
		self.ckb_fix_dict = {}
		self.wave_range_dict = {}


		# Pulsanti Get, Load e Save
		btn_getPar = QPushButton('Get from fit')
		btn_getPar.clicked.connect(get_parameters_from_fit)
		btn_load = QPushButton('Load')
		btn_load.clicked.connect(load_parameters_from_file)
		btn_save = QPushButton('Save')
		btn_save.clicked.connect(save_parameters)
		hbox_buttons = QHBoxLayout()
		hbox_buttons.addWidget(btn_getPar)
		hbox_buttons.addWidget(btn_load)
		hbox_buttons.addWidget(btn_save)
		widg_buttons = QWidget()
		widg_buttons.setLayout(hbox_buttons)
		#grid.addWidget(widg_buttons)

		def changePlotScale(xlog=None, ylog=None):
			global xscale, yscale
			if not xlog is None:
				xscale = 'log' if xlog else 'linear'
			if not ylog is None:
				yscale = 'log' if ylog else 'linear'
			update_plots()

		# Checkbox plot
		ckb_xlog = QCheckBox("X log")
		ckb_ylog = QCheckBox("Y log")
		ckb_xlog.setChecked(True if xscale == 'log' else False)
		ckb_ylog.setChecked(True if yscale == 'log' else False)
		ckb_xlog.stateChanged.connect(lambda: changePlotScale(xlog=ckb_xlog.isChecked()))
		ckb_ylog.stateChanged.connect(lambda: changePlotScale(ylog=ckb_ylog.isChecked()))
		btn_toggle_view = QPushButton('E/λ')
		btn_toggle_view.clicked.connect(toggle_view)
		#grid = QGridLayout()
		#grid.addWidget(btn_toggle_fix)
		
		widg_scale = QWidget()
		#hbox_scale = QHBoxLayout()
		hbox_scale = QHBoxLayout(widg_scale)
		hbox_scale.addStretch(1)
		hbox_scale.addWidget(ckb_xlog)
		hbox_scale.addWidget(ckb_ylog)
		hbox_scale.addWidget(btn_toggle_view)
		hbox_scale.addStretch(1)
		hbox_scale.setContentsMargins(0, 0, 0, 0)
		#widg_scale.setLayout(hbox_scale)
		#grid.addWidget(widg_scale)

		# Intervallo di fit e pulsante di fit
		lbl_wav_min = QLabel("E max" if viewE else "λ min")
		self.lbl_wav_min = lbl_wav_min
		txt_wav_min = QLineEdit()
		txt_wav_min.setFixedWidth(60)
		txt_wav_min.setText(str(round(EW/wav_min,3) if viewE else wav_min))
		txt_wav_min.editingFinished.connect(lambda: update_fit_range(float(txt_wav_min.text()),float(txt_wav_max.text())))
		self.wave_range_dict['txt_wav_min'] = txt_wav_min
		lbl_wav_max = QLabel("E min" if viewE else "λ max")
		self.lbl_wav_max = lbl_wav_max
		txt_wav_max = QLineEdit()
		txt_wav_max.setFixedWidth(60)
		txt_wav_max.setText(str(round(EW/wav_max,3) if viewE else wav_max))
		txt_wav_max.editingFinished.connect(lambda: update_fit_range(float(txt_wav_min.text()),float(txt_wav_max.text())))
		self.wave_range_dict['txt_wav_max'] = txt_wav_max
		#txt_wav_max.editingFinished.connect(lambda: (wav_max:=float(txt_wav_max.text())))
		btn_fit = QPushButton('Fit')
		btn_fit.clicked.connect(fit)
		btn_addTerm = QPushButton('Add El term')
		btn_addTerm.clicked.connect(add_fit_term)
		widg_wav = QWidget()
		#hbox_wav = QHBoxLayout()
		hbox_wav = QHBoxLayout(widg_wav)
		hbox_wav.addWidget(lbl_wav_min)
		hbox_wav.addWidget(txt_wav_min)
		hbox_wav.addStretch(1)
		hbox_wav.addWidget(lbl_wav_max)
		hbox_wav.addWidget(txt_wav_max)
		hbox_wav.addWidget(btn_fit)
		hbox_wav.addWidget(btn_addTerm)
		#widg_wav.setLayout(hbox_wav)
		#grid.addWidget(widg_wav)
		
		widg_commands = QWidget()
		vbox_commands = QVBoxLayout(widg_commands)
		vbox_commands.addWidget(widg_buttons)
		vbox_commands.addWidget(widg_scale)
		vbox_commands.addWidget(widg_wav)
		#widg_commands.setLayout(vbox_commands)
		
		btn_toggle_fix = QPushButton('Toggle Fix')
		btn_toggle_fix.clicked.connect(toggleFix)
		
		btn_fix_all = QPushButton('Fix all')
		btn_fix_all.clicked.connect(lambda: fix(True))
		
		btn_fix_none = QPushButton('Unfix all')
		btn_fix_none.clicked.connect(lambda: fix(False))
		
		widg_fix = QWidget()
		hbox_fix = QHBoxLayout(widg_fix)
		hbox_fix.addWidget(btn_toggle_fix)
		hbox_fix.addWidget(btn_fix_all)
		hbox_fix.addWidget(btn_fix_none)
		
		grid = QGridLayout()
		grid.addWidget(widg_fix)
		
		# Create controls for all the fit parameters
		#i=2
		for key,val in self.parvals.items():
			#grid.addWidget(self.createParameterGroup(key=key), i, 0)
			grid.addWidget(self.createParameterGroup(key=key))
		#	i += 1
		grid.setSpacing(0)
		
		widg_scroll = QWidget()
		widg_scroll.setLayout(grid)
		#scrollarea = QScrollArea(widg_scroll)
		scrollarea = QScrollArea()
		scrollarea.setWidget(widg_scroll)
		scrollarea.setWidgetResizable(True)
		#scrollarea.setFixedHeight(200)
		
		widg_central = QWidget()
		vbox_central = QVBoxLayout(widg_central)
		#vbox_central.addLayout(hbox_buttons)
		vbox_central.addWidget(widg_commands)
		vbox_central.addWidget(scrollarea)
		#vbox_central.addStretch(1)
		
		self.setCentralWidget(widg_central)
		#self.setCentralWidget(scrollarea)
		

	def createParameterGroup(self,key=''):
		vals = (params[key].min, params[key].max, params[key].value)
		#groupBox = QGroupBox(key)
		groupBox = QFrame()
		lbl_key = QLabel(f"{key:5s}")
		lbl_key.mousePressEvent = lambda event: plot_vertical_lines(key)
		slider = DoubleSlider(Qt.Horizontal)
		slider.setFocusPolicy(Qt.StrongFocus)
		#slider.setTickPosition(QSlider.TicksAbove)
		#slider.setTickInterval(10)
		#slider.setSingleStep(1)
		slider.setMaximum(vals[1])
		slider.setMinimum(vals[0])
		slider.setValue(vals[2])
		txt_min = QLineEdit()				# minimo dello slider
		txt_min.setText(f"{vals[0]:.2e}")
		txt_min.setFixedWidth(80)
		txt_max = QLineEdit()				# massimo dello slider
		txt_max.setText(f"{vals[1]:.2e}")
		txt_max.setFixedWidth(80)
		txt_val = QLineEdit()	# indicatore di valore
		txt_val.setText(f"{vals[2]:.3e}")
		txt_val.setFixedWidth(90)
		ckb_fix = QCheckBox("Fix")
		ckb_fix.setChecked(not params[key].vary)
		
		hbox = QHBoxLayout()
		hbox.addWidget(lbl_key)
		hbox.addWidget(txt_min)
		#hbox.addStretch(1)
		hbox.addWidget(slider)
		hbox.addWidget(txt_val)
		hbox.addWidget(txt_max)
		hbox.addWidget(ckb_fix)
		groupBox.setLayout(hbox)
		
		self.lbl_dict[key] = lbl_key
		self.slide_dict[key] = slider
		self.txt_min_dict[key] = txt_min
		self.txt_max_dict[key] = txt_max
		self.txt_val_dict[key] = txt_val
		self.ckb_fix_dict[key] = ckb_fix
		
		slider.valueChanged.connect(lambda: update_from_gui(key))
		txt_min.editingFinished.connect(lambda: update_from_gui(key))
		txt_max.editingFinished.connect(lambda: update_from_gui(key))
		txt_val.editingFinished.connect(lambda: update_from_gui(key,src='txt'))
		
		def change_fix(key,val):
			params[key].vary=val
			#print(f"{key} changed to {val}")
			#print(params[key])
			#print(params[key].vary)
			
			
		ckb_fix.stateChanged.connect(lambda: change_fix(key,(not ckb_fix.isChecked())))
		return groupBox

def printOutParam(key = None):
	if not key:
		print(out.params)
	else:
		print(out.params[key])
def printParam(key = None):
	if not key:
		print(params)
	else:
		print(params[key])


if __name__ == '__main__':
	app = QtCore.QCoreApplication.instance()
	if app is None:
		app = QtWidgets.QApplication(sys.argv)
	#app = QApplication(sys.argv)
	load_data_file()
	do_plots()
	wind = Window()
	wind.show()
	#sys.exit(app.exec_())
	
'''
PROFILING
import cProfile, pstats
cProfile.run("update_from_gui('El0')",'cProfile.stats')
s = pstats.Stats('cProfile.stats').sort_stats('cumulative') # .reverse_order()
s.print_stats(0.1)	# stampa il 10% della lista
s.print_stats('n-fit')	# stampa solo le chiamate a funzioni che contengono n-fit nel nome
'''