Complex number assignment

Question

Here is my MATLAB code which calculates the fresnel diffraction.And i wanna implement the same code in C++.

clc;clear all;close all;

N=512;
M=512;
lambda=632e-9;
X = 12.1e-6;
k=2*pi/lambda;
z0 = (N*X^2)/lambda;
a = (exp(j*k*z0)/(j*lambda))
hX=(exp(j*k*z0)/(j*lambda))*exp(j*pi/(z0*lambda)*X^2*[0:N-1].^2);
hY=(exp(j*k*z0)/(j*lambda))*exp(j*pi/(z0*lambda)*X^2*[0:M-1].^2);
h = hX.'*hY;
figure; imshow(real(h), []);

And while im trying to implement the same code in C++ I tried this:

int main (){
std::complex<double> com_two(0,1);
double mycomplex = 1;
double pi = 3.1415926535897;
double N = 512;
double M = 512;
double lambda = 632e-9;
double X = 12.1e-6;
double k = (2*pi)/lambda;
double z0=(N*pow(X,2))/lambda;
//std::complex<double> hx; // Both definitions don't work
//double hy;
/*for (int j=0; j < N; j++){
    hy = (exp(mycomplex*k*z0) / (mycomplex*lambda))*exp(mycomplex*pi/(z0*lambda)*pow(X,2)*pow(j,2));
    cout << hy <<endl;
}*/
system("pause");
return 0;

}

But the thing is that while calculation performing hx and hy values return complex values.

So How should i define the "mycomplex" like i,j as in the MATLAB code. Also I found some asssignments like std::complex<double> complex;

But i guess that won't work.In my opinion i have to get hx and hy values as complex numbers.But I couldn't find a proper solution.

I'll be very glad if someone helps about that.

Answer 1

the problem is not complex number itself, but calculations that you perform. First you do different things in Matlab alg (the version that you showed at least) and in C++. Here is how you should port this code to C++.

#include <cstdlib>
#include <iostream>
#include <complex>
using namespace std;
/*
 * 
 */
int main(int argc, char** argv) {
    double pi = 3.1415926535897;
    double N = 512;
    double M = 512;
    double lambda = 632e-9;
    double X = 12.1e-6;
    double k = (2 * pi) / lambda;
    double z0 = (N * pow(X, 2)) / lambda;
    std::complex<double> hx, hy; // Both definitions don't work
    //double hy;
    int j = 1;
    //for (int j=0; j < N; j++){
        hy = (exp(j*k*z0) / (j*lambda))*exp(j*pi/(z0*lambda)*pow(X,2)*pow(j,2));
        cout << hy <<endl;
    //}
    system("pause");
    return 0;
}

now the problem is in your lambda = 632e-9 which is close to zero, then k is very big k = (2 * pi) / lambda; as the inverse of something which approaches 0, and in turn

exp(j*k*z0)

is very big, and this

(exp(j*k*z0) / (j*lambda))

even more. Then this

(exp(j*k*z0) / (j*lambda))*exp(j*pi/(z0*lambda) 

is huge, and the remaining factor pow(X,2)*pow(j,2)) is meaningful. So then

cout << hy <<endl; 

will print (inf,0) even for lambda = 0.0001; (how about 632e-9 !).

However you can still print it if you use small enough j value to absorb huge k value in exp (j is also in denominator but de l'Hopital rule will assure this will be smaller).

As follows from your comments (not your code) you wanted j to be complex number z=1i. Then the solution is to change

double j =1;

to

std::complex<double> j(0,1); 

---

SOLUTION:

int main(int argc, char** argv) {
    double pi = 3.1415926535897;
    double N = 512;
    double M = 512;
    double lambda = 0.00001;//632e-9;
    double X = 12.1e-6;
    double k = (2 * pi) / lambda;
    double z0 = (N * pow(X, 2)) / lambda;
    std::complex<double> hx, hy;
    std::complex<double> j(0,1);
        hy = (exp(j*k*z0) / (j*lambda))*exp(j*pi/(z0*lambda)*pow(X,2)*pow(j,2));
        cout << hy <<endl;
    return 0;
}

Answer 2

I have changed a few lines in your code and it gives results different from 1#INF

1) Change mycomplex from double to std::complex<double> and assign (0,1) to make it imaginary 1. Actually com_two has the same type and value you need for mycomplex and you are not using it.

2) Remove comments. Change hy type from double to std::complex<double>. Remove unused hx declaration.

3) My compiler finds an error at pow(j,2) because of ambiguity, so i fixed it too.

#include <iostream>
#include <complex>

using namespace std;

int main (){
  complex<double> mycomplex(0,1);
  double pi = 3.1415926535897;
  double N = 512;
  double M = 512;
  double lambda = 632e-9;
  double X = 12.1e-6;
  double k = (2*pi)/lambda;
  double z0=(N*pow(X,2))/lambda;
  complex<double> hy;
  for (int j=0; j < N; j++){
      hy = (exp(mycomplex*k*z0) / (mycomplex*lambda))*exp(mycomplex*pi/(z0*lambda)*pow(X,2)*pow((double)j,2));
      cout << hy <<endl;
  }
  system("pause");
  return 0;
}