Untitled
unknown
plain_text
2 months ago
3.7 kB
1
Indexable
% Constants a = 3e-10; % Lattice constant in meters b = 3e-10; % Potential width in meters V0 = 1.60218e-19; % Potential height in Joules hbar = 1.0545718e-34; % Reduced Planck's constant in J·s m = 9.10938356e-31; % Electron mass in kg % Define a "safe" sine function that avoids division by zero. % For very small arguments (|x| < tol), sin(x)/x is approximated by 1. tol = 1e-8; safe_sin = @(x) (abs(x) < tol).*1 + (abs(x) >= tol).*(sin(x)./x); % Define the range of k values (crystal momentum) k = linspace(-pi/a, pi/a, 1000); % Initialize arrays to store energies for two bands E1 = NaN(1, length(k)); % First band energy E2 = NaN(1, length(k)); % Second band energy % Loop over each k value for i = 1:length(k) % For the current k, define the function f(k1) that must be zero: % f(k1) = cos(k(i)*a) - (m*V0*b*a/hbar^2)*[safe_sin(k1*a)] - cos(k1*a) func = @(k1) cos(k(i)*a) - (m*V0*b*a/hbar^2)*safe_sin(k1*a) - cos(k1*a); % Define a range for k1 (avoid k1 = 0 to steer clear of division by zero) k1_range = linspace(0.1, 2*pi/a, 100); % 100 points from 0.1 to 2*pi/a k1_solutions = []; % To store found roots for this k % Bracket the root by checking consecutive intervals for j = 1:length(k1_range)-1 x1 = k1_range(j); x2 = k1_range(j+1); f1 = func(x1); f2 = func(x2); % Only proceed if both endpoint values are finite... if ~isfinite(f1) || ~isfinite(f2) continue; end % ...and there is a sign change (or one of the endpoints is zero) if f1 * f2 <= 0 try root = fzero(func, [x1, x2]); % Add the root if it is finite and not already in the list (within tolerance) if isfinite(root) if isempty(k1_solutions) || ~any(abs(k1_solutions - root) < 1e-8) k1_solutions = [k1_solutions, root]; end end catch continue; end end end % Sort the roots and, if at least two are found, record the two lowest-energy solutions. k1_solutions = sort(k1_solutions); if length(k1_solutions) >= 2 E1(i) = (hbar^2 * k1_solutions(1)^2) / (2*m); % First band energy E2(i) = (hbar^2 * k1_solutions(2)^2) / (2*m); % Second band energy end end % Calculate overall band statistics using only valid k values (where both bands were found) validIndices = ~isnan(E1) & ~isnan(E2); k_valid = k(validIndices); E1_valid = E1(validIndices); E2_valid = E2(validIndices); % Find the minimum and maximum energies for the first two bands (in Joules) E1_min = min(E1_valid); E1_max = max(E1_valid); E2_min = min(E2_valid); E2_max = max(E2_valid); % Calculate band widths (in Joules) band_width_1 = E1_max - E1_min; band_width_2 = E2_max - E2_min; % Calculate energy gap between the first and second bands energy_gap = max(0, E2_min - E1_max); % Display results in electron-volts (eV) fprintf('Width of the first energy band: %.4f eV\n', band_width_1 / 1.60218e-19); fprintf('Width of the second energy band: %.4f eV\n', band_width_2 / 1.60218e-19); fprintf('Energy gap between the first and second bands: %.4f eV\n', energy_gap / 1.60218e-19); % ------------------------- % Plot the energy bands vs. k % Convert energies from Joules to eV for plotting. E1_eV = E1_valid / 1.60218e-19; E2_eV = E2_valid / 1.60218e-19; figure; plot(k_valid, E1_eV, 'b-', 'LineWidth', 2); hold on; plot(k_valid, E2_eV, 'r-', 'LineWidth', 2); xlabel('Crystal Momentum k (1/m)'); ylabel('Energy (eV)'); title('Energy Bands'); legend('Band 1', 'Band 2', 'Location', 'Best'); grid on;
Editor is loading...
Leave a Comment