# -*- coding: utf-8 -*-
"""
Created on Fri Oct 21 23:10:09 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,ax1,ax1111,ax15) ,(ax222,ax2,ax2222,ax25),(ax333,ax3,ax3333,ax35)) = plt.subplots(nrows=3,ncols=4, figsize = (14,8))
fig, ((ax111,ax1,ax15) ,(ax222,ax2,ax25),(ax333,ax3,ax35),(ax111t,ax1t,ax15t) ,(ax222t,ax2t,ax25t),(ax333t,ax3t,ax35t)) = plt.subplots(nrows=6,ncols=3, figsize = (10,13))#,constrained_layout=True)
plt.subplots_adjust(wspace=0.55, hspace = 0.3)
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']



datam=np.genfromtxt('di_FIDs_3_crop.csv', dtype=float, delimiter=',', names=['af','bf','cf','df'],skip_header=2)


af=datam['af']
bf=datam['bf'] 
cf=datam['cf'] 
df=datam['df'] 

dataw1t=np.genfromtxt('tri_SE_spec_3.csv', dtype=float, delimiter=',',skip_header=1, names=['ppmt','bt','ct','dt'])
ppmt = dataw1t['ppmt']
bt = dataw1t['bt']
ct = dataw1t['ct']
dt = dataw1t['dt']
dataw2t=np.genfromtxt('tri_SE_M2_cumulative.csv', dtype=float, delimiter=',', names=['x11t','i11t', 'y11t','x22t','i22t','y22t','x33t','i33t','y33t'], skip_header=1)
x11t = dataw2t['x11t']
i11t = dataw2t['i11t']
y11t = dataw2t['y11t']
x22t = dataw2t['x22t']
i22t = dataw2t['i22t']
y22t = dataw2t['y22t']
x33t = dataw2t['x33t']
i33t = dataw2t['i33t']
y33t = dataw2t['y33t']


datamt=np.genfromtxt('tri_SE_FIDs.csv', dtype=float, delimiter=',', names=['aft','bft','cft','dft'],skip_header=2)
aft=datamt['aft']
bft=datamt['bft'] 
cft=datamt['cft'] 
dft=datamt['dft'] 

fs=14


ax5 = ax1.twinx()
ax6 = ax2.twinx()
ax7 = ax3.twinx()

ax9 = ax1.twinx()
ax10 = ax2.twinx()
ax11 = ax3.twinx()


ax9.yaxis.tick_left()
ax10.yaxis.tick_left()
ax11.yaxis.tick_left()


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 fitted
# 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')

parscale3 = fit_parameter('scale', 40000.)
parX3 = fit_parameter('X', 0., displayname='X / kHz')
parsigma3 = fit_parameter('sigma', 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(fitgauss1, [parscale3, parX3, parsigma3])
#fitobj = FitObject(fitgauss2, [parscale1, parX1, parsigma1,parscale2, parX2, parsigma2, parc])

r1=b
r2=c
r3=d

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

# Plot the results
# Data
ax1.plot(xdc, r1, 'k', label='Data')
ax2.plot(xdc, r2, 'k', label='Data')
ax3.plot(xdc, r3, 'k', label='Data')




r1=b
r2=c
r3=d

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

# Plot the results
# Data
ax111.plot(xdc, r1, 'k', label='Data')
ax222.plot(xdc, r2, 'k', label='Data')
ax333.plot(xdc, r3, 'k', label='Data')


#plt.plot(xd, b1, 'g', label='Data im')
# Fitted function
ax111.plot(xdc, fitdata1, 'r-', label='Fit')
ax222.plot(xdc, fitdata2, 'r-', label='Fit')
ax333.plot(xdc, fitdata3, 'r-', label='Fit')

#plt.plot(xdc, rawgauss(xdc, 1100772, 1.5, 8.0))
limp1=70
limp2=70
limp3=15
fs=14
ax111.set_xlim(-limp1,limp1)
ax222.set_xlim(-limp2,limp2)
ax333.set_xlim(-7,9)
#plt.xlim(-200,-100)
#plt.ylim(-5000,5000)
#plt.legend(frameon=False, prop={'size':12} )


#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)
#plt.savefig('Gaussian_fit_P3.png', dpi = 300)
#plt.plot((fit[0]+fit[2]/1.1774, fit[0]-fit[2]/1.1774), (fit[1]/2, fit[1]/2), 'k-')
#plt.axvline(x=0)


#plt.plot(xdc, rawgauss(xdc, scale1, x1, sigma1), 'm-', label='Fit')
#plt.plot(xdc, rawgauss(xdc, parscale1.value, parX1.value, parsigma1.value)+parc.value, 'k-', label='Fit')
#plt.plot(xdc, rawgauss(xdc, parscale2.value, parX2.value, parsigma2.value)+parc.value, 'k-', label='Fit')

#plt.plot(xdc, r1-parc.value-rawgauss(xdc, parscale1.value, parX1.value, parsigma1.value), 'g', label='Data')
#plt.plot(xdc, r1-parc.value-rawgauss(xdc, parscale2.value, parX2.value, parsigma2.value), 'g', label='Data')

#plt.plot((0,0), (-20,0), 'r-')

#plt.plot((-60,60), (0,0), 'k-')
"""
ax111.text(0.05, 0.8, u'$-$52 \u00B0C',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax1.transAxes,
        color='black', fontsize=fs)
