如何在MEEP中定义会聚的高斯光束

FDTD计算中，经常要用到高斯光束 ，MEEP中没有给出现成的高斯光束的定义。
按浙江大学出版的＜光电子学＞书中的定义，写了一段定义聚焦高斯光束的代码。

其中，
g_beam_width 是束腰的宽度。
g_beam_freq 是工作频率(归一化的)。
g_beam_xcen, g_beam_ycen 是束腰的中心位置。
调用的方式：
continuous_src_time src(g_beam_freq, 20);

geometric_volume twist(vec(g_pml, 0), vec(g_pml, g_y));
在y方向，源伸入到了pml里，这样可以保证边缘处不发生突变。

f.add_volume_source(Hz, src, twist, gaussian, 1e3);
f.add_volume_source(Ez, src, twist, gaussian, 1e3);


complex gaussian(const vec &p) {

double width = g_beam_width;
double lambda = 1./g_beam_freq;
double y = g_beam_ycen - g_y/2.;
double z = -g_beam_xcen;
double z0 = pi * width*width / lambda;
double wz = sqrt( width*width * ( 1.0 + pow(z/z0, 2) ) );
double phiz = atan(z/z0);
double r = sqrt( pow( p.x() - z, 2 ) + pow( p.y() - y, 2 ) );
double k = 2*pi/lambda;
double Rz = 0;

complex u(0.0, 0.0);

if ( abs(z) > 1e-6 ) {
Rz = z * ( 1.0 + pow(z0/z, 2) );
u = width/wz * exp( -r*r / pow(wz, 2) ) * polar(1.0, ( k * (z + r*r/(2*Rz)) - phiz) );
}
else {
Rz = 0;
u = width/wz * exp( -r*r / pow(wz, 2) );
}

return u;
}

如何判断一个点在不在三角形内部

用MEEP做计算时，初始化epsilon，想做一个三角形，这就需要判断这个点在不在三角形的内部。
写了一小段C++代码来判断一个点是不是在三角形内部。
基本思想：
这个点和三角形三个顶点构成了三个小三角形，
比较这三个小三角形的面积之和与原三角形的面积，如果相等，说明这个点落在三角形内部。
如果三个小三角形的面积之和大于原三角形的面积，则说明这个点落在三角形的外面。

技术细节：
1、已知三角形三个顶点的坐标，求三角形的面积。1/2*|axb|即是三角形的面积。
a和b分别是两条边。
2、判断大小时，要注意数值误差，也就是1e-6的作用。

double area(const vec p, const vec p2, const vec p3)
{
double x = p2.x() - p.x();
double y = p2.y() - p.y();

double x2 = p3.x() - p.x();
double y2 = p3.y() - p.y();

double s = abs(x*y2 - x2*y)/2.;
return s;
}

bool inside(const vec insert, const vec corner, const vec corner2, const vec corner3)
{
double s = area(corner, corner2, corner3);
double ss = area(insert, corner, corner2);
double ss2 = area(insert, corner2, corner3);
double ss3 = area(insert, corner, corner3);

if ( ss + ss2 + ss3 - s > 1e-6 ) return false;

return true;
}