function  [boat gen] = MainSexualProduction(boat, BOAT_NUM, ROWER, Pc, Pm, radius)

% Get the maximum number of germans 
gnum = max(sum(boat, 2));
% Loop until all-German is available
gen = 0;   % number of generation
while gnum<ROWER
      % pick a location L (from 1 to BOAT_NUM)
    L = round(rand* (BOAT_NUM-1)) + 1;
    % Find the range with respect to L in the wrap around boat
    range = FindRange(L, BOAT_NUM, radius);
    % Randomly select two individual indices A and B whose distances away from L
    % are within the radius 
    A = range(round(rand*2*radius)+1);
    B = range(round(rand*2*radius)+1);
    % Select a parent according to fitness
    fA = ComputeFitness(boat(A,:));
    fB = ComputeFitness(boat(B,:));
    % Select parent 1
    if   fA > fB
        P1 = boat(A,:);
    else
        P1 = boat(B,:);
    end
    % Randomly select another two individual indices A and B whose distances away from L
    % are within the radius 
    A = range(round(rand*2*radius)+1);
    B = range(round(rand*2*radius)+1);
    % Select a parent according to fitness
    fA = ComputeFitness(boat(A,:));
    fB = ComputeFitness(boat(B,:));
    % Select parent 2
    if   fA > fB
        P2 = boat(A,:);
    else
        P2 = boat(B,:);
    end
    % Pick another individual R nearby to L
    R = range(round(rand*2*radius)+1);
    % Produce a child and replace R with it
    boat(R,:) = SexuallyProduce(P1, P2, Pc, Pm);
    % Get the maximum number of germans 
    gnum = max(sum(boat, 2));
    gen = gen + 1;
end