3 % PLOT will, of course, plot your results. The only parameter that needs
4 % to be set is SAVEFIG, which determines whether to save your plot to the
7 % Remember, any existing Figures/x.png plot will be replaced!
10 % choose whether to save the output plot
16 addpath(
'../Setup/Various_functions');
18 if(~exist(
'FIGURE',
'var'))
22 % Resize figure based on how many plots it has
23 if(~strcmp(FIGURE,
'3') && ~strcmp(FIGURE,
'1') && ~strcmp(FIGURE,
'4'))
24 figure(
'Position',[10 200 900 300]);
26 figure(
'Position',[10 200 900 600]);
32 % This code will reproduce our diffusion coefficients vs. L-star plots.
36 plot(L, 1./RadialDiffCoeff(Kp,L))
37 set(gca,'ylim',[1e-2 1e6],'xlim',[2 8],'yscale','log');
39 plot(L,ones(size(L)),'r')
40 ylabel('1/DLL (Days/R_E^2)');
43 % This code will reproduce our phase space
density vs. L-star plots
47 DLL = RadialDiffCoeff(Kp, L);
50 set(gca,'ylim',[1e-2 1e1],'xlim',[2 8],'yscale','log');
51 ylabel('Phase Space Density');
54 % This code will reproduce our electron flux vs. L-star plots
59 plot(L, ElectronFlux(L, Kp, tau))
60 set(gca,'ylim',[10e-1 10e4],'xlim',[2 8],'yscale','log');
61 ylabel('Electron Flux');
65 if(strcmp(FIGURE,'1'))
67 saveas(
gcf,['../Figures/1.png']);
74 PSD_filename = ['../Execute/Output_' FIGURE '/OutPSD_rad.dat'];
75 Grid_filename = ['../Execute/Output_' FIGURE '/rad_grid.plt'];
78 [L, epc, alpha, pc] = load_plt(Grid_filename, 'squeeze', 'permute');
79 [
PSD] = load_plt(PSD_filename, 'permute', 'nosqueeze', 'skip_zones', 0, 'n_zones', 99);%
82 % This is where the actual pcolor plots of Figure 2 and 3 in the paper are
83 % housed. All ELECTRON FLUX plots are done here.
84 J1 = ElectronFlux(L.arr(:,2,2),Kp,tau,
PSD.arr(1,:,2,2)');
85 % J1 is the electron flux given the Phase space
density that we provide
86 J2 = ElectronFlux(L.arr(:,2,2),Kp,tau);
87 % J2 is the electron flux calculated straight with
steady_state.m
90 J = zeros(size(
PSD.arr(:,:,2,2)));
91 time = (4/numel(
PSD.arr(:,1,1,1)))*(1:numel(
PSD.arr(:,1,1,1)));
93 J(idx,:) = ElectronFlux(L.arr(:,2,2),Kp,tau,
PSD.arr(idx,:,2,2)');
96 % Produce subplots according to
# of plots in the figure
97 if(strcmp(FIGURE,
'3'))
99 elseif(strcmp(FIGURE,
'4'))
105 %plot(L.arr(:,2,2),J1, 'b',L.arr(:,2,2),J2, 'r')
107 plot(L.arr(:,2,2), J([1 3*24 end],:))
108 set(gca,'yscale','log', 'ylim',[1e0 1e4])
109 ylabel('Electron Flux');
112 if(strcmp(FIGURE,'3'))
114 elseif(strcmp(FIGURE,'4'))
120 pcolor(time,squeeze(L.arr(:,2,2)),J')
121 set(gca,'ylim',[1 6]);
123 xlabel('Time, days');
129 % Phase Space Density Plot (Middle Row)
130 set(gca,'Clim',[0 1e3]);
132 plot(L.arr(:,2,2),
PSD.arr(:,:,2,2))
133 ylabel('Phase Space Density');
136 % This is a pcolor plot representing the same thing as above
138 time = 1:numel(
PSD.arr(:,1,1,1));
139 pcolor(time,squeeze(L.arr(:,2,2)),squeeze(
PSD.arr(:,:,2,2))')
142 set(gca,'ylim',[1 6]);
144 xlabel('Time, days');
146 % Plot of time-dependent boundary conditions (Bottom Row)
147 x=[0 0.5 1 1.5 2 4 999];
148 y=[1 1 0.1 1 1 1 1 ];
151 set(gca,'xlim',[0 4], 'ylim', [0 1]);
152 ylabel('Boundary Fluxes');
153 xlabel('Time, days');
155 % Plot of time-dependent lifetimes
160 set(gca,'xlim',[0 4], 'ylim', [0 3]);
162 xlabel('Time, days');
163 case '3' % Second row of Figure 3
168 set(gca,'ylim',[0 6]);
170 xlabel('Time, days');
172 xtau = 6; % the VERB code use forumla for tau(time) = xtau / Kp.
176 set(gca,'ylim',[0 3]);
177 ylabel('\tau_L, days');
178 xlabel('Time, days');
181 if(~strcmp(FIGURE,'3') && ~strcmp(FIGURE,'4'))
184 set(gcf, 'renderer', 'zbuffer');
187 saveas(gcf,['../Figures/' FIGURE '.png']);
193 % To reproduce the electron flux pcolor plots in Figures 2:
194 % 1. Use create_PSD0.m. Instructions are clearly outlined at beginning.
196 % Some additional information for Parameters.ini:
197 % To reproduce row 1 of Figure 2, set:
202 % Further, in create_PSD0.m, set Kp = 2, tau = 1 to match the initial
205 % To reproduce row 2 of Figure 2, set:
210 % Further, in create_PSD0.m, set Kp = 2, tau = 10 to match the initial
213 % To reproduce row 1 of Figure 3, set:
218 % Further, in create_PSD0.m, set Kp = 2, tau = 3 to match the initial
221 % To reproduce row 1 of Figure 4, set:
226 % Further, in create_PSD0.m, set Kp = 2, tau = 3 to match the initial
229 %% THE FOLLOWING CANNOT BE REPRODUCED
231 % To reproduce row 1, first figure, of Figure 5, set:
236 % Edit kp.dat file in Input folder so that at t = 1, kp = 60
237 % Further, in create_PSD0.m, set Kp = 2, tau = 3 to match the initial
240 % To reproduce row 1, second figure, of Figure 5, set:
245 % Edit kp.dat file in Input folder so that at t = 1, kp = 12
246 % Further, in create_PSD0.m, set Kp = 2, tau = 3 to match the initial
void steady_state(Matrix1D< double > &f, double tau, double Kp, int nx, Matrix1D< double > &L, double f_bnd_out, double f_bnd_in)
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.
double density(double L)
Chorus density model.