1 addpath(
'./Various_functions')
8 movie_name = '../Plots/movie1';
13 %{
'../Execute/previous_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\'}
19 output_formats = [{'fig
'}, {'epsc
'}, {'pdf
'}];
28 path_file = path_file_arr{1};
30 gridFileName = [path_file, 'perp_grid.plt
'];
31 % gridFileName = [path_file, 'rad_grid.plt
'];
33 oneDimDataFileName = [path_file, 'out1d.dat
'];
35 if (gridFileName(length(gridFileName)-3:length(gridFileName)) == '.dat
')
36 [alpha, epc, L, pc] = load_plt(gridFileName, 'squeeze
', 'permute
');
38 [L, epc, alpha, pc] = load_plt(gridFileName, 'squeeze
', 'permute
');
40 [time, Kp, Bf, Lpp] = load_plt(oneDimDataFileName);
43 for i = 1:length(path_file_arr)
44 path_file = path_file_arr{i}
46 % compressions: Indeo5
47 if (save_movie) mov = avifile([path_file, movie_name, '.avi
'], 'fps
', 10, 'quality
', 0, 'compression
', 'none
'); end;
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]);%
53 dataFileName = [path_file, 'OutPSD.dat
'];
54 %dataFileName = [path_file, 'OutPSD_rad.dat
'];
56 PSD = load_plt(dataFileName, 'permute
', 'nosqueeze
', 'skip_zones
', 1, 'n_zones
', 999);%
57 PSD.arr = PSD.arr/3e7;
61 [val, L_slice] = min(abs(L.arr(:,1,1) - target_L));
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,:,:));
67 alpha_2d = alpha.arr*180/pi;
71 set(gcf, 'Name
', path_file);
73 % comment/uncomment loop if you need just one pitcture
74 time_slice = length(PSD.time);
77 for time_slice = 1:length(PSD.time)%
81 set(gca,'nextplot
','replacechildren
');
83 PSD_2d = squeeze(PSD.arr(time_slice,L_slice,:,:));
85 % PSD_min = 10^(log10(max(max(PSD_2d))) - 2);
86 % PSD_2d(PSD_2d<PSD_min)=PSD_min;
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;
93 contourf(alpha_2d, energy_2d, log_abs_PSD, n_colors);
94 set(gca, 'yscale
', 'log
');
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)]);
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]);
113 xlabel(
'Pitch angle, deg')
114 ylabel('Energy, MeV')
117 set(gca,'nextplot','replacechildren');
119 %contourf(alpha_2d, energy_2d, sign(PSD_2d));
120 pcolor(alpha_2d, energy_2d, sign(PSD_2d));
122 set(gca, 'yScale', 'log')
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)]);
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');
141 if (save_movie) mov = addframe(mov, getframe(f1)); end;
143 if (save_pictures) save_picture(
gcf, [path_file, picture_name, '=', num2str(PSD.time(time_slice))], 'jpg', 'pdf', 'eps'); end;
146 if (save_movie) mov = close(mov); end;
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...
Phase Space Density class.