Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
fopr_CloverTerm_imp.cpp
Go to the documentation of this file.
1 
14 #include "fopr_CloverTerm.h"
15 
16 #ifdef USE_PARAMETERS_FACTORY
17 #include "parameters_factory.h"
18 #endif
19 
20 using std::valarray;
21 using std::string;
22 
23 //====================================================================
24 //- parameter entries
25 namespace {
26  void append_entry(Parameters& param)
27  {
28  param.Register_double("hopping_parameter", 0.0);
29  param.Register_double("clover_coefficient", 0.0);
30  param.Register_int_vector("boundary_condition", std::valarray<int>());
31 
32  param.Register_string("verbose_level", "NULL");
33  }
34 
35 
36 #ifdef USE_PARAMETERS_FACTORY
37  bool init_param = ParametersFactory::Register("", append_entry);
38 #endif
39 }
40 //- end
41 
42 //- parameters class
43 Parameters_Fopr_CloverTerm::Parameters_Fopr_CloverTerm() { append_entry(*this); }
44 //- end
45 
46 //====================================================================
48 namespace {
49  inline double mult_uv_r(double *u, double *v)
50  {
51  return u[0] * v[0] - u[1] * v[1]
52  + u[2] * v[2] - u[3] * v[3]
53  + u[4] * v[4] - u[5] * v[5];
54  }
55 
56 
57  inline double mult_uv_i(double *u, double *v)
58  {
59  return u[0] * v[1] + u[1] * v[0]
60  + u[2] * v[3] + u[3] * v[2]
61  + u[4] * v[5] + u[5] * v[4];
62  }
63 }
64 
65 //====================================================================
66 void Fopr_CloverTerm::set_parameters(const Parameters& params)
67 {
68  const string str_vlevel = params.get_string("verbose_level");
69 
70  m_vl = vout.set_verbose_level(str_vlevel);
71 
72  //- fetch and check input parameters
73  double kappa, cSW;
74  valarray<int> bc;
75 
76  int err = 0;
77  err += params.fetch_double("hopping_parameter", kappa);
78  err += params.fetch_double("clover_coefficient", cSW);
79  err += params.fetch_int_vector("boundary_condition", bc);
80 
81  if (err) {
82  vout.crucial(m_vl, "Fopr_CloverTerm: fetch error, input parameter not found.\n");
83  abort();
84  }
85 
86 
87  set_parameters(kappa, cSW, bc);
88 }
89 
90 
91 //====================================================================
92 void Fopr_CloverTerm::set_parameters(double kappa, double cSW,
93  valarray<int> bc)
94 {
95  //- print input parameters
96  vout.general(m_vl, "Parameters of Fopr_CloverTerm:\n");
97  vout.general(m_vl, " kappa = %8.4f\n", kappa);
98  vout.general(m_vl, " cSW = %8.4f\n", cSW);
99  for (int mu = 0; mu < m_Ndim; ++mu) {
100  vout.general(m_vl, " boundary[%d] = %2d\n", mu, bc[mu]);
101  }
102 
103  //- range check
104  // NB. kappa,cSW == 0 is allowed.
105  assert(bc.size() == m_Ndim);
106 
107  //- store values
108  m_kappa = kappa;
109  m_cSW = cSW;
110 
111  // m_boundary.resize(m_Ndim); // already resized in init.
112  for (int mu = 0; mu < m_Ndim; ++mu) {
113  m_boundary[mu] = bc[mu];
114  }
115 }
116 
117 
118 //====================================================================
119 void Fopr_CloverTerm::set_config(Field *U)
120 {
121  m_U = (Field_G *)U;
122  set_csw();
123 }
124 
125 
126 //====================================================================
127 void Fopr_CloverTerm::init(string repr)
128 {
129  m_Nvol = CommonParameters::Nvol();
130  m_Ndim = CommonParameters::Ndim();
131  m_Nc = CommonParameters::Nc();
132  m_Nd = CommonParameters::Nd();
133  m_NinF = 2 * m_Nc * m_Nd;
134 
135  m_U = 0;
136 
137  m_repr = repr;
138 
139  m_boundary.resize(m_Ndim);
140  m_SG.resize(m_Ndim * m_Ndim);
141 
142  GammaMatrixSet *gmset = GammaMatrixSet::New(m_repr);
143 
144  m_GM5 = gmset->get_GM(gmset->GAMMA5);
145 
146  m_SG[sg_index(0, 1)] = gmset->get_GM(gmset->SIGMA12);
147  m_SG[sg_index(1, 2)] = gmset->get_GM(gmset->SIGMA23);
148  m_SG[sg_index(2, 0)] = gmset->get_GM(gmset->SIGMA31);
149  m_SG[sg_index(3, 0)] = gmset->get_GM(gmset->SIGMA41);
150  m_SG[sg_index(3, 1)] = gmset->get_GM(gmset->SIGMA42);
151  m_SG[sg_index(3, 2)] = gmset->get_GM(gmset->SIGMA43);
152 
153  m_SG[sg_index(1, 0)] = m_SG[sg_index(0, 1)].mult(-1);
154  m_SG[sg_index(2, 1)] = m_SG[sg_index(1, 2)].mult(-1);
155  m_SG[sg_index(0, 2)] = m_SG[sg_index(2, 0)].mult(-1);
156  m_SG[sg_index(0, 3)] = m_SG[sg_index(3, 0)].mult(-1);
157  m_SG[sg_index(1, 3)] = m_SG[sg_index(3, 1)].mult(-1);
158  m_SG[sg_index(2, 3)] = m_SG[sg_index(3, 2)].mult(-1);
159 
160  m_SG[sg_index(0, 0)] = gmset->get_GM(gmset->UNITY);
161  m_SG[sg_index(1, 1)] = gmset->get_GM(gmset->UNITY);
162  m_SG[sg_index(2, 2)] = gmset->get_GM(gmset->UNITY);
163  m_SG[sg_index(3, 3)] = gmset->get_GM(gmset->UNITY);
164  // these 4 gamma matrices are actually not used.
165 
166 
167  // m_fopr_w = new Fopr_Wilson(repr);
168  if (m_repr == "Dirac") {
169  m_csw = &Fopr_CloverTerm::mult_csw_dirac;
170  m_gm5 = &Fopr_CloverTerm::gm5_dirac;
171  } else if (m_repr == "Chiral") {
172  m_csw = &Fopr_CloverTerm::mult_csw_chiral;
173  m_gm5 = &Fopr_CloverTerm::gm5_chiral;
174  }
175 
176  delete gmset;
177 }
178 
179 
180 //====================================================================
181 void Fopr_CloverTerm::tidyup()
182 {
183  // nothing to do.
184 }
185 
186 
187 //====================================================================
188 void Fopr_CloverTerm::mult_gm5(Field& v, const Field& f)
189 {
190  (this->*m_gm5)(v, f);
191 }
192 
193 
194 //====================================================================
195 void Fopr_CloverTerm::gm5_dirac(Field& w, const Field& f)
196 {
197  int Nvc = 2 * CommonParameters::Nc();
198  int Nd = CommonParameters::Nd();
199 
200  double *v1;
201  double *v2;
202 
203  v1 = const_cast<Field *>(&f)->ptr(0);
204  v2 = w.ptr(0);
205 
206  int id1 = 0;
207  int id2 = Nvc;
208  int id3 = Nvc * 2;
209  int id4 = Nvc * 3;
210 
211  for (int site = 0; site < m_Nvol; ++site) {
212  for (int icc = 0; icc < Nvc; icc++) {
213  int in = Nvc * Nd * site;
214  v2[icc + id1 + in] = v1[icc + id3 + in];
215  v2[icc + id2 + in] = v1[icc + id4 + in];
216  v2[icc + id3 + in] = v1[icc + id1 + in];
217  v2[icc + id4 + in] = v1[icc + id2 + in];
218  }
219  }
220 }
221 
222 
223 //====================================================================
224 void Fopr_CloverTerm::gm5_chiral(Field& w, const Field& f)
225 {
226  int Nvc = 2 * CommonParameters::Nc();
227  int Nd = CommonParameters::Nd();
228 
229  double *v1;
230  double *v2;
231 
232  v1 = const_cast<Field *>(&f)->ptr(0);
233  v2 = w.ptr(0);
234 
235  int id1 = 0;
236  int id2 = Nvc;
237  int id3 = Nvc * 2;
238  int id4 = Nvc * 3;
239 
240  for (int site = 0; site < m_Nvol; ++site) {
241  for (int icc = 0; icc < Nvc; icc++) {
242  int in = Nvc * Nd * site;
243  v2[icc + id1 + in] = v1[icc + id1 + in];
244  v2[icc + id2 + in] = v1[icc + id2 + in];
245  v2[icc + id3 + in] = -v1[icc + id3 + in];
246  v2[icc + id4 + in] = -v1[icc + id4 + in];
247  }
248  }
249 }
250 
251 
252 //====================================================================
253 void Fopr_CloverTerm::mult_isigma(Field_F& v, const Field_F& w,
254  const int mu, const int nu)
255 {
256  assert(mu != nu);
257  v.mult_iGM(m_SG[sg_index(mu, nu)], w);
258 }
259 
260 
261 //====================================================================
262 void Fopr_CloverTerm::mult_sigmaF(Field& v, const Field& f)
263 {
264  mult_csw(v, f);
265 }
266 
267 
268 //====================================================================
269 void Fopr_CloverTerm::mult_csw(Field& v, const Field& w)
270 {
271  (this->*m_csw)(v, w);
272 }
273 
274 
275 //====================================================================
276 void Fopr_CloverTerm::mult_csw_chiral(Field& v, const Field& w)
277 {
278  assert(w.nex() == 1);
279 
280  int Nc = CommonParameters::Nc();
281  int Nvc = 2 * Nc;
282  int Ndf = 2 * Nc * Nc;
283  int Nd = CommonParameters::Nd();
284  int Nvol = w.nvol();
285 
286  int id1 = 0;
287  int id2 = Nvc;
288  int id3 = Nvc * 2;
289  int id4 = Nvc * 3;
290 
291  double *w2 = const_cast<Field *>(&w)->ptr(0);
292  double *v2 = v.ptr(0);
293 
294  double *Bx = m_Bx.ptr(0);
295  double *By = m_By.ptr(0);
296  double *Bz = m_Bz.ptr(0);
297  double *Ex = m_Ex.ptr(0);
298  double *Ey = m_Ey.ptr(0);
299  double *Ez = m_Ez.ptr(0);
300 
301  v = 0.0;
302 
303  for (int site = 0; site < Nvol; ++site) {
304  int iv = Nvc * Nd * site;
305  int ig = Ndf * site;
306 
307  for (int ic = 0; ic < Nc; ++ic) {
308  int icr = 2 * ic;
309  int ici = icr + 1;
310  int icg = ic * Nvc + ig;
311 
312  // isigma_23 * Bx
313  v2[icr + id1 + iv] -= mult_uv_i(&Bx[icg], &w2[id2 + iv]);
314  v2[ici + id1 + iv] += mult_uv_r(&Bx[icg], &w2[id2 + iv]);
315  v2[icr + id2 + iv] -= mult_uv_i(&Bx[icg], &w2[id1 + iv]);
316  v2[ici + id2 + iv] += mult_uv_r(&Bx[icg], &w2[id1 + iv]);
317 
318  v2[icr + id3 + iv] -= mult_uv_i(&Bx[icg], &w2[id4 + iv]);
319  v2[ici + id3 + iv] += mult_uv_r(&Bx[icg], &w2[id4 + iv]);
320  v2[icr + id4 + iv] -= mult_uv_i(&Bx[icg], &w2[id3 + iv]);
321  v2[ici + id4 + iv] += mult_uv_r(&Bx[icg], &w2[id3 + iv]);
322 
323  // isigma_31 * By
324  v2[icr + id1 + iv] += mult_uv_r(&By[icg], &w2[id2 + iv]);
325  v2[ici + id1 + iv] += mult_uv_i(&By[icg], &w2[id2 + iv]);
326  v2[icr + id2 + iv] -= mult_uv_r(&By[icg], &w2[id1 + iv]);
327  v2[ici + id2 + iv] -= mult_uv_i(&By[icg], &w2[id1 + iv]);
328 
329  v2[icr + id3 + iv] += mult_uv_r(&By[icg], &w2[id4 + iv]);
330  v2[ici + id3 + iv] += mult_uv_i(&By[icg], &w2[id4 + iv]);
331  v2[icr + id4 + iv] -= mult_uv_r(&By[icg], &w2[id3 + iv]);
332  v2[ici + id4 + iv] -= mult_uv_i(&By[icg], &w2[id3 + iv]);
333 
334  // isigma_12 * Bz
335  v2[icr + id1 + iv] -= mult_uv_i(&Bz[icg], &w2[id1 + iv]);
336  v2[ici + id1 + iv] += mult_uv_r(&Bz[icg], &w2[id1 + iv]);
337  v2[icr + id2 + iv] += mult_uv_i(&Bz[icg], &w2[id2 + iv]);
338  v2[ici + id2 + iv] -= mult_uv_r(&Bz[icg], &w2[id2 + iv]);
339 
340  v2[icr + id3 + iv] -= mult_uv_i(&Bz[icg], &w2[id3 + iv]);
341  v2[ici + id3 + iv] += mult_uv_r(&Bz[icg], &w2[id3 + iv]);
342  v2[icr + id4 + iv] += mult_uv_i(&Bz[icg], &w2[id4 + iv]);
343  v2[ici + id4 + iv] -= mult_uv_r(&Bz[icg], &w2[id4 + iv]);
344 
345  // isigma_41 * Ex
346  v2[icr + id1 + iv] += mult_uv_i(&Ex[icg], &w2[id2 + iv]);
347  v2[ici + id1 + iv] -= mult_uv_r(&Ex[icg], &w2[id2 + iv]);
348  v2[icr + id2 + iv] += mult_uv_i(&Ex[icg], &w2[id1 + iv]);
349  v2[ici + id2 + iv] -= mult_uv_r(&Ex[icg], &w2[id1 + iv]);
350 
351  v2[icr + id3 + iv] -= mult_uv_i(&Ex[icg], &w2[id4 + iv]);
352  v2[ici + id3 + iv] += mult_uv_r(&Ex[icg], &w2[id4 + iv]);
353  v2[icr + id4 + iv] -= mult_uv_i(&Ex[icg], &w2[id3 + iv]);
354  v2[ici + id4 + iv] += mult_uv_r(&Ex[icg], &w2[id3 + iv]);
355 
356  // isigma_42 * Ey
357  v2[icr + id1 + iv] -= mult_uv_r(&Ey[icg], &w2[id2 + iv]);
358  v2[ici + id1 + iv] -= mult_uv_i(&Ey[icg], &w2[id2 + iv]);
359  v2[icr + id2 + iv] += mult_uv_r(&Ey[icg], &w2[id1 + iv]);
360  v2[ici + id2 + iv] += mult_uv_i(&Ey[icg], &w2[id1 + iv]);
361 
362  v2[icr + id3 + iv] += mult_uv_r(&Ey[icg], &w2[id4 + iv]);
363  v2[ici + id3 + iv] += mult_uv_i(&Ey[icg], &w2[id4 + iv]);
364  v2[icr + id4 + iv] -= mult_uv_r(&Ey[icg], &w2[id3 + iv]);
365  v2[ici + id4 + iv] -= mult_uv_i(&Ey[icg], &w2[id3 + iv]);
366 
367  // isigma_43 * Ez
368  v2[icr + id1 + iv] += mult_uv_i(&Ez[icg], &w2[id1 + iv]);
369  v2[ici + id1 + iv] -= mult_uv_r(&Ez[icg], &w2[id1 + iv]);
370  v2[icr + id2 + iv] -= mult_uv_i(&Ez[icg], &w2[id2 + iv]);
371  v2[ici + id2 + iv] += mult_uv_r(&Ez[icg], &w2[id2 + iv]);
372 
373  v2[icr + id3 + iv] -= mult_uv_i(&Ez[icg], &w2[id3 + iv]);
374  v2[ici + id3 + iv] += mult_uv_r(&Ez[icg], &w2[id3 + iv]);
375  v2[icr + id4 + iv] += mult_uv_i(&Ez[icg], &w2[id4 + iv]);
376  v2[ici + id4 + iv] -= mult_uv_r(&Ez[icg], &w2[id4 + iv]);
377  }
378  }
379 
380  v *= m_kappa * m_cSW;
381 }
382 
383 
384 //====================================================================
385 void Fopr_CloverTerm::mult_csw_dirac(Field& v, const Field& w)
386 {
387  assert(w.nex() == 1);
388 
389  int Nc = CommonParameters::Nc();
390  int Nvc = 2 * Nc;
391  int Ndf = 2 * Nc * Nc;
392  int Nd = CommonParameters::Nd();
393  int Nvol = w.nvol();
394 
395  int id1 = 0;
396  int id2 = Nvc;
397  int id3 = Nvc * 2;
398  int id4 = Nvc * 3;
399 
400  double *w2 = const_cast<Field *>(&w)->ptr(0);
401  double *v2 = v.ptr(0);
402 
403  double *Bx = m_Bx.ptr(0);
404  double *By = m_By.ptr(0);
405  double *Bz = m_Bz.ptr(0);
406  double *Ex = m_Ex.ptr(0);
407  double *Ey = m_Ey.ptr(0);
408  double *Ez = m_Ez.ptr(0);
409 
410  v = 0.0;
411 
412  for (int site = 0; site < Nvol; ++site) {
413  int iv = Nvc * Nd * site;
414  int ig = Ndf * site;
415 
416  for (int ic = 0; ic < Nc; ++ic) {
417  int icr = 2 * ic;
418  int ici = icr + 1;
419  int icg = ic * Nvc + ig;
420 
421  // isigma_23 * Bx
422  v2[icr + id1 + iv] -= mult_uv_i(&Bx[icg], &w2[id2 + iv]);
423  v2[ici + id1 + iv] += mult_uv_r(&Bx[icg], &w2[id2 + iv]);
424  v2[icr + id2 + iv] -= mult_uv_i(&Bx[icg], &w2[id1 + iv]);
425  v2[ici + id2 + iv] += mult_uv_r(&Bx[icg], &w2[id1 + iv]);
426 
427  v2[icr + id3 + iv] -= mult_uv_i(&Bx[icg], &w2[id4 + iv]);
428  v2[ici + id3 + iv] += mult_uv_r(&Bx[icg], &w2[id4 + iv]);
429  v2[icr + id4 + iv] -= mult_uv_i(&Bx[icg], &w2[id3 + iv]);
430  v2[ici + id4 + iv] += mult_uv_r(&Bx[icg], &w2[id3 + iv]);
431 
432  // isigma_31 * By
433  v2[icr + id1 + iv] += mult_uv_r(&By[icg], &w2[id2 + iv]);
434  v2[ici + id1 + iv] += mult_uv_i(&By[icg], &w2[id2 + iv]);
435  v2[icr + id2 + iv] -= mult_uv_r(&By[icg], &w2[id1 + iv]);
436  v2[ici + id2 + iv] -= mult_uv_i(&By[icg], &w2[id1 + iv]);
437 
438  v2[icr + id3 + iv] += mult_uv_r(&By[icg], &w2[id4 + iv]);
439  v2[ici + id3 + iv] += mult_uv_i(&By[icg], &w2[id4 + iv]);
440  v2[icr + id4 + iv] -= mult_uv_r(&By[icg], &w2[id3 + iv]);
441  v2[ici + id4 + iv] -= mult_uv_i(&By[icg], &w2[id3 + iv]);
442 
443  // isigma_12 * Bz
444  v2[icr + id1 + iv] -= mult_uv_i(&Bz[icg], &w2[id1 + iv]);
445  v2[ici + id1 + iv] += mult_uv_r(&Bz[icg], &w2[id1 + iv]);
446  v2[icr + id2 + iv] += mult_uv_i(&Bz[icg], &w2[id2 + iv]);
447  v2[ici + id2 + iv] -= mult_uv_r(&Bz[icg], &w2[id2 + iv]);
448 
449  v2[icr + id3 + iv] -= mult_uv_i(&Bz[icg], &w2[id3 + iv]);
450  v2[ici + id3 + iv] += mult_uv_r(&Bz[icg], &w2[id3 + iv]);
451  v2[icr + id4 + iv] += mult_uv_i(&Bz[icg], &w2[id4 + iv]);
452  v2[ici + id4 + iv] -= mult_uv_r(&Bz[icg], &w2[id4 + iv]);
453 
454  // isigma_41 * Ex
455  v2[icr + id1 + iv] += mult_uv_i(&Ex[icg], &w2[id4 + iv]);
456  v2[ici + id1 + iv] -= mult_uv_r(&Ex[icg], &w2[id4 + iv]);
457  v2[icr + id2 + iv] += mult_uv_i(&Ex[icg], &w2[id3 + iv]);
458  v2[ici + id2 + iv] -= mult_uv_r(&Ex[icg], &w2[id3 + iv]);
459 
460  v2[icr + id3 + iv] += mult_uv_i(&Ex[icg], &w2[id2 + iv]);
461  v2[ici + id3 + iv] -= mult_uv_r(&Ex[icg], &w2[id2 + iv]);
462  v2[icr + id4 + iv] += mult_uv_i(&Ex[icg], &w2[id1 + iv]);
463  v2[ici + id4 + iv] -= mult_uv_r(&Ex[icg], &w2[id1 + iv]);
464 
465  // isigma_42 * Ey
466  v2[icr + id1 + iv] -= mult_uv_r(&Ey[icg], &w2[id4 + iv]);
467  v2[ici + id1 + iv] -= mult_uv_i(&Ey[icg], &w2[id4 + iv]);
468  v2[icr + id2 + iv] += mult_uv_r(&Ey[icg], &w2[id3 + iv]);
469  v2[ici + id2 + iv] += mult_uv_i(&Ey[icg], &w2[id3 + iv]);
470 
471  v2[icr + id3 + iv] -= mult_uv_r(&Ey[icg], &w2[id2 + iv]);
472  v2[ici + id3 + iv] -= mult_uv_i(&Ey[icg], &w2[id2 + iv]);
473  v2[icr + id4 + iv] += mult_uv_r(&Ey[icg], &w2[id1 + iv]);
474  v2[ici + id4 + iv] += mult_uv_i(&Ey[icg], &w2[id1 + iv]);
475 
476  // isigma_43 * Ez
477  v2[icr + id1 + iv] += mult_uv_i(&Ez[icg], &w2[id3 + iv]);
478  v2[ici + id1 + iv] -= mult_uv_r(&Ez[icg], &w2[id3 + iv]);
479  v2[icr + id2 + iv] -= mult_uv_i(&Ez[icg], &w2[id4 + iv]);
480  v2[ici + id2 + iv] += mult_uv_r(&Ez[icg], &w2[id4 + iv]);
481 
482  v2[icr + id3 + iv] += mult_uv_i(&Ez[icg], &w2[id1 + iv]);
483  v2[ici + id3 + iv] -= mult_uv_r(&Ez[icg], &w2[id1 + iv]);
484  v2[icr + id4 + iv] -= mult_uv_i(&Ez[icg], &w2[id2 + iv]);
485  v2[ici + id4 + iv] += mult_uv_r(&Ez[icg], &w2[id2 + iv]);
486  }
487  }
488 
489  v *= m_kappa * m_cSW;
490 }
491 
492 
493 //====================================================================
494 void Fopr_CloverTerm::set_csw()
495 {
496  set_fieldstrength(m_Bx, 1, 2);
497  set_fieldstrength(m_By, 2, 0);
498  set_fieldstrength(m_Bz, 0, 1);
499  set_fieldstrength(m_Ex, 3, 0);
500  set_fieldstrength(m_Ey, 3, 1);
501  set_fieldstrength(m_Ez, 3, 2);
502 }
503 
504 
505 //====================================================================
506 void Fopr_CloverTerm::set_fieldstrength(Field_G& Fst,
507  const int mu, const int nu)
508 {
509  int Nvol = CommonParameters::Nvol();
510 
511  Staples staple;
512 
513  Field_G Cup(Nvol, 1), Cdn(Nvol, 1);
514  Field_G Umu(Nvol, 1);
515  Field_G v(Nvol, 1), v2(Nvol, 1);
516 
517  Cup = staple.upper(*m_U, mu, nu);
518  Cdn = staple.lower(*m_U, mu, nu);
519  Umu.setpart_ex(0, *m_U, mu);
520 
521  Fst.mult_Field_Gnd(0, Umu, 0, Cup, 0);
522  Fst.multadd_Field_Gnd(0, Umu, 0, Cdn, 0, -1.0);
523 
524  v.mult_Field_Gdn(0, Cup, 0, Umu, 0);
525  v.multadd_Field_Gdn(0, Cdn, 0, Umu, 0, -1.0);
526 
527  m_shift.forward(v2, v, mu);
528 
529  Fst += v2;
530 
531  Fst.ah_Field_G(0);
532  Fst *= 0.25;
533 }
534 
535 
536 //====================================================================
537 //============================================================END=====