VERB_code_2.3
create_PSD0.m
1 %%
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
6 % README.md
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.
10 
11 
12 %%
13 
14 % Section 1
15 
16 user = 2;
17 
18 switch(user)
19  % User is Akshat
20  case 1
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';
24  % User is Alec
25  case 2
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';
29 end
30 
31 [L, epc, alpha, pc] = load_plt(Grid_filename, 'squeeze', 'permute');
32 
33 %%
34 
35 % Section 2
36 
37 % Run this section after Parameters.ini has been modified and VERB code has been run
38 
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.
43 
44 tau = 3;
45 Kp = 2;
46 
47 L1d = L.arr(:,1,1);
48 
49 PSD_prof = steady_state(L1d, RadialDiffCoeff(Kp,L1d),tau)';
50 
51 PSD_ongrid.arr = squeeze(repmat(PSD_prof,[1 1 3 3]));
52 
53 save_plt(PSD0_filename,PSD_ongrid);
54 
55 [PSD] = load_plt(PSD_filename, 'permute', 'nosqueeze', 'skip_zones', 0, 'n_zones', 99);%
56 
57 %%
58 
59 %Section 3
60 
61 % ALL PHASE SPACE DENSITY plots are done here.
62 
63 figure
64 % Top figure is PSD profile across all L and all times.
65 subplot(2,1,1)
66 plot(L.arr(:,2,2),PSD.arr(:,:,2,2))
67 
68 % This is a pcolor plot representing the same thing as above.
69 subplot(2,1,2)
70 time = 1:numel(PSD.arr(:,1,1,1));
71 pcolor(time,squeeze(L.arr(:,2,2)),squeeze(PSD.arr(:,:,2,2))')
72 colorbar
73 set(gcf, 'renderer', 'zbuffer');
74 
75 %%
76 
77 %Section 4
78 
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.
81 
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
86 
87 J = zeros(size(PSD.arr(:,:,2,2)));
88 time = (4/numel(PSD.arr(:,1,1,1)))*(1:numel(PSD.arr(:,1,1,1)));
89 for idx=1:numel(time)
90  J(idx,:) = ElectronFlux(L.arr(:,2,2),Kp,tau,PSD.arr(idx,:,2,2)');
91 end
92 
93 figure
94 subplot(2,1,1)
95 plot(L.arr(:,2,2),J1, 'b',L.arr(:,2,2),J2, 'r')
96 subplot(2,1,2)
97 pcolor(time,squeeze(L.arr(:,2,2)),J')
98 colorbar
99 ylim([1,6])
100 set(gcf, 'renderer', 'zbuffer');
101 %set(gca,'Clim',[0 1e6]);
102 
103 
104 %%
105 
106 % Section 5
107 
108 figure;
109 plot(L.arr(:,2,2), J([1 3*24 end],:))
110 set(gca,'yscale','log', 'ylim',[1e0 1e4])
111 
112 %%
void steady_state(Matrix1D< double > &f, double tau, double Kp, int nx, Matrix1D< double > &L, double f_bnd_out, double f_bnd_in)
Definition: PSD.cpp:2373
functions for write log and support files. Functions are defined in Output.h and descripted in Output...
Definition: Output.cpp:15
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...
Definition: erf.cpp:103
Phase Space Density class.
Definition: PSD.h:48
Parameters_structure parameters
Parameters structure, with all parameters from the parameters.ini file. The default parameters define...
Definition: Main.cpp:185
double density(double L)
Chorus density model.