Matlab绘图代码共享-硕士阶段

%% 基础信息
% 反复修改的代码简直屎山！
clc;clear
close all


% color1=[77,133,189]/255;
% color2=[247,144,61]/255;
% color3=[89,169,90]/255;

% 马卡龙配色
color1=[147 213 220]/255;
color2=[242 202 201]/255;
color3=[198 223 200]/255;

% color1=[255 230 215]/255;
% color2=[225 240 219]/255;
% color3=[222 234 248]/255;

% color1=[240   181   125]/255;
% color2=[211   225   174]/255;
% color3=[113   171   182]/255;
% 
% color1=[   76   139   192]/255;
% color2=[  232    49    51]/255;
% color3=[   94   183    91]/255;

% 大论文使用的配色
% color1=[  110 158 206]/255;
% color2=[  118 186 128]/255;
% color3=[  230 146 143]/255;

FontSize=13;
color3s={color1;color2;color3};
colors={color1;color1;color1;color2;color2;color2;color3;color3;color3};
colors9=[  
   199    92   100
  240   181   125
  211   225   174
  113   171   182
   75    90   161
   76   139   192
  232    49    51
  163    96   173
   94   183    91]/255;


ki=[ 1	0	0
    0.75	0.25	0
    0.5	0.5	0
    0.25	0.75	0
    0	1	0
    0	0.75	0.25
    0	0.5	0.5
    0	0.25	0.75
    0	0	1];
bar_colors=ki*[color1;color3;color2];

legends=["v=240mm/s  T=210°C  h=0.16mm"
    "v=250mm/s  T=215°C  h=0.20mm"
    "v=260mm/s  T=220°C  h=0.24mm"
    "v=240mm/s  T=215°C  h=0.24mm"
    "v=250mm/s  T=220°C  h=0.16mm"
    "v=260mm/s  T=210°C  h=0.20mm"
    "v=240mm/s  T=220°C  h=0.20mm"
    "v=250mm/s  T=210°C  h=0.24mm"
    "v=260mm/s  T=215°C  h=0.16mm"]';
legends9=[legends;repmat("",5,9)];
%%

hexc=[
'32aeec'
'E3863c'
'E9262e'];

r=hex2dec(hexc(:,[1,2]));
g=hex2dec(hexc(:,[3 4]));
b=hex2dec(hexc(:,[5 6]));
rgb=[r,g,b];
%% 拉伸


Tensile_x=cell(9,6);       % read strain to Tensile_x
Tensile_y=cell(9,6);       % read stress to Tensile_y

% 读取原始数据并储存

