

Lat = 5;
Lon = 5;
NumRows = 20038;

%% Create Grid

BaseBin(1) = 1;
for Row = 1:NumRows
   LatBin(Row) = ((Row-.5)*180/NumRows) - 90; 
   NumBin(Row) = floor((2*NumRows*cosd(LatBin(Row)))+.5); 
   if Row>1
       BaseBin(Row) = BaseBin(Row-1)+NumBin(Row-1);
   end
end
TotBins = BaseBin(NumRows)+NumBin(NumRows)-1;

%% Find bin IDX from Lat Lon

Row = floor((90+Lat)*(NumRows/180))+1;
Row = min(Row,NumRows);
LonP = Lon+180;
Col = floor(LonP.*NumBin(Row)'/360)+1;
Col = min(Col,NumBin(Row)');
Idx = BaseBin(Row)' + Col-1;

%% Find bin IDX center Lat Lon

Row = NumRows;
while Idx<BaseBin(Row)
    Row = Row-1;
end
CenLat = LatBin(Row);
CenLon = 360*((Idx-BaseBin(Row)+.5)/NumBin(Row));
CenLon = CenLon - 180;

%% Get bin IDX corner Lat Lon

LatBinDeg = (180/NumRows)/2;
LonBinDeg = (360/NumBin(Row))/2;

% bottom left, top left, top right, bottom right
LonCor = [CenLon-LonBinDeg CenLon-LonBinDeg CenLon+LonBinDeg CenLon+LonBinDeg];
LatCor = [CenLat-LatBinDeg CenLat+LatBinDeg CenLat+LatBinDeg CenLat-LatBinDeg];

%% plot bin
plot(Lon,Lat,'+k')
hold on
plot([LonCor LonCor(1)],[LatCor LatCor(1)],'r')