Cumtrapz with cell arrays in a double loop

clear 

%g = 1; %Coupling strength
k = 20; %Decay linewidth for cavity
no = 0; %input photon number for cavity
nm = 40; %input phonon number for mechanical oscillator
Del = -10^2; %Delta, laser detuning
Ohm = 10^2; %Omega, natural frequency of mechanical oscillator
gam = 0.0032; %Decay linewidth for mechanical oscillator

w = (-1199:3:1200);
g = (0:1.25:999);

xat = 1 - (1i.*(w + Del).*(2./k));
xbt = 1 - (1i.*(w - Ohm).*(2./gam));

ca_c = cell(length(g),length(w));
da_c = cell(length(g),length(w));
Sbb_c = cell(length(g),length(w));

for u=1:length(g)
    for v=1:length(w)
        C = sqrt((4.*(g(u).^2))./(k.*gam)); %Cooperativity
        Ca(u) = C;

          c = (2.*1i.*Ca(u))./(sqrt(gam).*(xat(v).*xbt(v)+(Ca(u).^2))); %power spectrum coefficient
          ca_c{u,v}= c; %c-cell

          d =  (2.*xat(v))./(sqrt(gam).*(xat(v).*xbt(v)+(Ca(u).^2))); %power spectrum coefficient
          da_c{u,v} = d; %d-cell

          Sbb = (abs(ca_c{u,v}).^2).*(no + 0.5)+(abs(da_c{u,v}).^2).*(nm + 0.5);
          Sbb_c{u,v} = Sbb; %Sbb-cell

          A = cumtrapz(w,Sbb_c{u,v})
    end   
end

The code above wouldn't compile because of the cumtrapz function. The error says ORDER contains an invalid permutation index. I've ask something somewhat similar on MathWorks a while ago but it was not answered. https://www.mathworks.com/matlabcentral/answers/350535-on-double-looping-and-confusion-with-cell-arrays-and-cumtrapz-function

Maybe my explanation was too detailed that people didn't want to get into the details. In this case I am simply trying to find the area under the curve using cumtrapz for Sbb against w. cumptrapz has to stay in the loop in order for me to obtain an array of areas for different g values. Much help would be appreciated. Thank you.

Answer 1

[dpb]

Well, I've no idea what you're really trying to do, either, but setting a breakpoint at the line containing cumtrapz we find the arguments to it are:

K>> whos w
  Name      Size             Bytes  Class     Attributes

    w         1x800             6400  double              

K>> whos Sbb
  Name      Size            Bytes  Class     Attributes

    Sbb       1x1                 8  double              

K>> Sbb_c{u,v}
ans =
   7.6805e-08
K>>

I'll admit the error message is more than a little confusing as to what it means, but the problem is you're trying to integrate 800 points in the X vector with a single element for the corresponding Y.

Not sure why all the huge cell arrays; looks like you could simply compute an array of Y of the length of X (w) to integrate.

But, that's why it crashes; there's a big mismatch in the two arrays that must have 1:1 consonance in size in order to work...

ADDENDUM

Try

...  
Lg=length(g); Lw=length(w);
Sbb = zeros(Lg,Lw);
A=Sbb;
C=sqrt((4*(g.^2))/(k*gam));
for u=1:Lg
  for v=1:Lw
    Ca=(2*1i*C(u))/(sqrt(gam)*(xat(v)*xbt(v)+C(u)^2));
    Da=(2*xat(v))/(sqrt(gam)*(xat(v)*xbt(v)+(C(u)^2)));
    Sbb(u,v) = (abs(Ca)^2)*(no + 0.5)+(abs(Da)^2)*(nm + 0.5);
  end   
  A(u,:) = cumtrapz(w,Sbb(u,:));
end

that eliminates the cell arrays and the unneeded storage for the intermediaries--I think it is what you were looking for. Note your C is independent of v so doesn't need to be in the innermost loop on v and I went ahead and moved it out and precomputed the vector, computing the two coefficients dependent upon it on the fly.

Then the integral can be computed on each vector at the end of the inner loop.

Unless there's other need for it later, there's no reason to save the full 2D array, either; simply the 1D vector.

The inner loops could also be vectorized and the Sbb array computed without any loops at all via expanding the two x vectors to a grid via meshgrid.

I believe this addresses the comments/questions raised in the various responses to me and Walter.

Walter, thnx for the digging on the error from cumtrapz -- I was somewhat confused; didn't realize it would take a 2nd-position single value as the DIM argument without the [] placeholders...not one of the more prevalent functions I've called over the years! :)

ADDENDUM 2:

It's not that hard to vectorize over v altho following is air code; check thoroughly...

sqrtgam=sqrt(gam);   % take constant outside loops at initialization section
xabt=xat.*xbt;       % ditto--this is a constant vector, too...
CSq=C(u).^2;         % and so is this...

...
for u=1:Lg
  Ca=2*1i*C(u)./(sqrtgam*(xabt+Csq(u)));
  Da=2*xat./(sqrtgam*(xabt+Csq(u)));
  Sbb(u,:) = (abs(Ca).^2)*(no + 0.5)+(abs(Da).^2)*(nm + 0.5);
  A(u,:) = cumtrapz(w,Sbb(u,:));
end

[Alvin]

Hi, first of all thanks for the response.

What do you mean by computing an array of Y of length X(w) to integrate? You seem to notice something that I'm not noticing and it may just help trivialize the problem. Would you mind clarifying ?

Edit: If you look at Sbb_c, it is a 800x800 cell array with rows representing g and columns representing w, i need to find the area for each g value (each row in Sbb_c) and I'm using cumtrapz to do it. Hope that clarifies things.

Edit2:

I added

   Y = Sbb_c(1,:);

outside the loop and did

figure(1)
plot(w,Y)

And it comes up with an error saying invalid second data argument ?! Why is it not able to plot when Y is a 1x800 array and w is also a 1x800 array? This is truly baffling and depressing

Edit3: I fixed it, seems like Y is still a cell-array and plot doesn't like that so I converted it to a regular array using cell2mat() function. My question still stands: How should I loop this in such a way that cumtrapz takes each row in Sbb_c and compute the area with respect to w ?

[Walter Roberson]

You have

cumtrapz(w,Sbb_c{u,v})

but your Sbb_c{u,v} is a scalar.

When you call cumptraz() with two arguments, cumtrapz() looks at the size of the second argument. In the case where it is a scalar, it assumes that it is the dimension number. So you are effectively calling

cumtrapz(w, [], Sbb_c{u,v})

but 7.6805e-08 is not a valid dimension number.

[Alvin]

Thank you for the response! I think I have remedied that, my current issue is with the looping, and it is of the following:

Y = Sbb_c(100,:);
Z = cell2mat(Y);
A = cumtrapz(w,Z);

I added the three lines above and being outside the double loops and

plot(w,Z)

clearly works. What I'm doing here is picking the 100th row out of my 800x800 Sbb_c cell array and doing a cumtrapz() of that row with w. My wish is to loop it so that it includes all 800 rows (not just the 100th).