code: def lower_brkpts(a,b,c,e,f,p_min,p_max,n):
      
      EPS = 1.e-12        # for avoiding round-off errors
      if f == 0: f = math.pi/(p_max-p_min)
      brk = []
      nvalve = int(math.ceil(f*(p_max-p_min)/math.pi))
      for i in range(nvalve+1):
          p0 = p_min + i*math.pi/f
          if p0 >= p_max-EPS:
              brk.append((p_max,cost(a,b,c,e,f,p_min,p_max)))
              break
          for j in range(n):
              p = p0 + j*math.pi/f/n
              if p >= p_max:
                  break
              brk.append((p,cost(a,b,c,e,f,p_min,p)))
      return brk
docs: lower_brkpts: lower approximation of the cost function
      Parameters:
          - a,...,p_max: cost parameters
          - n: number of breakpoints' intervals to insert between valve points
      Returns: list of breakpoints in the form [(x0,y0),...,(xK,yK)]