00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include <cmath>
00014 #include <tut.h>
00015 #include "CircleZ.h"
00016 #include "LineAtoB.h"
00017 #include "romberg.h"
00018
00019 using mesh::Vector3;
00020
00021 template <typename Curve>
00022 struct Curve_length {
00023 explicit Curve_length(Curve curve)
00024 : curve_(curve)
00025 {
00026 }
00027 double operator()(double t) {
00028 Vector3 d(curve_.deriv(t));
00029 return sqrt(d.x*d.x + d.y*d.y + d.z*d.z);
00030 }
00031 Curve curve_;
00032 };
00033
00036 template <typename Curve>
00037 class MyEfieldDs {
00038 public:
00046 MyEfieldDs(Curve curve, Vector3 E)
00047 : curve_(curve),
00048 E_(E)
00049 { }
00051 double operator()(const double t) {
00052 Vector3 field(E_);
00053 return field.dot_product(curve_.deriv(t));
00054 }
00055 private:
00057 Curve curve_;
00059 const Vector3 E_;
00060 };
00061
00062 namespace
00063 {
00064
00065 struct testdata
00066 {
00067 testdata()
00068 : pi(std::acos(-1.0)),
00069 tol(1e-10)
00070 {
00071 }
00072 const double pi;
00073 const double tol;
00074 };
00075
00076 typedef tut::test_group<testdata> testgroup;
00077 typedef testgroup::object testobject;
00078 testgroup romberg_group("line_integral");
00079
00081 template<>
00082 template<>
00083 void testobject::test<1>()
00084 {
00085 int info;
00086 const double radius = 2.0;
00087 CircleZ curve(radius, 0.0);
00088 Curve_length<CircleZ> f_len(curve);
00089 double I = NR::qromb<Curve_length<CircleZ> >(f_len, 0.0, 2.0*pi, 1e-12, info);
00090 double expected = 2.0*pi*radius;
00091 ensure(fabs(I - expected) <= tol);
00092 }
00093
00095 template<>
00096 template<>
00097 void testobject::test<2>()
00098 {
00099 int info;
00100 Vector3 a(2.3, 7.8, -1.2);
00101 Vector3 b(-2.0, 5.1, 3.7);
00102 Vector3 e(-2.5, -8.9, 4.5);
00103 LineAtoB curve(a, b);
00104 MyEfieldDs<LineAtoB> e_field_integrand(curve, e);
00105 double I = NR::qromb<MyEfieldDs<LineAtoB> >(e_field_integrand, 0.0, 1.0, 1e-12, info);
00106 Vector3 d(b-a);
00107 double Eproj = e.dot_product(d) / d.length();
00108 double expected = Eproj*d.length();
00109 ensure(fabs(I - expected) <= tol);
00110 }
00111
00112 };