VERB_code_2.3
Interp_Dxx_2d_2d.m
1 clear all
2 %close all
3 
4 do_corr = true
5 do_smooth = false
6 
7 % for hiss
8 Kp = 4
9 
10 % Hiss scale plasmasphere
11 BwInFile = 100;
12 BwForDxxKp = 40;
13 DxxKp = 4/(BwForDxxKp/BwInFile);
14 %HissScale = (BwForDxxKp/BwInFile)^2 * (Kp/DxxKp)^2
15 HissScale_pp = (Kp/DxxKp)^2
16 
17 %Chorus scale
18 BwInFile = 50;
19 BwForDxxKp = 50;
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
25 
26 %%% name, MLT, Scale:
27 Dxx_file_arr = {
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},
33  };
34 
35 output_arr = [
36  {'./DiffCoeff/'},
37  ];
38 
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/';
44 
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}
49  try
50  postfix = ['_', Dxx_file_arr{i}{4}];
51  catch
52  postfix = '';
53  end
54 
55  tic
56 
57  'loading'
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;
61 
62  L.comments = '';
63  epc.comments = '';
64  alpha.comments = '';
65 
66  % for VLF
67  Daa.comments = [Daa.comments, ['# MLT = ', num2str(MLT)]];
68  Dpp.comments = [Dpp.comments, ['# MLT = ', num2str(MLT)]];
69  Dpa.comments = [Dpa.comments, ['# MLT = ', num2str(MLT)]];
70  if (Scale~=1)
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)]];
74  end
75 
76  for output_idx = 1:size(output_arr,1)
77  output = output_arr{output_idx};
78 
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);
84 
85 
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))];
91  L_slice = 1;
92 
93 % new_L.arr = L.arr(epc_slice, :);
94 % new_epc.arr = epc.arr(epc_slice, :);
95 % new_alpha.arr = alpha.arr(epc_slice, :);
96 
97  Daa_new = Daa;
98  Daa_new.name = [Daa_new.name, postfix];
99  Daa_new.name
100  Daa_new.fname = [strrep(Daa_new.name, ' ', '_'), '_Full_LPA'];
101  Daa_new.comments = [Daa_new.comments, ['# File: ', pwd, Dxx_file]];
102 
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);
107 
108  Dpp_new = Dpp;
109  Dpp_new.name = [Dpp_new.name, postfix];
110  Dpp_new.name
111  Dpp_new.fname = [strrep(Dpp_new.name, ' ', '_'), '_Full_LPA'];
112  Dpp_new.comments = [Dpp_new.comments, ['# File: ', pwd, Dxx_file]];
113 
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);
118 
119  Dpa_new = Dpa;
120  Dpa_new.name = [Dpa_new.name, postfix];
121  Dpa_new.name
122  Dpa_new.fname = [strrep(Dpa_new.name, ' ', '_'), '_Full_LPA'];
123  Dpa_new.comments = [Dpa_new.comments, ['# File: ', pwd, Dxx_file]];
124 
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;
128 
129  % correct Dpa where too big
130  if (do_corr)
131  total_correction = 0;
132  tmp = Dpa_new.arr;
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']];
137  end;
138 
139  save_plt_lpa([output, Dpa_new.fname, '.plt'], new_L, new_epc, new_alpha, Dpa_new);
140 
141  toc
142 
143  new_pc.arr = pfunc(new_epc.arr);
144  new_L.name = 'L';
145  new_epc.name = 'Energy, MeV';
146  new_alpha.name = 'Pitch angle, deg';
147  new_pc.name = 'pc';
148  save_plt_lpa([output, 'grid1.plt'], new_L, new_epc, new_alpha, new_pc);
149 
150  end
151 
152 
153  %return
154  % XIAO compare
155 
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]);%
159 
160  subplot1(2,2, 'Max', [0.90 0.90], 'Min', [0.12 0.10]);
161 
162  units = 1
163  c_min = -6
164  c_max = 2
165  %units = 1/60/60/24
166  %c_min = -12
167  %c_max = -4
168 
169  subplot1(1);
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);
174  shading flat
175  caxis([c_min c_max])
176  axis tight
177  set(gca, 'yScale', 'log')
178  ylim([1e-1 1e1])
179 
180  subplot1(2);
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);
185  shading flat
186  set(gca, 'yScale', 'log')
187  caxis([c_min c_max])
188  axis tight
189  ylim([1e-1 1e1])
190 
191  subplot1(3);
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);
196  shading flat
197  set(gca, 'yScale', 'log')
198  caxis([-1 1])
199  axis tight
200  ylim([1e-1 1e1])
201 
202  subplot1(4);
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);
207  shading flat
208  set(gca, 'yScale', 'log')
209  caxis([c_min c_max])
210  axis tight
211  ylim([1e-1 1e1])
212 
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)]);
218 
219 end
220 
221 
222 return
223 
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]);%
227 
228 subplot1(1,3, 'Max', [0.90 0.90], 'Min', [0.12 0.10]);
229 
230 units = 1;%/60/60/24
231 c_min = -6
232 c_max = 2
233 
234 subplot1(1);
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);
239 shading flat
240 caxis([c_min c_max])
241 axis tight
242 set(gca, 'yScale', 'log')
243 ylim([2e-1 1e1])
244 
245 subplot1(2);
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);
250 shading flat
251 set(gca, 'yScale', 'log')
252 caxis([c_min c_max])
253 axis tight
254 ylim([2e-1 1e1])
255 
256 subplot1(3);
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);
261 shading flat
262 set(gca, 'yScale', 'log')
263 caxis([c_min c_max])
264 axis tight
265 ylim([2e-1 1e1])
266 
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)]);
272 
273 return
274 
275 gridFileName = [output, 'perp_grid.plt'];
276 [new_L2, new_epc2, new_alpha2, new_pc2] = load_plt(gridFileName, 'squeeze', 'permute');
277 
278 
279 figure
280 hold on
281 plot(squeeze(new_epc.arr(1,:,91)))
282 plot(squeeze(new_epc2.arr(:,91)))
283 
284 
285 figure
286 hold on
287 plot(squeeze(new_alpha.arr(1,1,:)), 'o')
288 plot(squeeze(new_alpha2.arr(1,:)), 'o')
289 
290 
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]);%
294 
295 min_c = -6;
296 max_c = 2;
297 
298 
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]);
304 colorbar
305 
306 return
307 
308 % day side
309 
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))
313 
314 new_epc_plot = [new_epc_0_5, new_epc_1_0, new_epc_2_0]
315 
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]);%
319 
320 subplot1(3,1);
321 
322 subplot1(1);
323 plot(squeeze(new_alpha.arr(1,:))*180/pi, squeeze(Daa_new.arr(new_epc_plot,:)/60/60/24));
324 set(gca, 'yScale', 'log')
325 ylim([5e-7 1e-4])
326 xlim([0 90])
327 text(max(xlim), max(ylim), [num2str(new_epc.arr(new_epc_plot)), ' MeV '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'FontSize', 14, 'Color', 'black');
328 
329 
330 subplot1(2);
331 plot(squeeze(new_alpha.arr(1,:))*180/pi, squeeze(Dpp_new.arr(new_epc_plot,:)/60/60/24));
332 set(gca, 'yScale', 'log')
333 ylim([1e-8 3e-6])
334 xlim([0 90])
335 
336 subplot1(3);
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')
339 ylim([1e-9 1e-5])
340 xlim([0 90])
341 
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]
346 
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]);%
350 
351 subplot1(3,1);
352 
353 subplot1(1);
354 plot(squeeze(alpha.arr(1,:))*180/pi, squeeze(Daa.arr(epc_plot,:)/60/60/24));
355 set(gca, 'yScale', 'log')
356 ylim([5e-7 1e-4])
357 xlim([0 90])
358 text(max(xlim), max(ylim), [num2str(epc.arr(epc_plot)), ' MeV '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'FontSize', 14, 'Color', 'black');
359 
360 
361 subplot1(2);
362 plot(squeeze(alpha.arr(1,:))*180/pi, squeeze(Dpp.arr(epc_plot,:)/60/60/24));
363 set(gca, 'yScale', 'log')
364 ylim([1e-8 3e-6])
365 xlim([0 90])
366 
367 subplot1(3);
368 plot(squeeze(alpha.arr(1,:))*180/pi, squeeze(abs(Dpa.arr(epc_plot,:))/60/60/24));
369 set(gca, 'yScale', 'log')
370 ylim([1e-9 1e-5])
371 xlim([0 90])
372 
373 
374 % night side
375 
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))
379 
380 new_epc_plot = [new_epc_0_5, new_epc_1_0, new_epc_2_0]
381 
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]);%
385 
386 subplot1(3,1);
387 
388 subplot1(1);
389 plot(squeeze(new_alpha.arr(1,:))*180/pi, squeeze(Daa_new.arr(new_epc_plot,:)/60/60/24));
390 set(gca, 'yScale', 'log')
391 ylim([1e-7 2e-4])
392 xlim([0 90])
393 text(max(xlim), max(ylim), [num2str(new_epc.arr(new_epc_plot)), ' MeV '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'FontSize', 14, 'Color', 'black');
394 
395 
396 subplot1(2);
397 plot(squeeze(new_alpha.arr(1,:))*180/pi, squeeze(Dpp_new.arr(new_epc_plot,:)/60/60/24));
398 set(gca, 'yScale', 'log')
399 ylim([1e-9 1e-5])
400 xlim([0 90])
401 
402 subplot1(3);
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')
405 ylim([1e-9 4e-5])
406 xlim([0 90])
407 
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]
412 
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]);%
416 
417 subplot1(3,1);
418 
419 subplot1(1);
420 plot(squeeze(alpha.arr(1,:))*180/pi, squeeze(Daa.arr(epc_plot,:)/60/60/24));
421 set(gca, 'yScale', 'log')
422 ylim([1e-7 2e-4])
423 xlim([0 90])
424 text(max(xlim), max(ylim), [num2str(epc.arr(epc_plot)), ' MeV '], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'FontSize', 14, 'Color', 'black');
425 
426 
427 subplot1(2);
428 plot(squeeze(alpha.arr(1,:))*180/pi, squeeze(Dpp.arr(epc_plot,:)/60/60/24));
429 set(gca, 'yScale', 'log')
430 ylim([1e-9 1e-5])
431 xlim([0 90])
432 
433 subplot1(3);
434 plot(squeeze(alpha.arr(1,:))*180/pi, squeeze(abs(Dpa.arr(epc_plot,:))/60/60/24));
435 set(gca, 'yScale', 'log')
436 ylim([1e-9 4e-5])
437 xlim([0 90])
438 
439 return
440 
441 gridFileName = [output, 'perp_grid.plt'];
442 [my_L, my_epc, my_alpha, my_pc] = load_plt(gridFileName, 'squeeze', 'permute');
443 
444 
445 figure
446 plot(new_epc.arr(:,15));
447 hold on
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 .