% 设置导入选项并导入数据
opts = delimitedTextImportOptions("NumVariables", 8);
% 指定范围和分隔符
opts.DataLines = [8, Inf];
opts.Delimiter = "\t";
% 指定列名称和类型
opts.VariableNames = ["Load", "LoadY", "Position", "Displacement", "Extension", "ExtensionY", "StepIndex", "Time"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double"];
% 指定文件级属性
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
T_L0=50;              % 标距L0
T_S0=4*13;     % 横截面积S0
for Gi=1:9
    for ti=1:6
        filesname=sprintf("/Users/mawenqian/Documents/硕士/VSL/科研/小论文1/基本试件打印参数正交试验/G%d/N%d.txt",Gi,ti);
        readdata = readtable(filesname, opts);
        % 转换为输出类型
        readdata = table2array(readdata);
        Tensile_x{Gi,ti}=readdata(:,5)/T_L0*100;   % 应变 6sh/L^2
        Tensile_y{Gi,ti}=readdata(:,1)/T_S0;        % 应力 3Fl/2bh^2
        clear readdata filesname
    end
end
clear opts
save Tensile_raw_data Tensile_x Tensile_y
%%

% 通过下方原始图像人工判断获得，亦可绘制后展示。
Tensile_start_line=[
    133	137	137	144	145	139
    151	151	147	157	160	155
    176	178	177	179	191	174
168	205	143	194	162	192
157	162	163	154	222	170
196	242	175	206	246	174
    143	141	145	130	138	146
    133	129	135	126	130	127
145	150	163	145	138	142];
Tensile_end_line=[
    2022	1764	3568	2073	2596	3338
    2222	1912	2363	3740	1927	1832
    2613	1821	2390	2353	2041	2040
4000	3000	2449	2667	2557	1937
2000	1943	2000	1776	2000	2000
1719	1719	1719	1719	1719	1719
  2312	1975	2537	2161	1730	1979
1407	1238	1345	1260	1593	1684
1656	2419	2000	2000	2000	2000];
Tensile_maxy_line=[
740	709	744	730	755	716
1066	1103	1064	1167	1058	1157
834	709	735	779	779	723
1377	1441	1622	1512	1698	1399
1043	1074	1069	1140	1072	1036
1214	1142	1164	1214	1142	1164
1178	1099	1213	1252	1191	1182
1097	1091	1081	1035	1100	1063
853	903	920	896	891	922];

% ttt=1;
% for Gi=[ 7 8]
% for ti=1:6
%     subplot(4,6,ttt)
% 
% 
% x=Tensile_x{Gi,ti};
% y=Tensile_y{Gi,ti};
% [~,Tensile_maxy_line(Gi,ti)]=max(y);
% plot(y,MarkerIndices=Tensile_end_line(Gi,ti),MarkerSize=6,Marker="*",MarkerFaceColor='r')
% title(sprintf("%d-%d",Gi,ti))
%     ttt=ttt+1;
% end
% 
% end

%%
% % 绘制原始数据图像总图3x3
% figure()
% set(gcf,'unit','centimeters','Position',[5 5 30 30]);
% for Gi=1:9
%     subplot(3,3,Gi)
%     hold on
%     for ti=1:6
%         x=Tensile_x{Gi,ti}(Tensile_start_line(Gi,ti):Tensile_end_line(Gi,ti));x=x-x(1);
%         y=Tensile_y{Gi,ti}(Tensile_start_line(Gi,ti):Tensile_end_line(Gi,ti));y=y-y(1);
% Tensile_x{Gi,ti}=x;
% Tensile_y{Gi,ti}=y;
% plot(Tensile_y{Gi,ti},LineWidth=1.5)
% %         plot(Tensile_x{Gi,ti},Tensile_y{Gi,ti},LineWidth=1.5,Marker='o',MarkerEdgeColor='red',MarkerIndices=[Tensile_start_line(Gi,ti)])
%         %         plot(Tensile_y{Gi,ti},LineWidth=1.5)
%     end
%     xlabel("Strain (%)")
%     ylabel("Stress (MPa)")
% %     xlim([0 0.1])
%     % ylim([0,70])
%     % legend(Location="northwest",Box="off");
%     box('on');
%     title(['Tensile ',num2str(Gi)])
%     set(gca,'FontName','Arial','FontSize',FontSize)
% end
%%
%
% % 分组绘制原始图像9张单独图（应力）
% for Gi=[1 2 3 7 8]
%     figure()
%     set(gcf,'unit','centimeters','Position',[5 5 40 28]);
%     for ti=1:6
%         subplot(2,3,ti)
% %         [~,maxy]=max(Tensile_y{Gi,ti}(1:Tensile_end_line(Gi,ti)));
% %         Tensil_maxy_line(Gi,ti)=maxy;
%         plot(Tensile_y{Gi,ti},LineWidth=1.5,Marker='o',MarkerEdgeColor='red',MarkerIndices=[Tensile_start_line(Gi,ti),Tensil_maxy_line(Gi,ti),Tensile_end_line(Gi,ti)])
%         xlabel("N")
%         ylabel("Stress (MPa)")
%         % xlim([0 10])
%     % ylim([0,70])
% %     legend(Location="northwest",Box="off");
%     box('on');
%
%     title(['Tensile ',num2str(Gi), '-',num2str(ti)])
%     set(gca,'FontName','Arial','FontSize',FontSize)
%     end
% end

% 绘制受拉总图
%%
Tensile_maxy=zeros(9,6);

clear fh1
fh1=figure();
set(gcf,'unit','centimeters','Position',[15 15 15 15]);
axes('parent',fh1,'unit','centimeters','position',[1.5 1.5 12.5 12.5]);
hold on
LineStyles=repmat(["-",":","--"],1,3);
for Gi=1:9
    for ti=1:6
        x=Tensile_x{Gi,ti}(Tensile_start_line(Gi,ti):Tensile_end_line(Gi,ti));x=x-x(1);
        y=Tensile_y{Gi,ti}(Tensile_start_line(Gi,ti):Tensile_end_line(Gi,ti));y=y-y(1);
        Tensile_maxy(Gi,ti)=max(y);
        plot(x,y,LineWidth=1.5,Color=colors{Gi},LineStyle=LineStyles(Gi))
        clear x y
    end
    
    xlabel("Strain (%)")
    ylabel("Stress (MPa)")
    xlim([0 8])
    ylim([0,60])
 legend(legends9,Location="northeast",Box="off");
    box('on');
 
    set(gca,'FontName','Arial','FontSize',FontSize)  
end



%% 压缩1屈服强度
compressive1_x=cell(9,6);       % read strain to compressive1_x
compressive1_y=cell(9,6);       % read stress to compressive1_y
Com1_L0=25.4;                   % 标距L0
Com1_S0=12.7*12.7*pi*0.25;      % 横截面积S0


% 读取原始数据并储存
% 设置导入选项并导入数据
opts = spreadsheetImportOptions("NumVariables", 2);
% 指定工作表和范围
opts.Sheet = "Sheet1";
opts.DataRange = "A2:B2000";
% 指定列名称和类型
opts.VariableNames = ["LoadValue", "PositionValue", "PlayTime", "ExtendValue"];
opts.VariableTypes = ["double", "double", "double", "double"];
for Gi=19:27
    for ti=1:6
        filesname=sprintf("/Users/mawenqian/Documents/硕士/VSL/科研/小论文1/基本试件打印参数正交试验/G%d/N%d.xls",Gi,ti);
        readdata = readtable(filesname, opts);
        % 转换为输出类型
        readdata = table2array(readdata);
        compressive1_x{Gi-18,ti}=readdata(:,2)/Com1_L0*100;  % 应变
        compressive1_y{Gi-18,ti}=readdata(:,1)/Com1_S0;      % 应力
        clear readdata filesname
    end
end
clear opts
% 通过下方原始图像人工判断获得，亦可绘制后展示。
compressive1_start_line=[
    73	78	74	77	12	1
    77	87	79	78	75	77
    82	82	80	82	86	83
    369	70	64	72	72	74
    71	67	66	63	64	66
    70	70	68	67	69	69
    62	6	6	5	4	4
    3	3	3	181	1	7
    65	12	5	4	5	64];
compressive1_end_line=[
    471	413	432	390	303	322
    422	432	413	380	380	378
    471	442	480	390	303	322
    827	529	493	497	503	482
    511	526	453	475	496	496
    513	509	528	486	473	494
    528	516	475	414	463	470
    463	454	519	605	463	457
    619	497	496	365	420	478];
compressive1_maxy_line=[310	307	318	314	234	219
    354	377	344	374	352	345
    321	325	333	324	303	322
    728	412	381	403	404	404
    357	350	340	344	345	347
    349	330	332	341	329	339
    331	295	293	279	286	314
    306	285	300	499	313	293
    436	274	300	289	298	323];
% % 绘制原始数据图像总图3x3
% figure()
% set(gcf,'unit','centimeters','Position',[5 5 30 30]);
% for Gi=19:27
%     subplot(3,3,Gi-18)
%     hold on
%     for ti=1:6
%         [~,maxy]=max(compressive1_y{Gi-18,ti}(1:compressive1_end_line(Gi-18,ti)));
%         plot(compressive1_x{Gi-18,ti},compressive1_y{Gi-18,ti},LineWidth=1.5,Marker='o',MarkerEdgeColor='red',MarkerIndices=[maxy])
%         %         plot(compressive1_y{Gi-18,ti},LineWidth=1.5)
%     end
%     xlabel("Strain (%)")
%     ylabel("Stress (MPa)")
%     % xlim([0 10])
%     % ylim([0,70])
%     %     legend(Location="northwest",Box="off");
%     box('on');
%     title(['compressive1 ',num2str(Gi-18)])
%     set(gca,'FontName','Arial','FontSize',FontSize)
% end
% % 分组绘制原始图像9张单独图（应力）
% for Gi=27:-1:19
%     figure()
%     set(gcf,'unit','centimeters','Position',[5 5 40 28]);
%     for ti=1:6
%         subplot(2,3,ti)
%         [~,maxy]=max(compressive1_y{Gi-18,ti}(1:compressive1_end_line(Gi-18,ti)));
% compressive1_maxy_line(Gi,ti)=maxy;
%         plot(compressive1_y{Gi-18,ti},LineWidth=1.5,Marker='o',MarkerEdgeColor='red',MarkerIndices=[compressive1_start_line(Gi-18,ti),compressive1_maxy_line(Gi,ti),compressive1_end_line(Gi-18,ti)])
%         xlabel("N")
%         ylabel("Stress (MPa)")
%         % xlim([0 10])
%         % ylim([0,70])
%         %     legend(Location="northwest",Box="off");
%         box('on');
%
%         title(['compressive1 ',num2str(Gi-18), '-',num2str(ti)])
%         set(gca,'FontName','Arial','FontSize',FontSize)
%     end
% end



% 绘制受压总图
compressive1_maxy=zeros(9,6);

% figure()
% set(gcf,'unit','centimeters','Position',[5 5 30 30]);
% for Gi=19:27
%     Gi=Gi-18;
%     subplot(3,3,Gi)
% 
%     hold on
%     for ti=1:6
%         compressive1_maxy(Gi,ti)=compressive1_y{Gi,ti}(compressive1_maxy_line(Gi,ti))-compressive1_x{Gi,ti}(compressive1_start_line(Gi,ti));
%         x=compressive1_x{Gi,ti}(compressive1_start_line(Gi,ti):compressive1_end_line(Gi,ti));x=x-x(1);
%         y=compressive1_y{Gi,ti}(compressive1_start_line(Gi,ti):compressive1_end_line(Gi,ti));y=y-y(1);
%         compressive1_maxy(Gi,ti)=max(y);
%         plot(x,y,LineWidth=1.5,Color=colors{Gi})
%     end
%     xlabel("Strain (%)")
%     ylabel("Stress (MPa)")
%     xlim([0 8])
%     ylim([0,60])
% 
%     legend(legends(Gi),Location="north",Box="off");
%     box('on');
%     title(['Compressive ',num2str(Gi)])
%     set(gca,'FontName','Arial','FontSize',FontSize)
% end

clear fh1
fh1=figure();
set(gcf,'unit','centimeters','Position',[15 15 15 15]);
axes('parent',fh1,'unit','centimeters','position',[1.5 1.5 12.5 12.5]);
hold on
LineStyles=repmat(["-",":","--"],1,3);
for Gi=19:27
    Gi=Gi-18;
    hold on
    for ti=1:6
        x=compressive1_x{Gi,ti}(compressive1_start_line(Gi,ti):compressive1_end_line(Gi,ti));x=x-x(1);
        y=compressive1_y{Gi,ti}(compressive1_start_line(Gi,ti):compressive1_end_line(Gi,ti));y=y-y(1);
        compressive1_maxy(Gi,ti)=max(y);
        plot(x,y,LineWidth=1.5,Color=colors{Gi},LineStyle=LineStyles(Gi))
    end
    xlabel("Strain (%)")
    ylabel("Stress (MPa)")
    xlim([0 8])
    ylim([0,60])
 legend(legends9,Location="southeast",Box="off");
    box('on');
%     title('Compressive')
    set(gca,'FontName','Arial','FontSize',FontSize)
end
%% 弯曲


Flexural_x=cell(9,6);       % read strain to Flexural_x
Flexural_y=cell(9,6);       % read stress to Flexural_y

% 读取原始数据并储存

% 设置导入选项并导入数据
opts = delimitedTextImportOptions("NumVariables", 8);
% 指定范围和分隔符
opts.DataLines = [8, Inf];
opts.Delimiter = "\t";
% 指定列名称和类型
opts.VariableNames = ["Load", "LoadY", "Position", "Displacement", "Extension", "ExtensionY", "StepIndex", "Time"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double"];
% 指定文件级属性
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
L0=6*3.2/52^2;              % 标距L0
S0=3*52/(2*12.7*3.2^2);     % 横截面积S0
for Gi=10:18
    for ti=1:6
        filesname=sprintf("/Users/mawenqian/Documents/硕士/VSL/科研/小论文1/基本试件打印参数正交试验/G%d/N%d.txt",Gi,ti);
        readdata = readtable(filesname, opts);
        % 转换为输出类型
        readdata = table2array(readdata);
        Flexural_x{Gi-9,ti}=readdata(:,4).*L0*100;   % 应变 6sh/L^2
        Flexural_y{Gi-9,ti}=readdata(:,1).*S0;        % 应力 3Fl/2bh^2
        clear readdata filesname
    end
end
clear opts
% 通过下方原始图像人工判断获得，亦可绘制后展示。
Flexural_start_line=[
    467	518	522	779	368	484
    519	337	270	390	415	252
    517	357	717	515	553	654
    770	553	453	344	498	318
    498	521	569	667	510	417
    360	431	407	364	472	358
    308	493	677	617	648	471
    379	529	641	781	915	439
    693	701	745	696	734	744];
Flexural_maxy_line=[
    4554	4543	4534	4729	4459	4457
    4391	4252	4080	4298	4264	4114
    4496	4306	4691	4533	4453	4458
    4431	4205	4134	3991	4242	4116
    4529	4633	4710	4794	4624	4509
    4078	4153	4095	4148	4172	4140
    4059	4267	4356	4334	4289	4146
    4165	4235	4366	4562	4611	4207
    4725	4725	4725	4725	4725	4725];

% % 绘制原始数据图像总图3x3
% figure()
% set(gcf,'unit','centimeters','Position',[5 5 30 30]);
% for Gi=10:18
%     subplot(3,3,Gi-9)
%     hold on
%     for ti=1:6
%         plot(Flexural_x{Gi-9,ti},Flexural_y{Gi-9,ti},LineWidth=1.5,Marker='o',MarkerEdgeColor='red',MarkerIndices=[Flexural_start_line(Gi-9,ti)])
%         %         plot(Flexural_y{Gi-9,ti},LineWidth=1.5)
%     end
%     xlabel("Strain (%)")
%     ylabel("Stress (MPa)")
%     % xlim([0 10])
%     % ylim([0,70])
%     % legend(Location="northwest",Box="off");
%     box('on');
%     title(['Flexural ',num2str(Gi)])
%     set(gca,'FontName','Arial','FontSize',FontSize)
% end
%
% % 分组绘制原始图像9张单独图（应力）
% for Gi=18:-1:10
%     figure()
%     set(gcf,'unit','centimeters','Position',[5 5 40 28]);
%     for ti=1:6
%         subplot(2,3,ti)
%         [~,maxy]=max(Flexural_y{Gi-9,ti});
%         Flexural_maxy_line(Gi-9,ti)=maxy;
%         plot(Flexural_y{Gi-9,ti},LineWidth=1.5,Marker='o',MarkerEdgeColor='red',MarkerIndices=[Flexural_start_line(Gi-9,ti),Flexural_maxy_line(Gi-9,ti)])
%         xlabel("N")
%         ylabel("Stress (MPa)")
%         % xlim([0 10])
%     % ylim([0,70])
% %     legend(Location="northwest",Box="off");
%     box('on');
%
%     title(['Flexural ',num2str(Gi-9), '-',num2str(ti)])
%     set(gca,'FontName','Arial','FontSize',FontSize)
%     end
% end

% 绘制受弯总图
Flexural_maxy=zeros(9,6);

% figure()
% set(gcf,'unit','centimeters','Position',[5 5 30 30]);
% for Gi=10:18
%     Gi=Gi-9;
%  subplot(3,3,Gi)
% 
%     hold on
%     for ti=1:6
%         Flexural_maxy(Gi,ti)=Flexural_y{Gi,ti}(Flexural_maxy_line(Gi,ti))-Flexural_x{Gi,ti}(Flexural_start_line(Gi,ti));
%         x=Flexural_x{Gi,ti}(Flexural_start_line(Gi,ti):end);x=x-x(1);
%         y=Flexural_y{Gi,ti}(Flexural_start_line(Gi,ti):end);y=y-y(1);
%         Flexural_maxy(Gi,ti)=max(y);
%         plot(x,y,LineWidth=1.5,Color=colors{Gi})
%     end
%     xlabel("Strain (%)")
%     ylabel("Stress (MPa)")
%     xlim([0 8])
%     ylim([0,70])
% 
%     legend(legends(Gi),Location="north",Box="off");
%     box('on');
%     title(['Flexural ',num2str(Gi)])
%     set(gca,'FontName','Arial','FontSize',FontSize)
% end
% 
% 
clear fh1
fh1=figure();
set(gcf,'unit','centimeters','Position',[15 15 15 15]);
axes('parent',fh1,'unit','centimeters','position',[1.5 1.5 12.5 12.5]);
hold on
LineStyles=repmat(["-",":","--"],1,3);
for Gi=10:18
    Gi=Gi-9;

    hold on
    for ti=1:6
        x=Flexural_x{Gi,ti}(Flexural_start_line(Gi,ti):end);x=x-x(1);
        y=Flexural_y{Gi,ti}(Flexural_start_line(Gi,ti):end);y=y-y(1);
        Flexural_maxy(Gi,ti)=max(y);
        plot(x,y,LineWidth=1.5,Color=colors{Gi},LineStyle=LineStyles(Gi))
    end

end

xlabel("Strain (%)")
    ylabel("Stress (MPa)")
    xlim([0 8])
    ylim([0,70])

    legend(legends9,Location="southeast",Box="off");
    box('on');
%     title(['Flexural'])
    set(gca,'FontName','Arial','FontSize',FontSize)

%% 汇总柱状图
% close all
com_bar=figure();
set(gcf,'unit','centimeters','Position',[5 5 28 28]);
axes('parent',com_bar,'unit','centimeters','position',[1.5 1.5 25 25],'fontsize',20);

x=1:9;
all_mean=reshape(mean([Tensile_maxy;compressive1_maxy;Flexural_maxy],2),9,3)';
all_mean_std=reshape(std([Tensile_maxy;compressive1_maxy;Flexural_maxy],0,2),9,3)';

bar_fig= bar(x,all_mean,FaceColor="flat",BarWidth=0.8);
hold on


for typei=1:3

    set(bar_fig(typei),'FaceColor',color3s{typei});

    errorbar(x+0.225*(typei-2),all_mean(typei,:),all_mean_std(typei,:),'LineStyle','none',LineWidth=1.5,Color=[0.5 0.5 0.5]);

    xtips = bar_fig(typei).XEndPoints;
    ytips = bar_fig(typei).YEndPoints+all_mean_std(typei,:);
    label = string(vpa(bar_fig(typei).YData,3));
    text(xtips,ytips,label,'HorizontalAlignment','center','VerticalAlignment','bottom','FontSize',12)
end
legend(["Tensile Yield Strength" "Compressive Yield Strength" "Flexural Yield Strength" ],Location="northwest",box="on")

%     xlabel("Type")
ylabel("Stress (MPa)")
set(gca,'FontName','Arial','FontSize',14)
%% 单独柱状图
figure();
ti=1;
set(gcf,'unit','centimeters','Position',[3 3 15 15]);
axes('unit','centimeters','position',[1.5 1.5 12.5 12.5]);
bar_fig_i= bar(x,all_mean(ti,:),BarWidth=0.8,FaceColor=color2);
hold on
errorbar(x,all_mean(ti,:),all_mean_std(ti,:),'LineStyle','none',LineWidth=1.5,Color=[0.5 0.5 0.5]);

xtips = bar_fig_i.XEndPoints;
ytips = bar_fig_i.YEndPoints+all_mean_std(ti,:);
label = string(vpa(bar_fig_i.YData,3));
text(xtips,ytips,label,'HorizontalAlignment','center','VerticalAlignment','bottom')


ylabel("Stress (MPa)")
set(gca,'FontName','Arial','FontSize',FontSize)
%% 保存计算获得的应力矩阵
all_sigma=[Tensile_maxy(:),compressive1_maxy(:),Flexural_maxy(:)];

save sigma_mat all_sigma all_mean all_mean_std  Tensile_maxy compressive1_maxy Flexural_maxy


%% 加载应力矩阵
load('sigma_mat.mat')
A=repmat([1	1	1	2	2	2	3	3	3]',6,1);
B=repmat([1	2	3	1	2	3	1	2	3]',6,1);
C=repmat([1	2	3	2	3	1	3	1	2]',6,1);
D=repmat([1	2	3	3	1	2	2	3	1]',6,1);
save taguchi9_mat A B C D % 保存相关变量
%% SN方差分析
sn=-10*log10(1./all_sigma.^2);
anova_soure_tbl= [{'A','B','C','D','Tensile','Compressive','Flexural'};[num2cell([A B C D]), num2cell(sn)]];
all_anova=cell(3,1);
for ni=1:3
    [~,tbli]=anovan(sn(:,ni),{A,B,C,D},"varnames",{'A','B','C','D'});
    tbli{1,1}='Factor';
    percent_source=cell2mat(tbli(2:6,5));
    percent=[{'Contribution'};num2cell(round(100*[percent_source;sum(percent_source)]./sum(percent_source),2))];


    all_anova{ni}=[tbli(:,[1 2 3 5 6 7]),percent];

end
save anovan_mat  all_anova anova_soure_tbl

%% 信噪比分析（采用越大越好计算公式）
all_sigma2square=1./all_sigma.^2;
all_sn=cell(3,1);
% all_ki=cell(3,1);
for ni=1:3
    [tbl_A,grp_A]=grpstats(all_sigma2square(:,ni),A,["mean","gname"]);
    [tbl_B,grp_B]=grpstats(all_sigma2square(:,ni),B,["mean","gname"]);
    [tbl_C,grp_C]=grpstats(all_sigma2square(:,ni),C,["mean","gname"]);
    [tbl_D,grp_D]=grpstats(all_sigma2square(:,ni),D,["mean","gname"]);

    sni=-10*log10([tbl_A,tbl_B,tbl_C,tbl_D]);
    sni_delta=max(sni)-min(sni);
    [~,sni_sort]=sort(sni_delta,'descend'); % 降序排列
    [~,sni_sort]=sort(sni_sort);           % 获得排位
    all_sn{ni}=cell(size(sni)+[3 1]);
    all_sn{ni}(1,:)={'Level' 'A','B','C','D'};
    all_sn{ni}(2:end,1)={'1' '2' '3' 'Delta' 'Rank'};
    all_sn{ni}(2:end,2:end)=[num2cell([sni;sni_delta]);num2cell(sni_sort)];

    % 直观分析计算过程（不采用）
    % tbl_A=grpstats(all_sigma(:,ni),A,@(x)sum(x,1));
    % tbl_B=grpstats(all_sigma(:,ni),B,@(x)sum(x,1));
    % tbl_C=grpstats(all_sigma(:,ni),C,@(x)sum(x,1));
    % tbl_D=grpstats(all_sigma(:,ni),D,@(x)sum(x,1));
    % all_ki{ni}=[tbl_A,tbl_B,tbl_C,tbl_D];

end
save all_sn_mat all_sn
%% 绘制信噪比SN曲线
close all
which_sigma=1; % 1:Tensile 2:Compressive1 3:Flexural
y1=cell2mat(all_sn{which_sigma}(2:4,2));
y2=cell2mat(all_sn{which_sigma}(2:4,3));
y3=cell2mat(all_sn{which_sigma}(2:4,4));
y4=cell2mat(all_sn{which_sigma}(2:4,5));
% y1=[0 0 0 ];
% y2=[0 0 0 ];
% y3=[0 0 0 ];
% y4=[0 0 0 ];
plot([y1;y2;y3;y4])
ylim_range=[18 32]; % set ylim range from the upper figure


FontSize=13;
fh1=figure();
set(gcf,'unit','centimeters','Position',[2 2 23 11]);

axes('parent',fh1,'unit','centimeters','position',[2 2 5 8],'fontsize',10);
plot([1 2 3],y1,'linewidth',3,'Marker','.','MarkerSize',30);
xlim([0.5 3.5])
ylim(ylim_range)
set(gca,'XTick',[1 2 3],'yTickLabel',num2str(get(gca,'yTick')','%.1f'));
xlabel("Filament","FontSize",FontSize,'FontName','Arial')
ylabel("Mean of SN Radio","FontSize",FontSize,'FontName','Arial')
box off

axes('parent',fh1,'unit','centimeters','position',[7 2 5 8],'fontsize',10);
plot([240 250 260],y2,'linewidth',3,'Marker','.','MarkerSize',30);
xlim([235 265])
ylim(ylim_range)
set(gca,'XTick',[240 250 260],'YTickLabel','')
xlabel("Speed(mm/s)","FontSize",FontSize,'FontName','Arial')
box off

axes('parent',fh1,'unit','centimeters','position',[12 2 5 8],'fontsize',10);
plot([210 215 220],y3,'linewidth',3,'Marker','.','MarkerSize',30)
xlim([205 225])
ylim(ylim_range)
set(gca,'XTick',[210 215 220],'YTickLabel','')
xlabel("Temperature(\circC)","FontSize",FontSize,'FontName','Arial')
box off

axes('parent',fh1,'unit','centimeters','position',[17 2 5 8],'fontsize',10);
plot([0.16 0.20 0.24],y4,'linewidth',3,'Marker','.','MarkerSize',30)
xlim([0.14 0.26])
ylim(ylim_range)
set(gca,'XTick',[0.16 0.20 0.24],'YTickLabel','')
xlabel("Layer thickness(mm)","FontSize",FontSize,'FontName','Arial')
box off

ax2 = axes('parent',fh1,'unit','centimeters','position',[2 2 20 8],...  
           'XAxisLocation','top',...  
           'YAxisLocation','right',...  
           'Color','none',...  
           'XColor','k','YColor','k');  
set(ax2,'YTick', []);  
set(ax2,'XTick', []); 
%% 线性回归预测
X=[ones(size(A)) A B C D];
all_regress=cell(3,1);
for ni=1:3
    [b,~,~,~,stats] = regress(all_sigma(:,ni),X);
    all_regress{ni}=[{'' 'A' 'B' 'C' 'D' 'R^2' 'p-value'}',num2cell([b;stats([1 3])'])];
end

save all_regress_mat all_regress

clc;clear
close all


% color1=[77,133,189]/255;
% color2=[247,144,61]/255;
% color3=[89,169,90]/255;
color1=[147 213 220]/255;
color2=[242 202 201]/255;
color3=[198 223 200]/255;
FontSize=12;
%% 拉压弯标准对比图
% 试验时间：2022年10月19-24日
% 3个子图

fh1=figure(1);

% Tensile
set(gcf,'unit','centimeters','Position',[5 5 33 13]);
axes('parent',fh1,'unit','centimeters','position',[3,3,8,8],'fontsize',10);
hold on

% 设置导入选项并导入数据
opts = delimitedTextImportOptions("NumVariables", 8);

% 指定范围和分隔符
opts.DataLines = [8, Inf];
opts.Delimiter = "\t";

% 指定列名称和类型
opts.VariableNames = ["Load", "LoadY", "Position", "Displacement", "Extension", "ExtensionY", "StepIndex", "Time"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double"];

% 指定文件级属性
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";

% 标距L0
L0=[ones(1,10).*50,ones(1,5)*7.62];

% 横截面积S
S=[ones(1,5).*40,ones(1,5).*4*13,ones(1,5).*4*3.18];


% 曲线颜色
% colorrgb=[ones(5,1)*color1;ones(5,1)*color2;ones(5,1)*color3];
colorrgb=[repmat(color1,5,1);repmat(color2,5,1);repmat(color3,5,1)];

maxf_T=zeros(5,3);
for filesnum=1:15
    filesname=sprintf("/Users/wenqian/Documents/硕士/VSL/科研/小论文1/各种标准对比试验/实验数据/整理/Tensile/N%d.txt",filesnum);
    readdata = readtable(filesname, opts);
    % 转换为输出类型
    readdata = table2array(readdata);
    x=readdata(:,4)./L0(filesnum)*100;   % 应变
    y=readdata(:,1)./S(filesnum);        % 应力
    plot(x,y,LineWidth=1.5,Color=colorrgb(filesnum,:))
    maxf_T(filesnum)=max(y);
end
ave_stress_T=sum(maxf_T,1)/5;

legendtext=["ISO 527 1B",repmat([""],1,4) ,"ASTM D638 Type I",repmat([""],1,4),"ASTM D638 Type V"];

xlabel("\epsilon (%)")
ylabel("\sigma (MPa)")
xlim([0 10])
ylim([0,50])
legend(legendtext,Location="northwest",Box="off");
box('on');
  title('Tensile')
set(gca,'FontName','Arial','FontSize',FontSize)

clear opts
clear colorrgb

% Compressive
axes('parent',fh1,'unit','centimeters','position',[13,3,8,8],'fontsize',10);
hold on
% 设置导入选项并导入数据
opts = spreadsheetImportOptions("NumVariables", 4);

% 指定工作表和范围
opts.Sheet = "Sheet1";
opts.DataRange = "A2:D5860";

% 指定列名称和类型
opts.VariableNames = ["LoadValue", "PositionValue", "PlayTime", "ExtendValue"];
opts.VariableTypes = ["double", "double", "double", "double"];

% 标距L0
L0=ones(1,10).*50.8;

% 横截面积S
S=[ones(1,5).*12.7*12.7,ones(1,5).*12.7*12.7*pi*0.25];

colorrgb=[repmat(color2,5,1);repmat(color3,5,1)];
maxf_C=zeros(5,2);
for filesnum=1:10
    filesname=sprintf("/Users/wenqian/Documents/硕士/VSL/科研/小论文1/各种标准对比试验/实验数据/整理/Compressive/N%d.xls",filesnum);
    readdata = readtable(filesname, opts, "UseExcel", false);
    % 转换为输出类型
    readdata = table2array(readdata);
    x=readdata(:,2)./L0(filesnum)*100;   % 应变
    y=readdata(:,1)./S(filesnum);        % 应力
    plot(x,y,LineWidth=1.5,Color=colorrgb(filesnum,:))
    maxf_C(filesnum)=max(y);
end
ave_stress_C=sum(maxf_C,1)/5;

legendtext=["ASTM D695 Prise",repmat([""],1,4) ,"ASTM D695 Cylinder"];

xlabel("\epsilon (%)")
ylabel("\sigma (MPa)")
xlim([0 10])
ylim([0,50])
legend(legendtext,Location="northwest",Box="off");
box('on');
  title('Compressive')
set(gca,'FontName','Arial','FontSize',FontSize)

clear opts
clear colorrgb

% Flexural

axes('parent',fh1,'unit','centimeters','position',[23,3,8,8],'fontsize',10);
hold on

% 设置导入选项并导入数据
opts = delimitedTextImportOptions("NumVariables", 8);

% 指定范围和分隔符
opts.DataLines = [8, Inf];
opts.Delimiter = "\t";

% 指定列名称和类型
opts.VariableNames = ["Load", "LoadY", "Position", "Displacement", "Extension", "ExtensionY", "StepIndex", "Time"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double"];

% 指定文件级属性
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";

% 标距L0
L0=[ones(1,5).*6*4/56^2,ones(1,5)*6*3.2/52^2];

% 横截面积S
S=[ones(1,5).*3*56/(2*10*4^2),ones(1,5).*3*52/(2*12.7*3.2^2)];


% 曲线颜色
% colorrgb=[ones(5,1)*color1;ones(5,1)*color2;ones(5,1)*color3];
colorrgb=[repmat(color1,5,1);repmat(color2,5,1)];

maxf_F=zeros(5,2);
for filesnum=1:10
    filesname=sprintf("/Users/wenqian/Documents/硕士/VSL/科研/小论文1/各种标准对比试验/实验数据/整理/Flexural/N-%d.txt",filesnum);
    readdata = readtable(filesname, opts);
    % 转换为输出类型
    readdata = table2array(readdata);
    x=readdata(:,4).*L0(filesnum)*100;   % 应变 6sh/L^2
    y=readdata(:,1).*S(filesnum);        % 应力 3Fl/2bh^2
    plot(x,y,LineWidth=1.5,Color=colorrgb(filesnum,:))
    maxf_F(filesnum)=max(y);
end

ave_stress_F=sum(maxf_F,1)/5;

legendtext=["ISO 178",repmat([""],1,4) ,"ASTM D790"];

xlabel("\epsilon (%)")
ylabel("\sigma (MPa)")
xlim([0 10])
ylim([0,50])
legend(legendtext,Location="northwest",Box="off");
box('on');
  title('Flexural')
set(gca,'FontName','Arial','FontSize',FontSize)

clear opts
clear colorrgb

%% ASTM标准下压、弯不同打印方向对比图
% 试验时间：2022年10月19-24日
% 3个子图

fh2=figure(2);

% Compressive-1
set(gcf,'unit','centimeters','Position',[5 5 33 13]);
axes('parent',fh2,'unit','centimeters','position',[3,3,8,8],'fontsize',10);
hold on

% 设置导入选项并导入数据
opts = spreadsheetImportOptions("NumVariables", 4);

% 指定工作表和范围
opts.Sheet = "Sheet1";
opts.DataRange = "A2:D5860";

% 指定列名称和类型
opts.VariableNames = ["LoadValue", "PositionValue", "PlayTime", "ExtendValue"];
opts.VariableTypes = ["double", "double", "double", "double"];

% 标距L0
L0=[ones(1,10).*50.8];

% 横截面积S
S=[ones(1,10).*12.7*12.7];

colorrgb=[repmat(color2,5,1);repmat(color3,5,1)];
maxf_CN5=zeros(5,2);
for filesnum=1:10
    filesname=sprintf("/Users/wenqian/Documents/硕士/VSL/科研/小论文1/各种标准对比试验/实验数据/整理/compressive-diff-directions/N5-%d.xls",filesnum);
    readdata = readtable(filesname, opts, "UseExcel", false);
    % 转换为输出类型
    readdata = table2array(readdata);
    x=readdata(:,2)./L0(filesnum)*100;   % 应变
    y=readdata(:,1)./S(filesnum);        % 应力
    plot(x,y,LineWidth=1.5,Color=colorrgb(filesnum,:))
    maxf_CN5(filesnum)=max(y);
end
ave_stress_CN5=sum(maxf_CN5,1)/5;

legendtext=["水平打印",repmat([""],1,4) ,"竖直打印"];

xlabel("\epsilon (%)")
ylabel("\sigma (MPa)")
xlim([0 10])
ylim([0,50])
legend(legendtext,Location="northwest",Box="off");
box('on');
  title('ASTM D695 Prise')
set(gca,'FontName','Arial','FontSize',FontSize)

clear opts
clear colorrgb






% Compressive-2
axes('parent',fh2,'unit','centimeters','position',[13,3,8,8],'fontsize',10);
hold on

% 设置导入选项并导入数据
opts = spreadsheetImportOptions("NumVariables", 4);

% 指定工作表和范围
opts.Sheet = "Sheet1";
opts.DataRange = "A2:D5860";

% 指定列名称和类型
opts.VariableNames = ["LoadValue", "PositionValue", "PlayTime", "ExtendValue"];
opts.VariableTypes = ["double", "double", "double", "double"];

% 标距L0
L0=[ones(1,10).*50.8];

% 横截面积S
S=[ones(1,10).*12.7*12.7*pi*0.25];

colorrgb=[repmat(color2,5,1);repmat(color3,5,1)];
maxf_CN6=zeros(5,2);
for filesnum=1:10
    filesname=sprintf("/Users/wenqian/Documents/硕士/VSL/科研/小论文1/各种标准对比试验/实验数据/整理/compressive-diff-directions/N6-%d.xls",filesnum);
    readdata = readtable(filesname, opts, "UseExcel", false);
    % 转换为输出类型
    readdata = table2array(readdata);
    x=readdata(:,2)./L0(filesnum)*100;   % 应变
    y=readdata(:,1)./S(filesnum);        % 应力
    plot(x,y,LineWidth=1.5,Color=colorrgb(filesnum,:))
    maxf_CN6(filesnum)=max(y);
end
ave_stress_CN6=sum(maxf_CN6,1)/5;

legendtext=["水平打印",repmat([""],1,4) ,"竖直打印"];

xlabel("\epsilon (%)")
ylabel("\sigma (MPa)")
xlim([0 10])
ylim([0,50])
legend(legendtext,Location="northwest",Box="off");
box('on');
  title('ASTM D695 Cylinder')
set(gca,'FontName','Arial','FontSize',FontSize)

clear opts
clear colorrgb



% Flexural
axes('parent',fh2,'unit','centimeters','position',[23,3,8,8],'fontsize',10);
hold on

% 设置导入选项并导入数据
opts = delimitedTextImportOptions("NumVariables", 8);

% 指定范围和分隔符
opts.DataLines = [8, Inf];
opts.Delimiter = "\t";

% 指定列名称和类型
opts.VariableNames = ["Load", "LoadY", "Position", "Displacement", "Extension", "ExtensionY", "StepIndex", "Time"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double"];

% 指定文件级属性
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";

% 标距L0
L0=[ones(1,5).*6*4/56^2,ones(1,5)*6*3.2/52^2];

% 横截面积S
S=[ones(1,5).*3*56/(2*10*4^2),ones(1,5).*3*52/(2*12.7*3.2^2)];


% 曲线颜色
% colorrgb=[ones(5,1)*color1;ones(5,1)*color2;ones(5,1)*color3];
colorrgb=[repmat(color2,5,1);repmat(color3,5,1)];

maxf_FN8=zeros(5,2);
for filesnum=1:10
    filesname=sprintf("/Users/wenqian/Documents/硕士/VSL/科研/小论文1/各种标准对比试验/实验数据/整理/Flexural-diff-direction/N-%d.txt",filesnum);
    readdata = readtable(filesname, opts);
    % 转换为输出类型
    readdata = table2array(readdata);
    x=readdata(:,4).*L0(filesnum)*100;   % 应变 6sh/L^2
    y=readdata(:,1).*S(filesnum);        % 应力 3Fl/2bh^2
    plot(x,y,LineWidth=1.5,Color=colorrgb(filesnum,:))
    maxf_FN8(filesnum)=max(y);
end

ave_stress_FN8=sum(maxf_FN8,1)/5;

legendtext=["水平打印",repmat([""],1,4) ,"竖直打印"];

xlabel("\epsilon (%)")
ylabel("\sigma (MPa)")
xlim([0 10])
ylim([0,50])
legend(legendtext,Location="northwest",Box="off");
box('on');
  title('ASTM D790')
set(gca,'FontName','Arial','FontSize',FontSize)

clear opts
clear colorrgb

clc;clear
close all

out_data_file=[...
    4	5	6
    7	8	9
    10	11	12
    13	14	15
    16	19	20
    21	22	23
    27	28	29
    41	30	31
    32	33	34
    24	25	26
    35	36	37
    38	39	40];
group_name={'ITIL-折3';'ITIL-折6';'ITIL-折9';'ITIL-直-3';'ITIL-直-6';'ITIL-直-9';'Zebra-3';'Zebra-6';'Zebra-9';'Theta-15';'复合ITIL-直-9';'复合ITIL-折-9'};
data_size=size(out_data_file);
stress=zeros(data_size);
maxf=[];
ave_xys={};



h1=figure(1);
colorrgb=[77,133,189;247,144,61;89,169,90]/255;
for rows=1:data_size(1)
    subplot(4,3,rows)
    hold on
    legendtext=[];
    xys={};
    for cols=1:data_size(2)
        
        readdata = readmatrix(sprintf("/Users/wenqian/Documents/硕士/VSL/2022.5/PLA- PLA拉伸试验/20220510.is_tens_Exports/20220510_%d.csv",out_data_file(rows,cols)));
        [x,y,stress(rows,cols)]=filter_date(readdata(:,2),readdata(:,3));
%         N->kN
        y=1000*y;
        stress(rows,cols)=stress(rows,cols)*1000;
        plot(x,y,LineWidth=1.5,Color=colorrgb(cols,:))
        legendtext=[legendtext,sprintf("%s-%d",cell2mat(group_name(rows)),cols)];
        maxf(rows,cols)=max(y);
        xys=[xys;{x y}];
    end
    [ax,ay]=average_date(xys);


    ave_xys=[ave_xys;{ax ay}];
    plot(ax,smooth(ay),LineWidth=2.5);
    legendtext=[legendtext,"平均曲线"];
    xlabel("Displacement (mm)")
    ylabel("Force (N)")
    xlim([0 15])
    ylim([0,1250])
    legend(legendtext,Location="southeast",Box="off");
    box('on');
    title(sprintf('(%s)',char(96+rows)),'position',[7.5,-600])
    set(gca,'FontName','Arial','FontSize',14)
end

ave_force=mean(maxf,2);
eb=std(stress,0,2);

fig2_date_index={[1 2 3] [4 5 6] [7 8 9 ] [11 12]};
draw_from_index(ave_xys,group_name,fig2_date_index,2,4)
draw_from_index_bar(ave_force,eb,group_name,fig2_date_index,2,4)

fig3_date_index={[1 4 7] [2 5 8] [3 6 9 11 12 ]};
draw_from_index(ave_xys,group_name,fig3_date_index,2,4)
draw_from_index_bar(ave_force,eb,group_name,fig3_date_index,2,4)


function draw_from_index(fig_date_xys,fig_titles,fig_date_index,m,n)
fig_size=length(fig_date_index);
h=figure();

for i=1:fig_size
    subplot(m,n,i)

    hold on
    for j=1:length(fig_date_index{i})
        subfig_index=fig_date_index{i};
        index_num=subfig_index(j);
        x=fig_date_xys{index_num,1};
        y=fig_date_xys{index_num,2};
        plot(x,y,LineWidth=1.5)
        legend_title=fig_titles(subfig_index);
    end
    xlabel("Displacement (mm)")
    ylabel("Force (N)")
    xlim([0 15])
    ylim([0,1200])
    legend(legend_title,Location="southeast",Box="off");
    box('on');
    title(sprintf('(%s)',char(96+i)),'position',[7.5,-300])
    set(gca,'FontName','Arial','FontSize',14)
end
end

function draw_from_index_bar(bar_date,ebbar_date,fig_titles,fig_date_index,m,n)
fig_size=length(fig_date_index);
h=figure();

for i=1:fig_size
    subplot(m,n,i)

    hold on

    subfig_index=fig_date_index{i};
    x=categorical(fig_titles(subfig_index));
    x=reordercats(x,fig_titles(subfig_index));
    y=bar_date(subfig_index);
    eb=ebbar_date(subfig_index);
    bar_fig= bar(x,y,1.5*length(y)/10,FaceColor="flat");
    errorbar(x,y,eb,'LineStyle','none',LineWidth=1.5);
    %         legend_title=fig_titles(subfig_index);

    %     xlabel("Type")
    ylabel("Force (N)")
    %     xlim([0 60])
    ylim([0,1300])
    %     legend(legend_title,Location="southeast",Box="off");
    date=y;
    colorrgb=[77,133,189;247,144,61;89,169,90;77,133,189;247,144,61;89,169,90]/255;
    for num=1:size(date)
        bar_fig.CData(num,:)=colorrgb(num,:);
        text_date=string(round(date(num),2));
        text(num,date(num)+eb(num)+0.1,text_date,'VerticalAlignment','bottom','HorizontalAlignment','center')

    end
    box('on');
    title(sprintf('(%s)',char(96+i)),Position=[length(y)/2+0.5 -200])
    set(gca,'FontName','Arial','FontSize',14)
end
end