function tennissimulation
   np = 50	% number of point-winning probabilities tested
   ngames = 30000
   p1_win_fraction = [];
   p1point = linspace(0,1,np); % probability of player 1 winning any given point
   for j = 1:np
      gws = zeros(1,ngames);
      for i=1:ngames 
         gw(i) = game_winner(p1point(j));
      end
      gw;
      p1wins = (gw == 1);
      p1_win_fraction(j) = norm(double(p1wins), 1)/ngames;
	  fprintf('%i \n',j);
   end
   fprintf('(total games played: %i)\n',np*ngames);
   plot(p1point,p1_win_fraction,'-or');
   axis equal
   xlabel('probability player 1 wins a point')
   ylabel('probability player 1 wins a game')
   title(sprintf('%i trials at each probability',ngames));
   axis([0 1 0 1])
end
 
function pw = point_winner(player_1_probability_of_winning_point)
   if( rand(1,1) < player_1_probability_of_winning_point )
     pw = 1;
   else
     pw = 2;
   end 
 end

function gw = game_winner(player_1_probability_of_winning_point)
   s = [0,0];
   while ( norm(s,inf) < 3 )
      pw = point_winner(player_1_probability_of_winning_point);
      s( pw ) = s( pw ) + 1;
   end

   while ( s ~= [3,3] )
      pw = point_winner(player_1_probability_of_winning_point);
      if ( s( pw ) == 3 )
         gw = pw;
         s;
         return;
      else
         s( pw ) = s( pw ) + 1;
      end
   end

   % "deuce"
   while ( 1 )
%      fprintf('Deuce\n');
      pw1 = point_winner(player_1_probability_of_winning_point); 
      pw2 = point_winner(player_1_probability_of_winning_point); 
      if ( pw2 == pw1 ) 
         gw = pw2;
         s( pw2 ) = s( pw2 ) + 1;
         return;
      end
   end
end