function export_comsol_normE_from_dat(varargin)
%EXPORT_COMSOL_NORME_FROM_DAT  按 dat 中 Vertex 顺序从 COMSOL 插值 normE
%
%   export_comsol_normE_from_dat
%   export_comsol_normE_from_dat('normExpr', 'ewfd.Ex*conj(ewfd.Ex)+...')
%   export_comsol_normE_from_dat('exExpr', 'ewfd.Ex', 'eyExpr', 'ewfd.Ey', 'ezExpr', 'ewfd.Ez')
%
% 若 normExpr 未指定或为空，将自动探测 Ex/Ey/Ez 并组合成 normE。
% 若模型没有预定义 normE 变量（emw.normE 报错即属此类），请用自动探测或手动指定分量。

opts = parseOpts(varargin);

datPath = opts.datFile;
if ~isfile(datPath)
    error('找不到网格文件: %s（请先 getMesh + export 生成 SBCmesh.dat）', datPath);
end
if ~isfile(opts.mphFile)
    error('找不到 COMSOL 模型: %s', opts.mphFile);
end

[V, nbr] = load_dat_vertices(datPath);
fprintf('已从 %s 读取 %d 个顶点（dat 顺序）\n', datPath, nbr);

model = mphload(opts.mphFile);
coord = V.';  % mphinterp: 3×N

probeCoord = mean(V, 1).';
if opts.autoDetect
    info = resolve_comsol_e_field(model, probeCoord);
    if isempty(opts.exExpr)
        opts.exExpr = info.exExpr;
        opts.eyExpr = info.eyExpr;
        opts.ezExpr = info.ezExpr;
    end
    if isempty(opts.normExpr)
        opts.normExpr = info.normExpr;
    end
end

if isempty(opts.normExpr)
    if isempty(opts.exExpr)
        error(['未指定 normExpr，且自动探测失败。\n', ...
            '请先运行 probe_comsol_fields，或手动传入 ex/ey/ez 表达式。']);
    end
    opts.normExpr = sprintf('sqrt(%s*conj(%s)+%s*conj(%s)+%s*conj(%s))', ...
        opts.exExpr, opts.exExpr, opts.eyExpr, opts.eyExpr, opts.ezExpr, opts.ezExpr);
end

fprintf('正在插值 normE 表达式:\n  %s\n', opts.normExpr);
normE = mphinterp(model, opts.normExpr, 'coord', coord);
normE = real(normE(:));

if numel(normE) ~= nbr
    error('插值结果长度 %d 与顶点数 %d 不一致', numel(normE), nbr);
end

nanMask = isnan(normE) | isinf(normE);
if any(nanMask)
    warning('export_comsol_normE:NaN', ...
        '%d 个顶点插值为 NaN/Inf（可能落在网格外或未求解）', nnz(nanMask));
end

if opts.exportComponents
    fprintf('正在插值 Ex/Ey/Ez ...\n');
    Ex = mphinterp(model, opts.exExpr, 'coord', coord);
    Ey = mphinterp(model, opts.eyExpr, 'coord', coord);
    Ez = mphinterp(model, opts.ezExpr, 'coord', coord);
end

if ~exist(opts.outDir, 'dir')
    mkdir(opts.outDir);
end

writeNormE(fullfile(opts.outDir, 'normE'), normE);

if opts.exportComponents
    writeComplexField(fullfile(opts.outDir, 'Ex'), Ex);
    writeComplexField(fullfile(opts.outDir, 'Ey'), Ey);
    writeComplexField(fullfile(opts.outDir, 'Ez'), Ez);
end

if opts.exportTable
    T = [(1:nbr).', V, normE];
    writematrix(T, fullfile(opts.outDir, 'vertex_normE.dat'), 'Delimiter', 'tab');
end

metaFile = fullfile(opts.outDir, 'comsol_expr.txt');
fid = fopen(metaFile, 'w');
if fid > 0
    fprintf(fid, 'normExpr=%s\n', opts.normExpr);
    fprintf(fid, 'exExpr=%s\n', opts.exExpr);
    fprintf(fid, 'eyExpr=%s\n', opts.eyExpr);
    fprintf(fid, 'ezExpr=%s\n', opts.ezExpr);
    fclose(fid);
end

fprintf('已写出 %s\n', fullfile(opts.outDir, 'normE'));
fprintf('  max(normE) = %.12g\n', max(normE(~nanMask)));
fprintf('  行顺序与 %s 的 Vertex 第 1..%d 行一一对应\n', datPath, nbr);
end

function opts = parseOpts(args)
p = inputParser;
p.addParameter('mphFile', 'SBC.mph', @ischar);
p.addParameter('datFile', 'SBCmesh.dat', @ischar);
p.addParameter('normExpr', '', @ischar);
p.addParameter('exExpr', '', @ischar);
p.addParameter('eyExpr', '', @ischar);
p.addParameter('ezExpr', '', @ischar);
p.addParameter('outDir', 'OutFile_comsol', @ischar);
p.addParameter('exportComponents', false, @islogical);
p.addParameter('exportTable', true, @islogical);
p.addParameter('autoDetect', true, @islogical);
p.parse(args{:});
opts = p.Results;
end

function writeNormE(path, normE)
fid = fopen(path, 'w');
if fid < 0
    error('无法写入: %s', path);
end
c = onCleanup(@() fclose(fid));
fprintf(fid, '//\n');
fprintf(fid, '%.12g\n', normE);
end

function writeComplexField(path, field)
fid = fopen(path, 'w');
if fid < 0
    error('无法写入: %s', path);
end
c = onCleanup(@() fclose(fid));
fprintf(fid, '//\n');
fprintf(fid, '%.12g %.12g\n', [real(field(:)), imag(field(:))].');
end