/******************************************************************** * * * 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 10 /* 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]; gsl_integration_workspace *w_theta; gsl_integration_workspace *w_sphere; /* 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)); } /* Set up parameters for integration param[0] = coefficient of exponential in density distribution param[1] = coefficient of r^2 in density distribution param[2] = average temperature coefficient for potential param[3] = radius r param[4] = angle from x, y plane (theta) param[5] = angle from x axis (phi) param[6] = radius r' param[7] = angle from x, y plane (theta') param[8] = angle from x axis (phi') Note that last 3 parameters are being integrated over so I need 3 different subroutines. All end up at same place. */ /* Subroutine to find zone 0 equivelent position. Mapping is if theta < 0, theta -> -theta if phi > PI/2, map phi to first 0 - PI/2 it fist well with. then find correct sextant. Order of inputs and outpupts is r[0] = r, r[1] = theta, r[2] = phi. */ void map( double *r, double *rp) { double theta, phi, test1, test2; if(r[1] < 0.0) theta = -r[1]; else theta = r[1]; test1 = M_PI/2.0; if( r[2] >= 3.0*test1) phi = 2.0*M_PI - r[2]; else if( r[2] >= 2.0*test1) phi = r[2] - M_PI; else if( r[2] >= test1) phi = M_PI - r[2]; else phi = r[2]; rp[0] = r[0]; test2 = M_PI/4.0; if(theta > test2) { if( phi > test2) { rp[1] = test1 - phi; rp[2] = test1 - theta; return; } else { rp[1] = phi; rp[2] = test1 - theta; return; } } if( phi > test2) { if( theta > test1 - phi) { rp[1] = test1 - phi; rp[2] = theta; return; } else { rp[1] = theta; rp[2] = test1 - phi; return; } } if (theta > phi) { rp[1] = phi; rp[2] = theta; return; } rp[1] = theta; rp[2] = phi; } /* given a position sphere, return parameters which allow interpolation in either potential or vector table. input is r, theta, phi in pointer *r output is index as the set 1,1,1 ; 0,0,0 1,1,0 ; 0,0,1 1,0,1 ; 0,1,0 1,0,0 ; 0,1,1 which is 0 or 1 plus i, j, or k index from the closest corner to the input point. the pointer *delta contains the difference computed for a line going across each corner: dr + dtheta + dphi dr + dtheta - dphi dr - dtheta + dphi dr - dtheta - dphi If the input point is close to a plane so that 1 dimension is lost, or close to a line so that 2 dimensions are lost, we get a special case. The return will be different indexes, but the operations on the data returned will be the same. Some of the *delta values will be zero so they will not be included. In the rare instance when an exact value is the closest match, all the delta's will be zero so that the [i, j, k] index will be the correct one to use. This will always be returned in the same offset. */ void interp( double *r, int *index, double *delta) { double rp[3], dr, dp1, dt1, da, dp2, dt2, dt, dp; int i, j, k, i0, i1, j0, j1, l; dr = MAXWALL/MAXSTEPS; map(r, rp); // convert to 0th sextant dr = rp[0]/dr; i = floor(dr); dt = 4.0/M_PI*rp[1]*i; j = floor(dt); dp = 4.0/M_PI*rp[2]*i; k = floor(dp); dr -= i; dt -= j; dp -= k; i0 = i*(i+1)*(i+2)/6; i1 = (i+1)*(i+2)*(i+3)/6; j0 = j*(j+1)/2; j1 = (j+1)*(j+2)/2; for(l=0; l<8; l++) index[l] = 0; index[1] = i0 + j0 + k; if( dr < 1e-5) { // dr = 0, check if dtheta or dphi are also if( dt < 1e-5) { if( dp < 1e-5) // dtheta and dphi also small, return point { for(l=0; l<4; l++) delta[l] = 0.0; return; } index[0] = i0 + j0 + k + 1; delta[0] = dp; return; } if( dp < 1e-5) { index[0] = i0 + j1 + k; delta[0] = dt; return; } index[0] = i0 + j1 + k + 1; index[2] = i0 + j0 + k + 1; index[3] = i0 + j1 + k; delta[0] = dt + dp; delta[1] = dt + dp; return; } da = M_PI/4.0/i/(i+1); dt1 = (i-j)*da; dp1 = (i-k)*da; dt2 = (i+j+1)*da; dp2 = (i+k+1)*da; dt = rp[1] - j*M_PI/4.0/i; dp = rp[2] - k*M_PI/4.0/i; if( dt < 1e-5) { if( dp < 1e-5) { index[0] = i1 + j0 + k; delta[0] = dr; return; } index[0] = i1 + j0 + k + 1; index[2] = i1 + j0 + k; index[3] = i0 + j0 + k + 1; delta[0] = dr + dp/dp1; delta[1] = dr - dp/dp2; return; } if( dp < 1e-5) { index[0] = i1 + j1 + k; index[2] = i1 + j0 + k; index[3] = i0 + j1 + k; delta[0] = dr + dt/dt1; delta[1] = dr - dt/dt2; return; } delta[0] = dr + dt/dt1 + dp/dp1; delta[1] = dr + dt/dt1 - dp/dp2; delta[2] = dr - dt/dt2 + dp/dp1; delta[3] = dr - dt/dt2 - dp/dp2; index[0] = i1 + j1 + k+1; index[1] = i0 + j0 + k; index[2] = i1 + j1 + k; index[3] = i0 + j0 + k+1; index[4] = i1 + j0 + k+1; index[5] = i0 + j1 + k; index[6] = i1 + j0 + k; index[7] = i0 + j1 + k+1; } /* main integral subroutine. Enter with parameters as above. Assumes GridPot = grid potental loaded externally is global and that ElectPot = electron potential computed is also global. */ double electron_core( double *param) { double temp1, temp2, psi, fi, delta[4], r2, c1, c2, c3; double test1, test2, test3; int index[8]; /* find index into potential tables and compute exponential term */ interp( ¶m[6], index, delta); psi = (GridPot[index[0]] - GridPot[index[1]])*delta[0]; psi += (GridPot[index[2]] - GridPot[index[3]])*delta[1]; psi += (GridPot[index[4]] - GridPot[index[5]])*delta[2]; psi += (GridPot[index[6]] - GridPot[index[7]])*delta[3]; psi += GridPot[index[1]]; fi = (ElecPot0[index[0]] - ElecPot0[index[1]])*delta[0]; fi += (ElecPot0[index[2]] - ElecPot0[index[3]])*delta[1]; fi += (ElecPot0[index[4]] - ElecPot0[index[5]])*delta[2]; fi += (ElecPot0[index[6]] - ElecPot0[index[7]])*delta[3]; fi += ElecPot0[index[1]]; temp1 = exp(-(psi+fi)/param[2]); /* next compute radial density function from parameters and input radius */ r2 = param[6]*param[6]; /* temp2 = exp(-param[0]*r2) + param[1]*r2; temp1 /= temp2; /* finally do geometric terms */ c1 = cos(param[4] - param[7]); c2 = cos(param[4] + param[7]); c3 = cos(param[5] - param[8]); temp2 = param[3]* param[6]*(c1*(1.0 + c3) - c2*(1.0 - c3)); temp2 = sqrt(param[3]*param[3] + r2 - temp2); temp1 /= temp2; return temp1; } /* subroutines for potential distribution integration Because this is a triple integral, each subroutine calls an integral over another dimension. */ double distribution_phi( double phi, void *params) { double *p; p = (double*)params; p[8] = phi; return electron_core(p); } double theta_integral( double theta, void *params) { double ring, err, *p; gsl_function F; int i; p = (double*)params; p[7] = theta; F.function = &distribution_phi; F.params = params; gsl_integration_qag(&F, 0.0, M_PI*2.0, 1e-4, 1e-4, 500, GSL_INTEG_GAUSS15, w_theta, &ring, &err); ring *= cos(theta); return ring; } double sphere( double r, void *params) { double sphr, err, *p; double r2, temp2; int i; gsl_function S; if(r < MAXWALL/MAXSTEPS) return 0.0; p = (double*)params; p[6] = r; S.function = &theta_integral; S.params = params; gsl_integration_qag(&S, -M_PI/2.0, M_PI/2.0, 1e-4, 1e-4, 500, GSL_INTEG_GAUSS15, w_sphere, &sphr, &err); /* next compute radial density function from parameters and input radius */ r2 = p[6]*p[6]; temp2 = exp(-p[0]*r2) + p[1]*r2; sphr *= r2/temp2; return sphr; } main() { FILE *readmag, *fi_pot; int numread; char filename[64], input[8]; double rtlde, rmax, abserr, densparms[9]; double across, up, fi, r, phi, theta; int err, numval, i, j, k, idex, jdex, dex; /* verify that map subroutines does what I want for(i=0; i<10; i++) { theta = M_PI/2.0 - (i+1)*M_PI/11.0; for(j=0; j<20; j++) { phi = 2.0*M_PI - (j+1)*M_PI/11.0; densparms[1] = theta; densparms[2] = phi; map( densparms, &densparms[6]); printf("%lf %lf -> %lf %lf\n", theta, phi, densparms[7], densparms[8]); } } exit(0); */ gsl_integration_workspace *w; gsl_function F; /* set up global work spaces so allocation is only done once */ w_theta = gsl_integration_workspace_alloc(500); w_sphere = gsl_integration_workspace_alloc(500); /* compute normalization coefficient for density distribution chosen */ densparms[0] = 9.0; densparms[1] = 10.0; F.function = &density; F.params = densparms; rmax = MAXWALL; err = gsl_integration_qng( &F, 0.0, rmax, 1e-6, 1e-6, &rtlde, &abserr, (size_t*)&numval); printf("Density integral result = %lf\n", rtlde); printf("absolute error = %lf\n", abserr); /* add electron temperature parameter to list of parameters */ densparms[2] = 0.01; // start small! /* load grid potential from previous calculation */ sprintf(filename, "potential.dat"); readmag = fopen(filename, "rb"); if( !readmag) { printf("can't read in mag data file %s.\n", filename); exit(0); } printf("opening data files and reading them in.\n"); numread = fread(GridPot, sizeof(double), NUMELMNTS, readmag); if(numread < NUMELMNTS) printf("\nWARNING: not all field file read in. \n\n"); fclose(readmag); /* create file to hold computed output from this program */ fi_pot = fopen("electron_potential.dat", "w"); if (!fi_pot) { printf("can't create electric potential data file. \n"); exit(0); } /* integrate over the entire volume */ w = gsl_integration_workspace_alloc(500); across = M_PI/4.0; // x, y plane angle up = across; // z up from xy plane angle F.function = &sphere; F.params = densparms; /* I need to ping pong solution, so start with nothing and see how fast it stabilizes */ for( i=0; i