1 addpath(
'./Various_functions')
7 {
'../Execute/Output/'},
15 output_formats = [{'fig'}, {'epsc'}, {'pdf'}];
41 mc0 = 0.511;%/10^8/10^8;
45 target_K = target_J./sqrt(8.*mc0.*target_mu);
47 %target_J = target_K*sqrt(8*m0*target_mu);
51 units_constant = 1;%3*10^-8*(c/MeV/cm)^3 after that
53 %-------------------------------------------------------------------------
55 path_file = path_file_arr{1};
57 % gridFileName = [path_file, 'perp_grid.plt'];
58 gridFileName = [path_file, 'rad_grid.plt'];
60 oneDimDataFileName = [path_file, 'out1d.dat'];
62 if (gridFileName(length(gridFileName)-3:length(gridFileName)) ==
'.dat')
63 [alpha, epc, L, pc] = load_plt(gridFileName,
'squeeze',
'permute');
65 [L, epc, alpha, pc] = load_plt(gridFileName,
'squeeze',
'permute');
69 for il = 1:size(new_L)
70 % wrong RHS = (target_J*speed_of_light/earth_rad) * sqrt(new_L(il)) / sqrt(8.0 * B0(1) * mc2 * target_mu);
71 RHS = (target_J) * sqrt(new_L(il)) / sqrt(8.0 * B0(1) * mc2 * target_mu);
72 %%% ->> RHS = target_K * sqrt(new_L(il)) / sqrt(B0(1));
73 new_alpha(il) = fzero(@(alpha)
Y(alpha)/sin(alpha) - RHS, 1.5);
74 new_pc(il) = sqrt(2.0 * target_mu * mc2 * B0(new_L(il))) / sin(new_alpha(il));
76 %%%% RHS = (target_J*speed_of_light/earth_rad) * sqrt(new_L(il)) / sqrt(8.0 * B0(1) * mc2 * target_mu);
79 new_epc =
Kfunc(new_pc);
81 %f1 = figure(
'Units',
'inches',
'Position', [0.5 0.5 8 7], ...
82 %
'DefaultAxesFontName',
'times new roman',
'DefaultAxesFontWeight',
'bold',
'DefaultAxesFontSize', 12,...
83 %
'PaperOrientation',
'portrait',
'PaperPositionMode',
'auto',
'DefaultLineLineWidth', 2);%,
'PaperPosition', [1 1 3.8 9]);%
87 %subplot1(2, 1,
'Gap', [0.07, 0.03],
'Max', [0.935, 0.95]);
89 %-------------------------------------------------------------------------
91 for i = 1:length(path_file_arr)
92 path_file = path_file_arr{i}
94 f1 = figure(
'Units',
'inches',
'Position', [0.5 0.5 6 4], ...
95 'DefaultAxesFontName',
'times new roman',
'DefaultAxesFontWeight',
'bold',
'DefaultAxesFontSize', 12,...
96 'PaperOrientation',
'portrait',
'PaperPositionMode',
'auto',
'DefaultLineLineWidth', 2);%,
'PaperPosition', [1 1 3.8 9]);%
99 cmap(1, :) = [0, 0, 0];
102 % dataFileName = [path_file, 'OutPSD.dat'];
103 dataFileName = [path_file, 'OutPSD_rad.dat'];
105 [
PSD] = load_plt(dataFileName,
'permute',
'nosqueeze',
'skip_zones', 0,
'n_zones', 9999);%
108 [time, Kp, Bf, Lpp] = load_plt(oneDimDataFileName);
112 [new_PSD] = fast_interpolate(L.arr, epc.arr, alpha.arr,
PSD.
arr, new_L, new_epc, new_alpha,
'log');
113 % [new_PSD] = interpolate(L.arr, epc.arr, alpha.arr,
PSD.
arr, new_L, new_epc, new_alpha,
'log');
115 eval([
'save ', [path_file,
'PSD_mu=', num2str(target_mu),
'_J', num2str(target_J)],
' new_L time_arr new_PSD'])
117 %--------------------------------------------------------------------------
125 %PSD_max = 10^
max(
max(
max(
max(new_PSD([3:time.size1],:,:,:))))) - 0.2;
127 PSD_min = 10^c_min;%PSD_max - 10^delta_color;
128 new_PSD(find(new_PSD<PSD_min))=PSD_min;
129 log_PSD = log10(new_PSD);
130 % cmap = colormap(jet);
132 % cmap(1,:) = [0 0 0];
134 contourf(time_arr, new_L, log_PSD', n_colors);
139 caxis([c_min c_max]);
141 if (pic_n == 1 || pic_n == 3 || pic_n == 5) ylabel('L'); end;
142 if (pic_n == 5 || pic_n == 6) xlabel('Time, days'); end;
143 if (pic_n == 1 || pic_n == 3 || pic_n == 5)
144 ax_pos = get(gca, 'Position');
145 cl = colorbar('West');
146 cl_pos = get(cl, 'Position');
147 cl_x_pos = ax_pos(1) + ax_pos(3) + 0.025;%0.015; %0.925
148 set(cl, 'Position', [cl_x_pos, cl_pos(2), cl_pos(3), cl_pos(4)]);
149 text(
max(xlim), (min(ylim) +
max(ylim))/2, ' 3 \cdot 10^{-8} (c/cm/MeV)^3', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold', 'FontSize', 10);
153 plot(time.arr, Lpp.arr,
':',
'color',
'white',
'Linewidth', 1);
154 plot(time.arr, Kp.arr,
':',
'color',
'white',
'Linewidth', 1);
156 set(
gcf,
'Name', path_file);
162 save_picture(gca, [path_file, 'PSD_', num2str(target_mu), '_J', num2str(target_J), '_K', num2str(target_K)], 'jpg', 'pdf', 'epsc', 'fig');
167 %subplot1(2); delete(gca);
168 %subplot1(3); delete(gca);
169 %subplot1(4); delete(gca);
177 plot(squeeze(epc.arr(
PSD.size1, :,
PSD.size3)), squeeze(
PSD.arr(1,
PSD.size1, :,
PSD.size3)));
178 set(gca, 'xscale', 'log');
179 set(gca, 'yscale', 'log');
184 contourf(time_arr, squeeze(epc.arr(
PSD.size1, :,
PSD.size3)), log10(squeeze(
PSD.arr(:,
PSD.size1, :,
PSD.size3)))', 50);
186 set(gca, 'yscale', 'log');
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...
double pc2mu(double L, double pc, double alpha)
Convert pc to μ.
Matrix3D< double > arr
array of PSD values
Phase Space Density class.
double Y(double alpha)
Y(α) function.
Parameters_structure parameters
Parameters structure, with all parameters from the parameters.ini file. The default parameters define...
double Kfunc(double pc)
Computation of Kinetic energy from given momentum .