function [] = shells(jsize, dr, rmax, graphradius, A0, A1, nmajor, nminorpermajor)

testrun = 0;            % 0 means normal operation which takes much time and writes files to disk.  
                        % 1 means quick operation (no projected densities) and does not write files to disk.

% Example command line:
% ripples(500, .1, 100000, 60000, 1, -1, 5, 10);

nviews = 30;             % This is the number of points of view of the dark matter density.
                        % The plots of the dark matter density and its potential are jsize x jsize.
jjsize = 1000;          % jjsize is the number of pixels for the potential along the x-axis.

mu0 = 8.7e-13;            % not important - makes no difference - carryover from previous code.
patternperiod = 50000000; % not important - makes no difference for t=0 - carryover from previous code.

% dr is the radial stepsize, and rmax is the radius at which the dark
% matter is cutoff.

% wavelength0 and wavelength1 are the wavelengths of the zeroth and first
% spherical harmonic functions f0 and f1.

% The first peak in the rotation curve (for the case of a spiral galaxy) 
% due to dark matter is at 4.5/2pi = 0.7 times wavelength0.

% patternperiod is the time for the rotating potential which results from
% the interference pattern to go around once.  Positive values give counterclockwise 
% rotation in the xy plane, and negative values give clockwise rotation.  

% mu0 is the dark matter density corresponding to the standard scalar wave with amplitude one.

% graphradius is the radius of visible region in the figures.

% Begin Program
   
%medianr = initialradius;
%graphradius = medianr*nmedians4graphs;
%nfigs = size(dotsizes,2);

% Compute k0, k1, wavelength0, wavelength1

majorspacing = rmax/nmajor;
kdiff = 2*pi/majorspacing;
ksum = nminorpermajor*2*kdiff;
k0 = (ksum+kdiff)/2
k1 = (ksum-kdiff)/2

wavelength0 = 2*pi/k0
wavelength1 = 2*pi/k1

% Computer k0, k1, w0, w1, Omega.

%k0 = 2*pi / wavelength0;
%k1 = 2*pi / wavelength1;

w0 = 2*patternperiod/8/pi * (k1^2 - k0^2) - 2*pi/patternperiod;     %inserted first 2
w1 = 2*patternperiod/8/pi * (k1^2 - k0^2) + 2*pi/patternperiod;     %inserted first 2
dw = w1-w0;

Omega0Squared = (2*pi/(2*patternperiod))^2 + (2*patternperiod/8/pi)^2 * (k1^2 - k0^2)^2 - (k1^2 + k0^2)/2;  %doubled patternperiod everywhere

if Omega0Squared < 0
    
    Omega0Squared
    ErrorMessage = 'Omega0Squared is negative!'
    
end

Omega0 = sqrt(Omega0Squared);
OneOverOmega0 = 1/Omega0
        
% Initialize f0, f1, U0, U2, UU1, W0, W2, and WW1 (which stands for W tilde sub 1).

nsteps = 1 + round(rmax/dr);
dr = rmax / (nsteps - 1);
rvalues = 0:dr:rmax;

f0 = zeros(1, nsteps);
f1 = zeros(1, nsteps);

U0 = zeros(1, nsteps);
U2 = zeros(1, nsteps);
UU1 = zeros(1, nsteps);

W0 = zeros(1, nsteps);
W0p = zeros(1, nsteps);

W2 = zeros(1, nsteps);
W2p = zeros(1, nsteps);

WW1 = zeros(1, nsteps);
WW1p = zeros(1, nsteps);

% Compute f0 and f1.

f0(1) = 1;
f0p = 0;
f0pp = -k0^2 * f0(1) / 3;

f1(1) = 1;
f1p = 0;
f1pp = -k1^2 * f1(1) / 5;   %Changed 7 to 5

for i = 2:nsteps
    
    r = (i-1)*dr;
    
    f0(i) = f0(i-1) + f0p * dr + f0pp * dr^2 / 2;
    f0p = f0p + f0pp * dr;
    f0pp = -k0^2 * f0(i) - 2 * f0p / r;
    
    f1(i) = f1(i-1) + f1p * dr + f1pp * dr^2 / 2;
    f1p = f1p + f1pp * dr;
    f1pp = -k1^2 * f1(i) - 4 * f1p / r;    %Changed 6 to 4
    
