a

calcEBV_FFT.cpp

@@ -426,7 +426,7 @@ int calcEBV(fields *fi, par *par)
     float Tcyclotron = 2.0 * pi * mp[0] / (e_charge_mass * (par->Bmax + 1e-5f));
     float acc_e = par->Emax * e_charge_mass;
     float vel_e = sqrt(kb * Temp_e / e_mass);
-    float TE = (sqrt(vel_e * vel_e / (acc_e * acc_e) + 2 * a0 / acc_e) - vel_e / acc_e)*500;
+    float TE = (sqrt(vel_e * vel_e / (acc_e * acc_e) + 2 * a0 / acc_e) - vel_e / acc_e)*1000;
     float TE1= a0/par->Emax *(par->Bmax+.00001);
     TE=max(TE,TE1);
     // cout << "Tcyclotron=" << Tcyclotron << ",Bmax= " << par->Bmax << ", TE=" << TE << ",Emax= " << par->Emax << endl;

generaterandp.cpp

@@ -4,7 +4,7 @@ void generate_rand_sphere(particles *pt, par *par)
     // spherical plasma set plasma parameters
     float Temp[2] = {Temp_e, Temp_d}; // in K convert to eV divide by 1.160451812e4
     // initial bulk electron, ion velocity
-    float v0[2][3] = {{0, 0, -vz0}, {0, 0, vz0/60}};
+    float v0[2][3] = {{0, 0, -vz0}, {0, 0, vz0 / 60}};
 
     float r0 = r0_f * a0; // if sphere this is the radius
     float area = 4 * pi * r0 * r0;
@@ -61,6 +61,7 @@ void generate_rand_sphere(particles *pt, par *par)
                         pt->pos0x[p][na] = ((float)(i - n_space_divx / 2) + (float)rand() / RAND_MAX) * a0;
                         pt->pos0y[p][na] = ((float)(j - n_space_divy / 2) + (float)rand() / RAND_MAX) * a0;
                         pt->pos0z[p][na] = ((float)(k - n_space_divz / 2) + (float)rand() / RAND_MAX) * a0;
+
                         pt->pos1x[p][na] = pt->pos0x[p][na];
                         pt->pos1y[p][na] = pt->pos0y[p][na];
                         pt->pos1z[p][na] = pt->pos0z[p][na];
@@ -75,7 +76,8 @@ void generate_rand_sphere(particles *pt, par *par)
 #pragma omp parallel for ordered
         for (int n = na; n < n_partd; n++)
         {
-            float r = r0 * pow(gsl_ran_flat(rng, 0, 1), 0.3333333333);
+            // float r = r0 * pow(gsl_ran_flat(rng, 0, 1), 0.3333333333);
+            float r = gsl_ran_gaussian(rng, r0) ;
             // float r = r0 * pow(gsl_ran_flat(rng, 0, 1), 0.5);
             double x, y, z;
             gsl_ran_dir_3d(rng, &x, &y, &z);

include/traj_physics.h

@@ -10,24 +10,24 @@ constexpr int f2 = f1 * 1.2;
 constexpr float incf = 1.2f;        // increment
 constexpr float decf = 1.0f / incf; // decrement factor
 
-constexpr int n_space = 200;                                      // should be 2 to power of n for sater FFT
+constexpr int n_space = 128;                                      // should be 2 to power of n for sater FFT
 constexpr float nback = 1;                                       // background particles per cell - improves stability
-constexpr int n_partd = n_space * n_space * n_space * nback * 1; // must be 2 to power of n
+constexpr int n_partd = n_space * n_space * n_space * nback * 8; // must be 2 to power of n
 constexpr int n_parte = n_partd;
 
 constexpr float R_s = n_space / 1;  // LPF smoothing radius
 constexpr float r0_f = 4; //  radius of sphere or cylinder
 
 // The maximum expected E and B fields. If fields go beyond this, the the time step, cell size etc will be wrong. Should adjust and recalculate.
 //  maximum expected magnetic field
-constexpr float Bmax0 = 0.01;    // in T earth's magnetic field is of the order of 1e-4 T 
+constexpr float Bmax0 = 0.01;    // in T earth's magnetic field is of the order of ~ 1e-4 T DPF ~ 100T
 constexpr float Emax0 = 1e5; // 1e11V/m is approximately interatomic E field -extremely large fields implies poor numerical stability
 
 constexpr float Bz0 = 0.01; // in T, static constant fields
 constexpr float Ez0 = 0.0f;//in V/m
 constexpr float vz0 = 0.0f;
-constexpr float a0 = 1.0e-3;      // typical dimensions of a cell in m This needs to be smaller than debye length otherwise energy is not conserved if a particle moves across a cell
-constexpr float target_part = 1e11; // 3.5e22 particles per m^3 per torr of ideal gas. 7e22 electrons for 1 torr of deuterium
+constexpr float a0 = 1.0e-4;      // typical dimensions of a cell in m This needs to be smaller than debye length otherwise energy is not conserved if a particle moves across a cell
+constexpr float target_part = 1e10; // 3.5e22 particles per m^3 per torr of ideal gas. 7e22 electrons for 1 torr of deuterium
 
 // technical parameters
 
@@ -40,7 +40,7 @@ constexpr int n_output_part = (n_partd > 9369) ? 9369 : n_partd; // maximum numb
 // const int nprtd=floor(n_partd/n_output_part);
 
 constexpr int ndatapoints = 300; // total number of time steps to calculate
-constexpr int nc1 = 2;            // f1 * 1;      // number of times to calculate E and B between printouts
+constexpr int nc1 = 1;            // f1 * 1;      // number of times to calculate E and B between printouts
 constexpr int md_me = 60;        // ratio of electron speed/deuteron speed at the same KE. Used to calculate electron motion more often than deuteron motion
 
 #define Hist_n 512