Initial commit with the newest version of the code; to be opensourced

+akiles2d/+electrons/+parabolic/+semimaxwellian/moment.m

@@ -0,0 +1,119 @@
+%{
+Computes the (evz,evr,evtheta) moment of the electron distribution
+function at the requested points. 
+
+Trivial cases are forced to be zero. Intermediate variables are stored
+as persistent variables to speed up subsequent calls with the same
+solution and ipoints. 
+
+INPUT
+* data: structure with simulation data
+* solution: current solution structure with fields h,phi,ne00p,npoints
+* evz,evr,evtheta: indices of the moment to compute
+* ipoints: points at which the moment will be computed. Omit or give an
+  empty array to compute at all points 
+
+OUTPUT
+* moment: total (evz,evr,evtheta) moment at the requested points
+* moment1,moment2,moment4: moments for the free, reflected and
+  doubly-trapped subpoppulations
+%}
+function [moment,moment1,moment2,moment4] = moment(data,solution,evz,evr,evtheta,ipoints)
+
+%% Unpack
+nintegrationpoints = data.electrons.nintegrationpoints(:);
+alpha = data.electrons.alpha;
+
+npoints = solution.npoints; 
+h = solution.h(:); % force column
+r = solution.r(:);
+phi = solution.phi(:);
+ne00p = solution.ne00p;
+
+if ~exist('ipoints','var') || isempty(ipoints)
+    ipoints = 1:length(solution.h);
+end
+ipoints = ipoints(:); % force column
+
+nipoints = length(ipoints);
+nE_ = sum(nintegrationpoints);
+
+%% Allocate and zero
+moment1 = ipoints.*0;
+moment2 = moment1;
+moment4 = moment1;
+moment = moment1;
+
+%% Trivial cases
+if mod(evr,2) == 1 || mod(evtheta,2) == 1 
+    return;
+end
+
+%% Obtain integration points and weights; find limits of Pperp' (see Eq. A13)
+% save for future calls with same solution and ipoints
+persistent previous_solution previous_ipoints Pperp_limbwd Pperp_limfwd E_
+if isempty(previous_solution) || ~isequal(previous_solution,solution) || ~isequal(previous_ipoints,ipoints)   
+    Pperp_limbwd = zeros(nipoints,nE_); % allocate
+    Pperp_limfwd = Pperp_limbwd;    
+    E_ = Pperp_limbwd;
+    % prepare E' integration nodes
+    for iipoints = 1:nipoints 
+        i = ipoints(iipoints); 
+        E_transition = 1.5*(phi(i)-phi(end));
+        if isinf(phi(end)) || E_transition <= 0
+            E_transition = 5;
+        end        
+        sep = E_transition/(nintegrationpoints(1)-1); % basic separation between linspaced points
+        f = 1.05; % grow factor of this separation for the second segment
+        E_(iipoints,1:nintegrationpoints(1)) = linspace(0,E_transition,nintegrationpoints(1));
+        E_(iipoints,nintegrationpoints(1)+1:nE_) = E_transition + sep*(f.^[1:nintegrationpoints(2)]-1)/(f-1);
+        % Compute Pperp' limits        
+        if i == npoints % the point is h = infty 
+            continue;
+        end
+        for iE_ = 1:nE_ 
+            Pperp_ = (E_(iipoints,iE_) + phi - phi(i))./(h(i)^2 ./h.^2 - 1) - (r(i)^2/h(i)^4);
+            Pperp_limbwd(iipoints,iE_) = max(0,min(Pperp_(1:i)));
+            Pperp_limfwd(iipoints,iE_) = max(0,max(Pperp_(i+1:npoints))); 
+        end
+        if i == 1 % the point is h = 1
+            Pperp_limbwd(iipoints,1) = Inf;
+        end        
+    end  
+    % Store values in persistent variables for future comparisons
+    previous_solution = solution;
+    previous_ipoints = ipoints;
+end
+
+%% Prepare Integrals
+Hijk = @(Pperp_a,Pperp_b) ... % first integral in Pperp'
+        max(0,gamma((2+evr+evtheta)/2)*(gammainc(Pperp_b,(2+evr+evtheta)/2) - gammainc(Pperp_a,(2+evr+evtheta)/2)));
+dGijk = @(E_a,E_b,Gijka,Gijkb) ... % second integral in E', assuming Gijk varies linearly in each subinterval
+       (gamma((1+evz)/2)*(Gijkb.*E_a-Gijka.*E_b).*(gammainc(E_b,(1+evz)/2) - gammainc(E_a,(1+evz)/2)) + ...
+        gamma((3+evz)/2)*(Gijka-Gijkb).*(gammainc(E_b,(3+evz)/2) - gammainc(E_a,(3+evz)/2))) ./ (E_a-E_b);
+
+%% Integrate moment (see Eq. A10)
+factor = ne00p*2^((evz+evr+evtheta)/2)* ...
+         gamma((1+evr)/2)*gamma((1+evtheta)/2)/gamma(1+(evr+evtheta)/2)/pi^(3/2)* ...
+         exp(phi(ipoints)-r(ipoints)./h(ipoints).^4);
+
+Hijk1 = Hijk(Pperp_limfwd,Pperp_limbwd);
+dGijk1 = dGijk(E_(:,1:nE_-1),E_(:,2:nE_),Hijk1(:,1:nE_-1),Hijk1(:,2:nE_));
+moment1 = factor.*(sum(dGijk1,2) + Hijk1(:,nE_).*gamma((1+evz)/2).*gammainc(E_(:,nE_),(1+evz)/2,'upper')); % adds tail to infinity, assuming Gijk1 == Gijk1(end) from the last E' value onward
+
+if mod(evz,2) ~= 1    
+    Hijk2 = Hijk(0,min(Pperp_limbwd,Pperp_limfwd));
+    dGijk2 = dGijk(E_(:,1:nE_-1),E_(:,2:nE_),Hijk2(:,1:nE_-1),Hijk2(:,2:nE_));
+    moment2 = 2*factor.*sum(dGijk2,2);    
+    Hijk4 = Hijk(Pperp_limbwd,Pperp_limfwd);
+    dGijk4 = dGijk(E_(:,1:nE_-1),E_(:,2:nE_),Hijk4(:,1:nE_-1),Hijk4(:,2:nE_));
+    moment4 = 2*alpha*factor.*sum(dGijk4,2);
+else
+    moment2 = moment1.*0;
+    moment4 = moment1.*0;
+end
+moment = moment1 + moment2 + moment4;
+
+
+
+