end

% Normalize f.

f0 = -f0 / min(f0 .* rvalues);
f1 = -f1 / min(f1 .* rvalues .^ 2);    %Changed ^3 to ^2

fac = max(f0);

f0 = f0 / fac;
f1 = f1 / fac;

% Plot f.

fmin = min(min(f0),min(f1 .* rvalues));
fmax = max(max(f0),max(f1 .* rvalues));

graphsteps = min(1+round(graphradius/dr),size(rvalues,2));
Gradius = dr*(graphsteps-1);
    
% Setup graphics names for initial plots.  

fn = [num2str(jsize) '_' num2str(dr) '_' num2str(rmax) '_' num2str(graphradius) '_' num2str(A0) '_' num2str(A1) '_' num2str(nmajor) '_' num2str(nminorpermajor) '_' num2str(wavelength0) '_' num2str(wavelength1)];
fn = nodecimals(fn);

if 1==1
    figure(1)
    plot(rvalues(1:graphsteps), f0(1:graphsteps))
    hold on
    plot(rvalues(1:graphsteps), f1(1:graphsteps) .* rvalues(1:graphsteps) .^ 1)  %Changed ^2 to ^1
    hold off
    title('f_0(r) and f_1(r) * r^1')  %Changed ^2 to ^1
    axis([0 Gradius fmin fmax])
    fname = ['f1_' fn '.jpg'];     
    if testrun == 0
        saveas(1,fname)
    end
end

% Compute U0, U2, UU1.

for i = 1:nsteps
    
    r = (i-1)*dr;
    
    U0(i) = A0^2 * f0(i)^2 + 2/3 * A1^2 * r^2 * f1(i)^2;   %updated
    U2(i) = - 1/3 * A1^2 * f1(i)^2;                        %updated 
    UU1(i) = 2*A0*A1 * f0(i) * f1(i);                      %no changed needed
       
end

%------------------------------------------------------------------------

% Choose points of view.

xaxislist = zeros(nviews,3);
yaxislist = zeros(nviews,3);
zaxislist = zeros(nviews,3);
labellist = cell(nviews,1);

for i = 1:nviews
    
    xaxis = randomunitvector;
    yaxis = cross(xaxis,randomunitvector);
    yaxis = yaxis / norm(yaxis);
    zaxis = cross(xaxis,yaxis);

    xaxislist(i,:) = xaxis;
    yaxislist(i,:) = yaxis;
    zaxislist(i,:) = zaxis;

    labellist(i) = cellstr('');
end

xaxislist(1,:) = [1 0 0];
yaxislist(1,:) = [0 1 0];
zaxislist(1,:) = [0 0 1];
labellist(1) = cellstr(' in XY Plane');

xaxislist(2,:) = [1  0  0];
yaxislist(2,:) = [0  0  1];
zaxislist(2,:) = [0 -1  0];
labellist(2) = cellstr(' in XZ Plane');
  
xaxislist(3,:) = [0 1 0];
yaxislist(3,:) = [0 0 1];
zaxislist(3,:) = [1 0 0];
labellist(3) = cellstr(' in YZ Plane');
  
% Compute cross sectional densities and 3D density projected to plane (picture if DM radiated slightly).

