VERB_code_2.3
fast_interpolate.m
1 function [new_arr] = interpolate(L_arr, epc_arr, alpha_arr, f_arr, new_L, new_epc, new_alpha, varargin)
2 
3  loginterp = false;
4  if (size(varargin) > 0)
5  if (strcmp(varargin(1), 'log'))
6  loginterp = true;
7  for it = 1:size(f_arr, 1)
8  f_arr(it, :, :, :) = log10(f_arr(it, :, :, :));
9  end
10  end
11  end
12 
13  for it = 1:size(f_arr, 1)
14  ['interpolation step = ', num2str(it)]
15 
16  if (loginterp)
17  new_arr(it, 1) = -22; %Trick - PSD=0 for L=1
18  else
19  new_arr(it, 1) = 1e-22; %Trick - PSD=0 for L=1
20  end
21  for L_ind = 2:size(f_arr, 2)
22 
23  try
24 
25  % using the fact, that alpha is undependent of epc!
26  alpha_l_ind = max(find(alpha_arr(L_ind,1,:)<new_alpha(L_ind)));
27  alpha_r_ind = min(find(alpha_arr(L_ind,1,:)>new_alpha(L_ind)));
28 
29  % finding bottom and top indexes of epc for each alpha
30  epc_lb_ind = max(find(epc_arr(L_ind,:,alpha_l_ind)<new_epc(L_ind)));
31  epc_lt_ind = min(find(epc_arr(L_ind,:,alpha_l_ind)>new_epc(L_ind)));
32  epc_rb_ind = max(find(epc_arr(L_ind,:,alpha_r_ind)<new_epc(L_ind)));
33  epc_rt_ind = min(find(epc_arr(L_ind,:,alpha_r_ind)>new_epc(L_ind)));
34 
35  %['L=', num2str(L_arr(L_ind,1,1)), '; al=', num2str(alpha_l_ind), '; ar=', num2str(alpha_r_ind),'; epc_lb=', num2str(epc_lb_ind), '; epc_lt=', num2str(epc_lt_ind), '; epc_rb=', num2str(epc_rb_ind), '; epc_rt=', num2str(epc_rt_ind)]
36 
37  % interpolating
38 % double a = (x - m_x[i - 1]) / (m_x[i] - m_x[i - 1]);
39 % return m_y[i - 1] + a * (m_y[i] - m_y[i - 1]);
40  % left side
41  fb = f_arr(it, L_ind, epc_lb_ind, alpha_l_ind);
42  ft = f_arr(it, L_ind, epc_lt_ind, alpha_l_ind);
43  epcb = epc_arr(L_ind, epc_lb_ind, alpha_l_ind);
44  epct = epc_arr(L_ind, epc_lt_ind, alpha_l_ind);
45 
46  a = ( new_epc(L_ind) - epcb ) / ( epct - epcb );
47  new_arr_l = fb + a * ( ft - fb );
48 
49  % right side
50  fb = f_arr(it, L_ind, epc_rb_ind, alpha_r_ind);
51  ft = f_arr(it, L_ind, epc_rt_ind, alpha_r_ind);
52  epcb = epc_arr(L_ind, epc_rb_ind, alpha_r_ind);
53  epct = epc_arr(L_ind, epc_rt_ind, alpha_r_ind);
54 
55  a = ( new_epc(L_ind) - epcb ) / ( epct - epcb );
56  new_arr_r = fb + a * ( ft - fb );
57 
58  % between sides
59  alphab = alpha_arr(L_ind, epc_lb_ind, alpha_l_ind);
60  alphat = alpha_arr(L_ind, epc_rb_ind, alpha_r_ind);
61 
62  a = ( new_alpha(L_ind) - alphab ) / ( alphat - alphab );
63  new_arr(it, L_ind) = new_arr_l + a * ( new_arr_r - new_arr_l );
64 
65  % temp!!!
66  %new_arr(it, L_ind) = f_arr(it, L_ind, epc_rt_ind, alpha_r_ind);
67  catch
68  err = lasterror;
69  [err.message, '; L_ind=', num2str(L_ind)]
70  if (loginterp)
71  new_arr(it, L_ind) = -22; %Trick - PSD=0 for L=1
72  else
73  new_arr(it, L_ind) = 1e-22; %Trick - PSD=0 for L=1
74  end
75  end
76  end
77 
78  end
79 
80  if (loginterp)
81  new_arr = 10.^new_arr;
82  end
83 
84 end
double max(double v1, double v2)
Return maximum.
Phase Space Density class.
Definition: PSD.h:48