VERB_code_2.3
Plot.m
1 %% Plot
2 %
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
5 % Figures folder.
6 %
7 % Remember, any existing Figures/x.png plot will be replaced!
8 
9 %% Parameter
10 % choose whether to save the output plot
11 SAVEFIG = false;
12 
13 %% Plotting code
14 
15 addpath('../Setup/');
16 addpath('../Setup/Various_functions');
17 
18 if(~exist('FIGURE','var'))
19  FIGURE = '4';
20 end
21 
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]);
25 else
26  figure('Position',[10 200 900 600]);
27 end
28 
29 switch (FIGURE)
30  % Figure 1 plots
31  case '1'
32  % This code will reproduce our diffusion coefficients vs. L-star plots.
33  subplot(2,2,1);
34  L = 2:0.1:8;
35  Kp = 6;
36  plot(L, 1./RadialDiffCoeff(Kp,L))
37  set(gca,'ylim',[1e-2 1e6],'xlim',[2 8],'yscale','log');
38  hold on
39  plot(L,ones(size(L)),'r')
40  ylabel('1/DLL (Days/R_E^2)');
41  xlabel('L value');
42 
43  % This code will reproduce our phase space density vs. L-star plots
44  subplot(2,2,2);
45  L = 2:0.1:8;
46  Kp = 6;
47  DLL = RadialDiffCoeff(Kp, L);
48  tau = 10;
49  plot(L, steady_state(L, DLL, tau))
50  set(gca,'ylim',[1e-2 1e1],'xlim',[2 8],'yscale','log');
51  ylabel('Phase Space Density');
52  xlabel('L value');
53 
54  % This code will reproduce our electron flux vs. L-star plots
55  subplot(2,2,3);
56  L = [2:0.1:8]';
57  Kp = 6;
58  tau = 10;
59  plot(L, ElectronFlux(L, Kp, tau))
60  set(gca,'ylim',[10e-1 10e4],'xlim',[2 8],'yscale','log');
61  ylabel('Electron Flux');
62  xlabel('L value');
63 end
64 
65 if(strcmp(FIGURE,'1'))
66  if(SAVEFIG)
67  saveas(gcf,['../Figures/1.png']);
68  close(gcf);
69  end
70  return;
71 end
72 
73 % Figures 2 - 4
74 PSD_filename = ['../Execute/Output_' FIGURE '/OutPSD_rad.dat'];
75 Grid_filename = ['../Execute/Output_' FIGURE '/rad_grid.plt'];
76 
77 
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);%
80 
81 %Section 4
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
88 
89 % Generate time scale
90 J = zeros(size(PSD.arr(:,:,2,2)));
91 time = (4/numel(PSD.arr(:,1,1,1)))*(1:numel(PSD.arr(:,1,1,1)));
92 for idx=1:numel(time)
93  J(idx,:) = ElectronFlux(L.arr(:,2,2),Kp,tau,PSD.arr(idx,:,2,2)');
94 end
95 
96 % Produce subplots according to # of plots in the figure
97 if(strcmp(FIGURE,'3'))
98  subplot(2,2,1);
99 elseif(strcmp(FIGURE,'4'))
100  subplot(3,2,1);
101 else
102  subplot(1,2,1);
103 end
104 
105 %plot(L.arr(:,2,2),J1, 'b',L.arr(:,2,2),J2, 'r')
106 
107 plot(L.arr(:,2,2), J([1 3*24 end],:))
108 set(gca,'yscale','log', 'ylim',[1e0 1e4])
109 ylabel('Electron Flux');
110 xlabel('L');
111 
112 if(strcmp(FIGURE,'3'))
113  subplot(2,2,2);
114 elseif(strcmp(FIGURE,'4'))
115  subplot(3,2,2);
116 else
117  subplot(1,2,2);
118 end
119 
120 pcolor(time,squeeze(L.arr(:,2,2)),J')
121 set(gca,'ylim',[1 6]);
122 ylabel('L');
123 xlabel('Time, days');
124 colorbar
125 colormap jet
126 
127 switch (FIGURE)
128  case '4'
129  % Phase Space Density Plot (Middle Row)
130  set(gca,'Clim',[0 1e3]);
131  subplot(3,2,3)
132  plot(L.arr(:,2,2),PSD.arr(:,:,2,2))
133  ylabel('Phase Space Density');
134  xlabel('L');
135 
136  % This is a pcolor plot representing the same thing as above
137  subplot(3,2,4)
138  time = 1:numel(PSD.arr(:,1,1,1));
139  pcolor(time,squeeze(L.arr(:,2,2)),squeeze(PSD.arr(:,:,2,2))')
140  colorbar
141  colormap jet
142  set(gca,'ylim',[1 6]);
143  ylabel('L');
144  xlabel('Time, days');
145 
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 ];
149  subplot(3,2,5)
150  plot(x,y);
151  set(gca,'xlim',[0 4], 'ylim', [0 1]);
152  ylabel('Boundary Fluxes');
153  xlabel('Time, days');
154 
155  % Plot of time-dependent lifetimes
156  x=[0 1 4 999];
157  y=[3 1 3 3 ];
158  subplot(3,2,6)
159  plot(x,y)
160  set(gca,'xlim',[0 4], 'ylim', [0 3]);
161  ylabel('Lifetimes');
162  xlabel('Time, days');
163  case '3' % Second row of Figure 3
164  x=[0,1,4];
165  y=[2,6,2];
166  subplot(2,2,3);
167  plot(x,y);
168  set(gca,'ylim',[0 6]);
169  ylabel('K_p');
170  xlabel('Time, days');
171 
172  xtau = 6; % the VERB code use forumla for tau(time) = xtau / Kp.
173  z=xtau./y;
174  subplot(2,2,4);
175  plot(x,z);
176  set(gca,'ylim',[0 3]);
177  ylabel('\tau_L, days');
178  xlabel('Time, days');
179 end
180 
181 if(~strcmp(FIGURE,'3') && ~strcmp(FIGURE,'4'))
182 ylim([1,6])
183 end
184 set(gcf, 'renderer', 'zbuffer');
185 
186 if(SAVEFIG)
187 saveas(gcf,['../Figures/' FIGURE '.png']);
188 close(gcf);
189 end
190 
191 %%
192 
193 % To reproduce the electron flux pcolor plots in Figures 2:
194 % 1. Use create_PSD0.m. Instructions are clearly outlined at beginning.
195 
196 % Some additional information for Parameters.ini:
197 % To reproduce row 1 of Figure 2, set:
198 % usetau = constant,
199 % tau = 1,
200 % useKp = file,
201 % useBf = constant
202 % Further, in create_PSD0.m, set Kp = 2, tau = 1 to match the initial
203 % conditions.
204 %
205 % To reproduce row 2 of Figure 2, set:
206 % usetau = constant,
207 % tau = 10,
208 % useKp = file,
209 % useBf = constant
210 % Further, in create_PSD0.m, set Kp = 2, tau = 10 to match the initial
211 % conditions.
212 %
213 % To reproduce row 1 of Figure 3, set:
214 % usetau = x/Kp,
215 % tau = 6,
216 % useKp = file,
217 % useBf = constant,
218 % Further, in create_PSD0.m, set Kp = 2, tau = 3 to match the initial
219 % conditions.
220 %
221 % To reproduce row 1 of Figure 4, set:
222 % usetau = x/Kp,
223 % tau = 6,
224 % useKp = file,
225 % useBf = file,
226 % Further, in create_PSD0.m, set Kp = 2, tau = 3 to match the initial
227 % conditions.
228 
229 %% THE FOLLOWING CANNOT BE REPRODUCED
230 
231 % To reproduce row 1, first figure, of Figure 5, set:
232 % usetau = x/Kp,
233 % tau = 6,
234 % useKp = file,
235 % useBf = file,
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
238 % conditions.
239 %
240 % To reproduce row 1, second figure, of Figure 5, set:
241 % usetau = x/Kp,
242 % tau = 6,
243 % useKp = file,
244 % useBf = file,
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
247 % conditions.
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
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
double density(double L)
Chorus density model.