for i = 1:nviews    
    
    xaxis = xaxislist(i,:);
    yaxis = yaxislist(i,:);
    zaxis = zaxislist(i,:);
    
    CrossDensity = zeros(jsize);
    CrossDensityR2 = zeros(jsize);
    ProjDensity = zeros(jsize);
    ProjDensityR1 = zeros(jsize);
    
    for jrow = 1:jsize
        for jcol = 1:jsize
        
            v1 = (jcol-1)/(jsize-1) * 2 * graphradius - graphradius;
            v2 = (jrow-1)/(jsize-1) * 2 * graphradius - graphradius;
        
            xyz = v1*xaxis + v2*yaxis;
            CrossDensity(jsize+1-jrow,jcol) = DMdensity(xyz(1), xyz(2), xyz(3), 0);
            CrossDensityR2(jsize+1-jrow,jcol) = DMdensity(xyz(1), xyz(2), xyz(3), 0) * (v1^2 + v2^2);
        
            totaldensity = 0;
            dz = 2*graphradius/(jsize-1);
            nz = 1 + round(rmax/dz); 

            if testrun == 0
                totaldensity = 0;
                for k = -nz:nz
                    v3 = k*dz;
                    newxyz = xyz + v3*zaxis;
                    totaldensity = totaldensity + DMdensity(newxyz(1), newxyz(2), newxyz(3), 0);
                end
                ProjDensity(jsize+1-jrow,jcol) = totaldensity;
                ProjDensityR1(jsize+1-jrow,jcol) = totaldensity * sqrt(v1^2 + v2^2);            
            end
            
        end
    end
    
    % Plot cross sectional densities - 1.
    
    figurenumber = 100*i+21;
    figure(figurenumber)
    colormap(gray);
    maximum = 1 * sum(sum(CrossDensity))/jsize/jsize;
    imagesc(CrossDensity, [0 maximum]);
    title(['Cross Sectional Dark Matter Density' char(labellist(i))])

    redvector = [1 0 0];
    xxx = dot(redvector, xaxis);
    yyy = dot(redvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'r', 'LineWidth', 4)

    greenvector = [0 1 0];
    xxx = dot(greenvector, xaxis);
    yyy = dot(greenvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'g', 'LineWidth', 4)

    bluevector = [0 0 1];
    xxx = dot(bluevector, xaxis);
    yyy = dot(bluevector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'b', 'LineWidth', 4)
    
    drawnow
    fname = ['f' num2str(figurenumber) '_' fn '.jpg'];     
    if testrun == 0
        saveas(figurenumber,fname)
    end    
         
    % Plot cross sectional densities - 2.
    
    figurenumber = 100*i+22;
    figure(figurenumber)
    colormap(gray);
    maximum = 2 * sum(sum(CrossDensity))/jsize/jsize;
    imagesc(CrossDensity, [0  maximum]);
    title(['Cross Sectional Dark Matter Density' char(labellist(i))])

    redvector = [1 0 0];
    xxx = dot(redvector, xaxis);
    yyy = dot(redvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'r', 'LineWidth', 4)

    greenvector = [0 1 0];
    xxx = dot(greenvector, xaxis);
    yyy = dot(greenvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'g', 'LineWidth', 4)

    bluevector = [0 0 1];
    xxx = dot(bluevector, xaxis);
    yyy = dot(bluevector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'b', 'LineWidth', 4)
    
    drawnow
    fname = ['f' num2str(figurenumber) '_' fn '.jpg'];     
    if testrun == 0
        saveas(figurenumber,fname)
    end 
    
    % Plot cross sectional densities - 3.
    
    figurenumber = 100*i+23;
    figure(figurenumber)
    colormap(gray);
    maximum = 4 * sum(sum(CrossDensity))/jsize/jsize;
    imagesc(CrossDensity, [0  maximum]);
    title(['Cross Sectional Dark Matter Density' char(labellist(i))])

    redvector = [1 0 0];
    xxx = dot(redvector, xaxis);
    yyy = dot(redvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'r', 'LineWidth', 4)

    greenvector = [0 1 0];
    xxx = dot(greenvector, xaxis);
    yyy = dot(greenvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'g', 'LineWidth', 4)

    bluevector = [0 0 1];
    xxx = dot(bluevector, xaxis);
    yyy = dot(bluevector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'b', 'LineWidth', 4)
    
    drawnow
    fname = ['f' num2str(figurenumber) '_' fn '.jpg'];     
    if testrun == 0
        saveas(figurenumber,fname)
    end 
        
    % Plot cross sectional densities - 4.
    
    figurenumber = 100*i+24;
    figure(figurenumber)
    colormap(gray);
    maximum = 8 * sum(sum(CrossDensity))/jsize/jsize;
    imagesc(CrossDensity, [0  maximum]);
    title(['Cross Sectional Dark Matter Density' char(labellist(i))])

    redvector = [1 0 0];
    xxx = dot(redvector, xaxis);
    yyy = dot(redvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'r', 'LineWidth', 4)

    greenvector = [0 1 0];
    xxx = dot(greenvector, xaxis);
    yyy = dot(greenvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'g', 'LineWidth', 4)

    bluevector = [0 0 1];
    xxx = dot(bluevector, xaxis);
    yyy = dot(bluevector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'b', 'LineWidth', 4)
    
    drawnow
    fname = ['f' num2str(figurenumber) '_' fn '.jpg'];     
    if testrun == 0
        saveas(figurenumber,fname)
    end 
    
    % Plot cross sectional densities - 5.
    
    figurenumber = 100*i+25;
    figure(figurenumber)
    colormap(gray);
    maximum = 16 * sum(sum(CrossDensity))/jsize/jsize;
    imagesc(CrossDensity, [0  maximum]);
    title(['Cross Sectional Dark Matter Density' char(labellist(i))])

    redvector = [1 0 0];
    xxx = dot(redvector, xaxis);
    yyy = dot(redvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'r', 'LineWidth', 4)

    greenvector = [0 1 0];
    xxx = dot(greenvector, xaxis);
    yyy = dot(greenvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'g', 'LineWidth', 4)

    bluevector = [0 0 1];
    xxx = dot(bluevector, xaxis);
    yyy = dot(bluevector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'b', 'LineWidth', 4)
    
    drawnow
    fname = ['f' num2str(figurenumber) '_' fn '.jpg'];     
    if testrun == 0
        saveas(figurenumber,fname)
    end 
    
    % Plot cross sectional densities times r^2.
    
    figurenumber = 100*i+20;
    figure(figurenumber)
    colormap(gray);
    maximum = max(max(CrossDensityR2));
    imagesc(CrossDensityR2, [0  maximum]);
    title(['Cross Sectional Dark Matter Density' char(labellist(i)) ' (Times r^2)'])
   
    redvector = [1 0 0];
    xxx = dot(redvector, xaxis);
    yyy = dot(redvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'r', 'LineWidth', 4)

    greenvector = [0 1 0];
    xxx = dot(greenvector, xaxis);
    yyy = dot(greenvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'g', 'LineWidth', 4)

    bluevector = [0 0 1];
    xxx = dot(bluevector, xaxis);
    yyy = dot(bluevector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'b', 'LineWidth', 4)
    
    drawnow
    fname = ['f' num2str(figurenumber) '_' fn '.jpg'];     
    if testrun == 0
        saveas(figurenumber,fname)
    end 
    
    % Plot projected densities - 1.
    
    if testrun == 0
    
    figurenumber = 100*i+31;
    figure(figurenumber)
    colormap(gray);
    maximum = 1 * sum(sum(ProjDensity))/jsize/jsize;
    imagesc(ProjDensity, [0  maximum]);
    title(['Projected Dark Matter Density' char(labellist(i))])
    
    redvector = [1 0 0];
    xxx = dot(redvector, xaxis);
    yyy = dot(redvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'r', 'LineWidth', 4)

    greenvector = [0 1 0];
    xxx = dot(greenvector, xaxis);
    yyy = dot(greenvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'g', 'LineWidth', 4)

    bluevector = [0 0 1];
    xxx = dot(bluevector, xaxis);
    yyy = dot(bluevector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'b', 'LineWidth', 4)
    
    drawnow
    fname = ['f' num2str(figurenumber) '_' fn '.jpg'];     
    if testrun == 0
        saveas(figurenumber,fname)
    end 
    
    % Plot projected densities - 2.
    
    figurenumber = 100*i+32;
    figure(figurenumber)
    colormap(gray);
    maximum = 2 * sum(sum(ProjDensity))/jsize/jsize;
    imagesc(ProjDensity, [0  maximum]);
    title(['Projected Dark Matter Density' char(labellist(i))])
    
    redvector = [1 0 0];
    xxx = dot(redvector, xaxis);
    yyy = dot(redvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'r', 'LineWidth', 4)

    greenvector = [0 1 0];
    xxx = dot(greenvector, xaxis);
    yyy = dot(greenvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'g', 'LineWidth', 4)

    bluevector = [0 0 1];
    xxx = dot(bluevector, xaxis);
    yyy = dot(bluevector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'b', 'LineWidth', 4)
    
    drawnow
    fname = ['f' num2str(figurenumber) '_' fn '.jpg'];     
    if testrun == 0
        saveas(figurenumber,fname)
    end 
    
    % Plot projected densities - 3.
    
    figurenumber = 100*i+33;
    figure(figurenumber)
    colormap(gray);
    maximum = 4 * sum(sum(ProjDensity))/jsize/jsize;
    imagesc(ProjDensity, [0  maximum]);
    title(['Projected Dark Matter Density' char(labellist(i))])
    
    redvector = [1 0 0];
    xxx = dot(redvector, xaxis);
    yyy = dot(redvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'r', 'LineWidth', 4)

    greenvector = [0 1 0];
    xxx = dot(greenvector, xaxis);
    yyy = dot(greenvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'g', 'LineWidth', 4)

    bluevector = [0 0 1];
    xxx = dot(bluevector, xaxis);
    yyy = dot(bluevector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'b', 'LineWidth', 4)
    
    drawnow
    fname = ['f' num2str(figurenumber) '_' fn '.jpg'];     
    if testrun == 0
        saveas(figurenumber,fname)
    end 
    
    % Plot projected densities - 4.
    
    figurenumber = 100*i+34;
    figure(figurenumber)
    colormap(gray);
    maximum = 8 * sum(sum(ProjDensity))/jsize/jsize;
    imagesc(ProjDensity, [0  maximum]);
    title(['Projected Dark Matter Density' char(labellist(i))])
    
    redvector = [1 0 0];
    xxx = dot(redvector, xaxis);
    yyy = dot(redvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'r', 'LineWidth', 4)

    greenvector = [0 1 0];
    xxx = dot(greenvector, xaxis);
    yyy = dot(greenvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'g', 'LineWidth', 4)

    bluevector = [0 0 1];
    xxx = dot(bluevector, xaxis);
    yyy = dot(bluevector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'b', 'LineWidth', 4)
    
    drawnow
    fname = ['f' num2str(figurenumber) '_' fn '.jpg'];     
    if testrun == 0
        saveas(figurenumber,fname)
    end 
    
    % Plot projected densities - 5.
    
    figurenumber = 100*i+35;
    figure(figurenumber)
    colormap(gray);
    maximum = 16 * sum(sum(ProjDensity))/jsize/jsize;
    imagesc(ProjDensity, [0  maximum]);
    title(['Projected Dark Matter Density' char(labellist(i))])
    
    redvector = [1 0 0];
    xxx = dot(redvector, xaxis);
    yyy = dot(redvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'r', 'LineWidth', 4)

    greenvector = [0 1 0];
    xxx = dot(greenvector, xaxis);
    yyy = dot(greenvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'g', 'LineWidth', 4)

    bluevector = [0 0 1];
    xxx = dot(bluevector, xaxis);
    yyy = dot(bluevector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'b', 'LineWidth', 4)
    
    drawnow
    fname = ['f' num2str(figurenumber) '_' fn '.jpg'];     
    if testrun == 0
        saveas(figurenumber,fname)
    end 
    
    % Plot projected densities times r^1.
    
    figurenumber = 100*i+30;
    figure(figurenumber)
    colormap(gray);
    maximum = max(max(ProjDensityR1));
    imagesc(ProjDensityR1, [0  maximum]);
    title(['Projected Dark Matter Density' char(labellist(i)) ' (Times r)'])

    redvector = [1 0 0];
    xxx = dot(redvector, xaxis);
    yyy = dot(redvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'r', 'LineWidth', 4)

    greenvector = [0 1 0];
    xxx = dot(greenvector, xaxis);
    yyy = dot(greenvector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'g', 'LineWidth', 4)

    bluevector = [0 0 1];
    xxx = dot(bluevector, xaxis);
    yyy = dot(bluevector, yaxis);
    fac = (jsize+1)/2;
    xend = fac + fac*xxx;
    yend = fac - fac*yyy;
    line([fac xend], [fac yend], 'Color', 'b', 'LineWidth', 4)
    
    drawnow
    fname = ['f' num2str(figurenumber) '_' fn '.jpg'];     
    if testrun == 0
        saveas(figurenumber,fname)
    end 
    
    end
    
