Jovian
Sign In
Learn practical skills, build real-world projects, and advance your career
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
Notebook Image
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()
Notebook Image