+akiles2d/+electrons/+parabolic/Ueff.m

@@ -0,0 +1,15 @@
+%{
+Effective potential for the electron axial motion in the parabolic case.
+
+INPUT
+* h: effective plume radius at the axial position of the point(s)
+* phiz: potential at the axis at those point(s)
+* Jr, ptheta: momenta at those point(s)
+
+OUTPUT
+* phi: potential at the requested points
+%}
+function Ueff = Ueff(h,phiz,Jr,ptheta)
+
+%% Compute Ueff
+Ueff = - phiz + sqrt(2)*(Jr/pi+abs(ptheta))./h;

+akiles2d/+electrons/+parabolic/getbetar.m

@@ -0,0 +1,18 @@
+
+%{
+This function returns betar at a point (h,r), given Jr and ptheta.
+
+INPUT
+* h,r: position of the point(s)
+* Jr, ptheta: momenta at those point(s)
+
+OUTPUT
+* betar: the value of the radial angle coordinate betar at those
+  point(s)
+%}
+function betar = getbetar(h,r,Jr,ptheta)
+
+%% Compute betar
+betar = acos((Jr/pi+abs(ptheta) - sqrt(2)*r.^2 ./h.^2)./ ...
+         sqrt(Jr/pi.*(Jr/pi+2*abs(ptheta))))/(2*pi);
+betar(isnan(betar)) = 0; % For Jr = ptheta = r = 0

+akiles2d/+electrons/+parabolic/getmomenta.m

@@ -0,0 +1,25 @@
+%{
+This function returns the mechanical energy E, and the momenta Jr,
+ptheta from velocity variables at points (h,r).
+
+INPUT
+* h: effective plume radius at the axial position of the point(s)
+* r: radius of the point(s)
+* phiz: potential at the axis at those point(s)
+* vz,vr,vtheta: velocity components at those point(s)
+
+OUTPUT
+* E, Jr, ptheta: mechanical energy and momenta at those point(s)
+%}
+function [E,Jr,ptheta] = getmomenta(h,r,phiz,vz,vr,vtheta)
+
+%% Compute energies and momenta
+% Definition of |ptheta| 
+ptheta = abs(r.*vtheta);
+
+% Energy
+E = (vz.^2+vr.^2+vtheta.^2)/2 - phiz + r.^2 ./h.^4;
+
+% Jr
+Jr = (sqrt(1/2)*(E+phiz-vz.^2/2).*h.^2 - abs(ptheta))*pi;
+

