Skip to content

Aenglish #7

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
humdrumcomet wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Aenglish #7

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
96388ef
Added files via upload
Apr 20, 2016
6b186fc
Added files via upload
Apr 20, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
84 changes: 51 additions & 33 deletions MDCode/InitVStream.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,33 +1,51 @@
doPlot = 1;
dt = 5e-15;
TStop = 3000 * dt;
InitDist = 0.0;
Method = 'VE'; % VE -- verlot; FD -- Forward Difference

Mass0 = 14 * C.am; % Silicon
Mass1 = 100 * C.am; % Argon

AtomSpacing = 0.5430710e-9;
LJSigma = AtomSpacing / 2^(1 / 6);
LJEpsilon = 1e-21;

PhiCutoff = 3 * AtomSpacing * 1.1;

T = 30;

AddRectAtomicArray(10, 10, 0, 0, 0, 0, 0, T, 0);
% vy0 = -sqrt(0.02*Ep/Mass1);
% AddRectAtomicArray(4,4,0,12*AtomSpacing,0,vy0,0,T,1);
Ep = 2;
AddParticleStream(5, 0.1, 10, -pi / 2, 1, Ep * C.q_0, 5);

Size = 10*AtomSpacing;
Limits = [-Size +Size -Size +Size]; % square is good
PlDelt = 5 * dt;

PlotFile = 'BlockSt.gif';
doPlotImage = 1;
PlotSize = [100, 100, 1049, 1049];

ScaleV = .02e-11;
ScaleF = 10;
doPlot = 1;
dt = 5e-15;
TStop = 3000 * dt;
InitDist = 0.0;
Method = 'VE'; % VE -- verlot; FD -- Forward Difference

Mass0 = 14 * C.am; % Silicon
Mass1 = 100 * C.am; % Argon

AtomSpacing = 0.5430710e-9;
LJSigma = AtomSpacing / 2^(1 / 6);
LJEpsilon = 1e-21;

PhiCutoff = 3 * AtomSpacing * 1.1;

T = 30;

LAtoms=10;
WAtoms=10;
X0=0;
Y0=0;
VX0=0;
VY0=0;
InitDist=0;
Temp=T;
Type=0;
AddRectAtomicArray(LAtoms, WAtoms, X0, Y0, VX0, VY0, InitDist, Temp, Type);

% vy0 = -sqrt(0.02*Ep/Mass1);
% AddRectAtomicArray(4,4,0,12*AtomSpacing,0,vy0,0,T,1);

Ep = 2;
numStream=5;
x0=0.1;
y0=10;
PartAng=-pi/2;
Type=1;
Ep=Ep*C.q_0;
Seper=5;
AddParticleStream(numStream, x0, y0, PartAng, Type, Ep, Seper);

Size = 10*AtomSpacing;
Limits = [-Size +Size -Size +Size]; % square is good
PlDelt = 5 * dt;

PlotFile = 'BlockSt.gif';
doPlotImage = 1;
PlotSize = [100, 100, 1049, 1049];

ScaleV = .02e-11;
ScaleF = 10;
216 changes: 111 additions & 105 deletions MDCode/PlotVars.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,105 +1,111 @@
function [ output_args ] = PlotVars(c, Limits)
global Vx Vy L W x y Fx Fy C
global Phi nAtoms time Mass0 Mass1 Pty0in Pty1in
global LJEpsilon Phi0 PhiTot KETot MinX MaxX MinY MaxY
global T T0 T1 ScaleV MarkerSize doPlotImage PlotCount
global PlotFig map im PlotSize ScaleF

if isempty(Limits)
Limits = Limits;
end

if doPlotImage == 1
if PlotCount == 0
close all;
set(0,'DefaultFigureWindowStyle','normal')
PlotFig = figure('Position', PlotSize);
PlotCount = 1;
else
PlotCount = PlotCount+1;
end
end

V2 = Vx.*Vx + Vy.*Vy;
MaxV = max(sqrt(V2));
dx = MaxX - MinX;
dy = MaxY - MinY;

