OSHUN  beta
Arbitrary Order Spherical-Harmonic 1D-3P Vlasov-Fokker-Planck-Maxwell code
vlasov.cpp
Go to the documentation of this file.
1 
8 //--------------------------------------------------------------
9 // Standard libraries
10 #include <iostream>
11 #include <vector>
12 #include <valarray>
13 #include <complex>
14 #include <algorithm>
15 #include <cstdlib>
16 
17 #include <math.h>
18 #include <map>
19 
20 // My libraries
21 #include "lib-array.h"
22 #include "lib-algorithms.h"
23 #include "nmethods.h"
24 
25 // Declerations
26 #include "state.h"
27 #include "formulary.h"
28 #include "input.h"
29 #include "fluid.h"
30 #include "vlasov.h"
31 
32 
33 
34 //**************************************************************
35 //--------------------------------------------------------------
36 // Current
37 
38 Current_1D::Current_1D( double pmin, double pmax, size_t Np,
39  size_t Nx )
40  : Jx(Nx), Jy(Nx), Jz(Nx), small(0.5*pmax/Np) {
41 // pr(Algorithms::MakeAxis(complex<double>(pmin),
42 // complex<double>(pmax),
43 // Np)),
44 // invg(pr) {
45  small *= small; small *= small; small *= 0.5*pmax/Np;
46  small *= 0.2; small *= 1.0/(1.5*pmax/Np);
47 // for (size_t i(0); i < pr.size(); ++i) {
48 // invg[i] = 1.0 / (sqrt(1.0+pr[i]*pr[i]));
49 // }
50 };
51 //--------------------------------------------------------------
52 
53 void Current_1D::operator()(const DistFunc1D& Din, Field1D& Exh, Field1D& Eyh, Field1D& Ezh) {
54 
55  Array2D<double> temp(3,Din(0).numx());
56 
57  temp = Din.getcurrent();
58 
59  for (size_t i(0); i < Jx.numx(); ++i) {
60  Jx(i) = complex<double >(temp(0,i));
61  Jy(i) = complex<double >(temp(1,i));
62  Jz(i) = complex<double >(temp(2,i));
63  }
64 
65  Exh += Jx;
66  Eyh += Jy;
67  Ezh += Jz;
68 
69 }
70 
71 void Current_1D::es1d(const DistFunc1D& Din, Field1D& Exh) {
72 
73  valarray<double> temp(Din(0).numx());
74 
75  temp = Din.getcurrent(0);
76 
77  for (size_t i(0); i < Jx.numx(); ++i) {
78  Jx(i) = complex<double >(temp[i]);//+4.0/3.0*3.1415926*small*Din(1,0)(1,i);
79  }
80 
81  Exh += Jx;
82 
83 }
84 //--------------------------------------------------------------
85 //**************************************************************
86 
87 
88 //**************************************************************
89 //--------------------------------------------------------------
90 Electric_Field_1D::Electric_Field_1D(size_t Nl, size_t Nm,
91  double pmin, double pmax, size_t Np,
92  double xmin, double xmax, size_t Nx)
93 //--------------------------------------------------------------
94 // Constructor
95 //--------------------------------------------------------------
96  : A1(Nl+1,Nm+1), A2(Nl+1,Nm+1),
97  B1(Nl+1), B2(Nl+1),
98  C1(Nl+1), C2(Nl+1,Nm+1), C3(Nl+1), C4(Nl+1,Nm+1),
99  Hp0(Nl+1),
100  H(Np,Nx), G(Np,Nx), TMP(Np,Nx),
101  // pr(Algorithms::MakeAxis(complex<double>(complex<double>(pmax/2.0/Np)),
102  // complex<double>(pmax),
103  // Np)),
104 
105  pr(Algorithms::MakeCAxis(complex<double>(0.0),
106  complex<double>(pmax),
107  Np)),
108  invpr(pr)
109 {
110 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
111  complex<double> lc, mc;
112 
113 // Inverted momentum axis
114  for (size_t i(0); i < pr.size(); ++i) {
115  invpr[i] = 1.0/pr[i];
116  }
117  double idp = (-1.0)/ (2.0*(pmax-pmin)/double(Np)); // -1/(2dp)
118  // double idp = (-1.0)/ (2.0*(pmax-pmax/2.0/double(Np))/double(Np-1));
119 
120 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
121 // Calculate the A1 * l/(l+1) * 1/(2Dp), A2 * (-1)/(2Dp) parameters
122  for (size_t l(1); l < Nl+1; ++l){
123  for (size_t m(0); m<((Nm<l)?Nm:l)+1; ++m){
124  lc = complex<double>(l);
125  mc = complex<double>(m);
126  A1(l,m) = idp *(-1.0) * (lc+1.0-mc)/(2.0*lc+1.0) * lc/(lc+1.0);
127  A2(l,m) = idp * (lc+mc)/(2.0*lc+1.0);
128  }
129  }
130  A1(0,0) = idp;
131  // A2[0] is not used
132 
133 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
134 
135 // Calculate the "B1, B2" parameters
136 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
137  for (size_t l(1); l<Nl+1; ++l){
138  lc = complex<double>(l);
139  B1[l] = idp * lc* lc /(2.0*lc+1.0);
140  B2[l] = idp * lc*(lc+1.0)/(2.0*lc+1.0);
141  }
142  B2[0] = 1.0;
143 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
144 
145 // Calculate the "C1, C3" parameters
146 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
147  for (size_t l(1); l<Nl+1; ++l){
148  lc = complex<double>(l);
149  C1[l] = (-0.5) * idp * lc /((2.0*lc+1.0)*(lc+1.0));
150  C3[l] = (-0.5) * idp /(2.0*lc+1.0);
151  }
152  C1[0] = 0.5 * idp;
153 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
154 
155 // Calculate the "C2, C4" parameters
156 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
157  for (size_t l(2); l<Nl+1; ++l){
158  for (size_t m(2); m<((Nm<l)?Nm:l)+1; ++m){
159  lc = complex<double>(l);
160  mc = complex<double>(m);
161  C2(l,m) = (0.5) * idp * lc * (lc-mc+2.0)*(lc-mc+1.0)/((2.0*lc+1.0)*(lc+1.0));
162  C4(l,m) = (0.5) * idp * (lc+mc-1.0)*(lc+mc)/(2.0*lc+1.0);
163  }
164  }
165 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
166 
167 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
168 // -2*Dp*H at the 0 momentum cell
169  Hp0[0] = 1.0 / pr[0];
170  for (size_t l(1); l < Nl+1; ++l) {
171  double ld(l);
172  Hp0[l] = Hp0[l-1] * (pr[0]/pr[1]) * (2.0*ld+1.0)/(2.0*ld-1.0);
173  }
174  Hp0 *= (-2.0)*(pr[1]-pr[0]);
175 
176  A100 = complex<double>(idp);
177  C100 = complex<double>(0.5*idp);
178  A210 = complex<double>(1.0/3.0*idp);
179  B211 = complex<double>(2.0/3.0);
180  A310 = complex<double>(2.0/5.0*idp);
181  C311 = complex<double>(-0.5*idp/5.0);
182 
183 }
184 //--------------------------------------------------------------
185 
186 //--------------------------------------------------------------
187 void Electric_Field_1D::operator()(const DistFunc1D& Din,
188  const Field1D& FEx, const Field1D& FEy, const Field1D& FEz,
189  DistFunc1D& Dh) {
190 //--------------------------------------------------------------
191 // This is the core calculation for the electric field
192 //--------------------------------------------------------------
193 
194  complex<double> ii(0.0,1.0);
195 
196  valarray<complex<double> > Ex(FEx.array());
197  valarray<complex<double> > Em(FEz.array());
198  Em *= (-1.0)*ii;
199  Em += FEy.array();
200  valarray<complex<double> > Ep(FEz.array());
201  Ep *= ii;
202  Ep += FEy.array();
203 
204  Ex *= Din.q();
205  Em *= Din.q();
206  Ep *= Din.q();
207 
208  size_t l0(Din.l0());
209  size_t m0(Din.m0());
210 
211 
212 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
213 // m = 0, l = 0
214 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215  MakeG00(Din(0,0));
216  Ex *= A1(0,0); TMP = G; Dh(1,0) += G.mxaxis(Ex);
217  Em *= C1[0]; Dh(1,1) += TMP.mxaxis(Em);
218 
219 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
220 // m = 0, l = 1
221 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
222  MakeGH(Din(1,0),1);
223  Ex *= A2(1,0) / A1(0,0); Dh(0,0) += H.mxaxis(Ex);
224  Ex *= A1(1,0) / A2(1,0); TMP = G; Dh(2,0) += G.mxaxis(Ex);
225  Em *= C1[1] / C1[0] ; Dh(2,1) += TMP.mxaxis(Em);
226 
227 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
228 // m = 0, 1 < l < l0
229 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
230  for (size_t l(2); l < l0; ++l){
231  MakeGH(Din(l,0),l);
232  Ex *= A2(l,0) / A1(l-1,0); TMP = H; Dh(l-1,0) += H.mxaxis(Ex);
233  Em *= C3[l] / C1[l-1]; Dh(l-1,1) += TMP.mxaxis(Em);
234  Ex *= A1(l,0) / A2(l,0); TMP = G; Dh(l+1,0) += G.mxaxis(Ex);
235  Em *= C1[l] / C3[l]; Dh(l+1,1) += TMP.mxaxis(Em);
236  }
237 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
238 // m = 0, l = l0
239 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
240  MakeGH(Din(l0,0),l0);
241  Ex *= A2(l0,0) / A1(l0-1,0); TMP = H; Dh(l0-1,0) += H.mxaxis(Ex);
242  Em *= C3[l0] / C1[l0-1]; Dh(l0-1,1) += TMP.mxaxis(Em);
243 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
244 
245 
246 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
247 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
248 // m = 1, l = 1
249 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
250  MakeGH(Din(1,1),1);
251  Ep *= B2[1]; H = H.mxaxis(Ep); Dh(0,0) += H.Re();
252  Ex *= A1(1,1) / A2(l0,0); TMP = G; Dh(2,1) += G.mxaxis(Ex);
253  Em *= C1[1] / C3[l0]; G = TMP; Dh(2,2) += TMP.mxaxis(Em);
254  Ep *= B1[1] / B2[1]; G = G.mxaxis(Ep); Dh(2,0) += G.Re();
255 
256 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
257 // m = 1, l = 2
258 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
259  MakeGH(Din(2,1),2);
260  Ex *= A2(2,1) / A1(1,1); TMP = H; Dh(1,1) += TMP.mxaxis(Ex);
261  Ep *= B2[2] / B1[1]; H = H.mxaxis(Ep); Dh(1,0) += H.Re();
262  Ex *= A1(2,1) / A2(2,1); TMP = G; Dh(3,1) += G.mxaxis(Ex);
263  Em *= C1[2] / C1[1]; G = TMP; Dh(3,2) += TMP.mxaxis(Em);
264  Ep *= B1[2] / B2[2]; G = G.mxaxis(Ep); Dh(3,0) += G.Re();
265 
266 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
267 // m = 1, 1 < l < l0
268 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
269  for (size_t l(3); l < l0; ++l){
270  MakeGH(Din(l,1),l);
271  Ex *= A2(l,1) / A1(l-1,1); TMP = H; Dh(l-1,1) += H.mxaxis(Ex);
272  Em *= C3[l] / C1[l-1]; H = TMP; Dh(l-1,2) += TMP.mxaxis(Em);
273  Ep *= B2[l] / B1[l-1]; H = H.mxaxis(Ep); Dh(l-1,0) += H.Re();
274  Ex *= A1(l,1) / A2(l,1); TMP = G; Dh(l+1,1) += G.mxaxis(Ex);
275  Em *= C1[l] / C3[l]; G = TMP; Dh(l+1,2) += TMP.mxaxis(Em);
276  Ep *= B1[l] / B2[l]; G = G.mxaxis(Ep); Dh(l+1,0) += G.Re();
277  }
278 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
279 // m = 1, l = l0
280 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
281  MakeGH(Din(l0,1),l0);
282  Ex *= A2(l0,1) / A1(l0-1,1); TMP = H; Dh(l0-1,1) += H.mxaxis(Ex);
283  Em *= C3[l0] / C1[l0-1]; H = TMP; Dh(l0-1,2) += TMP.mxaxis(Em);
284  Ep *= B2[l0] / B1[l0-1]; H = H.mxaxis(Ep); Dh(l0-1,0) += H.Re();
285 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
286 
287 
288 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
289  C4(l0,1) = B2[l0];
290  for (size_t m(2); m < m0; ++m){
291 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
292 // m > 1 , l = m
293 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
294  MakeGH(Din(m,m),m);
295  Ep *= C4(m,m) / C4(l0,m-1); Dh(m-1,m-1) += H.mxaxis(Ep);
296  Ex *= A1(m,m) / A2(l0,m-1); TMP = G; Dh(m+1,m ) += G.mxaxis(Ex);
297  Em *= C1[m] / C3[l0]; G = TMP; Dh(m+1,m+1) += TMP.mxaxis(Em);
298  Ep *= C2(m,m) / C4(m,m); Dh(m+1,m-1) += G.mxaxis(Ep);
299 
300 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
301 // m > 1 , l = m+1
302 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
303  MakeGH(Din(m+1,m),m+1);
304  Ex *= A2(m+1,m) / A1(m,m); TMP = H; Dh(m ,m ) += TMP.mxaxis(Ex);
305  Ep *= C4(m+1,m) / C2(m,m); Dh(m ,m-1) += H.mxaxis(Ep);
306  if ( m+1 < l0) { //always true except when m = m0-1 = l0-1
307  Ex *= A1(m+1,m) / A2(m+1,m); TMP = G; Dh(m+2,m ) += G.mxaxis(Ex);
308  Em *= C1[m+1] / C1[m]; G = TMP; Dh(m+2,m+1) += TMP.mxaxis(Em);
309  Ep *= C2(m+1,m) / C4(m+1,m); Dh(m+2,m-1) += G.mxaxis(Ep);
310 
311 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
312 // m > 1, 1 < l < l0
313 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
314  for (size_t l(m+2); l < l0; ++l){
315  MakeGH(Din(l,m),l);
316  Ex *= A2(l,m) / A1(l-1,m); TMP = H; Dh(l-1,m ) += H.mxaxis(Ex);
317  Em *= C3[l] / C1[l-1]; H = TMP; Dh(l-1,m+1) += TMP.mxaxis(Em);
318  Ep *= C4(l,m) / C2(l-1,m); Dh(l-1,m-1) += H.mxaxis(Ep);
319  Ex *= A1(l,m) / A2(l,m); TMP = G; Dh(l+1,m ) += G.mxaxis(Ex);
320  Em *= C1[l] / C3[l]; G = TMP; Dh(l+1,m+1) += TMP.mxaxis(Em);
321  Ep *= C2(l,m) / C4(l,m); Dh(l+1,m-1) += G.mxaxis(Ep);
322  }
323 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
324 // m > 1, l = l0
325 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
326  MakeGH(Din(l0,m),l0);
327  Ex *= A2(l0,m) / A1(l0-1,m); TMP = H; Dh(l0-1,m ) += H.mxaxis(Ex);
328  Em *= C3[l0] / C1[l0-1]; H = TMP; Dh(l0-1,m+1) += TMP.mxaxis(Em);
329  Ep *= C4(l0,m) / C2(l0-1,m); Dh(l0-1,m-1) += H.mxaxis(Ep);
330  }
331  }
332 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
333 
334 
335 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
336  MakeGH(Din(m0,m0),m0);
337  Ep *= C4(m0,m0) / C4(l0,m0-1); Dh(m0-1,m0-1) += H.mxaxis(Ep);
338 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
339 // m = m0, l0 > l > m0
340 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
341  if ( m0 < l0) {
342  Ex *= A1(m0,m0) / A2(l0,m0-1); TMP = G; Dh(m0+1,m0 ) += TMP.mxaxis(Ex);
343  Ep *= C2(m0,m0) / C4(m0,m0); Dh(m0+1,m0-1) += G.mxaxis(Ep);
344 
345 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
346 // m = m0 , l = m0+1
347 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
348  MakeGH(Din(m0+1,m0),m0+1);
349  Ex *= A2(m0+1,m0) / A1(m0,m0); TMP = H; Dh(m0,m0 ) += TMP.mxaxis(Ex);
350  Ep *= C4(m0+1,m0) / C2(m0,m0); Dh(m0,m0-1) += H.mxaxis(Ep);
351  if ( m0+1 < l0) {
352  Ex *= A1(m0+1,m0) / A2(m0+1,m0); TMP = G; Dh(m0+2,m0 )+= TMP.mxaxis(Ex);
353  Ep *= C2(m0+1,m0) / C4(m0+1,m0); Dh(m0+2,m0-1)+= G.mxaxis(Ep);
354 
355 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
356 // m = m0, m0+2 < l < l0
357 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
358  for (size_t l(m0+2); l < l0; ++l){
359  MakeGH(Din(l,m0),l);
360  Ex *= A2(l,m0) / A1(l-1,m0); TMP = H; Dh(l-1,m0) += TMP.mxaxis(Ex);
361  Ep *= C4(l,m0) / C2(l-1,m0); Dh(l-1,m0-1) += H.mxaxis(Ep);
362  Ex *= A1(l,m0) / A2(l, m0); TMP = G; Dh(l+1,m0 ) += TMP.mxaxis(Ex);
363  Ep *= C2(l,m0) / C4(l, m0); Dh(l+1,m0-1) += G.mxaxis(Ep);
364  }
365 
366 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
367 // m > 1, l = l0
368 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
369  MakeGH(Din(l0,m0),l0);
370  Ex *= A2(l0,m0) / A1(l0-1,m0); TMP = H; Dh(l0-1,m0) += TMP.mxaxis(Ex);
371  Ep *= C4(l0,m0) / C2(l0-1,m0); Dh(l0-1,m0-1) += H.mxaxis(Ep);
372  }
373  }
374 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
375 
376 }
377 //--------------------------------------------------------------
378 //--------------------------------------------------------------
379 void Electric_Field_1D::es1d(const DistFunc1D& Din,
380  const Field1D& FEx, const Field1D& FEy, const Field1D& FEz,
381  DistFunc1D& Dh) {
382 //--------------------------------------------------------------
383 // This is the core calculation for the electric field
384 //--------------------------------------------------------------
385 
386  complex<double> ii(0.0,1.0);
387 
388  valarray<complex<double> > Ex(FEx.array());
389  // valarray<complex<double> > Em(FEz.array());
390  // Em *= (-1.0)*ii;
391  // Em += FEy.array();
392  // valarray<complex<double> > Ep(FEz.array());
393  // Ep *= ii;
394  // Ep += FEy.array();
395 
396  Ex *= Din.q();
397  // Em *= Din.q();
398  // Ep *= Din.q();
399 
400  size_t l0(Din.l0());
401  // size_t m0(Din.m0());
402 
403 
404 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
405 // m = 0, l = 0
406 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
407  MakeG00(Din(0,0));
408  Ex *= A1(0,0); TMP = G; Dh(1,0) += G.mxaxis(Ex);
409 
410 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
411 // m = 0, l = 1
412 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
413  MakeGH(Din(1,0),1);
414  Ex *= A2(1,0) / A1(0,0); Dh(0,0) += H.mxaxis(Ex);
415  Ex *= A1(1,0) / A2(1,0); TMP = G; Dh(2,0) += G.mxaxis(Ex);
416 
417 
418 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
419 // m = 0, 1 < l < l0
420 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
421  for (size_t l(2); l < l0; ++l){
422  MakeGH(Din(l,0),l);
423  Ex *= A2(l,0) / A1(l-1,0); TMP = H; Dh(l-1,0) += H.mxaxis(Ex);
424 
425  Ex *= A1(l,0) / A2(l,0); TMP = G; Dh(l+1,0) += G.mxaxis(Ex);
426 
427  }
428 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
429 // m = 0, l = l0
430 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
431  MakeGH(Din(l0,0),l0);
432  Ex *= A2(l0,0) / A1(l0-1,0); TMP = H; Dh(l0-1,0) += H.mxaxis(Ex);
433 
434 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
435 
436 
437 }
438 //--------------------------------------------------------------
439 
440 
441 
442 //--------------------------------------------------------------
443 void Electric_Field_1D::Implicit_Ex(const DistFunc1D& Din, const Field1D& FEx, DistFunc1D& Dh) {
444 //--------------------------------------------------------------
445 // This is the calculation for the implicit electric field in x,
446 // which is to say Ex as it acts on the first few harmonics.
447 //--------------------------------------------------------------
448 
449  valarray<complex<double> > Ex(FEx.array());
450  Ex *= Din.q();
451 
452 // m = 0, l = 0
453  MakeG00(Din(0,0));
454  Ex *= A1(0,0); Dh(1,0) += G.mxaxis(Ex);
455 
456 // m = 0, l = 1
457  MakeGH(Din(1,0),1);
458  Ex *= A2(1,0) / A1(0,0); Dh(0,0) += H.mxaxis(Ex);
459  Ex *= A1(1,0) / A2(1,0); Dh(2,0) += G.mxaxis(Ex);
460 
461 // m = 0, l = 2
462  MakeGH(Din(2,0),2);
463  Ex *= A2(2,0) / A1(1,0); Dh(1,0) += H.mxaxis(Ex);
464 
465 // m = 0, l = 3
466  MakeGH(Din(3,0),3);
467  Ex *= A2(3,0) / A2(2,0); Dh(2,0) += H.mxaxis(Ex);
468 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
469 
470 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
471 // m = 1, l = 1
472  MakeGH(Din(1,1),1);
473  Ex *= A1(1,1) / A2(3,0); Dh(2,1) += G.mxaxis(Ex);
474 
475 // m = 1, l = 2
476  MakeGH(Din(2,1),2);
477  Ex *= A2(2,1) / A1(1,1); Dh(1,1) += H.mxaxis(Ex);
478 
479 // m = 1, l = 3
480  MakeGH(Din(3,1),3);
481  Ex *= A2(3,1) / A2(2,1); Dh(2,1) += H.mxaxis(Ex);
482 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
483 
484 
485 }
486 //--------------------------------------------------------------
487 //--------------------------------------------------------------
488 void Electric_Field_1D::Implicit_Ey(const DistFunc1D& Din, const Field1D& FEy, DistFunc1D& Dh) {
489 //--------------------------------------------------------------
490 // This is the calculation for the implicit electric field in y
491 // and z, which is to say Ex as it acts on the first few harmonics.
492 //--------------------------------------------------------------
493 
494  valarray<complex<double> > Em(FEy.array());
495  valarray<complex<double> > Ep(FEy.array());
496 
497 
498  Em *= Din.q();;
499  Ep *= Din.q();;
500 
501 // m = 0, l = 0
502  MakeG00(Din(0,0));
503  Em *= C1[0]; Dh(1,1) += G.mxaxis(Em);
504 
505 // m = 0, l = 1
506  MakeGH(Din(1,0),1);
507  Em *= C1[1] / C1[0]; Dh(2,1) += G.mxaxis(Em);
508 
509 // m = 0, 1 < l < l0
510  MakeGH(Din(2,0),2);
511  Em *= C3[2] / C1[1]; Dh(1,1) += H.mxaxis(Em);
512 
513 // m = 0, l = 3
514  MakeGH(Din(3,0),3);
515  Em *= C3[3] / C3[2]; Dh(2,1) += H.mxaxis(Em);
516 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
517 
518 
519 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
520 // m = 1, l = 1
521  MakeGH(Din(1,1),1);
522  Ep *= B2[1]; H = H.mxaxis(Ep); Dh(0,0) += H.Re();
523  Em *= C1[1] / C3[3]; TMP = G; Dh(2,2) += TMP.mxaxis(Em);
524  Ep *= B1[1] / B2[1]; G = G.mxaxis(Ep); Dh(2,0) += G.Re();
525 
526 // m = 1, l = 2
527  MakeGH(Din(2,1),2);
528  Ep *= B2[2] / B1[1]; H = H.mxaxis(Ep); Dh(1,0) += H.Re();
529 
530 // m = 1, l = 3
531  MakeGH(Din(3,1),3);
532  Em *= C3[3] / C1[1]; TMP = H; Dh(2,2) += TMP.mxaxis(Em);
533  Ep *= B2[3] / B2[2]; H = H.mxaxis(Ep); Dh(2,0) += H.Re();
534 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
535 
536 
537 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
538 // m = 2, l = 2
539  MakeGH(Din(2,2),2);
540  Ep *= C4(2,2) / B2[3]; Dh(1,1) += H.mxaxis(Ep);
541 
542 // m = 2, l = 3
543  MakeGH(Din(3,2),3);
544  Ep *= C4(3,2) / C4(2,2); Dh(2,1) += H.mxaxis(Ep);
545 
546  size_t m0(Din.m0());
547  if ( m0 > 2) {
548 // m = 3, l = 3
549  MakeGH(Din(3,3),3);
550  Ep *= C4(3,3) / C4(3,2); Dh(2,2) += H.mxaxis(Ep);
551  }
552 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
553 }
554 //
556 void Electric_Field_1D::Implicit_Ez(const DistFunc1D& Din, const Field1D& FEz, DistFunc1D& Dh) {
557 //--------------------------------------------------------------
558 // This is the calculation for the implicit electric field in y
559 // and z, which is to say Ex as it acts on the first few harmonics.
560 //--------------------------------------------------------------
561  complex<double> ii(0.0,1.0);
562 
563 
564  valarray<complex<double> > Em(FEz.array());
565  Em *= (-1.0)*ii;
566 
567 
568  valarray<complex<double> > Ep(FEz.array());
569  Ep *= ii;
570  // Ep += FEy.array();
571  //
572  // Ex *= Din.q();;
573  Em *= Din.q();
574  Ep *= Din.q();
575 
576 
577 
578 // m = 0, l = 0
579  MakeG00(Din(0,0));
580  Em *= C1[0]; Dh(1,1) += G.mxaxis(Em);
581 
582 // m = 0, l = 1
583  MakeGH(Din(1,0),1);
584  Em *= C1[1] / C1[0]; Dh(2,1) += G.mxaxis(Em);
585 
586 // m = 0, 1 < l < l0
587  MakeGH(Din(2,0),2);
588  Em *= C3[2] / C1[1]; Dh(1,1) += H.mxaxis(Em);
589 
590 // m = 0, l = 3
591  MakeGH(Din(3,0),3);
592  Em *= C3[3] / C3[2]; Dh(2,1) += H.mxaxis(Em);
593 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
594 
595 
596 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
597 // m = 1, l = 1
598  MakeGH(Din(1,1),1);
599  Ep *= B2[1]; H = H.mxaxis(Ep); Dh(0,0) += H.Re();
600  Em *= C1[1] / C3[3]; TMP = G; Dh(2,2) += TMP.mxaxis(Em);
601  Ep *= B1[1] / B2[1]; G = G.mxaxis(Ep); Dh(2,0) += G.Re();
602 
603 // m = 1, l = 2
604  MakeGH(Din(2,1),2);
605  Ep *= B2[2] / B1[1]; H = H.mxaxis(Ep); Dh(1,0) += H.Re();
606 
607 // m = 1, l = 3
608  MakeGH(Din(3,1),3);
609  Em *= C3[3] / C1[1]; TMP = H; Dh(2,2) += TMP.mxaxis(Em);
610  Ep *= B2[3] / B2[2]; H = H.mxaxis(Ep); Dh(2,0) += H.Re();
611 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
612 
613 
614 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
615 // m = 2, l = 2
616  MakeGH(Din(2,2),2);
617  Ep *= C4(2,2) / B2[3]; Dh(1,1) += H.mxaxis(Ep);
618 
619 // m = 2, l = 3
620  MakeGH(Din(3,2),3);
621  Ep *= C4(3,2) / C4(2,2); Dh(2,1) += H.mxaxis(Ep);
622 
623  size_t m0(Din.m0());
624  if ( m0 > 2) {
625 // m = 3, l = 3
626  MakeGH(Din(3,3),3);
627  Ep *= C4(3,3) / C4(3,2); Dh(2,2) += H.mxaxis(Ep);
628  }
629 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
630 }
631 //--------------------------------------------------------------
632 //--------------------------------------------------------------
633 void Electric_Field_1D::f1only(const DistFunc1D& Din,
634  const Field1D& FEx, const Field1D& FEy, const Field1D& FEz,
635  DistFunc1D& Dh) {
636 //--------------------------------------------------------------
637 // This is the core calculation for the electric field
638 //--------------------------------------------------------------
639 
640  complex<double> ii(0.0,1.0);
641 
642  valarray<complex<double> > Ex(FEx.array());
643  valarray<complex<double> > Em(FEz.array());
644  Em *= (-1.0)*ii;
645  Em += FEy.array();
646  valarray<complex<double> > Ep(FEz.array());
647  Ep *= ii;
648  Ep += FEy.array();
649 
650  Ex *= Din.q();
651  Em *= Din.q();
652  Ep *= Din.q();
653 
654  size_t l0(Din.l0());
655  size_t m0(Din.m0());
656 
657  SHarmonic1D f20(Din(0,0));
658  f20 *= complex<double>(1.0/3.0);
659 
660 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
661 // m = 0, l = 0
662 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
663  MakeG00(Din(0,0));
664  Ex *= A100; TMP = G; Dh(1,0) += G.mxaxis(Ex);
665  Em *= C100; Dh(1,1) += TMP.mxaxis(Em);
666 
667 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
668 // m = 0, l = 1
669 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
670  MakeGH(Din(1,0),1);
671  Ex *= A210 / A100; Dh(0,0) += H.mxaxis(Ex);
672 
673 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
674 // m = 0, 1 < l < l0
675 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
676 
677 // MakeGH(f20,2);
678 // Ex *= A310 / A210; TMP = H; Dh(1,0) += H.mxaxis(Ex);
679 // Em *= C311 / C100; Dh(1,1) += TMP.mxaxis(Em);
680 
681 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
682 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
683 // m = 1, l = 1
684 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
685  MakeGH(Din(1,1),1);
686  Ep *= B211; H = H.mxaxis(Ep); Dh(0,0) += H.Re();
687 
688 }
689 //--------------------------------------------------------------
690 
691 
692 //--------------------------------------------------------------
693 void Electric_Field_1D::Implicit_Ex_f1only(const DistFunc1D& Din, const Field1D& FEx, DistFunc1D& Dh) {
694 //--------------------------------------------------------------
695 // This is the calculation for the implicit electric field in x,
696 // which is to say Ex as it acts on the first few harmonics.
697 //--------------------------------------------------------------
698 
699  valarray<complex<double> > Ex(FEx.array());
700  Ex *= Din.q();
701 
702  SHarmonic1D f20(Din(0,0));
703  f20 *= complex<double>(1.0/3.0);
704 
705 // m = 0, l = 0
706  MakeG00(Din(0,0));
707  Ex *= A100; Dh(1,0) += G.mxaxis(Ex);
708 
709 // m = 0, l = 1
710  MakeGH(Din(1,0),1);
711  Ex *= A210 / A100; Dh(0,0) += H.mxaxis(Ex);
712  // Ex *= A1(1,0) / A2(1,0); Dh(2,0) += G.mxaxis(Ex);
713 
714 // m = 0, l = 2
715 // MakeGH(f20,2);
716 // Ex *= A310 / A210; TMP = H; Dh(1,0) += H.mxaxis(Ex);
717 
718 
719 }
720 //--------------------------------------------------------------
721 //--------------------------------------------------------------
722 void Electric_Field_1D::Implicit_Ey_f1only(const DistFunc1D& Din, const Field1D& FEy, DistFunc1D& Dh) {
723 //--------------------------------------------------------------
724 // This is the calculation for the implicit electric field in y
725 // and z, which is to say Ex as it acts on the first few harmonics.
726 //--------------------------------------------------------------
727 
728  valarray<complex<double> > Em(FEy.array());
729  valarray<complex<double> > Ep(FEy.array());
730 
731  SHarmonic1D f20(Din(0,0));
732  f20 *= complex<double>(1.0/3.0);
733 
734  Em *= Din.q();
735  Ep *= Din.q();
736 
737 // m = 0, l = 0
738  MakeG00(Din(0,0));
739  Em *= C100; Dh(1,1) += G.mxaxis(Em);
740 
741 // m = 0, l = 1
742  // MakeGH(Din(1,0),1);
743  // Em *= C1[1] / C1[0]; Dh(2,1) += G.mxaxis(Em);
744 
745 // m = 0, 1 < l < l0
746 // MakeGH(f20,2);
747 // Em *= C311 / C100; Dh(1,1) += TMP.mxaxis(Em);
748 
749 // m = 0, l = 3
750  // MakeGH(Din(3,0),3);
751  // Em *= C3[3] / C3[2]; Dh(2,1) += H.mxaxis(Em);
752 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
753 
754 
755 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
756 // m = 1, l = 1
757  MakeGH(Din(1,1),1);
758  Ep *= B211; H = H.mxaxis(Ep); Dh(0,0) += H.Re();
759 }
760 //
762 void Electric_Field_1D::Implicit_Ez_f1only(const DistFunc1D& Din, const Field1D& FEz, DistFunc1D& Dh) {
763 //--------------------------------------------------------------
764 // This is the calculation for the implicit electric field in y
765 // and z, which is to say Ex as it acts on the first few harmonics.
766 //--------------------------------------------------------------
767  complex<double> ii(0.0,1.0);
768  SHarmonic1D f20(Din(0,0));
769  f20 *= complex<double>(1.0/3.0);
770 
771  valarray<complex<double> > Em(FEz.array());
772  Em *= (-1.0)*ii;
773 
774 
775  valarray<complex<double> > Ep(FEz.array());
776  Ep *= ii;
777  // m = 0, l = 0
778  MakeG00(Din(0,0));
779  Em *= C100; Dh(1,1) += G.mxaxis(Em);
780 
781 // m = 0, l = 1
782  // MakeGH(Din(1,0),1);
783  // Em *= C1[1] / C1[0]; Dh(2,1) += G.mxaxis(Em);
784 
785 // m = 0, 1 < l < l0
786 // Din(2,0) = Din(0,0);
787 // Din(2,0) *= complex<double>(1.0/3.0);
788 //
789 // MakeGH(f20,2);
790 // Em *= C311 / C100; Dh(1,1) += TMP.mxaxis(Em);
791 
792 // // m = 0, l = 3
793 // MakeGH(Din(3,0),3);
794 // Em *= C3[3] / C3[2]; Dh(2,1) += H.mxaxis(Em);
795 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
796 
797 
798 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
799 // m = 1, l = 1
800  MakeGH(Din(1,1),1);
801  Ep *= B211; H = H.mxaxis(Ep); Dh(0,0) += H.Re();
802 }
803 //--------------------------------------------------------------
804 
805 //--------------------------------------------------------------
806 // Make derivatives 2*Dp*(l+1/l)*G and -2*Dp*H for a given f
807 void Electric_Field_1D::MakeGH(SHarmonic1D& f, size_t el){
808 //--------------------------------------------------------------
809  valarray<complex<double> > invpax(invpr);
810  complex<double> ld(el);
811 
812  invpax *= (-2.0)*(ld+1.0) * (pr[1]-pr[0]);
813 
814  G = f; H = f;
815  G = G.Dp();
816  H = H.mpaxis(invpax);
817  H += G;
818  G *= -(2.0*ld+1.0)/ld;
819  G += H;
820 
821  for (size_t i(0); i < G.numx(); ++i) G(0,i) = 0.0;
822  for (size_t i(0); i < H.numx(); ++i) H(0,i) = f(1,i) * Hp0[el];
823 }
824 //--------------------------------------------------------------
825 
826 //--------------------------------------------------------------
827 // Calculation of G00 = -2*Dp* df/dp(p0)
828 void Electric_Field_1D::MakeG00(SHarmonic1D& f) {
829 //--------------------------------------------------------------
830  G = f; G = G.Dp();
831 
832  complex<double> p0p1_sq( pr[0]*pr[0]/(pr[1]*pr[1]) ),
833  inv_mp0p1_sq( 1.0/(1.0-p0p1_sq) ),
834  g_r = -4.0*(pr[1]-pr[0]) * pr[0]/(pr[1]*pr[1]),
835  // g_r = 2.0*(pr[0]/pr[1]/pr[1]),
836  f00;
837 
838  for (size_t i(0); i < f.numx(); ++i) {
839  f00 = ( f(0,i) - f(1,i) * p0p1_sq) * inv_mp0p1_sq;
840  G(0,i) = ( f(1,i) - f00) * g_r;
841  }
842 }
843 //--------------------------------------------------------------
844 //**************************************************************
845 
846 
847 //**************************************************************
848 //--------------------------------------------------------------
849 Magnetic_Field_1D::Magnetic_Field_1D(size_t Nl, size_t Nm,
850  double pmin, double pmax, size_t Np,
851  double xmin, double xmax, size_t Nx)
852 //--------------------------------------------------------------
853 // Constructor
854 //--------------------------------------------------------------
855  : A1(Nm+1), B1(Nl+1), A2(Nl+1,Nm+1), A3(0.5),
856  FLM(Np,Nx)
857 {
858 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
859  complex<double> lc, mc;
860  complex<double> c01(0.0,1.0);
861 
862 // Calculate the "A1" parameters
863 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
864  for (size_t m(0); m < Nm+1; ++m){
865  mc = complex<double>(m);
866  A1[m] = (-1.0)*c01*mc;
867  }
868 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
869 
870 // Calculate the "A2" parameters
871 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
872  for (size_t l(1); l < Nl+1; ++l)
873  for (size_t m=0; m<((Nm<l)?Nm:l)+1; ++m){
874  lc = complex<double>(l);
875  mc = complex<double>(m);
876  A2(l,m) = (-0.5)*(lc+1.0-mc)*(lc+mc);
877  }
878 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
879 
880 // Calculate the "A3" parameters
881 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
882  A3 = 0.5;
883 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
884 
885 // Calculate the "B1" parameters
886 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
887  for (size_t l(0); l < Nl+1; ++l){
888  lc = complex<double>(l);
889  B1[l] = (-1.0)*lc*(lc+1.0);
890  }
891 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
892 
893 }
894 //--------------------------------------------------------------
895 
896 
897 //--------------------------------------------------------------
898 void Magnetic_Field_1D::operator()(const DistFunc1D& Din,
899  const Field1D& FBx, const Field1D& FBy, const Field1D& FBz,
900  DistFunc1D& Dh) {
901 //--------------------------------------------------------------
902 // This is the core calculation for the magnetic field
903 //--------------------------------------------------------------
904 
905  complex<double> ii(0.0,1.0);
906 
907  valarray<complex<double> > Bx(FBx.array());
908  valarray<complex<double> > Bm(FBy.array());
909  Bm *= (-1.0)*ii;
910  Bm += FBz.array();
911  valarray<complex<double> > Bp(FBy.array());
912  Bp *= ii;
913  Bp += FBz.array();
914 
915  Bx *= Din.q();
916  Bm *= Din.q();
917  Bp *= Din.q();
918 
919  size_t l0(B1.size()-1);
920  size_t m0(A1.size()-1);
921 
922 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
923 // m = 0, 1 < l < l0+1
924 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
925  Bp *= A3;
926  for (size_t l(1); l < l0+1; ++l){
927  FLM = Din(l,0); Dh(l,1) += FLM.mxaxis(Bp);
928  }
929 
930 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
931 // m = 1, l = 1
932 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
933  FLM = Din(1,1); Bx *= A1[1]; Dh(1,1) += FLM.mxaxis(Bx);
934  FLM = Din(1,1); Bm *= B1[1]; FLM = FLM.mxaxis(Bm); Dh(1,0) += FLM.Re();
935 
936 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
937 // m = 1, l > 1
938 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
939  for (size_t l(2); l < l0+1; ++l){
940  FLM = Din(l,1); Dh(l,2) += FLM.mxaxis(Bp);
941  FLM = Din(l,1); Dh(l,1) += FLM.mxaxis(Bx);
942  FLM = Din(l,1); Bm *= B1[l]/B1[l-1]; FLM = FLM.mxaxis(Bm); Dh(l,0) += FLM.Re();
943  }
944  Bm *= 1.0/B1[l0];
945 
946 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
947 // m > 1, l = m
948 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
949  for (size_t m(2); m < m0; ++m){
950  FLM = Din(m,m); Bx *= A1[m]/A1[m-1]; Dh(m,m ) += FLM.mxaxis(Bx);
951  FLM = Din(m,m); Bm *= A2(m,m); Dh(m,m-1) += FLM.mxaxis(Bm);
952  for (size_t l(m+1); l < l0+1; ++l){
953  FLM = Din(l,m); Dh(l,m+1) += FLM.mxaxis(Bp);
954  FLM = Din(l,m); Dh(l,m ) += FLM.mxaxis(Bx);
955  FLM = Din(l,m); Bm *= A2(l,m)/A2(l-1,m); Dh(l,m-1) += FLM.mxaxis(Bm);
956  }
957  Bm *= 1.0/A2(l0,m);
958  }
959 
960 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
961 // m = m0, l >= m0
962 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
963  FLM = Din(m0,m0); Bx *= A1[m0]/A1[m0-1]; Dh(m0,m0) += FLM.mxaxis(Bx);
964  FLM = Din(m0,m0); Bm *= A2(m0,m0)/*/A2(l0,m0-1)*/; Dh(m0,m0-1) += FLM.mxaxis(Bm);
965  for (size_t l(m0+1); l < l0+1; ++l){
966  FLM = Din(l,m0); Dh(l,m0 ) += FLM.mxaxis(Bx);
967  FLM = Din(l,m0); Bm *= A2(l,m0)/A2(l-1,m0); Dh(l,m0-1) += FLM.mxaxis(Bm);
968  }
969 
970 }
971 //--------------------------------------------------------------
972 //--------------------------------------------------------------
973 void Magnetic_Field_1D::implicit(DistFunc1D& Din,
974  const Field1D& FBx, const Field1D& FBy, const Field1D& FBz,
975  double dt) {
976 //--------------------------------------------------------------
977 // This is the core calculation for the magnetic field
978 //--------------------------------------------------------------
979 
980  complex<double> ii(0.0,1.0);
981 
982  valarray<complex<double> > Bx(FBx.array());
983  valarray<complex<double> > Bm(FBy.array());
984  Bm *= (-1.0)*ii;
985  Bm += FBz.array();
986  valarray<complex<double> > Bp(FBy.array());
987  Bp *= ii;
988  Bp += FBz.array();
989 
990  Bx *= Din.q();//Din.mass();
991  Bm *= Din.q();//Din.mass();
992  Bp *= Din.q();//Din.mass();
993 
994  size_t l0(1);
995  size_t m0(1);
996 
997  int Nm(0);
998  size_t ind(0);
999  size_t indp(0);
1000 
1001  complex<double> temp(0.0,0.0);
1002 
1003  for (size_t i(0); i<Din(0,0).numx(); ++i)
1004  {
1005  for (size_t l(1); l < l0+1; ++l)
1006  {
1007  Nm = (l<m0?l:m0);
1008  valarray<complex<double>> fin((0.0,0.0),2*Nm+1);
1009  valarray<complex<double>> RHSv((0.0,0.0),2*Nm+1);
1010  Array2D<complex<double>> LHS(2*Nm+1,2*Nm+1);
1011  Array2D<complex<double>> RHS(2*Nm+1,2*Nm+1);
1012 
1013  //* The matrices corresponding to RHS and LHS are computed
1014  // They apply to the vector
1015  // [ f_l^m ] -> 0 index but applies to Nm th harmonic
1016  // ...
1017  // [ f_l^2 ]
1018  // [ f_l^1 ]
1019  // [ f_l^0 ] -> Nm index but applies to 0th harmonic
1020  // [ f_l^-1 ]
1021  // [ f_l^-2 ]
1022  // ...
1023  // [ f_l^-m ] *//
1024  RHS(0,0) = 1.0-ii*double(Nm)*dt/2.0*Bx[i];
1025  LHS(0,0) = conj(RHS(0,0));
1026  RHS(0,1) = 0.25*dt*Bp[i];
1027  LHS(0,1) = -0.25*dt*Bp[i];
1028 
1029  for (int m(Nm-1);m>0;--m)
1030  {
1031  ind = Nm-m;
1032 
1033  RHS(ind,ind) = 1.0-ii*double(m)*dt/2.0*Bx[i];
1034  LHS(ind,ind) = conj(RHS(ind,ind));
1035 
1036  RHS(ind,ind+1) = 0.25*dt*Bp[i];
1037  LHS(ind,ind+1) = -0.25*dt*Bp[i];
1038 
1039  RHS(ind,ind-1) = -0.25*dt*(l-m)*(l+m+1)*Bm[i];
1040  LHS(ind,ind-1) = 0.25*dt*(l-m)*(l+m+1)*Bm[i];
1041  }
1042 
1043  RHS(Nm,Nm) = 1.0;
1044  LHS(Nm,Nm) = 1.0;
1045  RHS(Nm,Nm-1) = -0.25*dt*(l)*(l+1)*Bm[i];
1046  LHS(Nm,Nm-1) = 0.25*dt*(l)*(l+1)*Bm[i];
1047 
1048  //* Fill in bottom half which is a series of reflections and complex conjugates
1049  RHS(Nm,Nm+1) = conj(RHS(Nm,Nm-1));
1050  LHS(Nm,Nm+1) = conj(LHS(Nm,Nm-1));
1051 
1052  for (int m(-1);m>-Nm;--m)
1053  {
1054  ind = Nm-m;
1055  indp = Nm-abs(m);
1056 
1057  RHS(ind,ind) = conj(RHS(indp,indp));
1058  LHS(ind,ind) = conj(LHS(indp,indp));
1059 
1060  RHS(ind,ind+1) = conj(RHS(indp,indp-1));
1061  LHS(ind,ind+1) = conj(LHS(indp,indp-1));
1062 
1063  RHS(ind,ind-1) = conj(RHS(indp,indp+1));
1064  LHS(ind,ind-1) = conj(LHS(indp,indp+1));
1065  }
1066 
1067  RHS(2*Nm,2*Nm) = conj(RHS(0,0));
1068  LHS(2*Nm,2*Nm) = conj(LHS(0,0));
1069  RHS(2*Nm,2*Nm-1) = conj(RHS(0,1));
1070  LHS(2*Nm,2*Nm-1) = conj(LHS(0,1));
1071 
1072  for (size_t k(0); k < Din(0,0).nump(); ++k) {
1074  fin[0] = Din(l, Nm)(k, i);
1075  for (int m(Nm - 1); m > 1; --m) {
1076  ind = Nm - m;
1077  fin[ind] = Din(l, m)(k, i);
1078  }
1079  fin[Nm] = Din(l, 0)(k, i);
1080  for (int m(-1); m > -Nm; --m) {
1081  ind = Nm - m;
1082  fin[ind] = conj(Din(l, abs(m))(k, i));
1083  }
1084  fin[2 * Nm] = conj(Din(l, Nm)(k, i));
1085 
1087  RHSv[0] = fin[0] * RHS(0, 0) + fin[1] * RHS(0, 1);
1088  for (size_t mm(1); mm < 2 * Nm; ++mm) {
1089  RHSv[mm] = fin[mm - 1] * RHS(mm, mm - 1) + fin[mm] * RHS(mm, mm) + fin[mm + 1] * RHS(mm, mm + 1);
1090  }
1091  RHSv[2 * Nm] = fin[2 * Nm - 1] * RHS(2 * Nm, 2 * Nm - 1) + fin[2 * Nm] * RHS(2 * Nm, 2 * Nm);
1092 
1093 // std::cout << "\n\n LHS = \n";
1094 // for (size_t i(0); i < LHS.dim1(); ++i) {
1095 // std::cout << "i = " << i << " :::: ";
1096 // for (size_t j(0); j < LHS.dim2(); ++j) {
1097 // std::cout << LHS(i, j) << " ";
1098 // }
1099 // std::cout << "\n";
1100 // }
1101 
1102  Thomas_Tridiagonal(LHS,RHSv,fin);
1103 
1105  Din(l,Nm)(k,i) = fin[0];
1106  for (int m(Nm-1);m>1;--m)
1107  {
1108  ind = Nm-m;
1109  Din(l,m)(k,i) = fin[ind];
1110  }
1111  Din(l,0)(k,i) = fin[Nm];
1112 
1113  }
1114 
1115  }
1116  }
1117 
1118 }
1119 //--------------------------------------------------------------
1120 void Magnetic_Field_1D::f1only(const DistFunc1D& Din,
1121  const Field1D& FBx, const Field1D& FBy, const Field1D& FBz,
1122  DistFunc1D& Dh) {
1123 //--------------------------------------------------------------
1124 // This is the core calculation for the magnetic field
1125 //--------------------------------------------------------------
1126 
1127  complex<double> ii(0.0,1.0);
1128 
1129  valarray<complex<double> > Bx(FBx.array());
1130  valarray<complex<double> > Bm(FBy.array());
1131  Bm *= (-1.0)*ii;
1132  Bm += FBz.array();
1133  valarray<complex<double> > Bp(FBy.array());
1134  Bp *= ii;
1135  Bp += FBz.array();
1136 
1137  Bx *= Din.q();
1138  Bm *= Din.q();
1139  Bp *= Din.q();
1140 
1141  size_t l0(B1.size()-1);
1142  size_t m0(A1.size()-1);
1143 
1144 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
1145 // m = 0, 1 < l < l0+1
1146 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
1147  Bp *= A3;
1148  for (size_t l(1); l < l0+1; ++l){
1149  FLM = Din(l,0); Dh(l,1) += FLM.mxaxis(Bp);
1150  }
1151 
1152 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
1153 // m = 1, l = 1
1154 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
1155  FLM = Din(1,1); Bx *= A1[1]; Dh(1,1) += FLM.mxaxis(Bx);
1156  FLM = Din(1,1); Bm *= B1[1]; FLM = FLM.mxaxis(Bm); Dh(1,0) += FLM.Re();
1157 
1158 // // - - - - - - - - - - - - - - - - - - - - - - - - - - -
1159 // // m = 1, l > 1
1160 // // - - - - - - - - - - - - - - - - - - - - - - - - - - -
1161 // for (size_t l(2); l < l0+1; ++l){
1162 // FLM = Din(l,1); Dh(l,2) += FLM.mxaxis(Bp);
1163 // FLM = Din(l,1); Dh(l,1) += FLM.mxaxis(Bx);
1164 // FLM = Din(l,1); Bm *= B1[l]/B1[l-1]; FLM = FLM.mxaxis(Bm); Dh(l,0) += FLM.Re();
1165 // }
1166 // Bm *= 1.0/B1[l0];
1167 
1168 // // - - - - - - - - - - - - - - - - - - - - - - - - - - -
1169 // // m > 1, l = m
1170 // // - - - - - - - - - - - - - - - - - - - - - - - - - - -
1171 // for (size_t m(2); m < m0; ++m){
1172 // FLM = Din(m,m); Bx *= A1[m]/A1[m-1]; Dh(m,m ) += FLM.mxaxis(Bx);
1173 // FLM = Din(m,m); Bm *= A2(m,m); Dh(m,m-1) += FLM.mxaxis(Bm);
1174 // for (size_t l(m+1); l < l0+1; ++l){
1175 // FLM = Din(l,m); Dh(l,m+1) += FLM.mxaxis(Bp);
1176 // FLM = Din(l,m); Dh(l,m ) += FLM.mxaxis(Bx);
1177 // FLM = Din(l,m); Bm *= A2(l,m)/A2(l-1,m); Dh(l,m-1) += FLM.mxaxis(Bm);
1178 // }
1179 // Bm *= 1.0/A2(l0,m);
1180 // }
1181 
1182 // // - - - - - - - - - - - - - - - - - - - - - - - - - - -
1183 // // m = m0, l >= m0
1184 // // - - - - - - - - - - - - - - - - - - - - - - - - - - -
1185 // FLM = Din(m0,m0); Bx *= A1[m0]/A1[m0-1]; Dh(m0,m0) += FLM.mxaxis(Bx);
1186 // FLM = Din(m0,m0); Bm *= A2(m0,m0)/*/A2(l0,m0-1)*/; Dh(m0,m0-1) += FLM.mxaxis(Bm);
1187 // for (size_t l(m0+1); l < l0+1; ++l){
1188 // FLM = Din(l,m0); Dh(l,m0 ) += FLM.mxaxis(Bx);
1189 // FLM = Din(l,m0); Bm *= A2(l,m0)/A2(l-1,m0); Dh(l,m0-1) += FLM.mxaxis(Bm);
1190 // }
1191 
1192 }
1193 //--------------------------------------------------------------
1194 //**************************************************************
1195 
1196 
1197 //**************************************************************
1198 //--------------------------------------------------------------
1199 Spatial_Advection_1D::Spatial_Advection_1D(size_t Nl, size_t Nm,
1200  double pmin, double pmax, size_t Np,
1201  double xmin, double xmax, size_t Nx)
1202 //--------------------------------------------------------------
1203 // Constructor
1204 //--------------------------------------------------------------
1205  : A1(Nl+1,Nm+1), A2(Nl+1,Nm+1),
1206  fd1(Np,Nx), fd2(Np,Nx),
1207  // vr(Algorithms::MakeAxis(complex<double>(pmax/2.0/Np),
1208  // complex<double>(pmax),
1209  // Np)) {
1210 
1211  vr(Algorithms::MakeCAxis(complex<double>(0.0),
1212  complex<double>(pmax),
1213  Np)){
1214 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
1215 
1216  complex<double> lc, mc;
1217 
1218  for (size_t i(0); i < vr.size(); ++i) {
1219  vr[i] = vr[i]/(sqrt(1.0+vr[i]*vr[i]));
1220  }
1221 
1222  double idx = (-1.0) / (2.0*(xmax-xmin)/double(Nx)); // -1/(2dx)
1223 
1224 
1225 // - - - - - - - - - - - - - - - - - - - - - - - - - - -
1226 // Calculate the "A1, A2" parameters
1227  for (size_t l(0); l < Nl+1; ++l){
1228  for (size_t m=0; m<((Nm<l)?Nm:l)+1; ++m){
1229  lc = complex<double>(l);
1230  mc = complex<double>(m);
1231  A1(l,m) = idx *(-1.0) * (lc-mc+1.0) / (2.0*lc+1.0);
1232  A2(l,m) = idx *(-1.0) * (lc+mc) / (2.0*lc+1.0);
1233  }
1234  }
1235  A2(0,0) = 1.0;
1236 
1237  // Calculate the "A1, A2" parameters
1238  A00 = complex<double>(idx);
1239  A10 = complex<double>(idx/3.0);
1240  A20 = complex<double>(idx*2.0/5.0);
1241 
1242 
1243 
1244 }
1245 //--------------------------------------------------------------
1246 
1247 //--------------------------------------------------------------
1248 // Advection in x
1249 void Spatial_Advection_1D::operator()(const DistFunc1D& Din, DistFunc1D& Dh) {
1250 //--------------------------------------------------------------
1251 
1252  valarray<complex<double> > vt(vr); vt *= 1.0/Din.mass();
1253 
1254  size_t l0(Din.l0());
1255  size_t m0(Din.m0());
1256 
1257  for (size_t m(0); m < ((m0<l0)?(m0+1):m0); ++m){
1258 
1259  fd1 = Din(m,m); fd1 = fd1.Dx();
1260 
1261  vt *= A1(m,m); Dh(m+1,m) += fd1.mpaxis(vt);
1262 
1263  for (size_t l(m+1); l < l0; ++l) {
1264  fd1 = Din(l,m); fd1 = fd1.Dx();
1265 
1266  vt *= A2(l,m)/A1(l-1,m); fd2 = fd1; Dh(l-1,m) += fd1.mpaxis(vt);
1267  vt *= A1(l,m)/A2(l ,m); Dh(l+1,m) += fd2.mpaxis(vt);
1268  }
1269 
1270  fd1 = Din(l0,m); fd1 = fd1.Dx();
1271  vt *= A2(l0,m)/A1(l0-1,m); Dh(l0-1,m) += fd1.mpaxis(vt);
1272  vt *= 1.0/A2(l0,m);
1273  }
1274 
1275 
1276 }
1277 //--------------------------------------------------------------
1278 //--------------------------------------------------------------
1279 // Advection in x
1280 void Spatial_Advection_1D::es1d(const DistFunc1D& Din, DistFunc1D& Dh) {
1281 //--------------------------------------------------------------
1282 
1283  valarray<complex<double> > vt(vr); vt *= 1.0/Din.mass();
1284 
1285  size_t l0(Din.l0());
1286 
1287  fd1 = Din(0,0); fd1 = fd1.Dx();
1288 
1289  vt *= A1(0,0); Dh(1,0) += fd1.mpaxis(vt);
1290 
1291  for (size_t l(1); l < l0; ++l) {
1292  fd1 = Din(l,0); //std::cout << "\n \n before dx, l = " << l << " \n";
1293  fd1 = fd1.Dx(); //std::cout << " \n after dx\n";
1294 
1295  vt *= A2(l,0)/A1(l-1,0); fd2 = fd1; Dh(l-1,0) += fd1.mpaxis(vt);
1296  vt *= A1(l,0)/A2(l ,0); Dh(l+1,0) += fd2.mpaxis(vt);
1297  }
1298 
1299  fd1 = Din(l0,0); fd1 = fd1.Dx();
1300  vt *= A2(l0,0)/A1(l0-1,0); Dh(l0-1,0) += fd1.mpaxis(vt);
1301  // vt *= 1.0/A2(l0,0);
1302 }
1303 //--------------------------------------------------------------
1304 //--------------------------------------------------------------
1305 // Advection in x
1306 void Spatial_Advection_1D::f1only(const DistFunc1D& Din, DistFunc1D& Dh) {
1307 //--------------------------------------------------------------
1308 
1309  valarray<complex<double> > vt(vr); vt *= 1.0/Din.mass();
1310 
1311  fd1 = Din(0,0);
1312  fd1 = fd1.Dx(); vt *= A00;
1313  Dh(1,0) += (fd1.mpaxis(vt));
1314 
1315  fd1 = Din(1,0);
1316  fd1 = fd1.Dx(); vt *= A10/A00;
1317  Dh(0,0) += (fd1.mpaxis(vt));
1318 
1319 // fd1 = Din(0,0); fd1 *= complex<double>(1.0/3.0);
1320 // fd1 = fd1.Dx(); vt *= A20/A10;
1321 // Dh(1,0) += (fd1.mpaxis(vt));
1322 
1323 }
1324 //--------------------------------------------------------------
1325 //--------------------------------------------------------------
1326 //**************************************************************
1327 //--------------------------------------------------------------
1328 // Update B with E term from Faraday's Law
1329 //--------------------------------------------------------------
1330 Faraday_1D::Faraday_1D( double xmin, double xmax, size_t Nx)
1331 //--------------------------------------------------------------
1332 // Constructor
1333 //--------------------------------------------------------------
1334  : tmpE(Nx), numx(Nx) {
1335 
1336  idx = complex<double>((-1.0)/ (2.0*(xmax-xmin)/double(Nx))); // -1/(2dx)
1337 
1338 }
1339 //--------------------------------------------------------------
1340 
1341 //--------------------------------------------------------------
1342 void Faraday_1D::operator()(EMF1D& EMFin, EMF1D& EMFh) {
1343 //--------------------------------------------------------------
1344 // This is the core calculation for Faraday's Law
1345 //--------------------------------------------------------------
1346 
1347 // dBx/dt += - dEz/dy
1348 // tmpE = Fin.Ez();
1349 // tmpE *= (-1.0) * idy;
1350 // Yh.EMF().Bx() += tmpE.Dy();
1351 
1352 // dBy/dt += dEz/dx
1353  tmpE = EMFin.Ez();
1354  tmpE *= idx;
1355  EMFh.By() += tmpE.Dx();
1356  // EMFh.By()(numx-1) = 0.0;
1357 
1358 // dBz/dt += dEx/dy
1359 // tmpE = Ein.Ex();
1360 // tmpE *= idy;
1361 // Yh.EMF().Bz() += tmpE.Dy();
1362 
1363 // dBz/dt += - dEy/dx
1364  tmpE = EMFin.Ey();
1365  tmpE *= (-1.0) * idx;
1366  EMFh.Bz() += tmpE.Dx();
1367 
1368  // EMFh.Bz()(numx-1) = 0.0;
1369 
1370 }
1371 //--------------------------------------------------------------
1372 //--------------------------------------------------------------
1373 //**************************************************************
1374 
1375 //**************************************************************
1376 //--------------------------------------------------------------
1377 // Update E with B term from Faraday's Law
1378 //--------------------------------------------------------------
1379 Ampere_1D::Ampere_1D( double xmin, double xmax, size_t Nx)
1380 //--------------------------------------------------------------
1381 // Constructor
1382 //--------------------------------------------------------------
1383  : tmpB(Nx), numx(Nx) {
1384 
1385  idx = complex<double>((-1.0)/ (2.0*(xmax-xmin)/double(Nx))); // -1/(2dx)
1386 
1387 }
1388 //--------------------------------------------------------------
1389 
1390 //--------------------------------------------------------------
1391 void Ampere_1D::operator()(EMF1D& EMFin, EMF1D& EMFh) {
1392 //--------------------------------------------------------------
1393 // This is the core calculation for Ampere's Law
1394 //--------------------------------------------------------------
1395 
1396 // dEx/dt += dBz/dy
1397 // tmpB = Fin.Bz();
1398 // tmpB *= idy;
1399 // Fh.Ex() += tmpB.Dy();
1400 
1401 // dEy/dt += - dBz/dx
1402  tmpB = EMFin.Bz();
1403  tmpB *= (-1.0) * idx;
1404  EMFh.Ey() += tmpB.Dx();
1405 
1406  //.Dx();
1407  // EMFh.Ey()(numx-1) = 0.0;
1408 
1409 // dEz/dt += - dBx/dy
1410 // tmpB = Fin.Bx();
1411 // tmpB *= (-1.0) * idy;
1412 // Fh.Ez() += tmpB.Dy();
1413 
1414 // dEz/dt += dBy/dx
1415  tmpB = EMFin.By();
1416  tmpB *= idx;
1417  EMFh.Ez() += tmpB.Dx();
1418  // EMFh.Ez()(numx-1) = 0.0;
1419 
1420 }
1421 //--------------------------------------------------------------
1422 //--------------------------------------------------------------
1423 //**************************************************************
lib-algorithms.h
Algorithms
formulary.h
Plasma Formulary and Units - Declarations.
Electric_Field_1D::Implicit_Ez_f1only
void Implicit_Ez_f1only(const DistFunc1D &Din, const Field1D &FEz, DistFunc1D &Dh)
Definition: vlasov.cpp:762
Field1D::numx
size_t numx() const
Definition: state.h:197
Electric_Field_1D::C3
valarray< complex< double > > C3
Definition: vlasov.h:76
Field1D::Dx
Field1D & Dx()
Definition: state.cpp:450
Thomas_Tridiagonal
bool Thomas_Tridiagonal(Array2D< double > &A, valarray< double > &d, valarray< double > &xk)
The tridiagonal solver for implicit collisions.
Definition: nmethods.cpp:216
Magnetic_Field_1D::B1
valarray< complex< double > > B1
Definition: vlasov.h:106
EMF1D::By
Field1D & By()
Definition: state.h:294
Magnetic_Field_1D::A3
complex< double > A3
Definition: vlasov.h:108
Electric_Field_1D::A100
complex< double > A100
Definition: vlasov.h:72
Faraday_1D::idx
complex< double > idx
Definition: vlasov.h:146
input.h
Input reader - Declarations.
Electric_Field_1D::MakeG00
void MakeG00(SHarmonic1D &f)
Definition: vlasov.cpp:828
Ampere_1D::Ampere_1D
Ampere_1D(double xmin, double xmax, size_t Nx)
Definition: vlasov.cpp:1379
Field1D
A 1D Field.
Definition: state.h:184
lib-array.h
Underlying data structures.
Electric_Field_1D::C2
Array2D< complex< double > > C2
Definition: vlasov.h:77
Ampere_1D::tmpB
Field1D tmpB
Definition: vlasov.h:163
Electric_Field_1D::H
SHarmonic1D H
Definition: vlasov.h:70
Electric_Field_1D::Electric_Field_1D
Electric_Field_1D(size_t Nl, size_t Nm, double pmin, double pmax, size_t Np, double xmin, double xmax, size_t Nx)
Definition: vlasov.cpp:90
Magnetic_Field_1D::FLM
SHarmonic1D FLM
Definition: vlasov.h:104
Spatial_Advection_1D::fd2
SHarmonic1D fd2
Definition: vlasov.h:31
Spatial_Advection_1D::A00
complex< double > A00
Definition: vlasov.h:32
Magnetic_Field_1D::A1
valarray< complex< double > > A1
Definition: vlasov.h:106
DistFunc1D::mass
double mass() const
Definition: state.h:400
Magnetic_Field_1D::implicit
void implicit(DistFunc1D &Din, const Field1D &FBx, const Field1D &FBy, const Field1D &FBz, double dt)
Definition: vlasov.cpp:973
Electric_Field_1D::Implicit_Ex
void Implicit_Ex(const DistFunc1D &Din, const Field1D &FEx, DistFunc1D &Dh)
Definition: vlasov.cpp:443
Electric_Field_1D::C100
complex< double > C100
Definition: vlasov.h:72
Ampere_1D::idx
complex< double > idx
Definition: vlasov.h:164
Magnetic_Field_1D::Magnetic_Field_1D
Magnetic_Field_1D(size_t Nl, size_t Nm, double pmin, double pmax, size_t Np, double xmin, double xmax, size_t Nx)
Definition: vlasov.cpp:849
Magnetic_Field_1D::A2
Array2D< complex< double > > A2
Definition: vlasov.h:107
SHarmonic1D::mxaxis
SHarmonic1D & mxaxis(const valarray< complex< double > > &shmulti)
Definition: state.cpp:126
Spatial_Advection_1D::A2
Array2D< complex< double > > A2
Definition: vlasov.h:33
Current_1D::Jz
Field1D Jz
Definition: vlasov.h:128
Spatial_Advection_1D::vr
valarray< complex< double > > vr
Definition: vlasov.h:34
Current_1D::es1d
void es1d(const DistFunc1D &Din, Field1D &FExh)
Definition: vlasov.cpp:71
Spatial_Advection_1D::Spatial_Advection_1D
Spatial_Advection_1D(size_t Nl, size_t Nm, double pmin, double pmax, size_t Np, double xmin, double xmax, size_t Nx)
Definition: vlasov.cpp:1199
Spatial_Advection_1D::A20
complex< double > A20
Definition: vlasov.h:32
Faraday_1D::operator()
void operator()(EMF1D &EMFin, EMF1D &EMFh)
Definition: vlasov.cpp:1342
Electric_Field_1D::Implicit_Ey_f1only
void Implicit_Ey_f1only(const DistFunc1D &Din, const Field1D &FEy, DistFunc1D &Dh)
Definition: vlasov.cpp:722
Electric_Field_1D::Implicit_Ez
void Implicit_Ez(const DistFunc1D &Din, const Field1D &FEz, DistFunc1D &Dh)
Definition: vlasov.cpp:556
Electric_Field_1D::A210
complex< double > A210
Definition: vlasov.h:72
SHarmonic1D::Re
SHarmonic1D & Re()
Definition: state.cpp:130
vlasov.h
Functors for various time-integration methodds - Declarations.
SHarmonic1D::numx
size_t numx() const
Definition: state.h:72
EMF1D::Ez
Field1D & Ez()
Definition: state.h:292
Electric_Field_1D::A310
complex< double > A310
Definition: vlasov.h:72
SHarmonic1D
A 1D Spherical Harmonic.
Definition: state.h:57
SHarmonic1D::Dp
SHarmonic1D & Dp()
Definition: state.cpp:139
Field1D::array
valarray< complex< double > > & array() const
Definition: state.h:196
Electric_Field_1D::C311
complex< double > C311
Definition: vlasov.h:72
Magnetic_Field_1D::operator()
void operator()(const DistFunc1D &Din, const Field1D &FBx, const Field1D &FBy, const Field1D &FBz, DistFunc1D &Dh)
Definition: vlasov.cpp:898
Faraday_1D::numx
size_t numx
Definition: vlasov.h:147
Electric_Field_1D::B211
complex< double > B211
Definition: vlasov.h:72
Electric_Field_1D::f1only
void f1only(const DistFunc1D &Din, const Field1D &FEx, const Field1D &FEy, const Field1D &FEz, DistFunc1D &Dh)
Definition: vlasov.cpp:633
Magnetic_Field_1D::f1only
void f1only(const DistFunc1D &Din, const Field1D &FBx, const Field1D &FBy, const Field1D &FBz, DistFunc1D &Dh)
Definition: vlasov.cpp:1120
Spatial_Advection_1D::f1only
void f1only(const DistFunc1D &Din, DistFunc1D &Dh)
Definition: vlasov.cpp:1306
Electric_Field_1D::Implicit_Ey
void Implicit_Ey(const DistFunc1D &Din, const Field1D &FEy, DistFunc1D &Dh)
Definition: vlasov.cpp:488
Current_1D::Jy
Field1D Jy
Definition: vlasov.h:128
Current_1D::small
double small
Definition: vlasov.h:129
Spatial_Advection_1D::A1
Array2D< complex< double > > A1
Definition: vlasov.h:33
Faraday_1D::tmpE
Field1D tmpE
Definition: vlasov.h:145
DistFunc1D::m0
size_t m0() const
Definition: state.h:397
DistFunc1D::q
double q() const
Definition: state.h:399
Electric_Field_1D::A1
Array2D< complex< double > > A1
Definition: vlasov.h:74
EMF1D::Bz
Field1D & Bz()
Definition: state.h:295
Current_1D::operator()
void operator()(const DistFunc1D &Din, Field1D &FExh, Field1D &FEyh, Field1D &FEzh)
Definition: vlasov.cpp:53
Electric_Field_1D::B2
valarray< complex< double > > B2
Definition: vlasov.h:75
Ampere_1D::operator()
void operator()(EMF1D &EMFin, EMF1D &EMFh)
Definition: vlasov.cpp:1391
state.h
Fields, Distributions, Harmonics, States - Declarations.
Electric_Field_1D::pr
valarray< complex< double > > pr
Definition: vlasov.h:78
EMF1D
Definition: state.h:271
Spatial_Advection_1D::fd1
SHarmonic1D fd1
Definition: vlasov.h:31
Electric_Field_1D::TMP
SHarmonic1D TMP
Definition: vlasov.h:70
Electric_Field_1D::MakeGH
void MakeGH(SHarmonic1D &f, size_t l)
Definition: vlasov.cpp:807
Algorithms::MakeCAxis
valarray< T > MakeCAxis(const T min, const T max, const size_t N)
Definition: lib-algorithms.h:36
SHarmonic1D::Dx
SHarmonic1D & Dx()
Definition: state.cpp:185
Spatial_Advection_1D::operator()
void operator()(const DistFunc1D &Din, DistFunc1D &Dh)
Definition: vlasov.cpp:1249
Electric_Field_1D::G
SHarmonic1D G
Definition: vlasov.h:70
Electric_Field_1D::C1
valarray< complex< double > > C1
Definition: vlasov.h:76
EMF1D::Ey
Field1D & Ey()
Definition: state.h:291
Electric_Field_1D::B1
valarray< complex< double > > B1
Definition: vlasov.h:75
Spatial_Advection_1D::A10
complex< double > A10
Definition: vlasov.h:32
Electric_Field_1D::C4
Array2D< complex< double > > C4
Definition: vlasov.h:77
Current_1D::Current_1D
Current_1D(double pmin, double pmax, size_t Np, size_t Nx)
Definition: vlasov.cpp:38
Faraday_1D::Faraday_1D
Faraday_1D(double xmin, double xmax, size_t Nx)
Definition: vlasov.cpp:1330
DistFunc1D::l0
size_t l0() const
Definition: state.h:396
Electric_Field_1D::operator()
void operator()(const DistFunc1D &Din, const Field1D &FEx, const Field1D &FEy, const Field1D &FEz, DistFunc1D &Dh)
Definition: vlasov.cpp:187
Electric_Field_1D::es1d
void es1d(const DistFunc1D &Din, const Field1D &FEx, const Field1D &FEy, const Field1D &FEz, DistFunc1D &Dh)
Definition: vlasov.cpp:379
Electric_Field_1D::invpr
valarray< complex< double > > invpr
Definition: vlasov.h:78
Electric_Field_1D::Hp0
valarray< complex< double > > Hp0
Definition: vlasov.h:78
DistFunc1D::getcurrent
valarray< double > getcurrent(size_t dir)
Definition: state.cpp:965
Electric_Field_1D::A2
Array2D< complex< double > > A2
Definition: vlasov.h:74
DistFunc1D
The 1D distribution function is the container for all SHarmonic1D per species.
Definition: state.h:376
Electric_Field_1D::Implicit_Ex_f1only
void Implicit_Ex_f1only(const DistFunc1D &Din, const Field1D &FEx, DistFunc1D &Dh)
Definition: vlasov.cpp:693
Spatial_Advection_1D::es1d
void es1d(const DistFunc1D &Din, DistFunc1D &Dh)
Definition: vlasov.cpp:1280
Current_1D::Jx
Field1D Jx
Definition: vlasov.h:128
SHarmonic1D::mpaxis
SHarmonic1D & mpaxis(const valarray< complex< double > > &shmulti)
Definition: state.cpp:122
nmethods.h
Numerical Methods - Declarations.