Home| Inhalt| Theorie| Experimente| Versionen| PDF
| de

Source-Code

C++ Objekt zur numerischen Berechnung elektrischer Felder

Header-File

Von hier kann der Code heruntergeladen werden. Anwendungsbeispiele finden sich hier, hier und hier.
0/*
1 Dr. Steffen Kuehn / 2015 / steffen.kuehn@quantino-theory.de
2
3 This program is free software: you can redistribute it and/or modify
4 it under the terms of the GNU General Public License as published by
5 the Free Software Foundation, either version 3 of the License, or
6 (at your option) any later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program. If not, see <http://www.gnu.org/licenses/>.
15 */
16
17#ifndef CHARGE_H
18#define CHARGE_H
19
20typedef long double sim_double;
21
22class TVector
23{
24public:
25    TVector(void);
26    TVector(sim_double x, sim_double y, sim_double z);
27    TVector(const TVector& v);
28    ~TVector();
29    
30    friend TVector operator+(const TVector &v1, const TVector &v2);
31    friend TVector operator-(const TVector &v1, const TVector &v2);    
32    friend TVector operator-(const TVector &v);
33    friend sim_double operator*(const TVector &v1, const TVector &v2);
34    friend TVector operator*(sim_double v1, const TVector &v2);
35    friend TVector operator*(const TVector &v1, sim_double v2);
36    friend TVector operator/(const TVector &v1, sim_double v2);
37    
38    friend sim_double nrm(const TVector &v);
39    
40    sim_double x;
41    sim_double y;
42    sim_double z;
43};
44
45class TOptParam
46{
47public:
48    TOptParam(void);
49    
50    // buffer size (should be great for speeds near the light speed!)
51    int n;
52    // for partial interval calculation
53    sim_double re;
54    // maximal partial integration interval number
55    int lim;
56};
57
58class TCharge
59{
60public:
61    TCharge(void);
62    ~TCharge();
63    
64    // returns the electric field strength at position r for the current time.
65    // v is the receiver speed.
66    TVector Efield(const TVector& r, const TVector& v);
67    
68    // classically calculated electric field for calibration
69    TVector CoulombEfield(const TVector& r, const TVector& v);
70    
71    // has to be called after the constructor
72    // r - initial position in m
73    // v - initial speed in m/s
74    // q - charge in elementary charges
75    // m - inertial mass in kg
76    // op - calculation parameters
77    void Initialise(const TVector& r, const TVector& v,
78        int q, sim_double m, const TOptParam& op);
79    
80    // time evolution. new position and speed is given
81    void TimeStep(const TVector& r, const TVector& v, sim_double dt);
82    
83    // time evolution under influence of a force
84    void TimeStep(const TVector& efield, sim_double dt);
85    
86    // getter functions
87    TVector GetPosition(void) const {return rh[rbn];}
88    sim_double GetInertialMass(void) const {return m;}
89    TVector GetSpeed(void) const {return vh[rbn];}
90    sim_double GetCharge(void) const {return q;}
91    sim_double GetTime(void) const {if(th) return th[rbn]; else return 0;}
92    sim_double GetBufferState(void) const {return ((sim_double)rbs)/op.n;}
93    
94    const static sim_double c = 299792458;     // light speed in m/s
95    const static sim_double mu0 = 12.566370614e-7; // N/A^2
96    const static sim_double e = 1.602176565e-19; // elementary charge in C
97    
98private:    
99    // position history (ring buffer)
100    TVector* rh;
101    // speed history (ring buffer)
102    TVector* vh;
103    // time point history (ring buffer)
104    sim_double* th;
105    // position of the newest element in the ring buffer
106    int rbn;
107    // number of valid elements in the ring buffers
108    int rbs;
109    // charge in C    
110    sim_double q;
111    // inertial mass in m
112    sim_double m;
113    // for error message output
114    int oerrors;
115    int errors;
116    
117    // determines the numerical error and the speed of the calculation
118    TOptParam op;
119    
120    // for internal calculations
121    TVector Ai,Bi,Ci,Di;
122    TVector R,V;
123    sim_double Tn;
124    int lim;
125    
126    // helper functions
127    sim_double Gamma(sim_double w);    
128    sim_double sign(sim_double v);
129    void TimeStep(void);    
130    TVector wi(sim_double T);
131    TVector ui(sim_double T);
132    TVector Ki(sim_double T);
133    TVector simpson_rule(sim_double a, sim_double b);
134    TVector integrate(sim_double a, sim_double b, int hr);
135    void Index(int k, int& i, int& ip);
136    sim_double Tc(void);
137    TVector Core(int k);
138    void print_errors(void);
139};
140
141#endif

Implementierung