if c==0 && ~ScaleV
ScaleV = sqrt(dx*dx + dy*dy) / MaxV * 0.15;
end

subplot(3, 2, 1), plot(x(Pty0in), y(Pty0in), 'bo', 'markers',...
MarkerSize,'MarkerFaceColor', 'b');
hold on
subplot(3, 2, 1), plot(x(Pty1in), y(Pty1in), 'go', 'markers',...
MarkerSize,'MarkerFaceColor', 'g');
subplot(3, 2, 1),quiver(x, y, Fx * ScaleF, Fy * ScaleF, 0, 'm', 'linewidth', 2);
hold off
axis(Limits);
title('Atoms and Forces')
xlabel('X')
ylabel('Y')

subplot(3, 2, 2), plot(time, KETot, 'linewidth', 2);
hold on
subplot(3, 2, 2), plot(time, PhiTot + LJEpsilon * nAtoms, 'g', 'linewidth',2);
subplot(3, 2, 2), plot(time, PhiTot + KETot + LJEpsilon * nAtoms, 'k', 'linewidth',2);
% axis([0 TStop 0e-22 3e-23]);
xlabel('time')
ylabel('Energy')
title('Total, Kinetic and Potential Energy')

subplot(3,2,3), plot(x(Pty0in), y(Pty0in), 'bo', 'markers',...
MarkerSize,'MarkerFaceColor','b');
hold on
subplot(3,2,3), plot(x(Pty1in), y(Pty1in), 'go', 'markers',...
MarkerSize,'MarkerFaceColor','g');
subplot(3,2,3),quiver(x, y, Vx * ScaleV, Vy * ScaleV, 0, 'r', 'linewidth', 2);
hold off
axis(Limits);
title('Atoms and Velocities')
xlabel('X')
ylabel('Y')

AE(Pty0in) = 1 / 2 * Mass0 * V2(Pty0in) + Phi(Pty0in) - Phi0;
AE(Pty1in) = 1 / 2 * Mass1 * V2(Pty1in) + Phi(Pty1in) - Phi0;
% if AddParticle, AE(end) = 0;end
% AE = Phi;
subplot(3, 2, 4),scatter3(x,y,AE,ones(1,nAtoms)*300,AE,'fill');
axis(Limits);
view(2);
title('Atom Energy (KE + Pot)')

subplot(3, 2, 5), plot(x(Pty0in), y(Pty0in), 'bo', 'markers',...
MarkerSize, 'MarkerFaceColor', 'b');
hold on
subplot(3, 2, 5), plot(x(Pty1in), y(Pty1in), 'go', 'markers',...
MarkerSize, 'MarkerFaceColor', 'g');
hold off
title('Atoms (All)')
xlabel('X')
ylabel('Y')

subplot(3, 2, 6), plot(time, T, 'k', 'linewidth', 2);
hold on
subplot(3, 2, 6), plot(time, T0, 'b', 'linewidth', 2);
subplot(3, 2, 6), plot(time, T1, 'g', 'linewidth', 2);
%
xlabel('time')
ylabel('Temperature')
title('Temperature: Each type and all')

if doPlotImage
if PlotCount == 1
f = getframe(PlotFig);
[im, map] = rgb2ind(f.cdata, 256, 'nodither');
im(1, 1, 1, 2) = 0;

else
f = getframe(PlotFig);
im(:, :, 1, PlotCount) = rgb2ind(f.cdata, map, 'nodither');
end
end

end
function [ output_args ] = PlotVars(c, Limits)
global Vx Vy L W x y Fx Fy C
global Phi nAtoms time Mass0 Mass1 Pty0in Pty1in
global LJEpsilon Phi0 PhiTot KETot MinX MaxX MinY MaxY
global T T0 T1 ScaleV MarkerSize doPlotImage PlotCount
global PlotFig map im PlotSize ScaleF

if isempty(Limits)
Limits = Limits;
end

