/******************************************************************** * * * Plot potential from a Maxwellian electron density function. * * Experimentally based approximate electron distribution in * * polywell. See notes on "Electron Fluid in a Polywell" * * * * R ~ 1/(exp(-a*r*r) + b*r*r) * * * * Author = Mike Rosing * * date = Jan 23, 2008 * * * *******************************************************************/ #include #include #include //#define MAXSTEPS 100 #define MAXSTEPS 25 /* MAXWALL is the ratio of outer containment sphere radius to L (center of sphere to center of coil distance) */ #define MAXWALL (3.0) /* COILRADIUS is the ratio of coil radius to L */ #define COILRADIUS (0.9) #define NUMELMNTS ((MAXSTEPS+1)*(MAXSTEPS+2)*(MAXSTEPS+3)/6) /* Create global buffers for data storage so integrals are easy */ double GridPot[NUMELMNTS]; double ElecPot0[NUMELMNTS], ElecPot1[NUMELMNTS]; double dens[MAXSTEPS*(MAXSTEPS+1)/2], energy[NUMELMNTS]; double angles[2*48*NUMELMNTS]; int numangles[NUMELMNTS]; double test_integral[NUMELMNTS]; /* ELECTRON_TEMP is the temperature to MaGrid voltage ratio. It is really an energy ratio defined by ELECTRON_TEMP*e*V = k*T where V is the MaGrid voltage, k is Boltzmann's constant and T is the electron temperature. */ #define ELECTRON_TEMP 0.3 /* ANGLEPS is the smallest allowed angle for a match. I set it fairly large, but it could also be a function of MAXSTEPS. */ #define ANGLEPS 1e-4 /* params[0] = a in density distribution, params[1] = b */ double density( double p, void *params) { double p2, *b; b = (double*) params; p2 = p*p; return (p2/(exp(-b[0]*p2) + b[1]*p2)); } /* Find all symmetry points in every sextant and octant to a given point. The sextant split line is along theta=phi = PI/4, and the octants are split along the (x,y), (x,z) and (y,z) planes. All points are at same radius so only two angles are needed for each point. Input is phi, theta in 0,0 (octant, sextant) Output is vector list with return value of number of points in the list. First entry includes input value followed by N-1 points in phi, theta order. */ int symmetry( double *source, double *vector) { double phi, theta, pi2, phicheck, thcheck, vec[96]; int i, j, k, n, check; pi2 = M_PI/2.0; phi = source[0]; theta = source[1]; for(i=0; i<6; i++) // outer loop over sextants { k = 16*i; for(j=0; j<16; j+=2) // inner loop over octants { vec[k+j] = phi; vec[k+j+1] = theta; if(j==6 || j==14) theta = - theta; phi += pi2; if(phi > 2.0*M_PI) phi -= 2.0*M_PI; phicheck = phi - 2.0*M_PI; if(phicheck < 0.0) phicheck = -phicheck; if(phicheck < ANGLEPS) phi = 0.0; } switch (i) { case 0: phi = pi2 - source[0]; break; case 1: phi = pi2 - source[1]; theta = source[0]; break; case 2: phi = pi2 - source[1]; theta = pi2 - source[0]; break; case 3: phi = source[1]; theta = pi2 - source[0]; break; case 4: phi = source[1]; theta = source[0]; break; default: break; } } /* see if any of these points are on special case boundary */ check = 0; phi = source[0]; theta = source[1]; phicheck = phi - pi2/2.0; if (phicheck < 0.0) phicheck = -phicheck; if (phicheck < ANGLEPS) check = 1; if (theta < ANGLEPS) check = 1; thcheck = phi - theta; if (thcheck < 0.0) thcheck = -thcheck; if( thcheck < ANGLEPS) check = 1; if (check) { /* brute force simple check to remove duplicate points */ vector[0] = phi; vector[1] = theta; n = 1; for( i=2; i<96; i+=2) { check = 1; j = 0; while( (check > 0) && (j < n)) { k = 2*j; phicheck = vec[i] - vector[k]; if (phicheck < 0.0) phicheck = -phicheck; thcheck = vec[i+1] - vector[k+1]; if (thcheck < 0.0) thcheck = -thcheck; if ((phicheck < ANGLEPS) && (thcheck < ANGLEPS)) check = 0; j++; } if(check) { vector[2*n] = vec[i]; vector[2*n+1] = vec[i+1]; n++; } } return n; } else { for(i=0; i<96; i+=2) { vector[i] = vec[i]; vector[i+1] = vec[i+1]; } return 48; } } /* General integral over 1/48 volume segment. Input is pointer to table of NUMELMNTS which is summed over the volume segment. Note that r^2 dr cos(theta) dtheta dphi is computed here, not in input table. Returns integral value by accounting for duplicate faces. */ double integrate_sphere( double *func) { int i, j, k, idex, kdex; double dr, theta, acp, sum, dsum; sum = func[0]; dr = pow(MAXWALL/MAXSTEPS, 3.0); dr *= M_PI*M_PI/16.0; /* sum over lines which have very limited numbers */ /* x, y, z axies first */ dsum = 0.0; for(i=1; i