Nested loops and function output into a matrix inconsistency

I am attempting to fit a set of data simultaneously with two varying inputs: the temperature and h-constant. Both of these inputs are in single row, 400-element vector. The fundamental piece of my code is the function B that takes two inputs (scalars x and h) and outputs a column vector of 43 elements -- the 43 elements are based off of the "nu" vector (frequency) input. From my understanding, to achieve the best fit of the .txt file I implemented, I have to consider the initial temperature evaluated with ALL 400 h-constants and repeat that for the other 399 temperatures. As such, it seems like I need 160,000 outputted column vectors. So that's why I have a 43x160000 matrix.

The plan is to take each temperature and cycle through each h-constant before continuing on with the next T(i) value. Each column of my "md" matrix will be the output of the B function. The problems I'm encountering are that: the columns repeat themselves and that the output in the matrix doesn't match the numbers (example like B(T(1),H(1)) in the command window --> this should be the first column in "md" but it's not). I think the crux of the problem lies in the loop statement. Where did I misstep? How can I ensure the output replaces each column in the large matrix, accounting for every single T and h? I've attached the .txt file if the need arises. Thank you.

Dprime=readtable('FIRASdata.txt'); 
spints=Dprime.Var2; invwav=Dprime.Var1; % extracting the two columns
c=2.9979*10^(8); k=1.3806504*10^(-23) ; J=2/c^2;% constants
nu=c./(invwav.^(-1)*0.01); % CONVERSION; frequency of FIRAS data Hz^-1
T=linspace(1,5,400); % range of temperatures
H=linspace(2*10^(-34),1*10^(-33),400); %range of h-constants
B= @(x,h) 10^(20)*J*h*(nu.^3./(exp((h*nu)./(k*x))-1)); % anonymous function to return MJy/sr intensities
%% 
md=zeros(43,160000);
for q=1:400
    for r=1:400
        for w=1:160000
    md(:,w)=B(T(q),H(r));
        end
    end
end

Answer 1

[Sindar]

It seems like your 'w' columns are supposed to correspond to q-r pairs. Currently, this isn't happening, so each q-r pair will fill every column, resulting in the last one overwriting all the others. One way around this is to directly connect w to q-r. sub2ind is a convenient way to do this, since you can reverse it (extract q-r from w) with ind2sub:

qN = 400;
rN = 400;
for q=1:qN
    for r=1:rN
        w = sub2ind([rN qN],r,q);
        md(:,w)=B(T(q),H(r));
    end
end

Alternatively, you could create a 3-dimensional matrix:

qN = 400;
rN = 400;
for q=1:qN
    for r=1:rN
        md(:,r,q)=B(T(q),H(r));
    end
end

[Sean Rapp]

Thank you, Sindar. The sub2ind method is one I ought to implement more. The 3-dimensional matrix technique is apparent to me now that I see it executed as well. I appreciate the help.