MoM矩量法（三）：激励矩阵以及方程求解

MoM矩量法（三）：激励矩阵以及方程求解

填充完阻抗矩阵之后就是激励矩阵的填充以及方程的求解了。激励矩阵比较简单，后面矩阵方程的求解也可以采用库函数来完成。
这里采用的公式来自于聂在平主编的《目标与环境 电磁散射特性建模——理论、方法与实现（应用篇）》，所有数据均采用国际标准单位。

完整代码已贴在文章最后。文章和代码思路借鉴了https://blog.csdn.net/u014411646/article/details/99689460，谢谢原作者！

一、激励矩阵

上式中各符号的含义：

1. $E_{\theta }\theta$和$E_{\varphi }\varphi$代表$\theta$方向极化电场和$\varphi$方向极化电场，在这里我们只考虑二者其一，而不考虑其他极化方向的入射波。所以，当确定了极化方向时，上式中的$E_{\theta }\theta +E_{\varphi }\varphi$只需要保留一个即可；
2. $k_i$为入射方向单位向量；
3. 其他变量与MoM矩量法（二）：阻抗矩阵的填充定义的变量一致。

核心代码：

cx_vec setVm(rowvec ei) {
	cx_vec vm(EdgeTotal);
	for (int i = 0;i < EdgeTotal;i++) {
		Complex temp = 0;
		for (int k = 0;k < 4;k++) {
			rowvec rlmk(3);
			rlmk(0) = xx[k] * pdata(Edge_list(i, 2), 0) + yy[k] * pdata(Edge_list(i, 0), 0) + zz[k] * pdata(Edge_list(i, 1), 0);
			rlmk(1) = xx[k] * pdata(Edge_list(i, 2), 1) + yy[k] * pdata(Edge_list(i, 0), 1) + zz[k] * pdata(Edge_list(i, 1), 1);
			rlmk(2) = xx[k] * pdata(Edge_list(i, 2), 2) + yy[k] * pdata(Edge_list(i, 0), 2) + zz[k] * pdata(Edge_list(i, 1), 2);
			rowvec rrmk(3);
			rrmk(0) = xx[k] * pdata(Edge_list(i, 3), 0) + yy[k] * pdata(Edge_list(i, 1), 0) + zz[k] * pdata(Edge_list(i, 0), 0);
			rrmk(1) = xx[k] * pdata(Edge_list(i, 3), 1) + yy[k] * pdata(Edge_list(i, 1), 1) + zz[k] * pdata(Edge_list(i, 0), 1);
			rrmk(2) = xx[k] * pdata(Edge_list(i, 3), 2) + yy[k] * pdata(Edge_list(i, 1), 2) + zz[k] * pdata(Edge_list(i, 0), 2);
			temp += ww[k] * (exp(-J * k0 * getDot(ki, rlmk)) * getDot(ei, rlmk - pdata.row(Edge_list(i, 2))) -
				exp(-J * k0 * getDot(ki, rrmk)) * getDot(ei, rrmk - pdata.row(Edge_list(i, 3))));
		}
		vm(i) = EdgeLength(i) * temp;
	}
	return vm;
}

与前一节描述的阻抗矩阵类似，通过对数据运算将其填充进激励矩阵。之后就是对矩阵方程的求解，从而得到表面电流。

核心代码：

cx_vec solveMatrix(cx_mat zmn, cx_vec v) {
	int i, j;
	int ii = -1;
	Complex sum;
	for (i = 0;i < EdgeTotal;i++) {
		sum = v(i);
		if (ii > -1) {
			for (j = ii;j <= i - 1;j++) {
				sum -= zmn(i, j) * v(j);
			}
		}
		else if (sum != 0.) {
			ii = i;
		}
		v(i) = sum;
	}
	for (i = EdgeTotal - 1; i >= 0;i--) {
		sum = v(i);
		for (j = i + 1;j < EdgeTotal;j++) {
			sum -= zmn(i, j) * v(j);
		}
		v(i) = sum / zmn(i, i);
	}
	return v;
}

需要注意的是，在上述解矩阵方程的函数中，参数$Zmn$需要是上一节经过LU分解的矩阵。

下面贴上所有代码：

/***********************************************************/
/*              该函数对激励矩阵进行填充                    */
/*                   并求解方程组                          */
/***********************************************************/
/*                                                         */
/*              Author：Tianshaowenwe                      */
/*                                                         */
/***********************************************************/
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <complex>
#include <math.h>
#include <time.h>
#include <windows.h>
#include <cmath>
#include "rwg3.h"
#include "armadillo"

using namespace arma;
using namespace std;

