If statement with incorrect output

[JamieAl]

I have an if statement that works only for one case. For example I have a 0 and 1 case and the 1 case does not return any error but doesn't return the correct output. I am converting this function from MATLAB to c++.

My if statement looks like:

double dx = Lx/nx; //test
    double dy = Ly/ny; //test 
    for(int i = 0; i < nx; i++){
        for(int j = 0; j < nyk; j++){
            if ( FourierMeshType == 0.){    
                ksqu[j + nyk*i] = KX[j + nyk*i] * KX[j + nyk*i] + KY[j + nyk*i] * KY[j + nyk*i];                        
                ninvksqu[j + nyk*i] = -1 / (KX[j + nyk*i] * KX[j + nyk*i] + KY[j + nyk*i] * KY[j + nyk*i]);
            }
            if (FourierMeshType == 1.){
                //dx = Lx/nx; // calculate dx 
                //dy = Ly/ny; // calculates dy
                ksqu[j + nyk*i] = pow((sin(KX[j + nyk*i] * dx/2))/dx/2,2) + pow((sin(KY[j + nyk*i] * dy/2))/dy/2,2) ; // Use approximations for Kx, Ky, K^2
                // pow((sin(KX[j + nyk*i] * dx/2))/dx/2,2) + pow((sin(KY[j + nyk*i] * dy/2))/dy/2,2);
                ninvksqu[j + nyk*i] = -1 /(pow((sin(KX[j + nyk*i] * dx/2))/dx/2,2) + pow((sin(KY[j + nyk*i] * dy/2))/dy/2,2));
                // 1 /(pow((sin(KX[j + nyk*i] * dx/2))/dx/2,2) + pow((sin(KY[j + nyk*i] * dy/2))/dy/2,2))
                
                
            }
        }   
    }

This if statement is from my MATLAB code:

if meshType == 0
    % Calculate ksquared
    ksqu = KX.^2 + KY.^2;
end

% If meshType~= 1, we are using the exact k
% If meshType==1, we are using the numerical approximation
if meshType==1    
    
    % Need to calculate dx
    dx = L./n;    
    
    % Use approximations for kx, ky, and k^2. These come from Birdsall and Langdon. Find page number and put it here at some point,
    ksqu = (sin( KX * dx(1)/2)/(dx(1)/2)).^2 + (sin( KY * dx(2)/2)/(dx(2)/2)).^2 ;
    
    KXvec = sin(kx * dx(1)) / dx(1);
    KYvec = sin(ky * dx(2)) / dx(2);    
    
    [KX, KY] = meshgrid(KXvec, KYvec);
    
    % Take the tranpose since we want the first index to be for x and the second index to be for y
    KX = KX';
    KY = KY';
end

I will include the c++ code as well as an attachment for you to take a look at. What I believe the problem is in my syntax or the way I am writing the statements for the case == 1. I am still learning pretty basic stuff on c++ so my apologies for the confusion. Thanks.

[Ganado]

Specifically concerning your 0/1 switch:
What type is FourierMeshType? It appears you are comparing it against a double. This should generally be avoided. Make it be an int (or enum) and compare it against ints. You can also use else-if instead of just if (in both the C++ and MATLAB).

int FourierMeshType = 1;
if (FourierMeshType == 0)
{

}
else if (FourierMeshType == 1)
{

}

But I don't think the if statements are the problem here. Your 'ninvksqu' variable only exists in your C++ code and not in the matlab code, so it's hard to see the 1:1 mapping between them. What specifically is the "correct result" vs. actual result that you are calculating?

Edit:
Your C++ code:
pow((sin(KX[j + nyk*i] * dx/2))/dx/2,2)

Your MATLAB code:
(sin( KX * dx(1)/2)/(dx(1)/2)).^2

These are not equivalent.
x/y/z != x/(y/z)

If a line becomes too complicated and hard to read, I suggest breaking into parts to prevent mistakes like these.

PS: In general, you may want to remove if-statements within loops. Sometimes, for simple cases, the compiler can detect it and prevent unnecessary branching, but not always. This is unique to your C++ code, since your MATLAB code avoids the user-facing for loops through use of vector-based operations.

[JamieAl]

Hey, Ganado!

Yes, unfortunately I was using double for my FourierMeshType, which I just changed. Also, in the MATLAB Code "ninvksqu" is:

% When finding the inverse of ksqu for the potential equation, we need to get rid of the ksqu=0 values.
% We can do this because we don't care about any sort of constant potential.
% First, set equal to ksqu
ksqu_4inv = ksqu;

% Replace all 0 values with 1
ksqu_4inv(ksqu_4inv==0) = 1;

% Take the inverse
inv_ksqu = 1 ./ ksqu_4inv; 

Sorry about that!

These are not equivalent.
x/y/z != x/(y/z)

oh I see. Okay, would it be equivalent this way?

pow(sin(KX[j + nyk*i] * dx/2)/dx/2,2)

Thanks for the help!

[Ganado]

pow(sin(KX[j + nyk*i] * dx/2)/dx/2, 2)
is equivalent to what you already had; you simply removed a redundant parentheses set.

pow((sin(KX[j + nyk*i] * dx/2))/dx/2, 2)
pow( sin(KX[j + nyk*i] * dx/2) /dx/2, 2)

The problem that I see is with your division.
Look at the parentheses in your matlab code during the division part, and how you don't have these parentheses in the C++ code.
<...>/(dx/2)
<...>/dx/2

[JamieAl]

Ohhhh I SEE!! Yes, you're absolutely right. I have no idea how I could not see that. Thanks!!