2D快速多极子代码: nero2d.zip
2D快速多极子代码:
- /***************************************************************************
- * Copyright (C) 2006 by Jan Fostier *
- * jan.fostier@intec.ugent.be *
- * *
- * This program is free software; you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation; either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program; if not, write to the *
- * Free Software Foundation, Inc., *
- * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
- ***************************************************************************/
- #include "excitation.h"
- #include "scene2d.h"
- #include "mathlib.h"
- // ========================================================================
- // POINT SOURCES
- // ========================================================================
- dcomplex Excitation::pointSourceExField(const Point& r) const
- {
- dcomplex H2;
- fastHankel0(medium->gamma*(orig-r).norm(), &H2);
- return -0.25*2*pi*medium->frequency*medium->mu*I0*H2;
- }
- dcomplex Excitation::pointSourceHtField(const Point& r, const Point& t) const
- {
- Point arg = orig-r;
- double narg = arg.norm();
- dcomplex result[2];
- fastHankel01(medium->gamma*narg, result);
- dcomplex H2 = result[1];
- Point tmp(arg.getY(), -arg.getX());
- double dt = t.dot(tmp)/narg;
- return -0.25*I*I0*medium->gamma*dt*H2;
- }
- // ========================================================================
- // TM-WAVES
- // ========================================================================
- dcomplex Excitation::planeWaveTMExField(const Point& r) const
- {
- Point u(cos(angle), sin(angle));
- return Ex0*exp(-I*medium->gamma*r.dot(u));
- }
- dcomplex Excitation::planeWaveTMHtField(const Point& r, const Point& t) const
- {
- Point u(cos(angle), sin(angle));
- Point tmp(u.getY(), -u.getX());
- return Ex0/medium->Z*t.dot(tmp)*exp(-I*medium->gamma*r.dot(u));
- }
- // ========================================================================
- // GAUSSIAN BUNDLES
- // ========================================================================
- dcomplex Excitation::gaussBundleExField(const Point& r) const
- {
- // perform coordinate transformation to gauss bundle intrinsic coords
- double cphi = cos(angle);
- double sphi = sin(angle);
- Point arg(r.getX()*cphi+r.getY()*sphi, -r.getX()*sphi+r.getY()*cphi);
- arg = arg-orig;
- dcomplex RR = 0.5*medium->gamma*gaussW*gaussW;
- dcomplex q = arg.getX()+I*RR;
- return Ex0/sqrt(q)*exp(-I*medium->gamma*(0.5*arg.getY()*
- arg.getY()/q+arg.getX()));
- }
- dcomplex Excitation::gaussBundleHtField(const Point& r, const Point& t) const
- {
- // perform coordinate transformation to gauss bundle intrinsic coords
- double cphi = cos(angle);
- double sphi = sin(angle);
- Point arg(r.getX()*cphi+r.getY()*sphi, -r.getX()*sphi+r.getY()*cphi);
- arg = arg-orig;
- Point ta(t.getX()*cphi+t.getY()*sphi, -t.getX()*sphi+t.getY()*cphi);
- dcomplex RR = 0.5*medium->gamma*gaussW*gaussW;
- dcomplex q = arg.getX()+I*RR;
- dcomplex GB = Ex0/sqrt(q)*exp(-I*medium->gamma*(0.5*arg.getY()*
- arg.getY()/q+arg.getX()));
- dcomplex Htx = -I*medium->gamma*arg.getY()/q*GB;
- dcomplex Hty = -(0.5*I*medium->gamma*arg.getY()*arg.getY()/q/q
- -0.5/q-medium->gamma*I)*GB;
- return Ex0*I/2.0/pi/medium->frequency/medium->mu*
- (ta.getX()*Htx+ta.getY()*Hty);
- }
- // ========================================================================
- // PUBLIC FUNCTIONS
- // ========================================================================
- dcomplex Excitation::calcExField(const Point& r) const
- {
- switch (source) {
- case POINTSOURCE:
- return pointSourceExField(r);
- case GAUSSBUNDLE:
- return gaussBundleExField(r);
- case PLANEWAVE:
- return planeWaveTMExField(r);
- }
- }
- dcomplex Excitation::calcHtField(const Point& r, const Point& t) const
- {
- switch (source) {
- case POINTSOURCE:
- return pointSourceHtField(r, t);
- case GAUSSBUNDLE:
- return gaussBundleHtField(r, t);
- case PLANEWAVE:
- return planeWaveTMHtField(r, t);
- }
- }
- std::ostream& operator<<(std::ostream& os, const Excitation& E)
- {
- // shortcuts
- Point r = E.orig;
- dcomplex E0 = E.Ex0; dcomplex I0 = E.I0;
- double w = E.gaussW; double angle = E.angle;
- switch(E.source) {
- case POINTSOURCE:
- os << "PointSrc at (" << r << "), I = " << I0;
- break;
- case GAUSSBUNDLE:
- os << "GaussB at (" << r << "), E0 = " << E0
- << " angle = " << angle << " width = " << w;
- break;
- case PLANEWAVE:
- os << "PlaneWave E0 = " << E0 << ", angle = " << angle;
- break;
- }
- return os;
- }
支持代码分享
很好的东西，许多人都在寻找！
好东西啊，元柱辛苦了
谢谢！正在学习。
谢谢分享，支持源代码
缺文件了，楼主有完整版的程序吗？
