# -*- coding: utf-8 -*-
"""
Created on Fri Oct 14 10:24:08 2022

@author: wicki
"""

# -*- coding: utf-8 -*-
"""
Created on Fri Oct 14 00:14:05 2022

@author: wicki
"""

# -*- coding: utf-8 -*-
"""
Created on Wed Jun 15 15:52:13 2022

@author: wicki
"""


import numpy as np
import matplotlib.pyplot as plt
#from scipy.optimize import curve_fit
#Need to import at least FitObject and fit_parameter
from fitutils_P3 import FitObject, fit_parameter
from matplotlib import rc
rc('font', **{'family': 'sans-serif', 'sans-serif': ['arial']})

plt.rcParams['mathtext.it'] = 'Arial:italic'
plt.rcParams['mathtext.rm'] = 'Arial'
plt.rcParams['svg.fonttype'] = 'none'



fig, ((ax111,ax222,ax333) ) = plt.subplots(nrows=3,ncols=1, figsize = (2.6,8))
#fig, ((ax111,ax1,ax15) ,(ax222,ax2,ax25),(ax333,ax3,ax35)) = plt.subplots(nrows=3,ncols=3, figsize = (10,8))
plt.subplots_adjust(hspace=0.4)
dataw1=np.genfromtxt('di_spectra_3.csv', dtype=float, delimiter=',', skip_header=1,names=['ppm','khz','b','c','d'])
ppm = dataw1['ppm']
khz = dataw1['khz']
b = dataw1['b']
c = dataw1['c']
d = dataw1['d']



dataw3=np.genfromtxt('tri_SE_spec_3.csv', dtype=float, delimiter=',',skip_header=1, names=['ppmt','bt','ct','dt'])
ppmt = dataw3['ppmt']
bt = dataw3['bt']
ct = dataw3['ct']
dt = dataw3['dt']

fs=10


def rawgauss(x, scale, X, sigma):
    #zvalues = (x - X)/(1.4142*sigma)
    #return scale* np.exp(-(zvalues**2))
   # return scale*1/(sigma*(2*3.1415926)**0.5)*np.exp(-((x-X)**2/(2*sigma**2)))
    return scale*np.exp(-((x-X)**2/(2*(sigma**2))))
# Functions to be objted
# These have the "signature" of xvalues, dictionary of parameters, optional additional (fixed) arguments

def fitgauss1(x, pardict, args):
    X = pardict['X'].value
    sigma = pardict['sigma'].value
    scale = pardict['scale'].value
    return rawgauss(x, scale, X, sigma)

def fitgauss2(x, pardict, args):
    X1 = pardict['X1'].value
    sigma1 = pardict['sigma1'].value
    scale1 = pardict['scale1'].value
    X2 = pardict['X2'].value
    sigma2 = pardict['sigma2'].value
    scale2 = pardict['scale2'].value
    return rawgauss(x, scale1, X1, sigma1) + rawgauss(x, scale2, X2, sigma2)

# Initialization parameters
xdc=khz  #convert ppm to kHz

parscale1 = fit_parameter('scale', 20000.)
parX1 = fit_parameter('X', 1., displayname='X / kHz')
parsigma1 = fit_parameter('sigma', 1., displayname='sigma / kHz')

parscale2 = fit_parameter('scale', 400000.)
parX2 = fit_parameter('X', 0., displayname='X / kHz')
parsigma2 = fit_parameter('sigma', 1., displayname='sigma / kHz')

parscale12 = fit_parameter('scale', 20000.)
parX12 = fit_parameter('X', 1., displayname='X / kHz')
parsigma12 = fit_parameter('sigma', 1., displayname='sigma / kHz')

parscale22 = fit_parameter('scale', 400000.)
parX22 = fit_parameter('X', 0., displayname='X / kHz')
parsigma22 = fit_parameter('sigma', 1., displayname='sigma / kHz')

parscale3 = fit_parameter('scale1', 40000.)
parX3 = fit_parameter('X1', 0., displayname='X / kHz')
parsigma3 = fit_parameter('sigma1', 1., displayname='sigma / kHz')


