% Roots of Bessel functions, from
% http://mathworld.wolfram.com/BesselFunctionZeros.html
% For more roots see:
% http://www.mathworks.com/matlabcentral/fileexchange/6794-bessel-function-zeros

BesselRoots = ...  % first five roots of the first 5 Bessel functions of the 1st kind
[2.4048,  3.8317,  5.1356,  6.3802,  7.5883,  8.7715; ...
 5.5201,  7.0156,  8.4172,  9.7610, 11.0647, 12.3386; ...
 8.6537, 10.1735, 11.6198, 13.0152, 14.3725, 15.7002; ...
11.7915, 13.3237, 14.7960, 16.2235, 17.6160, 18.9801; ...
14.9309, 16.4706, 17.9598, 19.4094, 20.8269, 22.2178];

BesselExtrema = ...  % first five extrema of the first 5 Bessel functions of the 1st kind
[     0,  1.8412,  3.0542,  4.2012,  5.3175,  6.4156; ...
 3.8317,  5.3314,  6.7061,  8.0152,  9.2824, 10.5199; ...
 7.0156,  8.5363,  9.9695, 11.3459, 12.6819, 13.9872; ...
10.1735, 11.7060, 13.1704, 14.5858, 15.9641, 17.3128; ...
13.3237, 14.8636, 16.3475, 17.7887, 19.1960, 20.5755];


% VALUE cylindrical eigenmodes

% Linear plot as a function of x (theta = 0 and pi)
x = -1:0.01:1;
besselx = zeros(length(x),5);
for n = 0:4
    figure(n+1);
    for i = 1:5
        besselx(:,i) = besselj(n, x*BesselRoots(i,n+1));
    end
    plot(x, besselx);
    hold on
    plot([-1,1],[0,0],'--k')
    hold off
    axis([-1.1 1.1 -1.1 1.1])
    title(['n = ', num2str(n)])
    xlabel('x (units R)')
    ylabel(['J_{', num2str(n), '}((\lambda_{', num2str(n), 'i})^{1/2}\rho)'])
    print('-depsc', ['cylindrical_x_', num2str(n), '.eps'])
end

% Image plot as a function of rho and theta in polar coordinates
lin = -0.99:0.02:0.99; % cartesian linear grid
x = lin' * ones(1, length(lin));
y = x';
rho = abs(x + 1i*y);
theta = angle(x + 1i*y);
outside = rho > 1;  % mask points outside the unit radius to be zeroed
coldhot = 1 - [(31:-1:0)'*[1 1 0]; (0:31)'*[0 1 1]]/32;  % color map
figure(8)
for n = 0:4
    for i = 1:5
        subplot(5, 5, 5*n+i)
        besselni = besselj(n, rho*BesselRoots(i,n+1)) .* cos(n*theta);
        besselni(outside) = 0; % zero outside unit radius
        besselni = besselni / max(max(abs(besselni))); % normalize maximum amplitude
        imagesc(besselni, [-1,1])
        colormap(coldhot)
        axis square
        axis off
%        title(['u_{', num2str(n), num2str(i), '}(\rho, \theta)'])
        title(['\lambda_{', num2str(n), num2str(i), '} = ', num2str(BesselRoots(i,n+1)^2,4), '/R^2'])
    end
end
suptitle('{}^{Cylindrical Value Boundary Conditions}')
print('-depsc', 'cylindrical.eps')

% Same, cartesian reference (box rather than cylindrical boundary
% conditions)
xn = 0.5 * (x+1); % decentered, [0,1]-normalized x
yn = 0.5 * (y+1); % decentered, [0,1]-normalized y
figure(9)
for n = 1:5
    for i = 1:5
        subplot(5, 5, 5*(n-1)+i)
        imagesc(sin(n*pi*xn) .* sin(i*pi*yn), [-1,1])
        colormap(coldhot)
        axis square
        axis off
%        title(['u_{', num2str(n), num2str(i), '}(x, y)'])
        title(['\lambda_{', num2str(n), num2str(i), '} = ', num2str((n^2+i^2)*pi^2,4), '/L^2'])
    end
end
suptitle('{}^{Cartesian Box Value Boundary Conditions}')
print('-depsc', 'cartesian.eps')



% FLUX cylindrical eigenmodes

% Linear plot as a function of x (theta = 0 and pi)
x = -1:0.01:1;
besselx = zeros(length(x),5);
for n = 0:4
    figure(10+n+1);
    for i = 1:5
        besselx(:,i) = besselj(n, x*BesselExtrema(i,n+1));
    end
    plot(x, besselx);
    hold on
    plot([-1,1],[0,0],'--k')
    hold off
    axis([-1.1 1.1 -1.1 1.1])
    title(['n = ', num2str(n)])
    xlabel('x (units R)')
    ylabel(['J_{', num2str(n), '}((\lambda_{', num2str(n), 'i})^{1/2}\rho)'])
    print('-depsc', ['cyl_flux_x_', num2str(n), '.eps'])
end

% Image plot as a function of rho and theta in polar coordinates
lin = -0.99:0.02:0.99; % cartesian linear grid
x = lin' * ones(1, length(lin));
y = x';
rho = abs(x + 1i*y);
theta = angle(x + 1i*y);
outside = rho > 1;  % mask points outside the unit radius to be zeroed
coldhot = 1 - [(31:-1:0)'*[1 1 0]; (0:31)'*[0 1 1]]/32;  % color map
figure(18)
for n = 0:4
    for i = 1:5
        subplot(5, 5, 5*n+i)
        besselni = besselj(n, rho*BesselExtrema(i,n+1)) .* cos(n*theta);
        besselni(outside) = 0; % zero outside unit radius
        besselni = besselni / max(max(abs(besselni))); % normalize maximum amplitude
        imagesc(besselni, [-1,1])
        colormap(coldhot)
        axis square
        axis off
%        title(['u_{', num2str(n), num2str(i), '}(\rho, \theta)'])
        title(['\lambda_{', num2str(n), num2str(i), '} = ', num2str(BesselExtrema(i,n+1)^2,4), '/R^2'])
    end
end
suptitle('{}^{Cylindrical Flux Boundary Conditions}')
print('-depsc', 'cyl_flux.eps')

% Same, cartesian reference (box rather than cylindrical boundary
% conditions)
xn = 0.5 * (x+1); % decentered, [0,1]-normalized x
yn = 0.5 * (y+1); % decentered, [0,1]-normalized y
figure(19)
for n = 0:4
    for i = 0:4
        subplot(5, 5, 5*n+i+1)
        imagesc(cos(n*pi*xn) .* cos(i*pi*yn), [-1,1])
        colormap(coldhot)
        axis square
        axis off
%        title(['u_{', num2str(n), num2str(i), '}(x, y)'])
        title(['\lambda_{', num2str(n), num2str(i), '} = ', num2str((n^2+i^2)*pi^2,4), '/L^2'])
    end
end
suptitle('{}^{Cartesian Box Flux Boundary Conditions}')
print('-depsc', 'cart_flux.eps')
