波束在自由空间传播快速算法（平面到平面）

瑞利 - 索末菲衍射积分公式 被广泛的应用于这一问题的分析当中。但由于该公式计算过程中所涉及的巨大的计算量，瑞利-索末菲衍射积分可以理解为两个式子的卷积，可以利用先做FFT，再做IFFT来简化，以下是实现代码，用到了fftw库。

#ifndef CALCUlATION_H
#define CALCULATION_H

#include<cmath>
#include<complex>
#include "Tools/GEMS_Memory.h"
#include "Tools/Constant_Var.h"
#include"calculation/include/FFT.h"

class Calculation
{
public:
    Calculation(double f = 10.65e9, double theta = 0, double theta_h0 = 0, double z0 = 1,int N = 361, int M = 361);
    ~Calculation();

    void Initialization();
    void FreeCal();
    void Setds(double ds1); //设置ds

    //设置输入并补0
    void SetInput(std::complex<double> ** ExIn, std::complex<double> ** EyIn);

    void StartCal();

    //输出Ex Ey Ez Hx Hy Hz 并调用FreeCal
    void output(std::complex<double> ** Ex, std::complex<double> ** Ey, std::complex<double> ** Ez,
        std::complex<double> ** Hx, std::complex<double> ** Hy, std::complex<double> ** Hz );

private:
    //定义复数的乘
    void MulComp(const double r1, const double i1, const double r2, const double i2, double & outR, double & outI);
    std::complex<double> MulComp(std::complex<double> & In, double r2, double i2);
    //插值函数
    double InserVal(const double x0, const double yUp, const double y0, const double yDown);

    //计算传达函数
    void Calh0();

    //对输入场进行旋转
    void Rotate();

    //对输入场进行插值
    //void InserExEy();

private:
    double f; // 频率 默认 10.65e9
    double lamda; // 波长
    double k;
    //double theta;
    //double theta_h0;
    double ds;
    double z0; //目标距离

    int N, M; //计算的点数 一般设为奇数 默认360
    int N2, M2; // N2 = 2 * N -1

    //传递函数
    std::complex<double> ** HExy, ** HEz_x0, **HEz_y0;
    std::complex<double> ** HHx_x0, ** HHx_y0,
        ** HHy_x0, ** HHy_y0,
        ** HHz_x0, ** HHz_y0;

    //补0后的输入
    std::complex<double> ** Ex0, ** Ey0;

    // 帕斯维尔定理
    double     PowerHExy, PowerFFTHExy;
    double     PowerHEz_x0, PowerFFTHEz_x0;
    double     PowerHEz_y0, PowerFFTHEz_y0;
    double     PowerHHx_x0, PowerFFTHHx_x0;
    double     PowerHHx_y0, PowerFFTHHx_y0;
    double     PowerHHy_x0, PowerFFTHHy_x0;
    double     PowerHHy_y0, PowerFFTHHy_y0;
    double     PowerHHz_x0, PowerFFTHHz_x0;
    double     PowerHHz_y0, PowerFFTHHz_y0;

    //Ex1 Ey1 Ez1 Hx1 Hy1 Hz1
    std::complex<double> ** Ex1, ** Ey1, **Ez1;
    std::complex<double> ** Hx1, ** Hy1, **Hz1;
};

#include"calculation/include/Calculation.h"
#include <fstream>

Calculation::Calculation(double f, double theta, double theta_h0, double z0,int N, int M)
	:f(f),z0(z0),N(N),M(M)
{
	lamda = C_Speed / f;
	k = 2 * Pi * f / C_Speed;
	ds = lamda / 3.5;
	N2 = 2 * N - 1;
	M2 = 2 * M - 1;

	PowerHExy = 0;
	PowerHEz_x0 = 0;
	PowerHEz_y0 = 0;
	PowerHHx_x0 = 0;
	PowerHHx_y0 = 0;
	PowerHHy_x0 = 0;
	PowerHHy_y0 = 0;
	PowerHHz_x0 = 0;
	PowerHHz_y0 = 0;

	PowerFFTHExy = 0;
	PowerFFTHEz_x0 = 0;
	PowerFFTHEz_y0 = 0;
	PowerFFTHHx_x0 = 0;
	PowerFFTHHx_y0 = 0;
	PowerFFTHHy_x0 = 0;
	PowerFFTHHy_y0 = 0;
	PowerFFTHHz_x0 = 0;
	PowerFFTHHz_y0 = 0;

	Initialization();
}

Calculation::~Calculation(