MATLAB: Subscript indices must be real positive or logicals

subscript indices

      %Astro Final: Problem 2b Competition of Proton Capture and Beta
  %Decay
  %Carrie Elliott - 2015
T = linspace(0.1, 1, 100);
                        %Establishes range(T) in units of 10^9 Kelvin
p1B = -4.538;           %only parameter for Beta decay
Rbeta = exp(p1B);        %The rate of Beta decay is constant, and has only one 
                        %parameter independent of temperature.
p1_p1 = 17.29;          %These are the seven parameters

p2_p1 = 4.341*10^-4;    %for the first component of proton capture;
p3_p1 = -18.10;         % _p1 denotes this.
p4_p1 = 0.1471;         %







p5_p1 = -0.1550;        %
p6_p1 = 0.0167;         %
p7_p1 = -0.7411;        %
Rp1 = exp(p1_p1 + (p2_p1\T) + (p3_p1\((T).^(1\3))) + p4_p1(T.^1\3) + p5_p1*T + p6_p1*(T.^(5\3)) + p7_p1(ln(T)));
p1_p2 = 3.428;          %These are the seven parameters
p2_p2 = -4.654;         % for the second component of proton capture;
p3_p2 = -3.999;         % _p2 denotes this.    
p4_p2 = 8.947;          %       
p5_p2 = -0.6332;        %
p6_p2 = 0.0402;         %
p7_p2 = -5.113;         %
Rp2 = exp(p1_p2 + (p2_p2/T) + (p3_p2/((T)^(1/3))) + p4_p2(T^1/3) + p5_p2*T + p6_p2*(T^(5/3)) + p7_p2(ln(T)));
Rproton = Rp1 + Rp2;
hold on
plot(Rbeta)
plot(Rproton)
title('Competition between Beta+ Decay and Proton Capture as Function of Temperature');
xlabel('Temp, 10^9 Kelvin');
ylabel('Rate');
hold off %
What's going on here? I'm simply trying to plot two rates of reactions over a range of temperatures (0.1*10^9 – 1.0*10^9) on the same plot, and it says the error is in line 21 but I have no idea where or why.

Best Answer

This is probably the first instance that throws the error:
p4_p1(T.^1\3)
but p4_p1 is a scalar constant, not an array:
p4_p1 = 0.1471;         %
What do you want to do in that line? My guess is that you want to multiply where you specified no operator, so see if substituting this line solves the problem:
Rp1 = exp(p1_p1 + (p2_p1\T) + (p3_p1\((T).^(1\3))) + p4_p1*(T.^1\3) + p5_p1*T + p6_p1*(T.^(5\3)) + p7_p1*(log(T)));
and:
Rp2 = exp(p1_p2 + (p2_p2./T) + (p3_p2/((T).^(1/3))) + p4_p2(T.^1/3) + p5_p2*T + p6_p2*(T.^(5/3)) + p7_p2*(log(T)));
Also, the natural logarithm function is log not ln, and the two division operators ‘/’ and ‘\’ are not the same, but you seem to have used one in Rp1 in the same place you used the other in Rp2. I did not change the division operators because I don’t know which one you want.
Related Question