-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathtotcurrent combine.py
155 lines (122 loc) · 5.75 KB
/
totcurrent combine.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
# This is the script for plotting the spin
# texture and spin currents.
# Here we will add some noise to the order
# parameter
#
# Tian-Qi Chen, 20/11/2016
import numpy as np
from scipy.integrate import dblquad, quad
from scipy.stats import linregress
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 10, 8
# Parameter define
gamma = 0 # the strength of the noise
sigma_2 = 10000 # the sigma squared in the function of noise
sigma_1 = 10000 # the sigma squared in the exponential of all Chi's
sigma_3 = 10000 # the sigma squared in Chi_1 and Chi_{-1}
sigma_4 = 10000 # the sigma squared in Chi_0
C_1 = 1/np.sqrt(2*np.pi*sigma_3)
C_2 = 1/np.sqrt(2*np.pi*sigma_4)
x_0 = np.array([2, -2])
y_0 = np.array([0, 0]) # the coordinates of the centers of the noise
N = len(x_0) # the length of the coordinates array
#N_f = 6 # the winding number of f(x,y)
#N_g = 3 # the winding number of g(x,y)
#NN = (N_f**2-N_g**2)/N_g # the big N
L = 200 # the length of the density window
yl_i = 0 # the initial value of yl
yl_f = L/2 # the final value of yl
yl = np.linspace(yl_i, yl_f, 10) # the y position of placing single vortex when measuring
Nl = len(yl) # the length of yl
h_const = 1 # h_bar over m
N_rho = 100 # the N for the density function rho
delta = 0 # the length of the laser beam
# Define functions
def devx(x, y, t): # Derivative with respect to x
tmp_1 = 0
tmp_2 = 0 # initialization
rslt = 0
if t == 'f':
rslt = -y/(x**2+y**2)
elif t == 'g':
for i in range(N):
tmp_1 = tmp_1+np.exp(-((x-x_0[i])**2+(y-y_0[i])**2)/(2*sigma_2))
for j in range(N):
tmp_2 = tmp_2+(-1)*(x-x_0[j])/sigma_2 *np.exp(-((x-x_0[j])**2+(y-y_0[j])**2)/(2*sigma_2))
rslt = (-y/(x**2+y**2))*(1+gamma*tmp_1)+np.arctan(y/x)*gamma*tmp_2
return rslt
def rho(y, yk):
return np.sqrt((2*C_1**2+C_2**2)*np.exp(-(y+yk)**2/sigma_1))
def KdotF1(y, yk):
return 2*np.sqrt(2/3)*C_1*C_2*np.exp(-(y+yk)**2/(1*sigma_1))
def KdotF2(y, yk):
return 2/3 * (2*C_1**2+C_2**2) * np.exp(-(y+yk)**2/(1*sigma_1))
def integralResult(yl, N_f, N_g):
jc1, err_jc1 = dblquad(lambda x,y: rho(y, yl) * (devx(x, y, 'f') + devx(x, y, 'g') * KdotF1(y, yl)), -np.Inf, 0, lambda x: -np.Inf, lambda x: np.Inf)
jc2, err_jc2 = dblquad(lambda x,y: rho(y, yl) * (devx(x, y, 'f') + devx(x, y, 'g') * KdotF1(y, yl)), 0, np.Inf, lambda x: -np.Inf, lambda x: np.Inf)
j11, err_j11 = dblquad(lambda x,y: rho(y, yl) * (devx(x, y, 'f') * KdotF1(y, yl) + devx(x, y, 'g') * KdotF2(y, yl)), -np.Inf, 0, lambda x: -np.Inf, lambda x: np.Inf)
j12, err_j12 = dblquad(lambda x,y: rho(y, yl) * (devx(x, y, 'f') * KdotF1(y, yl) + devx(x, y, 'g') * KdotF2(y, yl)), 0, np.Inf, lambda x: -np.Inf, lambda x: np.Inf)
j21, err_j21 = dblquad(lambda x,y: rho(y, yl) * (devx(x, y, 'f') * KdotF2(y, yl) + devx(x, y, 'g') * KdotF2(y, yl)), -np.Inf, 0, lambda x: -np.Inf, lambda x: np.Inf)
j22, err_j22 = dblquad(lambda x,y: rho(y, yl) * (devx(x, y, 'f') * KdotF2(y, yl) + devx(x, y, 'g') * KdotF2(y, yl)), 0, np.Inf, lambda x: -np.Inf, lambda x: np.Inf)
nn = (N_f**2-N_g**2)/N_g
return (nn * (jc1+jc2) - (N_f * (j11+j12) - N_g * (j21+j22)))/np.pi
def integralAY(yl):
val1, abserr1 = quad(lambda y: rho(y, 0), -np.Inf, yl)
val2, abserr2 = quad(lambda y: rho(y, 0), yl, np.Inf)
return val1 - val2
# Plot the Jx-Ay relation
Jtot3 = [] # the total current measured
Jtot6 = []
Jtot9 = []
Ay = [] # the A_y
print('Start the integration')
for j in range(Nl):
#print('Start the integration')
Jtot3.append(integralResult(yl[j], 3, 2))
Jtot6.append(integralResult(yl[j], 4, 3))
Jtot9.append(integralResult(yl[j], 5, 4))
print('Integration completed for yl=', yl[j])
Ay.append(integralAY(yl[j]))
# Linear fitting
z3 = np.polyfit(Ay, Jtot3, 1)
z6 = np.polyfit(Ay, Jtot6, 1)
z9 = np.polyfit(Ay, Jtot9, 1)
# Plot the discrete data
plt.plot(Ay, Jtot3, 'ro')
plt.plot(Ay, Jtot6, 'go')
plt.plot(Ay, Jtot9, 'bo')
# Plot the fitted relation
xx = np.linspace(min(Ay), max(Ay), 50)
yy3 = z3[0]*xx+z3[1]
yy6 = z6[0]*xx+z6[1]
yy9 = z9[0]*xx+z9[1]
plt.plot(xx, yy3, 'r-')
plt.plot(xx, yy6, 'g-')
plt.plot(xx, yy9, 'b-')
#print(xx)
# Correlation coefficients calculation
#def rsquared(x, y):
# slope, intercept, r_value, p_value, std_err = linregress(x, y)
# return slope, r_value**2
#rsl1 = rsquared(Ay, Jtot)
#print('The correlation coefficient is: ', rsl1[1])
#print('The slope is: ', rsl1[0])
# Show the fits to data and fitted relation
#SIZE = 20
#plt.rc('axes', titlesize=SIZE) # fontsize of the axes title
#plt.rc('axes', labelsize=SIZE) # fontsize of the x and y labels
#plt.rc('xtick', labelsize=SIZE) # fontsize of the tick labels
#plt.rc('ytick', labelsize=SIZE) # fontsize of the tick labels
plt.xticks([0.0, 0.4, 0.8, 1.2], fontsize=20)
plt.yticks([1.0, 2.0, 3.0], fontsize=20)
plt.xlabel(r'$A_y$', fontsize=20)
plt.ylabel(r'$J_{tot}/G_{0}$', fontsize=20)
#plt.text(40, 150, 'N=3', color='blue')
#plt.text(60, 100, 'N=2', color='green')
#plt.text(60, 45, 'N=1', color='red')
#plt.title(r'$\sigma_1=$'+str(sigma_1)+', '+r'$\sigma_4=$'+str(sigma_4)+', '+r'$\sigma_3=$'+str(sigma_3))
#plt.grid(True)
#plt.figure(figsize=(60,60))
plt.savefig('JtotCombinedFractional.eps', format='eps')
plt.show()