% Display the geometry of the gt170 gyrotron gun for the first prototype.
% The magnetic field is the one of the asg magnet
% This also creates the h5 files of the geometry and of the magnetic field
% to be used in FENNECS

%% Define and calculate the magnetic field
z=linspace(-0.05,0.2,100);
r=linspace(0.02,0.09,100);
[Z,R] = meshgrid(z,r);

% I = [ 9.7 7.7 96.2499 91.0124 87.4723]; 
% I = [0  0  6.79541  1.967  88.2799];
% I = [0  0  5  1.967  88.2799];
% Current in coils
I = 0.6*[61.56  66.45  113.43  108.09 ]; 
res=zeros([size(Z),3]);


for i=1:size(R,1)
res(i,:,:) = B_Ellip_Cryogenic_170({'aphi','bz','br'},'cryogenic',I,R(i,:),Z(i,:));
end
Aphi=res(:,:,1);
Bz=res(:,:,2);
Br=res(:,:,3);


%% Define the individual boundaries
%read from input file
geom=importproto('geometry_proto_modif.data',[2,Inf]);
geomcells={};
j=1;
i=1;
n=1;
% clean and separate the curves
while i<=size(geom.Z,1)
    if isnan(geom.Z(i))
        j=j+1;
        while isnan(geom.Z(i))
            i=i+1;
        end
        n=1;
    end
    geomcells{j}.Z(n)=geom.Z(i);
    geomcells{j}.R(n)=geom.R(i);
    i=i+1;
    n=n+1;
end

geo_old = geomcells{2};
%% Correcting the corners

% Body -----------------------------------------
N = 30;

% First corner (lone corner)
start1 = 32;
end1 = 36;
R_corner = cat(2,linspace(geomcells{2}.R(start1),0.0869994,30),linspace(0.0869994,geomcells{2}.R(end1),30));
Z_corner = cat(2,linspace(geomcells{2}.Z(start1),0.104844,30),linspace(0.104844,geomcells{2}.Z(end1),30));
% Incomplete square
start2 =510;
R_square = cat(2,linspace(geomcells{2}.R(start2),0.0819994,N),linspace(0.0819994,0.0899991,N),linspace(0.08999911,0.0900011,N),linspace(0.0899911,0.085,N));
Z_square = cat(2,linspace(geomcells{2}.Z(start2),0.0143061,N),linspace(0.0143061,0.0143061,N),linspace(0.0143061,0.00316498,N),linspace(0.00316498,0.00330665,N));

geomcells{2}.R = cat(2,geomcells{2}.R(1:start1),R_corner,geomcells{2}.R(end1:start2),R_square);
geomcells{2}.Z = cat(2,geomcells{2}.Z(1:start1),Z_corner,geomcells{2}.Z(end1:start2),Z_square);

% Cathode --------------------------------------

%Incomplete square
end1 = 320;
R_square = cat(2,linspace(0.085,0.0900011,N),linspace(0.0900011,0.0899991,N),linspace(0.089991,0.08,N));
Z_square = cat(2,linspace(-0.0829438,-0.0829483,N),linspace(-0.0829483,-0.0939433,N),linspace(-0.0939433,-0.094,N));

geomcells{1}.R = cat(2,R_square,geomcells{1}.R(end1:end));
geomcells{1}.Z = cat(2,Z_square,geomcells{1}.Z(end1:end));

%% Define the ceramic boundaries

geomcells{4}.Z = cat(2,geomcells{2}.Z,linspace(geomcells{2}.Z(end),-0.08294,50),geomcells{1}.Z(1:668),linspace(geomcells{1}.Z(668),geomcells{3}.Z(4),50),geomcells{3}.Z(4:end));
geomcells{4}.R = cat(2,geomcells{2}.R,0.085*ones(1,50),geomcells{1}.R(1:668),0.085*ones(1,50),geomcells{3}.R(4:end));

% create the geometry structure
geomcells{1}.name='Cathode';
geomcells{1}.Dval=-30000;
geomcells{2}.name='Body';
geomcells{2}.Dval=0;
geomcells{3}.name='Coaxial insert';
geomcells{3}.Dval=0;
geomcells{4}.name='Ceramics';
geomcells{4}.Dval=0;
% geomcells{4}.name='insulator left';
% geomcells{4}.Dval=0;
% geomcells{5}.name='insulator right';
% geomcells{5}.Dval=0;

for k=1:(length(geomcells)-1)
geomcells{k}.order=2;
geomcells{k}.dim=2;
geomcells{k}.epsce=1e-9;
geomcells{k}.epsge=1e-9;
geomcells{k}.type=0;
geomcells{k}.periodic=0;
end
% geomcells{end-1}.type=2;
% geomcells{end}.type=2;

k = length(geomcells);
geomcells{k}.order=2;
geomcells{k}.dim=2;
geomcells{k}.epsce=1e-9;
geomcells{k}.epsge=1e-9;
geomcells{k}.type=2;
geomcells{k}.periodic=0;


%% Plots
f=figure;
for k=1:length(geomcells)
    plothandle=plot(geomcells{k}.Z, geomcells{k}.R,'k-','linewidth',1.5);
    hold on
    geomcells{k}.points=[geomcells{k}.Z; geomcells{k}.R];
    order=geomcells{k}.order;
    knots=linspace(0,1,length(geomcells{k}.Z)-(order-2));
    knots=augknt(knots, order);
    sizec=size(geomcells{k}.Z);
    order=length(knots)-sizec(end);
    coeffs=[geomcells{k}.Z; geomcells{k}.R];
    pp=spmak(knots,coeffs);
    s=linspace(0,1,10000);
    fittedpos=fnval(pp,s);
    plot(fittedpos(1,:),fittedpos(2,:),'x-')
end
hold on
%end
%axis equal


%Plot the magnetic field lines
%[~,cont1]=contour(Z,R,R.*Aphi,15,'r:','linewidth',3);

% legend([plothandle,cont1],{'Gun geometry', 'Magnetic field lines','Wall approximation'},'location','southwest')
f.PaperUnits='centimeters';
f.PaperSize=[12,8];
xlabel('z [m]')
ylabel('r [m]')
print(f,'phiBprofile_protoasg','-dpdf','-fillpage')
savefig(f,'phiBprofile_protoasg')
hold off


%% Save magnetic field and geometry to disk
save=true;
overwrite=true;
if save
% idr=1:1:length(magnet.r);
% idz=1:1:length(magnet.z);
% Aphi=zeros(length(idr),length(idz));
% Bz=zeros(length(idr),length(idz));
% Br=zeros(length(idr),length(idz));
% 
% for i=1:size(magnet.subcoils,1)
%     Aphi=Aphi+Icoil(i)*magnet.subcoils{i,1}(idr,idz);
%     Bz=Bz+Icoil(i)*magnet.subcoils{i,2}(idr,idz);
%     Br=Br+Icoil(i)*magnet.subcoils{i,3}(idr,idz);
% end
%     savemagtoh5([name,'.h5'],magnet.r,magnet.z,Aphi,Br,Bz,overwrite);
    savegeomtoh5('geom_crmcs.h5',geomcells,1e-2,overwrite);
end