#間違っていても責任はとりません、野元 #Single parabric band modelを使ってLorentz数、Keを計算します。 # -*- coding: utf-8 -*- import os import numpy as np from numpy import sqrt import scipy.integrate as integrate import sys import csv import pandas as pd import openpyxl from matplotlib import pyplot as plt from matplotlib import gridspec from tklib.tkutils import terminate, pfloat, pint, getarg, getintarg, getfloatarg from tklib.tkapplication import tkApplication from tklib.tkvariousdata import tkVariousData from tklib.tktransport.tkTransport import LorentzNumber_FEA, FermiIntegral_x2 from tklib.tktransport.tkTransport import FermiIntegral_fast as FermiIntegral from tklib.tktransport.tkTransport import meff2NC_FEA, NC2meff_FEA, meff2DC0_FEA #=================================== # physical constants #=================================== pi = 3.14159265358979323846 e = 1.60218e-19 # C"; kB = 1.380658e-23 # JK-1"; kBe = kB / e #8.6173303e-5 # kB/e kpi05 = 2.0 / sqrt(pi) # Global variables infile = "20221219Ba3GeO.xlsx" outfile = None meff = 1.0 T_label = 0 S_label = 1 sigma_label = 2 kappa_label = 3 nmaxiter = 100 eps = 1.0e-10 a = 1000 #max b = -300 #min #=================================== # figure configuration #=================================== figsize = (12, 8) fontsize = 16 legend_fontsize = 12 markersize = 8.0 app = None #=================================== # Treat arguments #=================================== infile = getarg(1, defval = infile) T_label = getarg(2, defval = T_label) S_label = getarg(3, defval = S_label) sigma_label = getarg(4, defval = sigma_label) kappa_label = getarg(5, defval = kappa_label) meff = getfloatarg(6, defval = meff) header, ext = os.path.splitext(infile) filebody = os.path.basename(header) outfile = f'{header}-out.xlsx' #=================================== # Other functions #=================================== def usage(app): argv = sys.argv print("") print("Usage: python {} infile T_label S_label(S/cm) sigma_label(S/cm) kappa_label(W/m/K)".format(argv[0])) print(" ex: python {} single '{}' '{}' '{}' '{}' '{}'".format(argv[0], infile, T_label, S_label, sigma_label, kappa_label)) def savecsv(outfile, header, datalist): try: print("Write to [{}]".format(outfile)) f = open(outfile, 'w') except: # except IOError: print("line 98 Error: Can not write to [{}]".format(outfile)) else: fout = csv.writer(f, lineterminator='\n') fout.writerow(header) # fout.writerows(data) for i in range(0, len(datalist[0])): a = [] for j in range(len(datalist)): a.append(datalist[j][i]) fout.writerow(a) f.close() def read_file(fname): print("") datafile = tkVariousData(infile) labels, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = usage) return datafile, labels, datalist #Calculation method def cal_ke(sigma, L, T): return sigma * L * T def Fr(EF, r): result = integrate.quad(lambda E: E**r / (1 + np.exp(E - EF)), 0, np.inf) #print("felmi_integral = {}".format(result)) fel = float(result[0]) #print(fel) return fel def S_comparison(T, S, sigma, EF): diff = S - kBe * (2.0 * Fr(EF, r = 1.0) / Fr(EF, r = 0.0) - EF) #print("diff = {}".format(diff)) return diff def Lorenz(EF): answer = kBe**2 * (3.0 * Fr(EF, r = 2.0) / Fr(EF, r = 0.0) - (2.0 * Fr(EF, r = 1.0) / Fr(EF, r = 0.0))**2) return answer def cal_EF_from_S_bisection(T, S, sigma, max, min, err): xEF = (max + min) / 2.0 ymin = S_comparison(T, S, sigma, min) ymax = S_comparison(T, S, sigma, max) yh = S_comparison(T, S, sigma, xEF) print(f"xEF={xEF} ymin={ymin} ymax={ymax} yh={yh} min={min} max={max}") if ymin * ymax > 0.0: print("line Error: ymin * ymax > 0.0") app.terminate(pause = True) for i in range(nmaxiter): print(f"xEF={xEF:8.3f} d=({ymin:8.4g} {yh:8.4g} {ymax:8.4g}) min={min:8.3f} max={max:8.3f}") if abs(yh) < eps: # if abs(yh) < eps or abs(max - min) < eps: return xEF if ymin * yh > 0.0: min = xEF else: max = xEF xEF = (max + min) / 2.0 ymin = S_comparison(T, S, sigma, min) ymax = S_comparison(T, S, sigma, max) yh = S_comparison(T, S, sigma, xEF) else: print(f"Did not converge in {nmaxiter} iterations.") exit() def cal_Ne(T, xEF, NC): F05 = Fr(xEF, r = 0.5) N = NC * kpi05 * F05 return N def main(): global app app = tkApplication() logfile = app.replace_path(infile) print(f"Open logfile [{logfile}]") app.redirect(targets = ["stdout", logfile], mode = 'w') print("") print("Estimate Lorentz factor from Seebeck coefficient with r = 0.0") print("input file :", infile) print("output file:", outfile) print("T_label :", T_label) print("S_label :", S_label) print("sigma_label:", sigma_label) print("kappa_label:", kappa_label) print("") if '***' in T_label: app.terminate("Error: Choose T", usage = usage, pause = True) if '***' in S_label: app.terminate("Error: Choose S", usage = usage, pause = True) if '***' in sigma_label: app.terminate("Error: Choose sigma", usage = usage, pause = True) if '***' in kappa_label: app.terminate("Error: Choose kappa", usage = usage, pause = True) if 'S/m' in sigma_label or 'cm' not in sigma_label: app.terminate(f"Error: The unit of sigma must be 'S/cm'. The given label is [{sigma_label}]", usage = usage, pause = True) if 'uK' in S_label or 'microK' in S_label: app.terminate(f"Error: The unit of S must be 'V/K'. The given label is [{S_label}]", usage = usage, pause = True) print("Read data from [{}]".format(infile)) datafile, labels, datalist = read_file(infile) T_idx = datafile.FindLabelIndex(T_label, flag = 'i') S_idx = datafile.FindLabelIndex(S_label, flag = 'i') sigma_idx = datafile.FindLabelIndex(sigma_label, flag = 'i') kappa_idx = datafile.FindLabelIndex(kappa_label, flag = 'i') print("T index: ", T_idx) print("S index: ", S_idx) print("sigma index: ", sigma_idx) print("kappa index: ", kappa_idx) Tlabel, Tlist = datafile.FindDataArray(T_label, flag = 'i') Slabel, Slist = datafile.FindDataArray(S_label, flag = 'i') # V/K sigmalabel, sigmalist = datafile.FindDataArray(sigma_label, flag = 'i') # S/cm kappalabel, kappalist = datafile.FindDataArray(kappa_label, flag = 'i') print("") print(f"meff = {meff} me") print(f"{'T(K)':8} {'S(V/K)':10} {'sigma(S/cm)':10} {'kappa':10}") for i in range(len(Tlist)): print(f"{Tlist[i]:8.3g} {Slist[i]:10.4g} {sigmalist[i]:10.4g} {kappalist[i]:10.4g}") print("") print(f"meff = {meff} me") EFlist = [] NClist = [] Nelist = [] EFkTlist = [] Llist = [] kelist = [] for ic in range(len(Tlist)): T = Tlist[ic] S = Slist[ic] sigma = sigmalist[ic] kappa = kappalist[ic] print(f"{ic:3d}: T={T} K S={S} V/K sigma={sigma} S/cm kappa={kappa}") xEF = cal_EF_from_S_bisection(T, abs(S), sigma * 1.0e-2, max = a, min = b, err = 1e-10) EF = xEF * kB * T / e NC = meff2NC_FEA(meff, T) # cm-3 Ne = cal_Ne(T, xEF, NC) # cm-3 L = Lorenz(xEF) ke = cal_ke(sigma / 1.0e-2, L, T) # sigma to S/m print(f"T = {T} K") print(f"S = {S*1.0e6} uV/K") print(f"EF/kT = {xEF}") print(f"EF = {EF} eV") print(f"NC = {NC} cm^-3") print(f"Ne = {Ne} cm^-3") print(f"L = {L}") print(f"ke = {ke} W/m/K") print(f"kT/e = { kB * T / e}") EFlist.append(EF) EFkTlist.append(xEF) NClist.append(NC) Nelist.append(Ne) Llist.append(L) kelist.append(ke) df = pd.DataFrame(np.array([Tlist, Nelist, Slist, sigmalist, kappalist, EFkTlist, EFlist, Llist, kelist]).T, columns = ["T(K)", "Ne(cm^-3)", "S(V/K)", "simga(S/cm)", "kappa(W/m/K)", "EF/kBT", "EF(eV)", "L", "kappa_e(W/m/K)"]) print("") print("Save L data to [{}]".format(outfile)) df.to_excel(outfile, index = False, header = True) #============================= # グラフの表示 #============================= print("") fig = plt.figure(figsize = figsize) ax1 = fig.add_subplot(2, 2, 1) ax1b = ax1.twinx() ax2 = fig.add_subplot(2, 2, 3) ax3 = fig.add_subplot(2, 2, 2) ax3b = ax3.twinx() ax4 = fig.add_subplot(2, 2, 4) ax1.tick_params(labelsize = fontsize) ax2.tick_params(labelsize = fontsize) ax3.tick_params(labelsize = fontsize) ax4.tick_params(labelsize = fontsize) ins1 = ax1.plot(Tlist, Nelist, label = '$N_e$ (cm$^{-3}$)', linestyle = '', marker = 'o', markerfacecolor = 'black', markeredgecolor = 'black') ax1.set_xlabel("$T$ (K)", fontsize = fontsize) ax1.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize) ax1.set_yscale('log') S = [Slist[i] * 1.0e6 for i in range(len(Slist))] ins2 = ax1b.plot(Tlist, S, label = '$S$ ($\\mu$V/K)', linestyle = '', marker = 's', markerfacecolor = 'blue', markeredgecolor = 'blue') # ax1b.set_xlabel("$T$ (K)", fontsize = fontsize) ax1b.set_ylabel("$S$ ($\\mu$V/K)", fontsize = fontsize) ins = ins1 + ins2 ax1.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'upper center') #loc = 'best') ax2.plot(Tlist, Llist, label = '$L$', linestyle = '', marker = 's') ax2.set_xlabel("$T$ (K)", fontsize = fontsize) ax2.set_ylabel("$L$", fontsize = fontsize) ins1 = ax3.plot(Nelist, sigmalist, label = '$\\sigma$ (S/cm)', linestyle = '', marker = 'o', markerfacecolor = 'black', markeredgecolor = 'black') ax3.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize) ax3.set_ylabel("$\\sigma$ (S/cm)", fontsize = fontsize) ax3.set_xscale('log') ins2 = ax3b.plot(Nelist, S, label = '$S$ ($\\mu$V/K)', linestyle = '', marker = 's', markerfacecolor = 'blue', markeredgecolor = 'blue') ax3b.set_ylabel("$S$ ($\\mu$V/K)", fontsize = fontsize) ins = ins1 + ins2 ax3.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'upper center') #loc = 'best') ax4.plot(Nelist, Llist, label = '$L$', linestyle = '', marker = 's') ax4.set_xlabel("$N_e$ (K)", fontsize = fontsize) ax4.set_ylabel("$L$", fontsize = fontsize) plt.tight_layout() plt.pause(0.1) print() print("NOTE: r is limited to 0") print() app.terminate("", usage = usage, pause = True) if __name__ == '__main__': main()