%% Create the geometry 2 for the T-REX experiment
% uses the upper ellipse and constant radius center column.
% bottom axial position
zbottom=305; % mm
radialindent=8; %mm
ztop=305+230; %mm
% Outer vacuum vessel
splid=1;
geomcells(splid).Dval=0; % V
geomcells(splid).order=3;
geomcells(splid).dim=2;
geomcells(splid).name='vacuum vessel';
geomcells(splid).type=0;
geomcells(splid).periodic=0;
rmax=60; %mm
lw=linspace(40,rmax,20);
rw=linspace(zbottom,ztop,60);
uw=flip(linspace(0,rmax,30));
geomcells(splid).points=flip([zbottom*ones(size(lw')) lw';
rw' rmax*ones(size(rw'));
ztop*ones(size(uw')) uw';
linspace(ztop,ztop+3,10)' -1*ones(10,1);
(ztop+3)*ones(30,1) linspace(0,rmax+3,30)';
linspace(ztop+3,zbottom-3,60)' (rmax+3)*ones(60,1);
(zbottom-3)*ones(30,1) linspace(rmax+3,40,30)';
linspace(zbottom-3,zbottom,10)' 40*ones(10,1)]/1e3);
% Central electrode
splid=2;
geomcells(splid).Dval=-20000; % V
geomcells(splid).order=3;
geomcells(splid).dim=2;
%splines(1).epsge=1e-6;
%splines(1).epsce=1e-6;
geomcells(splid).name='central electrode';
geomcells(splid).type=0;
geomcells(splid).periodic=0;
lw=linspace(-1,10,15);
rw=linspace(zbottom+30.5,429.5,60);
uw=flip(linspace(-1,10,10));
%ellipse=
lowerring=linspace(zbottom,zbottom+29.5,40);
geomcells(splid).points=[
(zbottom+30)*ones(size(lw')) lw';
rw' 10*ones(size(rw'));
430*ones(size(uw')) uw';
flip(rw') -1*ones(size(rw'))]/1e3;
% Bottom plate
splid=3;
geomcells(splid).Dval=0; % V
geomcells(splid).order=3;
geomcells(splid).dim=2;
%splines(1).epsge=1e-6;
%splines(1).epsce=1e-6;
geomcells(splid).name='vacuum vessel';
geomcells(splid).type=0;
geomcells(splid).periodic=0;
lw=linspace(10,24,40);
geomcells(splid).points=flip([
zbottom*ones(size(lw')) lw';
linspace(zbottom,zbottom-3,10)' 24*ones(10,1);
(zbottom-3)*ones(size(lw')) flip(lw');
linspace(zbottom-3,zbottom,10)' 10*ones(10,1);]/1e3);
% Outer electrode
splid=4;
geomcells(splid).Dval=0; % V
geomcells(splid).order=3;
geomcells(splid).dim=2;
%splines(1).epsge=1e-6;
%splines(1).epsce=1e-6;
geomcells(splid).name='outer electrode';
geomcells(splid).type=0;
geomcells(splid).periodic=0;
lbaser=linspace(24+radialindent,24,20);
base=linspace(zbottom+30,349,80);
lowerbase=linspace(zbottom,zbottom+29.5,50);
%tiltedwall
alpha=0.1745;
tiltedz=linspace(350,180+250,500);
tiltedr=(tiltedz-tiltedz(1))*tan(alpha)+24;
%tilted ellipse
cosa=cos(alpha);
sina=sin(alpha);
r_c=28;
z_c=375;
Lz=25;
Lr=5;
deltar=(tiltedr-r_c);
deltaz=(tiltedz-z_c);
D=((deltaz*cosa+deltar*sina)/Lz).^2+((deltaz*sina-deltar*cosa)/Lr).^2;
w=1-D;
ellipsez=tiltedz(w>=0);
deltaz=ellipsez-z_c;
a=(Lr^2*sina^2+Lz^2*cosa^2);
b=2*deltaz*sina*cosa*(Lr^2-Lz^2);
c=-Lr^2*Lz^2 +deltaz.^2*(Lr^2*cosa^2+Lz^2*sina^2);
tiltedr(w>=0)=r_c+(-b-sqrt(b.^2-4*a.*c))/2/a;
axialend=linspace(tiltedr(end)+0.5,39.5,15);
upperend=flip(linspace(zbottom+30,180+250,60));
geomcells(splid).points=flip([
base' 24*ones(size(base))';
tiltedz' tiltedr';
(180+250)*ones(size(axialend')) axialend';
upperend' 40*ones(size(upperend'));
zbottom+30*ones(size(lbaser')) linspace(40,24,20)']/1e3);
% Dielectric ring
% splid=5;
% geomcells(splid).Dval=0; % V
% geomcells(splid).order=3;
% geomcells(splid).dim=2;
% %splines(1).epsge=1e-6;
% %splines(1).epsce=1e-6;
% geomcells(splid).name='dielectric rings';
% geomcells(splid).type=2;
% geomcells(splid).periodic=0;
% ow=flip(linspace(zbottom,zbottom+30,130));
% lw=flip(linspace(10.1,23.5,30));
% rw=linspace(zbottom,429.5,60);
% uw=flip(linspace(0,10,10));
% geomcells(splid).points=[flip([upperend', 40*ones(size(upperend'))]);
% flip([(180+250)*ones(size(axialend')) axialend']);
% flip([tiltedz' tiltedr']);
% flip([base' 24*ones(size(base))']);
% ow' 24*ones(size(ow'));
% zbottom*ones(size(lw')) lw';
% rw' 10*ones(size(rw'));
% 430*ones(size(uw')) uw']/1e3;
% Dielectric ring inner
splid=5;
geomcells(splid).Dval=0; % V
geomcells(splid).order=2;
geomcells(splid).dim=2;
%splines(1).epsge=1e-6;
%splines(1).epsce=1e-6;
geomcells(splid).name='dielectric ring center';
geomcells(splid).type=2;
geomcells(splid).periodic=0;
ow=flip(linspace(zbottom-3,zbottom+30,130));
lw=flip(linspace(10.1,23.5,30));
rw=linspace(zbottom-3,334.5,50);
uw=flip(linspace(-1,10,10));
geomcells(splid).points=[(zbottom-3)*ones(size(uw')) flip(uw');
rw' 10*ones(size(rw'));
335*ones(size(uw')) uw';
flip(rw') -1*ones(size(rw'))]/1e3;
% Dielectric ring outer
splid=6;
geomcells(splid).Dval=0; % V
geomcells(splid).order=2;
geomcells(splid).dim=2;
%splines(1).epsge=1e-6;
%splines(1).epsce=1e-6;
geomcells(splid).name='dielectric ring center';
geomcells(splid).type=2;
geomcells(splid).periodic=0;
ow=flip(linspace(zbottom,zbottom+30,130));
lw=flip(linspace(10.1,23.5,30));
rw=linspace(zbottom-3,334.5,50);
uw=linspace(40,24,20);
geomcells(splid).points=[rw' uw(1)*ones(size(rw'));
335*ones(size(uw')) uw';
flip(rw') uw(end)*ones(size(rw'));
(zbottom-3)*ones(size(uw')) flip(uw');
]/1e3;
%% Plots
f=figure;
for k=1:length(geomcells)
plothandle=plot(geomcells(k).points(:,1), geomcells(k).points(:,2),'k-x','linewidth',1.5);
hold on
%geomcells(k).points=[geomcells(k).Z; geomcells(k).R];
order=geomcells(k).order;
knots=linspace(0,1,size(geomcells(k).points,1)-(order-2));
knots=augknt(knots, order);
sizec=size(geomcells(k).points,1);
order=length(knots)-sizec(end);
coeffs=geomcells(k).points';
geomcells(k).periodic=0;
pp=spmak(knots,coeffs);
s=linspace(0,1,1000);
fittedpos=fnval(pp,s);
plot(fittedpos(1,:),fittedpos(2,:),'x-')
end
legend(plothandle,{'Gun geometry'},'location','southwest')
f.PaperUnits='centimeters';
f.PaperSize=[12,8];
xlabel('z [m]')
ylabel('r [m]')
axis equal
savegeomtoh5('exp_outer_with_vessel_geom_2.h5',geomcells,1e-2,true);
% print(f,name,'-dpdf','-fillpage')
% savefig(f,name)
% set(f, 'Color', 'w');
% export_fig(f,name,'-eps')
hold off