VERB_code_2.3
Plot_PSD_L_time.m
1 addpath('./Various_functions')
2 
3 clear all;
4 %close all
5 
6 path_file_arr = [
7  {'../Execute/Output/'},
8  ];
9 
10 pic_n = 1;
11 TextColor = 'white';
12 TextColor1 = 'white';
13 
14 % output parameters:
15 output_formats = [{'fig'}, {'epsc'}, {'pdf'}];
16 FontSize = 9;
17 n_colors = 50;
18 
19 target_mu = [850];
20 c_min = 0;
21 c_max = 3.5;
22 
23 % target_mu = [2000];
24 % c_min = -2.0;
25 % c_max = 1.5;
26 
27 % target_mu = [2000];
28 % c_min = 0;
29 % c_max = 1.3;
30 
31 
32 % target_mu = [2000];
33 % c_min = -2.5;
34 % c_max = 0.8;
35 
36 %target_K = [0.025];
37 
38 target_J = [0.01];
39 %target_J = [1];
40 
41 mc0 = 0.511;%/10^8/10^8;
42 speed_of_light = 3e8;
43 earth_rad = 6.4e6;
44 mc2 = 0.511;
45 target_K = target_J./sqrt(8.*mc0.*target_mu);
46 
47 %target_J = target_K*sqrt(8*m0*target_mu);
48 
49 
50 
51 units_constant = 1;%3*10^-8*(c/MeV/cm)^3 after that
52 
53 %-------------------------------------------------------------------------
54 
55  path_file = path_file_arr{1};
56 
57  gridFileName = [path_file, 'perp_grid.plt'];
58 % gridFileName = [path_file, 'rad_grid.plt'];
59 
60  oneDimDataFileName = [path_file, 'out1d.dat'];
61 
62  if (gridFileName(length(gridFileName)-3:length(gridFileName)) == '.dat')
63  [alpha, epc, L, pc] = load_plt(gridFileName, 'squeeze', 'permute');
64  else
65  [L, epc, alpha, pc] = load_plt(gridFileName, 'squeeze', 'permute');
66  end;
67 
68  new_L = L.arr(:,1,1);
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));
75 
76 %%%% RHS = (target_J*speed_of_light/earth_rad) * sqrt(new_L(il)) / sqrt(8.0 * B0(1) * mc2 * target_mu);
77 
78  end
79  new_epc = Kfunc(new_pc);
80 
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]);%
84 
85 
86  %subplot(3,2,1);
87  %subplot1(2, 1, 'Gap', [0.07, 0.03], 'Max', [0.935, 0.95]);
88 
89  %-------------------------------------------------------------------------
90 
91  for i = 1:length(path_file_arr)
92  path_file = path_file_arr{i}
93 
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]);%
97 
98  cmap = colormap;
99  cmap(1, :) = [0, 0, 0];
100  colormap(cmap);
101 
102  dataFileName = [path_file, 'OutPSD.dat'];
103 % dataFileName = [path_file, 'OutPSD_rad.dat'];
104 
105  [PSD] = load_plt(dataFileName, 'permute', 'nosqueeze', 'skip_zones', 3, 'n_zones', 9999);%
106  PSD.arr = PSD.arr*units_constant;
107 
108  [time, Kp, Bf, Lpp] = load_plt(oneDimDataFileName);
109 
110  time_arr = PSD.time;
111 
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');
114 
115  eval(['save ', [path_file, 'PSD_mu=', num2str(target_mu), '_J', num2str(target_J)], ' new_L time_arr new_PSD'])
116 
117  %--------------------------------------------------------------------------
118 
119  pic_n = 1;
120 
121  %subplot(3,2,pic_n);
122  %subplot1(pic_n);
123 
124  %PSD_max = 10^c_max
125  %PSD_max = 10^max(max(max(max(new_PSD([3:time.size1],:,:,:))))) - 0.2;
126  %c_max = log10(max(max(max(max(new_PSD(:,:,:,:))))));
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);
131 % cmap(1,:);
132 % cmap(1,:) = [0 0 0];
133 % colormap(cmap)
134  contourf(time_arr, new_L, log_PSD', n_colors);
135  shading flat
136  axis tight
137  ylim([1 7]);
138 
139  caxis([c_min c_max]);
140 
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);
150  end;
151 
152  hold on
153  plot(time.arr, Lpp.arr, ':', 'color', 'white', 'Linewidth', 1);
154  plot(time.arr, Kp.arr, ':', 'color', 'white', 'Linewidth', 1);
155 
156  set(gcf, 'Name', path_file);
157 
158  xlabel('Time, days')
159  ylabel('L')
160 
161 
162  save_picture(gca, [path_file, 'PSD_', num2str(target_mu), '_J', num2str(target_J), '_K', num2str(target_K)], 'jpg', 'pdf', 'epsc', 'fig');
163 
164  end
165 
166 
167 %subplot1(2); delete(gca);
168 %subplot1(3); delete(gca);
169 %subplot1(4); delete(gca);
170 
171 
172 
173 return
174 
175 
176 figure
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');
180 ylim([1e-5 1e8])
181 
182 
183 figure
184 contourf(time_arr, squeeze(epc.arr(PSD.size1, :, PSD.size3)), log10(squeeze(PSD.arr(:, PSD.size1, :, PSD.size3)))', 50);
185 shading flat
186 set(gca, 'yscale', 'log');
187 caxis([0 6])
188 colorbar
189 
190 
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
double Y(double alpha)
Y(α) function.
Parameters_structure parameters
Parameters structure, with all parameters from the parameters.ini file. The default parameters define...
Definition: Main.cpp:185
double Kfunc(double pc)
Computation of Kinetic energy from given momentum .