function [res, res1, Res, Res1, err, err_ft, pixel_shifts, pixel_shifts1, Raw, Raw1] = find_N1_N2(raw, raw1, N11, N2, Ns, band_num, oversampling_factor)
iterations = 20;
points_in_TR1 = Ns + N11 + N2;
dt1 = 1/points_in_TR1;
% norm_time1 = linspace(N11*dt1, (points_in_TR1 - 1 - N2)*dt1, Ns);
norm_time1 = linspace((N11 + 1)*dt1, (points_in_TR1 - N2)*dt1, Ns);

% correct for phase shifts between columns of raw.
% sub_pixel_phase = norm_time vertically, multiplied by [-pi to pi] horizontally.
% correction done with the assumed norm_time1, which may not be equal to the
% true norm_time which is not known
% ========================================================================
Nband_vec = (0:band_num - 1) - floor(band_num/2);
sub_pixel_phase = 2*pi*norm_time1(:).*Nband_vec/band_num;
sub_pixel_phase1 = 2*pi*norm_time1(:).*(Nband_vec + 0.5)/band_num;

Raw = raw.*exp(1j*sub_pixel_phase);       % phase correct sampled data.
Raw1 = raw1.*exp(1j*sub_pixel_phase1);       % phase correct sampled data.
dataout = zeros(points_in_TR1, size(raw, 2));
dataout1 = dataout;
res = dataout;
res1 = dataout;
pixel_shifts = dataout;
pixel_shifts1 = dataout;
Dval = dataout;
Dval1 = dataout;
residue = zeros(iterations, size(raw, 2));
residue1 = zeros(iterations, size(raw, 2));

k1 = 2;
k2 = 15;
% k1 = 2;
% k2 = 5;

% correct res, res1 by: a) sparse completion of points to TR. b) ring
% removal by the method of Kellner et. al.
for j = 1:size(raw, 2)
    [dataout(:, j), residue(:, j)] = find_sparse(Raw(:, j), N11, N2, iterations);
    [dataout1(:, j), residue1(:, j)] = find_sparse(Raw1(:, j), N11, N2, iterations);

    [~, res(:, j), pixel_shifts(:, j), Dval(:, j)] = ring_remove(dataout(:, j), k1, k2, 1/8);
    %     [~, res(:, j)] = ring_remove(fftc(res(:, j)), k1, k2, 1/8);
    [~, res1(:, j), pixel_shifts1(:, j), Dval1(:, j)] = ring_remove(dataout1(:, j), k1, k2, 1/8);
    %     [~, res1(:, j)] = ring_remove(fftc(res1(:, j)), k1, k2, 1/8);
end
Res = ifftc(dataout);
Res1 = ifftc(dataout1);

% find error from the sparse number of points.
% Choose the best column in "raw". W weight more for high frequencies. Use
% the strongest high frequecy column. The "dc" term takes care of
% lines with only a dc component.
dc = 0.15;
[M, N] = size(res);
W = fftshift(fermi_filt(0.07, 0.1, round(M/oversampling_factor)));
W = (1 - dc)*W + dc;
Tmp = W(:).*crop(res, round(M/oversampling_factor), N);
[~, K] = max(sum(abs(Tmp)));
Tmp = W(:).*crop(res1, round(M/oversampling_factor), N);
[~, K1] = max(sum(abs(Tmp)));

% err = mean(crop(Dval(:, K), round(M/oversampling_factor))) + mean(crop(Dval1(:, K1), round(M/oversampling_factor)));
err = max(crop(Dval(:, K), round(M/oversampling_factor))) + max(crop(Dval1(:, K1), round(M/oversampling_factor)));

err_ft = find_error(Res(:, K), oversampling_factor) + find_error(Res1(:, K1), oversampling_factor);

function Error = find_error(datain, oversampling_factor)
datain = crop(abs(datain(:)), round(length(datain(:))/oversampling_factor));
Error = mean(datain/max(datain));
