clc
clear all

%参考四面体
%顶点1-[1,0,0] 2-[0,1,0] 3-[0,0,1] 4-[0,0,1]

% 算例配置（改波长/材料请编辑 case_config.m）
[phy, cfg] = case_config();

%读取网格（与 C++ SBCmesh.dat 一致；无 .mat 时读 .dat）
rootDir = fileparts(mfilename('fullpath'));
matFile = fullfile(rootDir, cfg.meshMatFile);
datFile = fullfile(rootDir, cfg.meshDatFile);
if isfile(matFile)
    load(matFile, 'mesh');
elseif isfile(datFile)
    mesh = load_SBCmesh_dat(datFile);
else
    altDat = fullfile(rootDir, '..', '..', '..', '3D opticsfem-master', 'SBCmesh.dat');
    if isfile(altDat)
        mesh = load_SBCmesh_dat(altDat);
    else
        error('未找到 SBCmesh.mat 或 SBCmesh.dat，请先运行 getMesh.m 或复制 C++ 网格。');
    end
end

%矩阵组装
solver.Ai=[];solver.Aj=[];solver.Av=[];
solver.Bi=[];solver.Bj=[];solver.Bv=[];
solver=assembly_equ(phy,mesh,solver);
solver.A=sparse(solver.Ai,solver.Aj,solver.Av);
solver.B=sparse(solver.Bi,solver.Bj,solver.Bv);

% 导出 A/B 矩阵（与 C++ FemType=5 OutFile 格式一致）
export_ab_matrix(solver);

%求解设置
targetff0=phy.ff0;
targetlam0=phy.c_const/targetff0;
targetk0=2*pi/targetlam0;
sigma=targetk0*targetk0;
numberSolve=cfg.NbrMode;
%求解
[x,GAMA_sq]=eigs(solver.A,solver.B,numberSolve,sigma);
solverk0=diag(sqrt(GAMA_sq));
solverlam0=2*pi./solverk0;
solverff0=phy.c_const./solverlam0


%计算电场
Ex=zeros(mesh.NbrVertex,numberSolve);
Ey=zeros(mesh.NbrVertex,numberSolve);
Ez=zeros(mesh.NbrVertex,numberSolve);
normE=zeros(mesh.NbrVertex,numberSolve);
for i=1:numberSolve
    xx=x(:,i);
    [Ex(:,i),Ey(:,i),Ez(:,i),normE(:,i)]=get_ele(mesh,xx);
end

% 导出 freq / Ex / Ey / Ez / normE（与 C++ Post_3D_EigenFreq 格式一致）
export_eigen_fields(solverff0, Ex, Ey, Ez, fullfile(fileparts(mfilename('fullpath')), 'OutFile'));

%绘图
triIndex=findTri(3,mesh);
tri=mesh.Tri(triIndex,:);
plotE(mesh.Vertex,tri,normE(:,3),0);