VERB_code_2.3
Plot_PSD_Energy_Alpha.m
1 %addpath('./Various_functions')
2 
3 %clear all;
4 %close all
5 
6 save_pictures = SAVEFIG;
7 save_movie = SAVEMOV;
8 
9 switch(FIGURE)
10  case 'a'
11  picture_name = 'Figure_12a';
12  case 'b'
13  picture_name = 'Figure_12a';
14  case 'c'
15  picture_name = 'Figure_12a';
16  case 'd'
17  picture_name = 'Figure_12a';
18 end
19 
20 %folder_name1=['../Execute/Output_', FIGURE, '/'];
21 
22 path_file_arr = [
23  %{'../Execute/previous_Output/'},
24  {['../Execute/Output_12', FIGURE, '/']}
25  %{'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\'}
26  ];
27 
28 % output parameters:
29 %output_formats = [{'fig'}, {'epsc'}, {'pdf'}];
30 output_formats = [{'png'}];
31 
32 n_colors = 50;
33 
34 c_min = 0;
35 c_max = 6;
36 
37 target_L = 4.5;
38 
39 
40 path_file = path_file_arr{1};
41 
42  gridFileName = [path_file, 'perp_grid.plt'];
43 % gridFileName = [path_file, 'rad_grid.plt'];
44 
45 oneDimDataFileName = [path_file, 'out1d.dat'];
46 
47 if (gridFileName(length(gridFileName)-3:length(gridFileName)) == '.dat')
48  [alpha, epc, L, pc] = load_plt(gridFileName, 'squeeze', 'permute');
49 else
50  [L, epc, alpha, pc] = load_plt(gridFileName, 'squeeze', 'permute');
51 end;
52 [time, Kp, Bf, Lpp] = load_plt(oneDimDataFileName);
53 
54 
55 for i = 1:length(path_file_arr)
56  path_file = path_file_arr{i};
57 
58  % compressions: Indeo5
59  if (save_movie)
60  mov = avifile([path_file, movie_name, '.avi'], 'fps', 10, 'quality', 0, 'compression', 'none');
61  end;
62 
63  f1 = figure('Units', 'inches', 'Position', [0.5 0.5 12 5], ...
64  'DefaultAxesFontName', 'times new roman', 'DefaultAxesFontWeight', 'bold', 'DefaultAxesFontSize', 12, 'DefaultTextFontWeight', 'bold', 'DefaultTextFontSize', 12,...
65  'PaperOrientation', 'portrait', 'PaperPositionMode', 'auto', 'DefaultLineLineWidth', 2);%, 'PaperPosition', [1 1 3.8 9]);%
66 
67  dataFileName = [path_file, 'OutPSD.dat'];
68  %dataFileName = [path_file, 'OutPSD_rad.dat'];
69 
70  PSD = load_plt(dataFileName, 'permute', 'nosqueeze', 'skip_zones', 1, 'n_zones', 999);%
71 
72  time_arr = PSD.time;
73 
74  [val, L_slice] = min(abs(L.arr(:,1,1) - target_L));
75 
76  if (length(size(alpha.arr)) == 3)
77  alpha_2d = squeeze(alpha.arr(L_slice,:,:));
78  energy_2d = squeeze(epc.arr(L_slice,:,:));
79  else
80  alpha_2d = alpha.arr;
81  energy_2d = epc.arr;
82  end
83 
84  set(gcf, 'Name', path_file);
85 
86  % comment/uncomment loop if you need just one pitcture
87  time_slice = length(PSD.time);
88 
89  if(strcmp(FIGURE,'a'))
90  time_slice_arr = 1;
91  else
92  time_slice_arr = length(PSD.time);
93  end
94 
95  for time_slice = time_slice_arr% 1:length(PSD.time)%
96 
97  subplot(1, 2, 1);
98  set(gca,'nextplot','replacechildren');
99 
100  PSD_2d = squeeze(PSD.arr(time_slice,L_slice,:,:));
101 
102 % PSD_min = 10^(log10(max(max(PSD_2d))) - 2);
103 % PSD_2d(PSD_2d<PSD_min)=PSD_min;
104 
105  log_abs_PSD = log10(abs(PSD_2d));
106  log_abs_PSD(isinf(log_abs_PSD)) = 0;
107  log_abs_PSD(find(log_abs_PSD>c_max))=c_max;
108  log_abs_PSD(find(log_abs_PSD<c_min))=c_min;
109 
110  contourf(alpha_2d, energy_2d, log_abs_PSD, n_colors);
111  set(gca, 'yscale', 'log');
112  shading flat
113 
114  caxis([c_min c_max]);
115  colorbar
116 
117  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');
118  text(min(xlim), max(ylim), ['Time=', num2str(PSD.time(time_slice)), ', days'], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'FontWeight', 'bold', 'Color', 'white');
119 
120  xlabel('Pitch angle, deg')
121  ylabel('Energy, MeV')
122 
123  subplot(1, 2, 2);
124  set(gca,'nextplot','replacechildren');
125 
126  %contourf(alpha_2d, energy_2d, sign(PSD_2d));
127  pcolor(alpha_2d, energy_2d, sign(PSD_2d));
128  shading flat
129  set(gca, 'yScale', 'log')
130  caxis([-1 1])
131  colorbar
132  xlabel('Pitch angle, deg');
133  ylabel('Energy, MeV');
134  text(max(xlim), 10^((min(log10(ylim)) + max(log10(ylim)))/2), 'sign(PSD)', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
135 
136  if (save_movie)
137  mov = addframe(mov, getframe(f1)); end;
138  getframe(f1);
139  %if (save_pictures) save_picture(gcf, [path_file, picture_name, '=', num2str(PSD.time(time_slice))], 'jpg', 'pdf', 'eps'); end;
140  if (save_pictures)
141  save_picture(gcf, ['../Figures/', picture_name, '=', num2str(PSD.time(time_slice))], 'png');
142  end; % ../Figures... add FIGURE to the picture_name
143 
144  end
145 
146  if (save_movie)
147  mov = close(mov);
148  end
149 
150 end
151 
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
Parameters_structure parameters
Parameters structure, with all parameters from the parameters.ini file. The default parameters define...
Definition: Main.cpp:185