implement new temprature compensation

This commit is contained in:
Sgr A* VMT
2024-01-28 21:52:57 +08:00
parent d43819eac8
commit 4ab46b3f9b
4 changed files with 442 additions and 346 deletions

View File

@@ -2,49 +2,37 @@
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import pandas as pd
import warnings
warnings.filterwarnings("ignore")
class TempModel: class TempModel:
def __init__(self, amfg, tcc, tcfl, tctl, fmin, fmin_temp): def __init__(self, a_a, a_b, b_a, b_b, fmin, fmin_temp):
self.amfg = amfg self.a_a=a_a
self.tcc = tcc self.a_b=a_b
self.tcfl = tcfl self.b_a=b_a
self.tctl = tctl self.b_b=b_b
self.fmin = fmin self.fmin = fmin
self.fmin_temp = fmin_temp self.fmin_temp = fmin_temp
def _tcf(self, f, df, dt, tctl): def compensate(self, freq, temp_source, temp_target):
tctl = self.tctl if tctl is None else tctl if self.a_a == None or self.a_b == None or self.b_a == None or self.b_b == None:
tc = self.tcc + self.tcfl * df + tctl * df * df return freq
return f + self.amfg * tc * dt * f A=4*(temp_source*self.a_a)**2+4*temp_source*self.a_a*self.b_a+self.b_a**2+4*self.a_a
B=8*temp_source**2*self.a_a*self.a_b+4*temp_source*(self.a_a*self.b_b+self.a_b*self.b_a)+2*self.b_a*self.b_b+4*self.a_b-4*(freq-model.fmin)*self.a_a
def compensate(self, freq, temp_source, temp_target, tctl=None): C=4*(temp_source*self.a_b)**2+4*temp_source*self.a_b*self.b_b+self.b_b**2-4*(freq-model.fmin)*self.a_b
dt = temp_target - temp_source if(B**2-4*A*C<0):
dfmin = self.fmin * self.amfg * self.tcc * \ param_c=freq-param_linear(freq-model.fmin,self.a_a,self.a_b)*temp_source**2-param_linear(freq-model.fmin,self.b_a,self.b_b)*temp_source
(temp_source - self.fmin_temp) return param_linear(freq-model.fmin,self.a_a,self.a_b)*temp_target**2+param_linear(freq-model.fmin,self.b_a,self.b_b)*temp_target+param_c
df = freq - (self.fmin + dfmin) ax=(np.sqrt(B**2-4*A*C)-B)/2/A
if dt < 0.: param_a=param_linear(ax,self.a_a,self.a_b)
f2 = self._tcf(freq, df, dt, tctl) param_b=param_linear(ax,self.b_a,self.b_b)
dfmin2 = self.fmin * self.amfg * self.tcc * \ return param_a*(temp_target+param_b/2/param_a)**2+ax+model.fmin
(temp_target - self.fmin_temp) #print(-param_linear(ax,self.b_a,self.b_b)/2/param_linear(ax,self.a_a,self.a_b))
df2 = f2 - (self.fmin + dfmin2) #param_c=freq-param_linear(freq-model.fmin,self.a_a,self.a_b)*temp_source**2-param_linear(freq-model.fmin,self.b_a,self.b_b)*temp_source
f3 = self._tcf(f2, df2, -dt, tctl) #return param_linear(freq-model.fmin,self.a_a,self.a_b)*temp_target**2+param_linear(freq-model.fmin,self.b_a,self.b_b)*temp_target+param_c
ferror = freq - f3
freq = freq + ferror
df = freq - (self.fmin + dfmin)
return self._tcf(freq, df, dt, tctl)
def line_fit(x,a,b,c): def line_fit(x,a,b,c):
return a*x**2+b*x+c return a*x**2+b*x+c
def fit(data,tcc,tcfl,tctl): def line0(x,a,c):
result=[] return a*x**2+c
model.tcc=tcc def line120(x,a,c):
model.tcfl=tcfl return a*x**2-240*a*x+c
model.tctl-tctl
for j in range(len(datas)):
for i in range(len(data)):
result.append(model.compensate(datas[j][3000],35,data[i]))
return result
def area_find(temp,freq): def area_find(temp,freq):
middle=int(len(temp)/100/2)*100 middle=int(len(temp)/100/2)*100
i=j=100 i=j=100
@@ -54,7 +42,7 @@ def area_find(temp,freq):
if(i_flag): if(i_flag):
i=i+100 i=i+100
if middle-i>=0: if middle-i>=0:
linear_params, params_covariance = curve_fit(line_fit, temp[middle-i:middle+j],freq[middle-i:middle+j],maxfev=100000,ftol=1e-10,xtol=1e-20) linear_params, params_covariance = curve_fit(line_fit, temp[middle-i:middle+j],freq[middle-i:middle+j],maxfev=100000,ftol=1e-10,xtol=1e-10)
minus=line_fit(temp[middle-i:middle+j],linear_params[0],linear_params[1],linear_params[2])-freq[middle-i:middle+j] minus=line_fit(temp[middle-i:middle+j],linear_params[0],linear_params[1],linear_params[2])-freq[middle-i:middle+j]
if np.sum(np.square(minus))/len(minus)>threshold: if np.sum(np.square(minus))/len(minus)>threshold:
i=i-100 i=i-100
@@ -62,73 +50,111 @@ def area_find(temp,freq):
if(j_flag): if(j_flag):
j=j+100 j=j+100
if middle+j<=len(freq): if middle+j<=len(freq):
linear_params, params_covariance = curve_fit(line_fit, temp[middle-i:middle+j],freq[middle-i:middle+j],maxfev=100000,ftol=1e-10,xtol=1e-20) linear_params, params_covariance = curve_fit(line_fit, temp[middle-i:middle+j],freq[middle-i:middle+j],maxfev=100000,ftol=1e-10,xtol=1e-10)
minus=line_fit(temp[middle-i:middle+j],linear_params[0],linear_params[1],linear_params[2])-freq[middle-i:middle+j] minus=line_fit(temp[middle-i:middle+j],linear_params[0],linear_params[1],linear_params[2])-freq[middle-i:middle+j]
if np.sum(np.square(minus))/len(minus)>threshold: if np.sum(np.square(minus))/len(minus)>threshold:
j=j-100 j=j-100
j_flag=False j_flag=False
linear_params, params_covariance = curve_fit(line_fit, temp[middle-i:middle+j],freq[middle-i:middle+j],maxfev=100000,ftol=1e-10,xtol=1e-20) linear_params, params_covariance = curve_fit(line_fit, temp[middle-i:middle+j],freq[middle-i:middle+j],maxfev=100000,ftol=1e-10,xtol=1e-10)
return linear_params return linear_params
def data_process(path): def data_process(path):
data=[] freq=[]
file_path = path # 替换为你的文件路径 temp=[]
with open(file_path, 'r') as file: with open(path, 'r') as file:
# 逐行读取文件内容 # 逐行读取文件内容
lines = file.readlines() lines = file.readlines()
# 遍历每行内容 # 遍历每行内容
for line in lines: for line in lines:
data.append(line.split(',')) data=line.split(',')
file.close() try:
full_data=pd.DataFrame(data[1:-1],columns=data[0]) freq.append(float(data[3]))
temp=np.array(full_data['temp']).astype(np.float32) temp.append(float(data[5]))
freq=np.array(full_data['freq']).astype(np.float32) except:pass
freq=freq[::100] dv=int(len(temp)/1000)
temp=temp[::100] if dv>1:
plt.plot(temp[10:],freq[10:]) freq=np.array(freq[::dv])
linear_params=area_find(temp,freq) temp=np.array(temp[::dv])
plt.plot(temp,line_fit(temp,linear_params[0],linear_params[1],linear_params[2])) plt.plot(temp[20:],freq[20:])
data0=line_fit(np.arange(5,80,0.01),linear_params[0],linear_params[1],linear_params[2]) #linear_params=area_find(temp[20:],freq[20:])
return data0 param_bounds=([0,-np.inf,-np.inf],[np.inf,np.inf,np.inf])
linear_params, params_covariance = curve_fit(line_fit, temp[20:],freq[20:],bounds=param_bounds,maxfev=100000,ftol=1e-10,xtol=1e-10)
try:
plt.title("Range:"+str(int(np.max(freq[20:])-np.min(freq[20:]))))
except:
pass
axis=-1*linear_params[1]/2/linear_params[0]
if(axis>120):
linear_params1, params_covariance = curve_fit(line120, temp[20:],freq[20:],bounds=([0,-np.inf],[np.inf,np.inf]),maxfev=100000,ftol=1e-10,xtol=1e-10)
plt.plot(temp[20:],line120(temp[20:],linear_params1[0],linear_params1[1]))
return [linear_params1[0],-240*linear_params1[0],line120(120,linear_params1[0],linear_params1[1])]
elif(axis<0):
linear_params1, params_covariance = curve_fit(line0, temp[20:],freq[20:],bounds=([0,-np.inf],[np.inf,np.inf]),maxfev=100000,ftol=1e-10,xtol=1e-10)
plt.plot(temp[20:],line0(temp[20:],linear_params1[0],linear_params1[1]))
return [linear_params1[0],0,line0(axis,linear_params1[0],linear_params1[1])]
plt.plot(temp[20:],line_fit(temp[20:],linear_params[0],linear_params[1],linear_params[2]))
linear_params[2]=line_fit(axis,linear_params[0],linear_params[1],linear_params[2])
return linear_params
def param_linear(x,a,b):
return a*x+b
while(1): while(1):
plt.figure(figsize=(25, 15)) plt.figure(figsize=(25, 15))
paths=['./data1','./data2','./data3','./data4'] paths=['./data1','./data2','./data3','./data4']
datas=[] a=[]
b=[]
freqs=[]
num=241 num=241
threshold=int(input('threshold set(recommend start from 250):\n请输入阈值设置(默认推荐250):\n')) #threshold=int(input('threshold set(recommend start from 1000):\n请输入阈值设置(默认推荐1000):\n'))
try: try:
for path in paths: for path in paths:
plt.subplot(num) plt.subplot(num)
num+=1 num+=1
datas.append(data_process(path)) temp=data_process(path)
a.append(temp[0])
b.append(temp[1])
freqs.append(temp[2])
except: except:
print("please make sure you have move the 4 data file to IDM folder\n请确认你有把4个文件拷到IDM文件夹内") print("please make sure you have move the 4 data file to IDM folder\n请确认你有把4个文件拷到IDM文件夹内")
break break
#反向求值 #反向求值
model=TempModel(1,-2.1429828e-05,-1.8980091e-10,3.6738370e-16,2943053.84,20.33) model=TempModel(None,None,None,None,2943053.8415908813,23.33)
p0=[-2.1429828e-05,-1.8980091e-10,3.6738370e-16] linear_params, params_covariance = curve_fit(param_linear, np.array(freqs)-model.fmin,a,maxfev=100000,ftol=1e-10,xtol=1e-10)
params, params_covariance = curve_fit(fit,np.arange(5,80,0.01),np.hstack(datas),p0=p0,maxfev=1000000,ftol=1e-10,xtol=1e-10) model.a_a=linear_params[0]
model.a_b=linear_params[1]
linear_params1, params_covariance = curve_fit(param_linear, np.array(freqs)-model.fmin,b,maxfev=100000,ftol=1e-10,xtol=1e-10)
model.b_a=linear_params1[0]
model.b_b=linear_params1[1]
for path in paths: for path in paths:
plt.subplot(num) plt.subplot(num)
num+=1 num+=1
data=[] freq=[]
file_path = path # 替换为你的文件路径 temp=[]
with open(file_path, 'r') as file: with open(path, 'r') as file:
# 逐行读取文件内容 # 逐行读取文件内容
lines = file.readlines() lines = file.readlines()
# 遍历每行内容 # 遍历每行内容
for line in lines: for line in lines:
data.append(line.split(',')) data=line.split(',')
file.close() try:
full_data=pd.DataFrame(data[1:-1],columns=data[0]) freq.append(float(data[3]))
temp=np.array(full_data['temp']).astype(np.float32) temp.append(float(data[5]))
freq=np.array(full_data['freq']).astype(np.float32) except:pass
freq=freq[::100] dv=int(len(temp)/10000)
temp=temp[::100] if dv>1:
freq=np.array(freq[::dv])
temp=np.array(temp[::dv])
temp=temp[200:]
freq=freq[200:]
result0=[] result0=[]
for i in range(len(temp)): for i in range(len(temp)):
result0.append(model.compensate(freq[i],temp[i],20.66)) result0.append(model.compensate(freq[i],temp[i],50))
plt.plot(temp[10:],result0[10:]) plt.plot(temp,result0)
plt.savefig('fit.png') try:
plt.title("Range:"+str(int(np.max(result0)-np.min(result0))))
except:
pass
plt.savefig('fit_output.png')
print('fit result:') print('fit result:')
print('tc_tcc:'+str(params[0])+'\ntc_tcfl:'+str(params[1])+'\ntc_tctl:'+str(params[2])) print('tc_a_a:'+str(model.a_a)+'\ntc_a_b:'+str(model.a_b)+'\ntc_b_a:'+str(model.b_a)+'\ntc_b_b:'+str(model.b_b))
break break