Der Code zum Download.
0/*
1 Dr. Steffen Kuehn / 2015 / steffen.kuehn@quantino-theory.de
2
3 This program is free software: you can redistribute it and/or modify
4 it under the terms of the GNU General Public License as published by
5 the Free Software Foundation, either version 3 of the License, or
6 (at your option) any later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program. If not, see <http://www.gnu.org/licenses/>.
15 */
16
17#include <math.h>
18#include <stdlib.h>
19#include <stdio.h>
20#include "charge.h"
21
22TVector::TVector(void)
23{
24    x = 0;
25    y = 0;
26    z = 0;
27}
28
29TVector::TVector(sim_double x, sim_double y, sim_double z)
30{
31    this->x = x;
32    this->y = y;
33    this->z = z;
34}
35
36TVector::TVector(const TVector& v)
37{
38    x = v.x;
39    y = v.y;
40    z = v.z;
41}
42
43TVector::~TVector()
44{
45    
46}
47
48TVector operator+(const TVector &v1, const TVector &v2)
49{
50    TVector o = v1;
51    
52    o.x += v2.x;
53    o.y += v2.y;
54    o.z += v2.z;
55    
56    return o;
57}
58
59TVector operator-(const TVector &v1, const TVector &v2)
60{
61    TVector o = v1;
62    
63    o.x -= v2.x;
64    o.y -= v2.y;
65    o.z -= v2.z;
66    
67    return o;
68}
69
70sim_double operator*(const TVector &v1, const TVector &v2)
71{
72    sim_double o = 0;
73    
74    o += v1.x*v2.x;
75    o += v1.y*v2.y;
76    o += v1.z*v2.z;
77    
78    return o;
79}
80
81TVector operator*(sim_double v1, const TVector &v2)
82{
83    TVector o;
84    
85    o.x = v1*v2.x;
86    o.y = v1*v2.y;
87    o.z = v1*v2.z;
88    
89    return o;
90}
91
92TVector operator*(const TVector &v1, sim_double v2)
93{
94    TVector o;
95    
96    o.x = v1.x*v2;
97    o.y = v1.y*v2;
98    o.z = v1.z*v2;
99    
100    return o;
101}
102
103TVector operator/(const TVector &v1, sim_double v2)
104{
105    TVector o;
106    
107    o.x = v1.x/v2;
108    o.y = v1.y/v2;
109    o.z = v1.z/v2;
110    
111    return o;
112}
113
114TVector operator-(const TVector &v)
115{
116    TVector o = v;
117    
118    o.x = -v.x;
119    o.y = -v.y;
120    o.z = -v.z;
121    
122    return o;
123}
124
125sim_double nrm(const TVector &v)
126{
127    return sqrt(v*v);
128}
129
130TOptParam::TOptParam(void)
131{    
132    n = 1000;
133    re = 0.99;
134    lim = 1000;
135}
136
137TCharge::TCharge(void)
138{
139    rbn = 0;
140    rbs = 0;    
141    q = 0;
142    m = 0;    
143    rh = NULL;
144    vh = NULL;
145    th = NULL;
146}
147
148TCharge::~TCharge()
149{
150    if (rh != NULL)    delete [] rh;
151    if (vh != NULL) delete [] vh;
152    if (th != NULL) delete [] th;
153    rh = NULL;
154    vh = NULL;
155    th = NULL;
156}
157
158void TCharge::Initialise(const TVector& r, const TVector& v,
159        int q, sim_double m, const TOptParam& op)
160{
161    rbn = 0;
162    rbs = 1;
163    TCharge::q = q*e;
164    TCharge::m = m;
165    TCharge::op = op;
166    
167    rh = new TVector [op.n];
168    vh = new TVector [op.n];
169    th = new sim_double [op.n];
170    
171    rh[0] = r;
172    vh[0] = v;
173    th[0] = 0.0;
174    
175    oerrors = 0;
176    errors = 0;
177}
178
179void TCharge::TimeStep(void)
180{    
181    rbn++;
182    if (rbn == op.n) rbn=0;    
183    if (rbs < op.n) rbs += 1;
184}
185
186void TCharge::TimeStep(const TVector& efield, sim_double dt)
187{    
188    int i = rbn;
189    TimeStep();
190    rh[rbn] = rh[i] + vh[i]*dt;
191    vh[rbn] = vh[i] + q*efield/m*dt;
192    th[rbn] = th[i] + dt;
193}
194
195void TCharge::TimeStep(const TVector& r, const TVector& v, sim_double dt)
196{
197    int i = rbn;
198    TimeStep();
199    rh[rbn] = r;
200    vh[rbn] = v;
201    th[rbn] = th[i] + dt;
202    print_errors();
203}
204
205sim_double TCharge::Gamma(sim_double w)
206{
207    // the exact form of the distribution is still unknown, especially for fast velocities.
208    return 6*w/c;
209}
210
211sim_double TCharge::sign(sim_double v)
212{
213    if (v < 0) return -1.0;
214    if (v > 0) return +1.0;
215    return 0.0;
216}
217
218TVector TCharge::wi(sim_double T)
219{
220    return (R-Ai-T*Bi)/(Tn-T)-Ci-T*Di;
221}
222
223TVector TCharge::ui(sim_double T)
224{
225    return (R-Ai-T*Bi)/(Tn-T)-V;
226}
227
228TVector TCharge::Ki(sim_double T)
229{
230    return (ui(T)*ui(T)*wi(T)*sign(ui(T)*wi(T))*Gamma(nrm(wi(T))))/(pow(Tn-T,3)*pow(nrm(wi(T)),3));
231}
232
233TVector TCharge::simpson_rule(sim_double a, sim_double b)
234{
235    return ((b-a)/6)*(Ki(a) + 4*Ki((b+a)/2) + Ki(b));
236}
237
238void TCharge::print_errors(void)
239{
240    if (errors != oerrors)
241    {
242        int nerror = errors ^ oerrors;
243        if (nerror & (1 << 0))
244        {
245            printf("lim is to small\n");
246        }
247        if (nerror & (1 << 1))
248        {
249            printf("code should not be reached!\n");
250        }
251        if (nerror & (1 << 2))
252        {
253            printf("code should not be reached!\n");
254        }
255        oerrors = errors;
256    }
257}
258
259TVector TCharge::integrate(sim_double a, sim_double b, int hr)
260{
261    TVector r1,r2;
262    sim_double h = (b-a)/2;    
263    lim++;
264    
265    r1 = simpson_rule(a,b);    
266    r2 = simpson_rule(a,a+h) + simpson_rule(a+h,b);
267    
268    sim_double ar1 = nrm(r1);
269    sim_double ar2 = nrm(r2);
270    
271    sim_double re = ar1/ar2;
272    if (re > 1) re = 1.0/re;
273    
274    if (isnan(re)) return r1;
275    
276    if (re > op.re)
277    {
278        return r2;
279    }
280    else if (lim > op.lim)
281    {
282        errors |= (1<<0);
283        return r2;
284    }
285    else
286    {
287        return integrate(a,a+h,hr+1) + integrate(a+h,b,hr+1);
288    }
289}
290
291sim_double TCharge::Tc(void)
292{
293    TVector Ei = R - Ai - V*Tn;
294    TVector Fi = V - Bi;
295    sim_double EiFi = Ei*Fi;
296    sim_double FiFi = Fi*Fi;
297    sim_double EiEi = Ei*Ei;
298    
299    sim_double sr = pow((c*c*Tn + EiFi),2) + (c*c-FiFi)*(EiEi-c*c*Tn*Tn);
300    if (sr < 0)
301    {
302        // light is too slow to reach this point in such short a time
303        return NAN;
304    }
305    sim_double v1 = c*c*Tn + EiFi - sqrtl(sr);
306    sim_double v2 = c*c-FiFi;
307    return v1/v2;
308}
309
310TVector TCharge::CoulombEfield(const TVector& r, const TVector& v)
311{
312    TVector u = v-vh[rbn];
313    sim_double mt = 1.5*((r-rh[rbn])*u)/nrm(r-rh[rbn]);    
314    return (mu0*q)/(4*M_PI)*(r-rh[rbn])/pow(nrm(r-rh[rbn]),3)*((c*c + u*u) - mt*mt);
315}
316
317void TCharge::Index(int k, int& i, int& ip)
318{
319    i = rbn - k;
320    ip = i + 1;
321    if (i < 0)
322    {
323        i = i + op.n;
324    }
325    if (ip < 0)
326    {
327        ip = ip + op.n;
328    }
329}
330
331TVector TCharge::Core(int k)
332{
333    TVector Pe;    
334    int i,ip;
335    
336    Index(k,i,ip);
337
338    sim_double dt = (th[ip] - th[i]);
339    Ai = (rh[i] * th[ip] - rh[ip] * th[i]) / dt;
340    Bi = (rh[ip] - rh[i]) / dt;
341    Ci = (vh[i] * th[ip] - vh[ip] * th[i]) / dt;
342    Di = (vh[ip] - vh[i]) / dt;
343
344    sim_double abs_ui_a = nrm(ui(th[i])) / c;
345    sim_double abs_ui_b = nrm(ui(th[ip])) / c;
346    lim = 0;
347
348    if ((abs_ui_a <= 1) && (abs_ui_b <= 1))
349    {
350        Pe = integrate(th[i], th[ip],1);
351    }
352    else if ((abs_ui_a > 1) && (abs_ui_b <= 1))
353    {
354        sim_double a = Tc();
355        if (a > th[ip])
356        {
357            errors |= (1 << 1);
358        }
359        if ((a < th[ip]) && (!isnan(a)))
360        {
361            Pe = integrate(a, th[ip], 1);
362        }
363    }
364    else if ((abs_ui_a <= 1) && ((abs_ui_b > 1) || (isnan(abs_ui_b))))
365    {
366        sim_double b = Tc();
367        if (b < th[i])
368        {
369            errors |= (1 << 2);
370        }
371        if ((b > th[i]) && (!isnan(b)))
372        {
373            Pe = integrate(th[i], b,1);
374        }
375    }
376    
377    return Pe;
378}
379
380TVector TCharge::Efield(const TVector& r, const TVector& v)
381{
382    TVector Ef;
383    sim_double fc = q*mu0/(8*M_PI);
384    
385    R = r;
386    V = v;
387    Tn = th[rbn];
388    
389    for (int k=1; k<rbs; k++)
390    {
391        TVector p = fc*Core(k);
392        if (!isnan(nrm(p)))
393        {
394            Ef = Ef + p;
395        }
396    }
397    
398    print_errors();
399    
400    return Ef;
401}