function [cell_els,cells_out] = LEM_cell_els(cells,pt1,pt2)
% LEM_cell_els - find cell elements traversed from pt1 to pt2
% On input:
%     cells (grid struct): defines cells
%       .x_min (float): x min for cells
%       .y_min (float): y min for cells
%       .x_max (float): x max for cells
%       .y_max (float): y max for cells
%       .del_x (float): x axis spacing
%       .del_y (float): y axis spacing
%       .M (int): number of rows
%       .N (int): number of cols
%       .launch_sites (kx3 array): x,y,z launch locations
%       .land_sites (kx3 array): x,y,z land locations
%       .hd (float): headway distance
%       .cell_lists (array struct): 
%          (ind1,ind2).list (nx1 vector): [ind1,...,ind2] linear indexes
%     pt1 (1x3 vector): 3D point
%     pt2 (1x3 vector): 3D point
%     step (float): dist to move along segment pt1-pt2
% On output:
%     cell_els (nx1 vector): cell elements (linearized by column)
% Call:
%     flights(f).cell_els = LEM_cell_els(grid,traj(1,1:3),traj(end,4:6));
% Author:
%     T. Henderson
%     UU
%     Spring 2020
%

cell_els = [];
cells_out = cells;

x_min = cells.x_min;
x_max = cells.x_max;
x_spread = x_max - x_min;
y_min = cells.y_min;
y_max = cells.y_max;
y_spread = y_max - y_min;
del_x = cells.del_x;
del_y = cells.del_y;
M = cells.M;
N = cells.N;

% row col indexes

c1 = ceil(N*(pt1(1)-x_min)/x_spread);
r1 = M - ceil(M*(pt1(2)-y_min)/y_spread) + 1;
c2 = ceil(N*(pt2(1)-x_min)/x_spread);
r2 = M - ceil(M*(pt2(2)-y_min)/y_spread) + 1;
ind1 = min(max(1,M*(c1-1) + r1),M*N);
ind2 = min(max(1,M*(c2-1) + r2),M*N);
if isempty(cells.cell_lists(ind1,ind2).list)
    [rows,cols] = LEM_line_between([r1,c1],[r2,c2]);
    clist = zeros(length(rows),1);
    for k = 1:length(rows)
        clist(k) = M*(cols(k)-1)+rows(k);
    end
    cells_out.cell_lists(ind1,ind2).list = clist;
    cell_els = clist;
else
    cell_els = cells.cell_lists(ind1,ind2).list;
end
