diff --git a/matlab/DispPressure.m b/matlab/DispPressure.m
index 532e2d4..e65848f 100644
--- a/matlab/DispPressure.m
+++ b/matlab/DispPressure.m
@@ -1,35 +1,35 @@
function DispPressure(M,it)
dens=mean(M.N(:,:,it),3);
maxdens=max(max(dens));
-f=figure();
+f=figure('title',M.name);
plotsubplot(subplot(2,3,1),squeeze(mean(M.Presstens(1,:,:,it),4)),'P_{rr} [Pa]');
plotsubplot(subplot(2,3,2),squeeze(mean(M.Presstens(2,:,:,it),4)),'P_{r\theta} [Pa]');
plotsubplot(subplot(2,3,3),squeeze(mean(M.Presstens(3,:,:,it),4)),'P_{rz} [Pa]');
plotsubplot(subplot(2,3,4),squeeze(mean(M.Presstens(4,:,:,it),4)),'P_{\theta\theta} [Pa]');
plotsubplot(subplot(2,3,5),squeeze(mean(M.Presstens(5,:,:,it),4)),'P_{z\theta} [Pa]');
plotsubplot(subplot(2,3,6),squeeze(mean(M.Presstens(6,:,:,it),4)),'P_{zz} [Pa]');
sgtitle(sprintf('Pressure t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[16 14];
f.PaperSize=papsize;
print(f,sprintf('%sfluid_Pressure',name),'-dpdf','-fillpage')
function plotsubplot(ax,Val,clabel)
h=surface(ax,M.zgrid,M.rgrid,Val);
xlim(ax,[M.zgrid(1) M.zgrid(end)])
rgridend=sum(M.nnr(1:2));
ylim(ax,[M.rgrid(1) M.rgrid(rgridend)])
xlabel(ax,'z [m]')
ylabel(ax,'r [m]')
c = colorbar(ax);
c.Label.String= clabel;
view(ax,2)
set(h,'edgecolor','none');
end
end
diff --git a/matlab/PotentialWell.m b/matlab/PotentialWell.m
new file mode 100644
index 0000000..f7755f6
--- /dev/null
+++ b/matlab/PotentialWell.m
@@ -0,0 +1,35 @@
+function [pot] = PotentialWell(Phi,Atheta,r,z,creategraph)
+%PotentialWell Computes and display the potential well given the electrostatic potentila Phi
+% and the magnetic vector potential Atheta
+if nargin <5
+ creategraph=true;
+end
+[rgrid,~]=meshgrid(r,z);
+rAthet=(rgrid.*Atheta)';
+pot=zeros(size(Phi));
+for j=1:size(pot,2)
+ potmin=pot(:,j);
+ for i=1:size(pot,1)
+ for k=1:size(pot,2)
+ rpos=find(rAthet(i,j)>rAthet(:,k),1,'last');
+ if rpos==size(rAthet,1)
+ continue
+ end
+ if isempty(rpos)
+ continue
+ end
+ rinterp=(rAthet(i,j)-rAthet(rpos,k))./(rAthet(rpos+1,k)-rAthet(rpos,k));
+ potmin(i)=min(potmin(i),pot(rpos,k)+rinterp*(pot(rpos+1,k)-pot(rpos,k)));
+
+ end
+ end
+ pot(:,j)=Phi(:,j)-potmin;
+end
+if(creategraph)
+figure
+surface(z(2:end),r(2:end),-pot(2:end,2:end),'edgecolor','none')
+xlabel('r')
+ylabel('z')
+colorbar
+end
+end
diff --git a/matlab/analyse_espicfields.m b/matlab/analyse_espicfields.m
index fe2014f..9ec8d5a 100644
--- a/matlab/analyse_espicfields.m
+++ b/matlab/analyse_espicfields.m
@@ -1,372 +1,282 @@
% file='teststablegauss_13_traject2.h5';
% file='unigauss.h5';
% file='Dav_5e14fine_wr+.h5';
% file='stablegauss_8e19fine.h5'
% % %file='teststable_Dav1.h5';
% % %close all hidden;
% if (~exist('M','var'))
% M.file=file;
% end
% M=load_espic2d(file,M,'all');
t2dlength=size(M.t2d,1)
-fieldstart=max(1,t2dlength-100);%floor(0.95*t2dlength);
+fieldstart=max(1,t2dlength-500);%floor(0.95*t2dlength);
deltat=t2dlength-fieldstart;
%[~, rgridend]=min(abs(M.rgrid-0.045));
rgridend=sum(M.nnr(1:2));
[~, name, ~] = fileparts(M.file);
M.plotEnergy
%fieldstart=2300; fieldend=2500;
dens=mean(M.N(:,:,fieldstart:end),3);
pot=mean(M.pot(:,:,fieldstart:end),3);
brillratio=2*dens*M.me./(M.eps_0*M.Bz'.^2);
[R,Z]=meshgrid(M.rgrid,M.zgrid);
Rinv=1./R;
Rinv(:,1)=0;
VTHET=mean(M.fluidUTHET(:,:,fieldstart:end),3);
omegare=(VTHET.*Rinv');
maxdens=max(max(dens));
% f=figure();
% ax3=gca;
% surface(ax3,M.zgrid(1:end-1),M.rgrid(1:end-1),(squeeze(mean(M.Presstens(4,:,:,fieldstart:end),4)))*M.vlight)
% xlabel(ax3,'z [m]')
% ylabel(ax3,'r [m]')
% xlim(ax3,[M.zgrid(1) M.zgrid(end)])
% ylim(ax3,[M.rgrid(1) M.rgrid(50)])
% colormap(ax3,'jet')
% c = colorbar(ax3);
% c.Label.String= 'thermal v_\theta [m/s]';
% %c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))];
% %caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))])
% title(ax3,'Azimuthal velocity')
% %set(ax3,'colorscale','log')
% grid(ax3, 'on')
% view(ax3,2)
% f.PaperOrientation='landscape';
% [~, name, ~] = fileparts(M.file);
% f.PaperUnits='centimeters';
% papsize=[16 14];
% f.PaperSize=papsize;
% print(f,sprintf('%sfluid_thermtheta',name),'-dpdf','-fillpage')
%% Psi function
try
- displaypsi(M,deltat)
+ M.displaypsi(deltat)
catch
end
%% Density
f=figure('Name', sprintf('%sfluid_dens',name));
ax1=gca;
h=surface(ax1,M.zgrid,M.rgrid,dens);
-%set(h,'edgecolor','none');
+set(h,'edgecolor','none');
ylim(ax1,[M.rgrid(1) M.rgrid(rgridend)])
xlim(ax1,[M.zgrid(1) M.zgrid(end)])
xlabel(ax1,'z [m]')
ylabel(ax1,'r [m]')
title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens)))
c = colorbar(ax1);
c.Label.String= 'n[m^{-3}]';
view(ax1,2)
grid on;
hold on;
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[16 14];
f.PaperSize=papsize;
print(f,sprintf('%sfluid_dens',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sfluid_dens',name))
%% Fieldswide
f=figure('Name', sprintf('%s fields',name));
ax1=gca;
h=contourf(ax1,M.zgrid,M.rgrid,dens,'Displayname','n_e [m^{-3}]');
hold on;
%[c1,h]=contour(ax1,M.zgrid,M.rgrid,pot./1e3,5,'m--','ShowText','off','linewidth',1.5,'Displayname','Equipotentials [kV]');
%clabel(c1,h,'Color','white')
Blines=zeros(size(M.Athet));
for i=1:length(M.zgrid)
Blines(i,:)=M.Athet(i,:).*M.rgrid';
end
zindex=floor(length(M.zgrid/2));
levels=logspace( log10(min(min(Blines(:,2:end)))), log10(max(max(Blines(:,:)))),8);
contour(ax1,M.zgrid,M.rgrid,Blines',levels,'r-.','linewidth',1.5,'Displayname','Magnetic field lines');
xlim(ax1,[M.zgrid(1) M.zgrid(end)])
ylim(ax1,[M.rgrid(1) M.rgrid(end)])
legend
xlabel(ax1,'z [m]')
ylabel(ax1,'r [m]')
title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens)))
c = colorbar(ax1);
c.Label.String= 'n[m^{-3}]';
view(ax1,2)
%set(h,'edgecolor','none');
grid on;
%rectangle('position',[M.zgrid(5) M.rgrid(3) M.zgrid(end-5)-M.zgrid(5) (M.rgrid(rgridend)-M.rgrid(1))],'Edgecolor','white','linestyle','--','linewidth',2)
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[18 10];
f.PaperSize=papsize;
print(f,sprintf('%sFieldswide',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sFieldswide',name))
%% Fields
f=figure('Name', sprintf('%s fields',name));
ax1=gca;
h=contourf(ax1,M.zgrid,M.rgrid,dens,'Displayname','n_e [m^{-3}]');
hold on;
[c1,h]=contour(ax1,M.zgrid,M.rgrid,pot./1e3,'m--','ShowText','on','linewidth',1.5,'Displayname','Equipotentials [kV]');
clabel(c1,h,'Color','white')
Blines=zeros(size(M.Athet));
for i=1:length(M.zgrid)
Blines(i,:)=M.Athet(i,:).*M.rgrid';
end
zindex=floor(length(M.zgrid/2));
levels=logspace( log10(min(min(Blines(:,2:end)))), log10(max(max(Blines(:,:)))),10);
contour(ax1,M.zgrid,M.rgrid,Blines',levels,'r-.','linewidth',1.5,'Displayname','Magnetic field lines');
xlim(ax1,[M.zgrid(1) M.zgrid(end)])
ylim(ax1,[M.rgrid(1) M.rgrid(rgridend)])
legend
xlabel(ax1,'z [m]')
ylabel(ax1,'r [m]')
title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens)))
c = colorbar(ax1);
c.Label.String= 'n[m^{-3}]';
view(ax1,2)
%set(h,'edgecolor','none');
grid on;
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[18 10];
f.PaperSize=papsize;
print(f,sprintf('%sFields',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sFields',name))
-%% Wells
-Wells(M,pot,[fieldstart t2dlength],maxdens)
-Wells(M,M.pot(:,:,1),1, max(max(M.N(:,:,1))))
-
%% Brillouin
f=figure('Name', sprintf('%s Brillouin',name));
ax1=gca;
h=surface(ax1,M.zgrid,M.rgrid,brillratio);
set(h,'edgecolor','none');
xlim(ax1,[M.zgrid(1) M.zgrid(end)])
ylim(ax1,[M.rgrid(1) M.rgrid(rgridend)])
xlabel(ax1,'z [m]')
ylabel(ax1,'r [m]')
title(ax1,'Position')
c = colorbar(ax1);
c.Label.String= 'Brillouin Ratio';
view(ax1,2)
caxis([0 1.2])
title(ax1,sprintf('Brillouin Ratio t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens)))
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[16 14];
f.PaperSize=papsize;
print(f,sprintf('%sfluid_Brillouin',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sfluid_Brillouin',name))
%% Omegar
f=figure('Name', sprintf('%s Omegar',name));
ax3=gca;
surface(ax3,M.zgrid,M.rgrid,omegare,'edgecolor','None')
xlabel(ax3,'z [m]')
ylabel(ax3,'r [m]')
xlim(ax3,[M.zgrid(1) M.zgrid(end)])
ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)])
colormap(ax3,'jet')
c = colorbar(ax3);
c.Label.String= '|\omega_\theta| [1/s]';
%c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))];
%caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))])
title(ax3,sprintf('Azimuthal frequency t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens)))
%set(ax3,'colorscale','log')
grid(ax3, 'on')
view(ax3,2)
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[16 14];
f.PaperSize=papsize;
print(f,sprintf('%sfluid_omegar',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sfluid_omegar',name))
%% HP
-if(size(M.H,2)>=2)
+if(M.R.nt>=2)
taveragestart=max(floor(0.8*size(M.tpart,1)),2);
M.displayHP(taveragestart);
end
%% 0D diagnostics
BrillouinRatio=2*M.omepe^2/M.omece^2
NpartsTrapped=mean(M.nbparts(end-10:end))*M.msim/M.me
dblengthiter=fieldstart:length(M.t2d);
-Pr=squeeze(mean(M.Presstens(1,:,:,dblengthiter),4));
-Pz=squeeze(mean(M.Presstens(6,:,:,dblengthiter),4));
-enddens=dens;
-Debye_length=sqrt(abs(min(min(min(Pr(Pr>0)./enddens(Pr>0).^2)),min(min(Pz(Pz>0)./enddens(Pz>0).^2))))*M.eps_0/M.qe^2)
+Pr=squeeze(M.Presstens(1,:,:,dblengthiter));
+Pz=squeeze(M.Presstens(6,:,:,dblengthiter));
+enddens=M.N(:,:,dblengthiter);
+meanPr=mean(Pr,3);
+meanPz=mean(Pz,3);
+Tr=Pr./enddens;
+Tz=Pz./enddens;
+Debye_length=sqrt(abs(min(min(min(mean(Tr(Tr>0)./enddens(Tr>0),3))),min(min(mean(Tz(Tz>0)./enddens(Tz>0),3)))))*M.eps_0/M.qe^2)
%% Debye length
f=figure('Name', sprintf('%s Debye',name));
-surface(M.zgrid,M.rgrid,min(Pr./enddens.^2,Pz./enddens.^2)*M.eps_0/M.qe^2)
+subplot(2,1,1)
+surface(M.zgrid,M.rgrid,sqrt(mean(abs(Tr./enddens*M.eps_0/M.qe^2),3)))
colorbar
f.PaperOrientation='landscape';
+xlabel('z [m]')
+ylabel('r [m]')
+c = colorbar;
+c.Label.String= '\lambda_r [m]';
+
+subplot(2,1,2)
+surface(M.zgrid,M.rgrid,sqrt(mean(abs(Tz./enddens*M.eps_0/M.qe^2),3)))
+colorbar
+f.PaperOrientation='landscape';
+xlabel('z [m]')
+ylabel('r [m]')
+c = colorbar;
+c.Label.String= '\lambda_z [m]';
+
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[16 14];
f.PaperSize=papsize;
print(f,sprintf('%s_dblength',name),'-dpdf','-fillpage')
savefig(f,sprintf('%s_dblength',name))
%% Temperature
f=figure('Name', sprintf('%s Temperature',name));
ax1=subplot(2,1,1);
title('Radial temperature')
-surface(M.zgrid,M.rgrid,Pr./enddens/M.qe,'edgecolor','none')
+surface(M.zgrid,M.rgrid,mean(Tr,3)/M.qe,'edgecolor','none')
colormap('jet')
c = colorbar;
c.Label.String= 'T_er [eV]';
templimits1=caxis;
ax2=subplot(2,1,2);
title('Axial tempreature')
-surface(M.zgrid,M.rgrid,Pz./enddens/M.qe,'edgecolor','none')
+surface(M.zgrid,M.rgrid,mean(Tz,3)/M.qe,'edgecolor','none')
colormap('jet')
c = colorbar;
c.Label.String= 'T_ez [eV]';
maxtemp=max([templimits1,caxis]);
caxis(ax1,[0 maxtemp])
caxis(ax2,[0 maxtemp])
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[16 14];
f.PaperSize=papsize;
print(f,sprintf('%s_temp',name),'-dpdf','-fillpage')
savefig(f,sprintf('%s_temp',name))
-function displaypsi(M,deltat)
-%% Psi function
-
-[~, name, ~] = fileparts(M.file);
-f=figure('Name', sprintf('%s Psi',name));
-zpos=floor(length(M.zgrid)/2);
-tinit=1;
-tend=length(M.t2d);
-if(size(M.H,2)<2)
- H0=M.H0;
- P0=M.P0;
-else
- H0=mean(M.H(:,end));
- P0=mean(M.P(:,end));
-end
-lw=1.5;
-Mirrorratio=(M.Rcurv-1)/(M.Rcurv+1);
-locpot=mean(M.pot(:,zpos,tend-deltat:tend),3);
-psi=1+M.qe*locpot(:)/H0-1/(2*M.me*H0)*(P0./M.rgrid+M.qe*0.5*M.B0.*(M.rgrid-M.width/pi*Mirrorratio*cos(2*pi*M.zgrid(zpos)/M.width)*besseli(1,2*pi*M.rgrid/M.width))).^2;
-locdens=mean(M.N(:,zpos,tend-deltat:tend),3);
-[maxn,In]=max(locdens);%M.N(:,zpos,tinit));
-plot(M.rgrid,M.N(:,zpos,tinit),'bx-','DisplayName',sprintf('t=%1.2f[ns]',M.t2d(tinit)*1e9),'linewidth',lw)
-hold on
-plot(M.rgrid,locdens,'rx-','DisplayName',sprintf('t=[%1.2f-%1.2f] [ns] averaged',M.t2d(tend-deltat)*1e9,M.t2d(tend)*1e9),'linewidth',lw)
-plot(M.rgrid(In-2:end),1./M.rgrid(In-2:end)*maxn*M.rgrid(In),'DisplayName','N=c*1/r','linewidth',lw)
-xlabel('r [m]')
-ylabel('n_e [m^-3]')
-I=find(psi>0);
-if (length(I)>1)
-I=[I(1)-2; I(1)-1; I; I(end)+1; I(end)+2];
-else
- I=M.nnr(1):length(psi);
-end
-rq=linspace(M.rgrid(max(I(1),1)),M.rgrid(I(end)),500);
-psiinterp=interp1(M.rgrid(I),psi(I),rq,'pchip');
-zeroindices=find(diff(psiinterp>=0),2);
-maxpsiinterp=max(psiinterp);
-plot(rq,maxn*psiinterp/abs(maxpsiinterp),'Displayname','normalized \Psi [a.u.]','linewidth',lw)
-ylim([0 inf])
-ylimits=ylim;
-for i=1:length(zeroindices)
- border=plot([rq(zeroindices(i)) rq(zeroindices(i))],[0 M.rgrid(In)/M.rgrid(In-2)*maxn],'k--','linewidth',lw);
- set(get(get(border,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
-end
-legend
-xlim([0 0.02])
-grid on
-title(sprintf('Radial density profile at z=%1.2e[m]',M.zgrid(zpos)))
-f.PaperOrientation='landscape';
-[~, name, ~] = fileparts(M.file);
-f.PaperUnits='centimeters';
-papsize=[15 10];
-f.PaperSize=papsize;
-print(f,sprintf('%sPsi',name),'-dpdf','-fillpage')
-savefig(f,sprintf('%sPsi',name))
-end
+%% Wells
+M.Wells(pot,[fieldstart t2dlength],maxdens)
+M.Wells(M.pot(:,:,1),1, max(max(M.N(:,:,1))))
-function Wells(M,pot,fieldtime,maxdens)
-[~, name, ~] = fileparts(M.file);
-f=figure('Name', sprintf('%s Wells',name));
-ax1=gca;
-rpos=zeros(size(pot));
-[r,~]=meshgrid(M.rgrid,M.zgrid);
-rathet=(M.Athet.*r)';
-poteff=zeros(size(pot));
-for j=2 :size(rpos,2)-1
- for i=1:size(rpos,1)
- if(isempty(find(rathet(i,j)>rathet(:,1),1,'last')))
- rpos(i,j)=1;
- rinterp=0;
- else
- rpos(i,j)=find(rathet(i,j)>rathet(:,1),1,'last');
- rinterp=(rathet(i,j)-rathet(rpos(i,j),1))./(rathet(rpos(i,j)+1,1)-rathet(rpos(i,j),1));
- end
- bmin=M.B(1,rpos(i,j))+rinterp*(M.B(1,rpos(i,j)+1)-M.B(1,rpos(i,j)));
- potmin=pot(rpos(i,j),1)+rinterp*(pot(rpos(i,j)+1,1)-pot(rpos(i,j),1));
- poteff(i,j)=-(pot(i,j)-potmin)*bmin/(M.B(j,i)-bmin);
- end
-end
-h=surface(ax1,M.zgrid,M.rgrid,poteff,'edgecolor','none','facecolor','interp');
-Blines=rathet';
-zindex=floor(length(M.zgrid/2));
-levels=logspace( log10(min(min(Blines(:,2:end)))), log10(max(max(Blines(:,:)))),8);
-hold on
-%contour(ax1,M.zgrid,M.rgrid,Blines',levels,'r-.','linewidth',1.5,'Displayname','Magnetic field lines');
-xlim(ax1,[M.zgrid(1) M.zgrid(end)])
-ylim(ax1,[M.rgrid(1) M.rgrid(end)])
-legend
-xlabel(ax1,'z [m]')
-ylabel(ax1,'r [m]')
-if size(fieldtime)<2
- title(ax1,sprintf('Ion Well t=%1.2g s n_e=%1.2gm^{-3}',M.t2d(fieldtime),double(maxdens)))
- filename=sprintf('%sWell%d',name,fieldtime);
-else
- title(ax1,sprintf('Ion Well t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldtime(1)),M.t2d(fieldtime(2)),double(maxdens)))
- filename=sprintf('%sWell%d_%d',name,fieldtime);
-end
-c = colorbar(ax1);
-c.Label.String= 'Potential [eV]';
-view(ax1,2)
-%set(h,'edgecolor','none');
-grid on;
-f.PaperOrientation='landscape';
-f.PaperUnits='centimeters';
-papsize=[16 14];
-f.PaperSize=papsize;
-print(f,filename,'-dpdf','-fillpage')
-savefig(f,filename)
-end
\ No newline at end of file
diff --git a/matlab/coaxial_ion_loss_criteria.m b/matlab/coaxial_ion_loss_criteria.m
index 0a8e4b6..c186d09 100644
--- a/matlab/coaxial_ion_loss_criteria.m
+++ b/matlab/coaxial_ion_loss_criteria.m
@@ -1,230 +1,295 @@
%% wr+
H0=3.2e-14;
P0=8.66e-25;
r0=0.005;
qe=1.60217662E-19;
me=9.109383E-31;
eps_0=8.85418781762e-12;
vlight=299792458;
ne=5e15;
ne=5*logspace(10,19,5000);
kb=1.38e-23;
T=300;
P=1e-9*100000;
estimneutraldensity=P/(kb*T)
I=40;
V=90;
alpha=1.3;
B0=0.21;
L=0.48;
z=0.0;
% Davidson
prefix='Davidson';
phia=-60000;
phib=0;
R=1.5;
B0=0.21;
width=0.48;
b=0.16;
a=0.001;
-r=0.0075;
-z=0.22;
-deltar=(R-1)/(R+1)*besseli(1,2*pi*r/width)*(width/2/pi)*(cos(2*pi*z/width)+1);
-rb_=0.007;
-rbplus=0.009872;
+r=0.0079;
+zmin=0.0;
+zmax=-width/2;
+deltar=(R-1)/(R+1)*besseli(1,2*pi*r/width)*(width/2/pi)*(cos(2*pi*zmin/width)-cos(2*pi*zmax/width));
+rb_=0.00764;
+rbplus=0.0842;
Te=400; %eV
%1rbplus=0.020;
+% Davidson test
+prefix='Davidson';
+phia=-60000;
+phib=0;
+R=1.5;
+B0=0.21;
+width=0.48;
+b=0.16;
+a=0.001;
+rb_=0.0076;
+r=9.7e-3 ;
+zmin=0.0;
+zmax=-width/8;
+deltar=(R-1)/(R+1)*besseli(1,2*pi*r/width)*(width/2/pi)*(cos(2*pi*zmin/width)-cos(2*pi*zmax/width));
+rbplus=rb_ + max(deltar,0.0026);
+rcenter=r+deltar
+Te=400; %eV
+%1rbplus=0.020;
+
+
% rb_=0.01;
% rbplus=rb_+0.1;
% r=(rb_+rbplus)/2;
% r=0.015
% deltar=(R-1)/(R+1)*r;
% ITER
% prefix='ITER'
% R=1.0001;
% phia=-10000;
% phib=0;
% b=0.09;
% a=0.068;
% r=0.085;
% deltar=(R-1)/(R+1)*r;
% rb_=0.075;
% rbplus=0.088;
spcharge=@(ne)-qe*ne/2/eps_0;
%
% Apart=@(ne)spcharge(ne)/2*(2*r*deltar+deltar^2);
%
% Bpart=@(ne)log(1+deltar/r)/log(b/a)*(phib-phia+spcharge(ne)*rb_^2*log(rbplus/rb_));
%
% Cpart=@(ne)-log(1+deltar/r)/log(rbplus/rb_)*(phib+spcharge(ne)*rb_^2*log(rbplus/rb_));
%
% Dpart=@(ne)spcharge(ne)*(rbplus^2-rb_^2)/(log(b/a)*log(rbplus/rb_))*log(1+deltar/r)*...
% ( log(b/rbplus)*(log(rbplus/a)+1/2) - log(rb_/a)*(log(b/rbplus)-1/2) );
deltaphicloud=@(ne) Phi(r+deltar,rb_,rbplus,a,b,phia,phib, ne)-Phi(r,rb_,rbplus,a,b,phia,phib, ne);
-deltaphivacuum=(phib-phia)/log(b/a)*log(2*R/(R+1));
-
+deltaphivacuum=(phib-phia)/log(b/a)*log((r+deltar)/r);
+
+
+% figure
+% rforphi=linspace(0.001,0.16,1000);
+% phir=zeros(length(rforphi),1);
+% for i=1:length(rforphi)
+% phir(i)=Phi(rforphi(i),rb_,rbplus,a,b,phia,phib, 1e15);
+% end
+% plot(rforphi,phir)
+% figure
+% plot((rforphi(1:end-1)+rforphi(2:end))/2,diff(phir)./diff(rforphi))
+%
Emaxcloud=1/(R-1)*deltaphicloud(ne);
Emaxvacuum=1/(R-1)*deltaphivacuum
+rgrid=linspace(1e-3,16e-2,500);
+zgrid=linspace(-0.24,0.24,1000);
+GainedE=zeros(length(rgrid),length(zgrid));
+for i=1:length(zgrid)
+ drgrid=(R-1)/(R+1)*besseli(1,2*pi*rgrid/width)*(width/2/pi)*(cos(2*pi*0/width)-cos(2*pi*zgrid(i)/width));
+ GainedE(:,i)=-(Phivacuum(rgrid,a,b,phia,phib)-Phivacuum(rgrid+drgrid,a,b,phia,phib));
+end
+figure
+title(sprintf('\\phi_a=%3.2f kV \\phi_b=%3.2f kV R=%3.2f',phia/1e3,phib/1e3, R))
+surface(zgrid,rgrid,GainedE);
+xlabel('z')
+ylabel('r')
+c=colorbar;
+c.Label.String='E [eV]';
+meanE=mean(mean(abs(GainedE(20:375,ceil(length(zgrid)/2):floor(5*length(zgrid)/6))),2))/2
+
figure('name',sprintf('%s Phi_a=%3.2f kV Phi_b=%3.2f kV r=%3.2f mm R=%3.2f',prefix,phia/1e3,phib/1e3,r*1e3, R));
loglog(ne,abs(Emaxcloud),'displayname','flat top density cloud')
hold on
loglog([ne(1) ne(end)], Emaxvacuum*[1 1],'displayname','Vacuum')
xlabel('electron density [m^{-3}]')
ylabel('|E_{min}| [eV]')
title(sprintf('%s \\Phi_a=%3.2f kV \\Phi_b=%3.2f kV r=%3.2f mm R=%3.2f',prefix,phia/1e3,phib/1e3,r*1e3, R))
savefig(sprintf('%s_Phia%3.2f_Phib%3.2f_r%3.2f_R%3.2f.fig',prefix,phia/1e3,phib/1e3,r*1e3, R))
figure
subplot(2,2,1)
sgtitle(sprintf('%s \\Phi_a=%3.2f kV \\Phi_b=%3.2f kV r_0=%3.2f mm z_0=%3.2f mm R=%3.2f',prefix,phia/1e3,phib/1e3,(r+deltar)*1e3,z*1e3, R))
neinter=5e13;
vpar=linspace(1,vlight,1000);
vper=sqrt(-2*qe/me*deltaphicloud(neinter)/(R-1)+vpar.^2/(R-1))/vlight;
vpar=vpar/vlight;
vpar=vpar(real(vper)~=0);
vper=vper(real(vper)~=0);
%Davidson vperp
phi0=Phi(r,rb_,rbplus,a,b,phia,phib, neinter);
vrz=sqrt(2/me*H0+2/me*qe*phi0-1/me^2*(P0/r+qe*Athet(r,z,B0,width,R))^2);
theta=2*pi*linspace(0,1,1000);
vthet=(P0/r+qe*Athet(r,z,B0,width,R))/me;
vperdav=sqrt((vrz^2*cos(theta).^2+vthet^2))/vlight;
vpardav=vrz*sin(theta)/vlight;
vpargauss=sqrt(Te*qe/me)*linspace(-1,1,200);
vpergauss=sqrt(2*Te*qe/me)*sqrt(1-vpargauss.^2*me/Te/qe);
vpargauss=[flip(vpargauss(2:end)) vpargauss]/vlight;
vpergauss=[-flip(vpergauss(2:end)) vpergauss]/vlight;
plot(vpar,vper,'b-')
hold on
plot(vpar,-vper,'b-')
plot(-vpar,vper,'b-')
plot(-vpar,-vper,'b-')
xlabel('\beta_z [m/s]')
ylabel('\beta_\perp [m/s]')
title(sprintf('electrons n_e=%3.2e m^{-3}',neinter))
grid
-plot(vpardav,vperdav,'r--')
-plot(vpargauss,vpergauss,'g--')
+h(1)=plot(vpardav,vperdav,'r--','Displayname',sprintf('Dav H0=%#.2g eV P0=%#.2g kg m^2 s^{-1}',H0/qe,P0));
+h(2)=plot(vpargauss,vpergauss,'g--','Displayname',sprintf('Gauss Te=%#.2g eV',Te));
+legend(h)
+legend('location','northoutside')
xlim([-1 1])
ylim([-1 1])
-axis equal
subplot(2,2,3)
neinter=5e17;
vpar=linspace(1,vlight,1000);
vper=sqrt(-2*qe/me*deltaphicloud(neinter)/(R-1)+vpar.^2/(R-1))/vlight;
vpar=vpar/vlight;
vpar=vpar(real(vper)~=0);
vper=vper(real(vper)~=0);
%Davidson vperp
phi0=Phi(r,rb_,rbplus,a,b,phia,phib, neinter);
vrz=sqrt(2/me*H0+2/me*qe*phi0-1/me^2*(P0/r+qe*Athet(r,z,B0,width,R))^2);
theta=2*pi*linspace(0,1,1000);
vthet=(P0/r+qe*Athet(r,z,B0,width,R))/me;
vperdav=sqrt((vrz^2*cos(theta).^2+vthet^2))/vlight;
vpardav=vrz*sin(theta)/vlight;
vpargauss=sqrt(Te*qe/me)*linspace(-1,1,200);
vpergauss=sqrt(2*Te*qe/me)*sqrt(1-vpargauss.^2*me/Te/qe);
vpargauss=[flip(vpargauss(2:end)) vpargauss]/vlight;
vpergauss=[-flip(vpergauss(2:end)) vpergauss]/vlight;
plot(vpar,vper,'b-')
hold on
plot(vpar,-vper,'b-')
plot(-vpar,vper,'b-')
plot(-vpar,-vper,'b-')
xlabel('\beta_z [m/s]')
ylabel('\beta_\perp [m/s]')
title(sprintf('electrons n_e=%3.2e m^{-3}',neinter))
grid
-plot(vpardav,vperdav,'r--')
-plot(vpargauss,vpergauss,'g--')
+h(1)=plot(vpardav,vperdav,'r--','Displayname',sprintf('Dav H0=%#.2g eV P0=%#.2g kg m^2 s^{-1}',H0/qe,P0));
+h(2)=plot(vpargauss,vpergauss,'g--','Displayname',sprintf('Gauss Te=%#.2g eV',Te));
xlim([-1 1])
ylim([-1 1])
-axis equal
subplot(2,2,2)
neinter=5e13;
vpar=linspace(1,vlight,1000);
m=28.0134/1000*6.02214076e23;
vper=sqrt(2*qe/m*deltaphicloud(neinter)/(R-1)+vpar.^2/(R-1))/vlight;
vpar=vpar/vlight;
vpar=vpar(real(vper)~=0);
vper=vper(real(vper)~=0);
plot(vpar,vper,'b-')
hold on
plot(vpar,-vper,'b-')
plot(-vpar,vper,'b-')
plot(-vpar,-vper,'b-')
xlabel('\beta_z [m/s]')
ylabel('\beta_\perp [m/s]')
title(sprintf('N_2^+ ions n_e=%3.2e m^{-3}',neinter))
grid
xlim([-1 1])
ylim([-1 1])
subplot(2,2,4)
neinter=5e17;
vpar=linspace(1,vlight,1000);
vper=sqrt(2*qe/m*deltaphicloud(neinter)/(R-1)+vpar.^2/(R-1))/vlight;
vpar=vpar/vlight;
vpar=vpar(real(vper)~=0);
vper=vper(real(vper)~=0);
plot(vpar,vper,'b-')
hold on
plot(vpar,-vper,'b-')
plot(-vpar,vper,'b-')
plot(-vpar,-vper,'b-')
xlabel('\beta_z [m/s]')
ylabel('\beta_\perp [m/s]')
title(sprintf('N_2^+ ions n_e=%3.2e m^{-3}',neinter))
grid
savefig(sprintf('%s_Phia%3.2f_Phib%3.2f_r%3.2f_z%3.2f_R%3.2f_cone.fig',prefix,phia/1e3,phib/1e3,(r+deltar)*1e3,z*1e3, R))
xlim([-1 1])
ylim([-1 1])
+%%
+figure
+r=double(M.rgrid);
+r(r<1e-3)=1e-3;
+r(r>16e-2)=16e-2;
+phi=zeros(length(r),1);
+for i=1:length(r)
+ phi(i)=Phi(r(i),14e-2,15e-2,0.001,16e-2,-60000,0,0);
+end
+plot(r,phi)
+
+phim=phi*ones(1,length(M.zgrid));
+
function Atheta=Athet(r,z,B0,width,R)
Atheta=0.5*B0*(r-width/pi*(R-1)/(R+1)...
.*besseli(1,2*pi*r/width).*cos(2*pi*z/width));
end
function phi=Phi(r,rbm,rbp,a,b,phia,phib, ne)
qe=1.60217662E-19;
me=9.109383E-31;
eps_0=8.85418781762e-12;
+q=-qe;
-phibp=1/log(b/a)*(ne/2*qe/eps_0*rbm^2*log(b/rbp)*log(rbp/rbm)+(-qe*ne/2/eps_0*(rbp^2-rbm^2)*(log(rbp/a)+0.5)+phia-phib)*log(b/rbp));
-phibm=1/log(b/a)*(-qe/eps_0*ne/2*rbm^2*log(rbm/a)*log(rbp/rbm)+(-qe/eps_0*ne/2*(rbp^2-rbm^2)*(log(b/rbp)-0.5)+phib)*log(rbm/a)+phia*log(b/rbm));
-if(r<rbm && r>a)
+phibp=1/log(b/a)*(-ne/2*q/eps_0*rbm^2*log(b/rbp)*log(rbp/rbm)+(-q*ne/2/eps_0*(rbp^2-rbm^2)*(log(rbp/a)-0.5)+phia)*log(b/rbp)+phib*log(rbp/a));
+phibm=1/log(b/a)*(q/eps_0*ne/2*rbm^2*log(rbm/a)*log(rbp/rbm)+(-q/eps_0*ne/2*(rbp^2-rbm^2)*(log(b/rbp)+0.5)+phib)*log(rbm/a)+phia*log(b/rbm));
+if(r<rbm && r>=a)
phi=((phibm-phia)*log(r)+phia*log(rbm)-phibm*log(a))/log(rbm/a);
elseif(r>=rbm && r<=rbp)
- phi=-qe/eps_0*ne/4*(r^2-rbp^2)+phibp+(phibp-phibm+qe/eps-0*ne/4*(rbp^2-rbm^2))*log(r/rbp)/log(rbp/rbm);
-elseif(r>rbp && r< b)
+ phi=-qe/eps_0*ne/4*(r^2-rbp^2)+phibp+(phibp-phibm+qe/eps_0*ne/4*(rbp^2-rbm^2))*log(r/rbp)/log(rbp/rbm);
+elseif(r>rbp && r<= b)
phi=((phib-phibp)*log(r)+phibp*log(b)-phib*log(rbp))/log(b/rbp);
else
error('invalid radial position')
end
+end
+function phi=Phivacuum(r,a,b,phia,phib)
+ phi=((phib-phia)*log(r)+phia*log(b)-phib*log(a))/log(b/a);
end
\ No newline at end of file
diff --git a/matlab/dispespicFields.m b/matlab/dispespicFields.m
index 5d481d3..52e5821 100644
--- a/matlab/dispespicFields.m
+++ b/matlab/dispespicFields.m
@@ -1,192 +1,195 @@
function dispespicFields(M,logdensity,showgrid,fixed)
-%UNTITLED Summary of this function goes here
-% Detailed explanation goes here
+%dispespicFields Allows to display the time evolution of the density, electric potential and electric fields
+% M is of class espic2dhdf5 and contains the simulation results
fieldstep=1;
if nargin <2
logdensity=false;
showgrid=false;
fixed=false;
end
if nargin <3
showgrid=false;
fixed=false;
end
if nargin <4
fixed=false;
end
fixed=fi(fixed);
f=uifigure('Name',sprintf('Grid data %s',M.name));
mf=uipanel(f,'Position',[5 50 f.Position(3)-10 f.Position(4)-55]);
mf.AutoResizeChildren='off';
m=uipanel(f,'Position',[5 5 f.Position(3)-10 40]);
sgtitle(mf,sprintf('step=%d t=%0.5e s',fieldstep*M.it1,M.t2d(fieldstep)))
sld = uislider(m,'Position',[10 30 0.6*m.Position(3) 3]);
sld.Value=fieldstep;
sld.Limits=[1 size(M.t2d,1)];
edt = uieditfield(m,'numeric','Limits',[1 size(M.t2d,1)],'Value',1);
edt.Position=[sld.Position(1)+sld.Position(3)+25 5 40 20];
edt.RoundFractionalValues='on';
MaxN=0;
Printbt=uibutton(m,'Position',[edt.Position(1)+edt.Position(3)+10 5 40 20],'Text', 'Save');
Play=uibutton(m,'Position',[Printbt.Position(1)+Printbt.Position(3)+10 5 40 20],'Text', 'Play');
Pause=uibutton(m,'Position',[Play.Position(1)+Play.Position(3)+10 5 40 20],'Text', 'Pause');
%Playbt=uibutton(m,'Position',[Printbt.Position(1)+Printbt.Position(3)+10 5 40 20],'Text', 'Play/Pause');
stop=false;
sld.ValueChangingFcn={@updatefigdata,edt,mf};
edt.ValueChangedFcn={@updatefigdata,sld,mf};
Printbt.ButtonPushedFcn={@plotGridButtonPushed};
Play.ButtonPushedFcn={@plotPlayButtonPushed};
Pause.ButtonPushedFcn={@PauseButtonPushed};
PlotEspic2dgriddata(mf,M,fieldstep);
function plotPlayButtonPushed(btn,ax)
stop=false;
- for i=edt.Value:sld.Limits(2)
+ i=sld.Value;
+ while ~stop
edt.Value=i;
sld.Value=i;
updatesubplotsdata(i,mf);
pause(0.01)
- if stop
- stop=false;
- break;
+ i=sld.Value;
+ i=i+10;
+ if(i>sld.Limits(2))
+ stop=true;
end
-
end
end
function PauseButtonPushed(btn,ax)
stop = true;
end
function PlotEspic2dgriddata(fig,M,fieldstep)
%PlotEspic2dgriddata Plot the 2d data of espic2d at time step fieldstep
sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it1,M.t2d(fieldstep)))
ax1=subplot(2,2,1,'Parent',fig);
sf=surface(ax1,M.zgrid,M.rgrid,M.N(:,:,fieldstep),'edgecolor','none');
if(showgrid)
set(sf,'edgecolor','black')
end
xlim(ax1,[M.zgrid(1) M.zgrid(end)])
ylim(ax1,[M.rgrid(1) M.rgrid(end)])
xlabel(ax1,'z [m]')
ylabel(ax1,'r [m]')
title(ax1,'Position')
c = colorbar(ax1);
c.Label.String= 'n[m^{-3}]';
%c.Limits=[0 max(M.N(:))];
if(isboolean(fixed) && fixed)
climits=caxis(ax1);
MaxN=climits(2);
elseif(~isboolean(fixed))
MaxN=fixed.data;
caxis(ax1,[-Inf MaxN]);
end
view(ax1,2)
if logdensity
set(ax1,'zscale','log')
set(ax1,'colorscale','log')
end
ax2=subplot(2,2,2,'Parent',fig);
contourf(ax2,M.zgrid,M.rgrid,M.pot(:,:,fieldstep));
%plot(ax2,M.zgrid,data.pot(:,5))
xlabel(ax2,'z [m]')
ylabel(ax2,'r [m]')
colormap(ax2,'jet')
c = colorbar(ax2);
c.Label.String= '\Phi [V]';
title(ax2,'Es potential')
grid(ax2, 'on')
ax3=subplot(2,2,3,'Parent',fig);
contourf(ax3,M.zgrid,M.rgrid,M.Er(:,:,fieldstep))
xlabel(ax3,'z [m]')
ylabel(ax3,'r [m]')
c = colorbar(ax3);
c.Label.String= 'E_r [V/m]';
title(ax3,'Radial Electric field')
grid(ax3, 'on')
ax4=subplot(2,2,4,'Parent',fig);
contourf(ax4,M.zgrid,M.rgrid,M.Ez(:,:,fieldstep))
xlabel(ax4,'z [m]')
ylabel(ax4,'r [m]')
c = colorbar(ax4);
c.Label.String= 'E_z [V/m]';
title(ax4,'Axial Electric field')
grid(ax4, 'on')
+
+ linkaxes([ax1,ax2,ax3,ax4]);
end
function plotGridButtonPushed(btn,ax)
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
f=figure();
PlotEspic2dgriddata(f,M,sld.Value);
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
print(f,sprintf('%sGrid%d',name,sld.Value),'-dpdf','-fillpage')
end
function updatefigdata(control, event, Othercontrol, fig)
if strcmp(event.EventName,'ValueChanged')
fieldstep=floor(control.Value);
control.Value=fieldstep;
else
fieldstep=floor(event.Value);
end
Othercontrol.Value=fieldstep;
updatesubplotsdata(fieldstep, fig);
end
function updatesubplotsdata(fieldstep, fig)
- sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it1,double((fieldstep-1)*M.it1)*M.dt))
+ sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it1,M.t2d(floor(fieldstep))))
%% update Position histogram
ax1=fig.Children(end);
dens=M.N(:,:,fieldstep);
ax1.Children.ZData=dens;
ax1.Children.CData=dens;
if(isboolean(fixed) && fixed)
climits=caxis(ax1);
MaxN=climits(2);
nmax=max(dens(:));
if(nmax>MaxN)
MaxN=nmax;
end
caxis(ax1,[0 MaxN]);
elseif(~isboolean(fixed))
MaxN=fixed.data;
caxis(ax1,[-Inf MaxN]);
end
% view(ax1,2)
%% update ES Potential
fig.Children(end-2).Children(1).ZData=M.pot(:,:,fieldstep);
%% update Radial electric field
fig.Children(end-4).Children(1).ZData=M.Er(:,:,fieldstep);
%% update Axial electric field
fig.Children(end-6).Children(1).ZData=M.Ez(:,:,fieldstep);
drawnow limitrate
end
end
diff --git a/matlab/dispespicFluid.m b/matlab/dispespicFluid.m
index 6b48169..64039e7 100644
--- a/matlab/dispespicFluid.m
+++ b/matlab/dispespicFluid.m
@@ -1,148 +1,148 @@
function f=dispespicFluid(M)
fieldstep=1;
f=uifigure('Name','Fluid data');
mf=uipanel(f,'Position',[5 50 f.Position(3)-10 f.Position(4)-55]);
mf.AutoResizeChildren='off';
m=uipanel(f,'Position',[5 5 f.Position(3)-10 40]);
sgtitle(mf,sprintf('t=%0.5e s',M.t2d(fieldstep)))
sld = uislider(m,'Position',[10 30 0.6*m.Position(3) 3]);
sld.Value=fieldstep;
sld.Limits=[1 length(M.t2d)];
edt = uieditfield(m,'numeric','Limits',[1 length(M.t2d)],'Value',1);
edt.Position=[sld.Position(1)+sld.Position(3)+25 5 40 20];
edt.RoundFractionalValues='on';
Printbt=uibutton(m,'Position',[edt.Position(1)+edt.Position(3)+10 5 40 20],'Text', 'Save');
%Playbt=uibutton(m,'Position',[Printbt.Position(1)+Printbt.Position(3)+10 5 40 20],'Text', 'Play/Pause');
sld.ValueChangingFcn={@updatefigdata,edt,mf};
edt.ValueChangedFcn={@updatefigdata,sld,mf};
Printbt.ButtonPushedFcn={@plotGridButtonPushed};
[R,Z]=meshgrid(M.rgrid,M.zgrid);
Rinv=1./R;
Rinv(:,1)=0;
PlotEspic2dfluiddata(mf,M,fieldstep);
function PlotEspic2dfluiddata(fig,M,fieldstep)
%PlotEspic2dgriddata Plot the 2d data of espic2d at time step fieldstep
sgtitle(fig,sprintf('t=%0.5e s',M.t2d(fieldstep)))
ax1=subplot(2,2,1,'Parent',fig);
surface(ax1,M.zgrid,M.rgrid,M.N(:,:,fieldstep),'edgecolor','none');
xlim(ax1,[M.zgrid(1) M.zgrid(end)])
ylim(ax1,[M.rgrid(1) M.rgrid(end)])
xlabel(ax1,'z [m]')
ylabel(ax1,'r [m]')
title(ax1,'Position')
c = colorbar(ax1);
c.Label.String= 'n[m^{-3}]';
%c.Limits=[0 max(M.N(:))];
%caxis(ax1,[0 max(M.N(:))]);
view(ax1,2)
%set(ax1,'colorscale','log')
ax2=subplot(2,2,2,'Parent',fig);
-surface(ax2,M.zgrid,M.rgrid,M.fluidUR(:,:,fieldstep)./M.vlight,'edgecolor','none');
+surface(ax2,M.zgrid,M.rgrid,M.fluidUR(:,:,fieldstep),'edgecolor','none');
%plot(ax2,M.zgrid,data.pot(:,5))
xlabel(ax2,'z [m]')
ylabel(ax2,'r [m]')
colormap(ax2,'jet')
c = colorbar(ax2);
%c.Limits=[min(M.fluidUR(:,:,:)) max(M.fluidUR(:,:,:))];
c.Label.String= 'U_r [m/s]';
title(ax2,'Radial velocity')
grid(ax2, 'on')
%caxis(ax2,[min(M.fluidUR(:)) max(M.fluidUR(:))])
view(ax2,2)
ax3=subplot(2,2,3,'Parent',fig);
surface(ax3,M.zgrid,M.rgrid,abs(M.fluidUTHET(:,:,fieldstep).*Rinv'),'edgecolor','none')
%mean(M.fluidUTHET(:,:,end).*Rinv')
xlabel(ax3,'z [m]')
ylabel(ax3,'r [m]')
colormap(ax3,'jet')
c = colorbar(ax3);
c.Label.String= '\omega_\theta [1/s]';
%c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))];
%caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))])
title(ax3,'Azimuthal velocity')
set(ax3,'colorscale','log')
grid(ax3, 'on')
view(ax3,2)
ax4=subplot(2,2,4,'Parent',fig);
-surface(ax4,M.zgrid,M.rgrid,M.fluidUZ(:,:,fieldstep)./M.vlight,'edgecolor','none')
+surface(ax4,M.zgrid,M.rgrid,M.fluidUZ(:,:,fieldstep),'edgecolor','none')
xlabel(ax4,'z [m]')
ylabel(ax4,'r [m]')
colormap(ax4,'jet')
c = colorbar(ax4);
%c.Limits=[min(M.fluidUZ(:)) max(M.fluidUZ(:))];
c.Label.String= 'U_z [m/s]';
title(ax4,'Axial velocity')
grid(ax4, 'on')
%caxis(ax4,[min(M.fluidUZ(:)) max(M.fluidUZ(:))])
view(ax4,2)
linkaxes([ax1 ax2 ax3 ax4],'xy')
end
function plotGridButtonPushed(btn,ax)
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
f=figure();
PlotEspic2dgriddata(f,M,sld.Value);
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
print(f,sprintf('%sfluid%d',name,sld.Value),'-dpdf','-fillpage')
end
function updatefigdata(control, event, Othercontrol, fig)
if strcmp(event.EventName,'ValueChanged')
fieldstep=floor(control.Value);
control.Value=fieldstep;
else
fieldstep=floor(event.Value);
end
Othercontrol.Value=fieldstep;
sgtitle(fig,sprintf('t=%0.5e s',double((fieldstep-1)*M.it1)*M.dt))
%% update Position histogram
ax1=fig.Children(end);
ax1.Children(1).CData=M.N(:,:,fieldstep);
ax1.Children(1).ZData=M.N(:,:,fieldstep);
view(ax1,2)
%% update Radial velocity
- fig.Children(end-2).Children(1).CData=M.fluidUR(:,:,fieldstep)./M.vlight;
- fig.Children(end-2).Children(1).ZData=M.fluidUR(:,:,fieldstep)./M.vlight;
+ fig.Children(end-2).Children(1).CData=M.fluidUR(:,:,fieldstep);
+ fig.Children(end-2).Children(1).ZData=M.fluidUR(:,:,fieldstep);
%% update Azimuthal velocity
fig.Children(end-4).Children(1).CData=abs(M.fluidUTHET(:,:,fieldstep).*Rinv');
fig.Children(end-4).Children(1).ZData=abs(M.fluidUTHET(:,:,fieldstep).*Rinv');
%% update Axial velocity
- fig.Children(end-6).Children(1).CData=M.fluidUZ(:,:,fieldstep)./M.vlight;
- fig.Children(end-6).Children(1).ZData=M.fluidUZ(:,:,fieldstep)./M.vlight;
+ fig.Children(end-6).Children(1).CData=M.fluidUZ(:,:,fieldstep);
+ fig.Children(end-6).Children(1).ZData=M.fluidUZ(:,:,fieldstep);
drawnow limitrate
end
end
diff --git a/matlab/dispespicParts.m b/matlab/dispespicParts.m
index 245f13f..ead72ca 100644
--- a/matlab/dispespicParts.m
+++ b/matlab/dispespicParts.m
@@ -1,138 +1,138 @@
function dispespicParts(M)
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
fieldstep=1;
f2=uifigure('Name','Particles data');
mf2=uipanel(f2,'Position',[5 50 f2.Position(3)-10 f2.Position(4)-55]);
mf2.AutoResizeChildren='off';
m2=uipanel(f2,'Position',[5 5 f2.Position(3)-10 40]);
sgtitle(mf2,sprintf('step=%d t=%0.5e s',fieldstep*M.it2,M.tpart(fieldstep)))
sld2 = uislider(m2,'Position',[10 30 0.6*m2.Position(3) 3]);
sld2.Value=fieldstep;
-sld2.Limits=[1 size(M.VR,2)];
+sld2.Limits=[1 length(M.tpart)];
-edt2 = uieditfield(m2,'numeric','Limits',[1 size(M.VR,2)],'Value',1);
+edt2 = uieditfield(m2,'numeric','Limits',[1 length(M.tpart)],'Value',1);
edt2.Position=[sld2.Position(1)+sld2.Position(3)+25 5 40 20];
edt2.RoundFractionalValues='on';
Printbt2=uibutton(m2,'Position',[edt2.Position(1)+edt2.Position(3)+10 5 40 20],'Text', 'Save');
sld2.ValueChangingFcn={@updatefigpartdata,edt2,mf2};
edt2.ValueChangedFcn={@updatefigpartdata,sld2,mf2};
Printbt2.ButtonPushedFcn={@plotPartButtonPushed};
PlotEspic2dpartdata(mf2,M,fieldstep);
function PlotEspic2dpartdata(fig,M,fieldstep)
%PlotEspic2dgriddata Plot the 2d data of espic2d at time step fieldstep
sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it2,M.tpart(fieldstep)))
ax1=subplot(3,2,1,'Parent',fig);
plot(ax1,M.Z(:,fieldstep),M.R(:,fieldstep),'.');
xlim(ax1,[M.zgrid(1) M.zgrid(end)])
ylim(ax1,[M.rgrid(1) M.rgrid(end)])
xlabel(ax1,'Z [m]')
ylabel(ax1,'R [m]')
title(ax1,'Position')
grid(ax1, 'on')
ax2=subplot(3,2,2,'Parent',fig);
plot(ax2,M.VZ(:,fieldstep),M.VR(:,fieldstep),'.');
xlabel(ax2,'VZ [m/s]')
ylabel(ax2,'VR [m/s]')
axis(ax2, 'equal')
grid(ax2, 'on')
ax3=subplot(3,2,3,'Parent',fig);
plot(ax3,M.Z(:,fieldstep),M.VZ(:,fieldstep),'.');
xlim(ax3,[M.zgrid(1) M.zgrid(end)])
xlabel(ax3,'Z [m]')
ylabel(ax3,'VZ [m/s]')
grid(ax3, 'on')
ax4=subplot(3,2,4,'Parent',fig);
plot(ax4,M.VZ(:,fieldstep),M.R(:,fieldstep),'.')
ylabel(ax4,'R [m]')
xlabel(ax4,'VZ [m/s]')
ylim(ax4,[M.rgrid(1) M.rgrid(end)])
grid(ax4, 'on')
ax5=subplot(3,2,5,'Parent',fig);
plot(ax5,M.Z(:,fieldstep),M.VR(:,fieldstep),'.');
xlim(ax5,[M.zgrid(1) M.zgrid(end)])
xlabel(ax5,'Z [m]')
ylabel(ax5,'VR [m/s]')
grid(ax5, 'on')
ax6=subplot(3,2,6,'Parent',fig);
plot(ax6,M.VR(:,fieldstep),M.R(:,fieldstep),'.');
ylim(ax6,[M.rgrid(1) M.rgrid(end)])
ylabel(ax6,'R [m]')
xlabel(ax6,'VR [m/s]')
grid(ax6, 'on')
end
function updatefigpartdata(control, event, Othercontrol, fig)
if strcmp(event.EventName,'ValueChanged')
fieldstep=floor(control.Value);
control.Value=fieldstep;
else
fieldstep=floor(event.Value);
end
Othercontrol.Value=fieldstep;
- sgtitle(fig,sprintf('t=%0.5e s',M.tpart(fieldstep-1)))
+ sgtitle(fig,sprintf('t=%0.5e s',M.tpart(fieldstep)))
%% update Position plot
ax1=fig.Children(end);
ax1.Children.XData=M.Z(:,fieldstep);
ax1.Children.YData=M.R(:,fieldstep);
view(ax1,2)
%% update VZ VPERP
fig.Children(end-1).Children(1).XData=M.VZ(:,fieldstep);
fig.Children(end-1).Children(1).YData=M.VR(:,fieldstep);
axis(fig.Children(end-1),'equal')
%% update Z VZ
fig.Children(end-2).Children(1).XData=M.Z(:,fieldstep);
fig.Children(end-2).Children(1).YData=M.VZ(:,fieldstep);
%% update VZ VR
fig.Children(end-3).Children(1).XData=M.VZ(:,fieldstep);
fig.Children(end-3).Children(1).YData=M.R(:,fieldstep);
%% update R VR
fig.Children(end-4).Children(1).XData=M.Z(:,fieldstep);
fig.Children(end-4).Children(1).YData=M.VR(:,fieldstep);
%% update Z VR
fig.Children(end-5).Children(1).XData=M.VR(:,fieldstep);
fig.Children(end-5).Children(1).YData=M.R(:,fieldstep);
drawnow limitrate
end
function plotPartButtonPushed(btn,ax)
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
f=figure();
PlotEspic2dpartdata(f,M,sld.Value);
f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
print(f,sprintf('%sParts%d',name,sld.Value),'-dpdf','-fillpage')
end
end
diff --git a/matlab/dispespicFields.m b/matlab/dispespicPressure.m
similarity index 51%
copy from matlab/dispespicFields.m
copy to matlab/dispespicPressure.m
index 5d481d3..a1644ab 100644
--- a/matlab/dispespicFields.m
+++ b/matlab/dispespicPressure.m
@@ -1,192 +1,185 @@
-function dispespicFields(M,logdensity,showgrid,fixed)
-%UNTITLED Summary of this function goes here
-% Detailed explanation goes here
+function dispespicPressure(M,logdensity,showgrid,fixed,temperature)
+%dispespicPressure Allows to display the time evolution of the pressure tensor
+% M is of class espic2dhdf5 and contains the simulation results
fieldstep=1;
if nargin <2
logdensity=false;
- showgrid=false;
- fixed=false;
end
if nargin <3
showgrid=false;
- fixed=false;
end
if nargin <4
fixed=false;
end
+if nargin <5
+ temperature=false;
+end
fixed=fi(fixed);
-f=uifigure('Name',sprintf('Grid data %s',M.name));
+f=uifigure('Name',sprintf('Pressure data %s',M.name));
mf=uipanel(f,'Position',[5 50 f.Position(3)-10 f.Position(4)-55]);
mf.AutoResizeChildren='off';
m=uipanel(f,'Position',[5 5 f.Position(3)-10 40]);
sgtitle(mf,sprintf('step=%d t=%0.5e s',fieldstep*M.it1,M.t2d(fieldstep)))
sld = uislider(m,'Position',[10 30 0.6*m.Position(3) 3]);
sld.Value=fieldstep;
sld.Limits=[1 size(M.t2d,1)];
edt = uieditfield(m,'numeric','Limits',[1 size(M.t2d,1)],'Value',1);
edt.Position=[sld.Position(1)+sld.Position(3)+25 5 40 20];
edt.RoundFractionalValues='on';
-MaxN=0;
+MaxP=0;
+MinP=0;
+Paxes=gobjects(6);
+
Printbt=uibutton(m,'Position',[edt.Position(1)+edt.Position(3)+10 5 40 20],'Text', 'Save');
Play=uibutton(m,'Position',[Printbt.Position(1)+Printbt.Position(3)+10 5 40 20],'Text', 'Play');
Pause=uibutton(m,'Position',[Play.Position(1)+Play.Position(3)+10 5 40 20],'Text', 'Pause');
%Playbt=uibutton(m,'Position',[Printbt.Position(1)+Printbt.Position(3)+10 5 40 20],'Text', 'Play/Pause');
stop=false;
sld.ValueChangingFcn={@updatefigdata,edt,mf};
edt.ValueChangedFcn={@updatefigdata,sld,mf};
Printbt.ButtonPushedFcn={@plotGridButtonPushed};
Play.ButtonPushedFcn={@plotPlayButtonPushed};
Pause.ButtonPushedFcn={@PauseButtonPushed};
PlotEspic2dgriddata(mf,M,fieldstep);
function plotPlayButtonPushed(btn,ax)
stop=false;
for i=edt.Value:sld.Limits(2)
edt.Value=i;
sld.Value=i;
updatesubplotsdata(i,mf);
pause(0.01)
if stop
stop=false;
break;
end
end
end
function PauseButtonPushed(btn,ax)
stop = true;
end
+
function PlotEspic2dgriddata(fig,M,fieldstep)
%PlotEspic2dgriddata Plot the 2d data of espic2d at time step fieldstep
sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it1,M.t2d(fieldstep)))
-
-
- ax1=subplot(2,2,1,'Parent',fig);
- sf=surface(ax1,M.zgrid,M.rgrid,M.N(:,:,fieldstep),'edgecolor','none');
- if(showgrid)
- set(sf,'edgecolor','black')
+ if temperature
+ Clabels={'T_{rr} [eV]', 'T_{r\theta} [eV]', 'T_{rz} [eV]', 'T_{\theta\theta} [eV]', 'T_{\theta z} [eV]', 'T_{zz} [eV]'};
+ else
+ Clabels={'P_{rr} [Pa]', 'P_{r\theta} [Pa]', 'P_{rz} [Pa]', 'P_{\theta\theta} [Pa]', 'P_{\theta z} [Pa]', 'P_{zz} [Pa]'};
end
- xlim(ax1,[M.zgrid(1) M.zgrid(end)])
- ylim(ax1,[M.rgrid(1) M.rgrid(end)])
- xlabel(ax1,'z [m]')
- ylabel(ax1,'r [m]')
- title(ax1,'Position')
- c = colorbar(ax1);
- c.Label.String= 'n[m^{-3}]';
- %c.Limits=[0 max(M.N(:))];
- if(isboolean(fixed) && fixed)
- climits=caxis(ax1);
- MaxN=climits(2);
- elseif(~isboolean(fixed))
- MaxN=fixed.data;
- caxis(ax1,[-Inf MaxN]);
+ for i=1:6
+ Paxes(i)=subplot(2,3,i,'Parent',fig);
+ if temperature
+ N=M.N(:,:,fieldstep);
+ invN=1./N;
+ invN(isinf(invN))=0;
+ p=squeeze(M.Presstens(i,:,:,fieldstep)).*invN/M.qe;
+ else
+ p=squeeze(M.Presstens(i,:,:,fieldstep));
+ end
+ sf=surface(Paxes(i),M.zgrid,M.rgrid,p,'edgecolor','none');
+ if(showgrid)
+ set(sf,'edgecolor','black')
+ end
+ xlim(Paxes(i),[M.zgrid(1) M.zgrid(end)])
+ ylim(Paxes(i),[M.rgrid(1) M.rgrid(end)])
+ xlabel(Paxes(i),'z [m]')
+ ylabel(Paxes(i),'r [m]')
+ title(Paxes(i),'Position')
+ c = colorbar(Paxes(i));
+ c.Label.String= Clabels{i};
+ if(isboolean(fixed) && fixed)
+ climits=caxis(Paxes(i));
+ MaxP=max(MaxP,climits(2));
+ MinP=min(MinP,climits(1));
+ elseif(~isboolean(fixed))
+ MaxP=fixed.data;
+ MinP=-Inf;
+ caxis(ax1,[MinP MaxP]);
+ end
+ view(Paxes(i),2)
+
+ if logdensity
+ set(Paxes(i),'zscale','log')
+ set(Paxes(i),'colorscale','log')
+ end
end
- view(ax1,2)
-
- if logdensity
- set(ax1,'zscale','log')
- set(ax1,'colorscale','log')
+ if(isboolean(fixed) && fixed)
+ for i=1:6
+ caxis(Paxes(i),[MinP MaxP]);
+ end
end
-
- ax2=subplot(2,2,2,'Parent',fig);
- contourf(ax2,M.zgrid,M.rgrid,M.pot(:,:,fieldstep));
- %plot(ax2,M.zgrid,data.pot(:,5))
- xlabel(ax2,'z [m]')
- ylabel(ax2,'r [m]')
- colormap(ax2,'jet')
- c = colorbar(ax2);
- c.Label.String= '\Phi [V]';
- title(ax2,'Es potential')
- grid(ax2, 'on')
-
-
- ax3=subplot(2,2,3,'Parent',fig);
- contourf(ax3,M.zgrid,M.rgrid,M.Er(:,:,fieldstep))
- xlabel(ax3,'z [m]')
- ylabel(ax3,'r [m]')
- c = colorbar(ax3);
- c.Label.String= 'E_r [V/m]';
- title(ax3,'Radial Electric field')
- grid(ax3, 'on')
-
- ax4=subplot(2,2,4,'Parent',fig);
- contourf(ax4,M.zgrid,M.rgrid,M.Ez(:,:,fieldstep))
- xlabel(ax4,'z [m]')
- ylabel(ax4,'r [m]')
- c = colorbar(ax4);
- c.Label.String= 'E_z [V/m]';
- title(ax4,'Axial Electric field')
- grid(ax4, 'on')
+ linkaxes(Paxes);
end
function plotGridButtonPushed(btn,ax)
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
f=figure();
PlotEspic2dgriddata(f,M,sld.Value);
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
print(f,sprintf('%sGrid%d',name,sld.Value),'-dpdf','-fillpage')
end
function updatefigdata(control, event, Othercontrol, fig)
if strcmp(event.EventName,'ValueChanged')
fieldstep=floor(control.Value);
control.Value=fieldstep;
else
fieldstep=floor(event.Value);
end
Othercontrol.Value=fieldstep;
updatesubplotsdata(fieldstep, fig);
end
function updatesubplotsdata(fieldstep, fig)
sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it1,double((fieldstep-1)*M.it1)*M.dt))
- %% update Position histogram
- ax1=fig.Children(end);
- dens=M.N(:,:,fieldstep);
-
- ax1.Children.ZData=dens;
- ax1.Children.CData=dens;
+ %% update Pressure subplots
+ for i=1:6
+ if temperature
+ N=M.N(:,:,fieldstep);
+ invN=1./N;
+ invN(isinf(invN))=0;
+ P=squeeze(M.Presstens(i,:,:,fieldstep)).*invN/M.qe;
+ else
+ P=squeeze(M.Presstens(i,:,:,fieldstep));
+ end
+ Paxes(i).Children.ZData=P;
+ Paxes(i).Children.CData=P;
+ if(isboolean(fixed) && fixed)
+ climits=caxis(Paxes(i));
+ MaxP=max([P(:);climits(2)]);
+ MinP=min([P(:);climits(1)]);
+ caxis(Paxes(i),[0 MaxP]);
+ elseif(~isboolean(fixed))
+ MaxP=fixed.data;
+ caxis(Paxes(i),[-Inf MaxP]);
+ end
+ end
if(isboolean(fixed) && fixed)
- climits=caxis(ax1);
- MaxN=climits(2);
- nmax=max(dens(:));
- if(nmax>MaxN)
- MaxN=nmax;
+ for i=1:6
+ caxis(Paxes(i),[MinP MaxP]);
end
- caxis(ax1,[0 MaxN]);
- elseif(~isboolean(fixed))
- MaxN=fixed.data;
- caxis(ax1,[-Inf MaxN]);
end
-
- % view(ax1,2)
- %% update ES Potential
- fig.Children(end-2).Children(1).ZData=M.pot(:,:,fieldstep);
- %% update Radial electric field
- fig.Children(end-4).Children(1).ZData=M.Er(:,:,fieldstep);
- %% update Axial electric field
- fig.Children(end-6).Children(1).ZData=M.Ez(:,:,fieldstep);
drawnow limitrate
-
-
end
end
diff --git a/matlab/dispespicFields.m b/matlab/dispespicWell.m
similarity index 58%
copy from matlab/dispespicFields.m
copy to matlab/dispespicWell.m
index 5d481d3..b165db1 100644
--- a/matlab/dispespicFields.m
+++ b/matlab/dispespicWell.m
@@ -1,192 +1,196 @@
-function dispespicFields(M,logdensity,showgrid,fixed)
-%UNTITLED Summary of this function goes here
-% Detailed explanation goes here
+function dispespicWell(M,logdensity,showgrid,fixed)
+%dispespicFields Allows to display the time evolution of the density, electric potential and electric fields
+% M is of class espic2dhdf5 and contains the simulation results
fieldstep=1;
if nargin <2
logdensity=false;
showgrid=false;
fixed=false;
end
if nargin <3
showgrid=false;
fixed=false;
end
if nargin <4
fixed=false;
end
fixed=fi(fixed);
-f=uifigure('Name',sprintf('Grid data %s',M.name));
+f=uifigure('Name',sprintf('Well data %s',M.name));
mf=uipanel(f,'Position',[5 50 f.Position(3)-10 f.Position(4)-55]);
mf.AutoResizeChildren='off';
m=uipanel(f,'Position',[5 5 f.Position(3)-10 40]);
sgtitle(mf,sprintf('step=%d t=%0.5e s',fieldstep*M.it1,M.t2d(fieldstep)))
sld = uislider(m,'Position',[10 30 0.6*m.Position(3) 3]);
sld.Value=fieldstep;
sld.Limits=[1 size(M.t2d,1)];
edt = uieditfield(m,'numeric','Limits',[1 size(M.t2d,1)],'Value',1);
edt.Position=[sld.Position(1)+sld.Position(3)+25 5 40 20];
edt.RoundFractionalValues='on';
MaxN=0;
+Minwell=0;
+Maxwell=0;
+plotaxes=gobjects(2);
+
Printbt=uibutton(m,'Position',[edt.Position(1)+edt.Position(3)+10 5 40 20],'Text', 'Save');
Play=uibutton(m,'Position',[Printbt.Position(1)+Printbt.Position(3)+10 5 40 20],'Text', 'Play');
Pause=uibutton(m,'Position',[Play.Position(1)+Play.Position(3)+10 5 40 20],'Text', 'Pause');
%Playbt=uibutton(m,'Position',[Printbt.Position(1)+Printbt.Position(3)+10 5 40 20],'Text', 'Play/Pause');
stop=false;
sld.ValueChangingFcn={@updatefigdata,edt,mf};
edt.ValueChangedFcn={@updatefigdata,sld,mf};
Printbt.ButtonPushedFcn={@plotGridButtonPushed};
Play.ButtonPushedFcn={@plotPlayButtonPushed};
Pause.ButtonPushedFcn={@PauseButtonPushed};
PlotEspic2dgriddata(mf,M,fieldstep);
function plotPlayButtonPushed(btn,ax)
stop=false;
for i=edt.Value:sld.Limits(2)
edt.Value=i;
sld.Value=i;
updatesubplotsdata(i,mf);
pause(0.01)
if stop
stop=false;
break;
end
end
end
function PauseButtonPushed(btn,ax)
stop = true;
end
function PlotEspic2dgriddata(fig,M,fieldstep)
%PlotEspic2dgriddata Plot the 2d data of espic2d at time step fieldstep
sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it1,M.t2d(fieldstep)))
- ax1=subplot(2,2,1,'Parent',fig);
- sf=surface(ax1,M.zgrid,M.rgrid,M.N(:,:,fieldstep),'edgecolor','none');
+ plotaxes(1)=subplot(1,2,1,'Parent',fig);
+ sf=surface(plotaxes(1),M.zgrid,M.rgrid,M.N(:,:,fieldstep),'edgecolor','none');
if(showgrid)
set(sf,'edgecolor','black')
end
- xlim(ax1,[M.zgrid(1) M.zgrid(end)])
- ylim(ax1,[M.rgrid(1) M.rgrid(end)])
- xlabel(ax1,'z [m]')
- ylabel(ax1,'r [m]')
- title(ax1,'Position')
- c = colorbar(ax1);
+ xlim(plotaxes(1),[M.zgrid(1) M.zgrid(end)])
+ ylim(plotaxes(1),[M.rgrid(1) M.rgrid(end)])
+ xlabel(plotaxes(1),'z [m]')
+ ylabel(plotaxes(1),'r [m]')
+ title(plotaxes(1),'Density')
+ c = colorbar(plotaxes(1));
c.Label.String= 'n[m^{-3}]';
%c.Limits=[0 max(M.N(:))];
if(isboolean(fixed) && fixed)
- climits=caxis(ax1);
+ climits=caxis(plotaxes(1));
MaxN=climits(2);
elseif(~isboolean(fixed))
MaxN=fixed.data;
- caxis(ax1,[-Inf MaxN]);
+ caxis(plotaxes(1),[-Inf MaxN]);
end
- view(ax1,2)
+ view(plotaxes(1),2)
if logdensity
- set(ax1,'zscale','log')
- set(ax1,'colorscale','log')
+ set(plotaxes(1),'zscale','log')
+ set(plotaxes(1),'colorscale','log')
+ end
+
+
+ plotaxes(2)=subplot(1,2,2,'Parent',fig);
+ well=-PotentialWell(M.pot(:,:,fieldstep),M.Athet(:,:),M.rgrid,M.zgrid,false);
+ sf=surface(plotaxes(2),M.zgrid(2:end),M.rgrid(2:end),well(2:end,2:end),'edgecolor','none');
+ if(showgrid)
+ set(sf,'edgecolor','black')
end
+ xlim(plotaxes(2),[M.zgrid(1) M.zgrid(end)])
+ ylim(plotaxes(2),[M.rgrid(1) M.rgrid(end)])
+ xlabel(plotaxes(2),'z [m]')
+ ylabel(plotaxes(2),'r [m]')
+ title(plotaxes(2),'Well')
+ c = colorbar(plotaxes(2));
+ c.Label.String= 'depth [V]';
+ %c.Limits=[0 max(M.N(:))];
+ if(isboolean(fixed) && fixed)
+ climits=caxis(plotaxes(2));
+ Minwell=climits(1);
+ Maxwell=climits(2);
+ elseif(~isboolean(fixed))
+ climits=caxis(plotaxes(2));
+ Minwell=climits(1);
+ Maxwell=climits(2);
+ end
+ view(plotaxes(2),2)
+ colormap(plotaxes(2),'jet');
- ax2=subplot(2,2,2,'Parent',fig);
- contourf(ax2,M.zgrid,M.rgrid,M.pot(:,:,fieldstep));
- %plot(ax2,M.zgrid,data.pot(:,5))
- xlabel(ax2,'z [m]')
- ylabel(ax2,'r [m]')
- colormap(ax2,'jet')
- c = colorbar(ax2);
- c.Label.String= '\Phi [V]';
- title(ax2,'Es potential')
- grid(ax2, 'on')
-
-
- ax3=subplot(2,2,3,'Parent',fig);
- contourf(ax3,M.zgrid,M.rgrid,M.Er(:,:,fieldstep))
- xlabel(ax3,'z [m]')
- ylabel(ax3,'r [m]')
- c = colorbar(ax3);
- c.Label.String= 'E_r [V/m]';
- title(ax3,'Radial Electric field')
- grid(ax3, 'on')
-
- ax4=subplot(2,2,4,'Parent',fig);
- contourf(ax4,M.zgrid,M.rgrid,M.Ez(:,:,fieldstep))
- xlabel(ax4,'z [m]')
- ylabel(ax4,'r [m]')
- c = colorbar(ax4);
- c.Label.String= 'E_z [V/m]';
- title(ax4,'Axial Electric field')
- grid(ax4, 'on')
+ linkaxes(plotaxes);
end
function plotGridButtonPushed(btn,ax)
- %UNTITLED2 Summary of this function goes here
- % Detailed explanation goes here
f=figure();
PlotEspic2dgriddata(f,M,sld.Value);
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
print(f,sprintf('%sGrid%d',name,sld.Value),'-dpdf','-fillpage')
end
function updatefigdata(control, event, Othercontrol, fig)
if strcmp(event.EventName,'ValueChanged')
fieldstep=floor(control.Value);
control.Value=fieldstep;
else
fieldstep=floor(event.Value);
end
Othercontrol.Value=fieldstep;
updatesubplotsdata(fieldstep, fig);
end
function updatesubplotsdata(fieldstep, fig)
sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it1,double((fieldstep-1)*M.it1)*M.dt))
%% update Position histogram
- ax1=fig.Children(end);
dens=M.N(:,:,fieldstep);
- ax1.Children.ZData=dens;
- ax1.Children.CData=dens;
+ well=-PotentialWell(M.pot(:,:,fieldstep),M.Athet(:,:),M.rgrid,M.zgrid,false);
+ well=well(2:end,2:end);
+
+ plotaxes(1).Children.ZData=dens;
+ plotaxes(1).Children.CData=dens;
if(isboolean(fixed) && fixed)
- climits=caxis(ax1);
+ climits=caxis(plotaxes(1));
MaxN=climits(2);
nmax=max(dens(:));
if(nmax>MaxN)
MaxN=nmax;
end
- caxis(ax1,[0 MaxN]);
+ caxis(plotaxes(1),[0 MaxN]);
elseif(~isboolean(fixed))
MaxN=fixed.data;
- caxis(ax1,[-Inf MaxN]);
+ caxis(plotaxes(1),[-Inf MaxN]);
end
- % view(ax1,2)
- %% update ES Potential
- fig.Children(end-2).Children(1).ZData=M.pot(:,:,fieldstep);
- %% update Radial electric field
- fig.Children(end-4).Children(1).ZData=M.Er(:,:,fieldstep);
- %% update Axial electric field
- fig.Children(end-6).Children(1).ZData=M.Ez(:,:,fieldstep);
+ plotaxes(2).Children.ZData=well;
+ plotaxes(2).Children.CData=well;
+ if(isboolean(fixed) && fixed || ~isboolean(fixed))
+ climits=caxis(plotaxes(2));
+ Maxwell=max([well(:);climits(2)]);
+ Minwell=min([well(:);climits(1)]);
+ caxis(plotaxes(2),[Minwell Maxwell]);
+ end
drawnow limitrate
end
end
diff --git a/matlab/distribution_function.m b/matlab/distribution_function.m
index 9cdc897..aacc2f4 100644
--- a/matlab/distribution_function.m
+++ b/matlab/distribution_function.m
@@ -1,83 +1,95 @@
fig=figure;
ax1=gca;
Ndistrib=M.N(:,:,end);
h=pcolor(ax1,M.zgrid,M.rgrid,Ndistrib);
set(h, 'EdgeColor', 'none');
hold on
[r,z]=find(Ndistrib~=0);
xlim(ax1,[M.zgrid(min(z)) M.zgrid(max(z))])
ylim(ax1,[M.rgrid(min(r)) M.rgrid(max(r))])
xlabel(ax1,'Z [m]')
ylabel(ax1,'R [m]')
c = colorbar(ax1);
c.Label.String= 'n [m^{-3}]';
view(ax1,2)
%set(ax1,'colorscale','log')
stepmin=2000;
stepend=size(M.t2d,1);
papsize=[16 14];
[x,y]=ginput(1);
Zindex=find(x>M.zgrid,1,'last');
Rindex=find(y>M.rgrid,1,'last');
Z=(M.zgrid(Zindex));
R=(M.rgrid(Rindex));
plot(Z,R,'rx','Markersize',12);
% Rindex=16; Zindex=64;
% Rindex=28; Zindex=floor(size(M.zgrid,1)/2)+1;
-Indices=M.Rindex(:,1)==Rindex & M.Zindex(:,1)==Zindex;
+Indices=M.Rindex(:,1,false)==Rindex & M.Zindex(:,1,false)==Zindex;
-Vr=M.VR(Indices,1)/M.vlight;
-Vz=M.VZ(Indices,1)/M.vlight;
-Vthet=M.VTHET(Indices,1)/M.vlight;
+Vr=M.VR(Indices,1,false)/M.vlight;
+Vz=M.VZ(Indices,1,false)/M.vlight;
+Vthet=M.VTHET(Indices,1,false)/M.vlight;
-Indicesend=M.Rindex(:,end)==Rindex & M.Zindex(:,end)==Zindex;
+Indicesend=M.Rindex(:,end,false)==Rindex & M.Zindex(:,end,false)==Zindex;
-Vrend=M.VR(Indicesend,end)/M.vlight;
-Vzend=M.VZ(Indicesend,end)/M.vlight;
-Vthetend=M.VTHET(Indicesend,end)/M.vlight;
+Vrend=M.VR(Indicesend,end,false)/M.vlight;
+Vzend=M.VZ(Indicesend,end,false)/M.vlight;
+Vthetend=M.VTHET(Indicesend,end,false)/M.vlight;
binwidth=abs(max(Vrend)-min(Vrend))/10;
f=figure('Name',sprintf("%s v distrib",M.file));
tstudied=0;
legtext=sprintf("t=%2.1f [ns]",M.tpart(end)*1e9);
-subplot(1,3,1)
-h1=histogram(Vr,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
+ax1=subplot(1,3,1);
+%h1=histogram(Vr,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
+h1=histfit(ax1,Vr);
+set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
hold on
-h1=histogram(Vrend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9));
-ylabel('probability')
+%h1=histogram(Vrend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9));
+h1=histfit(ax1,Vrend);
+set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9));
+ylabel('counts')
xlabel('\beta_r')
grid on
ax2=subplot(1,3,2);
-h1=histogram(Vthet(:,1),'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
+%h1=histogram(Vthet(:,1),'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
+h1=histfit(ax2,Vthet(:,1));
+set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
hold on
-h1=histogram(Vthetend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9));
-ylabel('probability')
+%h1=histogram(Vthetend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9));
+h1=histfit(ax2,Vthetend);
+set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9));
+ylabel('counts')
xlabel('\beta_\theta')
legend(ax2,'location','northoutside','orientation','vertical')
grid on
-subplot(1,3,3)
-h1=histogram(Vz(:,1),'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
+ax3=subplot(1,3,3);
+%h1=histogram(Vz(:,1),'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
+h1=histfit(ax3,Vz(:,1));
+set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
hold on
-h1=histogram(Vzend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9));
-ylabel('probability')
+%h1=histogram(Vzend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9));
+h1=histfit(ax3,Vzend);
+set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9));
+ylabel('counts')
xlabel('\beta_z')
grid on
f=gcf;
sgtitle(sprintf('R=%1.2e[m] Z=%1.2e[m] dt=%1.2e[ns]',M.rgrid(Rindex+1),M.zgrid(Zindex+1),M.dt*1e9))
f.PaperOrientation='landscape';
f.PaperUnits='centimeters';
papsize=[16 10];
f.PaperSize=papsize;
[~, name, ~] = fileparts(M.file);
print(f,sprintf('%sParts_V_RZ',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sParts_V_RZ',name))
\ No newline at end of file
diff --git a/matlab/espic2dhdf5.m b/matlab/espic2dhdf5.m
index ffdf7a8..d25e488 100644
--- a/matlab/espic2dhdf5.m
+++ b/matlab/espic2dhdf5.m
@@ -1,499 +1,697 @@
classdef espic2dhdf5
%UNTITLED Summary of this class goes here
% Detailed explanation goes here
properties
filename
name
folder
fullpath
timestamp
info
t0d
t1d
t2d
tpart
it0
it1
it2
%% Physical constants
vlight=299792458;
qe=1.60217662E-19;
me=9.109383E-31;
eps_0=8.85418781762E-12;
kb=1.38064852E-23;
%% Run parameters
dt
nrun
nlres
nlsave
nlclassical
nlPhis
nz
nnr
lz
nplasma
potinn
potout
B0
Rcurv
width
n0
temp
femorder
ngauss
plasmadim
radii
H0
P0
normalized
nbspecies
%% Frequencies
omepe
omece
%% Normalizations
tnorm
rnorm
bnorm
enorm
phinorm
vnorm
%% Grid data
rgrid
zgrid
dz
dr
Vcell
%% Magnetic field
Br
Bz
Athet
B
%% Energies
epot
ekin
etot
etot0
eerr
%% 2D time data
N
fluidUR
fluidUZ
fluidUTHET
pot
Er
Ez
CellVol
Presstens
%% Splines
knotsr
knotsz
%% Particle parameters
weight
qsim
msim
nbparts
- H
- P
partepot
R
Z
Rindex
Zindex
partindex
VR
VZ
VTHET
THET
species
end
methods
function file=file(obj)
file=obj.filename;
end
function obj = espic2dhdf5(filename,readparts)
filedata=dir(filename);
if (isempty(filedata))
error("File: ""%s"" doesn't exist",filename)
end
obj.folder=filedata.folder;
obj.filename=filename;
[~, obj.name, ~] = fileparts(obj.filename);
obj.fullpath=[obj.folder,'/',obj.filename];
obj.timestamp=filedata.date;
if nargin==1
readparts=true;
end
%obj.info=h5info(filename);
obj.dt = h5readatt(obj.fullpath,'/data/input.00/','dt');
obj.nrun = h5readatt(obj.fullpath,'/data/input.00/','nrun');
obj.nlres = strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlres'),'y');
obj.nlsave = strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlsave'),'y');
obj.nlclassical =strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlclassical'),'y');
obj.nlPhis =strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlPhis'),'y');
obj.nz = h5readatt(obj.fullpath,'/data/input.00/','nz');
obj.nnr = h5read(obj.fullpath,'/data/input.00/nnr');
obj.lz = h5read(obj.fullpath,'/data/input.00/lz');
obj.qsim = h5readatt(obj.fullpath,'/data/input.00/','qsim');
obj.msim = h5readatt(obj.fullpath,'/data/input.00/','msim');
try
obj.weight=h5readatt(obj.fullpath,'/data/part/','weight');
catch
obj.weight=obj.msim/obj.me;
end
obj.nplasma = h5readatt(obj.fullpath,'/data/input.00/','nplasma');
obj.potinn = h5readatt(obj.fullpath,'/data/input.00/','potinn');
obj.potout = h5readatt(obj.fullpath,'/data/input.00/','potout');
obj.B0 = h5readatt(obj.fullpath,'/data/input.00/','B0');
obj.Rcurv = h5readatt(obj.fullpath,'/data/input.00/','Rcurv');
obj.width = h5readatt(obj.fullpath,'/data/input.00/','width');
obj.n0 = h5readatt(obj.fullpath,'/data/input.00/','n0');
obj.temp = h5readatt(obj.fullpath,'/data/input.00/','temp');
try
obj.it0 = h5readatt(obj.fullpath,'/data/input.00/','it0d');
obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it2d');
obj.it2 = h5readatt(obj.fullpath,'/data/input.00/','itparts');
catch
obj.it0 = h5readatt(obj.fullpath,'/data/input.00/','it0');
obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it1');
obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it2');
end
try
obj.nbspecies=h5readatt(obj.fullpath,'/data/input.00/','nbspecies');
obj.normalized=strcmp(h5readatt(obj.fullpath,'/data/input.00/','rawparts'),'y');
catch
obj.nbspecies=1;
obj.normalized=false;
end
obj.omepe=sqrt(abs(obj.n0)*obj.qe^2/(obj.me*obj.eps_0));
obj.omece=obj.qe*obj.B0/obj.me;
obj.nbparts= h5read(obj.fullpath, '/data/var0d/nbparts');
try
obj.H0 = h5read(obj.fullpath,'/data/input.00/H0');
obj.P0 = h5read(obj.fullpath,'/data/input.00/P0');
catch
obj.H0=3.2e-14;
obj.P0=8.66e-25;
end
% Normalizations
obj.tnorm=1/obj.omepe;
obj.rnorm=obj.vlight*obj.tnorm;
obj.bnorm=obj.B0;
obj.enorm=obj.vlight*obj.bnorm;
obj.phinorm=obj.enorm*obj.rnorm;
obj.vnorm=obj.vlight;
% Grid data
obj.rgrid= h5read(obj.fullpath, '/data/var1d/rgrid')*obj.rnorm;
obj.zgrid= h5read(obj.fullpath, '/data/var1d/zgrid')*obj.rnorm;
obj.dz=(obj.zgrid(end)-obj.zgrid(1))/double(obj.nz);
rid=1;
for i=1:length(obj.nnr)
obj.dr(i)=(obj.rgrid(sum(obj.nnr(1:i))+1)-obj.rgrid(rid))/double(obj.nnr(i));
rid=rid+obj.nnr(i);
end
obj.Vcell=(obj.zgrid(2:end)-obj.zgrid(1:end-1))*((obj.rgrid(2:end).^2-obj.rgrid(1:end-1).^2)*pi)';
Br = h5read(obj.fullpath,'/data/fields/Br')*obj.bnorm;
obj.Br= reshape(Br,length(obj.zgrid),length(obj.rgrid));
Bz = h5read(obj.fullpath,'/data/fields/Bz')*obj.bnorm;
obj.Bz= reshape(Bz,length(obj.zgrid),length(obj.rgrid));
try
Atheta = h5read(obj.fullpath,'/data/fields/Athet')*obj.bnorm;
obj.Athet= reshape(Atheta,length(obj.zgrid),length(obj.rgrid));
catch
end
obj.B=sqrt(obj.Bz.^2+obj.Br.^2);
clear Br Bz
try
obj.t0d=h5read(obj.fullpath,'/data/var0d/time');
catch
obj.t0d=obj.dt.*double(0:length(obj.epot)-1);
end
obj.femorder = h5read(obj.fullpath,'/data/input.00/femorder');
obj.ngauss = h5read(obj.fullpath,'/data/input.00/ngauss');
obj.plasmadim = h5read(obj.fullpath,'/data/input.00/plasmadim');
obj.radii = h5read(obj.fullpath,'/data/input.00/radii');
obj.epot = h5read(obj.fullpath,'/data/var0d/epot');
obj.ekin = h5read(obj.fullpath,'/data/var0d/ekin');
obj.etot = h5read(obj.fullpath,'/data/var0d/etot');
try
obj.etot0 = h5read(obj.fullpath,'/data/var0d/etot0');
obj.eerr = obj.etot-obj.etot0;
catch
obj.eerr = obj.etot-obj.etot(2);
end
if(obj.normalized)
obj.pot=gridquantity(obj.fullpath,'/data/fields/pot',sum(obj.nnr)+1, obj.nz+1,1);
obj.Er=gridquantity(obj.fullpath,'/data/fields/Er',sum(obj.nnr)+1, obj.nz+1,1);
obj.Ez=gridquantity(obj.fullpath,'/data/fields/Ez',sum(obj.nnr)+1, obj.nz+1,1);
else
obj.pot=gridquantity(obj.fullpath,'/data/fields/pot',sum(obj.nnr)+1, obj.nz+1,obj.phinorm);
obj.Er=gridquantity(obj.fullpath,'/data/fields/Er',sum(obj.nnr)+1, obj.nz+1,obj.enorm);
obj.Ez=gridquantity(obj.fullpath,'/data/fields/Ez',sum(obj.nnr)+1, obj.nz+1,obj.enorm);
end
try
obj.t2d = h5read(obj.fullpath,'/data/fields/time');
catch
info=h5info(obj.fullpath,'/data/fields/partdensity');
obj.t2d=obj.dt*(0:info.objspace.Size(2)-1)*double(obj.it1);
end
try
info=h5info(obj.fullpath,'/data/fields/moments');
obj.femorder = h5read(obj.fullpath,'/data/input.00/femorder');
kr=obj.femorder(2)+1;
obj.knotsr=augknt(obj.rgrid,kr);
kz=obj.femorder(1)+1;
obj.knotsz=augknt(obj.zgrid,kz);
try
obj.CellVol= reshape(h5read(obj.fullpath,'/data/fields/volume'),length(obj.knotsz)-kz,length(obj.knotsr)-kr);
obj.CellVol=permute(obj.CellVol,[2,1,3])*obj.rnorm^3;
catch
zvol=fnder(spmak(obj.knotsz,ones(1,length(obj.knotsz)-kz)), -1 );
rvol=fnder(spmak(obj.knotsr,2*pi*[obj.rgrid' 2*obj.rgrid(end)-obj.rgrid(end-1)]), -1 );
ZVol=diff(fnval(zvol,obj.knotsz));
RVol=diff(fnval(rvol,obj.knotsr));
obj.CellVol=RVol(3:end-1)*ZVol(3:end-1)';
obj.CellVol=padarray(obj.CellVol,[1,1],'replicate','post');
end
if(obj.normalized)
obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, 1, 1);
else
obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, abs(obj.qsim/obj.qe), 1);
end
obj.fluidUR=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 2);
obj.fluidUTHET=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 3);
obj.fluidUZ=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 4);
if(obj.normalized)
obj.Presstens=splinepressure(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.me);
else
obj.Presstens=splinepressure(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.msim);
end
catch
obj.CellVol=(obj.zgrid(2:end)-obj.zgrid(1:end-1))*((obj.rgrid(2:end).^2-obj.rgrid(1:end-1).^2)*pi)';
obj.CellVol=obj.CellVol';
obj.N=griddensity(obj.fullpath, '/data/fields/partdensity', sum(obj.nnr)+1, obj.nz+1, obj.CellVol, abs(obj.qsim/obj.qe), true);
obj.fluidUR=gridquantity(obj.fullpath, '/data/fields/fluidur', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true);
obj.fluidUTHET=gridquantity(obj.fullpath, '/data/fields/fluiduthet', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true);
obj.fluidUZ=gridquantity(obj.fullpath, '/data/fields/fluiduz', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true);
end
if(readparts)
if(obj.normalized)
- obj.R = h5read(obj.fullpath,'/data/part/R');
- obj.Z = h5read(obj.fullpath,'/data/part/Z');
+ obj.R = h5partsquantity(obj.fullpath,'/data/part','R');
+ obj.Z = h5partsquantity(obj.fullpath,'/data/part','Z');
else
- obj.R = h5read(obj.fullpath,'/data/part/R')*obj.rnorm;
- obj.Z = h5read(obj.fullpath,'/data/part/Z')*obj.rnorm;
+ obj.R = h5partsquantity(obj.fullpath,'/data/part','R',obj.rnorm);
+ obj.Z = h5partsquantity(obj.fullpath,'/data/part','Z',obj.rnorm);
end
try
- obj.THET = h5read(obj.fullpath,'/data/part/THET');
+ obj.THET = h5partsquantity(obj.fullpath,'/data/part','THET');
catch
clear obj.THET
end
try
- obj.Rindex=h5read(obj.fullpath,'/data/part/Rindex');
- obj.Zindex=h5read(obj.fullpath,'/data/part/Zindex');
+ obj.Rindex=h5partsquantity(obj.fullpath,'/data/part','Rindex');
+ obj.Zindex=h5partsquantity(obj.fullpath,'/data/part','Zindex');
catch
clear obj.Rindex obj.Zindex
end
- UR = h5read(obj.fullpath,'/data/part/UR');
- UZ = h5read(obj.fullpath,'/data/part/UZ');
- UTHET= h5read(obj.fullpath,'/data/part/UTHET');
- if(obj.nlclassical)
- obj.VR = UR*obj.vnorm;
- clear UR
- obj.VZ = UZ*obj.vnorm;
- clear UZ
- obj.VTHET = UTHET*obj.vnorm;
- clear UTHET
- else
- GAMM = sqrt(1+UR.^2+UZ.^2+UTHET.^2);
- obj.VR = UR./GAMM*obj.vnorm;
- clear UR
- obj.VZ = UZ./GAMM*obj.vnorm;
- clear UZ
- obj.VTHET = UTHET./GAMM*obj.vnorm;
- clear UTHET
- clear GAMM
- end
- %obj.VPERP = sqrt(obj.VR.*obj.VR+obj.VTHET.*obj.VTHET);
+ obj.VR = h5partsquantity(obj.fullpath,'/data/part','UR',obj.vnorm);
+ obj.VZ = h5partsquantity(obj.fullpath,'/data/part','UZ',obj.vnorm);
+ obj.VTHET= h5partsquantity(obj.fullpath,'/data/part','UTHET',obj.vnorm);
+
if(obj.normalized)
- obj.partepot = sign(obj.qsim)*h5read(filename,'/data/part/pot')*obj.qe;
+ obj.partepot = h5partsquantity(obj.fullpath,'/data/part','pot',sign(obj.qsim)*obj.qe);
else
- obj.partepot = sign(obj.qsim)*h5read(filename,'/data/part/pot')*obj.phinorm*obj.qe;
+ obj.partepot = h5partsquantity(obj.fullpath,'/data/part','pot',sign(obj.qsim)*obj.qe*obj.phinorm);
end
- obj.H=0.5*obj.me*(obj.VR.^2+obj.VTHET.^2+obj.VZ.^2)+obj.partepot;
- %clear partepot
- obj.P=obj.R.*(obj.VTHET*obj.me+sign(obj.qsim)*obj.qe*Athet(obj.R,obj.Z));
+
try
- obj.partindex = h5read(obj.fullpath,'/data/part/partindex');
- partindex=obj.partindex;
- [~,Indices]=sort(partindex,'descend');
- for i=1:size(Indices,2)
- obj.H(:,i)=obj.H(Indices(:,i),i);
- obj.P(:,i)=obj.P(Indices(:,i),i);
- obj.R(:,i)=obj.R(Indices(:,i),i);
- obj.Z(:,i)=obj.Z(Indices(:,i),i);
- obj.Rindex(:,i)=obj.Rindex(Indices(:,i),i);
- obj.Zindex(:,i)=obj.Zindex(Indices(:,i),i);
- obj.VR(:,i)=obj.VR(Indices(:,i),i);
- obj.VZ(:,i)=obj.VZ(Indices(:,i),i);
- obj.VTHET(:,i)=obj.VTHET(Indices(:,i),i);
- obj.THET(:,i)=obj.THET(Indices(:,i),i);
- end
+ obj.partindex = h5partsquantity(obj.fullpath,'/data/part/','partindex');
catch
end
- clear partindex;
end
try
obj.tpart = h5read(obj.fullpath,'/data/part/time');
catch
obj.tpart=obj.dt*(0:size(obj.R,2)-1)*double(obj.it2);
end
- function Atheta=Athet(R,Z)
- Atheta=0.5*obj.B0*(R-obj.width/pi*(obj.Rcurv-1)/(obj.Rcurv+1)...
- .*besseli(1,2*pi*R/obj.width).*cos(2*pi*Z/obj.width));
- end
if(obj.nbspecies >1)
obj.species=h5parts.empty;
for i=2:obj.nbspecies
- obj.species(i-1)=h5parts(obj.fullpath,sprintf('/data/part/%2d',i),obj);
+ nbparts=h5read(obj.fullpath,sprintf('%s/Nparts',sprintf('/data/part/%2d',i)));
+ if (nbparts(1)>0)
+ obj.species(i-1)=h5parts(obj.fullpath,sprintf('/data/part/%2d',i),obj);
+ end
end
end
end
+ function Atheta=Atheta(obj,R,Z)
+ Atheta=0.5*obj.B0*(R-obj.width/pi*(obj.Rcurv-1)/(obj.Rcurv+1)...
+ .*besseli(1,2*pi*R/obj.width).*cos(2*pi*Z/obj.width));
+ end
+
+ function quantity=H(obj,indices)
+ %% computes the total energy for the main specie particle indices{1} at time indices{2}
+ if indices{1}==':'
+ p=1:obj.VR.nparts;
+ else
+ p=indices{1};
+ end
+ if indices{2}==':'
+ t=1:length(obj.tpart);
+ else
+ t=indices{2};
+ end
+ if size(indices,1)>2
+ track=indices{3};
+ else
+ track=false;
+ end
+ quantity=0.5*obj.me*(obj.VR(p,t,track).^2+obj.VTHET(p,t,track).^2+obj.VZ(p,t,track).^2)+obj.partepot(p,t,track);
+ end
+
+ function quantity=P(obj,indices)
+ %% computes the canonical angular momentum for the main specie particle indices{1} at time indices{2}
+ if indices{1}==':'
+ p=1:obj.R.nparts;
+ else
+ p=indices{1};
+ end
+ if indices{2}==':'
+ t=1:length(obj.tpart);
+ else
+ t=indices{2};
+ end
+ if size(indices,1)>2
+ track=indices{3};
+ else
+ track=false;
+ end
+ quantity=obj.R(p,t,track).*(obj.VTHET(p,t,track)*obj.me+sign(obj.qsim)*obj.qe*obj.Atheta(obj.R(p,t,track),obj.Z(p,t,track)));
+ end
+
+ function displaypsi(obj,deltat)
+ %% plot the initial and final radial profile at position z=0 and show the normalized enveloppe function Psi
+ % relevant for Davidson annular distribution function
+ f=figure('Name', sprintf('%s Psi',obj.name));
+ zpos=floor(length(obj.zgrid)/2);
+ tinit=1;
+ tend=length(obj.t2d);
+ if(obj.R.nt<2)
+ H0=obj.H0;
+ P0=obj.P0;
+ else
+ H0=mean(obj.H(:,end,false));
+ P0=mean(obj.P(:,end,false));
+ end
+ lw=1.5;
+ Mirrorratio=(obj.Rcurv-1)/(obj.Rcurv+1);
+ locpot=mean(obj.pot(:,zpos,tend-deltat:tend),3);
+ psi=1+obj.qe*locpot(:)/H0-1/(2*obj.me*H0)*(P0./obj.rgrid+obj.qe*0.5*obj.B0.*(obj.rgrid-obj.width/pi*Mirrorratio*cos(2*pi*obj.zgrid(zpos)/obj.width)*besseli(1,2*pi*obj.rgrid/obj.width))).^2;
+ locdens=mean(obj.N(:,zpos,tend-deltat:tend),3);
+ [maxn,In]=max(locdens);%M.N(:,zpos,tinit));
+ plot(obj.rgrid,obj.N(:,zpos,tinit),'bx-','DisplayName',sprintf('t=%1.2f[ns]',obj.t2d(tinit)*1e9),'linewidth',lw)
+ hold on
+ plot(obj.rgrid,locdens,'rx-','DisplayName',sprintf('t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(tend-deltat)*1e9,obj.t2d(tend)*1e9),'linewidth',lw)
+ plot(obj.rgrid(In-2:end),1./obj.rgrid(In-2:end)*maxn*obj.rgrid(In),'DisplayName','N=c*1/r','linewidth',lw)
+ xlabel('r [m]')
+ ylabel('n_e [m^-3]')
+ I=find(psi>0);
+ if (length(I)>1)
+ I=[I(1)-2; I(1)-1; I; I(end)+1; I(end)+2];
+ else
+ I=obj.nnr(1):length(psi);
+ end
+ rq=linspace(obj.rgrid(max(I(1),1)),obj.rgrid(I(end)),500);
+ psiinterp=interp1(obj.rgrid(I),psi(I),rq,'pchip');
+ zeroindices=find(diff(psiinterp>=0),2);
+ maxpsiinterp=max(psiinterp);
+ plot(rq,maxn*psiinterp/abs(maxpsiinterp),'Displayname','normalized \Psi [a.u.]','linewidth',lw)
+ ylim([0 inf])
+ for i=1:length(zeroindices)
+ border=plot([rq(zeroindices(i)) rq(zeroindices(i))],[0 obj.rgrid(In)/obj.rgrid(In-2)*maxn],'k--','linewidth',lw);
+ set(get(get(border,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
+ end
+ legend
+ xlim([0 0.02])
+ grid on
+ title(sprintf('Radial density profile at z=%1.2e[m]',obj.zgrid(zpos)))
+ obj.savegraph(f,sprintf('%sPsi',obj.name),[15 10]);
+ end
+
+ function sref = subsref(obj,s)
+ % obj(i) is equivalent to obj.Data(i)
+ switch s(1).type
+ case '.'
+ if(strcmp(s(1).subs,'H'))
+ sref=H(obj,s(2).subs);
+ elseif(strcmp(s(1).subs,'P'))
+ sref=P(obj,s(2).subs);
+ elseif(strcmp(s(1).subs,'plotEnergy'))
+ plotEnergy(obj);
+ elseif(strcmp(s(1).subs,'trappedParticles'))
+ trappedParticles(obj);
+ elseif(strcmp(s(1).subs,'displayHP'))
+ displayHP(obj,s(2).subs);
+ elseif(strcmp(s(1).subs,'displaypsi'))
+ displaypsi(obj,s(2).subs);
+ elseif(strcmp(s(1).subs,'Wells'))
+ Wells(obj,s(2).subs{1},s(2).subs{2},s(2).subs{3});
+ elseif(strcmp(s(1).subs,'displayaveragetemp'))
+ displayaveragetemp(obj);
+ else
+ sref=builtin('subsref',obj,s);
+ end
+ case '()'
+ sref=subsref(obj(s(1).subs{1}),s(2:end));
+ case '{}'
+ error('MYDataClass:subsref',...
+ 'Not a supported subscripted reference')
+ end
+ end
+
function changed=ischanged(obj)
%ischanged Check if the file has been changed and some data
%must be reloaded
try
filedata=dir(obj.fullpath);
checkedtimestamp=filedata.date;
if (max(checkedtimestamp > obj.timestamp) )
changed=true;
return
end
changed=false;
return
catch
changed=true;
return
end
end
function plotEnergy(obj)
+ %% Plot the time evolution of the system energy and number of simulated macro particles
tmin=2;
tmax=length(obj.ekin);
f=figure('Name', sprintf('%s Energy',obj.name));
subplot(2,1,1)
plot(obj.t0d(tmin:tmax),obj.ekin(tmin:tmax),'o-',...
obj.t0d(tmin:tmax),obj.epot(tmin:tmax),'d-',...
obj.t0d(tmin:tmax),obj.etot(tmin:tmax),'h-',...
obj.t0d(tmin:tmax),obj.eerr(tmin:tmax),'x--')
legend('ekin', 'epot', 'etot','eerr')
xlabel('Time [s]')
ylabel('Energies [J]')
grid on
xlimits=xlim();
subplot(2,1,2)
try
semilogy(obj.t0d(tmin:tmax),abs(obj.eerr(tmin:tmax)./obj.etot0(tmin:tmax)),'h-')
catch
semilogy(obj.t0d(tmin:tmax),abs(obj.eerr(tmin:tmax)/obj.etot(2)),'h-')
end
xlabel('t [s]')
ylabel('E_{err}/E_{tot}')
xlim(xlimits)
grid on
try
yyaxis right
plot(obj.t0d(tmin:tmax),abs(obj.nbparts(tmin:tmax)./obj.nbparts(1)*100),'d--')
ylabel('Nparts %')
%ylim([0 110])
catch
end
obj.savegraph(f,sprintf('%s/%sEnergy',obj.folder,obj.name));
end
function trappedParticles(obj)
+ %% Plot the time evolution of the number of physical trapped particles in the main specie
f=figure('Name', sprintf('%s Trapped particles',obj.name));
plot(obj.t0d,obj.nbparts*obj.weight)
xlabel('t [s]')
ylabel('N particles')
obj.savegraph(f,sprintf('%s/%sntrapped',obj.folder,obj.name));
end
function savegraph(obj, figure, name, papsize)
+ %% Saves the given figure as a pdf and a .fig in the
figure.PaperUnits='centimeters';
if (nargin < 4)
papsize=[14 16];
end
figure.PaperSize=papsize;
print(figure,name,'-dpdf','-fillpage')
savefig(figure,name)
end
function displayHP(obj,tstart)
- if(size(obj.H,2)>=2)
+ %% Plot the histogramm of the total energy and canonical angular momentum at the beginning and
+ % end of the simulation over the full simulation space for the main specie
+ if(iscell(tstart))
+ tstart=cell2mat(tstart);
+ end
+ if(obj.R.nt>=2)
f=figure('Name', sprintf('%s HP',obj.name));
legtext=sprintf("t=%2.1f - %2.1f [ns]",obj.tpart(tstart)*1e9,obj.tpart(end)*1e9);
subplot(1,2,1)
- partsmax=min(obj.nbparts(end),size(obj.H,1));
- H=obj.H(1:obj.nbparts(1),1);
- h1=histogram(H,10,'BinLimits',[min(H(:)) max(H(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(1)*1e9));
+ partsmax=min(obj.nbparts(end),obj.R.nparts);
+ Hloc=H(obj,{1:obj.nbparts(1),1,false});
+ h1=histogram(Hloc,20,'BinLimits',[min(Hloc(:)) max(Hloc(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(1)*1e9));
hold on
- H=mean(obj.H(1:partsmax,tstart:end),2);
+ Hloc=mean(H(obj,{1:partsmax,tstart:obj.R.nt,false}),2);
%,'Binwidth',h1.BinWidth
- h1=histogram(H,20,'BinLimits',[min(H(:)) max(H(:))],'DisplayName',legtext);
+ h1=histogram(Hloc,20,'BinLimits',[min(Hloc(:)) max(Hloc(:))],'DisplayName',legtext);
ylabel('counts')
xlabel('H [J]')
legend
subplot(1,2,2)
- P=obj.P(1:obj.nbparts(1),1);
- h2=histogram(P,10,'BinLimits',[min(P(:)) max(P(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(1)*1e9));
+ Ploc=P(obj,{1:obj.nbparts(1),1,false});
+ h2=histogram(Ploc,50,'BinLimits',[min(Ploc(:)) max(Ploc(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(1)*1e9));
hold on
- P=mean(obj.P(1:partsmax,tstart:end),2);
- histogram(P,50,'BinLimits',[min(P(:)) max(P(:))],'DisplayName',legtext);
+ Ploc=mean(P(obj,{1:partsmax,tstart:obj.R.nt,false}),2);
+ histogram(Ploc,50,'BinLimits',[min(Ploc(:)) max(Ploc(:))],'DisplayName',legtext);
ylabel('counts')
xlabel('P [kg\cdotm^2\cdots^{-1}]')
%clear P
%clear H
legend
%xlim([0.95*h2.BinLimits(1) 1.05*h2.BinLimits(2)])
obj.savegraph(f,sprintf('%s/%sParts_HP',obj.folder,obj.name));
end
end
+
+ function Wells(obj,pot,fieldtime,maxdens)
+ %% Plot the effective potential along the magnetic field lines given the potential pot at time fieldtime
+ rpos=zeros(size(pot));
+ [r,~]=meshgrid(obj.rgrid,obj.zgrid);
+ rathet=(obj.Athet.*r)';
+ poteff=zeros(size(pot));
+ potdepth=zeros(size(pot));
+ for j=1:size(rpos,2)
+ for i=1:size(rpos,1)
+ if(isempty(find(rathet(i,j)>rathet(:,1),1,'last')))
+ rpos(i,j)=1;
+ rinterp=0;
+ else
+ rpos(i,j)=find(rathet(i,j)>rathet(:,1),1,'last');
+ rinterp=(rathet(i,j)-rathet(rpos(i,j),1))./(rathet(rpos(i,j)+1,1)-rathet(rpos(i,j),1));
+ end
+ bmin=obj.B(1,rpos(i,j))+rinterp*(obj.B(1,rpos(i,j)+1)-obj.B(1,rpos(i,j)));
+ potmin=pot(rpos(i,j),1)+rinterp*(pot(rpos(i,j)+1,1)-pot(rpos(i,j),1));
+ poteff(i,j)=-(pot(i,j)-potmin)*bmin/(obj.B(j,i)-bmin);
+ end
+ end
+ zzero=ceil(length(obj.zgrid)/2);
+ for j=1:size(rpos,2)
+ for i=1:size(rpos,1)
+ if(isempty(find(rathet(i,j)<rathet(:,zzero),1,'first')))
+ rpos(i,j)=length(obj.rgrid);
+ else
+ rpos(i,j)=max(find(rathet(i,j)<rathet(:,zzero),1,'first'),2);
+ end
+ rinterp=(rathet(i,j)-rathet(rpos(i,j)-1,zzero))./(rathet(rpos(i,j),zzero)-rathet(rpos(i,j)-1,zzero));
+ potmin=pot(rpos(i,j)-1,zzero)+rinterp*(pot(rpos(i,j),zzero)-pot(rpos(i,j)-1,zzero));
+ potdepth(i,j)=pot(i,j)-potmin;
+ end
+ end
+
+ f=figure('Name', sprintf('%s Wells',obj.name));
+ ax1=gca;
+ h=surface(ax1,obj.zgrid(2:end-1),obj.rgrid(2:end-1),poteff(2:end-1,2:end-1),'edgecolor','none','facecolor','interp');
+ xlim(ax1,[obj.zgrid(1) obj.zgrid(end)])
+ ylim(ax1,[obj.rgrid(1) obj.rgrid(end)])
+ legend
+ xlabel(ax1,'z [m]')
+ ylabel(ax1,'r [m]')
+ if size(fieldtime)<2
+ title(ax1,sprintf('Ion Well t=%1.2g s n_e=%1.2gm^{-3}',obj.t2d(fieldtime),double(maxdens)))
+ figurename=sprintf('%sWell%d',obj.name,fieldtime);
+ else
+ title(ax1,sprintf('Ion Well t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',obj.t2d(fieldtime(1)),obj.t2d(fieldtime(2)),double(maxdens)))
+ figurename=sprintf('%sWell%d_%d',obj.name,fieldtime);
+ end
+ c = colorbar(ax1);
+ c.Label.String= 'Effective Potential [eV]';
+ view(ax1,2)
+ %set(h,'edgecolor','none');
+ grid on;
+ obj.savegraph(f,sprintf('%s/%s',obj.folder,figurename),[16 14]);
+
+ f=figure('Name', sprintf('%s Well depth',obj.name));
+ ax1=gca;
+ h=surface(ax1,obj.zgrid(2:end-1),obj.rgrid(2:end-1),potdepth(2:end-1,2:end-1),'edgecolor','none','facecolor','interp');
+ xlim(ax1,[obj.zgrid(1) obj.zgrid(end)])
+ ylim(ax1,[obj.rgrid(1) obj.rgrid(end)])
+ legend
+ xlabel(ax1,'z [m]')
+ ylabel(ax1,'r [m]')
+ if size(fieldtime)<2
+ title(ax1,sprintf('Potential Well t=%1.2g s n_e=%1.2gm^{-3}',obj.t2d(fieldtime),double(maxdens)))
+ figurename=sprintf('%sPotWell%d',obj.name,fieldtime);
+ else
+ title(ax1,sprintf('Potential Well t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',obj.t2d(fieldtime(1)),obj.t2d(fieldtime(2)),double(maxdens)))
+ figurename=sprintf('%sPotWell%d_%d',obj.name,fieldtime);
+ end
+ c = colorbar(ax1);
+ c.Label.String= 'Potential [eV]';
+ view(ax1,2)
+ %set(h,'edgecolor','none');
+ grid on;
+
+ obj.savegraph(f,sprintf('%s/%s',obj.folder,figurename),[16 14]);
+ end
+
+ function displayaveragetemp(obj)
+ f=figure('Name',sprintf('%s potinn=%f part temperature',obj.name,obj.potinn*obj.phinorm));
+ vr2=obj.VR(:,:);
+ vr2=mean(vr2.^2,1)-mean(vr2,1).^2;
+ vz2=obj.VZ(:,:);
+ vz2=mean(vz2.^2,1)-mean(vz2,1).^2;
+ vthet2=obj.VTHET(:,:);
+ vthet2=mean(vthet2.^2,1)-mean(vthet2,1).^2;
+ plot(obj.tpart,0.5*obj.me*vr2/obj.qe,'displayname','T_r')
+ hold on
+ plot(obj.tpart,0.5*obj.me*vz2/obj.qe,'displayname','T_z')
+ plot(obj.tpart,0.5*obj.me*vthet2/obj.qe,'displayname','T_{thet}')
+ xlabel('time [s]')
+ ylabel('T [eV]')
+ title(sprintf('\\phi_a=%.1f kV \\phi_b=%.1f kV R=%.1f',obj.potinn*obj.phinorm/1e3,obj.potout*obj.phinorm/1e3,obj.Rcurv))
+ legend
+ grid
+ obj.savegraph(f,sprintf('%s/%s_partstemp',obj.folder,obj.name));
+
+ end
end
end
diff --git a/matlab/h5parts.m b/matlab/h5parts.m
index 323e3bf..131c923 100644
--- a/matlab/h5parts.m
+++ b/matlab/h5parts.m
@@ -1,98 +1,139 @@
classdef h5parts
properties
fullpath
R
Z
THET
VR
VZ
VTHET
Rindex
Zindex
partindex
- H
- P
partepot
nlclassical
q
m
weight
tpart
it2
nbparts
rgrid
zgrid
vnorm
+
+ parent
end
methods
function obj=h5parts(filename,hdf5group,parent)
obj.fullpath=filename;
obj.nlclassical=parent.nlclassical;
obj.it2=parent.it2;
obj.rgrid=parent.rgrid;
obj.zgrid=parent.zgrid;
obj.vnorm=parent.vnorm;
- obj.R = h5read(obj.fullpath,sprintf('%s/R',hdf5group));
- obj.Z = h5read(obj.fullpath,sprintf('%s/Z',hdf5group));
- obj.THET = h5read(obj.fullpath,sprintf('%s/THET',hdf5group));
- obj.Rindex=h5read(obj.fullpath,sprintf('%s/Rindex',hdf5group));
- obj.Zindex=h5read(obj.fullpath,sprintf('%s/Zindex',hdf5group));
- UR = h5read(obj.fullpath,sprintf('%s/UR',hdf5group));
- UZ = h5read(obj.fullpath,sprintf('%s/UZ',hdf5group));
- UTHET= h5read(obj.fullpath,sprintf('%s/UTHET',hdf5group));
+ obj.tpart = h5read(obj.fullpath,sprintf('%s/time',hdf5group));
+ obj.nbparts = h5read(obj.fullpath,sprintf('%s/Nparts',hdf5group));
+ obj.parent=parent;
+
+ obj.R = h5partsquantity(obj.fullpath,hdf5group,'R');
+ obj.Z = h5partsquantity(obj.fullpath,hdf5group,'Z');
+ obj.THET = h5partsquantity(obj.fullpath,hdf5group,'THET');
+ obj.Rindex=h5partsquantity(obj.fullpath,hdf5group,'Rindex');
+ obj.Zindex=h5partsquantity(obj.fullpath,hdf5group,'Zindex');
+ obj.VR = h5partsquantity(obj.fullpath,hdf5group,'UR',obj.vnorm);
+ obj.VZ = h5partsquantity(obj.fullpath,hdf5group,'UZ',obj.vnorm);
+ obj.VTHET= h5partsquantity(obj.fullpath,hdf5group,'UTHET',obj.vnorm);
obj.q= h5readatt(obj.fullpath,hdf5group,'q');
obj.m= h5readatt(obj.fullpath,hdf5group,'m');
obj.weight= h5readatt(obj.fullpath,hdf5group,'weight');
- if(obj.nlclassical)
- obj.VR = UR*obj.vnorm;
- clear UR
- obj.VZ = UZ*obj.vnorm;
- clear UZ
- obj.VTHET = UTHET*obj.vnorm;
- clear UTHET
+ obj.partepot = h5partsquantity(filename,hdf5group,'pot',obj.q);
+
+ obj.partindex=h5partsquantity(obj.fullpath,hdf5group,'partindex');
+% try
+% obj.partindex = h5read(obj.fullpath,sprintf('%s/partindex',hdf5group));
+% partindex=obj.partindex;
+% partindex(partindex==-1)=NaN;
+% %[~,Indices]=sort(partindex,'ascend');
+% Indices=obj.partindex;
+% for i=1:size(Indices,2)
+% obj.H(Indices(:,i),i)=obj.H(:,i);
+% obj.P(Indices(:,i),i)=obj.P(:,i);
+%
+% obj.R(Indices(:,i),i)=obj.R(:,i);
+% obj.Z(Indices(:,i),i)=obj.Z(:,i);
+% obj.Rindex(Indices(:,i),i)=obj.Rindex(:,i);
+% obj.Zindex(Indices(:,i),i)=obj.Zindex(:,i);
+% obj.VR(Indices(:,i),i)=obj.VR(:,i);
+% obj.VZ(Indices(:,i),i)=obj.VZ(:,i);
+% obj.VTHET(Indices(:,i),i)=obj.VTHET(:,i);
+% obj.THET(Indices(:,i),i)=obj.THET(:,i);
+% end
+% catch
+% end
+% clear partindex;
+
+ end
+
+ function quantity=H(obj,indices)
+ if indices{1}==':'
+ p=1:obj.VR.nparts;
else
- GAMM = sqrt(1+UR.^2+UZ.^2+UTHET.^2);
- obj.VR = UR./GAMM*obj.vnorm;
- clear UR
- obj.VZ = UZ./GAMM*obj.vnorm;
- clear UZ
- obj.VTHET = UTHET./GAMM*obj.vnorm;
- clear UTHET
- clear GAMM
+ p=indices{1};
end
- obj.partepot = h5read(filename,sprintf('%s/pot',hdf5group))*obj.q;
- obj.H=0.5*obj.m*(obj.VR.^2+obj.VTHET.^2+obj.VZ.^2)+obj.partepot;
- %clear partepot
- obj.P=obj.R.*(obj.VTHET*obj.m+obj.q*Athet(obj.R,obj.Z));
- try
- obj.partindex = h5read(obj.fullpath,sprintf('%s/partindex',hdf5group));
- partindex=obj.partindex;
- [~,Indices]=sort(partindex,'descend');
- for i=1:size(Indices,2)
- obj.H(:,i)=obj.H(Indices(:,i),i);
- obj.P(:,i)=obj.P(Indices(:,i),i);
- obj.R(:,i)=obj.R(Indices(:,i),i);
- obj.Z(:,i)=obj.Z(Indices(:,i),i);
- obj.Rindex(:,i)=obj.Rindex(Indices(:,i),i);
- obj.Zindex(:,i)=obj.Zindex(Indices(:,i),i);
- obj.VR(:,i)=obj.VR(Indices(:,i),i);
- obj.VZ(:,i)=obj.VZ(Indices(:,i),i);
- obj.VTHET(:,i)=obj.VTHET(Indices(:,i),i);
- obj.THET(:,i)=obj.THET(Indices(:,i),i);
- end
- catch
+ if indices{2}==':'
+ t=1:length(obj.tpart);
+ else
+ t=indices{2};
end
- clear partindex;
- obj.tpart = h5read(obj.fullpath,sprintf('%s/time',hdf5group));
- obj.nbparts = h5read(obj.fullpath,sprintf('%s/Nparts',hdf5group));
-
- function Atheta=Athet(R,Z)
- Atheta=0.5*parent.B0*(R-parent.width/pi*(parent.Rcurv-1)/(parent.Rcurv+1)...
- .*besseli(1,2*pi*R/parent.width).*cos(2*pi*Z/parent.width));
+ if size(indices,1)>2
+ track=indices{3};
+ else
+ track=false;
end
+ quantity=0.5*obj.m*(obj.VR(p,t,track).^2+obj.VTHET(p,t,track).^2+obj.VZ(p,t,track).^2)+obj.partepot(p,t,track);
end
+
+ function quantity=P(obj,indices)
+ if indices{1}==':'
+ p=1:obj.R.nparts;
+ else
+ p=indices{1};
+ end
+ if indices{2}==':'
+ t=1:length(obj.tpart);
+ else
+ t=indices{2};
+ end
+ if size(indices,1)>2
+ track=indices{3};
+ else
+ track=false;
+ end
+ quantity=obj.R(p,t,track).*(obj.VTHET(p,t,track)*obj.m+obj.q*obj.parent.Atheta(obj.R(p,t,track),obj.Z(p,t,track)));
+ end
+
+ function sref = subsref(obj,s)
+ % obj(i) is equivalent to obj.Data(i)
+ switch s(1).type
+ case '.'
+ if(strcmp(s(1).subs,'H'))
+ sref=H(obj,s(2).subs);
+ elseif(strcmp(s(1).subs,'P'))
+ sref=P(obj,s(2).subs);
+ else
+ sref=builtin('subsref',obj,s);
+ end
+ case '()'
+ sref=subsref(obj(s(1).subs{1}),s(2:end));
+ case '{}'
+ error('MYDataClass:subsref',...
+ 'Not a supported subscripted reference')
+ end
+ end
+
end
end
\ No newline at end of file
diff --git a/matlab/h5partsquantity.m b/matlab/h5partsquantity.m
new file mode 100644
index 0000000..d545cff
--- /dev/null
+++ b/matlab/h5partsquantity.m
@@ -0,0 +1,85 @@
+classdef h5partsquantity
+
+ properties
+ filename
+ group
+ dataset
+ nparts
+ nt
+ scale
+ end
+
+ methods
+ function obj=h5partsquantity(filename, group, dataset, scale)
+ obj.filename=filename;
+ obj.dataset=dataset;
+ obj.group=group;
+ obj.nparts=h5info(filename,sprintf('%s/%s',group,dataset)).Dataspace.Size(1);
+ obj.nt=h5info(filename,sprintf('%s/%s',group,dataset)).Dataspace.Size(end);
+ if nargin < 4
+ obj.scale=1;
+ else
+ obj.scale=scale;
+ end
+ end
+
+ function ind=end(obj,k,n)
+ switch k
+ case 1
+ ind=obj.nparts;
+ case 2
+ ind=obj.nt;
+ case default
+ error('Invalid number of dimensions');
+ end
+ end
+
+ function sref = subsref(obj,s)
+ % obj(i) is equivalent to obj.Data(i)
+ switch s(1).type
+ case '.'
+ if(strcmp(s(1).subs,'val'))
+ sref = obj.(s(1).subs)(s(2).subs);
+ else
+ sref=builtin('subsref',obj,s);
+ end
+ case '()'
+ sref = val(obj,s(1).subs);
+ case '{}'
+ error('MYDataClass:subsref',...
+ 'Not a supported subscripted reference')
+ end
+ end
+
+ function quantity=val(obj,indices)
+ if indices{1}==':'
+ p=1:obj.nparts;
+ else
+ p=indices{1};
+ end
+ if indices{2}==':'
+ t=1:obj.nt;
+ else
+ t=indices{2};
+ end
+ if size(indices,1)>2
+ track=indices{3};
+ else
+ track=false;
+ end
+ quantity=zeros(obj.nparts,length(t));
+ if track
+ for i=1:length(t)
+ indices = h5read(obj.filename, sprintf('%s/%s',obj.group,'partindex'),[1 t(i)],[Inf 1]);
+ indices(indices<=0)=max(indices)+1;
+ quantity(indices,i) = h5read(obj.filename, sprintf('%s/%s',obj.group,obj.dataset),[1 t(i)],[Inf 1]);
+ end
+ else
+ for i=1:length(t)
+ quantity(:,i) = h5read(obj.filename, sprintf('%s/%s',obj.group,obj.dataset),[1 t(i)],[Inf 1]);
+ end
+ end
+ quantity=quantity(p,:)*obj.scale;
+ end
+ end
+end