function m20121024_01_solveFulopHelanderODE()
% This script solves the top line of equation (9) in
% Helander, Phys Plasmas 5, 3999 (1998)
% in the trace impurity limit alpha=0:
%
% d n/d theta = g * [n - b^2 + gamma * (n - ) * b^2]
%
% The same equation appears as equation (27) in
% Fulop & Helander, Phys Plasmas 8, 3305 (2001).
%
% The script works for an arbitrary normalized magnetic field strength
% B(theta), which is specified as a function at the end of this file.
%
% To solve the linear equation, we first isolate the inhomogeneous term on
% the right-hand side:
%
% n + gamma * (n - ) * b^2 - (1/g) * (d n/d theta) = b^2
%
% We discretize n on a uniformly spaced grid in theta, giving a column
% vector representing n on the grid points. Then notice the left-hand side
% of the above equation is a linear operation applied to n, in particular,
% an integro-differential operation. Since the operator is linear, we may
% represent it as a matrix M. For the derivative term, we use a
% spectral differentiation matrix. The integral term
% is performed by noting that integration is essentially
% multiplication of the n column vector by a row vector. The ODE is solved
% by solving the linear system M * n = b^2 with b^2 represented as a column
% vector.
% Physics parameters:
g = 0.7;
gamma = 0.4;
% Number of theta gridpoints for linear system
N = 15;
% Call an external function to make the Fourier spectral differentiation
% matrix for a uniform grid, as well as the grid itself:
[theta, ddtheta] = fourdif(N,1);
dtheta = theta(2)-theta(1);
% To integrate in theta, sum elements and multiply by dtheta.
Bs = B(theta);
fluxSurfaceAverageB2 = sum(Bs.^2)*dtheta/(2*pi);
bs = Bs ./ sqrt(fluxSurfaceAverageB2);
bSquared = bs.^2;
% Now begin construction of the matrix.
% Begin with the identity matrix, representing the n term:
matrix = eye(N);
% Add gamma * n * b^2 term:
matrix = matrix + gamma * diag(bSquared);
% Add -gamma * * b^2 term:
matrix = matrix - gamma * bSquared * (dtheta /(2*pi)* bSquared');
% Add -d/dtheta term:
matrix = matrix - (1/g)*ddtheta;
% We're now finished building the matrix.
% Now the key step: solve the linear system.
n = matrix \ bSquared;
% Plot results:
figure(1)
clf
numRows=2;
numCols=1;
% Number of point to use for interpolation:
NFine = 50;
boundary=1;
[thetaFine,bFine] = m20120528_03_FourierInterpolation(bs, 0, 2*pi, NFine, boundary);
subplot(numRows,numCols,1)
plot(thetaFine,bFine)
hold on
plot(theta,bs,'.r')
xlim([0,2*pi])
xlabel('\theta')
ylabel('b')
[thetaFine,nFine] = m20120528_03_FourierInterpolation(n, 0, 2*pi, NFine, boundary);
subplot(numRows,numCols,2)
plot(thetaFine,nFine)
hold on
plot(theta,n,'.r');
xlim([0,2*pi])
xlabel('\theta')
ylabel('n')
end
function BB = B(theta)
% User-supplied function for the field strength |B|.
% You can normalize B any way you like, since
% the output of this function will be divided by **^{-1/2}.
epsilon = 0.3;
BB = 1./(1 + epsilon * cos(theta));
end**