% Creates the geometry for a test case of an empty ring with sqare
% cross-section. The boundaries are defined using the curve b-splines
% representation.
% The magnetic field is the one created by the asg magnet 
% Both the magnet and geometry are definde in the Final report of the 
% Development of the european gyrotron CCCGDS6 


%% Define the individual boundaries
rmin=0.015;
rmax=0.065;
zmin=-0.025;
zmax=0.025;


npts=5;
geomcells={};
geomcells{1}.Z=zeros(1,4*npts-4);
geomcells{1}.R=zeros(1,4*npts-4);

geomcells{1}.Z(1)=zmin;
geomcells{1}.R(1)=rmin;

i=1;
dz=(zmax-zmin)/(npts-1);
dr=(rmax-rmin)/4;
geomcells{1}.Z(i+(1:npts-1))=zmin+dz*(1:(npts-1));
geomcells{1}.R(i+(1:npts-1))=rmin+[dr/2 dr dr/2 0];

i=i+npts-1;
dr=(rmax-rmin)/(npts-1);
dz=(zmax-zmin)/4;
geomcells{1}.R(i+(1:npts-1))=rmin+dr*(1:(npts-1));
geomcells{1}.Z(i+(1:npts-1))=zmax+[-dz/2 -dz -dz/2 0];

i=i+npts-1;
dz=(zmax-zmin)/(npts-1);
dr=(rmax-rmin)/4;
geomcells{1}.Z(i+(1:npts-1))=zmax-dz*(1:(npts-1));
geomcells{1}.R(i+(1:npts-1))=rmax-[dr/2 dr dr/2 0];

i=i+npts-1;
dr=(rmax-rmin)/(npts-1);
dz=(zmax-zmin)/4;
geomcells{1}.R(i+(1:npts-2))=rmax-dr*(1:(npts-2));
geomcells{1}.Z(i+(1:npts-2))=zmin+[dz/2 dz dz/2];

    
% geomcells{1}.Z=flip(geomcells{1}.Z);
% geomcells{1}.R=flip(geomcells{1}.R);


% geomcells{1}.Z=[geomcells{1}.Z(5:end) geomcells{1}.Z(2:5)];
% geomcells{1}.R=[geomcells{1}.R(5:end) geomcells{1}.R(2:5)];
for k=1:length(geomcells)
geomcells{k}.order=3;
geomcells{k}.dim=2;
geomcells{k}.epsce=1e-9;
geomcells{k}.epsge=1e-9;
geomcells{k}.Dval=0;
geomcells{k}.name=sprintf('Electrode_%i',k);
geomcells{k}.Z=geomcells{k}.Z;
geomcells{k}.R=geomcells{k}.R;
geomcells{k}.points=[geomcells{k}.Z; geomcells{k}.R];
geomcells{k}.type=0;
geomcells{k}.periodic=1;
end

%% Plots
f=figure;
for k=1:length(geomcells)
    plothandle=plot(geomcells{k}.Z, geomcells{k}.R,'k-x','linewidth',1.5);
    hold on
    order=geomcells{k}.order;
    
    sizec=size(geomcells{k}.Z);
    
    coeffs=[[geomcells{k}.Z(1:end-1) geomcells{k}.Z(1:geomcells{k}.order)]; geomcells{k}.R(1:end-1) geomcells{k}.R(1:geomcells{k}.order)];
    
    if ~geomcells{k}.periodic
        knots=linspace(0,1,length(geomcells{k}.Z)-(order-2));
        knots=augknt(knots, order);
    else
        knots=linspace(0,1,size(coeffs,2)-(order-2));
        dknots=knots(2)-knots(1);
        knots=[knots(1)+dknots*(-order+1:-1) knots knots(end)+dknots*(1:order-1)];
    end
    pp=spmak(knots,coeffs);
    %pp=csape(geomcells{k}.Z,geomcells{k}.R,'periodic');
    
    s=linspace(0,1,1000);
    fittedpos=fnval(pp,s);
    plot(fittedpos(1,:),fittedpos(2,:),'x-')
end
%axis equal
%[~,cont2]=contour(Zphi,Rphi,Phi,20,'b');
%rectangle('Position',[-0.011, 0.06375, 0.032+0.011, 0.081-0.06375],'EdgeColor','magenta','Linestyle','--')

legend([plothandle],{'Gun geometry', 'Magnetic field lines'},'location','southwest')
f.PaperUnits='centimeters';
f.PaperSize=[12,8];
xlabel('z [m]')
ylabel('r [m]')


% print(f,name,'-dpdf','-fillpage')
% savefig(f,name)
% set(f, 'Color', 'w');
% export_fig(f,name,'-eps')
hold off


%% Save magnetic field and geometry to disk
save=true;
overwrite=true;
if save
    savegeomtoh5('test_square_inb2.h5',geomcells,1e-2,overwrite);
end