+akiles2d/+electrons/+parabolic/getr.m

@@ -0,0 +1,17 @@
+%{
+This function returns r for a given h, betar, Jr and ptheta. It is the
+inverse function of getbetar
+
+INPUT
+* h: effective plume radius at the axial position of the point(s)
+* betar: radial angle coordinate at those point(s)
+* Jr, ptheta: momenta at those point(s)
+
+OUTPUT
+* r: radius of the point(s)
+%}
+function r = getr(h,betar,Jr,ptheta)
+
+%% Compute r
+r = sqrt((h.^2 ./sqrt(2)).* ...
+    (Jr/pi+abs(ptheta)-cos(2*pi*betar).*sqrt(Jr/pi.*(Jr/pi+2*abs(ptheta)))));

+akiles2d/+electrons/+parabolic/getvelocities.m

@@ -0,0 +1,26 @@
+%{
+This function returns |vz|, |vr|, |vtheta| at a point (h,r), given phiz,
+E, Jr and ptheta. 
+
+INPUT
+* h: effective plume radius at the axial position of the point(s)
+* r: radius of the point(s)
+* phiz: potential at the axis at those point(s)
+* E, Jr, ptheta: mechanical energy and momenta at those point(s)
+
+OUTPUT
+* absvz,absvr,absvtheta: absolute value of the velocities
+%} 
+function [absvz,absvr,absvtheta] = getvelocities(h,r,phiz,E,Jr,ptheta)
+
+%% Compute velocities
+% Axial velocity |vz|
+temp = sqrt(2)*(Jr/pi+abs(ptheta))./h.^2;
+absvz = sqrt((E+phiz-temp)*2);
+
+% Azimuthal velocity |vtheta| (definition of ptheta)
+absvtheta = abs(ptheta./r);
+absvtheta(isnan(absvtheta)) = 0; % vtheta = 0 always at the axis
+
+% Radial velocity |vr|
+absvr = sqrt((temp - (r.^2)./h.^4)*2 - absvtheta.^2);

+akiles2d/+ions/+parabolic/+cold/moment.m

@@ -0,0 +1,54 @@
+%{
+Computes the (evz,evr,evtheta) moment of the ion distribution function
+at the requested points, which must be at the plume axis.  
+
+This is a simple paraxial, cold ion model, where the axial ion velocity
+uzi follows energy conservation and density decays as 1/(uz*h^2).  
+All other ion moments are zero
+
+INPUT
+* data: structure with simulation data
+* solution: current solution structure with fields h,phi
+* evz,evr,evtheta: indices of the moment to compute
+* ipoints: points at which the moment will be computed. Omit or give an
+  empty array to compute at all points 
+
+OUTPUT
+* moment: total (evz,evr,evtheta) moment at the requested points
+* moment1,moment2,moment4: moments for the free, reflected and
+  doubly-trapped subpoppulations (there are only free in this case)
+%}
+function [moment,moment1,moment2,moment4] = moment(data,solution,evz,evr,evtheta,ipoints)
+
+%% Unpack
+chi = data.ions.chi; % dimensionless ion axial velocity at the origin
+mu = data.ions.mu; % dimensionless ion mass
+
+h = solution.h(:);
+r = solution.r(:);
+phi = solution.phi(:);
+
+if ~exist('ipoints','var') || isempty(ipoints) % An empty ipoints means all points
+    ipoints = 1:length(h);
+end
+ipoints = ipoints(:); % force column
+
+%% Allocate
+moment = ipoints.*0;
+moment1 = moment;
+moment2 = moment;
+moment4 = moment;
+
+%% Throw error if r ~= 0
+if any(r(ipoints)~=0)
+    error('Error: akiles2d.ions.parabolic.cold.moment cannot currently compute moments at points not on the plume axis')
+end
+
+%% Compute moment 
+uz = sqrt(chi^2-(2/mu)*phi(ipoints)); % ion velocity
+ni = chi./(uz.*h(ipoints).^2);
+
+if evr == 0 && evtheta == 0 
+    moment = ni.*uz.^evz; 
+end
+moment1 = moment; % all ions are free

