%This script will plot a given mesh file in .ncdf format generated from the
%acdtool meshconvert routine. The mesh file can be of type 'full', 'halfx',
%'halfy', or 'quarter' if the mesh contains planes of symmetry. In such
%cases the mesh is tiled via reflection across the yz-plane, xz-plane, or
%both for display purposes.
%
%Written by David Bizzozero
%March 31, 2021

mesh_file = 'SRF-GUN-cathR10mm.ncdf';  %Filename for .ncdf mesh
field_file = 'omega3p_results/omega3p.l0.m0000.1.9948170e+08.mod';  %Filename for .mod fields
type = 'halfx';       %Mesh symm. type: 'full', 'halfx', 'halfy', 'quarter'

%Define symmetry plane conditions: xbc is used for type 'quarter' or
%'halfx', ybc is used for type 'quarter' or 'halfy'
xbc = 'magnetic';	%Symmetry plane x=0 boundary condition
ybc = 'magnetic';	%Symmetry plane y=0 boundary condition

%Define plot color function based on E and B field values at vertices.
%Note: E and B are 3 x Nv arrays containing the x,y,z components of the Nv
%number of vertices in the mesh. For example, E(1,:) contains the values of
%Ex at each of the Nv vertices.
fcolor = @(E,B) (E(1,:).^2 + E(2,:).^2 + E(3,:).^2).^0.5;

%% End of user inputs section

%Load interior tetrahedra
t_ext = ncread(mesh_file,'tetrahedron_exterior');
C = ncread(mesh_file,'coords');
X = C(1,:)';    Y = C(2,:)';    Z = C(3,:)';

Tets = t_ext(2:5,:)'+1;
Tris = zeros(4*size(Tets,1),3);

%Load fields
E = ncread(field_file,'efield');
B = ncread(field_file,'bfield');

%Create 4 triangles from each boundary tetrahedron
for i = 1:size(Tets,1)
    Tris(4*i-3,:) = [Tets(i,1) Tets(i,3) Tets(i,2)];
    Tris(4*i-2,:) = [Tets(i,1) Tets(i,4) Tets(i,3)];
    Tris(4*i-1,:) = [Tets(i,1) Tets(i,2) Tets(i,4)];
    Tris(4*i,:)   = [Tets(i,2) Tets(i,3) Tets(i,4)];
end

%Locate boundary triangles on from boundary tetrahedra (boundary code 6)
Tflag = t_ext(6:9,:)==6;
Tris = Tris(Tflag(:),:);

%Mirror boundary triangles and fields using symmetry if needed
if strcmpi(type,'halfx')
    verts = size(X,1);
    X = [X;X];
    Y = [Y;-Y];
    Z = [Z;Z];
    if strcmpi(xbc,'magnetic')
        E = [E(1,:), E(1,:);...
             E(2,:),-E(2,:);...
             E(3,:), E(3,:)];
        B = [B(1,:),-B(1,:);...
             B(2,:), B(2,:);...
             B(3,:), B(3,:)];
    else
        E = [E(1,:),-E(1,:);...
             E(2,:), E(2,:);...
             E(3,:), E(3,:)];
        B = [B(1,:), B(1,:);...
             B(2,:),-B(2,:);...
             B(3,:), B(3,:)];
    end
    Tris = [Tris;Tris+verts];
elseif strcmpi(type,'halfy')
    verts = size(X,1);
    X = [X;-X];
    Y = [Y;Y];
    Z = [Z;Z];
    if strcmpi(xbc,'magnetic')
        E = [E(1,:),-E(1,:);...
             E(2,:), E(2,:);...
             E(3,:), E(3,:)];
        B = [B(1,:), B(1,:);...
             B(2,:),-B(2,:);...
             B(3,:), B(3,:)];
    else
        E = [E(1,:), E(1,:);...
             E(2,:),-E(2,:);...
             E(3,:), E(3,:)];
        B = [B(1,:),-B(1,:);...
             B(2,:), B(2,:);...
             B(3,:), B(3,:)];
    end
    Tris = [Tris;Tris+verts];
elseif strcmpi(type,'quarter')
    verts = size(X,1);
    X = [X;X;-X;-X];
    Y = [Y;-Y;-Y;Y];
    Z = [Z;Z;Z;Z];
    if strcmpi(xbc,'magnetic')
        E = [E(1,:), E(1,:),-E(1,:),-E(1,:);...
             E(2,:),-E(2,:),-E(2,:), E(2,:);...
             E(3,:), E(3,:), E(3,:), E(3,:)];
        B = [B(1,:),-B(1,:),-B(1,:), B(1,:);...
             B(2,:), B(2,:),-B(2,:),-B(2,:);...
             B(3,:), B(3,:), B(3,:), B(3,:)];
    else
        E = [E(1,:),-E(1,:),-E(1,:), E(1,:);...
             E(2,:), E(2,:),-E(2,:),-E(2,:);...
             E(3,:), E(3,:), E(3,:), E(3,:)];
        B = [B(1,:), B(1,:),-B(1,:),-B(1,:);...
             B(2,:),-B(2,:),-B(2,:), B(2,:);...
             B(3,:), B(3,:), B(3,:), B(3,:)];
    end
    Tris = [Tris;Tris+verts;Tris+2*verts;Tris+3*verts];
end

%Plot coloring by definition
if exist('fcolor','var')
    F = fcolor(E,B);
else
    F = ones(size(X));
end

%Plot tetrahedra and Cartesian grid near z=0 and view from front
fontsz = 40;
figure('position',get(0,'screensize'))
tm = trisurf(Tris,Z,X,Y,F,...
    'facealpha',1,'edgecolor','k','edgealpha',0.25);
hold on
axis equal
view(0,0)
colormap(parula(1001))
set(gca,'fontsize',fontsz,'fontname','times')
xlabel('$z\,{\rm [m]}$','interpreter','latex','fontsize',fontsz)
ylabel('$x\,{\rm [m]}$','interpreter','latex','fontsize',fontsz)
zlabel('$y\,{\rm [m]}$','interpreter','latex','fontsize',fontsz)