Finding values within a matrix

Asked by Tomás Costa on 11 Nov 2019

Latest activity Commented on by Tomás Costa on 11 Nov 2019

4 views (last 30 days)

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

0 Comments

2 Answers

Link

Answer by Guillaume on 11 Nov 2019

Accepted Answer

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 on 11 Nov 2019

You are right, im not a matlab genius and I picked up the interpolation part of the code from an answer that had been given to a question someone made. Im yet to clean the code. Basically what im trying to achieve with this code is a plot with P0 (x-axis) and T0 (y-axis). However the values that i want to put in, are not so simple; i want the values of P0 and T0 for which the values of T22 are equal to the corresponding value of Tc. So i basically have 3 inputs, P0,T0 and Tc and an output T22. For the values of T22 that are equal to the value of Tc (used for its calculation) save the values of T0 and P0 and make a plot thats why i was using the For cycles. I haven't reached this part yet and to be honest i really dont know how to do this so if you could help me would be amazing. Thank you for your help so far.

Guillaume on 11 Nov 2019

First issue: "For the values of T22 that are equal to", what does that mean exactly bearing in mind that:

>> (0.1 + 0.1 + 0.1) == 0.3
ans =
  logical
   0

is false.

Once you have defined this (typically you'd compare the absolute difference of two numbers against an arbitrary tolerance), then it's trivial to filter the result:

issame = abs(T22 - Tc) < 1e-8;  %No idea if that's what you want. T22 and Tc are equal when their difference is less than 1e-8
%define P0, Tc and T0 as 3d matrices to make them easier to search:
[P0_full, Tc_full, T0_full] = ndgrid(P0, Tc, T0);  %do not change the order of the arguments as they match the order of dimensions
selectedresults = table(P0_full(issame), Tc_full(issame), T0_full(issame), 'VariableNames', {'P0', 'Tc', 'T0'})

Tomás Costa on 11 Nov 2019

That is perfect! Thank you so much

Link

Answer by 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 on 11 Nov 2019

If the previous answer is correct, you can further clean-up your code by

removing the initialization of i,j,k to 0 (useless),
removing the semicolon at the end of the for statements, and
preallocating the arrays/matrices inside your loops

You are now following this question

Tomás Costa (view profile)