parscale4 = fit_parameter('scale2', 40000.)
parX4 = fit_parameter('X2', 0., displayname='X / kHz')
parsigma4 = fit_parameter('sigma2', 1., displayname='sigma / kHz')
#parc = fit_parameter('c', 10000., displayname='c ')

#Create fitting object, initialised with fitting function, list of parameters to fit

fitobj1 = FitObject(fitgauss1, [parscale1, parX1, parsigma1])
fitobj2 = FitObject(fitgauss1, [parscale2, parX2, parsigma2])
fitobj3 = FitObject(fitgauss2, [parscale3, parX3, parsigma3,parscale4, parX4, parsigma4])
#fitobj22 = FitObject(fitgauss2, [parscale1, parX1, parsigma1,parscale2, parX2, parsigma2])
#fitobj = FitObject(fitgauss2, [parscale1, parX1, parsigma1,parscale2, parX2, parsigma2, parc])

r1=b
r2=c
r3=d
xdct=ppmt*400.14/1000
fitdata1 = fitobj1.fit(xdc, r1)
fitdata2 = fitobj2.fit(xdc, r2)
fitdata3 = fitobj3.fit(xdc, r3)

# Plot the results
# Data


r1=b
r2=c
r3=d

fitdata1 = fitobj1.fit(xdc, r1)
fitdata2 = fitobj2.fit(xdc, r2)
fitdata3 = fitobj3.fit(xdct, ct)

# Plot the results
# Data
ax111.plot(xdc, r1, 'k', label='Data')
ax222.plot(xdc, r2, 'k', label='Data')
ax333.plot(xdct, ct, 'k', label='Data') #1#2
#ax333.plot(xdcc, r65, 'k', label='Data')#3

ax111.plot(xdc, fitdata1, 'r-', label='Fit')
ax222.plot(xdc, fitdata2, 'r-', label='Fit')
ax333.plot(xdct, fitdata3, 'r-', label='Fit',alpha=0)#1
ax333.plot(xdc, rawgauss (xdc,197778.961186,0.729471,10.871957), 'r', label='Broad signal',zorder=5)#2
ax333.plot(xdc, rawgauss (xdc,120057.350475,0.754236,0.686342), 'r', label='Narrow signal',alpha=0.5)#2
#ax333.plot(xdc, rawgauss (xdc,120057.350475,0.754236,0.686342)+rawgauss (xdc,197778.961186,0.729471,10.871957), 'g', label='Narrow signal')#2
#fitdata2 = fitobj2.fit(xdcc, r65)#3
#ax333.plot(xdcc, fitdata2, 'r-', label='Fit')#3

#ax333.plot(xdct, fitdata3, 'g-', label='Fit')



#plt.plot(xdc, rawgauss(xdc, 1100772, 1.5, 8.0))
limp1=70
limp2=70
limp3=70

ax111.set_xlim(-limp1,limp1)
ax222.set_xlim(-limp2,limp2)
ax333.set_xlim(-limp3,limp3)


ax333.set_xlabel('Frequency / kHz ', fontsize=fs)



ax111.axes.yaxis.set_visible(False)
ax222.axes.yaxis.set_visible(False)
ax333.axes.yaxis.set_visible(False)

ax111.tick_params(axis='both', which='major', labelsize=fs)
ax222.tick_params(axis='both', which='major', labelsize=fs)
ax333.tick_params(axis='both', which='major', labelsize=fs)



ax111.invert_xaxis()
ax222.invert_xaxis()
ax333.invert_xaxis()


ax111.text(0.5, 1.08, u'Diamantane, 221 K',
        verticalalignment='bottom', horizontalalignment='center',
        transform=ax111.transAxes,
        color='black', fontsize=fs)
ax222.text(0.5, 1.08, u'Diamantane, 411 K',
        verticalalignment='bottom', horizontalalignment='center',
        transform=ax222.transAxes,
        color='black', fontsize=fs)
ax333.text(0.5, 1.08, u'Triamantane, 360 K',
        verticalalignment='bottom', horizontalalignment='center',
        transform=ax333.transAxes,
        color='black', fontsize=fs)

#plt.tight_layout()
plt.savefig('GaussianFit_SE2.png', bbox_inches='tight',dpi=400)
plt.savefig('GaussianFit_SE2.svg', bbox_inches='tight',dpi=400)
plt.show()