function group_sel(numturns)
 % This is an attempt to make a crude model of two bird populations that will
 % exhibit group selection.  See David Sloan Wilson, _Darwin's Cathedral_.
 % Last modified 2-13-2021 by Peter Taylor.
 page_screen_output(0);
 %
 % This is a turn-based emigration model rather than a wipeout model.  Replacement of
 % the less successful flock membership is gradual rather than catastrophic.
 % Two flocks, 1 and 2.  Two alleles of a gene, "altruistic" (A) and "selfish" (B).
 %
 % Initial  populations:
 pop = [ 10 90;  0 100];   % = [ A1 S1; A2 S2 ]
 %
 % The probability of reproduction is constant.  I only care about one gene with two
 % alleles.  Male vs. female doesn't matter.
 % I am assuming that food is available, but a large flock requires birds to spend more time
 % and travel longer distances looking for food, which increases deaths from predators.  That
 % may not be consistent with flock size being a meaningful "trait-group" for purposes of an
 % "altruistic" bird being able to warn other members of the group.
 Pr = 0.3;
 %
 % Deaths occur at rate Pd(k,n) per turn, which is a function of the flock population (fp(n)),
 % the individual's allele (k), and proportion of altruistic alleles within the flock (ar).
 % Pd = c3*fp/(c4+fp)*(1+c1*(k==1))/(1+ar), where fp(n) = sum(pop(n,:)) and ar = pop(n,1)/fp.
 c1 = 0.01;     % How much does altruism increase an individual bird's death rate?
 c2 = 0.9;      % How much do a lot of altruists reduce the flock's death rate?
 c3 = 3.0*Pr;   % What is the maximum death rate for a large, selfish flock?
 % For equilibrium population fpeq = 100 (at zero altruism condition), I need (1+Pr)*(1-Pd) = 1.
 fpeq = 100;
 c4 = c3*fpeq*(1+Pr)/Pr-fpeq;   % Parameter to get a reasonable equilibrium population.
 %
 % The probability of a bird with allele k moving from flock n to flock m is
 % Pm = max(0,fp(m)-fp(n))*0.5/fp(n)*c5*(1-c6*(i1==1))
 c5 = 0.5;   % How rapidly does the probability of migration rise with population difference?
 c6 = 0.9;   % How much does the "altruistic" allele inhibit migration?
 %
 % Each turn has three phases: (1) death, (2) reproduction, and (3) migration.  My problem
 % now is to avoid introducing bias into the simulation by doing something dumb with the
 % order in which I process the different alleles.  Calculate the four Pd and Pm values at the
 % beginning of the turn and hold them constant for that turn.  If there is an instability
 % problem with the steps being too big, either reduce the rates of change or put in standard 
 % differential equation solution logic (e.g. Runge-Kutta).
 %
 % Current model (2-13-2021):
 % The results are surprisingly sensitive to chance.
 % I want to test the sensitivity to the following: c1 and c6.
 % Baseline: c1 = 0.01;   % How much does altruism increase an individual bird's death rate?
 %           c6 = 0.9;    % How much does the "altruistic" allele inhibit migration?
 % After about 1000 turns, I typically end up with something like [295, 0; 240, 35].
 %
 % Results:
 %  c1     c5    c6   typical result
 % 0.01   0.5   0.9   [295, 0; 240, 35]
 % 0.03   0.5   0.9   Erratic.  Altruists often go extinct.
 % 0.05   0.5   0.9   Altruists usually go extinct, but not always.
 % 0.01   0.5   0.7   [295, 0; 260, 20]  Altruists sometimes go extinct.  Groups sometimes switch places.
 % 0.01   0.5   0.5   [240, 20; 190, 60]  Selfish survive in both groups.  Altruists often go extinct.
 % 0.01   0.5   0.3   Erratic.  Altruists usually to go extinct.
 %
 % Try raising c1:
 % At about 0.03, I see inconsistent results.  Sometime the altruists go extinct,
 % but sometimes they still take over population 1.  At c1 = 0.05, the altruists usually 
 % go extinct.
 % I have eliminated c2 from the simulation.  Hiding always works if the predator doesn't gain
 % surprise, and squawking always ends the predator's attack.
 % Try lowering c6:
 % At c6 = 0.7, the two flocks sometimes switch places, and the altruists sometimes go extinct.
 % At c6 = 0.5, the selfish survive in both groups.  Altruists often go extinct.
 % At c6 = 0.3, the two flocks sometimes switch places, and the selfish survive in both flocks.
 %
 Pd = zeros(2,2);
 Pm = zeros(2,2);
 %
 disp('');
 disp('------------------------------     Begin simulation     ------------------------------');
 disp('');
 turn = 0
 pop
 for turn = 1:numturns
  disp('');
  turn
  %
  % Phase 1.  Calculate deaths.
  for flock = 1:2
   for allele = 1:2
    fp = sum(pop(flock,:));   % Flock population
    ar = pop(flock,1)/fp;     % Altruism ratio within flock
    altruistic = mod(allele,2);
    Pd(flock,allele) = c3*fp/(c4+fp) * (1+c1*altruistic)/(1+ar);
   end
  end
  % I want to see some occasional unlikely results with small numbers.
  %
  newpop = pop-pop.*Pd;
  excess = newpop-floor(newpop);
  newpop = floor(newpop)+(rand(2)<excess);
  death = pop-newpop;
  pop = newpop;
  %
  % Phase 2.  Calculate reproduction.
  newpop = pop+pop.*Pr;
  excess = newpop-floor(newpop);
  newpop = floor(newpop)+(rand(2)<excess);
  hatch = newpop-pop;
  pop = newpop;
  %
  % Phase 3.  Calculate migration.
  for flock = 1:2
   for allele = 1:2
    crowding = sum(pop(flock,:))/(sum(pop(flock,:))+sum(pop(3-flock,:)));
    altruistic = mod(allele,2);
    % Pm(flock,allele) = c5*0.5*(1-cos(pi*crowding))*(1-c6*altruistic);
    Pm(flock,allele) = c5*0.5*max(0,sum(pop(flock,:))/sum(pop(3-flock,:))-1)*(1-c6*altruistic);
   end
  end
  movers = pop.*Pm;
  excess = movers-floor(movers);
  movers = floor(movers)+(rand(2)<excess);
  newpop = pop-movers+[0 1;1 0]*movers;
  pop = newpop
 end
 % pop
end

