Updated 3 years ago
#! /usr/bin/env python3
import numpy as np
#from sincfit import *
from scipy.optimize import brentq
import matplotlib.pyplot as plt
import matplotlib as mpl
# Load the Pandas libraries with alias 'pd'
import pandas as pd
# Read data from file 'filename.csv'
# (in the same directory that your python process is based)
# Control delimiters, rows, column names with read_csv (see later)
data = pd.read_csv("Kappa_115mT_LK.csv")
# Preview the first 5 lines of the loaded data
data.head()
exp = pd.read_csv("Kappa_exp.csv")
exp.head()
#! /usr/bin/env python3
import numpy as np
#from sincfit import *
from scipy.optimize import brentq,root
import matplotlib.pyplot as plt
import matplotlib as mpl
# Memoization from http://code.activestate.com/recipes/578231-probably-the-fastest-memoization-decorator-in-the-/
# globals:
kB = 0.0861733 # Boltzmann constant in meV/K
mu_B = 5.7883818012e-2 # meV/T
def memoize(f):
""" Memoization decorator for functions taking one or more arguments. """
class memodict(dict):
def __init__(self, f):
self.f = f
def __call__(self, *args):
return self[args]
def __missing__(self, key):
ret = self[key] = self.f(*key)
return ret
return memodict(f)
@memoize
def fivegammaT(T,gamma):
return 5*gamma*T/(36*np.pi)
@memoize
def kappaGisquared(T,gamma):
return fivegammaT(T,gamma)**(2/3)
@memoize
def delta(T,TMF,T0,gamma):
return 3/(2**(1/3))*kappaGisquared(T,gamma)*(T-TMF)/T0
@memoize
def yalpha(alpha):
if alpha==0:
return 1
else:
svals=np.linspace(-100,100,100000,endpoint=True)
sqas=np.lib.scimath.sqrt(alpha/(1+svals**2))
toint=1-1/sqas*np.arctan(sqas)
integ=np.trapz(toint,svals)
return np.real(np.lib.scimath.sqrt(1+alpha)-integ/np.pi)
@memoize
def yaprime(alpha):
dalpha = 1e-6
if alpha < dalpha:
return 1/6
return np.real(0.5*(yalpha(alpha+dalpha)-yalpha(alpha-dalpha))/dalpha)
@memoize
def yadiff(alpha):
if alpha==0:
return 1/6
else:
svals=np.linspace(-100,100,100000,endpoint=True)
ass=alpha/(1+svals**2)
# sqas=np.lib.scimath.sqrt(alpha/(1+svals**2))
sqas=np.lib.scimath.sqrt(ass)
toint=1/sqas*np.arctan(sqas)-1/(1+ass)
integ=np.trapz(toint,svals)
return np.real(1/(2*np.lib.scimath.sqrt(1+alpha))-integ/(2*np.pi*alpha))
def flux0(B,fu,J,Q):
#phi_0 = (chi_0 * f.u.)/(mu_b ( mu_0)) * B for Bexternal
chi0 = 0.274 #external
#flux0 = (fu/(4*np.pi*9.274*0.1))*chi0/(1+chi0)*B
#flux0 = (fu/(4*np.pi*9.274*0.1))*chi0*B
flux0 = (B*mu_B/fu)*1/(J*Q**2)
print('flux0: {}'.format(flux0))
return flux0
def magdens0(B,fu,r): # this one makes no sense
rphi0 = mu_B*B/fu
flux0 = rphi0/r
#print('r.phi0/r: {}'.format(flux0))
return flux0
def kappa2out(kappa2in,uguess,T,TMF,T0,Q,J,flux0):
k2perp = kappa2in[0]
k2para = kappa2in[1]
alpha = (k2perp-k2para)/(k2perp+k2para)
#print('alpha = {:.4f}'.format(alpha))
gamma = uguess*kB*Q/J
y = yalpha(alpha)
yprime = yadiff(alpha)
k12 = (k2perp+k2para)**(-1/2)
k32 = (k2perp+k2para)**(-3/2)
k2perpline = delta(T,TMF,T0,gamma) + 1/6*uguess*flux0**2 + 1/6*1/(np.sqrt(2)*np.pi)*gamma*T*(3/2*k12*y+(4*k2para-2*k2perp)*k32*yprime )-k2perp
k2paraline = delta(T,TMF,T0,gamma) + 3/6*uguess*flux0**2 + 1/6*1/(np.sqrt(2)*np.pi)*gamma*T*(2*k12*y + (2*k2para-6*k2perp)*k32*yprime )-k2para
ksquared = np.array((k2perpline,k2paraline))
return ksquared
def kappa0out(kappain,uguess,T,TMF,T0,Q,J):
kB = 0.0861733 # Boltzmann constant in meV/K
gamma = uguess*kB*Q/J # using ~u
k2res = delta(T,TMF,T0,gamma) + fivegammaT(gamma,T)*1/kappain - kappain**2
return k2res
TMF = 30.6 # K
Q = 0.039 # A-1
J = 2.8 # meV/A
fu = 24.018 # A3
ucalc = 0.25/(J*Q**2)
# print(ucalc)
#T0 = 1.1 # K
#uguess = 58.7 # from u = 0.25 meV/A^3
T0 = 1.3 # K
uguess = 75.1 # from u = 0.32 meV/A^3
#uguess = 5964 # ~u, Jonas
#uguess = 13419 # ~u, Jonas, high field
#uguess = 127 # ~u made up to fit?
B = 0.104
klow=0.0001/Q
khigh=0.2000/Q
# for B in (0.104,0.240,0.340,0.350):
phi0 = flux0(B,fu,J,Q)
if False:
T = 30.0
r = J*Q*Q*(1+delta(T,TMF,T0,gamma))
phiO = magdens0(B,fu,r)
kperp=0.010 # A-1
kpara=0.011 # A-1
kappa2in = np.array(((kperp/Q)**2,(kpara/Q)**2))
print(kappa2in)
print(kappa2out(kappa2in,uguess,T,TMF,T0,Q,J,phi0))
print('# B\tT\tkperp\tkpara')
Tstart=27.0
Tstop=32.0
Tstep=0.1
Tnum=int((Tstop-Tstart)/Tstep)
Tvals = np.linspace(Tstart,Tstop,Tnum)
kparaplot = np.zeros_like(Tvals)
kperpplot = np.zeros_like(Tvals)
i = 0
for T in Tvals:
#r = J*Q*Q*(1+delta(T,TMF,T0,gamma))
#gamma = uguess*kB*Q/J
#phi0 = magdens0(B,fu,r)
#print(B,T)
#zero B fit to guess initial k values
soln,r = brentq(kappa0out,klow,khigh,args=(uguess,T,TMF,T0,Q,J),full_output=True)
if(r.converged):
kperp=soln*Q
kpara=kperp
if(B==0):
print('{:.3f}\t{:.1f}\t{:.6f}\t{:.6f}'.format(B, T, kperp, kpara))
kperpplot[i]=kperp
kparaplot[i]=kpara
else:
kappa2in = np.array(((kperp/Q)**2,(kpara/Q)**2))
#print(kappa2in)
huh = root(kappa2out,kappa2in,args=(uguess,T,TMF,T0,Q,J,phi0),method='lm',options={'factor': 0.1})
# print(huh)
if(huh.success):
# print(huh.x)
kout=np.real(np.lib.scimath.sqrt((huh.x)))*Q
print('{:.3f}\t{:.1f}\t{:.6f}\t{:.6f}'.format(B, T, *kout))
kperpplot[i]=kout[0]
kparaplot[i]=kout[1]
else:
print('{:.3f}\t{:.1f}\tfail\tfail'.format(B,T))
i = i + 1
plt.figure('B = '+str(B))
plt.plot(Tvals,kparaplot,label='$\kappa_{||}$')
plt.plot(Tvals,kperpplot,label='$\kappa_\perp$')
plt.xlabel('T [ K ]')
plt.ylabel('$\kappa$ [ Å$^{-1}$ ]')
plt.legend()
plt.tight_layout()
plt.show()
flux0: 0.058852704195906796
# B T kperp kpara
0.104 27.0 0.001829 0.004657
0.104 27.1 0.001870 0.004824
0.104 27.2 0.001913 0.005000
0.104 27.3 0.001961 0.005187
0.104 27.4 0.002012 0.005385
0.104 27.5 0.002067 0.005595
0.104 27.6 0.002128 0.005817
0.104 27.7 0.002195 0.006052
0.104 27.8 0.002270 0.006302
0.104 27.9 0.002355 0.006566
0.104 28.0 0.002450 0.006845
0.104 28.1 0.002559 0.007141
0.104 28.2 0.002684 0.007453
0.104 28.3 0.002829 0.007782
0.104 28.4 0.002997 0.008129
0.104 28.5 0.003194 0.008494
0.104 28.6 0.003423 0.008876
0.104 28.7 0.003693 0.009277
0.104 28.8 0.004009 0.009696
0.104 28.9 0.004379 0.010135
0.104 29.0 0.004810 0.010595
0.104 29.1 0.005309 0.011079
0.104 29.2 0.005883 0.011591
0.104 29.3 0.006539 0.012137
0.104 29.4 0.007280 0.012725
0.104 29.6 0.008108 0.013362
0.104 29.7 0.009025 0.014057
0.104 29.8 0.010027 0.014815
0.104 29.9 0.011108 0.015641
0.104 30.0 0.012258 0.016534
0.104 30.1 0.013463 0.017489
0.104 30.2 0.014710 0.018498
0.104 30.3 0.015984 0.019551
0.104 30.4 0.017271 0.020634
0.104 30.5 0.018561 0.021739
0.104 30.6 0.019843 0.022853
0.104 30.7 0.021110 0.023970
0.104 30.8 0.022359 0.025083
0.104 30.9 0.023587 0.026188
0.104 31.0 0.024790 0.027281
0.104 31.1 0.025969 0.028360
0.104 31.2 0.027124 0.029423
0.104 31.3 0.028254 0.030470
0.104 31.4 0.029360 0.031500
0.104 31.5 0.030444 0.032514
0.104 31.6 0.031505 0.033511
0.104 31.7 0.032546 0.034492
0.104 31.8 0.033566 0.035457
0.104 31.9 0.034566 0.036407
0.104 32.0 0.035549 0.037342
plt.figure('~u = {}'.format(uguess))
plt.plot(Tvals, kparaplot,label='new Para')
plt.plot(Tvals, kperpplot,label='new Perp')
plt.plot(data['Temp Para'], data['Kappa Para'], label='Laura Para')
plt.plot(data['Temp Perp'], data['Kappa Perp'], label='Laura Perp')
plt.legend()
plt.xlabel('T [ K ]')
plt.ylabel('$\kappa$ [ Å$^{-1}$ ]')
plt.title('Comparing the fit of Laura with ours, B = 115mT')
plt.savefig('Laura_compare_115mT.png', dpi=300, bbox_inches='tight')
plt.show()