ax222.text(0.05, 0.8, u'138 \u00B0C',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax2.transAxes,
        color='black', fontsize=fs)
ax333.text(0.05, 0.8, u'225 \u00B0C',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax3.transAxes,
        color='black', fontsize=fs)
"""
#to here

ax5.set_ylabel('M2 / kHz2', fontsize=fs)
ax6.set_ylabel('M2 / kHz2', fontsize=fs)
ax7.set_ylabel('M2 / kHz2', fontsize=fs)

ax9.set_ylabel('Integral / kHz ', fontsize=fs)
ax10.set_ylabel('Integral / kHz ', fontsize=fs)
ax11.set_ylabel('Integral / kHz ', fontsize=fs)

ax9.yaxis.set_label_position("left")
ax10.yaxis.set_label_position("left")
ax11.yaxis.set_label_position("left")


ax5.tick_params(axis='y', which='major', labelsize=fs)#
ax6.tick_params(axis='y', which='major', labelsize=fs)
ax7.tick_params(axis='y', which='major', labelsize=fs)
ax1.tick_params(axis='x', which='major', labelsize=fs)#
ax2.tick_params(axis='x', which='major', labelsize=fs)
ax3.tick_params(axis='x', which='major', labelsize=fs)
ax9.tick_params(axis='y', which='major', labelsize=fs)#
ax10.tick_params(axis='y', which='major', labelsize=fs)
ax11.tick_params(axis='y', which='major', labelsize=fs)
#ax11.set_yticks(0,1,0.5)

ax5.yaxis.label.set_color('blue')
ax6.yaxis.label.set_color('blue')
ax7.yaxis.label.set_color('blue')
ax9.yaxis.label.set_color('green')
ax10.yaxis.label.set_color('green')
ax11.yaxis.label.set_color('green')

ax5.tick_params(axis='y', colors='blue')
ax6.tick_params(axis='y', colors='blue')
ax7.tick_params(axis='y', colors='blue')
ax9.tick_params(axis='y', colors='green')
ax10.tick_params(axis='y', colors='green')
ax11.tick_params(axis='y', colors='green')

#plt.plot(xd, b1, 'g', label='Data im')
# Fitted function
#ax1.plot(xdc, fitdata1, 'r-', label='Fit')
#ax2.plot(xdc, fitdata2, 'r-', label='Fit')
#ax3.plot(xdc, fitdata3, 'r-', label='Fit')

#plt.plot(xdc, rawgauss(xdc, 1100772, 1.5, 8.0))
limp1=70
limp2=70
limp3=8
fs=14
ax1.set_xlim(-limp1,limp1)

ax2.set_xlim(-limp2,limp2)
ax3.set_xlim(-7,9)

ppp = -0.05
qqq = 1.1
ax9.set_ylim(ppp,qqq)
ax10.set_ylim(ppp,qqq)
ax11.set_ylim(ppp,qqq)
aaa=160
bbb=85
ccc=1.1
ax5.set_ylim(-aaa/22,aaa)
ax6.set_ylim(-bbb/22,bbb)
ax7.set_ylim(-ccc/22,ccc)

ddd=105000
eee=91000
fff=1500000

ax1.set_ylim(-ddd/22,ddd)
ax2.set_ylim(-eee/22,eee)
ax3.set_ylim(-fff/22,fff)


#plt.xlim(-200,-100)
#plt.ylim(-5000,5000)
#plt.legend(frameon=False, prop={'size':12} )

ax1.invert_xaxis()
ax2.invert_xaxis()
ax3.invert_xaxis()

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

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

#frame1 = plt.gca()
#frame1.axes.yaxis.set_ticklabels([])

#frame2 = plt.gca()
#frame2.axes.yaxis.set_ticklabels([])

#frame3 = plt.gca()
#frame3.axes.yaxis.set_ticklabels([])

ax1.axes.yaxis.set_visible(False)
ax2.axes.yaxis.set_visible(False)
ax3.axes.yaxis.set_visible(False)

ax1.tick_params(axis='both', which='major', labelsize=fs)
ax2.tick_params(axis='both', which='major', labelsize=fs)
ax3.tick_params(axis='both', which='major', labelsize=fs)
 
x0001=af
y0001=bf
x0002=af
y0002=cf
x0003=af
y0003=df

ulA=0.001
ulB=0.001
ulC=0.01

ll=0.000002

#t0s = t0[t0<=ul]


#FID0s = t0[t0<=ul]


x000aa = x0001[x0001<=ulA]
x000a = x000aa[x000aa>=ll] 

x000bb = x0002[x0002<=ulB]
x000b = x000bb[x000bb>=ll] 

x000cc = x0003[x0003<=ulC]
x000c = x000cc[x000cc>=ll] 

y000aa = y0001[x0001<=ulA]
y000a = y000aa[x000aa>=ll]
y000bb = y0002[x0002<=ulB]
y000b = y000bb[x000bb>=ll]
y000cc = y0003[x0003<=ulC]
y000c = y000cc[x000cc>=ll]
parscale1 = fit_parameter('scale', 4000.)
parX1 = fit_parameter('X', 0., displayname='X')
parsigma1 = fit_parameter('sigma', 0.000002, displayname='sigma')

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

parscale3 = fit_parameter('scale', 30000.)
parX3 = fit_parameter('X', 0., displayname='X / kHz')
parsigma3 = fit_parameter('sigma', 0.0005, displayname='sigma / kHz')

def rawgaussfixed(x, scale, 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)**2/(2*(sigma**2))))

def fitgaussfixed(x, pardict, args):
    sigma = pardict['sigma'].value
    scale = pardict['scale'].value
    return rawgaussfixed(x, scale, sigma)


fitobjfixed1 = FitObject(fitgaussfixed, [parscale1, parsigma1])
fitobjfixed2 = FitObject(fitgaussfixed, [parscale2, parsigma2])
fitobjfixed3 = FitObject(fitgaussfixed, [parscale3, parsigma3])

fitdata1 = fitobjfixed1.fit(x000a, y000a)
fitdata2 = fitobjfixed2.fit(x000b, y000b)
fitdata3 = fitobjfixed3.fit(x000c, y000c)

time=1000




def M2func(tscale, pardict, args):
    A = pardict['A'].value
    M2 = pardict['M2'].value
    M4 = pardict['M4'].value
    return A*(1-((tscale**2.0)/2.0)*M2 + ((tscale**4.0)/24.0)*M4 )

PA = 6000.
PM2 = 1.7E10
PM4 = 9E19

parA1 = fit_parameter('A', PA)
parM21 = fit_parameter('M2',PM2, displayname='M_2')
parM41 = fit_parameter('M4',PM4, displayname='M_4')

parA2 = fit_parameter('A', PA)
parM22 = fit_parameter('M2',PM2, displayname='M_2')
parM42 = fit_parameter('M4',PM4, displayname='M_4')

parA3 = fit_parameter('A', PA)
parM23 = fit_parameter('M2',PM2, displayname= 'M_2')
parM43 = fit_parameter('M4',PM4, displayname='M_4')

x0001=af
y0001=bf
x0002=af
y0002=cf
x0003=af
y0003=df

ulA1=0.000015 #pol
ulB1=0.000030
ulA2=0.000020
ulB2=0.000040
ulA3=0.000100
ulB3=0.000300

ll=0.000002

x000aa = x0001[x0001<=ulA1]
x000a = x000aa[x000aa>=ll] 
x000bb = x0001[x0001<=ulB1]
x000b = x000bb[x000bb>=ll] 
y000aa = y0001[x0001<=ulA1]
y000a = y000aa[x000aa>=ll]
y000bb = y0001[x0001<=ulB1]
y000b = y000bb[x000bb>=ll]

x000aa2 = x0002[x0002<=ulA2]
x000a2 = x000aa2[x000aa2>=ll] 
x000bb2 = x0002[x0002<=ulB2]
x000b2 = x000bb2[x000bb2>=ll] 
y000aa2 = y0002[x0002<=ulA2]
y000a2 = y000aa2[x000aa2>=ll]
y000bb2 = y0002[x0002<=ulB2]
y000b2 = y000bb2[x000bb2>=ll]

x000aa3 = x0003[x0003<=ulA3]
x000a3 = x000aa3[x000aa3>=ll] 
x000bb3 = x0003[x0003<=ulB3]
x000b3 = x000bb3[x000bb3>=ll] 
y000aa3 = y0003[x0003<=ulA3]
y000a3 = y000aa3[x000aa3>=ll]
y000bb3 = y0003[x0003<=ulB3]
y000b3 = y000bb3[x000bb3>=ll]

fitobj1 = FitObject(M2func, [parA1, parM21, parM41])
fitobj2 = FitObject(M2func, [parA2, parM22, parM42])
fitobj3 = FitObject(M2func, [parA3, parM23, parM43])

fitdata1 = fitobj1.fit(x000a, y000a)
fitdata2 = fitobj1.fit(x000b, y000b)

fitdata3 = fitobj2.fit(x000a2, y000a2)
fitdata4 = fitobj2.fit(x000b2, y000b2)

fitdata5 = fitobj3.fit(x000a3, y000a3)
fitdata6 = fitobj3.fit(x000b3, y000b3)

ax15.scatter(af*time, bf, label='Data', color = 'black', marker = ".", linewidth = 0.05,)
ax25.scatter(af*time, cf, label='Data', color = 'black', marker = ".", linewidth = 0.05,)
ax35.scatter(af*time, df, label='Data', color = 'black', marker = ".", linewidth = 0.0005,alpha=0.1)

ax15.plot(x000b*time, fitdata2, 'b-', label='1', alpha = 0.6, linewidth = 3)
ax15.plot(x000a*time, fitdata1, 'r-', label='0', alpha = 0.6, linewidth = 3)

ax25.plot(x000b2*time, fitdata4, 'b-', label='1', alpha = 0.6, linewidth = 3)
ax25.plot(x000a2*time, fitdata3, 'r-', label='0', alpha = 0.6, linewidth = 3)

ax35.plot(x000b3*time, fitdata6, 'b-', label='1', alpha = 0.8, linewidth = 3)
ax35.plot(x000a3*time, fitdata5, 'r-', label='0', alpha = 0.8, linewidth = 3)

ax15.set_xlim(0,0.00015*time)
ax25.set_xlim(0,0.00015*time)
ax35.set_xlim(0,0.0015*time)

#ax35.set_xlabel('Time / ms ', fontsize=fs)
"""
frame15 = plt.gca()
frame15.axes.yaxis.set_ticklabels([])

