VERB_code_2.3
steady_state.m
1 function PSD_SS = steady_state(L1d, D1d, tau)
2  %%
3  % STEADY_STATE solves d/dL DLL/L^2 df/dL + f/tau= 0 equation
4  %
5  % PSD_SS = steady_state(L1d, D1d, tau)
6  %
7  % Input:
8  % L1d - 1d L-star array
9  % D1d - 1d radial diffusion coefficient
10  % tau - electron lifetime
11  %
12  % Output:
13  % PSD_SS - computed 1d PSD
14 
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
21 
22  Am1 = zeros(size(L1d));
23  A0 = zeros(size(L1d));
24  Ap1 = zeros(size(L1d));
25 
26  if (length(tau) == 1)
27  tau = repmat(tau, size(L1d));
28  end
29 
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;
37  end
38  RHS = zeros(size(L1d));
39  RHS(1) = 1e-99;
40  RHS(end) = 1;
41 
42  A0(1) = 1;
43  A0(end) = 1;
44 
45  A = diag(A0) + diag(Am1(2:end), -1) + diag(Ap1(1:end-1), 1);
46 
47  PSD_SS = A\RHS(:);
48  %PSD_SS = inv(A) * RHS(:);
49 
50 end
void steady_state(Matrix1D< double > &f, double tau, double Kp, int nx, Matrix1D< double > &L, double f_bnd_out, double f_bnd_in)
Definition: PSD.cpp:2373