function Rx1=Jetfire1_1(di,Ct,Tf,Tj,alpha,Ma,Mf,RH,P0,P,k,Cd,elta,deltaHc,tau,directory)
% 一共15个独立变量，为了方便输入修改，所有变量存入'Jetfire1_1'excel表，
% dj为孔口直径,m；Ct为燃料空气混合摩尔系数，可以由浓度算出Ct=1*浓度/22.4；Tf为火焰温度,K；Tj为环境温度,K；alpha为燃料空气化学反应中每摩尔燃烧物所需反应物的摩尔数；
% Ma为空气相对分子量,g/mol；Mf为燃料相对分子量,g/mol；RH为相对湿度,%；P0为环境压力,Pa；P为容器内介质压力,Pa；k为燃料绝热指数,k=cp/cv；
% Cd为气体泄漏系数，与孔形状有关，圆形为1.0；Rg为理想气体常数，大小为8.314；elta为热辐射系数，泄压阀压力以上为0.3，以下为0.4；deltaHc为燃料的燃烧热，KJ/Kg；
% tau为大气传输率,若给出tau请直接输入，若没有给出，请输入0，程序可以用已有参数计算得出。
% -----------------------------------------------------------------------------------以上参数需要输入
% r距离热源半径，m；L为火焰长度，m；Pw为大气中水蒸汽压，Pa；m为质量流速，Kg/s；qr为热辐射量，KW/m^2；Fp为视觉因子；r3为对应辐射量计算出的半径值，m；%% 计算模块
% -----------------
dj=0.02;
Ct=0.095;
Tf=2200;
Tj=298.15;
alpha=1;
Ma=29;
Mf=16;
RH=0;
P0=100000;
P=10000000;
k=1.32;
Cd=1;
elta=0.2;
deltaHc=50000;
tau=0.8;
directory='D:\image\';
% ________________________调试参数，解释即可
x=[1.6 4 12.5 25 37.5]';
x0=length(x);
Z1=ones(x0,1);syms r;
% r = sym('r');Ta=Tj;A=pi/4.*dj.^2*Z1;M=Mf*Z1./10^3;T=Tj;Rg=8.314*Z1;dj=dj*Z1;Ct=Ct*Z1;Tf=Tf*Z1;Tj=Tj*Z1;alpha=alpha*Z1;Ma=Ma*Z1;
Mf=Mf*Z1;RH=RH*Z1;P0=P0*Z1;P=P*Z1;k=k*Z1;Cd=Cd*Z1;elta=elta*Z1;deltaHc=deltaHc*Z1;tau=tau*Z1;%变量：dj/Ct/Tf/Tj/alpha/Ma/Mf;
L=dj.*5.3./Ct.*sqrt(Tf./Tj./alpha.*(Ct+(1-Ct).*Ma./Mf));%计算火焰长度 %变量：RH
if tau==0
Pw=1-1325.*RH.*exp(14.4114-5328./Ta);
tau=2.02.*(Pw.*r).^-0.09;%计算大气传输率
end%变量：P0/P/k/Cd
P1=P0./P;P2=(2./(k+1)).^(1./(k-1));%P1，P2为中间变量
Y=(P1).^(1./k).*sqrt(1-(P1).^((k-1)./k)).*sqrt(2./(k-1).*((k+1)./2).^((k+1)./(k-1)));
if (P1<=P2.^k)m=Cd.*A.*P.*sqrt(M.*k./Rg./T.*P2.^(k+1));
else m=Y.*Cd.*A.*P.*sqrt(M.*k./Rg./T.*P2.^(k+1));%计算质量流速
endFp=1./(4.*pi.*r.^2);%计算视觉因子%变量：elta./deltaHc
qr=tau.*elta.*m.*deltaHc.*Fp;%计算热辐射强度
qr1=tau.*elta.*m.*deltaHc.*Fp-x;
qr2=vpa(qr,5);for i=1:x0
Z2=vpa(solve(qr1(i),r),4);%求解对应的半径
% Z2 = vpa(solve(qr1(i), r(1)), 4);
%qr1(i): 这是一个QR分解的结果，其中qr1是一个函数或变量，i是传递给该函数或变量的参数。QR分解是一种将矩阵分解为正交矩阵Q和上三角矩阵R的方法。
%solve函数用于求解线性方程组。在这里，它将QR分解的结果作为第一个参数，而r(1)可能是另一个矩阵或向量，表示线性方程组的右侧部分。
%vpa函数用于设置数值精度。在这里，它将solve函数的结果保留到小数点后四位有效数字。
r2(i)=double(Z2(find(Z2>0)));
end
r3=double(vpa(r2',4));
Lx=L(1);%% 离散取值模块for i=1:length(r2)
X1=linspace(-r2(i)-L(i)/2,-L(i)/2,100);
Y1=sqrt(r2(i).^2-(X1+L(i)/2).^2);
XX1(:,i)=[X1,X1];YY1(:,i)=[Y1,-Y1];X2=linspace(-L(i)/2,L(i)/2,100);
Y2=r2(i).*ones(1,length(X2));
XX2(:,i)=[X2,X2];YY2(:,i)=[Y2,-Y2];X3=linspace(L(i)/2,r2(i)+L(i)/2,100);
Y3=sqrt(r2(i).^2-(X3-L(i)/2).^2);
XX3(:,i)=[X3,X3];YY3(:,i)=[Y3,-Y3];
ZZ(:,i)=x(i).*ones(6*length(X1),1);
endp=6*length(x)*length(X1);
Yp=[XX1;XX2;XX3];Xp=[YY1;YY2;YY3];
x1=double(reshape(Xp,p,1));y1=double(reshape(Yp,p,1));z1=double(reshape(ZZ,p,1));%% 绘图模块
warning('off')%关闭警告信息。在MATLAB中，有时会出现一些警告信息，使用这个命令可以关闭这些警告，使输出更加简洁。
[Xq,Yq,Zq]=griddata(x1,y1,z1,linspace(min(x1),max(x1))',linspace(min(y1),max(y1)),'v4');%插值
% 这是调用MATLAB的griddata函数来进行插值操作。griddata函数接受散列数据点（x1, y1, z1）和网格点的坐标范围（通过linspace函数生成），并返回插值后的网格数据（Xq, Yq, Zq）。'v4'参数表示使用双三次插值方法。
figure('visible','off');%创建一个不可见的图形窗口。这通常用于在后台执行绘图操作，而不显示图形界面。
set(gcf,'color',[1 1 1]);%设置当前图形窗口的背景颜色为白色（RGB值为[1 1 1]）。gcf表示获取当前图形窗口的句柄。
mycolor=[0.6784 0.9216 1.00000.5088 0.9412 1.00000.3392 0.9608 1.00000.1696 0.9804 1.00000 1.0000 1.00000 1.0000 0.93330 1.0000 0.86670 1.0000 0.80000 1.0000 0.73330 1.0000 0.66670 1.0000 0.60000 1.0000 0.53330 1.0000 0.46670 1.0000 0.40000 1.0000 0.33330 1.0000 0.26670 1.0000 0.20000 1.0000 0.13330 1.0000 0.06670 1.0000 00.0455 1.0000 00.0909 1.0000 00.1364 1.0000 00.1818 1.0000 00.2273 1.0000 00.2727 1.0000 00.3182 1.0000 00.3636 1.0000 00.4091 1.0000 00.4545 1.0000 00.5000 1.0000 00.5455 1.0000 00.5909 1.0000 00.6364 1.0000 00.6818 1.0000 00.7273 1.0000 00.7727 1.0000 00.8182 1.0000 00.8636 1.0000 00.9091 1.0000 00.9545 1.0000 01.0000 1.0000 01.0000 0.9545 01.0000 0.9091 01.0000 0.8636 01.0000 0.8182 01.0000 0.7727 01.0000 0.7273 01.0000 0.6818 01.0000 0.6364 01.0000 0.5909 01.0000 0.5455 01.0000 0.5000 01.0000 0.4545 01.0000 0.4091 01.0000 0.3636 01.0000 0.3182 01.0000 0.2727 01.0000 0.2273 01.0000 0.1818 01.0000 0.1364 01.0000 0.0909 01.0000 0.0455 01.0000 0 0];ax1 = subplot(2,5,[1:3 6:8]);
% 创建一个2x5的子图布局，并在指定位置绘制图形
% plot(x,y);绘制x和y的数据点
plot(x,Y);
contourf(Xq,Yq,Zq,x,'ShowText','on'),xlim([1.1*min(x1),1.1*max(x1)]);ylim([0 1.1*max(y1)]); colorbar%等高线图
grid on;% 显示网格线
xlabel('x/m'),ylabel('y/m');% 设置坐标轴标签
h1=title('垂直喷射火热通量等值线图','fontsize',10);% 设置标题
set(h1,'FontName','微软雅黑');% 设置字体
colormap(ax1,mycolor);% 设置颜色映射ax2 = subplot(2,5,[4:5]);
pcolor(Xq,Yq,Zq);ylim([0 max(y1)]);shading interp,colorbar%伪彩色图
ylabel(colorbar,'热通量/kW/m^2');
h2=title('垂直喷射火热通量插值图','fontsize',10);
set(h2,'FontName','微软雅黑');
colormap(ax2,mycolor);ax3=subplot(2,5,[9:10]);
surfc(Xq,Yq,Zq);ylim([0 max(y1)]);shading interp,colorbar%三维图
view(-20,10);ylabel(colorbar,'热通量/kW/m^2');
h3=title(' 垂直喷射火热通量三维图','fontsize',10);
set(h3,'FontName','微软雅黑');
colormap(ax3,mycolor);file_name=mfilename;
time=datestr(now,'yyyymmddHHMMSSFFF');
ran=rand(1)*10^4;
name=strcat(strcat(file_name,time),num2str(ran));
name(find(name=='.'))=[];
mkdir(directory);
F=getframe(gcf);
imwrite(F.cdata,[directory,name,'.png']);Rx1=char(num2str(roundn(r3',-2)),num2str(roundn(Lx,-2)),name);该函数实现了以下功能：
计算相关物理参数
- 根据输入的孔口直径
dj、燃料空气混合摩尔系数Ct、火焰温度Tf、环境温度Tj、燃料空气化学反应中每摩尔燃烧物所需反应物的摩尔数alpha、空气相对分子量Ma、燃料相对分子量Mf，计算火焰长度L。 - 依据相对湿度
RH（若未给定大气传输率tau时）计算大气中水蒸汽压Pw，进而计算大气传输率tau。 - 根据环境压力
P0、容器内介质压力P、燃料绝热指数k、气体泄漏系数Cd等参数计算质量流速m。 - 基于计算出的质量流速
m、热辐射系数elta、燃料的燃烧热deltaHc等计算热辐射强度qr，并通过求解热辐射强度与给定值（x）的方程得到对应半径值r2（进而得到r3）。
离散取值
对于每个计算出的半径值r2，在特定区间内生成一系列离散的坐标点（X1、Y1、X2、Y2、X3、Y3等），并组合这些坐标点和对应热通量值（ZZ），然后将这些数据进行重塑，以便后续绘图使用。
绘图
- 绘制垂直喷射火热通量的等高线图，包括设置图形窗口不可见、背景颜色为白色，自定义颜色映射，添加坐标轴标签、标题（字体为微软雅黑），显示网格线、颜色条等。
- 绘制垂直喷射火热通量的伪彩色图，设置相关绘图参数，如
ylim、颜色插值、颜色条标签等，并添加标题（字体为微软雅黑）和设置颜色映射。 - 绘制垂直喷射火热通量的三维图，设置视角、
ylim、颜色插值、颜色条标签等，并添加标题（字体为微软雅黑）和设置颜色映射。
结果保存
以当前函数名、时间和随机数生成文件名，创建指定目录，获取当前图形窗口内容并保存为png图像文件，最后将处理后的半径值r3、火焰长度Lx和文件名组合成字符串返回。
总体而言，该函数实现了一个与垂直喷射火相关物理参数计算、数据可视化及结果保存的完整流程，适用于相关领域的模拟、分析和结果展示等任务。