3 % STEADY_STATE solves d/dL DLL/L^2 df/dL + f/tau= 0 equation
8 % L1d - 1d L-star array
9 % D1d - 1d radial diffusion coefficient
10 % tau - electron lifetime
13 % PSD_SS - computed 1d PSD
15 % Author: Yuri Shprits
16 % Email: shprits@mit.edu
17 % Last change: 20014-01-01 by Dmitriy Subbotin
18 % The work was supported by:
19 % LANL grant 12-LR-235337, PI Yuri Shprits
20 % NASA grant NNX09AF51G, PI Yuri Shprits
22 Am1 = zeros(size(L1d));
23 A0 = zeros(size(L1d));
24 Ap1 = zeros(size(L1d));
27 tau = repmat(tau, size(L1d));
30 for iL = 2:length(L1d)-1
31 Gdh = 1.0/((L1d(iL+1) - L1d(iL-1))/2) .* (L1d(iL).^2);
32 dfdL_L = ((D1d(iL) + D1d(iL-1))/2)/(L1d(iL) - L1d(iL-1))/((((L1d(iL) + L1d(iL-1))/2).^2));
33 dfdL_R = ((D1d(iL+1) + D1d(iL))/2)/(L1d(iL+1) - L1d(iL))/((((L1d(iL+1) + L1d(iL))/2).^2));
34 Am1(iL) = Gdh * dfdL_L;
35 A0(iL) = -Gdh * (dfdL_L + dfdL_R) - 1./tau(iL);
36 Ap1(iL) = Gdh * dfdL_R;
38 RHS = zeros(size(L1d));
45 A = diag(A0) + diag(Am1(2:end), -1) + diag(Ap1(1:end-1), 1);
48 %PSD_SS = inv(A) * RHS(:);
void steady_state(Matrix1D< double > &f, double tau, double Kp, int nx, Matrix1D< double > &L, double f_bnd_out, double f_bnd_in)