frame25 = plt.gca()
frame25.axes.yaxis.set_ticklabels([])

frame35 = plt.gca()
frame35.axes.yaxis.set_ticklabels([])
"""

ax15.axes.yaxis.set_visible(False)
ax25.axes.yaxis.set_visible(False)
ax35.axes.yaxis.set_visible(False)
"""
ax15.text(0.67, 0.8, u'$-$52 \u00B0C',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax15.transAxes,
        color='black', fontsize=fs)
ax25.text(0.7, 0.8, u'138 \u00B0C',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax25.transAxes,
        color='black', fontsize=fs)
ax35.text(0.7, 0.8, u'225 \u00B0C',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax35.transAxes,
        color='black', fontsize=fs)
"""
ax111.text(0.28, 1.08, 'Gaussian',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax111.transAxes,
        color='black', fontsize=fs)
ax1.text(0.29, 1.08, u'Integral',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax1.transAxes,
        color='black', fontsize=fs)
ax15.text(0.38, 1.08, u'M2M4',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax15.transAxes,
        color='black', fontsize=fs)


ax111.text(-0.15, 0.1, u'Low T: 221 K',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax111.transAxes,rotation='vertical',
        color='black', fontsize=fs)

ax222.text(-0.15, 0.1, u'Mid T: 411 K',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax222.transAxes,rotation='vertical',
        color='black', fontsize=fs)

ax333.text(-0.15, 0.07, u'High T: 498 K',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax333.transAxes,rotation='vertical',
        color='black', fontsize=fs)

ax111.text(-0.3, 0.16, u'Diamantane',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax111.transAxes,rotation='vertical',
        color='black', fontsize=fs)

ax222.text(-0.3, 0.16, u'Diamantane',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax222.transAxes,rotation='vertical',
        color='black', fontsize=fs)

ax333.text(-0.3, 0.16, u'Diamantane',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax333.transAxes,rotation='vertical',
        color='black', fontsize=fs)

ax15.tick_params(axis='both', which='major', labelsize=fs)
ax25.tick_params(axis='both', which='major', labelsize=fs)
ax35.tick_params(axis='both', which='major', labelsize=fs)

#plt.tight_layout()
#fig, ((ax111,ax1,ax15) ,(ax222,ax2,ax25)) = plt.subplots(nrows=2,ncols=3, figsize = (10,5.45))





######TRI


fs=14


ax5t = ax1t.twinx()
ax6t = ax2t.twinx()
ax7t = ax3t.twinx()


ax9t = ax1t.twinx()
ax10t = ax2t.twinx()
ax11t = ax3t.twinx()


ax9t.yaxis.tick_left()
ax10t.yaxis.tick_left()
ax11t.yaxis.tick_left()




# Initialization parameters
xdct=ppmt*400.14/1000  #convert ppm to kHz

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

parscale2 = fit_parameter('scale2', 40000.)
parX2 = fit_parameter('X2', 1., displayname='X / kHz')
parsigma2 = fit_parameter('sigma2', 10., displayname='sigma / kHz')

parscale3 = fit_parameter('scale', 40000.)
parX3 = fit_parameter('X', 0., displayname='X / kHz')
parsigma3 = fit_parameter('sigma', 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(fitgauss2, [parscale1, parX1, parsigma1,parscale2, parX2, parsigma2])
fitobj2 = FitObject(fitgauss2, [parscale1, parX1, parsigma1,parscale2, parX2, parsigma2])
fitobj3 = FitObject(fitgauss1, [parscale3, parX3, parsigma3])
#fitobj = FitObject(fitgauss2, [parscale1, parX1, parsigma1,parscale2, parX2, parsigma2, parc])

r1t=bt
r2t=ct
r3t=dt

fitdata1t = fitobj1.fit(xdct, r1t)
fitdata2t = fitobj2.fit(xdct, r2t)
fitdata3t = fitobj3.fit(xdct, r3t)

# Plot the results
# Data
ax1t.plot(xdct, r1t, 'k', label='Data')
ax2t.plot(xdct, r2t, 'k', label='Data')
ax3t.plot(xdct, r3t, 'k', label='Data')

ax5t.plot(-x11t+0.8, y11t, 'b', label='Data')
ax6t.plot(-x22t+0.8, y22t, 'b', label='Data')
ax7t.plot(-x33t+1.8, y33t, 'b', label='Data')

ax9t.plot(-x11t+0.8, i11t, 'g', label='Data')
ax10t.plot(-x22t+0.8, i22t, 'g', label='Data')
ax11t.plot(-x33t+1.8, i33t, 'g', label='Data')


r1t=bt
r2t=ct
r3t=dt

fitdata1t = fitobj1.fit(xdct, r1t)
fitdata2t = fitobj2.fit(xdct, r2t)
fitdata3t = fitobj3.fit(xdct, r3t)

# Plot the results
# Data
ax111t.plot(xdct, r1t, 'k', label='Data')
ax222t.plot(xdct, r2t, 'k', label='Data')
ax333t.plot(xdct, r3t, 'k', label='Data')


#plt.plot(xd, b1, 'g', label='Data im')
# Fitted function
ax111t.plot(xdct, fitdata1t, 'r-', label='Fit')
ax222t.plot(xdct, fitdata2t, 'r-', label='Fit')
ax333t.plot(xdct, fitdata3t, 'r-', label='Fit')

ax111t.plot(xdc, rawgauss (xdc,47776.222767,0.687384,25.106533), 'b', label='Broad signal', alpha = 0.5)
ax111t.plot(xdc, rawgauss (xdc,15812.791831,1.330435,2.804505), 'b', label='Narrow signal',alpha = 0.5)
ax222t.plot(xdc, rawgauss (xdc,197778.961186,0.729471,10.871957), 'b', label='Broad signal', alpha = 0.5)
ax222t.plot(xdc, rawgauss (xdc,120057.350475,0.754236,0.686342), 'b', label='Narrow signal',alpha = 0.5)

#plt.plot(xdc, rawgauss(xdc, 1100772, 1.5, 8.0))
limp1=165
limp2=60
limp3=8
fs=14
ax111t.set_xlim(-limp1,limp1)
ax222t.set_xlim(-limp2,limp2)
ax333t.set_xlim(-7,9)
#plt.xlim(-200,-100)
#plt.ylim(-5000,5000)
#plt.legend(frameon=False, prop={'size':12} )


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

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

ax111t.axes.yaxis.set_visible(False)
ax222t.axes.yaxis.set_visible(False)
ax333t.axes.yaxis.set_visible(False)

ax111t.tick_params(axis='both', which='major', labelsize=fs)
ax222t.tick_params(axis='both', which='major', labelsize=fs)
ax333t.tick_params(axis='both', which='major', labelsize=fs)
#plt.savefig('Gaussian_fit_P3.png', dpi = 300)
#plt.plot((fit[0]+fit[2]/1.1774, fit[0]-fit[2]/1.1774), (fit[1]/2, fit[1]/2), 'k-')
#plt.axvline(x=0)


#plt.plot(xdc, rawgauss(xdc, scale1, x1, sigma1), 'm-', label='Fit')
#plt.plot(xdc, rawgauss(xdc, parscale1.value, parX1.value, parsigma1.value)+parc.value, 'k-', label='Fit')
#plt.plot(xdc, rawgauss(xdc, parscale2.value, parX2.value, parsigma2.value)+parc.value, 'k-', label='Fit')

#plt.plot(xdc, r1-parc.value-rawgauss(xdc, parscale1.value, parX1.value, parsigma1.value), 'g', label='Data')
#plt.plot(xdc, r1-parc.value-rawgauss(xdc, parscale2.value, parX2.value, parsigma2.value), 'g', label='Data')

#plt.plot((0,0), (-20,0), 'r-')

#plt.plot((-60,60), (0,0), 'k-')
"""
ax111.text(0.05, 0.8, u'$-$52 \u00B0C',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax1.transAxes,
        color='black', fontsize=fs)
