MATLAB：我正在尝试在色度图上绘制普拉克曲线。有人可以帮我修复我的代码吗？

Question

我有一项任务是在 Matlab 中的色度图上绘制普朗克曲线。给定三个文本文件：cie_x2j.txt、cie_y2j.txt 和 cie_z2j.txt。

我以前从未使用过 Matlab，所以非常感谢任何帮助！

情节应该是这样的： 但是我的代码现在看起来像这样：.

这是我的代码：

函数 s_PlanckCurve()

`%%% -----  Code ----- %%%

clear all
close all

% --- Plot the chromaticity diagram as in the previous assignment --- %

% Load data from a file
cie_x2j = load('cie_x2j.txt');
cie_y2j = load('cie_y2j.txt');
cie_z2j = load('cie_z2j.txt');

% Define the wavelength range (380-780 nm)
wavelengths = linspace(380, 780, numel(cie_x2j));

% Compute the chromaticity coordinates
x = cie_x2j ./ (cie_x2j + cie_y2j + cie_z2j);
y = cie_y2j ./ (cie_x2j + cie_y2j + cie_z2j);


% Define the desired wavelengths for labels
label_wavelengths = [380 470 480 490 500 510 520 530 540 550 560 570 580 590 600 610 620 780];

% Find the corresponding indices in the wavelength array
label_indices = arrayfun(@(w) find(wavelengths == w, 1), label_wavelengths);


% Plot the chromaticity diagram
figure;
plot(x, y, 'LineWidth', 2);
hold on;
% Plot the triangle x + y + z = 1
triangle_x = [1, 0, 0, 1];
triangle_y = [0, 1, 0, 0];
plot(triangle_x, triangle_y, 'k-', 'LineWidth', 1);
% Add labels for each desired wavelength
for i = 1:numel(label_wavelengths)
    label_x = x(label_indices(i));
    label_y = y(label_indices(i));
    label_text = sprintf('%d nm', label_wavelengths(i));
    text(label_x, label_y, label_text, 'FontSize', 7, 'HorizontalAlignment', 'left', 'VerticalAlignment', 'bottom');
end

% Connect the coordinates for 380 nm and 780 nm with a purple line
plot([x(1), x(end)], [y(1), y(end)], 'Color', [0.5 0 0.5], 'LineWidth', 2);

% Set axis labels and title
xlabel('x');
ylabel('y');
title('Chromaticity Diagram');
axis equal;
grid on;




% --- Calculate the spectral distribution of the radiance of a black body (Le) --- %

lambda = 380:50:780;
T = 1000:1125:10000;

omega_0 = 1;

k = 1.3806488E-23;
c = 2.99792458E8;
h = 6.62606957E-34;
c_1 = 2*pi*h*c^2;

C_2 = (c*h) ./ k;

% Ensuring lambda and T are in the correct dimensions for the calculation
[Lambda, T] = meshgrid(lambda, T);

Le = (c_1 ./ (pi*Lambda.^5)) .* (1 ./ (exp((C_2)./(Lambda.*T)) - 1)) .* (1 ./ omega_0);




% --- Plot the planckian curve on the chromaticity diagram --- %

% Interpolate the color matching functions
cie_x_interp = interp1(wavelengths, cie_x2j, lambda);
cie_y_interp = interp1(wavelengths, cie_y2j, lambda);
cie_z_interp = interp1(wavelengths, cie_z2j, lambda);

% Calculate the chromaticity coordinates (x, y) for each temperature
planckian_x = zeros(size(T));
planckian_y = zeros(size(T));

for i = 1:length(T)
    % Calculate the blackbody spectrum for the current temperature
    bb_spectrum = (1 ./ Lambda.^5) .* (1 ./ (exp((C_2) ./ (Lambda .* T(i))) - 1));
    
    % Normalize the blackbody spectrum
    bb_spectrum = bb_spectrum / max(bb_spectrum(:));
    
    % Calculate the tristimulus values (XYZ) for the blackbody spectrum
    X_temp = sum(bb_spectrum .* cie_x_interp, 2);
    Y_temp = sum(bb_spectrum .* cie_y_interp, 2);
    Z_temp = sum(bb_spectrum .* cie_z_interp, 2);
    
    % Calculate the chromaticity coordinates (x, y) for the blackbody spectrum
    planckian_x(i) = X_temp(i) / (X_temp(i) + Y_temp(i) + Z_temp(i));
    planckian_y(i) = Y_temp(i) / (X_temp(i) + Y_temp(i) + Z_temp(i));
end

% Plot the Planckian curve on the chromaticity diagram
hold on;
plot(planckian_x, planckian_y, 'r-', 'LineWidth', 2);
legend('Spectral Locus', 'x + y + z = 1', 'Planckian Curve', 'Location', 'Best');The plot should look like this:
[enter image description here](https://i.stack.imgur.com/QT52j.png) But with my code now looks like this: [enter image description here](https://i.stack.imgur.com/Qvw40.png).

I have never worked with Matlab before, so any help is very much appreciated!

This is my code:

function s_PlanckCurve()

%%% -----  Code ----- %%%

clear all
close all

% --- Plot the chromaticity diagram as in the previous assignment --- %

% Load data from a file
cie_x2j = load('cie_x2j.txt');
cie_y2j = load('cie_y2j.txt');
cie_z2j = load('cie_z2j.txt');

% Define the wavelength range (380-780 nm)
wavelengths = linspace(380, 780, numel(cie_x2j));

% Compute the chromaticity coordinates
x = cie_x2j ./ (cie_x2j + cie_y2j + cie_z2j);
y = cie_y2j ./ (cie_x2j + cie_y2j + cie_z2j);


% Define the desired wavelengths for labels
label_wavelengths = [380 470 480 490 500 510 520 530 540 550 560 570 580 590 600 610 620 780];

% Find the corresponding indices in the wavelength array
label_indices = arrayfun(@(w) find(wavelengths == w, 1), label_wavelengths);


% Plot the chromaticity diagram
figure;
plot(x, y, 'LineWidth', 2);
hold on;
% Plot the triangle x + y + z = 1
triangle_x = [1, 0, 0, 1];
triangle_y = [0, 1, 0, 0];
plot(triangle_x, triangle_y, 'k-', 'LineWidth', 1);
% Add labels for each desired wavelength
for i = 1:numel(label_wavelengths)
    label_x = x(label_indices(i));
    label_y = y(label_indices(i));
    label_text = sprintf('%d nm', label_wavelengths(i));
    text(label_x, label_y, label_text, 'FontSize', 7, 'HorizontalAlignment', 'left', 'VerticalAlignment', 'bottom');
end

% Connect the coordinates for 380 nm and 780 nm with a purple line
plot([x(1), x(end)], [y(1), y(end)], 'Color', [0.5 0 0.5], 'LineWidth', 2);

% Set axis labels and title
xlabel('x');
ylabel('y');
title('Chromaticity Diagram');
axis equal;
grid on;




% --- Calculate the spectral distribution of the radiance of a black body (Le) --- %

lambda = 380:50:780;
T = 1000:1125:10000;

omega_0 = 1;

k = 1.3806488E-23;
c = 2.99792458E8;
h = 6.62606957E-34;
c_1 = 2*pi*h*c^2;

C_2 = (c*h) ./ k;

% Ensuring lambda and T are in the correct dimensions for the calculation
[Lambda, T] = meshgrid(lambda, T);

Le = (c_1 ./ (pi*Lambda.^5)) .* (1 ./ (exp((C_2)./(Lambda.*T)) - 1)) .* (1 ./ omega_0);




% --- Plot the planckian curve on the chromaticity diagram --- %

% Interpolate the color matching functions
cie_x_interp = interp1(wavelengths, cie_x2j, lambda);
cie_y_interp = interp1(wavelengths, cie_y2j, lambda);
cie_z_interp = interp1(wavelengths, cie_z2j, lambda);

% Calculate the chromaticity coordinates (x, y) for each temperature
planckian_x = zeros(size(T));
planckian_y = zeros(size(T));

for i = 1:length(T)
    % Calculate the blackbody spectrum for the current temperature
    bb_spectrum = (1 ./ Lambda.^5) .* (1 ./ (exp((C_2) ./ (Lambda .* T(i))) - 1));
    
    % Normalize the blackbody spectrum
    bb_spectrum = bb_spectrum / max(bb_spectrum(:));
    
    % Calculate the tristimulus values (XYZ) for the blackbody spectrum
    X_temp = sum(bb_spectrum .* cie_x_interp, 2);
    Y_temp = sum(bb_spectrum .* cie_y_interp, 2);
    Z_temp = sum(bb_spectrum .* cie_z_interp, 2);
    
    % Calculate the chromaticity coordinates (x, y) for the blackbody spectrum
    planckian_x(i) = X_temp(i) / (X_temp(i) + Y_temp(i) + Z_temp(i));
    planckian_y(i) = Y_temp(i) / (X_temp(i) + Y_temp(i) + Z_temp(i));
end

% Plot the Planckian curve on the chromaticity diagram
hold on;
plot(planckian_x, planckian_y, 'r-', 'LineWidth', 2);
legend('Spectral Locus', 'x + y + z = 1', 'Planckian Curve', 'Location', 'Best');`