596
idm.py

File diff suppressed because it is too large Load Diff

View File

@@ -12,12 +12,13 @@ fi
# install idm requirements to env # install idm requirements to env
echo "idm: installing python requirements to env, this may take 10+ minutes." echo "idm: installing python requirements to env, this may take 10+ minutes."
sudo apt-get install g++ gfortran sudo apt-get install g++
sudo apt-get install gfortran
sudo apt-get install libopenblas-dev
"${KENV}/bin/pip" install -r "${BKDIR}/requirements.txt" "${KENV}/bin/pip" install -r "${BKDIR}/requirements.txt"
# update link to idm.py # update link to idm.py
echo "idm: linking klippy to idm.py." echo "idm: linking klippy to idm.py."
sudo apt-get install g++ gfortran libopenblas-dev
if [ -e "${KDIR}/klippy/extras/idm.py" ]; then if [ -e "${KDIR}/klippy/extras/idm.py" ]; then
rm "${KDIR}/klippy/extras/idm.py" rm "${KDIR}/klippy/extras/idm.py"
fi fi

View File

@@ -2,4 +2,3 @@
numpy>=1.16.6 numpy>=1.16.6
scipy>=1.10.0 scipy>=1.10.0
matplotlib>=3.7.0 matplotlib>=3.7.0
pandas>=1.4.2