summaryrefslogtreecommitdiff
path: root/lpic/src/propagate_particles.C
diff options
context:
space:
mode:
Diffstat (limited to 'lpic/src/propagate_particles.C')
-rw-r--r--lpic/src/propagate_particles.C910
1 files changed, 910 insertions, 0 deletions
diff --git a/lpic/src/propagate_particles.C b/lpic/src/propagate_particles.C
new file mode 100644
index 0000000..c7f4755
--- /dev/null
+++ b/lpic/src/propagate_particles.C
@@ -0,0 +1,910 @@
+/*
+ This file is part of LPIC++, a particle-in-cell code for
+ simulating the interaction of laser light with plasma.
+
+ Copyright (C) 1994-1997 Roland Lichters
+
+ LPIC++ is free software; you can redistribute it and/or
+ modify it under the terms of the GNU General Public License
+ as published by the Free Software Foundation; either version 2
+ of the License, or (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+*/
+
+#include <propagate.h>
+
+void propagate::particles( domain &grid )
+// acceleration according Boris (in Birdsall, Langdon)
+{
+ static error_handler bob("propagate::particles",errname);
+
+ struct cell *cell;
+ struct particle *part;
+
+ // assumes fields of the following domain in cell rbuf
+
+ for( cell=grid.left; cell!=grid.rbuf; cell=cell->next ) // for all cells
+ {
+ if (cell->npart!=0)
+ {
+ part=cell->first;
+ do
+ {
+#ifdef DEBUG
+ if (!part) bob.error("segmentation");
+ if ((part->x < cell->x - TINY) || (part->x >= cell->x + dx + TINY) ) { // particle in cell ? --
+ bob.message( "particle at", part->x );
+ bob.message( " in cell", cell->number, "at", cell->x, " .. ", cell->x + dx );
+ bob.message( "is linked to wrong cell:" );
+ bob.message( " particle -> vx ", part->ux * part->igamma );
+ bob.message( "cell->x - part->x = ", cell->x - part->x );
+ bob.message( "cell->x + dx - part->x = ", cell->x + dx - part->x );
+ bob.error("");
+ }
+#endif
+
+ deposit_charge( cell, part ); // not necessary for the local algorithm
+ // charge distribution of the
+ // preceeding half time step
+ accelerate_1( cell, part );
+ }
+ while( (part=part->next) );
+
+ part=cell->first;
+ do part->igamma = 1.0/sqrt(part->igamma);
+ while( (part=part->next) );
+
+ part=cell->first;
+ do accelerate_2( cell, part );
+ while( (part=part->next) );
+
+ part=cell->first;
+ do part->igamma = 1.0/sqrt(part->igamma);
+ while( (part=part->next) );
+
+ part=cell->first;
+ do
+ {
+ move( part ); // move all particles
+
+ has_to_change_cell( cell, part ); // put particles on stack
+
+ deposit_current( cell, part ); // this step is necessary
+ }
+ while( (part=part->next) );
+ }
+ }
+
+ do_change_cell( grid ); // particles are removed from stack and linked to their
+ // new cells
+ // cells "Lbuf" and "Rbuf" remain empty
+
+ mask_current( grid ); // makes currents invisible near the box boundaries
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+inline void propagate::mask_current( domain &grid )
+// this routine sets the currents smoothly equal to zero in
+// a region of size 2*MASK at the left and right boundary of the simulation box
+// this inhibits artificial radiation at the boundaries
+
+{
+ static error_handler bob("propagate::mask_current",errname);
+
+ struct cell *cell;
+ int i;
+
+ if (domain_number==1) {
+ for( i=1, cell=grid.left; i<=2*MASK; i++, cell=cell->next ) { // the first 2MASK cells
+ cell->jx *= mask(i);
+ cell->jy *= mask(i);
+ cell->jz *= mask(i);
+ }
+ }
+
+ if (domain_number==n_domains) {
+ for( i=1, cell=grid.right; i<=2*MASK; i++, cell=cell->prev ) { // the last 2MASK cells
+ cell->jx *= mask(i);
+ cell->jy *= mask(i);
+ cell->jz *= mask(i);
+ }
+ }
+}
+
+
+inline double propagate::mask( int i )
+{
+ if (i<MASK) return 0;
+ else {
+ if (i<=2*MASK ) return pow( sin(0.5*PI*(i-MASK)/MASK), 2 );
+ else return 1.0;
+ }
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+inline void propagate::accelerate( struct cell *cell, struct particle *part )
+ //
+ // this function is currently not used
+ //
+ // full acceleration according Boris (in Birdsall, Langdon)
+{
+ static error_handler bob("propagate::accelerate",errname);
+
+ register double zmpidt = part->zm * PI * dt;
+
+ register double w = weighting(cell,part);
+ register double notw = 1.0-w;
+ register double ex = w * cell->ex + notw * cell->next->ex; // interpolate fields
+ register double ey = w * cell->ey + notw * cell->next->ey; // to particle position
+ register double ez = w * cell->ez + notw * cell->next->ez;
+ register double by = w * cell->by + notw * cell->next->by;
+ register double bz = w * cell->bz + notw * cell->next->bz;
+
+ register double ux = part->ux + ex * zmpidt; // half acceleration
+ register double uy = part->uy + ey * zmpidt;
+ register double uz = part->uz + ez * zmpidt;
+
+ register double igamma = (1.0 / sqrt(1.0 + (ux*ux + uy*uy + uz*uz)));
+
+
+ register double ty = by * zmpidt * igamma; // rotation
+ register double tz = bz * zmpidt * igamma;
+ register double t2 = ty*ty + tz*tz;
+ register double sy = 2.0 * ty / ( 1.0 + t2 );
+ register double sz = 2.0 * tz / ( 1.0 + t2 );
+
+ register double ux2 = ux * (1.0-tz*sz-ty*sy) + uy * sz - uz * sy;
+ register double uy2 = - ux * sz + uy * (1.0-tz*sz) + uz * ty*sz;
+ register double uz2 = ux * sy + uy * tz*sy + uz * (1.0-ty*sy);
+
+ part->ux = ux = ux2 + ex * zmpidt; // second half acceleration
+ part->uy = uy = uy2 + ey * zmpidt;
+ part->uz = uz = uz2 + ez * zmpidt;
+
+ part->igamma = (1.0 / sqrt(1.0 + (ux*ux + uy*uy + uz*uz)));
+
+ // we take the square roots for all particles per cell in a separate loop !
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+inline void propagate::accelerate_1( struct cell *cell, struct particle *part )
+// acceleration according Boris (in Birdsall, Langdon)
+{
+ static error_handler bob("propagate::accelerate_1",errname);
+
+ register double zmpidt = part->zm * PI * dt;
+
+ register double w = weighting(cell,part);
+ register double notw = 1.0-w;
+ // register double w0 = weighting_0(cell,part);
+ // register double notw0 = 1.0-w0;
+ register double ex = w * cell->ex + notw * cell->next->ex; // interpolate fields
+ register double ey = w * cell->ey + notw * cell->next->ey; // to particle position
+ register double ez = w * cell->ez + notw * cell->next->ez;
+
+ register double ux = part->ux + ex * zmpidt; // half acceleration
+ register double uy = part->uy + ey * zmpidt;
+ register double uz = part->uz + ez * zmpidt;
+
+ part->ux = ux;
+ part->uy = uy;
+ part->uz = uz;
+
+ part->igamma = 1.0 + ux*ux + uy*uy + uz*uz;
+
+ // we take the inverse square roots for all particles per cell in a separate loop !
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+inline void propagate::accelerate_2( struct cell *cell, struct particle *part )
+// acceleration according Boris (in Birdsall, Langdon)
+{
+ static error_handler bob("propagate::accelerate_2",errname);
+
+ register double zmpidt = part->zm * PI * dt;
+ register double igamma = part->igamma;
+
+ register double ux = part->ux;
+ register double uy = part->uy;
+ register double uz = part->uz;
+
+ register double w = weighting(cell,part);
+ register double notw = 1.0-w;
+ // register double w0 = weighting_0(cell,part);
+ // register double notw0 = 1.0-w0;
+ register double ex = w * cell->ex + notw * cell->next->ex; // interpolate fields
+ register double ey = w * cell->ey + notw * cell->next->ey; // to particle position
+ register double ez = w * cell->ez + notw * cell->next->ez;
+ register double by = w * cell->by + notw * cell->next->by;
+ register double bz = w * cell->bz + notw * cell->next->bz;
+
+ register double ty = by * zmpidt * igamma; // rotation
+ register double tz = bz * zmpidt * igamma;
+ register double t2 = ty*ty + tz*tz;
+ register double sy = 2.0 * ty / ( 1.0 + t2 );
+ register double sz = 2.0 * tz / ( 1.0 + t2 );
+
+ register double ux2 = ux * (1.0-tz*sz-ty*sy) + uy * sz - uz * sy;
+ register double uy2 = - ux * sz + uy * (1.0-tz*sz) + uz * ty*sz;
+ register double uz2 = ux * sy + uy * tz*sy + uz * (1.0-ty*sy);
+
+ part->ux = ux = ux2 + ex * zmpidt; // second half acceleration
+ part->uy = uy = uy2 + ey * zmpidt;
+ part->uz = uz = uz2 + ez * zmpidt;
+
+ part->igamma = 1.0 + ux*ux + uy*uy + uz*uz;
+
+ // we take the inverse square roots for all particles per cell in a separate loop !
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+inline void propagate::move( struct particle *part )
+{
+ static error_handler bob("propagate::move",errname);
+
+ if ( part->fix==1 ) part->dx = 0;
+ else {
+
+ part->dx = dx * part->ux * part->igamma; // dx = Gamma * dt, see constructor!
+ part->x += part->dx;
+
+#ifdef DEBUG
+ if ( fabs(part->dx) > dx )
+ bob.error( "particle displacement larger than grid spacing!" );
+#endif
+
+ }
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+inline void propagate::has_to_change_cell( struct cell *cell, struct particle *part )
+{
+ static error_handler bob("propagate::has_to_change_cell",errname);
+
+ if ( part->x < cell->x ) stk.put_on_stack( cell->prev, part );
+ else if ( part->x >= cell->x + dx ) stk.put_on_stack( cell->next, part );
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+inline void propagate::do_change_cell( domain &grid )
+{
+ static error_handler bob("propagate::do_change_cell",errname);
+
+ stack_member *stack, *stack_delete;
+ struct cell *cell;
+ struct cell *new_cell;
+ struct particle *part;
+
+ for( cell=grid.Lbuf; cell!=grid.dummy; cell=cell->next )
+ { // keep in mind the pointer to the first particle
+ cell->insert = cell->first; // in cell before inserting additional particles
+ } // from stack
+
+ stack = stk.zero->next;
+ while( stack != stk.hole ) // now insert particles from stack into new cells
+ { // and delete them from stack
+ part = stack->part;
+ new_cell = stack->new_cell;
+
+ stk.insert_particle( new_cell, part );
+
+ stack_delete = stack;
+
+ stk.remove_from_stack( stack_delete );
+
+ stack = stk.zero->next;
+ }
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+void propagate::reflect_particles( domain &grid )
+{
+ static error_handler bob("propagate::reflect_particles",errname);
+
+ struct particle *part;
+
+ // stack is empty after particles() has been called
+
+ if ( domain_number == 1 ) {
+
+ for( part=grid.lbuf->first; part!=NULL; part=part->next ) {
+ part->ux = - part->ux;
+ part->x -= part->dx;
+ bob.message( "re l", part->number, "time", time );
+
+ stk.put_on_stack( grid.left, part );
+ }
+
+ do_change_cell( grid );
+ }
+
+ if ( domain_number == n_domains ) {
+
+ for( part=grid.rbuf->first; part!=NULL; part=part->next ) {
+ part->ux = - part->ux;
+ part->x -= part->dx;
+ bob.message( "re r", part->number, "time", time );
+
+ stk.put_on_stack( grid.right, part );
+ }
+
+ do_change_cell( grid );
+ }
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+inline void propagate::deposit_charge( struct cell *cell, struct particle *part )
+{
+ double here, next, prev; // contributions to charge of this cell, ...
+ register double dist = (part->x - cell->x) * idx;
+
+ if ( dist <= 0.5 ) {
+ prev = part->zn * ( 0.5 - dist );
+ here = part->zn * ( 0.5 + dist );
+ cell->prev->charge += prev;
+ cell->prev->dens[part->species] += prev;
+ cell->charge += here;
+ cell->dens[part->species] += here;
+ }
+ else {
+ here = part->zn * ( 1.5 - dist );
+ next = part->zn * ( -0.5 + dist );
+ cell->charge += here;
+ cell->dens[part->species] += here;
+ cell->next->charge += next;
+ cell->next->dens[part->species] += next;
+ }
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+inline void propagate::deposit_current( struct cell *cell, struct particle *part )
+// We distinguish six cases:
+// first, distinguish former position in the first or second half of the cell
+// second, distinguish one one-boundary move and two two-boundary moves
+// #-boundary move means contributions to # boundary currents Jx
+//
+// the currents are calculated from the continuity equation,
+// assuming rectangular particle shape and area weighting
+// J.Villasenor and O.Buneman, Comp. Phys. Comm. 69 (1992) 306-316
+{
+ static error_handler bob("propagate::deposit_current",errname);
+
+ register double xm = part->x - part->dx; // before move
+ register double xp = part->x; // afterwards
+ register double x0 = part->cell->x; // former left hand cell boundary
+
+ register double x0p05dx = x0 + 0.5*dx;
+ register double x0m05dx = x0 - 0.5*dx;
+ register double x0p15dx = x0 + 1.5*dx;
+
+ if ( xm < x0p05dx ) { // former position in first half of the cell
+
+ if ( xp < x0m05dx ) left_two_left( cell, part ); // two boundary move to the left
+
+ else {
+ if ( xp >= x0p05dx ) left_two_right( cell, part ); // two-boundary move to the right
+ else left_one ( cell, part ); // one boundary move
+ }
+ }
+
+ else { // former position in the second half of the cell
+
+ if ( xp > x0p15dx ) right_two_right( cell, part );// two boundary move to the right
+
+ else {
+ if ( xp <= x0p05dx ) right_two_left( cell, part ); // two boundary move to the left
+ else right_one( cell, part ); // one boundary move
+ }
+ }
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+inline void propagate::left_one( struct cell *cell, struct particle *part )
+{
+ static error_handler bob("propagate::left_one",errname);
+
+ register double xm = part->x - part->dx; // before move
+ register double xp = part->x; // afterwards
+ register double x0 = part->cell->x; // former left hand cell boundary
+
+ register double jx0 = part->zn * (xp-xm)*idx;
+ /*
+ register double jx = cell->jx;
+ register double jy = cell->jy;
+ register double jz = cell->jz;
+ register double pjy = cell->prev->jy;
+ register double pjz = cell->prev->jz;
+ */
+
+ register double xpart = (xp+xm-2.0*x0)*idx;
+
+ register double r_1 = 0.5 * part->zn * ( 1.0 - xpart ) * part->igamma;
+// register double r0 = part->zn * part->igamma - r_1;
+ register double r0 = 0.5 * part->zn * ( 1.0 + xpart ) * part->igamma;
+
+ register double jy_1 = r_1 * part->uy;
+ register double jy0 = r0 * part->uy;
+ register double jz_1 = r_1 * part->uz;
+ register double jz0 = r0 * part->uz;
+ /*
+ cell->jx = jx + jx0;
+ cell->jy = jy + jy0;
+ cell->jz = jz + jz0;
+ cell->prev->jy = pjy + jy_1;
+ cell->prev->jz = pjz + jz_1;
+ */
+ cell->jx += jx0;
+ cell->jy += jy0;
+ cell->jz += jz0;
+ cell->prev->jy += jy_1;
+ cell->prev->jz += jz_1;
+
+#ifdef DEBUG
+// r_1 /= part->igamma;
+// r0 /= part->igamma;
+
+ if ( fabs(r_1+r0 - (part->zn * part->igamma)) > TINY || fabs(r_1) > fabs(part->zn * part->igamma) + TINY ||
+ fabs(r0) > fabs(part->zn * part->igamma) + TINY ) {
+ bob.message( "r_1 =", r_1 );
+ bob.message( "r0 =", r0 );
+ bob.message( "part->zn =", part->zn );
+ bob.message( "dif =", fabs(r_1+r0 - (part->zn * part->igamma)) );
+ bob.message( "part->N =", part->number );
+ bob.message( "part->igamma =", part->igamma );
+ bob.message( "TINY =", TINY );
+ bob.message( "xpart =", xpart );
+ bob.message( "xp =", xp );
+ bob.message( "xm =", xm );
+ bob.message( "x0 =", x0 );
+ bob.message( "idx =", idx );
+ bob.error( "density splitting wrong!" );
+ }
+#endif
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+inline void propagate::left_two_left( struct cell *cell, struct particle *part )
+{
+ static error_handler bob("propagate::left_two_left",errname);
+
+ register double xm = part->x - part->dx; // before move
+ register double xp = part->x; // afterwards
+ register double x0 = part->cell->x; // former left hand cell boundary
+ /*
+ register double jx = cell->jx;
+ register double jy = cell->jy;
+ register double jz = cell->jz;
+ register double pjx = cell->prev->jx;
+ register double pjy = cell->prev->jy;
+ register double pjz = cell->prev->jz;
+ register double ppjy = cell->prev->prev->jy;
+ register double ppjz = cell->prev->prev->jz;
+ */
+ register double jx_1 = - part->zn * ( 0.5 - (xp-x0+dx)*idx );
+ register double jx0 = - part->zn * ( 0.5 + (xm-x0)*idx );
+
+ register double eps = ( xm - (x0-0.5*dx) ) / ( xm - xp );
+ register double r_2 = part->zn * part->igamma * 0.5*(1.0-eps) * ( 0.5 - (xp-x0+dx)*idx );
+ register double r0 = part->zn * part->igamma * 0.5*eps * ( 0.5 + (xm-x0)*idx );
+ register double r_1 = part->zn * part->igamma - r0 - r_2;
+
+ register double jy_2 = r_2 * part->uy;
+ register double jy_1 = r_1 * part->uy;
+ register double jy0 = r0 * part->uy;
+ register double jz_2 = r_2 * part->uz;
+ register double jz_1 = r_1 * part->uz;
+ register double jz0 = r0 * part->uz;
+ /*
+ cell->jx = jx + jx0;
+ cell->jy = jy + jy0;
+ cell->jz = jz + jz0;
+ cell->prev->jx = pjx + jx_1;
+ cell->prev->jy = pjy + jy_1;
+ cell->prev->jz = pjz + jz_1;
+ cell->prev->prev->jy = ppjy + jy_2;
+ cell->prev->prev->jz = ppjz + jz_2;
+ */
+ cell->jx += jx0;
+ cell->jy += jy0;
+ cell->jz += jz0;
+ cell->prev->jx += jx_1;
+ cell->prev->jy += jy_1;
+ cell->prev->jz += jz_1;
+ cell->prev->prev->jy += jy_2;
+ cell->prev->prev->jz += jz_2;
+
+#ifdef DEBUG
+// r_2 /= part->igamma;
+// r_1 /= part->igamma;
+// r0 /= part->igamma;
+
+ if ( fabs(r_2+r_1+r0 - (part->zn*part->igamma)) > TINY || fabs(r_2) > fabs(part->zn*part->igamma) + TINY ||
+ fabs(r_1) > fabs(part->zn*part->igamma) + TINY || fabs(r0) > fabs(part->zn*part->igamma) + TINY ) {
+ bob.message( "r_2 =", r_2 );
+ bob.message( "r_1 =", r_1 );
+ bob.message( "r0 =", r0 );
+ bob.message( "dif =", fabs(r_2+r_1+r0 - (part->zn*part->igamma)) );
+ bob.message( "part->zn =", part->zn );
+ bob.message( "part->N =", part->number );
+ bob.message( "part->igamma =", part->igamma );
+ bob.error( "density splitting wrong!" );
+ }
+ if ( eps < 0 || eps > 1 )
+ bob.error( "eps!" );
+#endif
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+inline void propagate::left_two_right( struct cell *cell, struct particle *part )
+{
+ static error_handler bob("propagate::left_two_right",errname);
+
+ register double xm = part->x - part->dx; // before move
+ register double xp = part->x; // afterwards
+ register double x0 = part->cell->x; // former left hand cell boundary
+ /*
+ register double jx = cell->jx;
+ register double jy = cell->jy;
+ register double jz = cell->jz;
+ register double pjy = cell->prev->jy;
+ register double pjz = cell->prev->jz;
+ register double njx = cell->next->jx;
+ register double njy = cell->next->jy;
+ register double njz = cell->next->jz;
+ */
+ register double jx0 = part->zn * ( 0.5 - (xm-x0)*idx );
+ register double jx1 = part->zn * ( 0.5 + (xp-x0-dx)*idx );
+
+ register double eps = ( x0 + 0.5*dx - xm ) / ( xp - xm );
+ register double r_1 = 0.5*eps * part->zn * part->igamma * ( 0.5 - (xm-x0)*idx );
+ register double r1 = 0.5*(1.0-eps) * part->zn * part->igamma * ( 0.5 + (xp-x0-dx)*idx );
+ register double r0 = part->zn * part->igamma - r_1 - r1;
+
+ register double jy_1 = r_1 * part->uy;
+ register double jy0 = r0 * part->uy;
+ register double jy1 = r1 * part->uy;
+
+ register double jz_1 = r_1 * part->uz;
+ register double jz0 = r0 * part->uz;
+ register double jz1 = r1 * part->uz;
+ /*
+ cell->prev->jy = pjy + jy_1;
+ cell->prev->jz = pjz + jz_1;
+ cell->jx = jx + jx0;
+ cell->jy = jy + jy0;
+ cell->jz = jz + jz0;
+ cell->next->jx = njx + jx1;
+ cell->next->jy = njy + jy1;
+ cell->next->jz = njz + jz1;
+ */
+ cell->prev->jy += jy_1;
+ cell->prev->jz += jz_1;
+ cell->jx += jx0;
+ cell->jy += jy0;
+ cell->jz += jz0;
+ cell->next->jx += jx1;
+ cell->next->jy += jy1;
+ cell->next->jz += jz1;
+
+#ifdef DEBUG
+// r_1 /= part->igamma;
+// r0 /= part->igamma;
+// r1 /= part->igamma;
+
+ if ( fabs(r_1+r0+r1-(part->zn*part->igamma)) > TINY || fabs(r_1) > fabs(part->zn*part->igamma) + TINY ||
+ fabs(r0) > fabs(part->zn*part->igamma) + TINY || fabs(r1) > fabs(part->zn*part->igamma) + TINY ) {
+ bob.message( "r_1 =", r_1 );
+ bob.message( "r0 =", r0 );
+ bob.message( "r1 =", r1 );
+ bob.message( "part->zn =", part->zn );
+ bob.message( "dif =", fabs(r_1+r0+r1-(part->zn*part->igamma)) );
+ bob.message( "part->N =", part->number );
+ bob.message( "part->igamma =", part->igamma );
+ bob.error( "density splitting wrong!" );
+ }
+ if ( eps < 0 || eps > 1 )
+ bob.error( "eps!" );
+#endif
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+inline void propagate::right_one( struct cell *cell, struct particle *part )
+{
+ static error_handler bob("propagate::right_one",errname);
+
+ register double xm = part->x - part->dx; // before move
+ register double xp = part->x; // afterwards
+ register double x0 = part->cell->x; // former left hand cell boundary
+ /*
+ register double jy = cell->jy;
+ register double jz = cell->jz;
+ register double njx = cell->next->jx;
+ register double njy = cell->next->jy;
+ register double njz = cell->next->jz;
+ */
+ register double jx1 = part->zn * (xp-xm)*idx;
+
+ register double xpart = (xp+xm-2.0*(x0+dx))*idx;
+
+ register double r0 = 0.5 * part->zn * part->igamma * ( 1.0 - xpart );
+// register double r1 = part->zn * part->igamma - r0;
+ register double r1 = 0.5 * part->zn * part->igamma * ( 1.0 + xpart );
+
+ register double jy0 = r0 * part->uy;
+ register double jy1 = r1 * part->uy;
+
+ register double jz0 = r0 * part->uz;
+ register double jz1 = r1 * part->uz;
+ /*
+ cell->jy = jy + jy0;
+ cell->jz = jz + jz0;
+ cell->next->jx = njx + jx1;
+ cell->next->jy = njy + jy1;
+ cell->next->jz = njz + jz1;
+ */
+
+ cell->jy += jy0;
+ cell->jz += jz0;
+ cell->next->jx += jx1;
+ cell->next->jy += jy1;
+ cell->next->jz += jz1;
+
+#ifdef DEBUG
+// r0 /= part->igamma;
+// r1 /= part->igamma;
+
+ if ( fabs(r0+r1 - (part->zn * part->igamma)) > TINY || fabs(r0) > fabs(part->zn * part->igamma) + TINY ||
+ fabs(r1) > fabs(part->zn * part->igamma) + TINY ) {
+ bob.message( "r0 =", r0 );
+ bob.message( "r1 =", r1 );
+ bob.message( "part->zn =", part->zn );
+ bob.message( "dif =", fabs(r0+r1 - (part->zn * part->igamma)) );
+ bob.message( "part->N =", part->number );
+ bob.message( "part->igamma =", part->igamma );
+ bob.message( "TINY =", TINY );
+ bob.message( "xpart =", xpart );
+ bob.message( "xp =", xp );
+ bob.message( "xm =", xm );
+ bob.message( "x0 =", x0 );
+ bob.message( "idx =", idx );
+ bob.error( "density splitting wrong!" );
+ }
+#endif
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+inline void propagate::right_two_right( struct cell *cell, struct particle *part )
+{
+ static error_handler bob("propagate::right_two_right",errname);
+
+ register double xm = part->x - part->dx; // before move
+ register double xp = part->x; // afterwards
+ register double x0 = part->cell->x; // former left hand cell boundary
+ /*
+ register double jy = cell->jy;
+ register double jz = cell->jz;
+ register double njx = cell->next->jx;
+ register double njy = cell->next->jy;
+ register double njz = cell->next->jz;
+ register double nnjx = cell->next->next->jx;
+ register double nnjy = cell->next->next->jy;
+ register double nnjz = cell->next->next->jz;
+ */
+ register double jx1 = part->zn * ( 0.5 - (xm-x0-dx)*idx );
+ register double jx2 = part->zn * ( 0.5 + (xp-x0-2.0*dx)*idx );
+
+ register double eps = (x0+1.5*dx - xm) / (xp - xm);
+ register double r0 = 0.5*eps * part->zn * part->igamma * ( 0.5 - (xm-x0-dx)*idx );
+ register double r2 = 0.5*(1.0-eps) * part->zn * part->igamma * ( 0.5 + (xp-x0-2.0*dx)*idx );
+ register double r1 = part->zn * part->igamma - r0 - r2;
+
+ register double jy0 = r0 * part->uy;
+ register double jy1 = r1 * part->uy;
+ register double jy2 = r2 * part->uy;
+
+ register double jz0 = r0 * part->uz;
+ register double jz1 = r1 * part->uz;
+ register double jz2 = r2 * part->uz;
+ /*
+ cell->jy = jy + jy0;
+ cell->jz = jz + jz0;
+ cell->next->jx = njx + jx1;
+ cell->next->jy = njy + jy1;
+ cell->next->jz = njz + jz1;
+ cell->next->next->jx = nnjx + jx2;
+ cell->next->next->jy = nnjy + jy2;
+ cell->next->next->jz = nnjz + jz2;
+ */
+
+ cell->jy += jy0;
+ cell->jz += jz0;
+ cell->next->jx += jx1;
+ cell->next->jy += jy1;
+ cell->next->jz += jz1;
+ cell->next->next->jx += jx2;
+ cell->next->next->jy += jy2;
+ cell->next->next->jz += jz2;
+
+#ifdef DEBUG
+// r0 /= part->igamma;
+// r1 /= part->igamma;
+// r2 /= part->igamma;
+
+ if ( fabs(r0+r1+r2 - (part->zn*part->igamma)) > TINY || fabs(r0) > fabs(part->zn*part->igamma) + TINY ||
+ fabs(r1) > fabs(part->zn*part->igamma) + TINY || fabs(r2) > fabs(part->zn*part->igamma) + TINY ) {
+ bob.message( "r0 =", r0 );
+ bob.message( "r1 =", r1 );
+ bob.message( "r2 =", r2 );
+ bob.message( "part->zn =", part->zn );
+ bob.message( "dif =", fabs(r0+r1+r2 - (part->zn*part->igamma)) );
+ bob.message( "part->N =", part->number );
+ bob.message( "part->igamma =", part->igamma );
+ bob.error( "density splitting wrong!" );
+ }
+ if ( eps < 0 || eps > 1 )
+ bob.error( "eps!" );
+#endif
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+inline void propagate::right_two_left( struct cell *cell, struct particle *part )
+{
+ static error_handler bob("propagate::right_two_left",errname);
+
+ register double xm = part->x - part->dx; // before move
+ register double xp = part->x; // afterwards
+ register double x0 = part->cell->x; // former left hand cell boundary
+ /*
+ register double pjy = cell->prev->jy;
+ register double pjz = cell->prev->jz;
+ register double jx = cell->jx;
+ register double jy = cell->jy;
+ register double jz = cell->jz;
+ register double njx = cell->next->jx;
+ register double njy = cell->next->jy;
+ register double njz = cell->next->jz;
+ */
+ register double jx0 = - part->zn * ( 0.5 - (xp-x0)*idx );
+ register double jx1 = - part->zn * ( 0.5 + (xm-x0-dx)*idx );
+
+ register double eps = ( xm - (x0+0.5*dx) ) / ( xm - xp );
+ register double r_1 = 0.5*(1.0-eps) * part->zn * part->igamma * ( 0.5 - (xp-x0)*idx );
+ register double r1 = 0.5*eps * part->zn * part->igamma * ( 0.5 + (xm-x0-dx)*idx );
+ register double r0 = part->zn * part->igamma - r_1 - r1;
+
+ register double jy_1 = r_1 * part->uy;
+ register double jy0 = r0 * part->uy;
+ register double jy1 = r1 * part->uy;
+
+ register double jz_1 = r_1 * part->uz;
+ register double jz0 = r0 * part->uz;
+ register double jz1 = r1 * part->uz;
+ /*
+ cell->prev->jy = pjy + jy_1;
+ cell->prev->jz = pjz + jz_1;
+ cell->jx = jx + jx0;
+ cell->jy = jy + jy0;
+ cell->jz = jz + jz0;
+ cell->next->jx = njx + jx1;
+ cell->next->jy = njy + jy1;
+ cell->next->jz = njz + jz1;
+ */
+
+ cell->prev->jy += jy_1;
+ cell->prev->jz += jz_1;
+ cell->jx += jx0;
+ cell->jy += jy0;
+ cell->jz += jz0;
+ cell->next->jx += jx1;
+ cell->next->jy += jy1;
+ cell->next->jz += jz1;
+
+#ifdef DEBUG
+// r_1 /= part->igamma;
+// r0 /= part->igamma;
+// r1 /= part->igamma;
+
+ if ( fabs(r_1+r0+r1 - (part->zn*part->igamma)) > TINY || fabs(r_1) > fabs(part->zn*part->igamma) + TINY ||
+ fabs(r0) > fabs(part->zn*part->igamma) + TINY || fabs(r1) > fabs(part->zn*part->igamma) + TINY ) {
+ bob.message( "r_1 =", r_1 );
+ bob.message( "r0 =", r0 );
+ bob.message( "r1 =", r1 );
+ bob.message( "part->zn =", part->zn );
+ bob.message( "dif =", fabs(r_1+r0+r1 - (part->zn*part->igamma)) );
+ bob.message( "part->N =", part->number );
+ bob.message( "part->igamma =", part->igamma );
+ bob.error( "density splitting wrong!" );
+ }
+
+ if ( eps < 0 || eps > 1 )
+ bob.error( "eps!" );
+#endif
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+inline double propagate::weighting( struct cell *cell, struct particle *part )
+//
+// returns contribution to the grid point left of the particle's position
+// (1 - return value) is the contribution to the grid point on its right
+//
+{
+// zero order weighting
+// return ( 1.0 - floor( (part->x - cell->x) * idx + 0.5 ) );
+
+// linear weighting
+ return ( 1.0 - (part->x - cell->x) * idx );
+}
+
+inline double propagate::weighting_0( struct cell *cell, struct particle *part )
+//
+// returns contribution to the grid point left of the particle's position
+// (1 - return value) is the contribution to the grid point on its right
+//
+{
+ // zero order weighting
+ return ( 1.0 - floor( (part->x - cell->x) * idx + 0.5 ) );
+
+ // linear weighting
+ // return ( 1.0 - (part->x - cell->x) * idx );
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+//EOF
generated by cgit v1.2.3-54-g00ecf (git 2.44.0) at 2024-05-21 19:16:14 +0000