end
































%-----------------------------Next Section:  Compute and Plot Dark Matter Potential
if 1==0

% Compute W0, W2, and WW1. 

W0(1) = 0;
W0p(1) = 0;
W0pp = U0(1) / 3;

W2(1) = 0;
W2p(1) = 0;
W2pp = U2(1) / 7;

WW1(1) = 0;
WW1p(1) = 0;
WW1pp = UU1(1) / 5;    %Changed 7 to 5

for i = 2:nsteps
    
    r = (i-1)*dr;
    
    W0(i) = W0(i-1) + W0p(i-1) * dr + W0pp * dr^2 / 2;
    W0p(i) = W0p(i-1) + W0pp * dr;
    W0pp = U0(i) - 2 * W0p(i) / r;
    
    W2(i) = W2(i-1) + W2p(i-1) * dr + W2pp * dr^2 / 2;
    W2p(i) = W2p(i-1) + W2pp * dr;
    W2pp = U2(i) - 6 * W2p(i) / r;
    
    WW1(i) = WW1(i-1) + WW1p(i-1) * dr + WW1pp * dr^2 / 2;
    WW1p(i) = WW1p(i-1) + WW1pp * dr;
    WW1pp = UU1(i) - 4 * WW1p(i) / r;                      %Changed 6 to 4
    
