################################################# 
## Created by Robert Scherf and Matthew Weinzierl
#################################################

###################################################################################################
## Note: Following a referee's suggestion, we have transposed the axes relative to the conventional 
## depiction of Lindahl's analysis, putting tax shares on the vertical axis and public good 
##provision on the horizontal, to better fit with modern economists' graphical intuition. 
###################################################################################################

library("ggplot2")
library("data.table")
library(gridExtra)
library(ggpubr)

Graphs <- function(slope_H,intercept_H, slope_L,intercept_L) {
  
  eq1 <- function(x) {x*slope_H + intercept_H}
  eq2 <- function(x) {x*slope_L + intercept_L}
  
  ####################################
  ## Lindahl 
  ####################################
  
  row <- rbind(c(-1*slope_H, 1),
               c(-1*slope_L, 1))
  col <- c(intercept_H, intercept_L)
  
  intersect <- solve(row,col)  
  Lindahl_share <- round(intersect[1],2) 
  G_star <- round(intersect[2],2) 
  totalsurplus <- round(c( 0.5*(min(1,-intercept_H/slope_H)-Lindahl_share)*(G_star+max(0,intercept_H+slope_H)) + 0.5*(Lindahl_share-max(0,-intercept_L/slope_L))*(G_star+max(0,intercept_L))  ),2)
  Lindahl_surplus <- round(c( 0.5*(min(1,-intercept_H/slope_H)-Lindahl_share)*(G_star+max(0,intercept_H+slope_H))/totalsurplus ),2)
  
  dt.LindahlTrapezoid_H <- data.table(group = c(2,2,2,2), LindahlTrapezoid_H.x = c(Lindahl_share,Lindahl_share,min(1,-intercept_H/slope_H),min(1,-intercept_H/slope_H)), 
                                      LindahlTrapezoid_H.y = c(G_star,0,0,max(0,intercept_H+slope_H)))
  dt.LindahlTrapezoid_L <- data.table(group = c(2,2,2,2), LindahlTrapezoid_L.x = c(Lindahl_share,Lindahl_share,max(0,-intercept_L/slope_L),max(0,-intercept_L/slope_L)), 
                                      LindahlTrapezoid_L.y = c(G_star,0,0,max(0,intercept_L)))
  
  LindahlPlot <- ggplot(data.frame(x = c(0, 1)), aes(x = x)) + 
    xlab("H's tax share") + ylab("Public good") + 
    theme(axis.text=element_text(size=14), axis.title=element_text(size=14)) +
    
    coord_cartesian(ylim = c(0,3)) +
    coord_flip() +

    geom_polygon(data = dt.LindahlTrapezoid_H,aes(x=LindahlTrapezoid_H.x,y=LindahlTrapezoid_H.y,group=group), fill="#ffffbf", alpha=.75) + 
    geom_polygon(data = dt.LindahlTrapezoid_L,aes(x=LindahlTrapezoid_L.x,y=LindahlTrapezoid_L.y,group=group), fill="#d7191c", alpha=.5) +
    
    geom_segment(x=0, y=G_star, xend=Lindahl_share, yend=G_star, colour="gray", linetype="solid") + 
    geom_segment(x=Lindahl_share, y=0, xend=Lindahl_share, yend=G_star, colour="gray", linetype="solid") + 
    
    stat_function(fun=eq1, xlim = c(0,min(1,-intercept_H/slope_H))) + 
    stat_function(fun=eq2, xlim = c(max(0,-intercept_L/slope_L),1)) +
    
    annotate("text", x = 0.72, y = 4, size = 6, label = paste("tau[H]^Lindahl == " , Lindahl_share), parse = TRUE) +
    annotate("text", x = 0.675, y = 4, size = 6, label = paste("G* = ", G_star)) +
    annotate("text", x = 0.64, y = 4, size = 6, label = paste("H's surplus share = ", Lindahl_surplus))
    
  ####################################
  ## Brennan  
  ####################################
  
  Brennan_surplus <- round(c( (G_star*(Lindahl_share-.5) + 0.5*(min(1,-intercept_H/slope_H)-Lindahl_share)*(G_star+max(0,intercept_H+slope_H)))/totalsurplus ),2)
  
  dt.BrennanRectangle_H <- data.table(group = c(2,2,2,2), BrennanRectangle_H.x = c(Lindahl_share,Lindahl_share,.5,.5), 
                                      BrennanRectangle_H.y = c(0,G_star,G_star,0))
  dt.BrennanTrapezoid_H <- data.table(group = c(2,2,2,2), BrennanTrapezoid_H.x = c(min(1,intercept_H/(-slope_H)),Lindahl_share,Lindahl_share,min(1,intercept_H/(-slope_H))), 
                                      BrennanTrapezoid_H.y = c(0,0,G_star,max(0,intercept_H+slope_H)))
  
  dt.BrennanTrapezoid_Lpos <- data.table(group = c(2,2,2,2), BrennanTrapezoid_Lpos.x = c(.5,.5, min(.5,max(0,-intercept_L/slope_L)), min(.5,max(0,-intercept_L/slope_L)) ), 
                                         BrennanTrapezoid_Lpos.y = c(0,intercept_L + slope_L*.5,max(0,intercept_L),0))
  dt.BrennanTrapezoid_Lneg <- data.table(group = c(2,2,2,2), BrennanTrapezoid_Lneg.x = c(.5,.5, max(.5, -intercept_L/slope_L) ,Lindahl_share), 
                                         BrennanTrapezoid_Lneg.y = c(G_star, max(0,intercept_L + slope_L*.5), max(0,intercept_L + slope_L*.5) , G_star))
  
  BrennanPlot <- ggplot(data.frame(x = c(0, 1)), aes(x = x)) + xlab("H's tax share") + ylab("Public good") +
    theme(axis.text=element_text(size=14), axis.title=element_text(size=14)) +
    
    coord_cartesian(ylim = c(0,3)) +
    coord_flip() +
    
    geom_polygon(data = dt.BrennanRectangle_H,aes(x=BrennanRectangle_H.x,y=BrennanRectangle_H.y,group=group), fill="#ffffbf", alpha=.75) + 
    geom_polygon(data = dt.BrennanTrapezoid_H,aes(x=BrennanTrapezoid_H.x,y=BrennanTrapezoid_H.y,group=group), fill="#ffffbf", alpha=.75) +
    geom_polygon(data = dt.BrennanTrapezoid_Lpos,aes(x=BrennanTrapezoid_Lpos.x,y=BrennanTrapezoid_Lpos.y,group=group), fill="#d7191c", alpha=.5) +
    
    geom_segment(x=0, y=G_star, xend=Lindahl_share, yend=G_star, colour="gray", linetype="solid") + 
    geom_segment(x=.5, y=0, xend=.5, yend=G_star, colour="gray", linetype="solid") + 
    
    # draw demand for L
    stat_function(fun=eq1, xlim = c(0,min(1,-intercept_H/slope_H))) +
    # draw demand for R 
    stat_function(fun=eq2, xlim = c(max(0,-intercept_L/slope_L),1)) +
    
    annotate("text", x = 0.72, y = 4, size = 6, label = paste("tau[H]^Brennan == " , 0.5), parse = TRUE) + # add text box
    annotate("text", x = 0.64, y = 4, size = 6, label = paste("H's surplus share = ", Brennan_surplus)) +
    
    geom_polygon(data = dt.BrennanTrapezoid_Lneg,aes(x=BrennanTrapezoid_Lneg.x,y=BrennanTrapezoid_Lneg.y,group=group), colour="#d7191c", fill=NA, size=1.1)
  
  ####################################
  ## Moulin plot 1 of 2
  ####################################
  
  moulin_func <- function(g_star){ 0.5*(G_star + max(intercept_L,g_star) - 2*g_star)*( Lindahl_share - max(0,(g_star-intercept_L)/slope_L) ) + 0.5*(G_star + max(intercept_H + slope_H,g_star) - 2*g_star)*( min(1,(g_star-intercept_H)/slope_H) - Lindahl_share ) - g_star}
  moulin_Hist <- uniroot(moulin_func, c(0,G_star),maxiter=100) 
  moulin_delist <- unlist(lapply(moulin_Hist, '[[', 1)) 
  g_star <- round(unname(moulin_delist[1]),2) 
  
  Moulin_share <- round( c( (0.5*(G_star + max(intercept_H + slope_H,g_star) - 2*g_star)*( min(1,(g_star-intercept_H)/slope_H) - Lindahl_share ) + Lindahl_share*(G_star-g_star))/G_star ),2)
  dt.MoulinRectangle <- data.table(group = c(2,2,2,2), MoulinRectangle.x = c(0,1,1,0), 
                                   MoulinRectangle.y = c(0,0,g_star,g_star))
  
  dt.MoulinRectangle_H <- data.table(group = c(2,2,2,2), MoulinRectangle_H.x = c(0,0,Lindahl_share,Lindahl_share), 
                                     MoulinRectangle_H.y = c(g_star,G_star,G_star,g_star))
  dt.MoulinTrapezoid_H <- data.table(group = c(2,2,2,2), MoulinTrapezoid_H.x = c(Lindahl_share,Lindahl_share, min(1,(g_star-intercept_H)/slope_H),min(1,(g_star-intercept_H)/slope_H)), 
                                     MoulinTrapezoid_H.y = c(g_star,G_star,max(g_star,intercept_H+slope_H),g_star))
  dt.MoulinRectangle_L <- data.table(group = c(2,2,2,2), MoulinRectangle_L.x = c(Lindahl_share,Lindahl_share,1,1), 
                                     MoulinRectangle_L.y = c(g_star,G_star,G_star,g_star))
  dt.MoulinTrapezoid_L <- data.table(group = c(2,2,2,2), MoulinTrapezoid_L.x = c(max(0,(g_star-intercept_L)/slope_L),Lindahl_share,Lindahl_share,max(0,(g_star-intercept_L)/slope_L)), 
                                     MoulinTrapezoid_L.y = c(max(intercept_L,g_star),G_star,g_star,g_star))
  
  MoulinPlot <- ggplot(data.frame(x = c(0, 1)), aes(x = x)) + xlab("H's tax share") + ylab("Public good") +
    theme(axis.text=element_text(size=14), axis.title=element_text(size=14)) +
    
    coord_cartesian(ylim = c(0,3)) +
    coord_flip() +
    
    geom_polygon(data = dt.MoulinRectangle,aes(x=MoulinRectangle.x,y=MoulinRectangle.y,group=group), fill="#2b83ba", alpha=.9) + #
    
    geom_polygon(data = dt.MoulinRectangle_H,aes(x=MoulinRectangle_H.x,y=MoulinRectangle_H.y,group=group), fill="#ffffbf", alpha=.75) +
    geom_polygon(data = dt.MoulinTrapezoid_H,aes(x=MoulinTrapezoid_H.x,y=MoulinTrapezoid_H.y,group=group), fill="#ffffbf", alpha=.75) +
    geom_polygon(data = dt.MoulinRectangle_L,aes(x=MoulinRectangle_L.x,y=MoulinRectangle_L.y,group=group), fill="#d7191c", alpha=.5) + 
    geom_polygon(data = dt.MoulinTrapezoid_L,aes(x=MoulinTrapezoid_L.x,y=MoulinTrapezoid_L.y,group=group), fill="#d7191c", alpha=.5) +
    
    geom_segment(x=0, y=G_star, xend=Lindahl_share, yend=G_star, colour="gray", linetype="solid") + 
    geom_segment(x=0, y=g_star, xend=1, yend=g_star, colour="gray", linetype="solid") + 
    
    stat_function(fun=eq1, xlim = c(0,min(1,-intercept_H/slope_H))) + 
    stat_function(fun=eq2, xlim = c(max(0,-intercept_L/slope_L),1)) + 

    annotate("text", x = 0.67, y = 4, size = 6, label = paste("g* = ", g_star)) 
    
  ####################################
  ## Moulin plot 2 of 2  
  ####################################

  Moulin_surplus <- round(c( (G_star*(Lindahl_share-Moulin_share) + 0.5*(min(1,-intercept_H/slope_H)-Lindahl_share)*(G_star+max(0,intercept_H+slope_H)))/totalsurplus ),2)

  dt.Moulin2Rectangle_H <- data.table(group = c(2,2,2,2), Moulin2Rectangle_H.x = c(Lindahl_share,Lindahl_share,Moulin_share,Moulin_share), 
                                      Moulin2Rectangle_H.y = c(0,G_star,G_star,0))
  dt.Moulin2Trapezoid_H <- data.table(group = c(2,2,2,2), Moulin2Trapezoid_H.x = c(min(1,intercept_H/(-slope_H)),Lindahl_share,Lindahl_share,min(1,intercept_H/(-slope_H))), 
                                      Moulin2Trapezoid_H.y = c(0,0,G_star,max(0,intercept_H+slope_H)))
  
  dt.Moulin2Trapezoid_Lpos <- data.table(group = c(2,2,2,2), Moulin2Trapezoid_Lpos.x = c(Moulin_share,Moulin_share, min(Moulin_share,max(0,-intercept_L/slope_L)), min(Moulin_share,max(0,-intercept_L/slope_L)) ), 
                                         Moulin2Trapezoid_Lpos.y = c(0,intercept_L + slope_L*Moulin_share,max(0,intercept_L),0))
  dt.Moulin2Trapezoid_Lneg <- data.table(group = c(2,2,2,2), Moulin2Trapezoid_Lneg.x = c(Moulin_share,Moulin_share, max(Moulin_share, -intercept_L/slope_L) ,Lindahl_share), 
                                         Moulin2Trapezoid_Lneg.y = c(G_star, max(0,intercept_L + slope_L*Moulin_share), max(0,intercept_L + slope_L*Moulin_share) , G_star))
  
  Moulin2Plot <- ggplot(data.frame(x = c(0, 1)), aes(x = x)) + xlab("H's tax share") + ylab("Public good") +
    theme(axis.text=element_text(size=14), axis.title=element_text(size=14)) +
    
    coord_cartesian(ylim = c(0,3)) +
    coord_flip() +
    
    geom_polygon(data = dt.Moulin2Rectangle_H,aes(x=Moulin2Rectangle_H.x,y=Moulin2Rectangle_H.y,group=group), fill="#ffffbf", alpha=.75) + 
    geom_polygon(data = dt.Moulin2Trapezoid_H,aes(x=Moulin2Trapezoid_H.x,y=Moulin2Trapezoid_H.y,group=group), fill="#ffffbf", alpha=.75) + 
    geom_polygon(data = dt.Moulin2Trapezoid_Lpos,aes(x=Moulin2Trapezoid_Lpos.x,y=Moulin2Trapezoid_Lpos.y,group=group), fill="#d7191c", alpha=.5) +
    
    geom_segment(x=0, y=G_star, xend=Lindahl_share, yend=G_star, colour="gray", linetype="solid") + 
    geom_segment(x=Moulin_share, y=0, xend=Moulin_share, yend=G_star, colour="gray", linetype="solid") + 
    
    stat_function(fun=eq1, xlim = c(0,min(1,-intercept_H/slope_H))) + 
    stat_function(fun=eq2, xlim = c(max(0,-intercept_L/slope_L),1)) + 
    
    annotate("text", x = 0.72, y = 4, size = 6, label = paste("tau[H]^Moulin == " , Moulin_share), parse = TRUE) +
    annotate("text", x = 0.64, y = 4, size = 6, label = paste("H's surplus share = ", Moulin_surplus)) +  
    
    geom_polygon(data = dt.Moulin2Trapezoid_Lneg,aes(x=Moulin2Trapezoid_Lneg.x,y=Moulin2Trapezoid_Lneg.y,group=group), colour="#d7191c", fill=NA, size=1.1)
  
  ####################################
  ## Hines plot 1 of 2
  ####################################
  
  hines_func <- function(rho){ rho*G_star - .5*(Lindahl_share - rho)*(intercept_H + rho*slope_H - G_star) + (1-Lindahl_share)*G_star + 0.5*(Lindahl_share - max((-intercept_L/slope_L),1-rho))*(G_star + max(0,intercept_L + (1-rho)*slope_L)) - G_star }
  hines_Hist <- uniroot(hines_func, c(.5,1),maxiter=100) 
  hines_delist <- unlist(lapply(hines_Hist, '[[', 1))
  rho <- round(unname(hines_delist[1]),2) 
  
  Hines_share <- round(c( (rho*G_star - .5*(Lindahl_share - rho)*(intercept_H + rho*slope_H - G_star))/G_star ),2)
  
  dt.HinesRectangle_H <- data.table(group = c(2,2,2,2), HinesRectangle_H.x = c(0,rho,rho,0), 
                                    HinesRectangle_H.y = c(0,0,G_star,G_star))
  
  dt.HinesTriangle_H <- data.table(group = c(1,1,1), HinesTriangle_H.x = c(rho,Lindahl_share,rho), 
                                   HinesTriangle_H.y = c(G_star,G_star,intercept_H+slope_H*rho))

  dt.HinesRectangle_L <- data.table(group = c(2,2,2,2), HinesRectangle_L.x = c(Lindahl_share,1,1,Lindahl_share), 
                                    HinesRectangle_L.y = c(0,0,G_star,G_star))

  dt.HinesTrapezoid_L <- data.table(group = c(2,2,2,2), HinesTrapezoid_L.x = c(Lindahl_share,Lindahl_share, max(1-rho,-intercept_L/slope_L), max(1-rho,-intercept_L/slope_L)), 
                                    HinesTrapezoid_L.y = c(G_star,0,0, max(0,intercept_L + slope_L*(1-rho)) ))
  
  HinesPlot <- ggplot(data.frame(x = c(0, 1)), aes(x = x)) + xlab("H's tax share") + ylab("Public good") +
    theme(axis.text=element_text(size=14), axis.title=element_text(size=14)) +
    
    coord_cartesian(ylim = c(0,3)) +
    coord_flip() +

    geom_polygon(data = dt.HinesRectangle_H,aes(x=HinesRectangle_H.x,y=HinesRectangle_H.y,group=group), fill="#ffffbf", alpha=.75) + 
    geom_polygon(data = dt.HinesRectangle_L,aes(x=HinesRectangle_L.x,y=HinesRectangle_L.y,group=group), fill="#d7191c", alpha=.5) + 
    geom_polygon(data = dt.HinesTrapezoid_L,aes(x=HinesTrapezoid_L.x,y=HinesTrapezoid_L.y,group=group), fill="#d7191c", alpha=.5) + 
    
    geom_segment(x=0, y=G_star, xend=Lindahl_share, yend=G_star, colour="gray", linetype="solid") + 
    geom_segment(x=rho, y=0, xend=rho, yend=G_star, colour="gray", linetype="solid") +
    geom_segment(x=1-rho, y=0, xend=1-rho, yend=G_star, colour="gray", linetype="solid") + 
    
    stat_function(fun=eq1, xlim = c(0,min(1,-intercept_H/slope_H))) + 
    stat_function(fun=eq2, xlim = c(max(0,-intercept_L/slope_L),1)) + 
    
    annotate("text", x = 0.675, y = 4, size = 6, label = paste("rho == " , rho), parse = TRUE) +
      
    geom_polygon(data = dt.HinesTriangle_H,aes(x=HinesTriangle_H.x,y=HinesTriangle_H.y,group=group), colour="#d7191c", linetype="dashed", fill=NA, size=1.1) 

  ####################################
  ## Hines plot 2 of 2
  ####################################
  
  Hines_surplus <- round(c( (G_star*(Lindahl_share-Hines_share) + 0.5*(min(1,-intercept_H/slope_H)-Lindahl_share)*(G_star+max(0,intercept_H+slope_H)))/totalsurplus ),2)
  
  dt.Hines2Rectangle_H <- data.table(group = c(2,2,2,2), Hines2Rectangle_H.x = c(Lindahl_share,Lindahl_share,Hines_share,Hines_share), 
                                     Hines2Rectangle_H.y = c(0,G_star,G_star,0))
  dt.Hines2Trapezoid_H <- data.table(group = c(2,2,2,2), Hines2Trapezoid_H.x = c(min(1,intercept_H/(-slope_H)),Lindahl_share,Lindahl_share,min(1,intercept_H/(-slope_H))), 
                                     Hines2Trapezoid_H.y = c(0,0,G_star,max(0,intercept_H+slope_H)))
  
  dt.Hines2Trapezoid_Lpos <- data.table(group = c(2,2,2,2), Hines2Trapezoid_Lpos.x = c(Hines_share,Hines_share, min(Hines_share,max(0,-intercept_L/slope_L)), min(Hines_share,max(0,-intercept_L/slope_L)) ), 
                                        Hines2Trapezoid_Lpos.y = c(0,intercept_L + slope_L*Hines_share,max(0,intercept_L),0))
  dt.Hines2Trapezoid_Lneg <- data.table(group = c(2,2,2,2), Hines2Trapezoid_Lneg.x = c(Hines_share,Hines_share, max(Hines_share, -intercept_L/slope_L) ,Lindahl_share), 
                                        Hines2Trapezoid_Lneg.y = c(G_star, max(0,intercept_L + slope_L*Hines_share), max(0,intercept_L + slope_L*Hines_share) , G_star))
  
  Hines2Plot <- ggplot(data.frame(x = c(0, 1)), aes(x = x)) + xlab("H's tax share") + ylab("Public good") +
    theme(axis.text=element_text(size=14), axis.title=element_text(size=14)) +
    
    coord_cartesian(ylim = c(0,3)) +
    coord_flip() +
    
    geom_polygon(data = dt.Hines2Rectangle_H,aes(x=Hines2Rectangle_H.x,y=Hines2Rectangle_H.y,group=group), fill="#ffffbf", alpha=.75) + 
    geom_polygon(data = dt.Hines2Trapezoid_H,aes(x=Hines2Trapezoid_H.x,y=Hines2Trapezoid_H.y,group=group), fill="#ffffbf", alpha=.75) + 
    geom_polygon(data = dt.Hines2Trapezoid_Lpos,aes(x=Hines2Trapezoid_Lpos.x,y=Hines2Trapezoid_Lpos.y,group=group), fill="#d7191c", alpha=.5) +
    
    geom_segment(x=0, y=G_star, xend=Lindahl_share, yend=G_star, colour="gray", linetype="solid") + 
    geom_segment(x=Hines_share, y=0, xend=Hines_share, yend=G_star, colour="gray", linetype="solid") +  
    
    stat_function(fun=eq1, xlim = c(0,min(1,-intercept_H/slope_H))) + 
    stat_function(fun=eq2, xlim = c(max(0,-intercept_L/slope_L),1)) + 
    
    annotate("text", x = 0.72, y = 4, size = 6, label = paste("tau[H]^Hines == " , Hines_share), parse = TRUE) +
    annotate("text", x = 0.64, y = 4, size = 6, label = paste("H's surplus share = ", Hines_surplus)) +
    
    geom_polygon(data = dt.Hines2Trapezoid_Lneg,aes(x=Hines2Trapezoid_Lneg.x,y=Hines2Trapezoid_Lneg.y,group=group), colour="#d7191c", fill=NA, size=1.1)
  
  #####################################
  ## Arranging plots and saving as PDF
  #####################################
  
  #combined <- ggarrange(LindahlPlot, BrennanPlot, MoulinPlot, Moulin2Plot, HinesPlot, Hines2Plot + rremove("x.text"), ncol = 2, nrow = 3) 
  
  #ggexport(combined, filename = "graphs.pdf") #export to PDF
  
  #####################################
  ## Saving plots separately as JPG
  #####################################  
  
  ggsave(filename = "BWfig1new.jpg", plot = print(LindahlPlot),
         scale = 1, width = 7.47, height = 7.47, units = c("in"),
         dpi = 300, limitsize = TRUE)
  
  ggsave(filename = "BWfig2new.jpg", plot = print(BrennanPlot),
         scale = 1, width = 7.47, height = 7.47, units = c("in"),
         dpi = 300, limitsize = TRUE)
  
  ggsave(filename = "BWfig3new.jpg", plot = print(MoulinPlot),
         scale = 1, width = 7.47, height = 7.47, units = c("in"),
         dpi = 300, limitsize = TRUE)
  
  ggsave(filename = "BWfig4new.jpg", plot = print(Moulin2Plot),
         scale = 1, width = 7.47, height = 7.47, units = c("in"),
         dpi = 300, limitsize = TRUE)
  
  ggsave(filename = "BWfig5new.jpg", plot = print(HinesPlot),
         scale = 1, width = 7.47, height = 7.47, units = c("in"),
         dpi = 300, limitsize = TRUE)
  
  ggsave(filename = "BWfig6new.jpg", plot = print(Hines2Plot),
         scale = 1, width = 7.47, height = 7.47, units = c("in"),
         dpi = 300, limitsize = TRUE)
  
}

#############################################################################################################################
## Final Output: Graph Functions
## In order to create the graphs for different demand equations, enter new parameter values in the Graphs function below
## Please note that for the code to run smoothly, individual H should have higher demand for the public good than individual L
##############################################################################################################################

##Graphs(slope_H,intercept_H, slope_L,intercept_L) 
Graphs(-5,5,6,-2)

#############
# End of file
#############