diff --git a/matlab/CheckDiocStability.m b/matlab/CheckDiocStability.m
deleted file mode 100644
index f99145d..0000000
--- a/matlab/CheckDiocStability.m
+++ /dev/null
@@ -1,59 +0,0 @@
-t2dlength=size(M.t2d,1);
-fieldstart=max(1,t2dlength-500);%floor(0.95*t2dlength);
-deltat=t2dlength-fieldstart;
-maxl=50;
-nbsave=200;
-
-%[~, rgridend]=min(abs(M.rgrid-0.045));
-rgridend=sum(M.nnr(1:2));
-[~, name, ~] = fileparts(M.file);
-
-dens=mean(M.N(:,:,fieldstart:end),3);
-[R,Z]=meshgrid(M.rgrid,M.zgrid);
-Rinv=1./R;
-Rinv(isnan(Rinv))=0;
-VTHET=mean(M.fluidUTHET(:,:,fieldstart:end),3);
-omegare=(VTHET.*Rinv');
-Er=mean(M.Er(:,:,fieldstart:end),3);
-Ez=mean(M.Ez(:,:,fieldstart:end),3);
-vdrift=(M.Br'.*Ez-Er.*M.Bz')./M.B'.^2;
-omegadrift=(vdrift.*Rinv');
-nbzid=9;
-zindices=linspace(length(M.zgrid)/6,length(M.zgrid)*5/6,nbzid);
-
-for zid=1:nbzid%1:length(zindices)%;%floor(length(M.zgrid)/2);
- zindex=floor(zindices(zid));
- rindex=M.geomweight(:,zindex,1)>=0;
-
- n=dens(rindex,zindex);
- omegar=omegadrift(rindex,zindex);
- rgrid=M.rgrid(rindex);
- Bz=mean(M.Bz(zindex,rindex),2);
- zfolder=sprintf('diocz_%03d',zindex);
- if ~isfolder(zfolder)
- mkdir(zfolder)
- end
- wce=abs(M.qsim/M.msim*Bz);
- wpe2=M.qsim^2/(M.msim*M.eps_0*M.weight)*n;
- wnum=complex(zeros(nbsave,maxl));
- rnumlength=max(length(rgrid)+2,4003);
- phinum=complex(zeros(nbsave,rnumlength,maxl));
- rnum=zeros(rnumlength,maxl);
- parfor l=1:maxl
- [w, phin, r] = diocotronstab(l, rgrid, n, omegar, wce, 1024+512);
-
- % different profiles of dphi
- [~,ind]=sort(abs(imag(w)),'descend');
-
- wnum(:,l)=w(ind(1:nbsave));
- phinum(:,:,l)=padarray(phin(:,ind(1:nbsave))',[0 rnumlength-size(phin,1)], 0, 'post');
- rnum(:,l)=padarray(r,[0 rnumlength-length(r(:))],0, 'post')';
- sprintf('l=%02d, zid=%02d done',l,zid)
- end
- lsteps=1:maxl;
- savecase(zfolder,rnum,wnum,phinum,rgrid,wce,Er,Ez,n,omegare,vdrift,omegadrift,M.geomweight,rgrid(1),rgrid(end),lsteps);
-end
-
-function savecase(folder,r,w,phinum,rgrid,wce,Er,Ez,n,omegare,vdrift,omegadrift,geomweight,a,b,lsteps)
- save(sprintf('%s/results',folder),'r','w','phinum','a','b','lsteps','rgrid','wce','Er','Ez','n','omegare', 'vdrift','omegadrift','geomweight','phinum')
-end
\ No newline at end of file
diff --git a/matlab/CheckDiocStability_post.m b/matlab/CheckDiocStability_post.m
deleted file mode 100644
index 6d524f7..0000000
--- a/matlab/CheckDiocStability_post.m
+++ /dev/null
@@ -1,126 +0,0 @@
-nbzid=15;
-zindices=linspace(length(M.zgrid)/6,length(M.zgrid)*5/6,nbzid);
-zconsid=1:nbzid;
-
-for zid=zconsid%1:length(zindices)%;%floor(length(M.zgrid)/2);
- zindex=floor(zindices(zid));
- zfolder=sprintf('diocz_%03d',zindex);
- if ~isfolder(zfolder) || ~isfile(sprintf('%s/results.mat',zfolder))
- fprintf('Warning: result %s doesn''t exist\n',zfolder)
- continue
- end
- result=load(sprintf('%s/results.mat',zfolder));
- if ~isfield(result,'geomweight')
- result.geomweight=M.geomweight(:,:,1);
- end
- rindex=M.geomweight(:,zindex,1)>=0;
-
- for l=result.lsteps
-
- r=result.r(:,l);
- phin=result.phinum(:,:,l);
- w=result.w(:,l);
- n=result.n;
- omegar=result.omegadrift;
- rgrid=result.rgrid;
- omegar=omegar(rindex,zindex);
- dn=(n(3:end)-n(1:end-2))./(rgrid(3:end)-rgrid(1:end-2));
- a=min(rgrid);
- b=max(rgrid);
- idmin=find(r==a);
- idmax=find(r==b);
- wd=(M.qsim/M.weight)^2*n/M.eps_0/(M.msim/M.weight)/2/result.wce;
- f=figure;
- subplot(4,1,1)
- hold on
- xlabel('R [m]')
- ylabel('n [m^{-3}]')
- plots=gobjects(2,1);
- plots(1)=plot(rgrid,n,'displayname','n');
- yyaxis right
- plots(2)=plot(rgrid,omegar,'displayname','\omega_{re}');
- hold on
- %plots(3)=plot(rgrid,wd,'displayname','\omega_{d}');
- legend(plots,{'n','\omega_r','\omega_d'},'location','northwest')
- ylabel('\omega_{re} [rad\cdots^{-1}]')
- grid minor
-
- ax=gca;
- ax.LineStyleOrder={'-','--','-.',':'};
- subplot(4,1,[2,3])
- hold on
- p=gobjects();
- j=1;
- for i=1:size(phin,1)
- if imag(w(i))/abs(real(w(i)))>1e-4
- freq=w(i);
- disphinum=phin(i,:);
- disphinum=disphinum/max(abs(disphinum));
- yyaxis left
- p(j)=plot(r(idmin:idmax),[0 real(disphinum(idmin:idmax-2)) 0],'displayname',sprintf('\\omega=%1.3g%+1.3gi',real(freq),imag(freq)));
- j=j+1;
- yyaxis right
- plot(r(idmin:idmax),[0 imag(disphinum(idmin:idmax-2)) 0],'displayname',sprintf('\\omega=%1.3g%+1.3gi',real(freq),imag(freq)))
- end
- end
- subplot(4,1,1)
- if (abs(imag(w(1))/real(w(1)))>1e-10)
- title(sprintf('Unstable case l=%d, \\omega_d=%1.3g, z=%1.3f [mm]',l,max(wd),M.zgrid(zindex)*1e3))
- else
- title(sprintf('Stable case l=%d, \\omega_d=%1.3g z=%1.3f [mm]',l,max(wd),M.zgrid(zindex)*1e3))
- end
- subplot(4,1,[2,3])
- yyaxis left
- ylim([-1 1])
- xlabel('R [m]')
- ylabel('real(\delta\phi) [a.u.]')
- grid on
- yyaxis right
- ylim([-1 1])
- grid minor
- ylabel('imag(\delta\phi) [a.u.]')
- f.PaperUnits='centimeters';
- f.PaperSize=[16,20];
- legend('location','northwest')
- if length(p)>=1 && all(ishghandle(p))
- legend(p)
- end
- subplot(4,1,4)
- plot(real(w),imag(w),'x')
- xlabel('\omega_r [rad\cdots^{-1}]')
- ylabel('\omega_i [rad\cdots^{-1}]')
- subplot(4,1,1)
- ax1=gca;
- subplot(4,1,[2,3])
- ax2=gca;
- drawnow
- posit1=get(ax1,'position');
- posit2=get(ax2,'position');
- posit2(3)=posit1(3);
- set(ax2,'position',posit2);
- print(f,sprintf('%s/%s_dioc_l%d',zfolder,M.name,l),'-dpdf','-fillpage')
- savefig(f,sprintf('%s/%s_dioc_l%d',zfolder,M.name,l))
- close(f)
- f=figure;
- subplot(2,1,1)
- hold on
- plot(rgrid,n/max(n),'displayname','n')
- plot(rgrid(2:end-1),dn/max(abs(dn)),'displayname','dn')
- plot(rgrid,(real(w(1))-l*omegar)/max(abs((real(w(1))-l*omegar))),'displayname','w-lwre')
- legend
- grid on
- xlabel('r [m]')
- ylabel('Normalized quantity a.u.')
- subplot(2,1,2)
- disphinum=phin(1,:);
- disphinum=disphinum/max(abs(disphinum));
- p(j)=plot(r(idmin:idmax),[0 real(disphinum(idmin:idmax-2)) 0],'displayname',sprintf('\\omega=%1.3g%+1.3gi',real(freq),imag(freq)));
- xlabel('r [m]')
- ylabel('Re(\delta\phi) a.u.')
- f.PaperUnits='centimeters';
- f.PaperSize=[12,6];
- print(f,sprintf('%s/%s_maxw_l%d',zfolder,M.name,l),'-dpdf','-fillpage')
- savefig(f,sprintf('%s/%s_maxw_l%d',zfolder,M.name,l))
- close(f)
- end
-end
\ No newline at end of file
diff --git a/matlab/CheckDiocStability_post_brd.m b/matlab/CheckDiocStability_post_brd.m
deleted file mode 100644
index 67a9570..0000000
--- a/matlab/CheckDiocStability_post_brd.m
+++ /dev/null
@@ -1,92 +0,0 @@
-nbzid=15;
-zindices=linspace(length(M.zgrid)/6,length(M.zgrid)*5/6,nbzid);
-zconsid=2:2:nbzid-1;
-
-
-dens=mean(M.N(:,:,fieldstart:fieldend),3);
-wd=M.qe*dens./(2*M.eps_0*M.B');
-
-wd=max(wd(:));
-f=figure('Name', sprintf('%s dioc_stab',M.name));
-fieldstart=max(1,length(M.t2d)-500);
-fieldend=length(M.t2d);
-ax1=subplot(2,1,1);
-geomw=M.geomweight(:,:,1);
-dens=mean(M.N(:,:,fieldstart:fieldend),3);
-dens(geomw<=0)=0;
-[C,h]=contourf(ax1,M.zgrid*1e3,M.rgrid*1e3,dens,'LineColor','none');
-set(h,'LineColor','none')
-hold on
-contour(ax1,M.zgrid*1e3,M.rgrid*1e3,M.geomweight(:,:,1),[0 0],'w-.','linewidth',1.5,'Displayname','Boundaries');
-if(M.conformgeom)
- ylim(ax1,[M.rgrid(1) M.rgrid(sum(M.nnr(1:2)))]*1e3)
-else
- ylim(ax1,[M.rgrid(1) M.rgrid(end)]*1e3)
-end
-xlim(ax1,[M.zgrid(1) M.zgrid(end)]*1e3)
-xlabel(ax1,'z [mm]')
-ylabel(ax1,'r [mm]')
-c = colorbar(ax1);
-c.Label.String= 'n[m^{-3}]';
-view(ax1,2)
-grid on;
-hold on;
-
-zpos=zeros(size(zconsid));
-wrslt=zeros(size(zconsid));
-lmax=zeros(size(zconsid));
-for zid=1:length(zconsid)%1:length(zindices)%;%floor(length(M.zgrid)/2);
- zindex=floor(zindices(zconsid(zid)));
- plot(ax1,M.zgrid(zindex)*[1 1]*1e3,[M.rgrid(1) M.rgrid(end)]*1e3,'r--')
- zpos(zid)=M.zgrid(zindex)*1e3;
- zfolder=sprintf('diocz_%03d',zindex);
- if ~isfolder(zfolder) || ~isfile(sprintf('%s/results.mat',zfolder))
- fprintf('Warning: result %s doesn''t exist\n',zfolder)
- continue
- end
- result=load(sprintf('%s/results.mat',zfolder));
- if ~isfield(result,'geomweight')
- result.geomweight=M.geomweight(:,:,1);
- end
- rindex=M.geomweight(:,zindex,1)>=0;
-
- for l=result.lsteps
-
- r=result.r(:,l);
- phin=result.phinum(:,:,l);
- w=result.w(:,l);
- [~,idwmax]=max(imag(w));
- if(imag(w(idwmax))>imag(wrslt(zid)))
- wrslt(zid)=w(idwmax);
- lmax(zid)=l;
- end
- end
-end
-ax2=subplot(2,1,2);
-semilogy(zpos,real(wrslt)/wd,'bx','displayname','Real')
-xlabel('z [mm]')
-hold on
-semilogy(zpos,imag(wrslt)/wd,'b+','displayname','Imag')
-xlim([M.zgrid(1) M.zgrid(end)]*1e3)
-ylabel('\omega/\omega_d')
-ylim([1e-2 2])
-yyaxis right
-plot(zpos,lmax,'d','displayname','l')
-ylabel('most unstable l')
-grid on
-title('Most unstable modes')
-legend('location','east')
-
-
-drawnow
-posit1=get(ax1,'position');
-posit2=get(ax2,'position');
-posit2(3)=posit1(3);
-set(ax2,'position',posit2);
-
-f.PaperOrientation='landscape';
-f.PaperUnits='centimeters';
-papsize=[10 12];
-f.PaperSize=papsize;
-print(f,sprintf('%sdioc_stab',M.name),'-dpdf','-fillpage')
-savefig(f,sprintf('%sdioc_stab',M.name))
\ No newline at end of file
diff --git a/matlab/DiocotronStability.m b/matlab/DiocotronStability.m
deleted file mode 100644
index fcd054c..0000000
--- a/matlab/DiocotronStability.m
+++ /dev/null
@@ -1,41 +0,0 @@
-rbm=0.005;
-rbp=0.0128;
-%rbm=0.09;
-%rbp=0.113;
-%b=0.16;
-%a=0.001;
-a=M.rgrid(1);
-b=M.rgrid(end);
-
-Q=mean(M.Er(1,floor(length(M.zgrid)/2),end-30:end),3)/2*a;
-%Q=mean(M.Er(1,100,end-30:end),3)/2*a;
-omega_d=M.omepe^2/(2*M.omece);
-if( a~=0)
- omega_q=2*Q/(M.B0*(rbm^2));
-else
- omega_q=0;
-end
-
-Levaluated=10;
-stability=zeros(Levaluated,1);
-for l=1:Levaluated
- stability(l)=(-l*(1-omega_q/omega_d)*(1-(rbm/rbp)^2)*(1-(a/b)^(2*l))+ ...
- 2*(1+(a/b)^(2*l))-(1+(rbm/rbp)^(2*l))*((a/rbm)^(2*l)+(rbp/b)^(2*l)))^2 -...
- 4*(rbm/rbp)^(2*l)*(1-(rbp/b)^(2*l))^2*(1-(a/rbm)^(2*l))^2;
- if(stability(l) <0)
- fprintf(2,"Diocotron mode l=%d is potentially unstable\n",l);
- end
-end
-fig=figure;
-plot(1:Levaluated,stability)
-hold on
-plot([1 Levaluated],[0 0],'k')
-for l=1:Levaluated
- if stability(l)<0
- plot(l,stability(l),'rx')
- end
-end
-xlabel('azimuthal mode number')
-ylabel('stability condition')
-title(sprintf('r_b^-=%1.2g [m] r_b^+=%1.2g [m] a=%1.2g [m] b=%1.2g [m]',rbm,rbp,a,b))
-M.savegraph(fig,sprintf('%s/%s_diocstab',M.folder,M.name),[15,10])
\ No newline at end of file
diff --git a/matlab/diocotronstab.m b/matlab/diocotronstab.m
deleted file mode 100644
index 00e5477..0000000
--- a/matlab/diocotronstab.m
+++ /dev/null
@@ -1,84 +0,0 @@
-function [w, phin, r] = diocotronstab(l, R, n, wre, wce, nRmin)
-if nargin<6
- nRmin=512;
-end
-% Necessary physical constants
-% electron charge and mass
- e = 1.602176634e-19;
- me = 9.1093837015e-31;
-% vacuum permittivity
- eps_0 = 8.854188e-12;
-
-% declaration of wpe2 (vector for numerical solution, renormalised on wce)
-
- wpe2 = n*(e^2)/(me*eps_0);
-
-% scale factor for the frequency
- wscale=max(abs(wpe2/(2*wce)));
-
- % construction of r, new vector of distances (chooses the thinner mesh and
- % interpolate all variables on that mesh
-
- delta = diff(R);
- delta=delta(delta>0);
- x = unique(delta);
- N = numel(x);
- count = zeros(N,1);
- for k = 1:N
- count(k) = sum(delta==x(k));
- end
- x(count<2)=[];
- rmin = min(x);
- nR = floor((max(R)-min(R))/rmin) + 1;
- if nR<nRmin
- nR=nRmin+1;
- end
- r = linspace(R(1),R(end),nR);
- b=max(R);
- rmin=rmin/b;
- r=r/b;
- R=R/b;
-
- % interpolation
- % of n, wre and wpe2 --> nin, wrin, wpe2in
- wrin = interp1(R(2:end),wre(2:end),r,'spline');
- wpe2in = interp1(R,wpe2,r);
- wpe2in=wpe2in/wscale^2;
- wrin=wrin/wscale;
- wce=wce/wscale;
-
- % finite difference to compute d/dr(wpe^2)
- Dwpe2 = zeros([nR 1]);
- Dwpe2(2:(nR-1)) = (wpe2in(3:end)-wpe2in(1:(nR-2)))/(2*rmin);
-
-% construction of matrices A & B such that A*phi = w*B*phi
-
- A = zeros(nR,nR);
- B = zeros(nR,nR);
-
-
- % finite differences in A & B + other terms of the eigenvalue
- % equation
-
- for i = 2:1:(nR-1)
- A(i,i) = l*wrin(i)*(-2/(rmin)^2 - l^2/r(i)^2) - l*Dwpe2(i)/wce/r(i);
- A(i,i-1) = l*wrin(i)*(1/(rmin)^2 - 1/(2*r(i)*rmin));
- A(i,i+1) = l*wrin(i)*(1/(rmin)^2 + 1/(2*r(i)*rmin));
-
- B(i,i) = - 2/(rmin)^2 - (l^2/r(i)^2);
- B(i,i-1) = 1/(rmin)^2 - 1/(2*r(i)*rmin);
- B(i,i+1) = 1/(rmin)^2 + 1/(2*r(i)*rmin);
-
- end
- A=A(2:(end-1),2:(end-1));
- B=B(2:(end-1),2:(end-1));
-
-% Solution of the eigenvalue problem
- [phin,w] = eig(A, B);
-
-% renormalisation of w and r
- w=diag(w);
- w = w*wscale;
- r=r*b;
-end
-