% This version is for the first half of the game.  The boundary condition is that
% a specified team will kick off to start the second half.  One must run dynprog.m
% with T equal to the full second half before running this program.

% R(t,d,y,g) = prob that we win if in regular (scrimmage) mode
% K(t,d,g) = prob that we win if in kickoff mode
% E(t,d,g) = prob that we win if in extra-point mode
% t is time period, with t=1 the final period
% d indexes point differential; maxd is largest differential, and d runs 1:2*maxd+1
% y is yardage marker; there are numy of them
% g indicates possession; g=1 if we have the ball, g=2 if they have the ball
% if onside kick, success with probability ponside, and either way ball is at our 40.

% Copyright 2004 by William S. Krasker

% Set gkickoff=1 if we will kick off to start the second half; otherwise gkickoff=2.
gkickoff=1;

T=181;

% Initialize arrays

R=zeros(T,dtop,numy,numg);
K=zeros(T,dtop,numg);
E=zeros(T,dtop,numg);
Oneortwo=zeros(T,dtop,numg);

% Boundary conditions

for g=1:2
   for d=1:dtop
      K(1,d,g)=B(d,gkickoff);
      for y=1:numy
         R(1,d,y,g)=B(d,gkickoff);
      end
   end
end

for d=1:dtop
   gofor1=p1pt*B(min(d+1,dtop),gkickoff)+(1-p1pt)*B(d,gkickoff);
   gofor2=p2pt*B(min(d+2,dtop),gkickoff)+(1-p2pt)*B(d,gkickoff);
   E(1,d,1)=max(gofor1,gofor2);
   if gofor2>gofor1
      Oneortwo(1,d,1)=2;
   else
      Oneortwo(1,d,1)=1;
   end
   gofor1=p1pt*B(max(d-1,1),gkickoff)+(1-p1pt)*B(d,gkickoff);
   gofor2=p2pt*B(max(d-2,1),gkickoff)+(1-p2pt)*B(d,gkickoff);
   E(1,d,2)=min(gofor1,gofor2);
   if gofor2<gofor1
      Oneortwo(1,d,2)=2;
   else
      Oneortwo(1,d,2)=1;
   end
end

% backward induction

for t=2:T
  
   for d=1:dtop
      kickaway=ptranskickoff(numy+1)*E(t-1,max(d-6,1),2)+ptranskickoff(2*numy+2)*E(t-1,min(d+6,dtop),1);
      for y2=1:numy
         kickaway=kickaway+ptranskickoff(y2)*R(t-1,d,y2,2)+ptranskickoff(y2+numy+1)*R(t-1,d,y2,1);
      end
      kickonside=ponside*R(t-1,d,3,1)+(1-ponside)*R(t-1,d,4,2);
      K(t,d,1)=max(kickaway,kickonside);
   end

   for d=1:dtop
      gofor1=p1pt*K(t,min(d+1,dtop),1)+(1-p1pt)*K(t,d,1);
      gofor2=p2pt*K(t,min(d+2,dtop),1)+(1-p2pt)*K(t,d,1);
      E(t,d,1)=max(gofor1,gofor2);
      if gofor2>gofor1
         Oneortwo(t,d,1)=2;
      else
         Oneortwo(t,d,1)=1;
      end
   end
  
   for d=1:dtop
      for y=1:numy
                  
         rvalue1=0;
         for k=1:min(t-1,tdrivemax)
            rvalue1=rvalue1+ptrans(y,numy+1)*tprobs(y,numy+1,k)*E(t-k,min(d+6,dtop),1) ...
               +ptrans(y,numy+2)*tprobs(y,numy+2,k)*K(t-k,min(d+3,dtop),1) ...
               +ptrans(y,numy+3)*tprobs(y,numy+3,k)*E(t-k,max(d-6,1),2);
            for y2=1:numy
               rvalue1=rvalue1+ptrans(y,y2)*tprobs(y,y2,k)*R(t-k,d,y2,2);
            end
         end
         for k=min(t-1,tdrivemax)+1:tdrivemax
            rvalue1=rvalue1+ptrans(y,numy+1)*tprobs(y,numy+1,k)*R(1,d,y,2) ...
               +ptrans(y,numy+2)*tprobs(y,numy+2,k)*R(1,d,y,2) ...
               +ptrans(y,numy+3)*tprobs(y,numy+3,k)*R(1,d,y,2);
            for y2=1:numy
               rvalue1=rvalue1+ptrans(y,y2)*tprobs(y,y2,k)*R(1,d,y2,2);
            end
         end
         R(t,d,y,1)=rvalue1;
         
      end
   end  
   
   for d=1:dtop
      kickaway=ptranskickoff(numy+1)*E(t-1,min(d+6,dtop),1)+ptranskickoff(2*numy+2)*E(t-1,max(d-6,1),2);
      for y2=1:numy
         kickaway=kickaway+ptranskickoff(y2)*R(t-1,d,y2,1)+ptranskickoff(y2+numy+1)*R(t-1,d,y2,2);
      end
      kickonside=ponside*R(t-1,d,3,2)+(1-ponside)*R(t-1,d,4,1);
      K(t,d,2)=min(kickaway,kickonside);
   end

   for d=1:dtop
      gofor1=p1pt*K(t,max(d-1,1),2)+(1-p1pt)*K(t,d,2);
      gofor2=p2pt*K(t,max(d-2,1),2)+(1-p2pt)*K(t,d,2);
      E(t,d,2)=min(gofor1,gofor2);
      if gofor2<gofor1
         Oneortwo(t,d,2)=2;
      else
         Oneortwo(t,d,2)=1;
      end
   end
  
   for d=1:dtop
      for y=1:numy
                  
         rvalue1=0;
         for k=1:min(t-1,tdrivemax)
            rvalue1=rvalue1+ptrans(y,numy+1)*tprobs(y,numy+1,k)*E(t-k,max(d-6,1),2) ...
               +ptrans(y,numy+2)*tprobs(y,numy+2,k)*K(t-k,max(d-3,1),2) ...
               +ptrans(y,numy+3)*tprobs(y,numy+3,k)*E(t-k,min(d+6,dtop),1);
            for y2=1:numy
               rvalue1=rvalue1+ptrans(y,y2)*tprobs(y,y2,k)*R(t-k,d,y2,1);
            end
         end
         for k=min(t-1,tdrivemax)+1:tdrivemax
            rvalue1=rvalue1+ptrans(y,numy+1)*tprobs(y,numy+1,k)*R(1,d,y,1) ...
               +ptrans(y,numy+2)*tprobs(y,numy+2,k)*R(1,d,y,1) ...
               +ptrans(y,numy+3)*tprobs(y,numy+3,k)*R(1,d,y,1);
            for y2=1:numy
               rvalue1=rvalue1+ptrans(y,y2)*tprobs(y,y2,k)*R(1,d,y2,1);
            end
         end
         R(t,d,y,2)=rvalue1;
         
      end
   end   
   
end