ax222.text(0.05, 0.8, u'138 \u00B0C',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax2.transAxes,
        color='black', fontsize=fs)
"""
#to here

ax5t.set_ylabel('M2 / kHz2', fontsize=fs)
ax6t.set_ylabel('M2 / kHz2', fontsize=fs)
ax7t.set_ylabel('M2 / kHz2', fontsize=fs)


ax9t.set_ylabel('Integral / kHz ', fontsize=fs)
ax10t.set_ylabel('Integral / kHz ', fontsize=fs)
ax11t.set_ylabel('Integral / kHz ', fontsize=fs)

ax9t.yaxis.set_label_position("left")
ax10t.yaxis.set_label_position("left")
ax11t.yaxis.set_label_position("left")


ax5t.tick_params(axis='y', which='major', labelsize=fs)#
ax6t.tick_params(axis='y', which='major', labelsize=fs)
ax7t.tick_params(axis='y', which='major', labelsize=fs)
ax1t.tick_params(axis='x', which='major', labelsize=fs)#
ax2t.tick_params(axis='x', which='major', labelsize=fs)
ax3t.tick_params(axis='x', which='major', labelsize=fs)
ax9t.tick_params(axis='y', which='major', labelsize=fs)#
ax10t.tick_params(axis='y', which='major', labelsize=fs)
ax11t.tick_params(axis='y', which='major', labelsize=fs)
#ax11.set_yticks(0,1,0.5)

ax5t.yaxis.label.set_color('blue')
ax6t.yaxis.label.set_color('blue')
ax7t.yaxis.label.set_color('blue')

ax9t.yaxis.label.set_color('green')
ax10t.yaxis.label.set_color('green')
ax11t.yaxis.label.set_color('green')

ax5t.tick_params(axis='y', colors='blue')
ax6t.tick_params(axis='y', colors='blue')
ax7t.tick_params(axis='y', colors='blue')
ax9t.tick_params(axis='y', colors='green')
ax10t.tick_params(axis='y', colors='green')
ax11t.tick_params(axis='y', colors='green')

#plt.plot(xd, b1, 'g', label='Data im')
# Fitted function
#ax1.plot(xdc, fitdata1, 'r-', label='Fit')
#ax2.plot(xdc, fitdata2, 'r-', label='Fit')
#ax3.plot(xdc, fitdata3, 'r-', label='Fit')

#plt.plot(xdc, rawgauss(xdc, 1100772, 1.5, 8.0))
limp1=165
limp2=60
limp3=30
fs=14
ax1t.set_xlim(-limp1,limp1)
ax2t.set_xlim(-limp2,limp2)
ax3t.set_xlim(-limp3,limp3)

ppp = -0.05
qqq = 1.1
ax9t.set_ylim(ppp,qqq)
ax10t.set_ylim(ppp,qqq)
ax11t.set_ylim(ppp,qqq)
aaa=750
bbb=170
ccc=14
ax5t.set_ylim(-aaa/22,aaa)
ax6t.set_ylim(-bbb/22,bbb)
ax7t.set_ylim(-ccc/22,ccc)

ddd=70000
eee=370000
fff=1900000

ax1t.set_ylim(-ddd/22,ddd)
ax2t.set_ylim(-eee/22,eee)
ax3t.set_ylim(-fff/22,fff)


#plt.xlim(-200,-100)
#plt.ylim(-5000,5000)
#plt.legend(frameon=False, prop={'size':12} )

ax1t.invert_xaxis()
ax2t.invert_xaxis()
ax3t.invert_xaxis()

ax111t.invert_xaxis()
ax222t.invert_xaxis()
ax333t.invert_xaxis()

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

#frame1 = plt.gca()
#frame1.axes.yaxis.set_ticklabels([])

#frame2 = plt.gca()
#frame2.axes.yaxis.set_ticklabels([])

#frame3 = plt.gca()
#frame3.axes.yaxis.set_ticklabels([])

ax1t.axes.yaxis.set_visible(False)
ax2t.axes.yaxis.set_visible(False)
ax3t.axes.yaxis.set_visible(False)

ax1t.tick_params(axis='both', which='major', labelsize=fs)
ax2t.tick_params(axis='both', which='major', labelsize=fs)
ax3t.tick_params(axis='both', which='major', labelsize=fs)




#plt.savefig('Gaussian_fit_P3.png', dpi = 300)
#plt.plot((fit[0]+fit[2]/1.1774, fit[0]-fit[2]/1.1774), (fit[1]/2, fit[1]/2), 'k-')
#plt.axvline(x=0)


#plt.plot(xdc, rawgauss(xdc, scale1, x1, sigma1), 'm-', label='Fit')
#plt.plot(xdc, rawgauss(xdc, parscale1.value, parX1.value, parsigma1.value)+parc.value, 'k-', label='Fit')
#plt.plot(xdc, rawgauss(xdc, parscale2.value, parX2.value, parsigma2.value)+parc.value, 'k-', label='Fit')

#plt.plot(xdc, r1-parc.value-rawgauss(xdc, parscale1.value, parX1.value, parsigma1.value), 'g', label='Data')
#plt.plot(xdc, r1-parc.value-rawgauss(xdc, parscale2.value, parX2.value, parsigma2.value), 'g', label='Data')

#plt.plot((0,0), (-20,0), 'r-')

#plt.plot((-60,60), (0,0), 'k-')
"""
ax1.text(0.05, 0.8, u'$-$81 \u00B0C',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax1.transAxes,
        color='black', fontsize=fs)
