field.pdf

#include <math.h>

#include <time.h>

#include "nrutil.h"

#include "someutil.h"

#include "glob_streamer.h"

static void ProcBoundary (void); // calculate the potential at boundary points

// [Liu and Pasko, J. Phys. D: Appl. Phys., 2006]

static void ProcSource (void); // form the source term for SOR iteration

// debug the output time

// endTime = clock();

// elapsed = (endTime-startTime)/CLOCKS_PER_SEC;

// printf("Execution Time for ProcSource: %f",elapsed);

// id determines whether update boundary condition or not

if (!id) {

//if(currStep == imodV) {

ProcBoundary ();

if (currStep <= 20)

imodV += 1;

Sor();

CylindricalGradient(E.phi, pErnode, pEznode, gd.r, gd.z, nr, nz, 0, 0, pEr, pEz,

g.r, g.z, z0, Rsph);

for (iz = 1; iz <= nz; iz++)

if (g.r[ir]*g.r[ir]+(g.z[iz]-z0)*(g.z[iz]-z0)>Rsph*Rsph){

pErnode[iz][ir] = -pErnode[iz][ir];

pEnode [iz][ir] = sqrt (pErnode[iz][ir] * pErnode[iz][ir] +

nrerror ("Error in solvefield, pEnode is not numuber. \n");

else {

pEnode[iz][ir] = 0.;

pErnode[iz][ir] =0.;

pEznode[iz][ir] =0.;

// pEr[iz][ir] =0.;

// pEz[iz][ir] =0.;

}

}

}

static void ProcSource (void)

dnorms = 0;

for (iz = 2; iz <= nz - 1; iz++) {

deltaz2 = gd.z[iz] + gd.z[iz - 1];

for (ir = 1; ir <= nr - 1; ir++) {

if (g.r[ir]*g.r[ir]+(g.z[iz]-z0)*(g.z[iz]-z0)>Rsph*Rsph){

if (ir != 1) {

deltar2 = gd.r[ir] + gd.r[ir - 1];

2./deltaz2/gd.z[iz]-2./deltaz2/gd.z[iz-1];

} else

-4./gd.r[ir]/gd.r[ir]-2./deltaz2/gd.z[iz]-2./deltaz2/gd.z[iz-1];

nrerror ("Error in procboundary, pne is negative. \n");

nrerror ("Error in procboundary, pne is not a numuber. \n");

if (pnn[iz][ir]>=0)

pnn[iz][ir]=pnn[iz][ir];

else if (pnn[iz][ir]<0. )

nrerror ("Error in procboundary, pnn is negative. \n");

else

nrerror ("Error in procboundary, pnn is not a numuber. \n");

E.srs[iz][ir] = -(pnp[iz][ir] - pne[iz][ir] - pnn[iz][ir]) * (eleCharge

/ epsilon0);

if (E.srs[iz][ir]>=0.||E.srs[iz][ir]<0. )

E.srs[iz][ir] = E.srs[iz][ir]/tmp;

dnorms += E.srs[iz][ir]*E.srs[iz][ir];