How do I plot an integral which is repeated over an interval?

I am attempting to plot the absporptance (alpha_s_bar) of an uncontaminated spacecraft as a function of the incident light wavelenght. The absorptance is found by integrating two different functions of wavelength over an interval from zero to the wavelength value of interest, on an interval from 0 to 0.7 microns. The code I have attached below runs and returns different values of absorptance as the wavelenght of interest (L) changes, but does not save the values before iterating, so I end up with one scalar value as alpha_s_bar. I am interested in saving the values of alpha_s_bar as a vector so I can plot them against the wavelength.

%Note that I will be using x for the wavelength for ease of typing

% Start by defining the constants for the absorptance for typical

% outgassing products, called alpha_c in problem statement, where

% alpha_c=M10*exp(M11*L)

M10=34.898;

M11=-8.6481;

% Next define the constants for solar absorptivity, alpha_s, where

% aplha_s=M0+M1*L+M2*L^2+M3*L^3

M0=0.28779;

M1=0.06028;

M2=-0.046055;

M3=0.020211;

% Next, define the constants for the solar spectrum (adjusted so units are microns), S, where

% S=2*C1/L^5[exp(C2/L*T)-1]

T=5500;

C1=5.594E-7;

C2=1438.8;

for L=0:0.01:0.7;

fnum=@(x) (M0+M1.*x+M2.*(x.^2)+M3.*(x.^3))./((x.^5).*exp(C2./(x.*T)));

fden=@(x) 1./((x.^5).*exp(C2./(x.*T)));

alpha_bar_s=integral(fnum,0,L)./integral(fden,0,L)

plot(L,alpha_bar_s)

hold on

end

Answer 1

[Abhas]

Hi Mitchell,

To save the values of "alpha_bar_s" as a vector and plot them against the wavelength, you need to store each computed value in an array during the loop. After the loop, you can plot the results.

Here's the modified MATLAB code to achieve the same:

% Define constants
M10 = 34.898;
M11 = -8.6481;
M0 = 0.28779;
M1 = 0.06028;
M2 = -0.046055;
M3 = 0.020211;
T = 5500;
C1 = 5.594E-7;
C2 = 1438.8;

% Initialize vectors to store results
L_values = 0:0.01:0.7;
alpha_bar_s_values = zeros(size(L_values));

% Loop over wavelength values
for i = 1:length(L_values)
    L = L_values(i);
    
    if L == 0
        % Handle the case where L is zero
        alpha_bar_s_values(i) = 0;
    else
        % Define the integrand functions
        fnum = @(x) (M0 + M1.*x + M2.*(x.^2) + M3.*(x.^3)) ./ ((x.^5) .* exp(C2 ./ (x .* T)));
        fden = @(x) 1 ./ ((x.^5) .* exp(C2 ./ (x .* T)));
        
        % Calculate alpha_bar_s for each L
        alpha_bar_s_values(i) = integral(fnum, 0.001, L) / integral(fden, 0.001, L);
    end
end

% Plot the results
plot(L_values, alpha_bar_s_values, 'b-', 'LineWidth', 2);
xlabel('Wavelength (microns)');
ylabel('Absorptance (\alpha_{s\_bar})');
title('Absorptance as a Function of Wavelength');
grid on;

You can address the situation when "L" is zero according to your needs. If unhandled, this might trigger a warning indicating that the integrand functions "fnum" or "fden" are producing "Inf" or "NaN" values at certain points, likely due to division by zero or overflow. This issue generally arises when the wavelength "L" is zero, as the functions involve division by "x^5", which becomes problematic when "x" is zero.

You may refer to the below MathWorks documentation links to have a better understanding on numerical integration and function handles:

1. https://www.mathworks.com/help/matlab/ref/integral.html
2. https://www.mathworks.com/help/matlab/matlab_prog/creating-a-function-handle.html

Answer 2

[Voss]

Make L a vector, use one element of L on each iteration of the for loop, pre-allocate alpha_bar_s before the loop to be the same size as L, and store one element of alpha_bar_s on each iteration. Also, since the definitions of fnum and fden don't depend on the value of L, they can be set once before the loop.

M10=34.898;
M11=-8.6481;
% Next define the constants for solar absorptivity, alpha_s, where
% aplha_s=M0+M1*L+M2*L^2+M3*L^3
M0=0.28779;
M1=0.06028;
M2=-0.046055;
M3=0.020211;
% Next, define the constants for the solar spectrum (adjusted so units are microns), S, where
% S=2*C1/L^5[exp(C2/L*T)-1]
T=5500;
C1=5.594E-7;
C2=1438.8;
L=0:0.01:0.7; 
fnum=@(x) (M0+M1.*x+M2.*(x.^2)+M3.*(x.^3))./((x.^5).*exp(C2./(x.*T)));
fden=@(x) 1./((x.^5).*exp(C2./(x.*T)));
alpha_bar_s = zeros(size(L));
for ii = 1:numel(L) 
  alpha_bar_s(ii)=integral(fnum,0,L(ii))./integral(fden,0,L(ii));
end
plot(L,alpha_bar_s)

[Mitchell]

Thank You!

[Voss]

You're welcome! Any questions please let me know. Otherwise please Accept this answer. Thanks!