#!/usr/bin/python3
#exec(open("/home/sasa/Dropbox/MyPy/dwl.py").read())
from sympy import *

#Problema del monopolista: max q*D(q)-c(q) dove D è il prezzo di domanda

#quantità q
q=Symbol('q')
#funzione costo quadratica paramtetri a,b,c 
[a,b,c]=[1,4,1]
def C(q):
	return a+b*q+c*q**2
#costo marginale
mc=diff(C(q),q)

#domanda lineare parametri alpha beta
[alpha,beta]=[40,2]
def D(q):
	return alpha-beta*q
#ricavo marginale
mr=diff(q*D(q),q)

#FOC: costo marginale = ricavo marginale => quantità prodotta qmon
qmon=solve(mr-mc,q)
#prezzo 
pmon=D(qmon[0])
#costo marginale nel punto di ottimo
mcmon=mc.subs(q,qmon[0])
#quantità efficiente data da prezzo = costo marginale
qeff=solve(D(q)-mc,q)

#perdita secca = area fra domanda e costo marginale da qmon a qeff
dwl=(pmon-mcmon)*(qeff[0]-qmon[0])/2
#print(dwl)

#plot everything
from numpy import *
import math
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon

fig = plt.figure()
t = linspace(0, 20,400)
plt.plot(t, D(t), 'r')
MR=lambdify(q, mr, modules=['numpy'])
plt.plot(t, MR(t), 'b')
MC=lambdify(q, mc, modules=['numpy'])
plt.plot(t, MC(t), 'g') # plotting t, c separately 
plt.xlim([0, 25])
plt.ylim([0, 50])
T=plt.Polygon([[qmon[0],mcmon], [qmon[0],pmon], [qeff[0],D(qeff[0])]],closed=True,color='y')
plt.gca().add_patch(T)
#plt.xticks([])
#plt.yticks([])
plt.title('Monopolio: perdita secca nel triangolo verde',fontsize=18)
#plt.xlabel('blah')
#plt.ylabel('blah')
#add Latex text
matplotlib.rc('text', usetex=True)
plt.text(15, 31, r'\textit{Costo Marginale}', {'color': 'g', 'fontsize': 16})
plt.text(14.5, 12, r'\textit{Domanda} $D(q)$', {'color': 'r', 'fontsize': 16})
plt.text(9.5, 3.5, r'\textit{Ricavo Marginale}', {'color': 'b', 'fontsize': 16})
# save figure if you want
#fig.savefig('test.pdf')
plt.show()