end

% Add the correct constants to the W so that they go to zero at infinity.

W0shift = - rmax / 1 * W0p(nsteps) - W0(nsteps);
W2shift = - rmax / 5 * W2p(nsteps) - W2(nsteps);
WW1shift = - rmax / 3 * WW1p(nsteps) - WW1(nsteps);        %Changed 5 to 3

W0 = W0 + W0shift;
W2 = W2 + W2shift;
WW1 = WW1 + WW1shift;

% Plot W functions.

if 1==1
    
    xmin = rvalues(1);
    xmax = rvalues(graphsteps);
    ymin = min(min(min(W0(1:graphsteps)),min(-rvalues(1:graphsteps) .^ 2 .* W2(1:graphsteps))),min(WW1(1:graphsteps) .* rvalues(1:graphsteps) .^ 1));
    ymax = max(max(max(W0(1:graphsteps)),max(-rvalues(1:graphsteps) .^ 2 .* W2(1:graphsteps))),max(WW1(1:graphsteps) .* rvalues(1:graphsteps) .^ 1));

    figure(5)
    plot(rvalues(1:graphsteps), W0(1:graphsteps))
    axis([xmin,xmax,ymin,ymax])
    title('W0')
    fname = ['f5_' fn '.jpg'];     
    % saveas(5,fname)

    figure(6)
    plot(rvalues(1:graphsteps), -rvalues(1:graphsteps) .^ 2 .* W2(1:graphsteps))
    axis([xmin,xmax,ymin,ymax])
    title('-r^2 * W2')
    fname = ['f6_' fn '.jpg'];     
    % saveas(6,fname)