ax2.text(0.05, 0.8, u'87 \u00B0C',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax2.transAxes,
        color='black', fontsize=fs)


ax111.text(0.05, 0.8, u'$-$81 \u00B0C',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax111.transAxes,
        color='black', fontsize=fs)
ax222.text(0.05, 0.8, u'87 \u00B0C',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax222.transAxes,
        color='black', fontsize=fs)
"""
'''
scale=61966.917
X=1.5972
sigma=11.235761*(2*np.log(2))**0.5
c = 367.223
smc=scale-c
plt.axvline(x=X, color = 'k', alpha = 0.1)
plt.plot((X,sigma), (scale+c,scale+c), 'k-')
plt.plot((X,60), (c,c), 'k-')
plt.plot((X,X+sigma), (scale/2+c,scale/2+c), 'k-')
plt.plot((+X,+X-sigma), (scale/2+c+4000,scale/2+c+4000), 'g-')
plt.plot((0,0), (0+c,scale/2+c), 'r-')
plt.plot((15,15), (c+31000,scale/2+c+31000), 'm-')
'''
 #or

x0001=aft
y0001=bft
x0002=aft
y0002=cft
x0003=aft
y0003=dft

ulA=0.00010
ulB=0.001
ulC=0.01

ll=0.000001

#t0s = t0[t0<=ul]


#FID0s = t0[t0<=ul]


x000aa = x0001[x0001<=ulA]
x000a = x000aa[x000aa>=ll] 

x000bb = x0002[x0002<=ulB]
x000b = x000bb[x000bb>=ll] 

x000cc = x0003[x0003<=ulC]
x000c = x000cc[x000cc>=ll] 

y000aa = y0001[x0001<=ulA]
y000a = y000aa[x000aa>=ll]
y000bb = y0002[x0002<=ulB]
y000b = y000bb[x000bb>=ll]
y000cc = y0003[x0003<=ulC]
y000c = y000cc[x000cc>=ll]
parscale1 = fit_parameter('scale', 4000.)
parX1 = fit_parameter('X', 0., displayname='X')
parsigma1 = fit_parameter('sigma', 0.000002, displayname='sigma')

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

parscale3 = fit_parameter('scale', 30000.)
parX3 = fit_parameter('X', 0., displayname='X / kHz')
parsigma3 = fit_parameter('sigma', 0.0005, displayname='sigma / kHz')


fitobjfixed1 = FitObject(fitgaussfixed, [parscale1, parsigma1])
fitobjfixed2 = FitObject(fitgaussfixed, [parscale2, parsigma2])
fitobjfixed3 = FitObject(fitgaussfixed, [parscale3, parsigma3])

fitdata1 = fitobjfixed1.fit(x000a, y000a)
fitdata2 = fitobjfixed2.fit(x000b, y000b)
fitdata3 = fitobjfixed3.fit(x000c, y000c)

time=1000




def M2func(tscale, pardict, args):
    A = pardict['A'].value
    M2 = pardict['M2'].value
    M4 = pardict['M4'].value
    return A*(1-((tscale**2.0)/2.0)*M2 + ((tscale**4.0)/24.0)*M4 )

PA = 6000.
PM2 = 1.7E10
PM4 = 9E19

parA1 = fit_parameter('A', PA)
parM21 = fit_parameter('M2',PM2, displayname='M_2')
parM41 = fit_parameter('M4',PM4, displayname='M_4')

parA2 = fit_parameter('A', PA)
parM22 = fit_parameter('M2',PM2, displayname='M_2')
parM42 = fit_parameter('M4',PM4, displayname='M_4')

parA3 = fit_parameter('A', PA)
parM23 = fit_parameter('M2',PM2, displayname= 'M_2')
parM43 = fit_parameter('M4',PM4, displayname='M_4')

x0001=aft
y0001=bft
x0002=aft
y0002=cft
x0003=aft
y0003=dft

ulA1=0.000008
ulB1=0.000012
ulA2=0.000018 #pol
ulB2=0.000028
ulA3=0.000050
ulB3=0.000200

ll=0.000001

x000aa = x0001[x0001<=ulA1]
x000a = x000aa[x000aa>=ll] 
x000bb = x0001[x0001<=ulB1]
x000b = x000bb[x000bb>=ll] 
y000aa = y0001[x0001<=ulA1]
y000a = y000aa[x000aa>=ll]
y000bb = y0001[x0001<=ulB1]
y000b = y000bb[x000bb>=ll]

x000aa2 = x0002[x0002<=ulA2]
x000a2 = x000aa2[x000aa2>=ll] 
x000bb2 = x0002[x0002<=ulB2]
x000b2 = x000bb2[x000bb2>=ll] 
y000aa2 = y0002[x0002<=ulA2]
y000a2 = y000aa2[x000aa2>=ll]
y000bb2 = y0002[x0002<=ulB2]
y000b2 = y000bb2[x000bb2>=ll]

x000aa3 = x0003[x0003<=ulA3]
x000a3 = x000aa3[x000aa3>=ll] 
x000bb3 = x0003[x0003<=ulB3]
x000b3 = x000bb3[x000bb3>=ll] 
y000aa3 = y0003[x0003<=ulA3]
y000a3 = y000aa3[x000aa3>=ll]
y000bb3 = y0003[x0003<=ulB3]
y000b3 = y000bb3[x000bb3>=ll]

fitobj1 = FitObject(M2func, [parA1, parM21, parM41])
fitobj2 = FitObject(M2func, [parA2, parM22, parM42])
fitobj3 = FitObject(M2func, [parA3, parM23, parM43])

fitdata1 = fitobj1.fit(x000a, y000a)
fitdata2 = fitobj1.fit(x000b, y000b)


fitdata3 = fitobj2.fit(x000a2, y000a2)
fitdata4 = fitobj2.fit(x000b2, y000b2)

fitdata5 = fitobj3.fit(x000a3, y000a3)
fitdata6 = fitobj3.fit(x000b3, y000b3)

ax15t.scatter(aft*time, bft, label='Data', color = 'black', marker = ".", linewidth = 0.05,)
ax25t.scatter(aft*time, cft, label='Data', color = 'black', marker = ".", linewidth = 0.05,)
ax35t.scatter(aft*time, dft, label='Data', color = 'black', marker = ".", linewidth = 0.05,)

ax15t.plot(x000b*time, fitdata2, 'b-', label='1', alpha = 0.6, linewidth = 3)
ax15t.plot(x000a*time, fitdata1, 'r-', label='0', alpha = 0.6, linewidth = 3)

ax25t.plot(x000b2*time, fitdata4, 'b-', label='1', alpha = 0.6, linewidth = 3)
ax25t.plot(x000a2*time, fitdata3, 'r-', label='0', alpha = 0.6, linewidth = 3)

ax35t.plot(x000b3*time, fitdata6, 'b-', label='1', alpha = 0.8, linewidth = 3)
ax35t.plot(x000a3*time, fitdata5, 'r-', label='0', alpha = 0.8, linewidth = 3)

ax15t.set_xlim(0,0.00007*time)
ax25t.set_xlim(0,0.00012*time)
ax35t.set_xlim(0,0.0015*time)

ax35t.set_xlabel('Time / ms ', fontsize=fs)
"""
frame15 = plt.gca()
frame15.axes.yaxis.set_ticklabels([])

