Storing Numbers from Function in a Script

I have a function designed to output thermodynamic quantities. I created it as a function, as I felt that was easiest at the time. I am trying to call the function in a script, which works successfully. My function is outputting numbers, and the right numbers. However, I want to vary T and p over a range, and store these numbers, and when I tried to store the numbers, all I got was 1X a large number complex double, and all rows and columns were 0. I will put my current code below (the one that calls the function successfully but does not store the numbers) and below that, I can put what I was trying that gave me the complex double with all zeroes.

Script that successfully calls the function and produces accurate numbers that aren't stored;

z = 1;
p = 100000:1:10000000;
T = 153:1:423;
pc = 5050000;
Tc = 154.6;
w = 0.022;
k = 0;
cpig = [26 0.01 0.005 0.0000013];
DHf = 0;
DGf = 0;
ph = 'V';
flag = 'PR';

for T = 153:1:423
    for p = 100000:1:10000000
        [rho,Z,phi,A,S,H,U,G] = prsrk(z,p,T,pc,Tc,w,k,cpig,DHf,DGf,ph,flag);
    end
end

Code that gives me complex doubles full of zeroes;

z = 1;
p = 100000:1:10000000;
T = 153:1:423;
pc = 5050000;
Tc = 154.6;
w = 0.022;
k = 0;
cpig = [26 0.01 0.005 0.0000013];
DHf = 0;
DGf = 0;
ph = 'V';
flag = 'PR';

IntEng = [];
Gibbs = [];
Enthalpy = [];
Entropy = [];

for T = 153:1:423
    for p = 100000:1:10000000
        [rho,Z,phi,A,S,H,U,G] = prsrk(z,p,T,pc,Tc,w,k,cpig,DHf,DGf,ph,flag);
        Gibbs(T,p) = G;
        IntEng(T,p) = U;
        Enthalpy(T,p) = H;
        Entropy(T,p) = S;
    end
end

[Rik]

Take a look here: https://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

Correctly formatting your code makes it easier for other to see what happens and to spot potential errors.

What I'm noticing is that you are overwriting the values rho,Z,phi,A,S,H,U and G. So the thing you can try is changing that to rho(T_index,p_index) etc.

[Jerome Campbell]

Apologies. Here is what I have tried. I am still getting the problem where Gibbs, IntEng, Enthalpy, and Entropy are all 153x12000 complex doubles (I stop the code to see if it is working, otherwise it would be larger.)

z = 1;
p = 100000:1:10000000;
T = 153:1:423;
pc = 5050000;
Tc = 154.6;
w = 0.022;
k = 0;
cpig = [26 0.01 0.005 0.0000013];
DHf = 0;
DGf = 0;
ph = 'V';
flag = 'PR';

IntEng = [];
Gibbs = [];
Enthalpy = [];
Entropy = [];

for T = 153:1:423
    for p = 100000:1:10000000
        [rho(T,p),Z(T,p),phi(T,p),A(T,p),S(T,p),H(T,p),U(T,p),G(T,p)] ...
        = prsrk(z,p,T,pc,Tc,w,k,cpig,DHf,DGf,ph,flag)
        % 
        Gibbs(T,p) = G;
        IntEng(T,p) = U;
        Enthalpy(T,p) = H;
        Entropy(T,p) = S;
    end
end

[Rik]

You are assigning a large matrix (G) to a single position in Gibbs. Does the loop below yield the results you need?

for T = 153:1:423
  for p = 100000:1:10000000
      [rho(T,p),Z(T,p),phi(T,p),A(T,p),S(T,p),H(T,p),U(T,p),G(T,p)] ...
      = prsrk(z,p,T,pc,Tc,w,k,cpig,DHf,DGf,ph,flag)
      % 
      Gibbs(T,p) = G(T,p);
      IntEng(T,p) = U(T,p);
      Enthalpy(T,p) = H(T,p);
      Entropy(T,p) = S(T,p);
  end
end

If so, I'll suggest another modification and move this to the answer.

Answer 1

[Rik]

The solution Sonam Gupta offered is one way to adapt your code for indexing, but it does not support other step sizes than 1. The code below was what I hinted at in the comments.

T_list=153:1:423;%can be any sorted or unsorted list
p_list=100000:1:10000000;
for T_index = 1:numel(T_list)
    T=T_list(T_index)
    for p_index = 1:numel(p_list)
        p=p_list(p_index);
        [rho,Z,phi,A,S,H,U,G] = prsrk(z,p,T,pc,Tc,w,k,cpig,DHf,DGf,ph,flag);
        Gibbs(T_index,p_index) = G;
        IntEng(T_index,p_index) = U;
        Enthalpy(T_index,p_index) = H;
        Entropy(T_index,p_index) = S;
    end
end

Edit: replaced length by numel

I really recommend looking into that prsrk script to see if you can make it accept matrices. Removing for-loops usually saves massive amounts of calculation time.

[Guillaume]

Replace length with numel and the code will work for any shape of T_list and p_list.

[Rik]

Thanks for the suggestion, I edited my answer accordingly.

[Jerome Campbell]

I appreciate the fix. It is now storing the numbers. It does take a while, as you noted. I will try to adapt the function so matrices could be accepted to speed it up. Thank you for your help.

Answer 2

[Sonam Gupta]

Using the indices as T and p directly for storing IntEng, Gibbs, Enthalpy, Entropy array is what is creating a large array of zeros. Changing the indices such that they begin from 1 should work. I suggest you to make changes to your code as below:

 p_offset = 99999;   % 100000 - 1, so that after subtracting the offset, your index begins from 1
 T_offset = 152;     % 153 - 1 , similar thing for T also   
for T = 153:1:423
    for p = 100000:1:10000000
        [rho,Z,phi,A,S,H,U,G] = prsrk(z,p,T,pc,Tc,w,k,cpig,DHf,DGf,ph,flag);
        Gibbs(T - T_offset,p - p_offset) = G;
        IntEng(T - T_offset,p - p_offset) = U;
        Enthalpy(T - T_offset,p - p_offset) = H;
        Entropy(T - T_offset,p - p_offset) = S;
    end
end

I hope this helps in resolving the issue.

[Jerome Campbell]

Thank you for the help. This did allow me to store numbers. However, the problem I found was I couldn't take a step other than 1 for T. But it did allow me to store numbers. Thank you.