function [] = ellipticalgalaxy(dr, rmax, A0, A2, wavelength0, wavelength2, patternperiod, mu0, initialradius, nparticles, TotalTime, dt, binsize)

n1 = 5000;               % This is the number of loops in the simulation between updates on the computer screen.

jsize = 300;              % The plots of the dark matter density and its potential are jsize x jsize.
jjsize = 1000;            % The images of the galaxies, both on the computer screen and what is written as jpegs to disk, are jjsize x jjsize.
dotsizes = [3 5 10];   % These are radii of the translucent dots used to create the images on the computer screen.  The second one is used to created the jpegs to disk.                  

nmedians4graphs = 3;    % The radius of the graphs created will be approximately this many standard deviations of the particles.

% ellipticalgalaxy(.1,100000,1,-.5,2000,1990,50000000,8.7e-13,15000,1500,1e20,100000,1000);

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

% wavelength0 and wavelength2 are the wavelengths of the zeroth and second
% spherical harmonic functions f0 and f2.
% 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.

% Initially there are nparticles evenly distributed on a sphere of radius
% initialradius.  The simulation will run for TotalTime with stepsizes dt.

% binsize is used to compute the frequency of stars at each radius
% (projected into the xy plane).  Each radius (in the xy plane) is then
% rounded to the nearest integer multiple of binsize and then put into the
% corresponding bin.  The resulting xy-radial frequency distribution can then be
% compared to observations of brightness as a function of radius in actual
% elliptical galaxies (see Galactic Astronomy by Binney and Merrifield, p. 186.)

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

% Computer k0, k2, w0, w2, Omega.

k0 = 2*pi / wavelength0;
k2 = 2*pi / wavelength2;

w0 = patternperiod/8/pi * (k2^2 - k0^2) - 2*pi/patternperiod
w2 = patternperiod/8/pi * (k2^2 - k0^2) + 2*pi/patternperiod
dw = w2-w0;

Omega0Squared = (2*pi/patternperiod)^2 + (patternperiod/8/pi)^2 * (k2^2 - k0^2)^2 - (k2^2 + k0^2)/2;

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

Omega0 = sqrt(Omega0Squared)
OneOverOmega0 = 1/Omega0
        
% Initialize f0, f2, U0, U2, U4, UU2, W0, W2, W4, and WW2 (which stands for W tilde sub 2).

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

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

U0 = zeros(1, nsteps);
U2 = zeros(1, nsteps);
U4 = zeros(1, nsteps);
UU2 = zeros(1, nsteps);

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

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

W4 = zeros(1, nsteps);
W4p = zeros(1, nsteps);

WW2 = zeros(1, nsteps);
WW2p = zeros(1, nsteps);

% Compute f0 and f2.

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

f2(1) = 1;
f2p = 0;
f2pp = -k2^2 * f2(1) / 7;

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;
    
    f2(i) = f2(i-1) + f2p * dr + f2pp * dr^2 / 2;
    f2p = f2p + f2pp * dr;
    f2pp = -k2^2 * f2(i) - 6 * f2p / r;
    
end

% Normalize f.

f0 = -f0 / min(f0 .* rvalues);
f2 = -f2 / min(f2 .* rvalues .^ 3);

fac = max(f0);

f0 = f0 / fac;
f2 = f2 / fac;

% Plot f.

fmin = min(min(f0),min(f2));
fmax = max(max(f0),max(f2));

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

fn = [num2str(dr) '_' num2str(rmax) '_' num2str(A0) '_' num2str(A2) '_' num2str(wavelength0) '_' num2str(wavelength2)];
fn = [fn '_' num2str(patternperiod) '_' num2str(mu0) '_' num2str(initialradius) '_' num2str(nparticles) '_' num2str(0) '_' num2str(dt) '_' num2str(binsize)];
fn = [fn 'P' num2str(OneOverOmega0) '_' num2str(graphradius) 'P'];
fn = nodecimals(fn);

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

% Compute U0, U2, U4, UU2.

for i = 1:nsteps
    
    r = (i-1)*dr;
    
    U0(i) = A0^2 * f0(i)^2 + 56/105 * A2^2 * r^4 * f2(i)^2;
    U2(i) = -40/105 * A2^2 * r^2 * f2(i)^2;
    U4(i) = 3/105 * A2^2 * f2(i)^2;
    UU2(i) = 2*A0*A2 * f0(i) * f2(i);
       
end

% Compute and plot dark matter density in xy plane.

XYDensity = zeros(jsize);
XZDensity = zeros(jsize);
YZDensity = zeros(jsize);

