%% 执行MODTRAN函数
% 将输出的tape7文件保存在destFile中
function total_radiance = execute_mod()
%%
currentpath = pwd;
MODPath = 'D:\PcModWin5\Bin';
if(~exist(MODPath,'dir'))
    MODPath = 'E:\PcModWin5\Bin';
end
cd(MODPath);
sourceFile = fullfile(MODPath, 'tape7'); % 待输出的文件： tape7
delete(sourceFile)
system('Mod5.2.1.0.exe > NUL 2>&1'); % 不输出命令行信息
% system('Mod5.2.1.0.exe');

cd(currentpath);
% fprintf('已成功输出tape7文件\n');
%% 读取MODOUT文件
fid = fopen(sourceFile, 'r');
if fid < 0
    warndlg('打开文件失败!');
    return;
end
% 跳过11行
for i = 1:11
    fgetl(fid);
end
% 读取数据
data = fscanf(fid, '%f', [16, inf])';  % 按照一行16个数据的格式读取
fclose(fid); % 关闭文件
if(isempty(data))
    fprintf("tape7数据为空！\n");
    return;
end
data = data(any(data, 2), :);  % 删除非零行
% disp('文件读取成功！');
%% 转换total radiance单位
wavenumber = data(:, 1);
original_TR = data(:, 10); % MODTRAN运行单位结果：W/cm2/sr/cm-1
% changed_TR = original_TR .* 1e7; % 单位：mW/m2/str/cm-1
%% 使用原始单位输出
total_radiance = [wavenumber, original_TR]; % cm-1,  W/(cm2 sr cm-1)
%% 需要单位为：mW/m2/str/cm-1时
% total_radiance = [wavenumber, changed_TR];
%% 需要单位为：mW/m2/str/nm时
% wavelength_nm = 10e6 ./ wavenumber;
% total_radiance = wavenumber .* changed_TR ./ wavelength_nm; %单位：mW/m2/str/nm
% total_radiance = [wavelength_nm,total_radiance];
% total_radiance = flip(total_radiance);
% figure
% plot(total_radiance(:,1),total_radiance(:,2));
%
% title('MODTRAN')
% legend('XCO2 = 420ppm')
% xlabel('wavelength(nm)')
% ylabel('radiance(mW/m2/str/nm)')
end