2 % To run
this script, obey the following rules:
3 % 1. Run Section 1 first. It loads our various quantities like L, epc, alpha and
PSD
4 % 2. THEN change appropriate
parameters in Parameters.ini, and execute the
5 % VERB code. To see what
parameters you
'd want to change, check
7 % 3. Run the second section
8 % 4. Then run section 3,4,5 - it's upto you. Consult the README
9 % to know what you
're looking for.
21 PSD_filename = 'C:\Users\Akshat Mahajan\Desktop\Paper Reproduction\Time dependent radial diffusion modeling of relativistic electrons with realistic loss rates\1D Radial Diffusion\Execute\Output\OutPSD_rad.dat
';
22 Grid_filename = 'C:\Users\Akshat Mahajan\Desktop\Paper Reproduction\Time dependent radial diffusion modeling of relativistic electrons with realistic loss rates\1D Radial Diffusion\Execute\Output\rad_grid.plt
';
23 PSD0_filename = 'C:\Users\Akshat Mahajan\Desktop\Paper Reproduction\Time dependent radial diffusion modeling of relativistic electrons with realistic loss rates\1D Radial Diffusion\Execute\Input\myPSD0.plt
';
26 PSD_filename = '/Users/alecjen/Desktop/VERB2/1D Radial Diffusion (2004)/Execute/
Output/OutPSD_rad.dat
';
27 Grid_filename = '/Users/alecjen/Desktop/VERB2/1D Radial Diffusion (2004)/Execute/
Output/rad_grid.plt
';
28 PSD0_filename = '/Users/alecjen/Desktop/VERB2/1D Radial Diffusion (2004)/Execute/Input/myPSD0.plt
';
31 [L, epc, alpha, pc] = load_plt(Grid_filename, 'squeeze
', 'permute
');
37 % Run this section after Parameters.ini has been modified and VERB code has been run
39 % Make sure tau and kp below match the INITIAL conditions expected. For
40 % instance, if usetau = x/Kp in Parameters.ini, carefully study the plot of
41 % tau versus time to understand what the initial conditions are. You can
42 % find this plot in the paper.
49 PSD_prof = steady_state(L1d, RadialDiffCoeff(Kp,L1d),tau)';
51 PSD_ongrid.arr = squeeze(repmat(PSD_prof,[1 1 3 3]));
53 save_plt(PSD0_filename,PSD_ongrid);
55 [
PSD] = load_plt(PSD_filename,
'permute',
'nosqueeze',
'skip_zones', 0,
'n_zones', 99);%
61 % ALL PHASE SPACE DENSITY plots are done here.
64 % Top figure is
PSD profile across all L and all times.
66 plot(L.arr(:,2,2),
PSD.arr(:,:,2,2))
68 % This is a pcolor plot representing the same thing as above.
70 time = 1:numel(
PSD.arr(:,1,1,1));
71 pcolor(time,squeeze(L.arr(:,2,2)),squeeze(
PSD.arr(:,:,2,2))')
73 set(
gcf, 'renderer', 'zbuffer');
79 % This is where the actual pcolor plots of Figure 2 and 3 in the paper are
80 % housed. All ELECTRON FLUX plots are done here.
82 J1 = ElectronFlux(L.arr(:,2,2),Kp,tau,
PSD.arr(1,:,2,2)');
83 % J1 is the electron flux given the Phase space
density that we provide
84 J2 = ElectronFlux(L.arr(:,2,2),Kp,tau);
85 % J2 is the electron flux calculated straight with
steady_state.m
87 J = zeros(size(
PSD.arr(:,:,2,2)));
88 time = (4/numel(
PSD.arr(:,1,1,1)))*(1:numel(
PSD.arr(:,1,1,1)));
90 J(idx,:) = ElectronFlux(L.arr(:,2,2),Kp,tau,
PSD.arr(idx,:,2,2)');
95 plot(L.arr(:,2,2),J1, 'b',L.arr(:,2,2),J2, 'r')
97 pcolor(time,squeeze(L.arr(:,2,2)),J')
100 set(gcf, 'renderer', 'zbuffer');
101 %set(gca,'Clim',[0 1e6]);
109 plot(L.arr(:,2,2), J([1 3*24 end],:))
110 set(gca,'yscale','log', 'ylim',[1e0 1e4])
void steady_state(Matrix1D< double > &f, double tau, double Kp, int nx, Matrix1D< double > &L, double f_bnd_out, double f_bnd_in)
functions for write log and support files. Functions are defined in Output.h and descripted in Output...
void gcf(double *gammcf, double a, double x, double *gln)
Returns the incomplete gamma function Q(a, x) evaluated by its continued fraction representation as g...
Phase Space Density class.
Parameters_structure parameters
Parameters structure, with all parameters from the parameters.ini file. The default parameters define...
double density(double L)
Chorus density model.