for jrow = 1:jsize
    for jcol = 1:jsize
        
        v1 = (jcol-1)/(jsize-1) * 2 * graphradius - graphradius;
        v2 = graphradius - 2 * graphradius * (jrow-1)/(jsize-1);
        
        XYDensity(jrow,jcol) = DMdensity(v1,v2,0,0) * (v1^2 + v2^2);
        XZDensity(jrow,jcol) = DMdensity(v1,0,v2,0) * (v1^2 + v2^2);
        YZDensity(jrow,jcol) = DMdensity(0,v1,v2,0) * (v1^2 + v2^2);
        
    end
end

figure(2)
colormap(gray);
maximum = max(YZDensity(:));
imagesc(YZDensity, [0  maximum]);
title('Dark Matter Density (Times r^2) in the yz Plane at t=0')
fname = ['f2_' fn '.jpg']     
saveas(2,fname)

figure(3)
colormap(gray);
maximum = max(XZDensity(:));
imagesc(XZDensity, [0  maximum]);
title('Dark Matter Density (Times r^2) in the xz Plane at t=0')
fname = ['f3_' fn '.jpg']     
saveas(3,fname)

figure(4)
colormap(gray);
maximum = max(XYDensity(:));
imagesc(XYDensity, [0  maximum]);
title('Dark Matter Density (Times r^2) in the xy Plane at t=0')
fname = ['f4_' fn '.jpg']     
saveas(4,fname)

% Compute W0, W2, W4, and WW2. 

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

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

W4(1) = 0;
W4p(1) = 0;
W4pp = U4(1) / 11;

WW2(1) = 0;
WW2p(1) = 0;
WW2pp = UU2(1) / 7;

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;
    
    W4(i) = W4(i-1) + W4p(i-1) * dr + W4pp * dr^2 / 2;
    W4p(i) = W4p(i-1) + W4pp * dr;
    W4pp = U4(i) - 10 * W4p(i) / r;
    
    WW2(i) = WW2(i-1) + WW2p(i-1) * dr + WW2pp * dr^2 / 2;
    WW2p(i) = WW2p(i-1) + WW2pp * dr;
    WW2pp = UU2(i) - 6 * WW2p(i) / r;
    
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);
W4shift = - rmax / 9 * W4p(nsteps) - W4(nsteps);
WW2shift = - rmax / 5 * WW2p(nsteps) - WW2(nsteps);

W0 = W0 + W0shift;
W2 = W2 + W2shift;
W4 = W4 + W4shift;
WW2 = WW2 + WW2shift;

% Plot W functions.

if 1==1
    figure(5)
    plot(rvalues(1:graphsteps), W0(1:graphsteps))
    title('W0')
    fname = ['f5_' fn '.jpg']     
    saveas(5,fname)

    figure(6)
    plot(rvalues(1:graphsteps), -rvalues(1:graphsteps) .^ 2 .* W2(1:graphsteps))
    title('-r^2 * W2')
    fname = ['f6_' fn '.jpg']     
    saveas(6,fname)

    figure(7)
    plot(rvalues(1:graphsteps), 3 * W4(1:graphsteps) .* rvalues(1:graphsteps) .^ 4)
    title('3r^4 * W4')
    fname = ['f7_' fn '.jpg']     
    saveas(7,fname)

    figure(8)
    plot(rvalues(1:graphsteps), WW2(1:graphsteps) .* rvalues(1:graphsteps) .^ 2)
    title('r^2 * WW2')
    fname = ['f8_' fn '.jpg']     
    saveas(8,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)
    
%----------------Begin Simulation of Galaxy-------------------------------------------

m = ones(1,nparticles);
x = zeros(1,nparticles);
y = zeros(1,nparticles);
z = zeros(1,nparticles);
vx = zeros(1,nparticles);
vy = zeros(1,nparticles);
vz = zeros(1,nparticles);

% Determine initial positions and velocities of particles.

randomphi = [0];
kount = 1;
kround = 0;

while kount < nparticles

   kround = kround + 1;
   for i = 1:kount
       randomphi(kount + i) = randomphi(i) + 1/2^kround;
   end
   kount = kount * 2;
    
end

randomphi = randomphi * (2*pi);

for i = 1:nparticles,
    
    phi = randomphi(i);
    
    theta = acos( i/(nparticles+1) *2-1);
    
    x(i) = initialradius*sin(theta)*cos(phi);
    y(i) = initialradius*sin(theta)*sin(phi);
    z(i) = initialradius*cos(theta);
   
    vx(i) = 0;
    vy(i) = 0;
    vz(i) = 0;
    
