VERB_code_2.3
Figure_5.m
1 clear all
2 
3 f1 = figure('Units', 'inches', 'Position', [0.5 0.5 8.5 4], ...
4  'DefaultAxesFontName', 'times new roman', 'DefaultAxesFontWeight', 'bold', 'DefaultAxesFontSize', 7, 'DefaultTextFontWeight', 'bold', 'DefaultTextFontSize', 7,...
5  'PaperOrientation', 'portrait', 'PaperPositionMode', 'auto', 'DefaultLineLineWidth', 2);%, 'PaperPosition', [1 1 3.8 9]);%
6 
7 %colormap('gray')
8 
9 cmap = colormap;
10 cmap(1,:);
11 cmap(1,:) = [0 0 0];
12 colormap(cmap)
13 
14 subplot1(2, 3, 'Max', [0.97 0.95], 'Min', [0.05 0.1], 'Gap', [0.02 0.02]);
15 
16 
17 
18 path_arr1 = [
19  {'../Execute/3D_Without_0.1hr v1.28/Output/', ' No mixed terms'},
20  {'../Execute/3D_With_0.1hr v1.28/Output/', ' With mixed terms '},
21  ];
22 
23 path_arr2 = [
24  {'../Execute/3D_Without_0.1hr v1.28/Output/', ' No mixed terms'},
25  {'../Execute/3D_With_0.1hr v1.28/Output/', ' With mixed terms '},
26  ];
27 
28 
29 path1 = path_arr1{2,1};
30 path2 = path_arr1{1,1};
31 
32 name1 = path_arr1{2,2};
33 name2 = path_arr1{1,2};
34 
35 target_epc = [1];
36 target_alpha = [80];
37 
38 %c_min = -0.8;
39 %c_min = 1.2;
40 %c_max = 3.2;
41 c_min = 1;
42 c_max = 4;
43 n_colors = 100;
44 TextColor = 'white';
45 
46 
47 eval(['load ''', [path1, 'Flux_E=', num2str(target_epc), '_alpha', num2str(target_alpha)], ''''])
48 Flux1 = new_Flux;
49 L1 = new_L;
50 time1 = time_arr;
51 
52 eval(['load ''', [path2, 'Flux_E=', num2str(target_epc), '_alpha', num2str(target_alpha)], ''''])
53 Flux2 = new_Flux;
54 L2 = new_L;
55 time2 = time_arr;
56 
57 if (max(size(Flux1) ~= size(Flux2)))
58  Flux2 = interp2(time2, L2, Flux2', time1', L1, 'spline')';
59 end
60 
61 [ND, AND, MND] = NormDiff(Flux1, Flux2, 2);
62 
63 subplot1(4);
64 
65 Flux_min = 10^c_min;%Flux_max - 10^delta_color;
66 Flux1(find(Flux1<Flux_min))=Flux_min;
67 log_Flux = log10(Flux1);
68 
69 contourf(time_arr, L1, log_Flux', n_colors);
70 shading flat
71 axis tight
72 ylim([1.5 7]);
73 caxis([c_min c_max]);
74 colorbar
75 text(max(xlim)-0.05, (min(ylim) + max(ylim))/2, 'Flux, log_{10}(#/sr/s/cm^2/KeV)', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
76 %text(max(xlim), (min(ylim) + max(ylim))/2, '-3 \cdot 10^{-8} (c/MeV/cm)^3', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
77 %xlabel('Time, day of Apr. 2001');
78 xlabel('Time, days');
79 ylabel('L*');
80 text(min(xlim), max(ylim), [' (D) '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'FontWeight', 'bold', 'Color', 'white');
81 if (target_alpha >= 88)
82  text(max(xlim), min(ylim), [{name1}, {[' Energy = ', num2str(target_epc), ' MeV, equatorial pitch angle ']}], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'color', TextColor);
83 else
84  text(max(xlim), min(ylim), [{name1}, {[' Energy = ', num2str(target_epc), ' MeV, Pitch angle = ', num2str(target_alpha), '^o ']}], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'color', TextColor);
85 end
86 
87 
88 
89 subplot1(5);
90 
91 Flux_min = 10^c_min;%Flux_max - 10^delta_color;
92 Flux2(find(Flux2<Flux_min))=Flux_min;
93 log_Flux = log10(Flux2);
94 
95 contourf(time_arr, L1, log_Flux', n_colors);
96 shading flat
97 axis tight
98 ylim([1.5 7]);
99 caxis([c_min c_max]);
100 colorbar
101 text(max(xlim)-0.05, (min(ylim) + max(ylim))/2, 'Flux, log_{10}(#/sr/s/cm^2/KeV)', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
102 %text(max(xlim), (min(ylim) + max(ylim))/2, '-3 \cdot 10^{-8} (c/MeV/cm)^3', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
103 %xlabel('Time, day of Apr. 2001');
104 xlabel('Time, days');
105 %ylabel('L*');
106 text(min(xlim), max(ylim), [' (E) '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'FontWeight', 'bold', 'Color', 'white');
107 if (target_alpha >= 88)
108  text(max(xlim), min(ylim), [{name2}, {[' Energy = ', num2str(target_epc), ' MeV, equatorial pitch angle ']}], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'color', TextColor);
109 else
110  text(max(xlim), min(ylim), [{name2}, {[' Energy = ', num2str(target_epc), ' MeV, Pitch angle = ', num2str(target_alpha), '^o ']}], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'color', TextColor);
111 end
112 
113 
114 subplot1(6);
115 
116 contourf(time_arr, L1, ND', n_colors);
117 shading flat
118 axis tight
119 ylim([1.5 7]);
120 caxis([-75 75]);
121 colorbar
122 text(max(xlim), (min(ylim) + max(ylim))/2, 'Normalized difference, %', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
123 %text(max(xlim), (min(ylim) + max(ylim))/2, '-3 \cdot 10^{-8} (c/MeV/cm)^3', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
124 %xlabel('Time, day of Apr. 2001');
125 xlabel('Time, days');
126 %ylabel('L*');
127 text(max(xlim), min(ylim), {['(F_2 - F_1)/max_{at each t}(F_2/2 + F_1/2)']}, 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'color', 'black');
128 text(min(xlim), max(ylim), [' (F) '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'FontWeight', 'bold', 'Color', 'black');
129 
130 %%save_picture(gca, [path1, '../Plots/', 'ND_Flux_E=', num2str(target_epc), '_alpha', num2str(target_alpha)], 'jpg', 'pdf', 'epsc', 'fig');
131 
132 
133 
134 path1 = path_arr2{2,1};
135 path2 = path_arr2{1,1};
136 
137 name1 = path_arr2{2,2};
138 name2 = path_arr2{1,2};
139 
140 target_mu = [2000];
141 target_J = [0.01];
142 
143 c_min = 0;
144 c_max = 1.3;
145 TextColor = 'white';
146 
147 
148 eval(['load ''', path1, 'PSD_mu=', num2str(target_mu), '_J', num2str(target_J), ''' -mat'])
149 PSD1 = new_PSD;
150 L1 = new_L;
151 time1 = time_arr;
152 
153 eval(['load ''', path2, 'PSD_mu=', num2str(target_mu), '_J', num2str(target_J), ''' -mat'])
154 PSD2 = new_PSD;
155 L2 = new_L;
156 time2 = time_arr;
157 
158 if (max(size(PSD1) ~= size(PSD2)))
159  PSd2 = interp2(time2, L2, PSD2', time1', L1, 'spline')';
160 end
161 
162 [ND, AND, MND] = NormDiff(PSD1, PSD2, 2);
163 
164 subplot1(1);
165 
166 PSD_min = 10^c_min;
167 PSD1(find(PSD1<PSD_min))=PSD_min;
168 log_PSD = log10(PSD1);
169 
170 contourf(time_arr, L1, log_PSD', n_colors);
171 shading flat
172 axis tight
173 ylim([1.5 7]);
174 caxis([c_min c_max]);
175 colorbar
176 text(max(xlim)-0.05, (min(ylim) + max(ylim))/2, 'PSD, log_{10}(3 \cdot 10^{-8} (c/cm/MeV)^3)', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
177 %text(max(xlim), (min(ylim) + max(ylim))/2, '-3 \cdot 10^{-8} (c/MeV/cm)^3', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
178 %xlabel('Time, day of Apr. 2001');
179 xlabel('Time, days');
180 ylabel('L*');
181 text(max(xlim), min(ylim), [{name1}, {[' \mu = ', num2str(target_mu), ' MeV / G, J = ', num2str(target_J), 'MeV/c Re ']}], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'color', TextColor);
182 text(min(xlim), max(ylim), [' (A) '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'FontWeight', 'bold', 'Color', 'black');
183 
184 subplot1(2);
185 
186 PSD_min = 10^c_min;%PSD_max - 10^delta_color;
187 PSD2(find(PSD2<PSD_min))=PSD_min;
188 log_PSD = log10(PSD2);
189 
190 contourf(time_arr, L1, log_PSD', n_colors);
191 shading flat
192 axis tight
193 ylim([1.5 7]);
194 caxis([c_min c_max]);
195 colorbar
196 text(max(xlim)-0.05, (min(ylim) + max(ylim))/2, 'PSD, log_{10}(3 \cdot 10^{-8} (c/cm/MeV)^3)', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
197 %text(max(xlim), (min(ylim) + max(ylim))/2, '-3 \cdot 10^{-8} (c/MeV/cm)^3', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
198 %xlabel('Time, day of Apr. 2001');
199 xlabel('Time, days');
200 %ylabel('L*');
201 text(max(xlim), min(ylim), [{name2}, {[' \mu = ', num2str(target_mu), ' MeV / G, J = ', num2str(target_J), 'MeV/c Re ']}], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'color', TextColor);
202 text(min(xlim), max(ylim), [' (B) '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'FontWeight', 'bold', 'Color', 'black');
203 
204 subplot1(3);
205 
206 contourf(time_arr, L1, ND', n_colors);
207 shading flat
208 axis tight
209 ylim([1.5 7]);
210 caxis([-75 75]);
211 colorbar
212 text(max(xlim), (min(ylim) + max(ylim))/2, 'Normalized difference, %', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
213 %text(max(xlim), (min(ylim) + max(ylim))/2, '-3 \cdot 10^{-8} (c/MeV/cm)^3', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
214 %xlabel('Time, day of Apr. 2001');
215 xlabel('Time, days');
216 %ylabel('L*');
217 text(max(xlim), min(ylim), {['(F_2 - F_1)/max_{at each t}(F_2/2 + F_1/2)']}, 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'color', 'black');
218 text(min(xlim), max(ylim), [' (C) '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'FontWeight', 'bold', 'Color', 'black');
219 
220 %%save_picture(gca, [path1, '../Plots/', 'ND_PSD_mu=', num2str(target_mu), 'J=', num2str(target_J)], 'jpg', 'pdf', 'epsc', 'fig');
221 
222 save_picture(gcf, ['Figure_5'], 'jpg', 'pdf', 'epsc');
223 %save_picture(gcf, ['Figure_5_bw'], 'jpg', 'pdf', 'eps');
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
Phase Space Density class.
Definition: PSD.h:48
void spline(double *x, double *y, int n, double yp1, double ypn, double *y2)
Definition: spline.cpp:30