#simple model of farming/gathering transition in the southwest


#--------------------------------------------------------------------
# note - time has been rescaled in terms of the decay rate of seeds. 
# thus rs = rs/ds, rp=rp/ds, etc.                                   
#--------------------------------------------------------------------

par rs=1,rp=2,dc=0,A=1,As=1,gs=1,b=0.5,sw=1
par Kmin=0.25,Ac=1,al_c=1,N=0.01

#----------------------------------------------------------
# Define clearing function

f(x,a,b) = (x**b)/(a**b + x**b)

#----------------------------------------------------------
# Compute optimal, state-dependent labor allocation   
#                                                     
#   Recall, we are solving
#   min La+Ls
#   subject to
#   AcLc + As(s)Ls^b >= Kmin
#   Lc + Ls <= L
#
#  note, I am assuming that L=1.  
#
#----------------------------------------------------------

#----------------------------------------------------------
# Compute optimal labor allocation to gathering
#
# That is - if the demand is high enough, and labor is sufficient, this
#
# is how much labor should be devoted to seed gathering

Lso = (b*As*(s**gs)/(Ac*p))**(1/(1-b))
Hso = As*(s**gs)*Lso**b

#-----------------------------------------------------------

#----------------------------------------------------------
# Apply constraints: if demand is above the maximum producible, or 
#the optimal is above demand, then compute the "desired" harvest and labor. 

Hsmax = As*(s**gs)
Hsd = min(min(Hso,Kmin),Hsmax)
Lsd = (Hsd/(As*(s**gs)))**(1/b)

# Given the optimal gathering level, compute the level of agriculture
# necessary to compliment it - i.e. the optimal labor allocation to 
#agriculture

Hc = max(Kmin-Hsd,0)
Lco = max(min(Hc/(Ac*p),1-Lsd),0)
Hco = Ac*p*Lco

# Compute cleared area required for ag output and actual. Derate labor in 
# agriculture accordingly.

Aao= al_c*Lco
Aaa = (A-p)/N
drate=(Aaa+0.001)/(Aao+0.001)
Lca = min(Lco*drate,Lco)
Hca = Ac*p*Lca
Lsa = min(((Kmin-Hca)/(As*(s**gs)))**(1/b),1-Lca)
Hsa = As*(s**gs)*Lsa**b

s' = rs*p - s - Hsa*N
p'= rp*s*(1-p/A) - dc*f(1/drate,1,10)*(1+ p*sw-sw)*N


init s=0.5,p=0.5

@ xp=p,yp=s
@ xlo=0,ylo=0,xhi=1.1*g*a,yhi=1.1*g*a*rs
@ total=10,dt=0.02

aux A_a_o=Aao
aux A_a_a=Aaa
aux H_s_o = Hso
#aux H_smax = Hsmax
aux H_s_d = Hsd
aux L_s_o = Lso
aux L_s_d = Lsd
aux H_s_a = Hsa
aux H_c_o = Hco
aux H_c_a = Hca
aux L_c_o= Lco
aux L_c_a = Lca
aux L_s_a = Lsa
done