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++ Release/SBCmesh.dat 同一文件）
meshFile = fullfile(fileparts(mfilename('fullpath')), cfg.meshDatFile);
mesh = load_SBCmesh_dat(meshFile);mesh.incIndex=findTri(phy.inc,mesh);
mesh.outIndex=findTri(phy.out,mesh);

%矩阵组装
solver.Ai=[];solver.Aj=[];solver.Av=[];
solver.b=zeros(mesh.NbrEdge,1);
solver=assembly_equ(phy,mesh,solver);
%出射处理
solver=assembly_out(phy,mesh,solver);
%入射处理
solver=assembly_inc(phy,mesh,solver);
solver.A=sparse(solver.Ai,solver.Aj,solver.Av);

% 导出 A、b（与 C++ FemType=4 OutFile 格式一致，lda0=2 m）
export_ab_matrix(solver, fullfile(fileparts(mfilename('fullpath')), 'OutFile'));

%求解
solver.x=solver.A\solver.b;

%计算电场
triIndex=findTri(4,mesh);
[x,y,normE]=get_ele2(triIndex,mesh,solver);
patch(x,y,normE,'EdgeAlpha',0)
colormap jet;
colorbar
max(max(normE))

condest(solver.A)
% model=mphopen('SBC.mph');
% [A, b] = extract_mph_Matrice(model, 1);
% full_A=full(A);
% full_A2=full(solver.A);