const Complex J = Complex(0, 1);
extern mat Edge_list;//[端点1， 端点2， 对点1， 对点2， 左三角单元， 右三角单元]
extern mat TriCenter;
extern vec EdgeLength;
extern vec TriArea;
extern int TriTotal, EdgeTotal;
extern mat pdata, tdata;
extern double k0;
extern double thi, phi;
extern cx_mat Zmn;
extern int pol;
extern double xx[4];
extern double yy[4];
extern double zz[4];
extern double ww[4];
extern double xxx[7];
extern double yyy[7];
extern double zzz[7];
extern double www[7];

cx_vec I, V;
rowvec ei(3), ki(3);
void rwg3_main() {
	rwg2_main();
	ki(0) = -1. * sin(thi) * cos(phi);
	ki(1) = -1. * sin(thi) * sin(phi);
	ki(2) = -1. * cos(phi);
	if (pol == 0) {
		ei(0) = cos(thi) * cos(phi);
		ei(1) = cos(thi) * sin(phi);
		ei(2) = -1. * sin(thi);
	}
	else {
		ei(0) = -1. * sin(phi);
		ei(1) = cos(phi);
		ei(2) = 0;
	}
	cout << "\n" << "填充激励矩阵..." << endl;
	clock_t start = clock();
	V = setVm(ei);
	clock_t end = clock();
	cout << "\n" << "激励矩阵填充时间： " << (end - start) / 1000.0 << "s" << endl;
	/*
	V.save("Vph.txt", raw_ascii);
	*/
	cout << "\n" << "求解矩阵方程..." << endl;
	start = clock();
	I = solveMatrix(Zmn, V);
	end = clock();
	cout << "\n" << "矩阵方程求解时间： " << (end - start) / 1000.0 << "s" << endl;
	/*
	I.save("Ith.txt", raw_ascii);
	*/
}
/***********************************************************/
/*            Function：cx_vec setVm()                     */
/***********************************************************/
/* 1.填充激励矩阵                                           */
/***********************************************************/
cx_vec setVm(rowvec ei) {
	cx_vec vm(EdgeTotal);
	for (int i = 0;i < EdgeTotal;i++) {
		Complex temp = 0;
		for (int k = 0;k < 4;k++) {
			rowvec rlmk(3);
			rlmk(0) = xx[k] * pdata(Edge_list(i, 2), 0) + yy[k] * pdata(Edge_list(i, 0), 0) + zz[k] * pdata(Edge_list(i, 1), 0);
			rlmk(1) = xx[k] * pdata(Edge_list(i, 2), 1) + yy[k] * pdata(Edge_list(i, 0), 1) + zz[k] * pdata(Edge_list(i, 1), 1);
			rlmk(2) = xx[k] * pdata(Edge_list(i, 2), 2) + yy[k] * pdata(Edge_list(i, 0), 2) + zz[k] * pdata(Edge_list(i, 1), 2);
			rowvec rrmk(3);
			rrmk(0) = xx[k] * pdata(Edge_list(i, 3), 0) + yy[k] * pdata(Edge_list(i, 1), 0) + zz[k] * pdata(Edge_list(i, 0), 0);
			rrmk(1) = xx[k] * pdata(Edge_list(i, 3), 1) + yy[k] * pdata(Edge_list(i, 1), 1) + zz[k] * pdata(Edge_list(i, 0), 1);
			rrmk(2) = xx[k] * pdata(Edge_list(i, 3), 2) + yy[k] * pdata(Edge_list(i, 1), 2) + zz[k] * pdata(Edge_list(i, 0), 2);
			temp += ww[k] * (exp(-J * k0 * getDot(ki, rlmk)) * getDot(ei, rlmk - pdata.row(Edge_list(i, 2))) -
				exp(-J * k0 * getDot(ki, rrmk)) * getDot(ei, rrmk - pdata.row(Edge_list(i, 3))));
		}
		vm(i) = EdgeLength(i) * temp;
	}
	return vm;
}
/***********************************************************/
/*   Function：cx_vec solveMatrix(cx_mat zmn, cx_vec v)    */
/***********************************************************/
/* 1.求解矩阵方程                                          */
/***********************************************************/
cx_vec solveMatrix(cx_mat zmn, cx_vec v) {
	int i, j;
	int ii = -1;
	Complex sum;
	for (i = 0;i < EdgeTotal;i++) {
		sum = v(i);
		if (ii > -1) {
			for (j = ii;j <= i - 1;j++) {
				sum -= zmn(i, j) * v(j);
			}
		}
		else if (sum != 0.) {
			ii = i;
		}
		v(i) = sum;
	}
	for (i = EdgeTotal - 1; i >= 0;i--) {
		sum = v(i);
		for (j = i + 1;j < EdgeTotal;j++) {
			sum -= zmn(i, j) * v(j);
		}
		v(i) = sum / zmn(i, i);
	}
	return v;
}