VERB_code_2.3
Plot_PSD_Energy_Alpha.m
1 addpath('./Various_functions')
2 
3 clear all;
4 %close all
5 
6 save_movie = false;
7 save_pictures = false;
8 movie_name = '../Plots/movie1';
9 picture_name = '';
10 
11 
12 path_file_arr = [
13  %{'../Execute/previous_Output/'},
14  {'../Output/'},
15  %{'I:\YSSVNL_SVN\paper\Subbotin and Shprits 2009 Three-dimensional modeling of the radiation belts using the VERB code\Calculations\2009-06-27 3D paper recalculation\Fixed_parameters\Block and Split\Split\0.1\Output\'}
16  ];
17 
18 % output parameters:
19 output_formats = [{'fig'}, {'epsc'}, {'pdf'}];
20 n_colors = 50;
21 
22 c_min = -16;
23 c_max = -8;
24 
25 target_L = 4.2;
26 
27 
28 path_file = path_file_arr{1};
29 
30  gridFileName = [path_file, 'perp_grid.plt'];
31 % gridFileName = [path_file, 'rad_grid.plt'];
32 
33 oneDimDataFileName = [path_file, 'out1d.dat'];
34 
35 if (gridFileName(length(gridFileName)-3:length(gridFileName)) == '.dat')
36  [alpha, epc, L, pc] = load_plt(gridFileName, 'squeeze', 'permute');
37 else
38  [L, epc, alpha, pc] = load_plt(gridFileName, 'squeeze', 'permute');
39 end;
40 [time, Kp, Bf, Lpp] = load_plt(oneDimDataFileName);
41 
42 cl1 = 0; cl2 = 0;
43 for i = 1:length(path_file_arr)
44  path_file = path_file_arr{i}
45 
46  % compressions: Indeo5
47  if (save_movie) mov = avifile([path_file, movie_name, '.avi'], 'fps', 10, 'quality', 0, 'compression', 'none'); end;
48 
49  f1 = figure('Units', 'inches', 'Position', [0.5 0.5 12 5], ...
50  'DefaultAxesFontName', 'times new roman', 'DefaultAxesFontWeight', 'bold', 'DefaultAxesFontSize', 12, 'DefaultTextFontWeight', 'bold', 'DefaultTextFontSize', 12,...
51  'PaperOrientation', 'portrait', 'PaperPositionMode', 'auto', 'DefaultLineLineWidth', 2);%, 'PaperPosition', [1 1 3.8 9]);%
52 
53  dataFileName = [path_file, 'OutPSD.dat'];
54  %dataFileName = [path_file, 'OutPSD_rad.dat'];
55 
56  PSD = load_plt(dataFileName, 'permute', 'nosqueeze', 'skip_zones', 1, 'n_zones', 999);%
57  PSD.arr = PSD.arr/3e7;
58 
59  time_arr = PSD.time;
60 
61  [val, L_slice] = min(abs(L.arr(:,1,1) - target_L));
62 
63  if (length(size(alpha.arr)) == 3)
64  alpha_2d = squeeze(alpha.arr(L_slice,:,:))*180/pi;
65  energy_2d = squeeze(epc.arr(L_slice,:,:));
66  else
67  alpha_2d = alpha.arr*180/pi;
68  energy_2d = epc.arr;
69  end
70 
71  set(gcf, 'Name', path_file);
72 
73  % comment/uncomment loop if you need just one pitcture
74  time_slice = length(PSD.time);
75 
76 
77  for time_slice = 1:length(PSD.time)%
78 
79 
80  subplot(1, 2, 1);
81  set(gca,'nextplot','replacechildren');
82 
83  PSD_2d = squeeze(PSD.arr(time_slice,L_slice,:,:));
84 
85 % PSD_min = 10^(log10(max(max(PSD_2d))) - 2);
86 % PSD_2d(PSD_2d<PSD_min)=PSD_min;
87 
88  log_abs_PSD = log10(abs(PSD_2d));
89  log_abs_PSD(isinf(log_abs_PSD)) = 0;
90  log_abs_PSD(find(log_abs_PSD>c_max))=c_max;
91  log_abs_PSD(find(log_abs_PSD<c_min))=c_min;
92 
93  contourf(alpha_2d, energy_2d, log_abs_PSD, n_colors);
94  set(gca, 'yscale', 'log');
95  shading flat
96 
97  caxis([c_min c_max]);
98  if (cl1 == 0)
99  ax_pos = get(gca, 'Position');
100  cl1 = colorbar('West');
101  cl_pos = get(cl1, 'Position');
102  cl_x_pos = ax_pos(1) + ax_pos(3) + 0.025;
103  set(cl1, 'Position', [cl_x_pos, cl_pos(2), cl_pos(3), cl_pos(4)]);
104  end
105 
106 
107  xlim([0 90])
108  ylim([1e-1 1e1])
109 
110  text(max(xlim), 10^((min(log10(ylim)) + max(log10(ylim)))/2), 'log_{10}(|PSD|, 3 \cdot 10^{-8} c/MeV/cm^3)', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
111  text(min(xlim), max(ylim), ['Time=', num2str(PSD.time(time_slice)), ', days'], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'FontWeight', 'bold', 'Color', [0.5 0.5 0.5]);
112 
113  xlabel('Pitch angle, deg')
114  ylabel('Energy, MeV')
115 
116  subplot(1, 2, 2);
117  set(gca,'nextplot','replacechildren');
118 
119  %contourf(alpha_2d, energy_2d, sign(PSD_2d));
120  pcolor(alpha_2d, energy_2d, sign(PSD_2d));
121  shading flat
122  set(gca, 'yScale', 'log')
123  caxis([-1 1])
124 
125  if (cl2 == 0)
126  ax_pos = get(gca, 'Position');
127  cl2 = colorbar('West');
128  cl_pos = get(cl2, 'Position');
129  cl_x_pos = ax_pos(1) + ax_pos(3) + 0.025;%0.015; %0.925
130  set(cl2, 'Position', [cl_x_pos, cl_pos(2), cl_pos(3), cl_pos(4)]);
131 
132  end
133 
134  xlabel('Pitch angle, deg');
135  ylabel('Energy, MeV');
136  text(max(xlim), 10^((min(log10(ylim)) + max(log10(ylim)))/2), 'sign(PSD)', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
137 
138  xlim([0 90])
139  ylim([1e-1 1e1])
140 
141  if (save_movie) mov = addframe(mov, getframe(f1)); end;
142  getframe(f1);
143  if (save_pictures) save_picture(gcf, [path_file, picture_name, '=', num2str(PSD.time(time_slice))], 'jpg', 'pdf', 'eps'); end;
144  end
145 
146  if (save_movie) mov = close(mov); end;
147 
148 end
149 
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