I need help getting this code to work in the current version of matlab

2 views (last 30 days)
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
  6 Comments
Eric Pahlke
Eric Pahlke on 17 Jul 2015
Also, I highly recommend specifically pointing out the problem you're having (not running, error on line x, produces unexpected output). It helps people answering make sure they're seeing the same results you are.
Walter Roberson
Walter Roberson on 17 Jul 2015
I ran it for myself, and it seemed to go fine up until the point I interrupted it. The only "problem" I saw was that for any one series of cross-section plots, the plots were essentially indistinguishable except for the label (and sometimes a very small difference in the amount of white available where lines intersected.) But I have no reason to expect that wasn't the intended output.

Sign in to comment.

Answers (1)

Steve Miller
Steve Miller on 21 Jun 2021
Comments above indicate the code is working. If you can explain why this code is not performing as as expected, add some details and we can look at it again.
--Steve

Categories

Find more on Discrete Data Plots in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!