10 % Hiss scale plasmasphere
13 DxxKp = 4/(BwForDxxKp/BwInFile);
14 %HissScale = (BwForDxxKp/BwInFile)^2 * (Kp/DxxKp)^2
15 HissScale_pp = (Kp/DxxKp)^2
20 DxxKp = 4/(BwForDxxKp/BwInFile);
21 Bw2_0 = 2.0*10^(2.5 + 0.18 * DxxKp);
22 if (Kp <= 2.333) Bw2 = 2.0*10^(0.73 + 0.91 * Kp); end
23 if (Kp > 2.333) Bw2 = 2.0*10^(2.5 + 0.18 * Kp); end
24 ChorusScale = Bw2/Bw2_0
28 % {
'C:/Math/Projects/Calculations/2009-01-09 Full Dxx 2/2009-03-23 Outer 2d L=4.5/FullCode_DaysideChorus_BavD_Matrix_L45.txt', 0.25, ChorusScale},
29 % {
'C:/Math/Projects/Calculations/2009-01-09 Full Dxx 2/2009-03-23 Outer 2d L=4.5/FullCode_NightsideChorus_BavD_Matrix_L45.txt', 0.25, ChorusScale},
30 % {
'./2009-03-23 Outer 2d L=4.5/FullCode_PlasmasphericHiss_SmallWaveNormalAngle_BavD_Matrix_L45.txt', 0.6, HissScale_pp}
31 {
'C:/Math/Projects/Calculations/2009-01-09 Full Dxx 2/2009-06-02 diffusion rate Li et al 2007/DaysideChorus_OuterBelt_BavD_Matrix_L45_Lietal2007.txt', 0.25, ChorusScale},
32 {
'C:/Math/Projects/Calculations/2009-01-09 Full Dxx 2/2009-06-02 diffusion rate Li et al 2007/NightsideChorus_OuterBelt_BavD_Matrix_L45_Lietal2007.txt', 0.25, ChorusScale},
39 %output =
'../Interpolated/2009-03-19 1d from 3d en=1.0965/';
40 %output =
'../Interpolated/2009-03-19 1d from 3d en=2.63026/';
41 %output =
'../Interpolated/2009-03-19 1d from 3d en=2.7542/';
42 %output =
'../Interpolated/2009-03-19 1d from 3d en=4.7863/';
43 %output =
'../Interpolated/2009-03-19 1d from 3d en=5.01187/';
45 for i = 1:size(Dxx_file_arr,1)
46 Dxx_file = Dxx_file_arr{i}{1}
47 MLT = Dxx_file_arr{i}{2}
48 Scale = Dxx_file_arr{i}{3}
50 postfix = ['_', Dxx_file_arr{i}{4}];
58 [L,epc,alpha,Daa,Dpp,Dpa] = load_plt(Dxx_file,
'permute',
'squeeze');
59 epc.arr = epc.arr / 10^3;
60 alpha.arr = alpha.arr / 180*pi;
67 Daa.comments = [Daa.comments, ['# MLT = ', num2str(MLT)]];
68 Dpp.comments = [Dpp.comments, ['# MLT = ', num2str(MLT)]];
69 Dpa.comments = [Dpa.comments, ['# MLT = ', num2str(MLT)]];
71 Daa.comments = [Daa.comments, ['# !!! Multiplied by ', num2str(Scale)]];
72 Dpp.comments = [Dpp.comments, ['# !!! Multiplied by ', num2str(Scale)]];
73 Dpa.comments = [Dpa.comments, ['# !!! Multiplied by ', num2str(Scale)]];
76 for output_idx = 1:size(output_arr,1)
77 output = output_arr{output_idx};
79 gridFileName = [output, '../Output/perp_grid.plt'];
80 [new_L, new_epc, new_alpha, new_pc] = load_plt(gridFileName,
'squeeze',
'permute');
81 new_L_size = size(new_L, 1);
82 new_epc_size = size(new_epc, 2);
83 new_alpha_size = size(new_alpha, 3);
86 % [dif, epc_slice] = min(abs(epc.arr(1,:,1) - unique(new_epc.arr)));
87 %epc_slice = [min(find(epc.arr(:,1) >= 0.01)) : 2 : max(find(epc.arr(:,1) <= 10.0))];
88 % epc_slice = [min(find(epc.arr(:,1) >= 0.02)) : 1 : max(find(epc.arr(:,1) <= 10.0))];
89 % epc_slice = [min(find(epc.arr(:,1) >= 0.1)) : 1 : max(find(epc.arr(:,1) <= 10.0))];
90 % epc_slice = [min(find(epc.arr(:,1) >= 0.2)) : 1 : max(find(epc.arr(:,1) <= 10.0))];
93 % new_L.arr = L.arr(epc_slice, :);
94 % new_epc.arr = epc.arr(epc_slice, :);
95 % new_alpha.arr = alpha.arr(epc_slice, :);
98 Daa_new.name = [Daa_new.name, postfix];
100 Daa_new.fname = [strrep(Daa_new.name, ' ', '_'), '_Full_LPA'];
101 Daa_new.comments = [Daa_new.comments, ['# File: ', pwd, Dxx_file]];
103 % Daa_new.arr = Daa.arr(epc_slice, :)*MLT*Scale;
104 Daa_new.arr = griddata(epc.arr, alpha.arr, Daa.arr, new_epc.arr, new_alpha.arr)*MLT*Scale;
105 if (do_smooth) Daa_new.arr = smooth2(Daa_new.arr, 1, 10); end;
106 save_plt_lpa([output, Daa_new.fname,
'.plt'], new_L, new_epc, new_alpha, Daa_new);
109 Dpp_new.name = [Dpp_new.name, postfix];
111 Dpp_new.fname = [strrep(Dpp_new.name, ' ', '_'), '_Full_LPA'];
112 Dpp_new.comments = [Dpp_new.comments, ['# File: ', pwd, Dxx_file]];
114 % Dpp_new.arr = Dpp.arr(epc_slice, :)*MLT*Scale;
115 Dpp_new.arr = griddata(epc.arr, alpha.arr, Dpp.arr, new_epc.arr, new_alpha.arr)*MLT*Scale;
116 if (do_smooth) Dpp_new.arr = smooth2(Dpp_new.arr, 1, 10); end;
117 save_plt_lpa([output, Dpp_new.fname,
'.plt'], new_L, new_epc, new_alpha, Dpp_new);
120 Dpa_new.name = [Dpa_new.name, postfix];
122 Dpa_new.fname = [strrep(Dpa_new.name, ' ', '_'), '_Full_LPA'];
123 Dpa_new.comments = [Dpa_new.comments, ['# File: ', pwd, Dxx_file]];
125 % Dpa_new.arr = Dpa.arr(epc_slice, :)*MLT*Scale;
126 Dpa_new.arr = griddata(epc.arr, alpha.arr, Dpa.arr, new_epc.arr, new_alpha.arr)*MLT*Scale;
127 if (do_smooth) Dpa_new.arr = smooth2(Dpa_new.arr, 1, 10); end;
129 % correct Dpa where too big
131 total_correction = 0;
133 diff = find((Dpa_new.arr.*Dpa_new.arr) ./ (Daa_new.arr.*Dpp_new.arr) >= 1);
134 Dpa_new.arr(diff) = sign(Dpa_new.arr(diff)).*(sqrt(Daa_new.arr(diff).*Dpp_new.arr(diff))*0.99);
135 total_correction = sum(sum(abs(tmp - Dpa_new.arr)))/(size(Dpa_new.arr,1)*size(Dpa_new.arr,2))
136 Dpa_new.comments = [Dpa_new.comments, [
'# Corrected to satisfy Dpp*Daa>Dpa^2 (in average): ', num2str(total_correction),
' per point']];
139 save_plt_lpa([output, Dpa_new.fname,
'.plt'], new_L, new_epc, new_alpha, Dpa_new);
143 new_pc.arr =
pfunc(new_epc.arr);
145 new_epc.name =
'Energy, MeV';
146 new_alpha.name =
'Pitch angle, deg';
148 save_plt_lpa([output,
'grid1.plt'], new_L, new_epc, new_alpha, new_pc);
156 f1 = figure(
'Units',
'inches',
'Position', [0.5 0.5 9 8], ...
157 'DefaultAxesFontName',
'times new roman',
'DefaultAxesFontWeight',
'bold',
'DefaultAxesFontSize', 12,...
158 'PaperOrientation',
'portrait',
'PaperPositionMode',
'auto',
'DefaultLineLineWidth', 2);%,
'PaperPosition', [1 1 3.8 9]);%
160 subplot1(2,2,
'Max', [0.90 0.90],
'Min', [0.12 0.10]);
170 Dxx = log10(Daa_new.arr*units);
171 Dxx(isnan(Dxx)) = c_min;
172 Dxx(find(Dxx<c_min)) = c_min;
173 contourf(new_alpha.arr*180/pi, new_epc.arr, Dxx);
177 set(gca,
'yScale',
'log')
181 Dxx = log10(Dpp_new.arr*units);
182 Dxx(isnan(Dxx)) = c_min;
183 Dxx(find(Dxx<c_min)) = c_min;
184 contourf(new_alpha.arr*180/pi, new_epc.arr, Dxx);
186 set(gca, 'yScale', 'log')
192 Dxx = sign(Dpa_new.arr);
193 Dxx(isnan(Dxx)) = c_min;
194 Dxx(find(Dxx<c_min)) = c_min;
195 contourf(new_alpha.arr*180/pi, new_epc.arr, Dxx);
197 set(gca, 'yScale', 'log')
203 Dxx = log10(abs(Dpa_new.arr)*units);
204 Dxx(isnan(Dxx)) = c_min;
205 Dxx(find(Dxx<c_min)) = c_min;
206 contourf(new_alpha.arr*180/pi, new_epc.arr, Dxx);
208 set(gca, 'yScale', 'log')
213 ax_pos = get(gca, 'Position');
214 cl = colorbar('West');
215 cl_pos = get(cl, 'Position');
216 cl_x_pos = ax_pos(1) + ax_pos(3) + 0.035;%0.015; %0.925
217 set(cl, 'Position', [cl_x_pos, cl_pos(2), cl_pos(3), cl_pos(4)]);
224 f1 = figure('Units', 'inches', 'Position', [0.5 0.5 12 5], ...
225 'DefaultAxesFontName', 'times new roman', 'DefaultAxesFontWeight', 'bold', 'DefaultAxesFontSize', 12,...
226 'PaperOrientation', 'portrait', 'PaperPositionMode', 'auto', 'DefaultLineLineWidth', 2);%, 'PaperPosition', [1 1 3.8 9]);%
228 subplot1(1,3, 'Max', [0.90 0.90], 'Min', [0.12 0.10]);
235 Dxx = log10(Daa_new.arr*units);
236 Dxx(isnan(Dxx)) = c_min;
237 Dxx(find(Dxx<c_min)) = c_min;
238 contourf(new_alpha.arr*180/pi, new_epc.arr, Dxx);
242 set(gca, 'yScale', 'log')
246 Dxx = log10(Dpp_new.arr*units);
247 Dxx(isnan(Dxx)) = c_min;
248 Dxx(find(Dxx<c_min)) = c_min;
249 contourf(new_alpha.arr*180/pi, new_epc.arr, Dxx);
251 set(gca, 'yScale', 'log')
257 Dxx = log10(Dpa_new.arr*units);
258 Dxx(isnan(Dxx)) = c_min;
259 Dxx(find(Dxx<c_min)) = c_min;
260 contourf(new_alpha.arr*180/pi, new_epc.arr, Dxx);
262 set(gca, 'yScale', 'log')
267 ax_pos = get(gca, 'Position');
268 cl = colorbar('West');
269 cl_pos = get(cl, 'Position');
270 cl_x_pos = ax_pos(1) + ax_pos(3) + 0.035;%0.015; %0.925
271 set(cl, 'Position', [cl_x_pos, cl_pos(2), cl_pos(3), cl_pos(4)]);
275 gridFileName = [output, 'perp_grid.plt'];
276 [new_L2, new_epc2, new_alpha2, new_pc2] = load_plt(gridFileName, 'squeeze', 'permute');
281 plot(squeeze(new_epc.arr(1,:,91)))
282 plot(squeeze(new_epc2.arr(:,91)))
287 plot(squeeze(new_alpha.arr(1,1,:)), 'o')
288 plot(squeeze(new_alpha2.arr(1,:)), 'o')
291 f1 = figure('Units', 'inches','Position', [1. 1. 11 6], ...
292 'DefaultAxesFontName', 'times new roman', 'DefaultAxesFontWeight', 'bold', 'DefaultAxesFontSize', 12,...
293 'PaperOrientation', 'landscape', 'PaperPositionMode', 'auto', 'DefaultLineLineWidth', 2);%, 'PaperPosition', [1 1 3.8 9]);%
299 Dplot(:,:) = log10(squeeze(Daa_new.arr(1, :,:)));
300 Dplot(find(~isfinite(Dplot)))=min_c; Dplot(find(Dplot<min_c))=min_c; Dplot(find(Dplot>max_c))=max_c;
301 contourf(squeeze(new_alpha.arr(1, :,:).*180./pi), squeeze(new_epc.arr(1, :,:)), Dplot);
302 xlim([0 90]); %ylim([emin emax]);
303 shading flat; set(gca,'Yscale','log'); caxis([min_c max_c]);
310 new_epc_0_5 =
max(find(new_epc.arr(:,1)<0.5))
311 new_epc_1_0 =
max(find(new_epc.arr(:,1)<1))
312 new_epc_2_0 =
max(find(new_epc.arr(:,1)<2))
314 new_epc_plot = [new_epc_0_5, new_epc_1_0, new_epc_2_0]
316 f1 = figure('Units', 'inches', 'Position', [0.5 0.5 5 9], ...
317 'DefaultAxesFontName', 'times new roman', 'DefaultAxesFontWeight', 'bold', 'DefaultAxesFontSize', 12,...
318 'PaperOrientation', 'portrait', 'PaperPositionMode', 'auto', 'DefaultLineLineWidth', 2);%, 'PaperPosition', [1 1 3.8 9]);%
323 plot(squeeze(new_alpha.arr(1,:))*180/pi, squeeze(Daa_new.arr(new_epc_plot,:)/60/60/24));
324 set(gca, 'yScale', 'log')
327 text(
max(xlim),
max(ylim), [num2str(new_epc.arr(new_epc_plot)), ' MeV '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'FontSize', 14, 'Color', 'black');
331 plot(squeeze(new_alpha.arr(1,:))*180/pi, squeeze(Dpp_new.arr(new_epc_plot,:)/60/60/24));
332 set(gca, 'yScale', 'log')
337 plot(squeeze(new_alpha.arr(1,:))*180/pi, squeeze(abs(Dpa_new.arr(new_epc_plot,:))/60/60/24));
338 set(gca, 'yScale', 'log')
342 epc_0_5 =
max(find(epc.arr(:,1)<0.5))
343 epc_1_0 =
max(find(epc.arr(:,1)<1))
344 epc_2_0 =
max(find(epc.arr(:,1)<2))
345 epc_plot = [epc_0_5, epc_1_0, epc_2_0]
347 f1 = figure('Units', 'inches', 'Position', [0.5 0.5 5 9], ...
348 'DefaultAxesFontName', 'times new roman', 'DefaultAxesFontWeight', 'bold', 'DefaultAxesFontSize', 12,...
349 'PaperOrientation', 'portrait', 'PaperPositionMode', 'auto', 'DefaultLineLineWidth', 2);%, 'PaperPosition', [1 1 3.8 9]);%
354 plot(squeeze(alpha.arr(1,:))*180/pi, squeeze(Daa.arr(epc_plot,:)/60/60/24));
355 set(gca, 'yScale', 'log')
358 text(
max(xlim),
max(ylim), [num2str(epc.arr(epc_plot)), ' MeV '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'FontSize', 14, 'Color', 'black');
362 plot(squeeze(alpha.arr(1,:))*180/pi, squeeze(Dpp.arr(epc_plot,:)/60/60/24));
363 set(gca, 'yScale', 'log')
368 plot(squeeze(alpha.arr(1,:))*180/pi, squeeze(abs(Dpa.arr(epc_plot,:))/60/60/24));
369 set(gca, 'yScale', 'log')
376 new_epc_0_5 =
max(find(new_epc.arr(:,1)<0.5))
377 new_epc_1_0 =
max(find(new_epc.arr(:,1)<1))
378 new_epc_2_0 =
max(find(new_epc.arr(:,1)<2))
380 new_epc_plot = [new_epc_0_5, new_epc_1_0, new_epc_2_0]
382 f1 = figure('Units', 'inches', 'Position', [0.5 0.5 5 9], ...
383 'DefaultAxesFontName', 'times new roman', 'DefaultAxesFontWeight', 'bold', 'DefaultAxesFontSize', 12,...
384 'PaperOrientation', 'portrait', 'PaperPositionMode', 'auto', 'DefaultLineLineWidth', 2);%, 'PaperPosition', [1 1 3.8 9]);%
389 plot(squeeze(new_alpha.arr(1,:))*180/pi, squeeze(Daa_new.arr(new_epc_plot,:)/60/60/24));
390 set(gca, 'yScale', 'log')
393 text(
max(xlim),
max(ylim), [num2str(new_epc.arr(new_epc_plot)), ' MeV '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'FontSize', 14, 'Color', 'black');
397 plot(squeeze(new_alpha.arr(1,:))*180/pi, squeeze(Dpp_new.arr(new_epc_plot,:)/60/60/24));
398 set(gca, 'yScale', 'log')
403 plot(squeeze(new_alpha.arr(1,:))*180/pi, squeeze(abs(Dpa_new.arr(new_epc_plot,:))/60/60/24));
404 set(gca, 'yScale', 'log')
408 epc_0_5 =
max(find(epc.arr(:,1)<0.5))
409 epc_1_0 =
max(find(epc.arr(:,1)<1))
410 epc_2_0 =
max(find(epc.arr(:,1)<2))
411 epc_plot = [epc_0_5, epc_1_0, epc_2_0]
413 f1 = figure('Units', 'inches', 'Position', [0.5 0.5 5 9], ...
414 'DefaultAxesFontName', 'times new roman', 'DefaultAxesFontWeight', 'bold', 'DefaultAxesFontSize', 12,...
415 'PaperOrientation', 'portrait', 'PaperPositionMode', 'auto', 'DefaultLineLineWidth', 2);%, 'PaperPosition', [1 1 3.8 9]);%
420 plot(squeeze(alpha.arr(1,:))*180/pi, squeeze(Daa.arr(epc_plot,:)/60/60/24));
421 set(gca, 'yScale', 'log')
424 text(
max(xlim),
max(ylim), [num2str(epc.arr(epc_plot)), ' MeV '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'FontSize', 14, 'Color', 'black');
428 plot(squeeze(alpha.arr(1,:))*180/pi, squeeze(Dpp.arr(epc_plot,:)/60/60/24));
429 set(gca, 'yScale', 'log')
434 plot(squeeze(alpha.arr(1,:))*180/pi, squeeze(abs(Dpa.arr(epc_plot,:))/60/60/24));
435 set(gca, 'yScale', 'log')
441 gridFileName = [output, 'perp_grid.plt'];
442 [my_L, my_epc, my_alpha, my_pc] = load_plt(gridFileName, 'squeeze', 'permute');
446 plot(new_epc.arr(:,15));
448 plot(my_epc.arr(:,15), 'r');
449 set(gca, 'yscale', 'log');
double max(double v1, double v2)
Return maximum.
double pfunc(double K)
Computation of moumentum from Kinetic energy .