VERB_code_2.3
Figure_2.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 addpath 'F:\Dmitriy\VERB_Code_2.0_3D_WithMixedTerm_EXAMPLE\Plot\Various_functions';
11 path_arr = [
12  {'..\Calculations\201x201_Mixed_Term_Block_Full_LAPACK\'},
13  {'..\Calculations\201x201_NoMT_Split_Full\'},
14  ];
15 
16 
17 %c_min = 10^-10;
18 %c_max = 2*10^5;
19 
20 %c_min = 10^-4;
21 %c_max = 10^5;
22 n_colors = 100;
23 
24  f1 = figure('Units', 'inches', 'Position', [1 1 8.5 2.4], ...
25  'DefaultAxesFontName', 'times new roman', 'DefaultAxesFontWeight', 'bold', 'DefaultAxesFontSize', 7, 'DefaultTextFontWeight', 'bold', 'DefaultTextFontSize', 7,...
26  'PaperOrientation', 'portrait', 'PaperPositionMode', 'auto', 'DefaultLineLineWidth', 2);%, 'PaperPosition', [1 1 3.8 9]);%
27 
28 %colormap('gray')
29 
30 % xmin = 0.1; ymin = 0.15;
31 % xmax = 0.95; ymax = 0.75;
32  xmin = 0.07; ymin = 0.15;
33  xmax = 0.80; ymax = 0.90;
34  subplot1(1, 2, 'Min', [xmin ymin], 'Max', [xmax ymax], 'Gap', [0.05 0.01], 'YTickL', 'All')
35 
36  hold on
37 
38 % f2 = figure('Units', 'inches', 'Position', [1 1 5 3], ...
39 % 'DefaultAxesFontName', 'times new roman', 'DefaultAxesFontWeight', 'bold', 'DefaultAxesFontSize', 12, 'DefaultTextFontWeight', 'bold', 'DefaultTextFontSize', 12,...
40 % 'PaperOrientation', 'portrait', 'PaperPositionMode', 'auto', 'DefaultLineLineWidth', 2);%, 'PaperPosition', [1 1 3.8 9]);%
41 % hold on
42 
43  col0 = 'black';
44  col1 = 'blue';
45  col2 = [0 0.5 0];
46  col3 = 'red';
47 % col0 = [0 0.0 0];
48 % col1 = [0.2 0.2 0.2];
49 % col2 = [0.4 0.4 0.4];
50 % col3 = [0.6 0.6 0.6];
51 
52 for i = 1:length(path_arr)
53  path1 = path_arr{i};
54 
55 
56  [L, epc, alpha, pc] = load_plt([path1, 'Output/perp_grid.plt'], 'squeeze', 'permute');
57  [PSD] = load_plt([path1, 'Output/OutPSD.dat'], 'permute', 'nosqueeze');%, 'skip_zones', 7);%, 'n_zones', 25);%, 'skip_zones', 10, 'n_zones', 5);
58 
59  c_min = 10^-9;
60  c_max = 10^0;
61 
62 
63  alpha_2d = alpha.arr.*180./pi;
64  epc_2d = epc.arr;
65 
66  time0=1;
67  [dif, time8]=min(abs(PSD.time - 0.1));
68  [dif, time24]=min(abs(PSD.time - 1.0));
69  [dif, time48]=min(abs(PSD.time - 2.0));
70  PSD.time([time8,time24,time48])*24
71 
72  if (i == 1) line = ''; end;
73  if (i == 2) line = '--'; end;
74 
75 % figure(f1);
76  subplot1(1);
77  pos1 = get(gca,'Position');
78 
79  hold on
80  %energy_slice = max(find(epc_2d(:,1) <= 0.5))
81  %if (i == 2) energy_slice = energy_slice + 1; end;
82  [tmp, energy_slice] = min(abs(epc_2d(:,1) - 0.5));
83 
84  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
85  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);
86  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);
87  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);
88  set(gca, 'yscale', 'log');
89  xlim([0 90])
90  ylim([3e-3 6e-1])
91  xlabel('Pitch angle, deg');
92  ylabel('Flux, arbitrary units');
93  grid on
94  set(gca, 'yminorgrid', 'off')
95 % text(max(xlim), min(ylim), [num2str(round(epc.arr(energy_slice, 1)*10)/10), 'MeV'], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'Color', 'black');
96  text(min(xlim), max(ylim), [' (A)'], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'FontWeight', 'bold', 'Color', 'black');
97  title([num2str(round(epc.arr(energy_slice, 1)*10)/10), 'MeV'])
98 
99  [tmp, alpha30] = min(abs(alpha_2d(energy_slice,:) - 30));
100  [tmp, alpha60] = min(abs(alpha_2d(energy_slice,:) - 60));
101  [tmp, alpha90] = min(abs(alpha_2d(energy_slice,:) - 90));
102  if (i == 1)
103  F1_1_30 = PSD.arr(time48, 1, energy_slice, alpha30);
104  F1_1_60 = PSD.arr(time48, 1, energy_slice, alpha60);
105  F1_1_90 = PSD.arr(time48, 1, energy_slice, alpha90);
106  end;
107  if (i == 2)
108  F2_1_30 = PSD.arr(time48, 1, energy_slice, alpha30);
109  F2_1_60 = PSD.arr(time48, 1, energy_slice, alpha60);
110  F2_1_90 = PSD.arr(time48, 1, energy_slice, alpha90);
111  end
112 
113 % figure(f2);
114  subplot1(2);
115 
116  hold on
117  %energy_slice = max(find(epc_2d(:,1) <= 2))
118  %if (i == 2) energy_slice = energy_slice + 1; end;
119  [tmp, energy_slice] = min(abs(epc_2d(:,1) - 2));
120 
121  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
122  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);
123  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);
124  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);
125  set(gca, 'yscale', 'log');
126  xlim([0 90])
127  ylim([1e-10 2e-4])
128  xlabel('Pitch angle, deg');
129 % ylabel('Flux, arbitrary units');
130  grid on
131 % text(max(xlim), min(ylim), [num2str(round(epc.arr(energy_slice, 1)*10)/10), 'MeV'], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'Color', 'black');
132  text(min(xlim), max(ylim), [' (B)'], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'FontWeight', 'bold', 'Color', 'black');
133  title([num2str(round(epc.arr(energy_slice, 1)*10)/10), 'MeV']);
134 
135  [tmp, alpha30] = min(abs(alpha_2d(energy_slice,:) - 30));
136  [tmp, alpha60] = min(abs(alpha_2d(energy_slice,:) - 60));
137  [tmp, alpha90] = min(abs(alpha_2d(energy_slice,:) - 90));
138  if (i == 1)
139  F1_2_30 = PSD.arr(time48, 1, energy_slice, alpha30);
140  F1_2_60 = PSD.arr(time48, 1, energy_slice, alpha60);
141  F1_2_90 = PSD.arr(time48, 1, energy_slice, alpha90);
142  end;
143  if (i == 2)
144  F2_2_30 = PSD.arr(time48, 1, energy_slice, alpha30);
145  F2_2_60 = PSD.arr(time48, 1, energy_slice, alpha60);
146  F2_2_90 = PSD.arr(time48, 1, energy_slice, alpha90);
147  end
148 
149 end
150 
151 F2_1_30/F1_1_30
152 F2_1_60/F1_1_60
153 F2_1_90/F1_1_90
154 
155 F2_2_30/F1_2_30
156 F2_2_60/F1_2_60
157 F2_2_90/F1_2_90
158 
159 %lg = legend('Initial', [num2str(PSD.time(time8)*24), ' hrs, d/o'], [num2str(PSD.time(time24)*24), ' hrs, d/o'], [num2str(PSD.time(time48)*24), ' hrs, d/o'], [num2str(PSD.time(time8)*24), ' hrs, m.t.'], [num2str(PSD.time(time24)*24), ' hrs, m.t.'], [num2str(PSD.time(time48)*24), ' hrs, m.t.'], 'Location', [xmin ymax+0.1 (xmax-xmin) 0.04], 'Orientation', 'Horizontal', 0);
160 %set(lg, 'Box', 'off')
161 
162 lg = legend('Initial', ...
163  [num2str(PSD.time(time8)*24), ' hrs, mixed+diag.'], [num2str(PSD.time(time24)*24), ' hrs, mixed+diag.'], [num2str(PSD.time(time48)*24), ' hrs, mixed+diag.'], ...
164  [num2str(PSD.time(time8)*24), ' hrs, diag. only'], [num2str(PSD.time(time24)*24), ' hrs, diag. only'], [num2str(PSD.time(time48)*24), ' hrs, diag. only'], ...
165  'Location', [xmax ymin (0.95-xmax) (ymax-ymin)], 'Orientation', 'Vertical', 0);
166 
167  save_picture(gcf, ['Figure_2_bw'], 'jpg', 'pdf', 'epsc');
168 
double max(double v1, double v2)
Return maximum.
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