%COMPARE_NORME_ZSLICE  Compare |E| on z=const plane: COMSOL vs OpticsFEM.
%
% Usage (in MATLAB):
%   cd('.../3D opticsfem-master/compare')
%   compare_normE_zslice
%
% Before first run, export COMSOL slice once:
%   run_export_comsol        % 一键导出（需 LiveLink + SBC.mph）
%
% Or manually export from COMSOL GUI to data/comsol_z050.txt
% with columns: x  y  z  normE  (tab/space separated).

clear; clc; close all;

thisDir = fileparts(mfilename('fullpath'));
projDir = fileparts(thisDir);

cfg = struct();
cfg.zCut = 0.5;
cfg.meshFile = fullfile(projDir, 'build', 'Release', 'SBCmesh.dat');
cfg.outFileDir = fullfile(projDir, 'build', 'Release', 'OutFile');
cfg.comsolFile = fullfile(thisDir, 'data', 'comsol_z050.txt');
cfg.outPng = fullfile(thisDir, 'data', 'compare_normE_z050.png');
cfg.gridN = 121;

if ~isfile(cfg.meshFile)
    cfg.meshFile = fullfile(projDir, 'SBCmesh.dat');
end

assert(isfile(cfg.meshFile), 'Mesh not found: %s', cfg.meshFile);
assert(isfile(fullfile(cfg.outFileDir, 'normE')), ...
    'OpticsFEM normE not found. Run OpticsFEM.exe first.');

if ~isfile(cfg.comsolFile)
    mphPath = fullfile(projDir, '..', '三维matlab代码', 'matlab 3D一阶散射问题', 'SBC.mph');
    if isfile(mphPath) && exist('mphload', 'file') == 2
        fprintf('COMSOL data missing, exporting from mph ...\n');
        export_comsol_slice(mphPath, cfg.comsolFile, cfg.zCut, cfg.gridN, cfg.gridN);
    else
        error(['COMSOL 切面数据不存在: %s\n\n' ...
            '请先任选一种方式生成该文件：\n' ...
            '  1) 运行 run_export_comsol.m（需 LiveLink + SBC.mph）\n' ...
            '  2) COMSOL GUI: Cut Plane z=0.5 → Export → 保存到上述路径\n\n' ...
            '说明：项目里只有 SBC.mph.lock，没有 mph 本体，需在本机 COMSOL 里打开/保存模型。'], ...
            cfg.comsolFile);
    end
end

fprintf('Loading mesh: %s\n', cfg.meshFile);
mesh = load_SBCmesh_dat(cfg.meshFile);

fprintf('Loading OpticsFEM normE: %s\n', cfg.outFileDir);
normE = load_opticsfem_normE(cfg.outFileDir);
assert(numel(normE) == mesh.NbrVertex, ...
    'normE length (%d) != NbrVertex (%d)', numel(normE), mesh.NbrVertex);

fprintf('Building OpticsFEM z=%.2f slice ...\n', cfg.zCut);
sliceOpt = build_z_slice(mesh, normE, cfg.zCut, 3);

fprintf('Loading COMSOL slice: %s\n', cfg.comsolFile);
sliceCom = load_comsol_slice(cfg.comsolFile, cfg.zCut, 1e-3);

xg = linspace(-0.5, 0.5, cfg.gridN);
yg = linspace(-0.5, 0.5, cfg.gridN);
gridCom = interp_slice_to_grid(sliceCom, xg, yg);
gridOpt = interp_slice_to_grid(sliceOpt, xg, yg);

fprintf('Plotting comparison ...\n');
stats = plot_compare_heatmaps(gridCom, gridOpt, cfg.zCut, cfg.outPng);

assignin('base', 'compareStats', stats);
assignin('base', 'sliceOpt', sliceOpt);
assignin('base', 'sliceCom', sliceCom);
fprintf('Done.\n');