MATLAB编写的相场程序，模拟枝晶的凝固过程.zip

function X = Untitled3(sys_length, t_max, deposition_rate, rule) % ---------------------------------------------------------------------------------- % Created by Henry C. Nugroho % For MATH 489 -- Research Interaction in Mathematics % Program SimulateDeposition(sys_length, t_max, deposition_rate, rule) % Parameters: % sys_length = size of system (in integer only) % t_max = maximum units of time to simulate % deposition_rate = numbers of atoms deposited per unit time % rule = indicates which deposition rule to use % if rule = 1, Simple Rule % an atom coming always sticks % on top of the other atom regardless % the height of its neighbors. % h(i,t+1) = h(i,t) + 1 % if rule = 2, Ballistic Deposition Rule % we increase h(i,t) to h(i,t+1) by % MAX[h(i-1,t),h(i,t)+1,h(i+1,t)] % % Example: In the command window, type 'SimulateDeposition(30,100,30,2) % Meaning: The length/size of the system is 30, simulating within % 100 units of time at deposition rate of 30 atoms per unit time % following Ballistic Deposition rule. % % Results: The color printout of the graphs are stored in WidthFluc_1D.jpg and % Surface_Plot1D.jpg. You can alter this by changing line 95 & 111 % --------------------------------------------------------------------------------- clc close all % initialization format long; sys_length=30; t_max=100; deposition_rate=30; L = sys_length; rate = deposition_rate; rule=2; temp_max = 0; % dummy variable MaxHeight = 100; for i = 1:L h(i) = 0; end for i = 1:t_max width(i) = 0; end % generating an array of random numbers A = ceil(L*rand(t_max,rate)); % the simulation if((rule == 1) || (rule == 2)) if (rule == 1) rule_name = 'Simple'; else rule_name = 'Ballistic'; end for t = 1: t_max for k = 1 : rate i = A(t,k); if (rule == 1) h(i) = h(i) + 1; hit(i,h(i)) = 1; else if (i == 1) h(i) = max(h(i)+1,h(i+1)); hit(i,h(i)) = 1; elseif (i == L) h(i) = max(h(i-1),h(i)+1); hit(i,h(i)) = 1; else temp_max = max(h(i-1),h(i+1)); h(i) = max(temp_max,h(i)+1); hit(i,h(i)) = 1; end end end width(t) = std(h,1); % Printing h(i,t) on the screen to verify the results % (note: You can uncomment the following 3 lines if you want to check) % fprintf('Time : %d\n', t); % fprintf('Width Fluctuation : %d\n', width(t)); % h end % Plotting the width as a function of time using semilogx function figure('Name','Width Fluctuation as a Function of Time in 1 Dimensional Lattice'); semilogx(width(1:t_max),'s',... 'MarkerEdgeColor','r',... 'MarkerFaceColor','b',... 'MarkerSize',4) xlabel('Log10(Time)','FontSize',12) ylabel('W(L,t)','FontSize',12) title(['L = ',num2str(L),' Time Max = ',num2str(t_max), ' Rate = ', num2str(rate), ' Rule = ', rule_name],'FontSize',12) grid on print -djpeg WidthFluc_1D.jpg; % Plotting the surface figure for iatom=1:L for jatom=h(iatom):-1:1 if (hit(iatom,jatom)==1) plotatom(iatom,jatom,'b') end hold on end end axis([1 L 0 (max(h)+10)]) xlabel('column i','FontSize',12) ylabel('h(i)','FontSize',12) title(['Surface Plot of L = ', num2str(L)],'FontSize',12) print -djpeg Surface_Plot1D.jpg; else fprintf('\nInvalid Rule\n'); fprintf('\nRule = 1 for Simple Rule where h(i,t+1) = h(i,t) + 1\n'); fprintf('\nRule = 2 for Ballistic Rule where h(i,t+1) = MAX[h(i-1,t),h(i,t)+1,h(i+1,t)]\n'); end %-------------------------------------------------------------------------% function plotatom(x0,y0,color) % plot one atom in the specified color, centered at (x0,y0) with radius r r = 0.5; theta = 0:0.1:6.29; x = x0 + r*cos(theta); y = y0 + r*sin(theta); fill(x,y,color) plot(x,y,color)