MATLAB Answers

Finding values within a matrix

4 views (last 30 days)
Tomás Costa
Tomás Costa on 11 Nov 2019
Commented: Tomás Costa on 11 Nov 2019
Hi, im not sure if I'm asking the right question here but I'm having trouble finding an answer on why this happens. In the program bellow, the question lies on the "if". Basically the program gives me a matrix of 301x8x26 of temperatures for a give (Tc, P0, T0). When i do the if im trying to compare and find in the matrix, if there is a value that equals Tc(i) (for the same index) and it gives a matrix full of zeros. However, if i try to print, for example Tc(279)==T22(279,8,26), the answer gives me 1... What im trying to do here is create a matrix with all the values that correspond to one another and then plot P0,T0 and those temperatures as a function of the first 2. If you could help, it would mean a lot. Thank you very much
clc
clear all
tic
file = 'areatemp1';
file1 = 'machtemp1';
opts = detectImportOptions(file,'NumHeaderLines',0);
opts1 = detectImportOptions(file1,'NumHeaderLines',0);
T = readtable(file,opts);
T.Properties.VariableNames = {'Area','Mach','Temp'};
t = readtable(file1,opts1);
t.Properties.VariableNames = {'mach','temp'};
%condições iniciais do problema
Tc = 20:0.1:50;
P0 = 18:25;
T0 = 25:50;
D0 = input ('What is the size of the hole in microns? ');
D0 = D0*10^-6;
A0 = ((pi*D0)^2)/4;
i = 0;
j = 0;
k = 0;
Patm = 1;
R = 8.314;
A01 = (1/5)*A0;
gama = 1.4;
for i = 1:length(Tc);
% Cálculo do rho para a temperatura prevista de saida
rhoe(i) = Patm/(R*Tc(i));
for j = 1:length(P0);
% Calculo inicial de P01 para dar inicio á interpolação
P011(j) = 0.666*P0(j) + 0.333*Patm;
for k = 1:length(T0);
mass(j,k) = (0.6856*P011(j)*A01)/(R*T0(k))^(1/2);
rho0(j,k) = 0.9805 * (P011(j)/(R*T0(k)));
V0(j,k) = mass(j,k)/(rho0(j,k)*A0);
P01(j,k) = P0(j)-(1/4)*(rho0(j,k)*V0(j,k)^2);
T1(k) = T0(k) * 0.8316;
V1(k) = (gama*R*T1(k))^(1/2);
P1(k) = P01(j,k) * 0.5274;
Areas(i,j,k) = 1/((((R*T0(k))^(1/2))*rhoe(i)*V1(k))/(1.3712*P01(j,k))+sqrt((((R*T0(k))^(1/2)*rhoe(i)*V1(k))/(1.3712*P01(j,k)))^2+(R*T0(k)*rhoe(i)*(P1(k)-Patm))/(0.6856*P01(j,k))^2));
%interpolaton
dif = Areas(i,j,k); %obtained however you want. Can be an array of temperatures
Tarray = table2array(T);
tarray = table2array(t);
machs(i,j,k) = interp1(Tarray(:,1), Tarray(:,2), dif);
temps(i,j,k) = interp1(Tarray(:,1), Tarray(:,3), dif);
t1(i,j,k) = temps(i,j,k)*T0(k);
temperatura(i,j,k) = interp1(tarray(:,1), tarray(:,2), machs(i,j,k));
T2(i,j,k) = temperatura(i,j,k)*t1(i,j,k);
T22(i,j,k) = round(T2(i,j,k),1);
end
% end
end
end
if Tc(i) == T22(i,j,k)
Tfinal(i,j,k) = 1;
else
Tfinal(i,j,k) = 0;
end
% P0 = linspace(18,25,numel(T22));
% T0 = linspace(25,50,numel(T22));
% plot (P0,T22,'b--o',T0,T22,'r--*')
toc

Accepted Answer

Guillaume
Guillaume on 11 Nov 2019
It appears that you are building your code from bits you've asked or found without really understanding how it works. In particular the "%obtained however you want. Can be an array of temperatures" looks like a comment I could have written, have I helped you before?
It makes zero sense to perform the interpolations in the loop. It should be done afterward. Also, you need to learn to work with tables. There's never any reason to convert them to array. And don't differentiate variables just based on case. Using T and t is just asking for trouble, same for Tarray and tarray, you can be sure that at one point you'll mispell one for the other.
Your code doesn't need a single loop:
%all constants as defined in your code
PO = 18:25; %ROW vector. THE DIRECTION IS IMPORTANT.
Tc = (20:0.1:50).'; %make TC a COLUMN vector. THE DIRECTION IS IMPORTANT. MUST be different from P0
T0 = permute(25:50, [1 3 2]); %make T0 a vector in the 3RD dimension. THE DIRECTION IS IMPORTANT. Must be different from P0 and Tc
%calculate everything at once. No loop needed:
rhoe = Patm./(R*Tc); %note that ./ is required since we're working with matrices/vectors
P011 = 2*P0/3 + Patm/3;
mass = 0.6856*P011*A01./sqrt(R*T0); %sqrt is simpler than ^(1/2)
rho0 = 0.9805 * P011./(R*T0);
V0 = mass ./ (rho0*A0);
P01 = P0 - rho0.*V0^2 / 4;
T1 = T0 * 0.8316;
V1 = sqrt(gama*R*T1);
P1 = P01 * 0.5274;
Areas = 1./(sqrt(R*T0).*rhoe.*V1 ./ (1.3712*P01) + sqrt((sqrt(R*T0).*rhoe.*V1 ./ (1.3712*P01)).^2 + (R*T0.*rhoe.*(P1-Patm))./(0.6856*P01)^2));
machs = interp1(T.Area, T.Mach, Areas); %A lot clearer than converting tables to arrays
temps = interp1(T.Area, T.Temp, Areas);
t1 = temps .* T0; %t1 not to be confused with T1. Differentiating variables based on case is a BAD IDEA!
temperatura = interp1(t.mach, t.temp, machs);
T2 = temperatura .* t1;
T22 = round(T2, 1);
I may have made some typos simplifying your code. How about calling T MachTempvsArea, and t TempvsMach? Much better variable names that actually explain what is in the table.
As for your comparison, it's not clear what you're trying to do. At the point you're using i, you're no longer in the i loop, so Tc(i) is always the last value of Tc, and T(i,j,k) is always T(end, end, end). No idea what you're trying to achieve here.
  3 Comments
Tomás Costa
Tomás Costa on 11 Nov 2019
That is perfect! Thank you so much

Sign in to comment.

More Answers (1)

Fabio Freschi
Fabio Freschi on 11 Nov 2019
It is difficult to answer because we don't have the data files 'areatemp1' and 'machtemp1'.
In any case in your if statement you are comparing Tc(i) == T22(i,j,k), i,j,k are here fixed values and in particular they assume the values that they have at the end of the loops
i = length(Tc);
j = length(P0);
k = length(T0);
So you are actually comparing only the last value of your T22 matrix with the last value of your Tc array.I think you should put your if statement inside the innermost loop.
  1 Comment
Fabio Freschi
Fabio Freschi on 11 Nov 2019
If the previous answer is correct, you can further clean-up your code by
  1. removing the initialization of i,j,k to 0 (useless),
  2. removing the semicolon at the end of the for statements, and
  3. preallocating the arrays/matrices inside your loops

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!