VERB_code_2.3
Plot_Flux_L_time.m
1 %%
2 clear all;
3 %close all;
4 addpath('./Various_functions')
5 
6 DOY_start=0; % day of the year of the first day of simulation
7 
8 
9 path_arr = [
10  {'../Execute/Output/'},
11  ];
12 
13 pic_n = 3;
14 TextColor = 'white';
15 TextColor1 = 'white';
16 
17 % output parameters:
18 output_formats = [{'fig'}, {'epsc'}, {'pdf'}];
19 FontSize = 9;
20 n_colors = 90;
21 
22 
23 for pic = 1:length(path_arr)
24  path1 = path_arr{pic};
25 
26  %c_min = 1;
27  %c_max = 6;
28  %target_epc = [0.1];
29 
30  %c_min = 1;
31  %c_max = 7;
32  %target_epc = [0.5];
33 
34  c_min = -1;
35  c_max = 1;
36  target_epc = [1];
37 
38  target_alpha = [85];
39 
40  units_constant = 1;%3*10^-2; % units 10^-6*(c/MeV/cm)^3 after that
41 
42  %----------------------------------------------------------------------
43  %---
44 
45 
46 
47  gridFileName = [path1, 'perp_grid.plt'];
48  dataFileName = [path1, 'OutPSD.dat'];
49  %%%%%%%%%%%%%%%%%%%%%%%%%%%
50  % dataFileName = [path1, 'OutPSD_rad.dat'];
51  % gridFileName = [path1, 'rad_grid.plt'];
52 
53  if (gridFileName(length(gridFileName)-3:length(gridFileName)) == '.dat')
54  [alpha, epc, L, pc] = load_plt(gridFileName, 'squeeze', 'permute');
55  else
56  [L, epc, alpha, pc] = load_plt(gridFileName, 'squeeze', 'permute');
57  end;
58 
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);
62 
63  PSD = load_plt(dataFileName, 'permute', 'n_zones', 9999, 'skip_zones', 0);%
64 
65  PSD.arr = PSD.arr*units_constant;
66 
67  oneDimDataFileName = [path1, 'out1d.dat'];
68  [time, Kp, Bf, Lpp, Bw2] = load_plt(oneDimDataFileName);
69  time_arr = PSD.time+DOY_start;
70 return
71  for it = 1: size(PSD.arr,1) %time.size1
72  Flux(it, :, :, :) = squeeze(PSD.arr(it, :, :, :)) .* pc.arr(:, :, :) .* pc.arr(:, :, :);
73  end
74  new_Flux_arr{pic}.arr = fast_interpolate(L.arr, epc.arr, alpha.arr, Flux, new_L, new_epc, new_alpha, 'log');
75  new_Flux = new_Flux_arr{pic}.arr;
76  eval(['save ', [path1, 'Flux_E=', num2str(target_epc), '_alpha', num2str(target_alpha)], ' new_L time_arr new_Flux'])
77 
78 
79  %--------------------------------------------------------------------------
80 
81  %time_arr = time_arr + datenum('21-04-2001', 'dd-mm-yyyy') - datenum('00-04-2001', 'dd-mm-yyyy');
82 
83  %PSD_max = 10^c_max
84  %PSD_max = 10^max(max(max(max(new_PSD([3:time.size1],:,:,:))))) - 0.2;
85  %c_max = log10(max(max(max(max(new_PSD(:,:,:,:))))));
86  Flux_min = 10^c_min;%Flux_max - 10^delta_color;
87  new_Flux(new_Flux<Flux_min)=Flux_min;
88  log_Flux = log10(new_Flux);
89 %%
90  f1 = figure('Units', 'inches', 'Position', [0.5 0.5 7 4], ...
91  'DefaultAxesFontName', 'times new roman', 'DefaultAxesFontWeight', 'bold', 'DefaultAxesFontSize', 12, 'DefaultTextFontWeight', 'bold', 'DefaultTextFontSize', 12,...
92  'PaperOrientation', 'portrait', 'PaperPositionMode', 'auto', 'DefaultLineLineWidth', 2);%, 'PaperPosition', [1 1 3.8 9]);%
93 
94  cmap = colormap;
95  cmap(1,:);
96  cmap(1,:) = [0 0 0];
97  colormap(cmap)
98  contourf(time_arr, new_L, 10.^log_Flux', n_colors);
99  shading flat
100  axis tight
101  ylim([1 7]);
102  caxis(10.^[c_min c_max]);
103  colorbar
104  text(max(xlim), (min(ylim) + max(ylim))/2, 'Flux, log(#/sr/s/cm^2/KeV)', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
105  %text(max(xlim), (min(ylim) + max(ylim))/2, '-3 \cdot 10^{-8} (c/MeV/cm)^3', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
106  %xlabel('Time, day of Apr. 2001');
107  xlabel('Time, days');
108  ylabel('L*');
109  if (target_alpha >= 88)
110  text(max(xlim), min(ylim), {['Energy = ', num2str(target_epc), ' MeV, equatorial pitch angle ']}, 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'color', TextColor);
111  else
112  text(max(xlim), min(ylim), {['Energy = ', num2str(target_epc), ' MeV, Pitch angle = ', num2str(target_alpha), '^o ']}, 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'color', TextColor);
113  end
114 
115  hold on
116  plot(time.arr, Lpp.arr, 'w');
117  plot(time.arr, Kp.arr, '--w');
118 
119  save_picture(gca, [path1, 'Flux_E=', num2str(target_epc), '_alpha', num2str(target_alpha), '_', num2str(c_max)], 'jpg', 'pdf', 'epsc', 'fig');
120 
121  figure()
122  plot(time.arr, Lpp.arr, ':', 'color', 'red', 'Linewidth', 2);
123  hold on
124  plot(time.arr, Kp.arr, '--', 'color', 'blue', 'Linewidth', 2);
125  title([{path1}, {' '}]);
126  set(gcf, 'Name', path1);
127 
128 
129  %close(f1)
130  if (false)
131 
132  c_min = -2;
133  c_max = 7.4
134  Flux_L_top = squeeze(PSD.arr(:, PSD.size1, :, PSD.size3));
135  Flux_min = 10^c_min;
136  Flux_L_top(find(Flux_L_top<Flux_min))=Flux_min;
137  log_Flux = log10(Flux_L_top);
138 
139  f1 = figure('Units', 'inches', 'Position', [0.5 0.5 7 4], ...
140  'DefaultAxesFontName', 'times new roman', 'DefaultAxesFontWeight', 'bold', 'DefaultAxesFontSize', 12, 'DefaultTextFontWeight', 'bold', 'DefaultTextFontSize', 12,...
141  'PaperOrientation', 'portrait', 'PaperPositionMode', 'auto', 'DefaultLineLineWidth', 2);%, 'PaperPosition', [1 1 3.8 9]);%
142 
143  cmap = colormap;
144  cmap(1,:);
145  cmap(1,:) = [0 0 0];
146  colormap(cmap)
147  contourf(time_arr, squeeze(epc.arr(PSD.size1, :, PSD.size3)), log_Flux', n_colors);
148  shading flat
149  axis tight
150  set(gca, 'yscale', 'log');
151  ylim([1e-1 6e0]);
152  caxis([c_min c_max]);
153  colorbar
154  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');
155  %text(max(xlim), (min(ylim) + max(ylim))/2, '-3 \cdot 10^{-8} (c/MeV/cm)^3', 'Rotation', 90, 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
156  %xlabel('Time, day of Apr. 2001');
157  xlabel('Time, days');
158  ylabel('Energy, MeV');
159 
160  save_picture(gca, [path1, 'Flux_L_top'], 'jpg', 'pdf', 'epsc', 'fig');
161 
162  end
163 
164 end;
165 
166 
167 return
168 
169 figure
170 plot(new_L, new_Flux(1, :));
171 hold on;
172 plot(new_L, new_Flux(1, :));
173 
174 
175 return
176 
177 f1 = figure('Units', 'inches', 'Position', [0.5 0.5 2 2], ...
178 'DefaultAxesFontName', 'times new roman', 'DefaultAxesFontWeight', 'bold', 'DefaultAxesFontSize', 12, 'DefaultTextFontWeight', 'bold', 'DefaultTextFontSize', 12,...
179 'PaperOrientation', 'portrait', 'PaperPositionMode', 'auto', 'DefaultLineLineWidth', 2);%, 'PaperPosition', [1 1 3.8 9]);%
180 plot(time.arr, Bf.arr, 'k', 'LineWidth', 2);
181 ylim([0 1.1])
182 xlabel('Time, days')
183 ylabel('Boundary flux')
184 %save_picture(gca, [path1, 'Bf'], 'jpg', 'pdf', 'epsc', 'fig');
185 
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...
Definition: erf.cpp:103
Matrix3D< double > arr
array of PSD values
Definition: PSD.h:53
Phase Space Density class.
Definition: PSD.h:48
Parameters_structure parameters
Parameters structure, with all parameters from the parameters.ini file. The default parameters define...
Definition: Main.cpp:185