frame25 = plt.gca()
frame25.axes.yaxis.set_ticklabels([])

frame35 = plt.gca()
frame35.axes.yaxis.set_ticklabels([])
"""
ax15t.axes.yaxis.set_visible(False)
ax25t.axes.yaxis.set_visible(False)
ax35t.axes.yaxis.set_visible(False)
"""
ax15.text(0.67, 0.8, u'$-$81 \u00B0C',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax15.transAxes,
        color='black', fontsize=fs)
ax25.text(0.7, 0.8, u'87 \u00B0C',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax25.transAxes,
        color='black', fontsize=fs)
"""
#ax111t.text(0.22, 1.08, 'Spec, Gaussian',
#        verticalalignment='bottom', horizontalalignment='left',
#        transform=ax111t.transAxes,
#        color='black', fontsize=fs)
#ax1t.text(0.25, 1.08, u'Spec, integral',
#        verticalalignment='bottom', horizontalalignment='left',
#        transform=ax1t.transAxes,
#        color='black', fontsize=fs)
#ax1111t.text(0.25, 1.08, u'FID, Gaussian',
#        verticalalignment='bottom', horizontalalignment='left',
#        transform=ax1111t.transAxes,
#        color='black', fontsize=fs)
#ax15t.text(0.3, 1.08, u'FID, $M_2M_4$',
#        verticalalignment='bottom', horizontalalignment='left',
#        transform=ax15t.transAxes,
 #       color='black', fontsize=fs)


ax111t.text(-0.15, 0.1, u'Low T: 192 K',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax111t.transAxes,rotation='vertical',
        color='black', fontsize=fs)

ax222t.text(-0.15, 0.1, u'Mid T: 360 K',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax222t.transAxes,rotation='vertical',
        color='black', fontsize=fs)

ax333t.text(-0.15, 0.07, u'High T: 469 K',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax333t.transAxes,rotation='vertical',
        color='black', fontsize=fs)

ax111t.text(-0.3, 0.13, u'Triamantane',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax111t.transAxes,rotation='vertical',
        color='black', fontsize=fs)

ax222t.text(-0.3, 0.13, u'Triamantane',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax222t.transAxes,rotation='vertical',
        color='black', fontsize=fs)

ax333t.text(-0.3, 0.13, u'Triamantane',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax333t.transAxes,rotation='vertical',
        color='black', fontsize=fs)

ax15t.tick_params(axis='both', which='major', labelsize=fs)
ax25t.tick_params(axis='both', which='major', labelsize=fs)
ax35t.tick_params(axis='both', which='major', labelsize=fs)

#plt.tight_layout()


plt.savefig('M2_Di_tri_SE_all.svg', bbox_inches='tight')
plt.show()