Mirror 3 d

2020-10-mirror_3D_scripts/dpitch_long_i.txt

@@ -0,0 +1,15 @@
+ 7.952452671E-02 1.129492901E-01 5.941727353E-02 4.820582348E-02 4.205907478E-02 3.804727771E-02 3.518877981E-02 3.304121677E-02 3.137077637E-02 3.004056109E-02
+ 2.896447062E-02 2.808543543E-02 2.736411117E-02 2.677256624E-02 2.629055317E-02 2.590320725E-02 2.559957714E-02 2.537166455E-02 2.521378977E-02 2.512217608E-02
+ 2.509468969E-02 2.513069855E-02 2.523103142E-02 2.539803169E-02 2.563571299E-02 2.595003751E-02 2.634935658E-02 2.684508125E-02 2.745269801E-02 2.819332669E-02
+ 2.909617023E-02 3.020250455E-02 3.157247937E-02 3.329740395E-02 3.552366840E-02 1.852512870E-02 5.273452973E-04 1.051280553E-03 1.574664544E-03 2.093190533E-03
+ 2.606449645E-03 3.114289632E-03 3.612702340E-03 4.103848307E-03 4.584888596E-03 5.052951443E-03 5.509961183E-03 5.952774356E-03 6.380680873E-03 6.793063026E-03
+ 7.187307084E-03 7.563505927E-03 7.919970508E-03 8.254927667E-03 8.567450786E-03 8.855783065E-03 9.117702407E-03 9.350836405E-03 9.553584580E-03 9.720618113E-03
+ 9.849303078E-03 9.936221244E-03 9.972213075E-03 9.952859362E-03 9.867414270E-03 9.702892042E-03 9.449252096E-03 9.080992811E-03 8.580439098E-03 7.927108237E-03
+ 7.086565160E-03 6.044876865E-03 4.789057312E-03 3.328063965E-03 1.021491200E-01 1.021491200E-01 3.328063965E-03 4.789057312E-03 6.044876865E-03 7.086565160E-03
+ 7.927108237E-03 8.580439098E-03 9.080992811E-03 9.449252096E-03 9.702892042E-03 9.867414270E-03 9.952859362E-03 9.972213075E-03 9.936221244E-03 9.849303078E-03
+ 9.720618113E-03 9.553584580E-03 9.350836405E-03 9.117702407E-03 8.855783065E-03 8.567450786E-03 8.254927667E-03 7.919970508E-03 7.563505927E-03 7.187307084E-03
+ 6.793063026E-03 6.380680873E-03 5.952774356E-03 5.509961183E-03 5.052951443E-03 4.584888596E-03 4.103848307E-03 3.612702340E-03 3.114289632E-03 2.606449645E-03
+ 2.093190533E-03 1.574664544E-03 1.051280553E-03 5.273452973E-04 1.852512870E-02 3.552366840E-02 3.329740395E-02 3.157247937E-02 3.020250455E-02 2.909617023E-02
+ 2.819332669E-02 2.745269801E-02 2.684508125E-02 2.634935658E-02 2.595003751E-02 2.563571299E-02 2.539803169E-02 2.523103142E-02 2.513069855E-02 2.509468969E-02
+ 2.512217608E-02 2.521378977E-02 2.537166455E-02 2.559957714E-02 2.590320725E-02 2.629055317E-02 2.677256624E-02 2.736411117E-02 2.808543543E-02 2.896447062E-02
+ 3.004056109E-02 3.137077637E-02 3.304121677E-02 3.518877981E-02 3.804727771E-02 4.205907478E-02 4.820582348E-02 5.941727353E-02 1.129492901E-01 7.952452671E-02

2020-10-mirror_3D_scripts/interpolate_numerical_distribution.m

