VERB_code_2.3
Figure_3.m
1 try
2  mov = close(mov);
3 catch
4 end
5 
6 clear all
7 
8 save_movie = false;
9 movie_name = 'mt_fix3';
10 
11 %c_min = 10^-10;
12 %c_max = 2*10^5;
13 
14 %c_min = 10^-4;
15 %c_max = 10^5;
16 n_colors = 100;
17 
18  f1 = figure('Units', 'inches', 'Position', [1 1 8.5 4.6], ...
19  'DefaultAxesFontName', 'times new roman', 'DefaultAxesFontWeight', 'bold', 'DefaultAxesFontSize', 7, 'DefaultTextFontWeight', 'bold', 'DefaultTextFontSize', 7,...
20  'PaperOrientation', 'portrait', 'PaperPositionMode', 'auto', 'DefaultLineLineWidth', 2);%, 'PaperPosition', [1 1 3.8 9]);%
21 
22 
23  xmin = 0.07; ymin = 0.15;
24  xmax = 0.80; ymax = 0.95;
25  subplot1(2, 2, 'Min', [xmin ymin], 'Max', [xmax ymax], 'Gap', [0.05 0.01], 'YTickL', 'All')
26 
27  hold on
28 
29 
30 path_arr1 = [
31 
32  {'..\Calculations\201x201 NoMT Split Full V.1.28\'},
33  {'..\Calculations\201x201 NoMT Split night V.1.28\'},
34 
35  ];
36 
37 path_arr2 = [
38 
39  {'..\Calculations\201x201 MT Block Full V.1.28 LAPACK\'},
40  {'..\Calculations\201x201 MT Block night V.1.28 LAPACK\'},
41 
42  ];
43 
44 % col0 = 'black';
45 % col1 = 'blue';
46 % col2 = [0 0.5 0];
47 % col3 = 'red';
48  col0 = [0 0.0 0];
49  col1 = [0.2 0.2 0.2];
50  col2 = [0.4 0.4 0.4];
51  col3 = [0.6 0.6 0.6];
52 
53 
54 for i = 1:length(path_arr1)
55  path1 = path_arr1{i};
56 
57 
58  [L, epc, alpha, pc] = load_plt([path1, 'Output/perp_grid.plt'], 'squeeze', 'permute');
59  [PSD] = load_plt([path1, 'Output/OutPSD.dat'], 'permute', 'nosqueeze');%, 'skip_zones', 7);%, 'n_zones', 25);%, 'skip_zones', 10, 'n_zones', 5);
60 
61  c_min = 10^-9;
62  c_max = 10^0;
63 
64 
65  alpha_2d = alpha.arr.*180./pi;
66  epc_2d = epc.arr;
67 
68  time0=1;
69  [dif, time8]=min(abs(PSD.time - 0.1));
70  [dif, time24]=min(abs(PSD.time - 1.0));
71  [dif, time48]=min(abs(PSD.time - 2.0));
72  PSD.time([time8,time24,time48])*24
73 
74  if (i == 1) line = ''; end;
75  if (i == 2) line = '--'; end;
76 
77 % figure(f1);
78  subplot1(1);
79 
80  hold on
81  %energy_slice = max(find(epc_2d(:,1) <= 0.5))
82  %if (i == 2) energy_slice = energy_slice + 1; end;
83  [tmp, energy_slice] = min(abs(epc_2d(:,1) - 0.5));
84 
85  if (i == 1) plot(alpha_2d(energy_slice,:), squeeze(PSD.arr(time0, 1, energy_slice, :))'.*pc.arr(energy_slice, :).*pc.arr(energy_slice, :), 'LineWidth', 2, 'Color', col0); end;
86  plot(alpha_2d(energy_slice,:), squeeze(PSD.arr(time8, 1, energy_slice, :))'.*pc.arr(energy_slice, :).*pc.arr(energy_slice, :), line, 'LineWidth', 2, 'Color', col1);
87  plot(alpha_2d(energy_slice,:), squeeze(PSD.arr(time24, 1, energy_slice, :))'.*pc.arr(energy_slice, :).*pc.arr(energy_slice, :), line, 'LineWidth', 2, 'Color', col2);
88  plot(alpha_2d(energy_slice,:), squeeze(PSD.arr(time48, 1, energy_slice, :))'.*pc.arr(energy_slice, :).*pc.arr(energy_slice, :), line, 'LineWidth', 2, 'Color', col3);
89  set(gca, 'yscale', 'log');
90  xlim([0 90])
91  ylim([3e-3 6e-1])
92  xlabel('Pitch angle, deg');
93  ylabel({'Flux, arbitrary units, ', [num2str(round(epc.arr(energy_slice, 1)*10)/10), 'MeV']});
94  grid on
95  set(gca, 'yminorgrid', 'off')
96 % text(max(xlim), min(ylim), [num2str(round(epc.arr(energy_slice, 1)*10)/10), 'MeV'], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'Color', 'black');
97  text(min(xlim), max(ylim), [' (A) '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'FontWeight', 'bold', 'Color', 'black');
98  text(max(xlim), min(ylim), [{' Diagonal only '}, {[num2str(round(epc.arr(energy_slice, 1)*10)/10), ' MeV ']}], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'Color', 'black');
99  title('Diagonal only');
100 
101 %% save_picture(gcf, [path1, 'w_wo,epc=', num2str(epc.arr(energy_slice, 1))], 'jpg', 'pdf', 'eps');
102 
103 % figure(f2);
104  subplot1(3);
105 
106  hold on
107  %energy_slice = max(find(epc_2d(:,1) <= 2))
108  %if (i == 2) energy_slice = energy_slice + 1; end;
109  [tmp, energy_slice] = min(abs(epc_2d(:,1) - 2));
110 
111  if (i == 1) plot(alpha_2d(energy_slice,:), squeeze(PSD.arr(time0, 1, energy_slice, :))'.*pc.arr(energy_slice, :).*pc.arr(energy_slice, :), 'LineWidth', 2, 'Color', col0); end;
112  plot(alpha_2d(energy_slice,:), squeeze(PSD.arr(time8, 1, energy_slice, :))'.*pc.arr(energy_slice, :).*pc.arr(energy_slice, :), line, 'LineWidth', 2, 'Color', col1);
113  plot(alpha_2d(energy_slice,:), squeeze(PSD.arr(time24, 1, energy_slice, :))'.*pc.arr(energy_slice, :).*pc.arr(energy_slice, :), line, 'LineWidth', 2, 'Color', col2);
114  plot(alpha_2d(energy_slice,:), squeeze(PSD.arr(time48, 1, energy_slice, :))'.*pc.arr(energy_slice, :).*pc.arr(energy_slice, :), line, 'LineWidth', 2, 'Color', col3);
115  set(gca, 'yscale', 'log');
116  xlim([0 90])
117 % ylim([1e-10 2e-4])
118  ylim([1e-10 1e-2])
119  xlabel('Pitch angle, deg');
120  ylabel({'Flux, arbitrary units, ', [num2str(round(epc.arr(energy_slice, 1)*10)/10), 'MeV']});
121  grid on
122 % text(max(xlim), min(ylim), [num2str(round(epc.arr(energy_slice, 1)*10)/10), 'MeV'], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'Color', 'black');
123  text(min(xlim), max(ylim), [' (C) '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'FontWeight', 'bold', 'Color', 'black');
124  text(max(xlim), min(ylim), [{' Diagonal only '}, {[num2str(round(epc.arr(energy_slice, 1)*10)/10), ' MeV ']}], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'Color', 'black');
125 % title([num2str(round(epc.arr(energy_slice, 1)*10)/10), 'MeV']);
126 
127 %% save_picture(gcf, [path1, 'w_wo,epc=', num2str(epc.arr(energy_slice, 1))], 'jpg', 'pdf', 'eps');
128 
129 end
130 
131 
132 lg = legend('Initial', [num2str(PSD.time(time8)*24), ' hrs, day+night'], [num2str(PSD.time(time24)*24), ' hrs, day+night'], [num2str(PSD.time(time48)*24), ' hrs, day+night'], ...
133  [num2str(PSD.time(time8)*24), ' hrs, 2 \times night side'], [num2str(PSD.time(time24)*24), ' hrs, 2 \times night side'], [num2str(PSD.time(time48)*24), ' hrs, 2 \times night side'], ...
134  'Location', [xmax ymin+(ymax-ymin)/2 (0.95-xmax) (ymax-ymin)/2-0.005], 'Orientation', 'Vertical', 0);
135 
136 
137 
138 
139 for i = 1:length(path_arr2)
140  path1 = path_arr2{i};
141 
142 
143  [L, epc, alpha, pc] = load_plt([path1, 'Output/perp_grid.plt'], 'squeeze', 'permute');
144  [PSD] = load_plt([path1, 'Output/OutPSD.dat'], 'permute', 'nosqueeze');%, 'skip_zones', 7);%, 'n_zones', 25);%, 'skip_zones', 10, 'n_zones', 5);
145 
146  c_min = 10^-9;
147  c_max = 10^0;
148 
149 
150  alpha_2d = alpha.arr.*180./pi;
151  epc_2d = epc.arr;
152 
153  time0=1;
154  [dif, time8]=min(abs(PSD.time - 0.1));
155  [dif, time24]=min(abs(PSD.time - 1.0));
156  [dif, time48]=min(abs(PSD.time - 2.0));
157  PSD.time([time8,time24,time48])*24
158 
159  if (i == 1) line = ''; end;
160  if (i == 2) line = '--'; end;
161 
162 % figure(f1);
163  subplot1(2);
164 
165  hold on
166  %energy_slice = max(find(epc_2d(:,1) <= 0.5))
167  %if (i == 2) energy_slice = energy_slice + 1; end;
168  [tmp, energy_slice] = min(abs(epc_2d(:,1) - 0.5));
169 
170  plot(alpha_2d(energy_slice,:), squeeze(PSD.arr(time0, 1, energy_slice, :))'.*pc.arr(energy_slice, :).*pc.arr(energy_slice, :), 'LineWidth', 2, 'Color', col0);
171  plot(alpha_2d(energy_slice,:), squeeze(PSD.arr(time8, 1, energy_slice, :))'.*pc.arr(energy_slice, :).*pc.arr(energy_slice, :), line, 'LineWidth', 2, 'Color', col1);
172  plot(alpha_2d(energy_slice,:), squeeze(PSD.arr(time24, 1, energy_slice, :))'.*pc.arr(energy_slice, :).*pc.arr(energy_slice, :), line, 'LineWidth', 2, 'Color', col2);
173  plot(alpha_2d(energy_slice,:), squeeze(PSD.arr(time48, 1, energy_slice, :))'.*pc.arr(energy_slice, :).*pc.arr(energy_slice, :), line, 'LineWidth', 2, 'Color', col3);
174  set(gca, 'yscale', 'log');
175  xlim([0 90])
176  ylim([3e-3 6e-1])
177  xlabel('Pitch angle, deg');
178 % ylabel('Flux, arbitrary units');
179  grid on
180  set(gca, 'yminorgrid', 'off')
181 % text(max(xlim), min(ylim), [num2str(round(epc.arr(energy_slice, 1)*10)/10), 'MeV'], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'Color', 'black');
182  text(min(xlim), max(ylim), [' (B) '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'FontWeight', 'bold', 'Color', 'black');
183  text(max(xlim), min(ylim), [{' Diagonal + Mixed '}, {[num2str(round(epc.arr(energy_slice, 1)*10)/10), ' MeV ']}], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'Color', 'black');
184  title('Diagonal + Mixed');
185 
186 %% save_picture(gcf, [path1, 'w_wo,epc=', num2str(epc.arr(energy_slice, 1))], 'jpg', 'pdf', 'eps');
187 
188 % figure(f2);
189  subplot1(4);
190 
191  hold on
192  %energy_slice = max(find(epc_2d(:,1) <= 2))
193  %if (i == 2) energy_slice = energy_slice + 1; end;
194  [tmp, energy_slice] = min(abs(epc_2d(:,1) - 2));
195 
196  plot(alpha_2d(energy_slice,:), squeeze(PSD.arr(time0, 1, energy_slice, :))'.*pc.arr(energy_slice, :).*pc.arr(energy_slice, :), 'LineWidth', 2, 'Color', col0);
197  plot(alpha_2d(energy_slice,:), squeeze(PSD.arr(time8, 1, energy_slice, :))'.*pc.arr(energy_slice, :).*pc.arr(energy_slice, :), line, 'LineWidth', 2, 'Color', col1);
198  plot(alpha_2d(energy_slice,:), squeeze(PSD.arr(time24, 1, energy_slice, :))'.*pc.arr(energy_slice, :).*pc.arr(energy_slice, :), line, 'LineWidth', 2, 'Color', col2);
199  plot(alpha_2d(energy_slice,:), squeeze(PSD.arr(time48, 1, energy_slice, :))'.*pc.arr(energy_slice, :).*pc.arr(energy_slice, :), line, 'LineWidth', 2, 'Color', col3);
200  set(gca, 'yscale', 'log');
201  xlim([0 90])
202 % ylim([1e-10 2e-4])
203  ylim([1e-10 1e-2])
204  xlabel('Pitch angle, deg');
205 % ylabel('Flux, arbitrary units');
206  grid on
207 % text(max(xlim), min(ylim), [num2str(round(epc.arr(energy_slice, 1)*10)/10), 'MeV'], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'Color', 'black');
208  text(min(xlim), max(ylim), [' (D) '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'FontWeight', 'bold', 'Color', 'black');
209  text(max(xlim), min(ylim), [{' Diagonal + Mixed '}, {[num2str(round(epc.arr(energy_slice, 1)*10)/10), ' MeV ']}], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'Color', 'black');
210 
211 %% save_picture(gcf, [path1, 'w_wo,epc=', num2str(epc.arr(energy_slice, 1))], 'jpg', 'pdf', 'eps');
212 
213 end
214 
215 
216 %save_picture(gcf, [path1, '../Figures/Figure_3'], 'jpg', 'pdf', 'epsc');
217 save_picture(gcf, ['Figure_3_bw'], 'jpg', 'pdf', 'epsc');
218 
double max(double v1, double v2)
Return maximum.
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
Matrix3D< double > arr
array of PSD values
Definition: PSD.h:53
Phase Space Density class.
Definition: PSD.h:48