diff --git a/geometries/experiment_upper/exp_outer_withvessel.m b/geometries/experiment_upper/exp_outer_withvessel.m index 88e2f68..8ca38af 100644 --- a/geometries/experiment_upper/exp_outer_withvessel.m +++ b/geometries/experiment_upper/exp_outer_withvessel.m @@ -1,170 +1,170 @@ %% Create the geometry 1 for the T-REX experiment % uses the upper ellipse and constant radius center column. % bottom axial position zbottom=305; % mm radialindent=8; %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=50; %mm lw=linspace(40,rmax,8); rw=linspace(zbottom,480,20); uw=flip(linspace(0,rmax,30)); geomcells(splid).points=flip([zbottom*ones(size(lw')) lw'; rw' rmax*ones(size(rw')); 480*ones(size(uw')) uw']/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(10-radialindent,10,10); rw=linspace(zbottom+30.5,429.5,60); uw=flip(linspace(0,10,10)); %ellipse= lowerring=linspace(zbottom,zbottom+29.5,40); geomcells(splid).points=[ lowerring' (10-radialindent)*ones(size(lowerring')); (zbottom+30)*ones(size(lw')) lw'; rw' 10*ones(size(rw')); 430*ones(size(uw')) uw']/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-radialindent,24+radialindent,50); geomcells(splid).points=flip([ zbottom*ones(size(lw')) lw';]/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.5,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,180+250,60)); geomcells(splid).points=flip([lowerbase' (24+radialindent)*ones(size(lowerbase')); zbottom+30*ones(size(lbaser')) lbaser'; base' 24*ones(size(base))'; tiltedz' tiltedr'; (180+250)*ones(size(axialend')) axialend'; upperend' 40*ones(size(upperend'))]/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,10)); +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; %% 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'; 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.h5',geomcells,1e-2,true); % print(f,name,'-dpdf','-fillpage') % savefig(f,name) % set(f, 'Color', 'w'); % export_fig(f,name,'-eps') hold off \ No newline at end of file diff --git a/geometries/gt170_proto/Draw_B_Ellip_170_geom_proto.m b/geometries/gt170_proto/Draw_B_Ellip_170_geom_proto.m index 612ceac..a9e7b9c 100644 --- a/geometries/gt170_proto/Draw_B_Ellip_170_geom_proto.m +++ b/geometries/gt170_proto/Draw_B_Ellip_170_geom_proto.m @@ -1,170 +1,196 @@ % 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 = 37; +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)); + +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),0.085*ones(1,10),linspace(0.08505,0.0899911,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),0.003165*ones(1,N),linspace(0.003165,0.0029,10),0.0029*ones(1,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); +geomcells{2}.R(284)=[]; +geomcells{2}.Z(284)=[]; +geomcells{2}.R(32)=[]; +geomcells{2}.Z(32)=[]; + % 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)); +R_square = cat(2,linspace(0.090,0.08505,N),0.085*ones(1,10),linspace(0.085,0.0900011,N),linspace(0.0900011,0.0899991,N),linspace(0.089991,0.08,N)); +Z_square = cat(2,-0.082*ones(1,N),linspace(-0.082,-0.0829,10),-0.0829438*ones(1,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)); +% adapt body + +geomcells{3}.R(2:4)=0.085; +geomcells{3}.R=[geomcells{3}.R(1:4) 0.0875 geomcells{3}.R(5:end)]; +geomcells{3}.Z=[geomcells{3}.Z(1:4) geomcells{3}.Z(4) geomcells{3}.Z(5:end)]; +geomcells{3}.Z(1:2)=-0.271; + %% 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)); +geomcells{4}.Z = cat(2,geomcells{2}.Z(1:end-10-N),linspace(geomcells{2}.Z(end-10-N),-0.08294,50),geomcells{1}.Z((1:668)+N+10),linspace(geomcells{1}.Z(668),geomcells{3}.Z(2),50),geomcells{3}.Z(3:end)); +geomcells{4}.R = cat(2,geomcells{2}.R(1:end-10-N),0.085*ones(1,50),geomcells{1}.R((1:668)+N+10),0.085*ones(1,50),geomcells{3}.R(3: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; dist{k}=sqrt(diff(geomcells{k}.Z).^2+diff(geomcells{k}.R).^2); +% remove the problematic points too close to each other +geomcells{k}.Z(dist{k}<1e-5)=[]; +geomcells{k}.R(dist{k}<1e-5)=[]; 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; +dist{k}=sqrt(diff(geomcells{k}.Z).^2+diff(geomcells{k}.R).^2); +% remove the problematic points too close to each other +geomcells{k}.Z(dist{k}<1e-5)=[]; +geomcells{k}.R(dist{k}<1e-5)=[]; %% Plots f=figure; +markers=true; for k=1:length(geomcells) plothandle=plot(geomcells{k}.Z, geomcells{k}.R,'k-','linewidth',1.5); hold on + if markers && k>3 + for j=1:length(geomcells{k}.Z) + text(geomcells{k}.Z(j),geomcells{k}.R(j),sprintf('%i',j),'fontsize',14) + end + end 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); + s=linspace(0,1,100000); 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