% figure(7) had been for W4

    figure(7)
    plot(rvalues(1:graphsteps), WW1(1:graphsteps) .* rvalues(1:graphsteps) .^ 1)    %Changed ^2 to ^1
    axis([xmin,xmax,ymin,ymax])
    title('r^1 * WW1')             %Changed ^2 to ^1
    fname = ['f7_' fn '.jpg'];     
    % saveas(7,fname)
end
    
% Create 3D graphs of cross sections of the potential function V.

XYPotential = zeros(jsize);
XZPotential = zeros(jsize);
YZPotential = zeros(jsize);

for jrow = 1:jsize
    for jcol = 1:jsize
        
        v1 = (jcol-1)/(jsize-1) * 2 * graphradius - graphradius;
        v2 = (jrow-1)/(jsize-1) * 2 * graphradius - graphradius;
        
        XYPotential(jrow,jcol) = Vfunct(v1,v2,0,0); 
        XZPotential(jrow,jcol) = Vfunct(v1,0,v2,0);
        YZPotential(jrow,jcol) = Vfunct(0,v1,v2,0);

    end
end

figure(9)
contour(XYPotential);
colormap('lines')
shading interp
title('The Potential Function in the xy Plane')
xlabel('x axis')
ylabel('y axis')
zlabel('Potential')
fname = ['f9_' fn '.jpg'];     
% saveas(9,fname)

figure(10)
contour(XZPotential);
colormap('lines')
shading interp
title('The Potential Function in the xz Plane')
xlabel('x axis')
ylabel('z axis')
zlabel('Potential')
fname = ['f10_' fn '.jpg'];     
% saveas(10,fname)

