1 addpath(
'./Various_functions')
7 {
'../Execute/Output/'},
15 output_formats = [{'fig'}, {'epsc'}, {'pdf'}];
37 mc0 = 0.511;%/10^8/10^8;
41 target_K = target_J./sqrt(8.*mc0.*target_mu);
43 %target_J = target_K*sqrt(8*m0*target_mu);
47 units_constant = 1;%3*10^-8*(c/MeV/cm)^3 after that
49 %-------------------------------------------------------------------------
51 path_file = path_file_arr{1};
53 gridFileName = [path_file, 'perp_grid.plt'];
54 % gridFileName = [path_file, 'rad_grid.plt'];
56 oneDimDataFileName = [path_file, 'out1d.dat'];
58 if (gridFileName(length(gridFileName)-3:length(gridFileName)) ==
'.dat')
59 [alpha, epc, L, pc] = load_plt(gridFileName,
'squeeze',
'permute');
61 [L, epc, alpha, pc] = load_plt(gridFileName,
'squeeze',
'permute');
65 for il = 1:size(new_L)
66 % wrong RHS = (target_J*speed_of_light/earth_rad) * sqrt(new_L(il)) / sqrt(8.0 * B0(1) * mc2 * target_mu);
67 RHS = (target_J) * sqrt(new_L(il)) / sqrt(8.0 * B0(1) * mc2 * target_mu);
68 %%% ->> RHS = target_K * sqrt(new_L(il)) / sqrt(B0(1));
69 new_alpha(il) = fzero(@(alpha)
Y(alpha)/sin(alpha) - RHS, 1.5);
70 new_pc(il) = sqrt(2.0 * target_mu * mc2 * B0(new_L(il))) / sin(new_alpha(il));
72 %%%% RHS = (target_J*speed_of_light/earth_rad) * sqrt(new_L(il)) / sqrt(8.0 * B0(1) * mc2 * target_mu);
75 new_epc =
Kfunc(new_pc);
77 %f1 = figure(
'Units',
'inches',
'Position', [0.5 0.5 8 7], ...
78 %
'DefaultAxesFontName',
'times new roman',
'DefaultAxesFontWeight',
'bold',
'DefaultAxesFontSize', 12,...
79 %
'PaperOrientation',
'portrait',
'PaperPositionMode',
'auto',
'DefaultLineLineWidth', 2);%,
'PaperPosition', [1 1 3.8 9]);%
83 %subplot1(2, 1,
'Gap', [0.07, 0.03],
'Max', [0.935, 0.95]);
85 %-------------------------------------------------------------------------
87 for i = 1:length(path_file_arr)
88 path_file = path_file_arr{i}
90 f1 = figure(
'Units',
'inches',
'Position', [0.5 0.5 6 4], ...
91 'DefaultAxesFontName',
'times new roman',
'DefaultAxesFontWeight',
'bold',
'DefaultAxesFontSize', 12,...
92 'PaperOrientation',
'portrait',
'PaperPositionMode',
'auto',
'DefaultLineLineWidth', 2);%,
'PaperPosition', [1 1 3.8 9]);%
94 dataFileName = [path_file, 'OutPSD.dat'];
95 % dataFileName = [path_file, 'OutPSD_rad.dat'];
97 [
PSD] = load_plt(dataFileName,
'permute',
'nosqueeze',
'skip_zones', 0,
'n_zones', 9999);%
100 [time, Kp, Bf, Lpp] = load_plt(oneDimDataFileName);
104 [new_PSD] = fast_interpolate(L.arr, epc.arr, alpha.arr,
PSD.
arr, new_L, new_epc, new_alpha,
'log');
105 % [new_PSD] = interpolate(L.arr, epc.arr, alpha.arr,
PSD.
arr, new_L, new_epc, new_alpha,
'log');
107 eval([
'save ', [path_file,
'PSD_mu=', num2str(target_mu),
'_J', num2str(target_J)],
' new_L time_arr new_PSD'])
109 %--------------------------------------------------------------------------
117 %PSD_max = 10^
max(
max(
max(
max(new_PSD([3:time.size1],:,:,:))))) - 0.2;
119 PSD_min = 10^c_min;%PSD_max - 10^delta_color;
120 new_PSD(find(new_PSD<PSD_min))=PSD_min;
121 log_PSD = log10(new_PSD);
122 % cmap = colormap(jet);
124 % cmap(1,:) = [0 0 0];
126 contourf(time_arr, new_L, log_PSD', n_colors);
131 caxis([c_min c_max]);
133 if (pic_n == 1 || pic_n == 3 || pic_n == 5) ylabel('L'); end;
134 if (pic_n == 5 || pic_n == 6) xlabel('Time, days'); end;
135 if (pic_n == 1 || pic_n == 3 || pic_n == 5)
136 ax_pos = get(gca, 'Position');
137 cl = colorbar('West');
138 cl_pos = get(cl, 'Position');
139 cl_x_pos = ax_pos(1) + ax_pos(3) + 0.025;%0.015; %0.925
140 set(cl, 'Position', [cl_x_pos, cl_pos(2), cl_pos(3), cl_pos(4)]);
141 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);
145 plot(time.arr, Lpp.arr,
':',
'color',
'white',
'Linewidth', 1);
146 plot(time.arr, Kp.arr,
':',
'color',
'white',
'Linewidth', 1);
148 set(
gcf,
'Name', path_file);
154 save_picture(gca, [path_file, 'PSD_', num2str(target_mu), '_J', num2str(target_J), '_K', num2str(target_K)], 'jpg', 'pdf', 'epsc', 'fig');
159 %subplot1(2); delete(gca);
160 %subplot1(3); delete(gca);
161 %subplot1(4); delete(gca);
169 plot(squeeze(epc.arr(
PSD.size1, :,
PSD.size3)), squeeze(
PSD.arr(1,
PSD.size1, :,
PSD.size3)));
170 set(gca, 'xscale', 'log');
171 set(gca, 'yscale', 'log');
176 contourf(time_arr, squeeze(epc.arr(
PSD.size1, :,
PSD.size3)), log10(squeeze(
PSD.arr(:,
PSD.size1, :,
PSD.size3)))', 50);
178 set(gca, 'yscale', 'log');
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...
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 .