On double looping and confusion with cell arrays and cumtrapz() function

clear 

g = 1000; 
k = 20; 
no = 0; 
nm = 40; 
Del = -(10.^2); 
Ohm = (10.^2); 
gam = 0.0032; 

C = sqrt((4.*(g.^2))./(k.*gam)); 

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

c =  (2.*1i.*C)./(sqrt(gam).*(xat.*xbt+(C.^2))); %power spectrum coefficient
d =  (2.*xat)./(sqrt(gam).*(xat.*xbt+(C.^2))); %power spectrum coefficient

Sbb = (abs(c).^2).*(no+(1/2))+(abs(d).^2).*(nm+(1/2)); %power spectrum

A = cumtrapz(w,Sbb);
LA = A(length(A));

figure(1)
plot(w,Sbb);

Should one run the code above, one will find a plot with two sharp peaks. At g = 0, there will be one peak; With increasing g, there will be a gradual splitting of the peak into the eventual two as shown in the setting above (g = 100). Moreover for each g setting, I wish to find the area under the curve of Sbb versus w (hence my use of the cumtrapz() function). And given how cumtrapz() works, I made LA just to obtain the last entry of my A (area) array (which would give the correct area value) for some g value (g = 1000 in the above case). As it stands, LA has only one value given a particular g value, whereas A is an array since it is dependent on w and w is an array.

My question is of the following: I wish to plot a graph of LA vs C (C is dependent on g) with C varying over several values (and hence g varies as well). So, clearly I will have to loop over g and make an array to store all my C values. But since c and d depends on xat, xbt, and C, they will have to be included in the loop as well since my C is an array (which is fine). The difficulty comes in introducing w into the loop (so a double loop with g and w). In short, I ended up having to make a cell array to store all the values but became confused over the long run. I would appreciate any help in resolving this issue (for all I know I might be overthinking it).

Thank you very much in advance!

[Alvin]

As a follow up, I attempted a solution below.

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));

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

    c =  (2.*1i.*Ca(u))./(sqrt(gam).*(xat(u).*xbt(u)+(Ca(u).^2))); %power spectrum coefficient
    ca(u)= c; %c-array
    d =  (2.*xat(u))./(sqrt(gam).*(xat(u).*xbt(u)+(Ca(u).^2))); %power spectrum coefficient
    da(u) = d; %d-array

    Sbb = (abs(ca(u)).^2).*(no + 0.5)+(abs(da(u)).^2).*(nm + 0.5);
    Sbba(u) = Sbb; %Sbb-array

    %Qa(u) = Q;
end

A = cumtrapz(w,Sbba); %phonon population
LA = A(length(A));
%Qa(u) = Q;

figure(1)
plot(Ca,A)
set(gca,'FontSize',13)
title('\Gamma_m = 0.2, \kappa = 1, \Omega_m = -1,\Delta = 1')
xlabel('C (cooperativity)')
ylabel('n_b (phonon population)')

The code above runs but does not make an accurate plot of LA vs Ca (I'm plotting A vs Ca afterall). Nonetheless I'm encountering issues in putting my cumtrapz function into my loop. I keep getting the error that says ORDER contains an invalid permutation index whenever I try to put my A in my loop. Moreover I feel like I'm overcomplicating this with all kinds of arrays. Please help!

Answer 1

[Saurabh Gupta]

Most of your original code is element-wise operation on vector/matrix elements, so I don't see any reason to use a loop to handle 'g' as a vector. If you set 'g' as follows, your original code should be able to evaluate all values up to the 'cumprapz' statement.

>> g = (0:1.25:999)';

Now, when you reach 'cumprapz', you will notice an error regarding dimensions mismatch because 'w' is 1x3001 and 'Sbb' is 800x3001, whereas it requires the first dimension of 'Sbb' to match the length of 'w' and it causes the error. You can loop through the first dimension of 'Sbb' and run 'cumprapz' command 800 times similar to following code, and plot corresponding results.

A = zeros(size(Sbb));
for i = 1:size(Sbb,1)
    A(i,:) = cumtrapz(w, Sbb(i,:)); 
end

Hope this helps!