figure(11)
contour(YZPotential);
colormap('lines')
shading interp
title('The Potential Function in the yz Plane')
xlabel('y axis')
ylabel('z axis')
zlabel('Potential')
fname = ['f11_' fn '.jpg'];     
% saveas(11,fname)
    

XaxisPotential = zeros(1,jjsize);
Xaxis = zeros(1,jjsize);

for j = 1:jjsize
        
    v1 = (j-1)/(jjsize-1) * 2 * graphradius - graphradius;
    Xaxis(j) = v1;            
    XaxisPotential(j) = Vfunct(v1,0,0,0); 
   
end

figure(8)
plot(Xaxis,XaxisPotential)
title('The Potential along the x Axis')

end












% Functions

    function ruv = randomunitvector()
        
        zrand = 2*rand - 1;
        rrand = sqrt(1 - zrand^2);
        
        thetarand = rand*2*pi;
        xrand = rrand*cos(thetarand);
        yrand = rrand*sin(thetarand);

        ruv = [xrand yrand zrand];

    end

    function output = DMdensity(x,y,z,t)
        
        r = sqrt(x^2 + y^2 + z^2);
        
        U0value = 0;
        U2value = 0;
        UU1value = 0;

        if r <= rmax
            ival = 1 + round(r/dr);
            U0value = U0(ival);
            U2value = U2(ival);
            UU1value = UU1(ival);
        end
       
        temp = U0value + U2value * (3*z^2 - r^2);                                   %updated
        temp2 = temp + UU1value * (cos(dw*t) * (x) + sin(dw*t) * (y));              %updated
        output = mu0 * temp2;
        
    end

    function output = Vfunct(x,y,z,t)
        
        r = sqrt(x^2 + y^2 + z^2);
        
        if r <= rmax
            ival = 1 + round(r/dr);
            W0value = W0(ival);
            W2value = W2(ival);
            WW1value = WW1(ival);
        else
            W0value = W0(nsteps) * (rmax/r) ^ 1;
            W2value = W2(nsteps) * (rmax/r) ^ 5;
            WW1value = WW1(nsteps) * (rmax/r) ^ 3;    %Changed 5 to 3
        end
        
        temp = W0value + W2value*(3*z^2 - r^2);                                      %updated
        temp2 = temp + WW1value * (cos(dw*t) * (x) + sin(dw*t) * (y));
        output = 4*pi*mu0 * temp2;
        
    end

    function output = GradVfunct(x,y,z,t)
        
        r = sqrt(x^2 + y^2 + z^2);
        if r == 0
            output = [0 0 0];
        else
        
            if r <= rmax
                ival = 1 + round(r/dr);
            
                W0value = W0(ival);
                W2value = W2(ival);
                WW1value = WW1(ival);
            
                W0pvalue = W0p(ival);
                W2pvalue = W2p(ival);
                WW1pvalue = WW1p(ival);
            
            else
                W0value = W0(nsteps) * (rmax/r) ^ 1;
                W2value = W2(nsteps) * (rmax/r) ^ 5;
                WW1value = WW1(nsteps) * (rmax/r) ^ 3;     %Changed 5 to 3
            
                W0pvalue = -W0value / r * 1;
                W2pvalue = -W2value / r * 5;
                WW1pvalue = -WW1value / r * 3;    %Changed 5 to 3   
            
            end
        
            term1 = W0pvalue + W2pvalue * (3*z^2 - r^2);                 %updated
            term2 = WW1pvalue * (cos(dw*t) * (x) + sin(dw*t) * (y));     %updated
            
            vector1 = (term1 + term2) * [x/r, y/r, z/r];
        
            vector2 = W2value * [-2*x, -2*y, 4*z];
            
            vector3 = WW1value * [cos(dw*t), sin(dw*t), 0];              %updated
            
            output = 4*pi*mu0*(vector1 + vector2 + vector3);             %updated
            
        end
        
    end 

    function [goodfilename] = nodecimals(badfilename)
        
        stringsize = size(badfilename);
        stringlength = stringsize(2);
        
        goodfilename = [];
        for k = 1:stringlength
            badchar = badfilename(k);
            if badchar == '.'
                goodchar = 'p';
            else
                goodchar = badchar;
            end
            goodfilename = [goodfilename goodchar];
        end
           
    end


end
