前言:
我真的是拖延死了,差不多从1月开始,就开始这个项目,到2月中旬才结束,感谢老师的不杀之恩,哎,老师脾气真好,估计已经快触及老师的底线了,哎,小说真好看,抖音真快乐,微博真有趣,睡觉真享受,哎,改正改正,从头开始,好好做人,好好学习,抱住老师的大腿。一定要守时啊!!!!
刚刚老师分享了一个Taichi,编程使我快乐,我要去扒一扒了
- 两个快捷的东西
polyfit 与 curve_fit 不解释,看代码Sec.1与Sec.2
注意:curve_fit不一定最优,检验以后有空把包下载下来,修改成自己需要的模式,polyfit比curve_fit要好一些 - eta vs T与pd vs T的批处理
eta vs T存在的问题是没有标注Ttr与Tc 及图例的位置
pd vs T存在的问题是图例的位置
Sec.1
clc
clear all
M=dir('D:\20191226_论文\论文\cs\paper_help\data_liquid\*.xlsx');
hang=0;
for nn=1:length(M)
[num, text, raw] = xlsread(M(nn).name);
[rowN, columnN]=size(raw);
T=[];
xnum=0;
sigma=[];
Tc=num(rowN-1,1);
tt=0;
for n=1:rowN-1
if strcmp(num2str(num(n,12)),'NaN')||num(n,1)==Tc
tt=1;
else
sigma=[sigma,num(n,12)/1000];
T=[T,num(n,1)];
end
end
M(nn).name
Tr=T./Tc;
AAD0=1000000;
f0=20;
n0=20;
b=20;
x=1./(T);
y=log(sigma);
AA=polyfit(x,y,1);
A=AA(2);
B=AA(1);
ycal=AA(2)+AA(1).*x;
sigmacal=exp(ycal);
AAD1=sum((abs((sigmacal-sigma))./sigma))./length(sigma);
if AAD1<AAD0
AAD0=AAD1;
n0=exp(A);
b=B;
end
M(nn).name
AAD0
n0
b
data=[AAD0,n0,b];
c1=num2str(hang+1);
c2=num2str(hang+1);
c=['A',c1,':C',c2]
d=['A',c1,':K',c2]
xlswrite('D:\20191226_论文\论文\cs\paper_help\paper_new_2\new.xlsx',data , 'Sheet1',c)
xlswrite('D:\20191226_论文\论\cs\paper_help\paper_new_2\new.xlsx',M(nn).name , 'Sheet2',d)
hang=hang+1;
end
Sec.2
import math
import os
import xlrd
import matplotlib.pyplot as plt
import math
import xlrd
from xlutils.copy import copy
from xlwt import *
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
global hang
hang=0
def file_name(file_dir):
for root, dirs, files in os.walk(file_dir):
for f in files:
cal(f)
def func(x,a,b,c):
return b/(a+c*np.log(1-x))
def cal(fileName):
x=[]
y=[]
t=[]
fileName_q="E:\\cs\\20191226_论文\\论文\\cs\\paper_help\\data_liquid\\"+fileName
bk=xlrd.open_workbook(fileName_q)
shxrange=range(bk.nsheets)
try:
sh=bk.sheet_by_name("Sheet1")
except:
print ("wrong")
nrows=int(sh.nrows) #获取行数
global tc
tc=0.1
for i in range(1,nrows):
row_data=sh.row_values(i)
try:
if(str(sh.cell_value(i,0))!='undefined'):
global tc
tc=sh.cell_value(i,0)
except:
aaa=1
for i in range(1,nrows):
row_data=sh.row_values(i)
try:
global tc
if(str(sh.cell_value(i,11))!='undefined'and str(sh.cell_value(i,0))!='undefined'and sh.cell_value(i,0)!=tc):
xx=sh.cell_value(i,11)/1000
y.append(xx)
x.append(sh.cell_value(i,0)/tc)
except:
aaa=1
res_1=10000000000000000000000000000000000000000000
y_res=[]
aad=0
mad=0
y_flag=[]
flag=1
if(len(x)>=1):
try:
popt, pcov = curve_fit(func, x, y,maxfev=5000000)#训练函数
except:
flag=0
if(flag==1):
a=popt[0]
b=popt[1]
c=popt[2]
yvals=[func(x[i],a,b,c) for i in range(len(x))]
for ii in range(0,len(x)):
t=yvals[ii]
aad=aad+np.abs((y[ii]-t)/y[ii])
y_flag.insert(ii,t)
try:
res_1=min(res_1,aad/len(x))
if(res_1==aad/len(x)):
y_res=y_flag
else:
y_res=y_res
except:
aaa=1
print("******1111111111***")
print(str(fileName)+" Tc:"+str(tc)+" aad:"+str(res_1)+" 参数a:"+str(a)+" 参数b:"+str(b)+" 参数c:"+str(c))
print("**************22222*******")
else:
a=2.2
b=2.2
res_1=2.2
res_2=2.2
print("******1111111111***")
print("2.2.2.2.2.2.2.2.2.2..2.2.2.2.2.2.2.2.2.2.2.2.2.2.2.2.2")
print("**************22222*******")
else:
a=1.1
b=1.1
res_1=1.1
res_2=1.1
print("******1111111111***")
print("1.1.1.1.1.1.1.1.1..1.1.1.1..1.1.1.1..1.1.1.1.1..11.")
print("**************22222*******")
book = xlrd.open_workbook("e:\\cs\\20191226_论文\\论文\\cs\\paper_help\\paper1_Eq8_logy\\paper_new.xls")
sheet = book.sheet_by_index(0)
rowNum = sheet.nrows
colNum = sheet.ncols
newbook = copy(book)
newsheet = newbook.get_sheet(0) # 在末尾增加新行
global hang
global t_m
global tc
newsheet.write(hang,0,fileName)
newsheet.write(hang,1,a)
newsheet.write(hang,2,b)
newsheet.write(hang,3,c)
newsheet.write(hang,4,res_1)
newsheet.write(hang,5,tc)
hang=hang+1
newbook.save("e:\\cs\\20191226_论文\\论文\\cs\\paper_help\\paper1_Eq8_logy\\paper_new.xls")
file_name("e:\\cs\\20191226_论文\\论文\\cs\\paper_help\\data_liquid\\")
#cal("1,1,1,2,3,3,3-Heptafluoropropane (R227ea)_Liquid.xlsx")
Sec.3
close all
clear all
clc
path = 'D:\20191226_论文\论文\cs\paper_help\data_liquid\';
list = dir([path,'*.xlsx']);
fileNames={list.name};
[temp,dt,raw]=xlsread('D:\20191226_论文\论文\cs\paper_help\20200203_总.xlsx','A1:F60');
[temp_para,dt_para,raw_para]=xlsread('D:\20191226_论文\论文\cs\paper_help\20200203_总.xlsx','A69:J171');
for i = 1:length(fileNames)
data_in = xlsread(fileNames{i});
[num, text, rawN] = xlsread(fileNames{i});
T_all = data_in(:, 1);
eta_all= data_in(:, 12);
Tc=T_all(end);
eta=[];
T=[];
Tf=0;
dt_name=dt(:,1);
for n=1:length(dt(:,1))
if 1==strcmp(dt_name{n},fileNames{i}(1:end-5))
Tf=cell2mat(raw(n,5));
end
end
for n = 2:length(eta_all)
k1 = cell2mat(rawN(n,1));
k2 = cell2mat(rawN(n,12));
if k1>Tf||strcmp(num2str(k2),'undefined')
tt=1;
else
eta=[eta,cell2mat(rawN(n,12))/1000];
T=[T,cell2mat(rawN(n,1))];
end
end
n=0;
if length(T)>=400
n=40;
fileNames{i}
else if length(T)>=350
n=35;
fileNames{i}
else if length(T)>=300
n=30;
fileNames{i}
else if length(T)>=250
n=25;
fileNames{i}
else if length(T)>=200
n=20;
fileNames{i}
else if length(T)>=150
n=15;
fileNames{i}
else if length(T)>=100
n=10;
else if length(T)>=60
n=6;
else if length(T)>=40
n=4;
else if length(T)>=20
n=2;
else
n=1
end
end
end
end
end
end
end
end
end
end
len=fix(length(T)/n);
eta_Selected=[];
T_Selected=[];
temp_k=0;
for k = 1:n:length(eta)-n*4+1
T_Selected=[T_Selected,T(k)];
eta_Selected=[eta_Selected,eta(k)];
temp_k=k;
end
temp_k
n=ceil(n/2);
if length(eta)>4
for k = temp_k+n*2:n:length(eta)
T_Selected=[T_Selected,T(k)];
eta_Selected=[eta_Selected,eta(k)];
end
end
if length(eta)<=4
T_Selected=T;
eta_Selected=eta;
end
eta_Selected;
T_Selected;
a1=0;
b1=0;
a2=0;
b2=0;
c2=0;
a3=0;
b3=0;
c3=0;
dt_name=dt_para(:,1);
for n=1:length(dt_name(:,1))
if 1==strcmp(dt_name{n},fileNames{i}(1:end-5))
a1=cell2mat(raw_para(n,2));
b1=cell2mat(raw_para(n,3));
a2=cell2mat(raw_para(n,4));
b2=cell2mat(raw_para(n,5));
c2=cell2mat(raw_para(n,6));
a3=cell2mat(raw_para(n,7));
b3=cell2mat(raw_para(n,8));
c3=cell2mat(raw_para(n,9));
end
end
Tr_Selected=T_Selected./Tc;
sigma_max=0;
sigma_min=0;
sigma1=a1.*exp(b1./T_Selected);
[sigma_max,index]=max(sigma1);
[sigma_min,index]=min(sigma1);
sigma2=a2.*exp((b2*c2)./(T_Selected-c2));
sigma22=[sigma2,sigma_max];
sigma23=[sigma2,sigma_min];
[sigma_max,index]=max(sigma22);
[sigma_min,index]=min(sigma23);
sigma3=b3./(a3+c3.*log(1-Tr_Selected));
sigma33=[sigma3,sigma_max];
sigma34=[sigma3,sigma_min];
[sigma_max,index]=max(sigma33);
[sigma_min,index]=min(sigma34);
gap=sigma_max-sigma_min;
gap_new=gap/10;
y1=sigma_min-gap_new:0.000001:sigma_max+gap_new;
x1=T(1)+0.*y1;
x2=T(end)+0.*y1;
set(0,'defaultfigurecolor','w') %设置边框背景为白色
%color = {'r', 'g', 'b', 'c', 'm', 'y', 'k', 'w', 'r', 'g', 'b'}; %颜色
color=[ 0 .75 .75 ] ;
linestyle = {'-.', '-', '--', ':', '-.', '-', '--', ':', '-.', '-', '--'}; %线型
markerstyle = {'+', 'o', '*', '.', 'x', 's', 'd', '^', 'v', '<', '>'}; %标记符
T=T_Selected;
sigma1;
plot(T,eta_Selected,'ko')
hold on
plot(T,sigma1,'bs-')
%plot(T,pd2,'color',color(1,:))
plot(T,sigma2,'gd-')
plot(T,sigma3,'r*-')
plot(x1,y1,'Color', [ .8 .8 .8 ],'Linestyle',char(linestyle(3)))
plot(x2,y1,'Color', [ .8 .8 .8 ],'Linestyle',char(linestyle(3)))
xlabel({'$T(K)$'},'interpreter','latex');
ylabel1=ylabel({'$\eta(cp)$'},'interpreter','latex');
set(ylabel1,'Rotation',0);%把Y坐标的标注翻转
legend('REFPROP Data','Eq. (1)','Eq. (2)','Eq. (6)')%添加图例
legend boxoff %去掉图例的边框
set(gca,'Fontname','times new Roman'); %设置当前画图的字体的格式
gap=sigma_max-sigma_min;
gap_new=gap/10;
ylim([sigma_min-gap_new,sigma_max+gap_new]);
gap=T_Selected(end)-T_Selected(1);
gap_new=gap/10;
xlim([T_Selected(1)-gap_new,T_Selected(end)+gap_new]);
saveas(gcf,['D:\20191226_论文\论文\cs\paper_help\picture_eta\',fileNames{i}(1:end-5),'.eps']);
pause(200);
% close;
% pause(1); % 进行延时操作
end
Sec.4
close all
clear all
clc
path = 'D:\20191226_论文\论文\cs\paper_help\data_liquid\';
list = dir([path,'*.xlsx']);
fileNames={list.name};
[temp,dt,raw]=xlsread('D:\20191226_论文\论文\cs\paper_help\20200203_总.xlsx','A1:F60');
[temp_para,dt_para,raw_para]=xlsread('D:\20191226_论文\论文\cs\paper_help\20200203_总.xlsx','A69:J171');
for i = 1:length(fileNames)
data_in = xlsread(fileNames{i});
[num, text, rawN] = xlsread(fileNames{i});
T_all = data_in(:, 1);
eta_all= data_in(:, 12);
Tc=T_all(end);
eta=[];
T=[];
Tf=0;
dt_name=dt(:,1);
for n=1:length(dt(:,1))
if 1==strcmp(dt_name{n},fileNames{i}(1:end-5))
Tf=cell2mat(raw(n,5));
end
end
for n = 2:length(eta_all)
k1 = cell2mat(rawN(n,1));
k2 = cell2mat(rawN(n,12));
if k1>Tf||strcmp(num2str(k2),'undefined')
tt=1;
else
eta=[eta,cell2mat(rawN(n,12))/1000];
T=[T,cell2mat(rawN(n,1))];
end
end
n=0;
if length(T)>=400
n=40;
fileNames{i}
else if length(T)>=350
n=35;
fileNames{i}
else if length(T)>=300
n=30;
fileNames{i}
else if length(T)>=250
n=25;
fileNames{i}
else if length(T)>=200
n=20;
fileNames{i}
else if length(T)>=150
n=15;
fileNames{i}
else if length(T)>=100
n=10;
else if length(T)>=60
n=6;
else if length(T)>=40
n=4;
else if length(T)>=20
n=2;
else
n=1
end
end
end
end
end
end
end
end
end
end
len=fix(length(T)/n);
eta_Selected=[];
T_Selected=[];
temp_k=0;
for k = 1:n:length(eta)-n*4+1
T_Selected=[T_Selected,T(k)];
eta_Selected=[eta_Selected,eta(k)];
temp_k=k;
end
temp_k
n=ceil(n/2);
if length(eta)>4
for k = temp_k+n*2:n:length(eta)
T_Selected=[T_Selected,T(k)];
eta_Selected=[eta_Selected,eta(k)];
end
end
if length(eta)<=4
T_Selected=T;
eta_Selected=eta;
end
eta_Selected;
T_Selected;
a1=0;
b1=0;
a2=0;
b2=0;
c2=0;
a3=0;
b3=0;
c3=0;
dt_name=dt_para(:,1);
for n=1:length(dt_name(:,1))
if 1==strcmp(dt_name{n},fileNames{i}(1:end-5))
a1=cell2mat(raw_para(n,2));
b1=cell2mat(raw_para(n,3));
a2=cell2mat(raw_para(n,4));
b2=cell2mat(raw_para(n,5));
c2=cell2mat(raw_para(n,6));
a3=cell2mat(raw_para(n,7));
b3=cell2mat(raw_para(n,8));
c3=cell2mat(raw_para(n,9));
end
end
Tr_Selected=T_Selected./Tc;
pd_max=0;
pd_min=0;
sigma1=a1.*exp(b1./T_Selected);
pd1=(sigma1-eta_Selected)./eta_Selected;
[pd_max,index]=max(pd1);
[pd_min,index]=min(pd1);
sigma2=a2.*exp((b2*c2)./(T_Selected-c2));
pd2=(sigma2-eta_Selected)./eta_Selected;
pd22=[pd2,pd_max];
pd23=[pd2,pd_min];
[pd_max,index]=max(pd22);
[pd_min,index]=min(pd23);
sigma3=b3./(a3+c3.*log(1-Tr_Selected));
pd3=(sigma3-eta_Selected)./eta_Selected;
pd33=[pd3,pd_max];
pd34=[pd3,pd_min];
[pd_max,index]=max(pd33);
[pd_min,index]=min(pd34);
pd_max;
pd_min;
set(0,'defaultfigurecolor','w') %设置边框背景为白色
T=T_Selected;
sigma1;
color=[ 0 .75 .75 ] ;
plot(T,pd1,'bs-')
hold on
%plot(T,pd2,'color',color(1,:))
plot(T,pd2,'gd-')
plot(T,pd3,'r*-')
xlabel({'$T(K)$'},'interpreter','latex');
ylabel1=ylabel({'$PD$'},'interpreter','latex');
set(ylabel1,'Rotation',0);%把Y坐标的标注翻转
legend('Eq. (1)','Eq. (2)','Eq. (6)')%添加图例
legend boxoff %去掉图例的边框
set(gca,'Fontname','times new Roman'); %设置当前画图的字体的格式
gap=pd_max-pd_min;
gap_new=gap/10;
ylim([pd_min-gap_new,pd_max+gap_new]);
gap=T_Selected(end)-T_Selected(1);
gap_new=gap/10;
xlim([T_Selected(1)-gap_new,T_Selected(end)+gap_new]);
saveas(gcf,['D:\20191226_论文\论文\cs\paper_help\picture\',fileNames{i}(1:end-5),'.eps']);
pause(120);
close;
pause(1); % 进行延时操作
end