end

ax = zeros(1,nparticles);
ay = zeros(1,nparticles);
az = zeros(1,nparticles);

if 1==0
    figure(12)
    quiver3(x,y,z,vx,vy,vz)
    axis([-graphradius graphradius -graphradius graphradius -graphradius graphradius]);
    title('Initial positions and velocities')
    xlabel('x axis')
    ylabel('y axis')
    zlabel('z axis')
    fname = ['f12_' fn '.jpg']     
    saveas(12,fname)
    drawnow 
end

% Begin Elliptical Galaxy Simulation ---------------------------------------

Maxloops = round(TotalTime/dt);

for nloops = 0:Maxloops
    
    t = nloops*dt;
    
    % Setup graphics names for repeated plots.

    fn = [num2str(dr) '_' num2str(rmax) '_' num2str(A0) '_' num2str(A2) '_' num2str(wavelength0) '_' num2str(wavelength2)];
    fn = [fn '_' num2str(patternperiod) '_' num2str(mu0) '_' num2str(initialradius) '_' num2str(nparticles) '_' num2str(t) '_' num2str(dt) '_' num2str(binsize)];
    fn = [fn 'P' num2str(OneOverOmega0) '_' num2str(graphradius) 'P'];
    fn = nodecimals(fn);
    
    % Define Kinetic Energy, Particle-Particle Potential Energy, 
    % Dark Matter-Particle Potential Energy, and Angular Momentum.
    
    KE = 0;
    DPE = 0;
    ANG = [0 0 0];
     
    % Note that KE + DPE and ANG would be conserved quantities if the dark
    % matter potential were static and spherically symmetric.  Since this
    % is not the case, their lack of conservation represents transfers of
    % energy and angular momentum to the test particles.
    
    clear radialcount;
    nbins = 1;
    radialcount(1:nbins) = 0;

    for j = 1:nparticles
        
        xyr = sqrt(x(j)^2 + y(j)^2);
        ibin = round(xyr/binsize) + 1;
        if ibin > nbins
            radialcount(nbins+1:ibin) = 0;
            nbins = ibin;
        end
        radialcount(ibin) = radialcount(ibin) + 1;
        
        KE = KE + (1/2)*m(j)*(vx(j)^2 + vy(j)^2 + vz(j)^2);        
        DPE = DPE + m(j)*Vfunct(x(j),y(j),z(j),t);
        ANG = ANG + m(j)*cross([x(j) y(j) z(j)], [vx(j) vy(j) vz(j)]);
        
        accel = -GradVfunct(x(j),y(j),z(j),t);
        ax(j) = accel(1);
        ay(j) = accel(2);
        az(j) = accel(3);
   
    end
    
    if nloops == 0;
        ANG0 = ANG;
        KE0 = KE;
        DPE0 = DPE;
        TE0 = KE+DPE;
    end 
    
    x = x + vx*dt + ax*dt*dt/2;
    y = y + vy*dt + ay*dt*dt/2;
    z = z + vz*dt + az*dt*dt/2;
    
    vx = vx + ax*dt;
    vy = vy + ay*dt;   
    vz = vz + az*dt;
    
    if nloops == round(nloops/n1)*n1
    
        nbins
        
        % Determine median radial value (to the nearest bin).
        
        ihalf = 0;
        halfcount = 0;
        while halfcount < nparticles/2
            ihalf = ihalf + 1; 
            halfcount = halfcount + radialcount(ihalf);
        end
        medianr = (ihalf - 1)*binsize
        if medianr == 0
            medianr = binsize/2
        end
        
        % Compute r^(1/4) radial distribution.
        % See Galactic Astronomy by Binney and Merrifield, p. 186.
        
        rrr = 0:binsize:(nbins-1)*binsize;
        Coef = nparticles / 7.22 / pi / medianr^2;
        radialmodel = binsize * Coef * 2*pi*(rrr .* exp( -7.67*((rrr/medianr).^(1/4) - 1)));
        
        maxmodel = 2*pi*Coef*(3);
        maxdata = max(radialcount);
        ymax = max(maxmodel, maxdata*1.05);
        xmax = medianr*10;
        
        % Plot xy radial distribution alongside r^(1/4) model.
        
        figure(13) 
        plot(rrr, radialcount, 'b')
        hold on
        plot(rrr, radialmodel, 'r')
        hold off
        title('Radial Frequency (in Blue) in the xy Plane Compared to 1/4 Model (in Red)')
        xlabel(['binsize = ' sprintf('%g', binsize)])
        axis([0 xmax 0 ymax])
        drawnow
        fname = ['f13_' fn '.jpg']     
        saveas(13,fname)
        
        % Plot star locations.
        
        figure(15)
        dsize = 1;
        %quiver3(x,y,z,vx,vy,vz)
        scatter3(x,y,z,dsize,'filled')
        angle = dw*t/2;
        xbar = graphradius*cos(angle);
        ybar = graphradius*sin(angle);
        hold on
        quiver3(0,0,0,xbar, ybar, 0)
        hold off
        axis([-graphradius graphradius -graphradius graphradius -graphradius graphradius]);
        title(sprintf('%g %g %g %g %g %g %g', [dr rmax A0 A2 wavelength0 wavelength2 patternperiod]))
        xlabel(sprintf('%g %g %g %g %g %g', [mu0 initialradius nparticles t dt binsize]))
        ylabel(sprintf('%g %g %g',[OneOverOmega0 graphradius dsize]))
        view([0,0,1])
        set(gca,'CameraUpVector',[0 1 0])       
        drawnow
        fname = ['f15_' fn '.jpg']     
        saveas(15,fname)
        
        for k = 1:nfigs
            
            [dd,maxdd] = xyellipticaldensity(jjsize, dotsizes(k)); 
            ddxy = dd/maxdd;
            figure(20 + k)
            colormap(gray);
            imagesc(ddxy, [0  1]);
            title(sprintf('%g %g %g %g %g %g %g', [dr rmax A0 A2 wavelength0 wavelength2 patternperiod]))
            xlabel(sprintf('%g %g %g %g %g %g', [mu0 initialradius nparticles t dt binsize]))
            ylabel(sprintf('%g %g %g',[OneOverOmega0 graphradius dotsizes(k)]))
            drawnow
            %fname = ['f' num2str(20+k) '_' fn '.jpg']     
            %saveas(20+k,fname)
        
            [dd,maxdd] = xzellipticaldensity(jjsize, dotsizes(k)); 
            ddxz = dd/maxdd;
            figure(30 + k)
            colormap(gray);
            imagesc(ddxz, [0  1]);
            title(sprintf('%g %g %g %g %g %g %g', [dr rmax A0 A2 wavelength0 wavelength2 patternperiod]))
            xlabel(sprintf('%g %g %g %g %g %g', [mu0 initialradius nparticles t dt binsize]))
            ylabel(sprintf('%g %g %g',[OneOverOmega0 graphradius dotsizes(k)]))
            drawnow
            %fname = ['f' num2str(30+k) '_' fn '.jpg']     
            %saveas(30+k,fname)
                          
            fn = [num2str(dr) '_' num2str(rmax) '_' num2str(A0) '_' num2str(A2) '_' num2str(wavelength0) '_' num2str(wavelength2)];
            fn = [fn '_' num2str(patternperiod) '_' num2str(mu0) '_' num2str(initialradius) '_' num2str(nparticles) '_' num2str(t) '_' num2str(dt) '_' num2str(binsize)];
            fn = [fn 'P' num2str(OneOverOmega0) '_' num2str(graphradius) '_' num2str(dotsizes(k)) 'P'];
            fn = nodecimals(fn);
            
            fname = ['xy_' num2str(k) '_' fn '.jpg']
            imwrite(ddxy, fname, 'jpg')
                 
            fname = ['xz_' num2str(k) '_' fn '.jpg']
            imwrite(ddxz, fname, 'jpg')
           
        end
        
        % Compute ellipticities and new graphradius, if necessary.
                 
        EllipticityMatrix = [x*x' x*y' x*z'; y*x' y*y' y*z'; z*x' z*y' z*z']/nparticles;
        a2b2c2 = eig(EllipticityMatrix);
        abc = sqrt(a2b2c2);
        abcnormalized = abc/abc(2);
                 
        %stddevradius = max(abc)*sqrt(3);
        graphradius = max([graphradius, medianr*nmedians4graphs]);
                 
        TE = KE + DPE;
        
        nloops
        ElapseTime = nloops*dt
        
        abc
        abcnormalized
        
        ANG0
        ANG
        TE0
        TE
        
        
    end
    
end

Omega0
OneOverOmega0


        
       


%-----------------------END OF PROGRAM-------------------------------






% Functions

    function [density, maxdensity] = xyellipticaldensity(jjsize, dotradius)
        
       
        density = zeros(jjsize);
        maxdensity = 0;

        for i = 1:nparticles
    
            jcol0 = round((x(i)+graphradius)/2/graphradius*(jjsize-1) + 1);
            jrow0 = (jjsize+1) - round((y(i)+graphradius)/2/graphradius*(jjsize-1) + 1);

            for djcol = -dotradius:dotradius
            for djrow = -dotradius:dotradius
        
                if djcol^2 + djrow^2 <= dotradius^2
            
                    jcol = jcol0 + djcol;
                    jrow = jrow0 + djrow;
                    if ((jcol > 0) & (jcol <= jjsize) & (jrow > 0) & (jrow <= jjsize))
                
                    density(jrow,jcol) = density(jrow,jcol)+1;
                    maxdensity = max(maxdensity,density(jrow,jcol));
                
                    end
            
                end
        
            end
            end
    
        end

    end


    function [density, maxdensity] = xzellipticaldensity(jjsize, dotradius)
        
       
        density = zeros(jjsize);
        maxdensity = 0;

        for i = 1:nparticles
    
            jcol0 = round((x(i)+graphradius)/2/graphradius*(jjsize-1) + 1);
            jrow0 = (jjsize+1) - round((z(i)+graphradius)/2/graphradius*(jjsize-1) + 1);

            for djcol = -dotradius:dotradius
            for djrow = -dotradius:dotradius
        
                if djcol^2 + djrow^2 <= dotradius^2
            
                    jcol = jcol0 + djcol;
                    jrow = jrow0 + djrow;
                    if ((jcol > 0) & (jcol <= jjsize) & (jrow > 0) & (jrow <= jjsize))
                
                    density(jrow,jcol) = density(jrow,jcol)+1;
                    maxdensity = max(maxdensity,density(jrow,jcol));
                
                    end
            
                end
        
            end
            end
    
        end

    end


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

        if r <= rmax
            ival = 1 + round(r/dr);
            U0value = U0(ival);
            U2value = U2(ival);
            U4value = U4(ival);
            UU2value = UU2(ival);
        end
       
        temp = U0value + U2value * (3*z^2 - r^2) + U4value * (35*z^4 - 30*r^2*z^2 + 3*r^4);
        temp2 = temp + UU2value * (cos(dw*t) * (x^2 - y^2) + sin(dw*t) * (2*x*y));
        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);
            W4value = W4(ival);
            WW2value = WW2(ival);
        else
            W0value = W0(nsteps) * (rmax/r) ^ 1;
            W2value = W2(nsteps) * (rmax/r) ^ 5;
            W4value = W4(nsteps) * (rmax/r) ^ 9;
            WW2value = WW2(nsteps) * (rmax/r) ^ 5;
        end
        
        temp = W0value + W2value*(3*z^2 - r^2) + W4value*(35*z^4 - 30*r^2*z^2 + 3*r^4);
        temp2 = temp + WW2value * (cos(dw*t) * (x^2 - y^2) + sin(dw*t) * (2*x*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);
                W4value = W4(ival);
                WW2value = WW2(ival);
            
                W0pvalue = W0p(ival);
                W2pvalue = W2p(ival);
                W4pvalue = W4p(ival);
                WW2pvalue = WW2p(ival);
            
            else
                W0value = W0(nsteps) * (rmax/r) ^ 1;
                W2value = W2(nsteps) * (rmax/r) ^ 5;
                W4value = W4(nsteps) * (rmax/r) ^ 9;
                WW2value = WW2(nsteps) * (rmax/r) ^ 5;
            
                W0pvalue = -W0value / r * 1;
                W2pvalue = -W2value / r * 5;
                W4pvalue = -W4value / r * 9;
                WW2pvalue = -WW2value / r * 5;   
            
            end
        
            term1 = W0pvalue + W2pvalue * (3*z^2 - r^2) + W4pvalue * (35*z^4 - 30*r^2*z^2 + 3*r^4);
            term2 = WW2pvalue * (cos(dw*t) * (x^2 - y^2) + sin(dw*t) * (2*x*y));
            vector1 = (term1 + term2) * [x/r, y/r, z/r];
        
            vector2 = W2value * [-2*x, -2*y, 4*z];
        
            vector3 = W4value * [12*x*(r^2 - 5*z^2), 12*y*(r^2 - 5*z^2), 16*z*(5*z^2 - 3*r^2)];
            
            vector4 = cos(dw*t) * WW2value * [2*x, -2*y, 0];
            
            vector5 = sin(dw*t) * WW2value * [2*y, 2*x, 0];
        
            output = 4*pi*mu0*(vector1 + vector2 + vector3 + vector4 + vector5);
            
        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