@@ -0,0 +1,61 @@
+function F_out = interpolate_numerical_distribution()
+
+    function out = F2(u,z)
+        out = ones(numel(u),numel(z));
+        if size(u) == size(z)
+            for i = 1:numel(u)
+                for j = 1:numel(z)
+                    out(i,j) = F(u(i),z(j));
+                end
+            end
+        elseif isscalar(u) 
+            for j = 1:numel(z)
+                out(j) = F(u,z(j));
+            end
+        end
+    end
+
+pitch_file = fopen('pitch_i.txt', 'r');
+vel_file = fopen('u_j.txt', 'r');
+dist_file = fopen('f_ij.txt', 'r');
+
+z = fscanf(pitch_file, '%f');
+u = fscanf(vel_file, '%f');
+f_vals = fscanf(dist_file, '%f');
+
+size_f = length(f_vals);
+num = sqrt(size_f);
+
+n_o = 8e14;
+temp = 1.6022e-10; %temperature in erg
+m_e = 9.109*10^-28;
+u_th = sqrt(2*temp/m_e);
+lIndex = [0 1 2];
+
+%[u,z] = meshgrid(vel_vals,pitch_vals);
+%maxwell = n_o.*exp(-u.^2./u_th^2).*(z.*0 + 1)./(u_th.^3.*pi^(3/2));
+
+analytic = @(x,y) n_o.*exp(-x.^2./u_th^2).*(y.*0 + 1)./(u_th.^3.*pi^(3/2));
+f_vals = reshape(f_vals, [num num]);
+
+% ax1 = subplot(1,3,1);
+% contour(vel_vals,pitch_vals,f_vals);
+% ax2 = subplot(1,3,2);
+% contour(u,z,analytic(u,z));
+
+%u = u';
+
+%z = z';
+
+F = griddedInterpolant({u,z},f_vals');
+F_out = @F2;
+
+%ax3 = subplot(1,3,3);
+%f_plot = F_out(u(:),z(:));
+%f_plot = reshape(f_plot, size(u));
+%contour(u,z,f_plot);
+% 
+
+
+
+end

2020-10-mirror_3D_scripts/pitch_i.txt

@@ -0,0 +1,15 @@
+ 0.00000E+00 1.59049E-01 2.25899E-01 2.77884E-01 3.22310E-01 3.62002E-01 3.98405E-01 4.32379E-01 4.64487E-01 4.95121E-01
+ 5.24568E-01 5.53050E-01 5.80739E-01 6.07778E-01 6.34284E-01 6.60359E-01 6.86091E-01 7.11558E-01 7.36834E-01 7.61986E-01
+ 7.87078E-01 8.12175E-01 8.37340E-01 8.62637E-01 8.88136E-01 9.13909E-01 9.40036E-01 9.66607E-01 9.93726E-01 1.02151E+00
+ 1.05011E+00 1.07971E+00 1.11052E+00 1.14285E+00 1.17711E+00 1.21390E+00 1.21416E+00 1.21495E+00 1.21627E+00 1.21810E+00
+ 1.22045E+00 1.22331E+00 1.22668E+00 1.23054E+00 1.23489E+00 1.23971E+00 1.24499E+00 1.25073E+00 1.25690E+00 1.26349E+00
+ 1.27049E+00 1.27787E+00 1.28561E+00 1.29371E+00 1.30212E+00 1.31084E+00 1.31983E+00 1.32908E+00 1.33854E+00 1.34818E+00
+ 1.35798E+00 1.36788E+00 1.37785E+00 1.38783E+00 1.39776E+00 1.40756E+00 1.41716E+00 1.42646E+00 1.43532E+00 1.44362E+00
+ 1.45118E+00 1.45779E+00 1.46327E+00 1.46737E+00 1.46992E+00 1.67167E+00 1.67422E+00 1.67833E+00 1.68380E+00 1.69042E+00
+ 1.69797E+00 1.70627E+00 1.71513E+00 1.72443E+00 1.73403E+00 1.74384E+00 1.75377E+00 1.76374E+00 1.77371E+00 1.78362E+00
+ 1.79341E+00 1.80306E+00 1.81252E+00 1.82176E+00 1.83075E+00 1.83947E+00 1.84789E+00 1.85598E+00 1.86373E+00 1.87111E+00
+ 1.87810E+00 1.88469E+00 1.89086E+00 1.89660E+00 1.90188E+00 1.90670E+00 1.91105E+00 1.91491E+00 1.91828E+00 1.92114E+00
+ 1.92349E+00 1.92533E+00 1.92664E+00 1.92743E+00 1.92770E+00 1.96448E+00 1.99874E+00 2.03107E+00 2.06189E+00 2.09148E+00
+ 2.12008E+00 2.14787E+00 2.17499E+00 2.20156E+00 2.22768E+00 2.25346E+00 2.27896E+00 2.30425E+00 2.32942E+00 2.35451E+00
+ 2.37961E+00 2.40476E+00 2.43003E+00 2.45550E+00 2.48123E+00 2.50731E+00 2.53381E+00 2.56085E+00 2.58854E+00 2.61702E+00
+ 2.64647E+00 2.67711E+00 2.70921E+00 2.74319E+00 2.77959E+00 2.81928E+00 2.86371E+00 2.91569E+00 2.98254E+00 3.14159E+00

2020-10-mirror_3D_scripts/pitch_long_i.txt

@@ -0,0 +1,15 @@
+ 0.000000000E+00 1.590490534E-01 2.258985802E-01 2.778836005E-01 3.223102271E-01 3.620017501E-01 3.984047825E-01 4.323793097E-01 4.644872161E-01 4.951208624E-01
+ 5.245683382E-01 5.530498037E-01 5.807392091E-01 6.077780260E-01 6.342843416E-01 6.603591323E-01 6.860907561E-01 7.115582866E-01 7.368340852E-01 7.619858662E-01
+ 7.870784373E-01 8.121752455E-01 8.373398344E-01 8.626373084E-01 8.881358978E-01 9.139087343E-01 9.400359728E-01 9.666074475E-01 9.937261353E-01 1.021512844E+00
+ 1.050112789E+00 1.079705184E+00 1.110517798E+00 1.142850143E+00 1.177112606E+00 1.213897480E+00 1.214162863E+00 1.214952170E+00 1.216265424E+00 1.218101499E+00
+ 1.220451805E+00 1.223314399E+00 1.226680385E+00 1.230539803E+00 1.234888081E+00 1.239709580E+00 1.244993984E+00 1.250729503E+00 1.256899533E+00 1.263490865E+00
+ 1.270485659E+00 1.277865479E+00 1.285612671E+00 1.293705420E+00 1.302122526E+00 1.310840321E+00 1.319834092E+00 1.329075726E+00 1.338535765E+00 1.348182895E+00
+ 1.357977001E+00 1.367881501E+00 1.377849444E+00 1.387825928E+00 1.397755162E+00 1.407560756E+00 1.417160946E+00 1.426459260E+00 1.435322932E+00 1.443620138E+00
+ 1.451177149E+00 1.457793269E+00 1.463266902E+00 1.467371383E+00 1.469923030E+00 1.671669623E+00 1.674221270E+00 1.678325751E+00 1.683799385E+00 1.690415505E+00
+ 1.697972515E+00 1.706269722E+00 1.715133393E+00 1.724431707E+00 1.734031898E+00 1.743837491E+00 1.753766726E+00 1.763743210E+00 1.773711152E+00 1.783615652E+00
+ 1.793409758E+00 1.803056889E+00 1.812516928E+00 1.821758561E+00 1.830752332E+00 1.839470128E+00 1.847887234E+00 1.855979983E+00 1.863727175E+00 1.871106995E+00
+ 1.878101789E+00 1.884693121E+00 1.890863151E+00 1.896598670E+00 1.901883073E+00 1.906704572E+00 1.911052850E+00 1.914912269E+00 1.918278255E+00 1.921140848E+00
+ 1.923491154E+00 1.925327229E+00 1.926640483E+00 1.927429790E+00 1.927695174E+00 1.964480048E+00 1.998742511E+00 2.031074856E+00 2.061887470E+00 2.091479865E+00
+ 2.120079810E+00 2.147866518E+00 2.174985206E+00 2.201556681E+00 2.227683919E+00 2.253456756E+00 2.278955345E+00 2.304252819E+00 2.329417408E+00 2.354514216E+00
+ 2.379606787E+00 2.404758568E+00 2.430034367E+00 2.455501898E+00 2.481233521E+00 2.507308312E+00 2.533814628E+00 2.560853444E+00 2.588542850E+00 2.617024315E+00
+ 2.646471791E+00 2.677105438E+00 2.709213344E+00 2.743187871E+00 2.779590904E+00 2.819282426E+00 2.863709053E+00 2.915694073E+00 2.982543600E+00 3.141592654E+00

2020-10-mirror_3D_scripts/plot_mirror_output.m

@@ -1,45 +1,196 @@
-function plot_mirror_output(nodes, f_nD, f_nD_analytic, pde)
+function plot_mirror_output(nodes, outputs, pde, opts)
 
-    x = nodes{1};
-    y = nodes{2};
-    z = nodes{3};
-
-    nx = numel(x);
-    ny = numel(y);
-    nz = numel(z);
-
-    %%
+      n_cgs = 8e14; %equilibrium density in cm.^-3
+      m_e_cgs = 9.109*10^-28; %electron mass in g
+      m_D_cgs = 3.3443*10^-24; %Deuterium mass in g
+      m_He_cgs = 6.7*10^-24; %helium 4 mass in g
+      m_B_cgs = 1.82*10^-23; %Boron 11 mass in g
+      m_Ne_cgs = 3.3509177*10^-23; %Neon mass in g
+      temp_cgs = 1.6022e-10; %temperature in erg
+      params_cgs.a.n = n_cgs;
+      params_cgs.b.n = n_cgs;
+      params_cgs.b2.n = n_cgs;
+      %         params_cgs.a.vth = sqrt(2*temp_cgs/m_e_cgs);
+      %         params_cgs.b.vth = sqrt(2*temp_cgs/m_e_cgs);
+      params_cgs.a.m = m_e_cgs; %beam is electrons
+      params_cgs.b.m = m_D_cgs; %background ions
+      params_cgs.b2.m = m_e_cgs; %background electrons
+      params_cgs.e = 4.803*10^-10; %charge in Fr
+      params_cgs.c = 3*10^10; %speed of light in cm/s
+      params_cgs.a.theta = pde.params.a.T_eV/(params_cgs.a.m*params_cgs.c^2);
+      nu_ee = 4*pi*params_cgs.a.n.*params_cgs.e.^4*pde.params.ln_delt/(params_cgs.a.theta^(3/2)...
+          *params_cgs.a.m^3*params_cgs.c^3);
+      franz_conduct = params_cgs.a.n*params_cgs.e^2/(pde.params.b.Z*params_cgs.a.m*nu_ee);
+      x = nodes{1};
+      y = nodes{2};
+      nx = numel(x);
+      ny = numel(y);
+      sx = max(1, floor(nx/2));
+      sy = max(1, floor(2*ny/8));
+      energy_func_v = @(x) 0.5*pde.params.a.m.*x.^2/(1.602*10^-19);
+      energy_func_z = @(z) z.*0 + 2*pi;
+      current_func_v = @(x) 2*pi*pde.params.a.Z.*pde.params.e.*x;
+      current_func_z = @(x) cos(x);
+      coord = get_realspace_coords(pde,nodes);
+      mass_func = @(x) x.*0 + sqrt(2*pi);
+      zero_func = @(x) x.*0;
+      num_dims = numel(pde.dimensions);
+      moment_func_nD{1} = energy_func_v;
+      moment_func_nD{2} = energy_func_z;
+      mass_func_nD{1} = mass_func;
+      mass_func_nD{2} = mass_func;
+      current_func_nD{1} = current_func_v;
+      current_func_nD{2} = current_func_z;
+      x_E = 0.5.*pde.params.a.m*x.^2/pde.params.e;
+      num_steps = length(outputs.time_array);
+      spitzer_conduct = (3/(4*sqrt(2*pi)))*(4*pi*pde.params.eps0)^2*(1.38e-23*pde.params.a.T_eV*11604)^(3/2)...
+          /(pde.params.b.Z*pde.params.e^2*pde.params.a.m^(1/2)*pde.params.ln_delt);
+      spitzer_conduct2 = 4*pi*pde.params.eps0^2*(pde.params.a.m*pde.params.a.vth^2)^(3/2) ... 
+          /(pde.params.a.m^(1/2)*pde.params.e^2*pde.params.ln_delt*pde.params.b2.Z);
+      f1d_analytic = outputs.f_realspace_analytic_nD_t{1,num_steps};
+      f1d_ic = outputs.f_realspace_analytic_nD_t{1,1};
+      %getting energy density
+      lev_vec = [opts.lev, opts.lev];
+     % data = struct2cell(outputs);
+     % f_data = cell2mat(data{6,1});
+     % f_data = reshape(f_data, [num_steps size(f1d_ic)]);
+%      [x_grid, y_grid] = meshgrid(outputs.time_array,x_E);
+      outputs.time_array(1) = 0;
+for i = 1:numel(outputs.time_array)
+      f_data(i,:) = outputs.f_realspace_nD_t{i}(ny/2,:);
+end
+      %contourf(outputs.time_array,x_E,real(f_data)',50);
+%      f_data = permute(f_data,[1 3 2]);
+
+%      contour(outputs.time_array',x',squeeze(f_data(:,:,sy))');
+for j = 1:num_steps
+    %plot(x,f1d);
+    hold on
+    fval_realspace = reshape(outputs.f_realspace_nD_t{j}, [numel(f1d_ic) 1]);
+    mass_vals(j) = moment_integral(lev_vec, opts.deg, coord, fval_realspace, mass_func_nD, pde.dimensions);
+    mass = moment_integral(lev_vec, opts.deg, coord, fval_realspace, mass_func_nD, pde.dimensions);
+    conduct_vals(j) = moment_integral(lev_vec,opts.deg,coord,fval_realspace, current_func_nD,pde.dimensions)/(pde.params.E);
+    conduct_vals(j) = conduct_vals(j)/spitzer_conduct;
+    energy_vals(j) = moment_integral(lev_vec, opts.deg, coord, fval_realspace, moment_func_nD, pde.dimensions)./mass;
+  %        vel_vals(i) = sum(outputs.f_realspace_nD_t{1,i}(:).*coord{:,:}*1e-4);
+%    x_hint_t = @(t) -x_hint(t,pde.params.a.vel_vals(j),pde.params.a,pde.params.b);
+     %e_hint = integral(x_hint_t, 0, outputs.time_array(j)); 
+    %hint_func(j) = energy_vals(1)*exp(-e_hint);
+     for i = 1:length(x)
+          f1d_new(i) = fval_realspace(i);
+          f1d_analytic_new(i) = f1d_analytic(i);
+     end
+     time = outputs.dt*j;
+     %loglog(x,f_data(j,:),'DisplayName',[num2str(time)]);
+     hold on
+end
+legend();
+hold off
+      %formulating the Hinton solution
+     % timespan = [0 1e-3]; %time span in seconds
+     % x0 = 3e3; %initial energy in eV
+     % y = @(x) sqrt(1.602e-19.*x./(0.5.*3.3443e-27)); %velocity function in terms of energy
+     % [t,E] = ode45(@(t,E) (-params.nu_E(y(E),params.a,params.b)-params.nu_E(y(E),params.a,params.b2)).*E,timespan,x0);
+     %plot(x,f1d_ic,'--');
+     %plotting Hinton solution alongside numerical values
+     %x = x./x0;
+     %hold on
+     %plot(t,E,'r','LineWidth', 2);
+     %plot(outputs.time_array,energy_vals,'-o','LineWidth',1);
+     %figure
+     %plot(outputs.time_array,mass_vals,'-o','LineWidth',2);
+     %hold off
+     %f_slice = outputs.f_realspace_nD_t{1,num_steps};
+     %v_norm = x./pde.params.b.vth;
+     %plot(v_norm,f_slice);
+     %figure
+     plot(outputs.time_array,conduct_vals,'-','LineWidth',2);
+     hold on
+     %loglog(x_E,f1d_analytic_new, '-o');
+     %semilogy(outputs.time_array,hint_func,'k');
+    % hold off
+%      hold off;
+%     y = nodes{2};
+%     z = nodes{3};
+%      nx = numel(x);
+%      ny = numel(y);
+%      ax1 = subplot(2,2,1);
+%      contourf(y,x,squeeze(f_nD{1,1}(:,:))');
+%      ax2 = subplot(2,2,2);
+%      contourf(y,x,squeeze(f_nD{1,num_steps}(:,:))');
+     %title('Pitch vs Velocity');
+
+%      sy = numel(f_nD{1,1}(:,1))/4;
+%      sx = 1;
+%      f1d = f_nD{1,num_steps}(:,sx);
+%      f1d_ic = f_nD{1,1}(:,sx);
+%      ax3 = subplot(2,2,3);
+%      plot(y,f1d,'-o');
+%      hold on
+%      plot(y,f1d_ic, '-');
+%      hold off
+%      f1d = f_nD{1,num_steps}(sy,:);
+%      f1d_ic = f_nD{1,1}(sy,:);
+%      ax4 = subplot(2,2,4);
+%      plot(x,f1d,'-o');
+%      hold on
+%      plot(x,f1d_ic, '-');
+
+%     nz = numel(z);     
+%     num_dims = 3;
+%     subset_dimensions = 2;
+%     deg = 6;
+%     space_func = @(x) x.*0 + 1;
+%     energy_func = @(x) x.^2;
+%     perp_func = @(x) sin(x).^2;
+%     par_func = @(x) cos(x).^2;
+%     coord = get_realspace_coords(pde,nodes);
+%     moment_func_nD = {energy_func,perp_func,space_func};
+%    % v_perp_temp = moment_integral(pde.get_lev_vec,deg,coord,f_nD,moment_func_nD,pde.dimensions,subset_dimensions);
+%     moment_func_nD = {energy_func,par_func,space_func};
+%    % v_par_temp = moment_integral(pde.get_lev_vec,deg,coord,f_nD,moment_func_nD,pde.dimensions,subset_dimensions);
+%     %%
     % Plot a 1D line through the solution
+    for i = 1:nx
+        for j = 1:ny
+            v_par(i,j) = x(i).*cos(y(j));
+            v_perp(i,j) = x(i).*sin(y(j));
+        end
+    end
+%     
+%     sz = numel(f_nD{1,1}(:,1,1))/2 + 1;
+
 
-    sz = numel(f_nD(:,1,1))/2;
-    sy = numel(f_nD(1,:,1))/2;
-    sx = numel(f_nD(1,1,:)) - 1;
-
+%    contourf(v_par,v_perp,outputs.f_realspace_nD_t{1,num_steps}(:,sx));
+%     hold on
+%plot(v_par,fval_realspace);
+%     plot(z,v_perp_temp, 'r');
+%     plot(z,v_par_temp, 'b');
+%     hold off
     %plotting 
-    f1d = f_nD(:,sy,sx);
-    f1d = f1d(1,:);
-    x = nodes{1};
-    y = nodes{2};
-    z = nodes{3};
-    x = 9.109*10^-31*x.^2/(1.602*10^-19);
-    plot(z,f1d,'-o');
+%     f1d = f1d(1,:);
+%     x = nodes{1};
+%     y = nodes{2};
+%     z = nodes{3};
+%     x = 9.109*10^-31*x.^2/(1.602*10^-19);
     %xlim([0.01,1e4]);
     %ylim([0.01,1e4]);
-    title('Mirror 2D Relaxation');
-
-    %%
-
-     %%
-    % Overplot analytic solution
-
-    if pde.checkAnalytic
-        f_1d_analytic = f_nD_analytic(:,sy,sx);
-        f_1d_analytic = f_1d_analytic(1,:);
-        hold on;
-        plot(z,f_1d_analytic,'-');
-        hold off;
-    end
-
+%     title('Mirror 2D Relaxation');
+%     
+%     %%
+%     
+%      %%
+%     % Overplot analytic solution
+%     
+% %     if pde.checkAnalytic
+% %         f_1d_analytic = f_nD_analytic(:,sy,sx);
+% %         f_1d_analytic = f_1d_analytic(1,:);
+%           hold on;
+%           f_1d_analytic = f_nD_analytic{1,2};
+%           plot(x,f_1d_analytic,'-');
+%           hold off;
+% %     end
+%     
 %     sx = max(1,floor(nx/2));
 %     if nx > 2
 %         sx = sx+2; % just to get off the exact middle
@@ -56,12 +207,22 @@ function plot_mirror_output(nodes, f_nD, f_nD_analytic, pde)
 %     %%
 %     % Plot 2D
 %     
-%     ax1 = subplot(2,2,3);
-%     f_nD_with_noise = f_nD;
-%     f_nD_with_noise(1,1) = f_nD_with_noise(1,1)*1.0001;
-%     contourf(x,y,f_nD_with_noise,'LineColor','none');
-%     title('numeric 2D solution');
+%     ax1 = subplot(2,3,1);
+%     contourf(z,y,squeeze(f_nD{1,num_steps}(:,:,sx))');
+%     title('Space vs. Pitch');
+%     ax2 = subplot(2,3,2);
+%     contourf(z,x,squeeze(f_nD{1,num_steps}(:,sy,:))');
+%     title('Space vs. Velocity');
+%     %sx = 46;
+%     ax3 = subplot(2,3,3);
+% 
 %     
-
-
+%     n_total = zeros(nz, num_steps);
+%     for i = 1:num_steps
+%          n_total(:,i) = squeeze(sum(f_nD{1,i}, [2 3]));
+%     end
+%     ax4 = subplot(2,3,5);
+%     contourf(z,time_array,squeeze(n_total(:,:))');
+%     title('Number Density vs. Time');
+
 end