VERB_code_2.3
Plot_Flux_L_time.m
1 %%
2 clear all;
3 %close all;
4 addpath('./Various_functions')
5 
6 DOY_start=293; % 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 = 80;
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 =0 ;
35  c_max = 5.5;
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', 3);%
64  PSD.arr = PSD.arr*units_constant;
65 
66  oneDimDataFileName = [path1, 'out1d.dat'];
67  [time, Kp, Bf, Lpp, Bw2] = load_plt(oneDimDataFileName);
68  time_arr = PSD.time+DOY_start;
69 
70  for it = 1: size(PSD.arr,1) %time.size1
71  Flux(it, :, :, :) = squeeze(PSD.arr(it, :, :, :)) .* pc.arr(:, :, :) .* pc.arr(:, :, :);
72  end
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'])
76 
77 
78  %--------------------------------------------------------------------------
79 
80  %time_arr = time_arr + datenum('21-04-2001', 'dd-mm-yyyy') - datenum('00-04-2001', 'dd-mm-yyyy');
81 
82  %PSD_max = 10^c_max
83  %PSD_max = 10^max(max(max(max(new_PSD([3:time.size1],:,:,:))))) - 0.2;
84  %c_max = log10(max(max(max(max(new_PSD(:,:,:,:))))));
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);
88 %%
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]);%
92 
93  cmap = colormap;
94  cmap(1,:);
95  cmap(1,:) = [0 0 0];
96  colormap(cmap)
97  contourf(time_arr, new_L, log_Flux', n_colors);
98  shading flat
99  axis tight
100  ylim([1 7]);
101  caxis([c_min c_max]);
102  colorbar
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');
107  ylabel('L*');
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);
110  else
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);
112  end
113 
114  hold on
115  plot(time.arr, Lpp.arr, 'w');
116  plot(time.arr, Kp.arr, '--w');
117 
118  save_picture(gca, [path1, 'Flux_E=', num2str(target_epc), '_alpha', num2str(target_alpha), '_', num2str(c_max)], 'jpg', 'pdf', 'epsc', 'fig');
119 
120  figure()
121  plot(time.arr, Lpp.arr, ':', 'color', 'red', 'Linewidth', 2);
122  hold on
123  plot(time.arr, Kp.arr, '--', 'color', 'blue', 'Linewidth', 2);
124  title([{path1}, {' '}]);
125  set(gcf, 'Name', path1);
126 
127 
128  %close(f1)
129  if (false)
130 
131  c_min = -2;
132  c_max = 7.4
133  Flux_L_top = squeeze(PSD.arr(:, PSD.size1, :, PSD.size3));
134  Flux_min = 10^c_min;
135  Flux_L_top(find(Flux_L_top<Flux_min))=Flux_min;
136  log_Flux = log10(Flux_L_top);
137 
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]);%
141 
142  cmap = colormap;
143  cmap(1,:);
144  cmap(1,:) = [0 0 0];
145  colormap(cmap)
146  contourf(time_arr, squeeze(epc.arr(PSD.size1, :, PSD.size3)), log_Flux', n_colors);
147  shading flat
148  axis tight
149  set(gca, 'yscale', 'log');
150  ylim([1e-1 6e0]);
151  caxis([c_min c_max]);
152  colorbar
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');
158 
159  save_picture(gca, [path1, 'Flux_L_top'], 'jpg', 'pdf', 'epsc', 'fig');
160 
161  end
162 
163 end;
164 
165 
166 return
167 
168 figure
169 plot(new_L, new_Flux(1, :));
170 hold on;
171 plot(new_L, new_Flux(1, :));
172 
173 
174 return
175 
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);
180 ylim([0 1.1])
181 xlabel('Time, days')
182 ylabel('Boundary flux')
183 save_picture(gca, [path1, 'Bf'], 'jpg', 'pdf', 'epsc', 'fig');
184 
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