4 addpath(
'./Various_functions')
6 DOY_start=293; % day of the year of the first day of simulation
10 {
'../Execute/Output/'},
18 output_formats = [{'fig'}, {'epsc'}, {'pdf'}];
23 for pic = 1:length(path_arr)
24 path1 = path_arr{pic};
40 units_constant = 1;%3*10^-2; % units 10^-6*(c/MeV/cm)^3 after that
42 %----------------------------------------------------------------------
47 gridFileName = [path1, 'perp_grid.plt'];
48 dataFileName = [path1, 'OutPSD.dat'];
49 %%%%%%%%%%%%%%%%%%%%%%%%%%%
50 % dataFileName = [path1, 'OutPSD_rad.dat'];
51 % gridFileName = [path1, 'rad_grid.plt'];
53 if (gridFileName(length(gridFileName)-3:length(gridFileName)) ==
'.dat')
54 [alpha, epc, L, pc] = load_plt(gridFileName,
'squeeze',
'permute');
56 [L, epc, alpha, pc] = load_plt(gridFileName,
'squeeze',
'permute');
59 target_L = L.arr(:,1,1);
60 [new_L, new_epc, new_alpha] = meshgrid(target_L, target_epc, target_alpha/180*pi);
61 new_pc =
pfunc(new_epc);
63 PSD = load_plt(dataFileName,
'permute',
'n_zones', 9999);%,
'skip_zones', 3);%
66 oneDimDataFileName = [path1, 'out1d.dat'];
67 [time, Kp, Bf, Lpp, Bw2] = load_plt(oneDimDataFileName);
68 time_arr =
PSD.time+DOY_start;
70 for it = 1: size(
PSD.
arr,1) %time.size1
71 Flux(it, :, :, :) = squeeze(
PSD.
arr(it, :, :, :)) .* pc.arr(:, :, :) .* pc.arr(:, :, :);
73 new_Flux_arr{pic}.arr = fast_interpolate(L.arr, epc.arr, alpha.arr, Flux, new_L, new_epc, new_alpha,
'log');
74 new_Flux = new_Flux_arr{pic}.arr;
75 eval([
'save ', [path1,
'Flux_E=', num2str(target_epc),
'_alpha', num2str(target_alpha)],
' new_L time_arr new_Flux'])
78 %--------------------------------------------------------------------------
80 %time_arr = time_arr + datenum('21-04-2001', 'dd-mm-yyyy') - datenum('00-04-2001', 'dd-mm-yyyy');
83 %PSD_max = 10^
max(
max(
max(
max(new_PSD([3:time.size1],:,:,:))))) - 0.2;
85 Flux_min = 10^c_min;%Flux_max - 10^delta_color;
86 new_Flux(find(new_Flux<Flux_min))=Flux_min;
87 log_Flux = log10(new_Flux);
89 f1 = figure(
'Units',
'inches',
'Position', [0.5 0.5 7 4], ...
90 'DefaultAxesFontName',
'times new roman',
'DefaultAxesFontWeight',
'bold',
'DefaultAxesFontSize', 12,
'DefaultTextFontWeight',
'bold',
'DefaultTextFontSize', 12,...
91 'PaperOrientation',
'portrait',
'PaperPositionMode',
'auto',
'DefaultLineLineWidth', 2);%,
'PaperPosition', [1 1 3.8 9]);%
97 contourf(time_arr, new_L, log_Flux', n_colors);
101 caxis([c_min c_max]);
103 text(
max(xlim), (min(ylim) +
max(ylim))/2, 'Flux, log(
#/sr/s/cm^2/KeV)', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
104 %text(
max(xlim), (min(ylim) +
max(ylim))/2,
'-3 \cdot 10^{-8} (c/MeV/cm)^3',
'Rotation', 90,
'VerticalAlignment',
'top',
'HorizontalAlignment',
'center',
'FontWeight',
'bold');
105 %xlabel(
'Time, day of Apr. 2001');
106 xlabel(
'Time, days');
108 if (target_alpha >= 88)
109 text(
max(xlim), min(ylim), {['Energy = ', num2str(target_epc), ' MeV, equatorial pitch angle ']},
'VerticalAlignment',
'bottom',
'HorizontalAlignment',
'right',
'FontWeight',
'bold',
'color', TextColor);
111 text(
max(xlim), min(ylim), {['Energy = ', num2str(target_epc), ' MeV, Pitch angle = ', num2str(target_alpha), '^o ']},
'VerticalAlignment',
'bottom',
'HorizontalAlignment',
'right',
'FontWeight',
'bold',
'color', TextColor);
115 plot(time.arr, Lpp.arr,
'w');
116 plot(time.arr, Kp.arr,
'--w');
118 save_picture(gca, [path1,
'Flux_E=', num2str(target_epc),
'_alpha', num2str(target_alpha),
'_', num2str(c_max)],
'jpg',
'pdf',
'epsc',
'fig');
121 plot(time.arr, Lpp.arr, ':', 'color', 'red', 'Linewidth', 2);
123 plot(time.arr, Kp.arr, '--', 'color', 'blue', 'Linewidth', 2);
124 title([{path1}, {
' '}]);
125 set(
gcf,
'Name', path1);
133 Flux_L_top = squeeze(
PSD.arr(:,
PSD.size1, :,
PSD.size3));
135 Flux_L_top(find(Flux_L_top<Flux_min))=Flux_min;
136 log_Flux = log10(Flux_L_top);
138 f1 = figure('Units', 'inches', 'Position', [0.5 0.5 7 4], ...
139 'DefaultAxesFontName', 'times new roman', 'DefaultAxesFontWeight', 'bold', 'DefaultAxesFontSize', 12, 'DefaultTextFontWeight', 'bold', 'DefaultTextFontSize', 12,...
140 'PaperOrientation', 'portrait', 'PaperPositionMode', 'auto', 'DefaultLineLineWidth', 2);%, 'PaperPosition', [1 1 3.8 9]);%
146 contourf(time_arr, squeeze(epc.arr(
PSD.size1, :,
PSD.size3)), log_Flux', n_colors);
149 set(gca, 'yscale', 'log');
151 caxis([c_min c_max]);
153 text(
max(xlim), 10.^((log10(min(ylim)) + log10(
max(ylim)))/2), 'Flux, log(
#/sr/s/cm^2/KeV)', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
154 %text(
max(xlim), (min(ylim) +
max(ylim))/2,
'-3 \cdot 10^{-8} (c/MeV/cm)^3',
'Rotation', 90,
'VerticalAlignment',
'top',
'HorizontalAlignment',
'center',
'FontWeight',
'bold');
155 %xlabel(
'Time, day of Apr. 2001');
156 xlabel(
'Time, days');
157 ylabel(
'Energy, MeV');
159 save_picture(gca, [path1,
'Flux_L_top'],
'jpg',
'pdf',
'epsc',
'fig');
169 plot(new_L, new_Flux(1, :));
171 plot(new_L, new_Flux(1, :));
176 f1 = figure(
'Units',
'inches',
'Position', [0.5 0.5 2 2], ...
177 'DefaultAxesFontName',
'times new roman',
'DefaultAxesFontWeight',
'bold',
'DefaultAxesFontSize', 12,
'DefaultTextFontWeight',
'bold',
'DefaultTextFontSize', 12,...
178 'PaperOrientation',
'portrait',
'PaperPositionMode',
'auto',
'DefaultLineLineWidth', 2);%,
'PaperPosition', [1 1 3.8 9]);%
179 plot(time.arr, Bf.arr,
'k',
'LineWidth', 2);
182 ylabel('Boundary flux')
183 save_picture(gca, [path1, 'Bf'], 'jpg', 'pdf', 'epsc', 'fig');
double max(double v1, double v2)
Return maximum.
double pfunc(double K)
Computation of moumentum from Kinetic energy .
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...
Matrix3D< double > arr
array of PSD values
Phase Space Density class.
Parameters_structure parameters
Parameters structure, with all parameters from the parameters.ini file. The default parameters define...