if doPlotImage == 1
if PlotCount == 0
close all;
set(0,'DefaultFigureWindowStyle','normal')
PlotFig = figure('Position', PlotSize);
PlotCount = 1;
else
PlotCount = PlotCount+1;
end
end

V2 = Vx.*Vx + Vy.*Vy;
MaxV = max(sqrt(V2));
dx = MaxX - MinX;
dy = MaxY - MinY;

if c==0 && ~ScaleV
ScaleV = sqrt(dx*dx + dy*dy) / MaxV * 0.15;
end

subplot(3, 2, 1), plot(x(Pty0in), y(Pty0in), 'bo', 'markers',...
MarkerSize,'MarkerFaceColor', 'b');
hold on
subplot(3, 2, 1), plot(x(Pty1in), y(Pty1in), 'go', 'markers',...
MarkerSize,'MarkerFaceColor', 'g');
subplot(3, 2, 1),quiver(x, y, Fx * ScaleF, Fy * ScaleF, 0, 'm', 'linewidth', 2);
hold off
axis(Limits);
title('Atoms and Forces')
xlabel('X')
ylabel('Y')

subplot(3, 2, 2), plot(time, KETot, 'linewidth', 2);
hold on
subplot(3, 2, 2), plot(time, PhiTot + LJEpsilon * nAtoms, 'g', 'linewidth',2);
subplot(3, 2, 2), plot(time, PhiTot + KETot + LJEpsilon * nAtoms, 'k', 'linewidth',2);
% axis([0 TStop 0e-22 3e-23]);
xlabel('time')
ylabel('Energy')
title('Total, Kinetic and Potential Energy')

subplot(3,2,3), plot(x(Pty0in), y(Pty0in), 'bo', 'markers',...
MarkerSize,'MarkerFaceColor','b');
hold on
subplot(3,2,3), plot(x(Pty1in), y(Pty1in), 'go', 'markers',...
MarkerSize,'MarkerFaceColor','g');
subplot(3,2,3),quiver(x, y, Vx * ScaleV, Vy * ScaleV, 0, 'r', 'linewidth', 2);
hold off
axis(Limits);
title('Atoms and Velocities')
xlabel('X')
ylabel('Y')

AE0(Pty0in) = 1 / 2 * Mass0 * V2(Pty0in) + Phi(Pty0in) - Phi0;
AE1(Pty1in) = 1 / 2 * Mass1 * V2(Pty1in) + Phi(Pty1in) - Phi0;
AE1=AE1(AE1~=0);
% if AddParticle, AE(end) = 0;end
% AE = Phi;
num0=length(AE0);
num1=length(AE1);
subplot(3, 2, 4),scatter3(x(1:num0),y(1:num0),AE0,ones(1,num0)*300,AE0,'fill');
hold on
subplot(3, 2, 4),scatter3(x(num0+1:end),y(num0+1:end),AE1,ones(1,num1)*100,AE1,'fill');
axis(Limits);
view(2);
title('Atom Energy (KE + Pot)')
hold off

subplot(3, 2, 5), plot(x(Pty0in), y(Pty0in), 'bo', 'markers',...
MarkerSize, 'MarkerFaceColor', 'b');
hold on
subplot(3, 2, 5), plot(x(Pty1in), y(Pty1in), 'go', 'markers',...
MarkerSize, 'MarkerFaceColor', 'g');
hold off
title('Atoms (All)')
xlabel('X')
ylabel('Y')

subplot(3, 2, 6), plot(time, T, 'k', 'linewidth', 2);
hold on
subplot(3, 2, 6), plot(time, T0, 'b', 'linewidth', 2);
subplot(3, 2, 6), plot(time, T1, 'g', 'linewidth', 2);

xlabel('time')
ylabel('Temperature')
title('Temperature: Each type and all')

if doPlotImage == 1
if PlotCount == 1
f = getframe(PlotFig);
[im, map] = rgb2ind(f.cdata, 256, 'nodither');
im(1, 1, 1, 2) = 0;

else
f = getframe(PlotFig);
im(:, :, 1, PlotCount) = rgb2ind(f.cdata, map, 'nodither');
end
end

end