Why do I get "complex" arrays?

I am not supposed to get complex values but unfortunatley something is wrong in my eigenvalue caclultions. When plotting I get the error: (which I shouldn't be getting technically)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

The part where the complex value show is in this specific line. However, everything above is a +ve number or value so I don't see how it turns complex:

   r=q(1,:); % this is +ve when printed out
   u=q(2,:)./r; E=q(3,:)./r;
   Bx = q(4,:)./r; By = q(5,:)./r; 
   p=(gamma-1)*r.*(E-0.5*(u.^2)-(Bx.^2+By.^2)/(2*mu0)); % this is +ve
   
    a=sqrt(gamma*p./r); % this is COMPLEX ???
    
    CA1 = sqrt((Bx.^2+By.^2)./(r*mu0)); 
    CAx1 = sqrt((Bx.^2)./(r*mu0)); 
    a_MHD=sqrt(0.5*(CA1.^2+a.^2+sqrt((CA1.^2+a.^2).^2-4*a.^2.*CAx1.^2))); 
    

I posted the entire code as well, main code is "MUSCLEULER_Test"

Answer 1

[Steven Lord]

If r contains only positive values and p contains only positive values, what about gamma?

areALLElementsOfRPositive = all(r > 0)
areALLElementsOfPNonnegative = all(p >= 0)
areALLElementsOfGammaNonnegative = all(gamma >= 0)

If any elements of r, p, or gamma would result in a complex value in the corresponding element of a, the array a will have the complex attribute.

[Jamie Al]

gamma is a constant = 1.4 so also a +ve value

[Jamie Al]

I see when I inspect every element of r there are only 2 values are -ve that I am not sure where they come from. I tried using something like:

areALLElementsOfRPositive = all(r > 0)

But it returns only 0

How can I change that?

[Steven Lord]

find where the elements of r that are less than or equal to 0 are.

Use breakpoints to determine where those elements in q(1, :) are set to a negative value. You might want to add assert statements to error as soon as any element becomes negative to help you locate the introduction of the negative numbers.

Work your way back to determine why those elements are negative.

[Jamie Al]

Thanks!