如何在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;
}