% slit_fn
% This script takes experimentally determined slit functions (meas_fns)
% and combines them to create a total slit function for our spectrometer.

% VARIABLES:
% Input variables:
% Read from slit function data files:
% coord array giving the first pixel # (ie lowest pixel #)
% of each slit funtion
% len array containing the width (in pixels) of each slit function
% meas_fn matrix containing the measured slit functions
% rows = measured function # (ie 1:n_fns)
% columns = function value (1:length(row #))
% Other input variables:
% n_fns # of measured slit functions
% n_pix # of pixels on the CCD
% noise # noise threshold
% h_width desired half width of slit functions
% max_iterations max # of iterations to do before quitting
% Output and other variables
% centre array giving the centre of "mass" of each slit function.
% After the script has run, this will be the same for
% all pixels.
% coord Okay, this is the same variable as mentioned above. But
% after the script has run, it contains the pixel
% coordinate of the centre of each measured slit function.
% s_fn output matrix for the slit functions
% rows = pixel #
% columns = function value (1:max_len)
% max_len maximum width of measured slit functions.
% max_cent centre coord of measured slit function closest to n_pix
% min_cent centre coord of measured slit function closest to 1
%
% --> So, just to summarize, the key output of this script is contained in the
% variables "s_fn", "centre" and "max_len".
%
% INPUT FILES:
% whatever files contain the slit function data
% CALLS:
% find_centre finds optimal centre pixel for each slit function

% INITIALIZE:
clear n_fns n_pix noise i coord centre meas_fn s_fn
clear max_len lhs_len rhs_len min_temp s_loc len
clear z1 z2 a1 a2 m n

clear global len meas_fn coord

global meas_fn
global len
global h_width
global max_iterations
global coord

%n_fns = input('Please enter # of functions to fit:');
%i = 1:1:n_fns;

n_pix = 2000;
noise = 1450 + 1500 % show user the noise thresh

h_width = 50; % desired half width, in pixels, of the slit functions
max_iterations = 1;

% Read in slit function data from CSV files
dir_path = 'c:\james\analyse\slit\slit_data\';
wv_fn = [-8.5631 4884] % polynomial for converty wave centre to pixel centre
n_fns = read_csv_slit_fns(dir_path, wv_fn) % see READ_CSV_SLIT_FNS.M
i = 1:1:n_fns;
if n_fns < 1
disp('SLIT_FN.M: Error, problem reading in slit functions.')
return
else
disp(strcat('SLIT_FN.M: No. of measured slit functions = ', num2str(n_fns)))
end

whos meas_fn
len
coord

%z1 = wk1read(strcat(dir_path,'s_260.wk1'));
%z2 = wk1read(strcat(dir_path,'s_185.wk1'));
%z1 = z1';
%z2 = z2';
%len(i) = [length(z1) length(z2)];
%meas_fn(1,1:len(1)) = z1;
%meas_fn(2,1:len(2)) = z2;
%coord(i) = [260 185];
%clear z1 z2 z3

%coord(i) = [4 7 12];

% READ IN SLIT FNS. temp fix:
%k = 1:10;
%temp(k) = 3;
%meas_fn(i,k) = 3;
%len(i) = 10;
%centre(i) = 2;

% note: use the fn "wk1read('filename')" to read in Lotus (wk1) spreadsheet data

%meas_fn(1,k) = 3;
%meas_fn(2,k) = 3;
%meas_fn(1,1) = 10;
%meas_fn(2,8) = 10;

% DEFINE NOISE THRESHOLD (AND SUBTRACT NOISE?)
for m = 1:1:n_fns
for k = 1:1:len(m)
meas_fn(m,k) = meas_fn(m,k) - noise;
end
end

centre = coord; % initialize centre to values that are definitely not the centre
lower(1:n_fns) = 1;

whos meas_fn

for m = 1:1:n_fns
centre(m) = find_centre(m, centre(m), lower(m), len(m), 0);
end

disp('printing from slit_fn.m, centre is...')
centre

% Shift coords of measured functions so that coord now stores the absolute pixel
% position of the centre of each one:
coord(i) = coord(i) + centre(i);
max_cent = max(coord);
min_cent = min(coord);
% note, the above line MUST be executed before the shifting code below runs!!

% get max slit function width
max_len = max(len(i));
% initialize output slit function array (not necessary in MatLab, but speeds up the code)
s_fn = zeros(n_pix, max_len);

% To account for the possibility of measured slit functions being of different
% pixel widths, we should find the maximum function width and then shift the
% other functions so that they are centred within that pixel range. That is,
% we shift the functions so that they are all centred at the same column index
% of meas_fn:
lhs_len = max(centre); % remember, centre = RELATIVE centre (in pixels)
% coord = ABSOLUTE centre (in pixels)
rhs_len = max(max_len - centre);
shift(i) = lhs_len - centre(i);
centre(i) = centre(i) + shift(i);
clear temp
temp = meas_fn;
for m = 1:1:n_fns
for k = 1:1:len(m)
if (k + shift(m) > 0) %bandaid fix
meas_fn(m,k + shift(m)) = temp(m,k);
end
end
for k = 1:1:shift(m)
meas_fn(i,k) = 0;
end
for k = (len(i) + shift(i) + 1):1:(centre(i) + rhs_len)
meas_fn(i,k) = 0;
end
end % note to self: this code redefines array lengths a lot,
% and so it's probably not as fast as it could be if the
% array lengths were preallocated in memory. FIX LATER

% Update centre and max_len:
disp('Updated centre and max_len:')
centre = centre(1)
max_len = centre + rhs_len

% Now sort meas_fn so that the measured slit functions are organized in order
% of increasing absolute pixel centre:
clear temp
p(i) = zeros(1,n_fns);
s_loc(i) = zeros(1, n_fns);
temp(i) = coord(i);
max_temp = max(temp);
for k = 1:1:n_fns
[p(k), s_loc(k)] = min(temp);
temp(s_loc(k)) = max_temp + 1;
end
clear temp
temp = meas_fn;
for k = 1:1:n_fns
meas_fn(k,1:max_len) = temp(s_loc(k),1:max_len);
end
coord = p;
clear temp

% Now we have coord and meas_fn as they were, except that the function indices
% (the rows of meas_fn) now increase along with the absolute centre pixel
% coordinates of the functions. That is, meas_fn(1,1:max_len) gives the
% measured slit function with lowest centre coordinate, meas_fn(2,1:max_len)
% is the one with next lowest centre coordinate, and so on.

% IGNORE THIS COMMENT:
% Now we have the vector "p" storing the absolute centre coordinates of every
% measured slit function, in ascending order. The vector "s_loc" contains the
% meas_fn indices of the function corresponding to each of these coordinates.

%interp_fns = 2; % # of functions to use in each interpolation

% Do interpolation to get the slit functions for each pixel:
method = 'linear'
for j = 1:1:n_pix
for k = 1:1:max_len
clear s_temp p_temp
for m = 1:1:n_fns
s_temp(m) = meas_fn(m,k);
p_temp(m) = coord(m) - (centre - k);
end
s_fn(j,k) = interp1(p_temp, s_temp, j - (centre - k), method);
end
end

% find the two closest measured functions to pixel j:
%clear temp
%temp = abs(p - j); % get absolute value of distance between pixel j and
% centres of measured slit functions.
%for k = 1:1:interp_fns
% [min_temp(k), p_loc(k)] = min();

% Now we have indices marking the interp_fns measured slit functions closest
% to pixel j (j is an absolute pixel coordinate) stored in s_loc(k).
%if ((j < coord(s_loc(1))) & (j < coord(s_loc(2)))) | ((j > coord(s_loc(1))) & (j > coord(s_loc(2)))) | (min_temp(1) == 0)
% s_fn(j,1:max_len) = meas_fn(s_loc(1),1:max_len);
%else
% interpolate to get the function values
% if (j < coord(s_loc(1)))
% x1 = coord(s_loc(1)) - j;
% x2 = j - coord(s_loc(2));
% else
% x2 = coord(s_loc(2)) - j;
% x1 = j - coord(s_loc(1));
% end
% Do linear interpolation:
% for k = 1:1:max_len
% s_fn(j,k) = (meas_fn(s_loc(1),k)*x2 + meas_fn(s_loc(2),k)*x1)/(x1 + x2);
% end
%end
%end % end of "for j = 1:1:n_pix"

% plot some results:
figure(1);
clf;
hold on;

%Plot 2 measure fns:
%m = coord(1); % coord of centre
%plot([(m - centre + 1):(m - centre + max_len)], s_fn(m,1:max_len), 'r');
%m = coord(2);
%plot([(m - centre + 1):(m - centre + max_len)], s_fn(m,1:max_len), 'b');

%m = abs(coord(2) - coord(1))/2;
%m = m - mod(m,1);
%m = min(coord(1),coord(2)) + m + 1;

% compare a measured fn to an interpolated one:
z1 = wk1read('s_230.wk1');
z2 = z1';
len(n_fns + 1) = length(z2)
k = 1:1:len(n_fns + 1);
z2(k) = z2(k) - noise;
meas_fn(n_fns + 1, k) = z2(k);
cent_z2 = 230;
cent_z2 = find_centre(n_fns + 1, cent_z2, 1, len(n_fns + 1), 0);
coord(n_fns + 1) = 230 + cent_z2
m = coord(n_fns + 1);
plot([(m - cent_z2 + 1):(m - cent_z2 + len(n_fns+1))], z2(1:len(n_fns+1)), 'g');

% Plot interpolated function:
n = m;
plot([(n - centre + 1):(n - centre + max_len)], s_fn(m,1:max_len), 'k');

title('BLACK = interpolated function, GREEN = measured function');

%n = coord(2);
%plot([(n - centre + 1):(n - centre + max_len)], s_fn(m,1:max_len), 'g');
%z1 = wk1read(strcat(dir_path,'s_230.wk1'));