rand('state',0);
randn('state',0);
numA=1000;
output=zeros(30,100);

zangle=randn(numA,1);
zdistance=randn(numA,1);
zbouncedistance=randn(numA,1);
udowned=rand(numA,1);

anglechange=zeros(numA,1);

% Ball kicked from left hash mark (3.0833 yards left of center), and 49 yards from own
% goal line. This corresponds approximately to a line of scrimmage at own 61 yard line,
% i.e. opponent's 39 yard line.
x0=-3.0833;
y0=49;

% opponents will fair catch outside their own 5 yard line, otherwise let the ball bounce.
faircatchlim=95;

kickthetasd=(pi/2)*0.11;
kickdistancesd=0.15;

bouncethetasd=(pi/2)*0.8;
bouncedistancemean=7;
bouncedistancesd=0.5;
bouncedistancehaircut=0.33;

for k=1:numA
   for m=1:1000
      anglechange0=randn*bouncethetasd;
      if abs(anglechange0) < pi
         anglechange(k)=anglechange0;
         break
      end
   end
end

bestsofar=0;

for i=1:30
   for j=1:100

%------------ start of simulation for given aim point -----------------
A=zeros(numA,1);

thetaaim=(pi/2)*(50-j)/100;
distanceaim=40+i;

for k=1:numA

spot=0;

theta=thetaaim+zangle(k)*kickthetasd;
while theta > pi
   theta=theta-2*pi;
end
while theta <= -pi
   theta=theta+2*pi;
end
distance=distanceaim*exp(kickdistancesd*zdistance(k)-0.5*kickdistancesd^2);

x=x0+distance*sin(theta);
y=y0+distance*cos(theta);

thetabounce=theta+anglechange(k);
while thetabounce > pi
   thetabounce=thetabounce-2*pi;
end
while thetabounce <= -pi
   thetabounce=thetabounce+2*pi;
end
b0=bouncedistancemean*cos(bouncedistancehaircut*abs(anglechange(k)));
bouncedistance=b0*exp(bouncedistancesd*zbouncedistance(k)-0.5*bouncedistancesd^2);

xbounce=x+bouncedistance*sin(thetabounce);
ybounce=y+bouncedistance*cos(thetabounce);

if theta > atan((26.67-x0)/(100-y0)) & x > 26.67
   spot=y0+(26.67-x0)/tan(theta);
elseif theta < atan((-26.67-x0)/(100-y0)) & x < -26.67
   spot=y0+(-26.67-x0)/tan(theta);
elseif y >= 100
   spot=80;
else
   if y < faircatchlim
      spot=y;
   elseif thetabounce > atan((26.67-x)/(100-y)) & xbounce > 26.67
      spot=y+(26.67-x)/tan(thetabounce);
   elseif thetabounce < atan((-26.67-x)/(100-y)) & xbounce < -26.67
      spot=y+(-26.67-x)/tan(thetabounce);
   elseif ybounce >= 100
      if udowned(k) < (100-y)/(ybounce-y)
         spot=max(98,y);
      else
         spot=80;
      end
   else
      spot=ybounce;
   end
end

A(k)=spot;

end
%------------ end of individual trial -----------------

      if mean(A) > bestsofar
         bestsofar=mean(A);
         istar=i;
         jstar=j;
      end

      output(i,j)=mean(A);
   end
end

bestsofar
xstar=x0+(40+istar)*sin((pi/2)*(50-jstar)/100)
ystar=y0+(40+istar)*cos((pi/2)*(50-jstar)/100)