+akiles2d/+postprocessor/moments.m

@@ -0,0 +1,96 @@
+%{
+Postprocesses a converged simulation to compute several moments electrons
+and ions, and standard quantities derived from them.
+
+INPUT
+* data: structure with simulation data
+* solution: converged solution structure
+
+OUTPUT
+* solution: updated structure with the moment fields (electrons, ions)
+%}
+function solution = moments(data,solution)
+
+%% Basic electron moments
+[solution.electrons.M000,solution.electrons.M0001,solution.electrons.M0002,solution.electrons.M0004] = ...
+    akiles2d.electrons.(data.potential.model).(data.electrons.model).moment(data,solution,0,0,0); 
+[solution.electrons.M100,solution.electrons.M1001,solution.electrons.M1002,solution.electrons.M1004] = ...
+    akiles2d.electrons.(data.potential.model).(data.electrons.model).moment(data,solution,1,0,0); 
+[solution.electrons.M200,solution.electrons.M2001,solution.electrons.M2002,solution.electrons.M2004] = ...
+    akiles2d.electrons.(data.potential.model).(data.electrons.model).moment(data,solution,2,0,0); 
+[solution.electrons.M020,solution.electrons.M0201,solution.electrons.M0202,solution.electrons.M0204] = ...
+    akiles2d.electrons.(data.potential.model).(data.electrons.model).moment(data,solution,0,2,0); 
+[solution.electrons.M002,solution.electrons.M0021,solution.electrons.M0022,solution.electrons.M0024] = ...
+    akiles2d.electrons.(data.potential.model).(data.electrons.model).moment(data,solution,0,0,2); 
+[solution.electrons.M300,solution.electrons.M3001,solution.electrons.M3002,solution.electrons.M3004] = ...
+    akiles2d.electrons.(data.potential.model).(data.electrons.model).moment(data,solution,3,0,0); 
+[solution.electrons.M120,solution.electrons.M1201,solution.electrons.M1202,solution.electrons.M1204] = ...
+    akiles2d.electrons.(data.potential.model).(data.electrons.model).moment(data,solution,1,2,0); 
+[solution.electrons.M102,solution.electrons.M1021,solution.electrons.M1022,solution.electrons.M1024] = ...
+    akiles2d.electrons.(data.potential.model).(data.electrons.model).moment(data,solution,1,0,2); 
+
+%% Derived electron quantities
+ll = {'','1','2','4'}; % Label for the whole population and each subpopulation
+
+for il = 1:length(ll)
+    l = ll{il};
+    % Density
+    solution.electrons.(['n',l]) = solution.electrons.(['M000',l]);  
+    % Axial velocity
+    solution.electrons.(['u',l]) = solution.electrons.(['M100',l])./solution.electrons.(['n',l]);  
+    % Axial temperature
+    solution.electrons.(['Tz',l]) = solution.electrons.(['M200',l])./solution.electrons.(['n',l]) - solution.electrons.(['u',l]).^2;  
+    % Radial temperature
+    solution.electrons.(['Tr',l]) = solution.electrons.(['M020',l])./solution.electrons.(['n',l]);
+    % Azimuthal temperature
+    solution.electrons.(['Ttheta',l]) = solution.electrons.(['M002',l])./solution.electrons.(['n',l]);
+    % Axial flux of axial heat
+    solution.electrons.(['qzz',l]) = 1/2*solution.electrons.(['M300',l]) - 3/2*solution.electrons.(['n',l]).*solution.electrons.(['u',l]).*solution.electrons.(['Tz',l]) - 1/2*solution.electrons.(['n',l]).*solution.electrons.(['u',l]).^3;
+    % Axial flux of radial heat
+    solution.electrons.(['qzr',l]) = 1/2*solution.electrons.(['M120',l]) - 3/2*solution.electrons.(['n',l]).*solution.electrons.(['u',l]).*solution.electrons.(['Tr',l]);
+    % Axial flux of azimuthal heat
+    solution.electrons.(['qztheta',l]) = 1/2*solution.electrons.(['M102',l]) - 3/2*solution.electrons.(['n',l]).*solution.electrons.(['u',l]).*solution.electrons.(['Ttheta',l]);
+end
+
+%% Basic ion moments
+[solution.ions.M000,solution.ions.M0001,solution.ions.M0002,solution.ions.M0004] = ...
+    akiles2d.ions.(data.potential.model).(data.ions.model).moment(data,solution,0,0,0); 
+[solution.ions.M100,solution.ions.M1001,solution.ions.M1002,solution.ions.M1004] = ...
+    akiles2d.ions.(data.potential.model).(data.ions.model).moment(data,solution,1,0,0); 
+[solution.ions.M200,solution.ions.M2001,solution.ions.M2002,solution.ions.M2004] = ...
+    akiles2d.ions.(data.potential.model).(data.ions.model).moment(data,solution,2,0,0); 
+[solution.ions.M020,solution.ions.M0201,solution.ions.M0202,solution.ions.M0204] = ...
+    akiles2d.ions.(data.potential.model).(data.ions.model).moment(data,solution,0,2,0); 
+[solution.ions.M002,solution.ions.M0021,solution.ions.M0022,solution.ions.M0024] = ...
+    akiles2d.ions.(data.potential.model).(data.ions.model).moment(data,solution,0,0,2); 
+[solution.ions.M300,solution.ions.M3001,solution.ions.M3002,solution.ions.M3004] = ...
+    akiles2d.ions.(data.potential.model).(data.ions.model).moment(data,solution,3,0,0); 
+[solution.ions.M120,solution.ions.M1201,solution.ions.M1202,solution.ions.M1204] = ...
+    akiles2d.ions.(data.potential.model).(data.ions.model).moment(data,solution,1,2,0); 
+[solution.ions.M102,solution.ions.M1021,solution.ions.M1022,solution.ions.M1024] = ...
+    akiles2d.ions.(data.potential.model).(data.ions.model).moment(data,solution,1,0,2); 
+
+%% Derived ion quantities
+ll = {'','1','2','4'}; % Label for the whole population and each subpopulation
+
+for il = 1:length(ll)
+    l = ll{il};
+    % Density
+    solution.ions.(['n',l]) = solution.ions.(['M000',l]);  
+    % Axial velocity
+    solution.ions.(['u',l]) = solution.ions.(['M100',l])./solution.ions.(['n',l]);  
+    % Axial temperature
+    solution.ions.(['Tz',l]) = solution.ions.(['M200',l])./solution.ions.(['n',l]) - solution.ions.(['u',l]).^2;  
+    % Radial temperature
+    solution.ions.(['Tr',l]) = solution.ions.(['M020',l])./solution.ions.(['n',l]);
+    % Azimuthal temperature
+    solution.ions.(['Ttheta',l]) = solution.ions.(['M002',l])./solution.ions.(['n',l]);
+    % Axial flux of axial heat
+    solution.ions.(['qzz',l]) = 1/2*solution.ions.(['M300',l]) - 3/2*solution.ions.(['n',l]).*solution.ions.(['u',l]).*solution.ions.(['Tz',l]) - 1/2*solution.ions.(['n',l]).*solution.ions.(['u',l]).^3;
+    % Axial flux of radial heat
+    solution.ions.(['qzr',l]) = 1/2*solution.ions.(['M120',l]) - 3/2*solution.ions.(['n',l]).*solution.ions.(['u',l]).*solution.ions.(['Tr',l]);
+    % Axial flux of azimuthal heat
+    solution.ions.(['qztheta',l]) = 1/2*solution.ions.(['M102',l]) - 3/2*solution.ions.(['n',l]).*solution.ions.(['u',l]).*solution.ions.(['Ttheta',l]);
+end
+
+

+akiles2d/+potential/+parabolic/phi.m

@@ -0,0 +1,15 @@
+%{
+Computes phi in the parabolic potential case.
+
+INPUT
+* h: effective plume radius at the axial position of the point(s)
+* r: radius of the point(s)
+* phiz: potential at the axis at those point(s)
+
+OUTPUT
+* phi: potential at the requested points
+%}
+function phi = phi(h,r,phiz)
+
+%% Compute phi
+phi = - r.^2 ./ h.^4 + phiz;