#!/usr/bin/env python
# -*- coding: utf-8 -*-

# ***************************************************************************
# *   Copyright (C) 2013, Paul Lutus                                        *
# *                                                                         *
# *   This program is free software; you can redistribute it and/or modify  *
# *   it under the terms of the GNU General Public License as published by  *
# *   the Free Software Foundation; either version 2 of the License, or     *
# *   (at your option) any later version.                                   *
# *                                                                         *
# *   This program is distributed in the hope that it will be useful,       *
# *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
# *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
# *   GNU General Public License for more details.                          *
# *                                                                         *
# *   You should have received a copy of the GNU General Public License     *
# *   along with this program; if not, write to the                         *
# *   Free Software Foundation, Inc.,                                       *
# *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
# ***************************************************************************

import re,sys,os,subprocess,optparse

from math import *

PROGRAM_VERSION = '1.4'

# a simple Cartesian data pair class

class Cart:
  def __init__(self,x,y):
    self.x = float(x)
    self.y = float(y)
  def __repr__(self):
    return '{%f,%f}' % (self.x,self.y)
    
# a class that creates an HTML or PDF page
    
class OutputFormatter:
  
  css_style_block = """
      <meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>
      <style type=\'text/css\'>
        body * {
          font-family: Verdana, Tahoma, Helvetica, Arial;
        }
        table {
          border-collapse:collapse;
        }
        td {
          border:1px solid #808080;
          padding-right:4px;
          padding-left:4px;
          white-space: nowrap;
          text-align:right;
        }
        th {
          border:1px solid #808080;
          padding-right:4px;
          padding-left:4px;
          background:#ffffe0;
          font-weight: bold;
          text-align: center;
        }
        .cell0 { background:#eff1f1; }
        .cell1 { background:#f8f8f8; }
      </style>"""
      
  xml_tab_str = '  '
    
  def wrap_tag(self,tag,data,mod = ''):
    if(len(mod) > 0): mod = ' ' + mod
    return '<%s%s>\n%s\n</%s>\n' % (tag,mod,data,tag)

  def beautify_xhtml(self,data,otab = 0):
    tab = otab
    xml = []
    data = re.sub('\n+','\n',data)
    for record in data.split('\n'):
      record = record.strip()
      outc = len(re.findall('</|/>',record))
      inc = len(re.findall('<\w',record))
      net = inc - outc
      tab += min(net,0)
      xml.append((self.xml_tab_str * tab) + record)
      tab += max(net,0)
    if(tab != otab):
      sys.stderr.write('Error: tag mismatch: %d\n' % tab)
    return '\n'.join(xml) + '\n'

  def encode_xhtml(self,title,extra,array):
    table = ''
    row_num = 0
    for record in array:
      row = ''
      for field in record:
        field = field.strip()
        # strip quotes
        field = re.sub('"','',field)
        # replace blank field with &nbsp;
        field = re.sub('^$','&nbsp;',field)
        row += self.wrap_tag(('td','th')[row_num == 0],field)
      row_num += 1
      mod = ('class="cell%d"' % (row_num % 2),'')[row_num == 0]
      table+= self.wrap_tag('tr',row,mod)
    table = self.wrap_tag('table',table, 'cellpadding="0" cellspacing="0" border="0"')
    footer = 'Created using '
    footer += self.wrap_tag('a','TankProfiler','href=\"http://arachnoid.com/TankProfiler\"')
    footer += ' &mdash; copyright &copy; 2012,  P. Lutus'
    
    if(extra):
      extra = re.sub('\n','<br/>\n',extra)
      table = self.wrap_tag('p',extra) + table
    if(title):
      table = self.wrap_tag('h3',title) + table
    table += self.wrap_tag('p',footer)
    page = self.wrap_tag('div',table,'align="center"')
    page = self.wrap_tag('body',page)
    head = ''
    if(title):
      head += self.wrap_tag('title',title)
    page = self.wrap_tag('head',head + self.css_style_block) + page
    page = self.wrap_tag('html',page)
    page = self.beautify_xhtml(page)
    return page
    
 # must have 'wkhtmlpdf' PDF formatter installed:
 # http://code.google.com/p/wkhtmltopdf/
 
  def encode_pdf(self,html_title,extra,array):
    try:
      out = self.encode_xhtml(html_title,extra,array)
      p = subprocess.Popen(
        ['wkhtmltopdf','-q','--footer-center', '[page]/[toPage]', '-', '-']
      ,stdin=subprocess.PIPE
      ,stdout=subprocess.PIPE
      ,stderr=subprocess.PIPE)
      stdoutdata, stderrdata = p.communicate(out)
      return stdoutdata
    except Exception as e:
      sys.stderr.write('Error: %s (possibly no pdf converter present)\n' % e)
      return ''

class TankProfiler:
  
  # read file/stream containing tank profile
  # in a robust, tolerant way
  
  def read_profile(self,path):
    profile = []
    if(path == '-'):
      data = sys.stdin.read()
    else:
      with open(path,'r') as f:
        data = f.read()
    n = 0
    a = None
    b = 0
    # strip comment lines
    data = re.sub('(?sm)^#.*?$','',data)
    # split on non-numeric chars
    array = re.split('(?is)[^\d\.e+-]+',data)
    for s in array:
      try:
        b = float(s)
        if(a != None and n % 2 != 0):
          profile.append(Cart(a,b))
        a = b
        n += 1
      except:
        None
    return profile
    
  # circle area partitioned by y
  # -r <= y <= r
    
  def circle_sec(self,y,r):
    if(r == 0): return 0
    else:
      if(y > r): return pi * r * r
      elif(y < -r): return 0
      else:
        y /= r
        return (y * sqrt(1-y*y) + acos(-y)) * r * r
        
  # full tank volume, closed-form
  # treat each section as a cylinder if a.y == b.y
  # or a conical frustum if not
  
  def full_volume(self):
    a = False
    s = 0
    for b in self.profile:
      # if paired and nonzero length
      if(a and a.x != b.x):
        length = b.x - a.x
        # if a cylindrical section
        if(a.y == b.y):
          s += pi * a.y*a.y * length
        else:
          # conical frustum
          s += pi * length * (a.y * a.y + a.y * b.y + b.y * b.y) / 3.0
      a = b
    return s
    
  # partial frustum volume partitioned by y
  
  def partitioned_frustum(self,r_a,r_b,h,y):
    return ((3*h**2*r_a**2*y
    + (r_a**2 - 2*r_a*r_b + r_b**2)*y**3
    - 3*(h*r_a**2 - h*r_a*r_b)*y**2)*pi/h**2) / 3.0
   
  # frustum surface area
   
  # from http://mathworld.wolfram.com/ConicalFrustum.html
  
  # this area doesn't include the top and bottom disc areas
  # but it does manage the extreme cases of no height (disk)
  # and same radii (cylinder)
  
  def frustum_area(self,ra,rb,h):
    return pi * (ra+rb) * sqrt((ra-rb)**2 + h**2)
    
  # tank surface area including possible unclosed initial and final radii
  
  def tank_surface_area(self):
    a = False
    area = 0
    for b in self.profile:
      if(a):
        area += self.frustum_area(b.y,a.y,b.x-a.x)
      else:
        # add initial nonzero radius 
        area += pi * b.y**2
      a = b
    # add final nonzero radius
    area += pi * b.y**2
    return area
      
  # numerical integral horizontal volume partitioned by y
    
  def num_vol_horiz(self,y):
    a = False
    s = 0
    for b in self.profile:
      # if paired and nonzero length
      if(a and a.x != b.x):
        length = b.x - a.x
        # if a cylindrical section
        if(a.y == b.y):
          s += self.circle_sec(y,a.y) * length
        # else model numerically
        else:
          ss = 0
          step = (b.y - a.y) / self.intsteps
          r = a.y + step * 0.5
          for i in range(int(self.intsteps)):
            ss += self.circle_sec(y,r)
            r += step
          s += ss * length / self.intsteps
      a = b
    return s
    
  # numerical integral vertical volume partitioned by y
    
  def num_vol_vert(self,y):
    a = False
    s = 0
    for b in self.profile:
      # only include volumes below y level
      if(a and a.x != b.x and (b.x < y or a.x < y)):
        span = ((y >= a.x and y <= b.x) or (y >= b.x and y <= a.x))
        if(span):
          if(b.x < a.x):
            height = b.x - y
          else:
            height = y - a.x
        else:
          height = b.x - a.x
        # if a cylindrical section
        if(a.y == b.y):
          s += pi * a.y * a.y * height
        # else model as conical frustum
        else:
          if(span):
            s += self.partitioned_frustum(a.y,b.y,b.x - a.x,height)
          else:  
            s += pi * height * (a.y * a.y + a.y * b.y + b.y * b.y) / 3
      a = b
    return s
    
  # root finder for y from volume argument
  
  def num_y(self,v,f):
    epsilon = 1e-8
    y = 0
    dy = self.ymax
    while(abs(dy) > epsilon):
      dv = f(y)
      dy *= 0.5
      y += (dy,-dy)[dv>v]
    return y
    
  def emit_prog(self):
    if(self.show_prog):
      sys.stderr.write('.')
      sys.stderr.flush()
  
  def print_format(self,a,b,m):
    return ('%.*f,%.*f,%.*f' %
    (self.decimals,a,self.decimals,b,self.decimals,b*100/m))
    
  def to_array(self,s):
    self.data_table.append(s.split(','))
    
  def height_to_volume(self,y,ff):
    v = ff(y) * self.vol_mult
    return self.print_format(y+self.soffset,
    v / self.conv+self.voffset,
    self.fv * self.vol_mult+self.voffset)
    
  def height_table(self,ff):
    fs = self.height_to_volume(self.ymax,ff)
    y = self.ymin
    while(y <= self.ymax):
      self.emit_prog()
      s = self.height_to_volume(y,ff)
      self.to_array(s)
      y += self.step
    if(fs != s):
      self.to_array(fs)
      
  def volume_to_height(self,v,ff):
    y = self.num_y(v * self.conv,ff)
    return self.print_format(
      v * self.vol_mult+self.voffset,
      y+self.soffset,
      self.ymax+self.soffset)
      
  def volume_table(self,ff):
    fs = self.volume_to_height(self.fv,ff)
    s = 0
    v = 0
    while(v <= self.fv):
      self.emit_prog()
      s = self.volume_to_height(v,ff)
      self.to_array(s)
      v += self.step
    if(fs != s):
      self.to_array(fs)

  def fill_plot(self,plt,px,py):
    plt.plot(py,px,'-',color='blue')
    plt.fill(py,px,color='#f0f0f0')
    if(self.plot_points):
      plt.plot(py,px,'o',color='red')
        
  # must have matplotlib installed
  
  def show_save_graph(self):
    import matplotlib.pyplot as plt
    # set 1:1 aspect ratio
    plt.figure().add_subplot(111).set_aspect(1)
    px = []
    py = []
    yl = 1e30
    yh = -1e30
    xl = 1e30
    xh = -1e30
    # draw positive-coordinate side
    for p in self.profile:
      xl = min(xl,p.x)
      xh = max(xh,p.x)
      yl = min(yl,-p.y)
      yh = max(yh,p.y)
      px.append(p.x)
      py.append(p.y)
    # draw opposite side
    for p in self.profile[::-1]:
      px.append(p.x)
      py.append(-p.y)
    # close the figure
    px.append(p.x)
    py.append(p.y)
    border = (xh-xl) / 10
    if(self.vmode):
      plt.ylim(xl - border,xh + border)
      plt.xlim(yl - border,yh + border)
      plt.plot([0,0],[xl-border,xh+border],'-',color='green')
      self.fill_plot(plt,px,py)
    else:
      plt.xlim(xl - border,xh + border)
      plt.ylim(yl - border + self.ymax,yh + border + self.ymax)
      plt.plot([xl-border,xh+border],[self.ymax,self.ymax],'-',color='green')
      py = [y + self.ymax for y in py]
      self.fill_plot(plt,py,px)
    plt.grid(True)
    plt.savefig('tank.png')
    plt.show()
      
  def output(self):
    if(self.show_prog):
      sys.stderr.write('\n')
    out = ''
    if(self.pdf_mode):
      out = OutputFormatter().encode_pdf(self.html_title,self.extra,self.data_table)
    elif(self.html_mode):
      out = OutputFormatter().encode_xhtml(self.html_title,self.extra,self.data_table)
    else:
      if(self.extra):
        out = self.extra
      for line in self.data_table:
        out += ','.join(line) + '\n'
    if(self.ofile):
      with open(self.ofile,'w') as f:
        f.write(out)
    else:
      sys.stdout.write(out)
      
  def version(self):
    print('  %s version %s' % (self.__class__.__name__,PROGRAM_VERSION))
    print('  %s is © 2012, P. Lutus' % self.__class__.__name__)
    print('  For full documentation and latest version, visit')
    print('  http://arachnoid.com/%s' % self.__class__.__name__)
    
      
  def __init__(self):
    
    # force help if no args provided
    
    if(len(sys.argv) < 2):
      sys.argv.append('-h')
    
    # a default example tank with radius 30 units,
    # a central, 100-unit cylindrical section,
    # and two end frusta with height 30 units
    
    self.profile = (Cart(0,10),Cart(30,30),Cart(130,30),Cart(160, 10))
 
    self.vmode = False
    self.revmode = False
    self.show_prog = True
    self.intsteps = 1000.0   
    self.step = 1.0
    self.vol_mult = 1.0
    self.soffset = 0.0
    self.voffset = 0.0
    self.conv = 1.0
    self.decimals = 4
    self.ofile = False
    self.html_mode = False
    self.html_title = False
    self.pdf_mode = False
    self.plot_points = False
    self.arg = False
    self.extra = False
    self.wall = 0
    
    self.data_table = []
    
    self.parser =  optparse.OptionParser()
    self.parser.add_option("-a", "--arg", dest="arg", help="use argument to compute a single height/volume result instead of a table")
    self.parser.add_option("-c", "--conv", dest="conv", help=u"volume conversion factor (cm³ to liters = 1000, in³ to gallons = 231), default %.4f" % self.conv)
    self.parser.add_option("-d", "--decimals", dest="dec", help="result decimal places, default %d" % self.decimals)
    self.parser.add_option("-D", "--draw", action="store_true",dest="draw", help="draw an image of the tank (and save 'tank.png')")
    self.parser.add_option("-f", "--ifile", dest="ifile", help="input file path containing tank profile (x,y value pairs) ('-' = stdin)")
    self.parser.add_option("-F", "--ofile", dest="ofile", help="output file path (default: stdout)")
    self.parser.add_option("-H", "--html", action="store_true",dest="html", help="use HTML output format (default: CSV)")
    self.parser.add_option("-i", "--intsteps", dest="intsteps", help="integration steps (only used in horizontal mode), default %.0f" % self.intsteps)
    self.parser.add_option("-m", "--mult", dest="mult", help="volume multiplier for oval tanks, default %.4f" % self.vol_mult)
    self.parser.add_option("-n", "--noprog", action="store_true",dest="noprog", help="don't print progress chars")
    self.parser.add_option("-P", "--pdf", action="store_true",dest="pdf", help="use PDF output format (default: CSV)")
    self.parser.add_option("-p", "--plot", action="store_true",dest="plot", help="plot coordinate points in graphic image")
    self.parser.add_option("-r", "--rev", action="store_true",dest="rev", help="create volume -> height table (default height -> volume) -- slower.")
    self.parser.add_option("-s", "--step", dest="step", help="table increment size, default %.4f" % self.step)
    self.parser.add_option("-S", "--soff", dest="soffset", help="sensor offset value, default %.4f" % self.soffset)
    self.parser.add_option("-t", "--title", dest="title", help="HTML/pdf page title (default none)")
    
    self.parser.add_option("-u", "--upright", action="store_true",dest="upright", help="upright (vertical) tank mode (x and y reversed)")
    self.parser.add_option("-v", "--version", action="store_true",dest="version", help="program version and additional information")
    self.parser.add_option("-V", "--voff", dest="voffset", help="volume offset value, default %.4f" % self.voffset)
    self.parser.add_option("-w", "--wall", dest="wall", help="tank wall thickness, default %.4f" % self.wall)
    self.parser.add_option("-x", "--extra", action="store_true",dest="extra", help="compute full volume, surface area, surface volume (if wall thickness provided)")
    
    (options, args) = self.parser.parse_args()
    if(options.upright):
      self.vmode = True
    if(options.arg):
      self.arg = float(options.arg)
    if(options.rev):
      self.revmode = True
    if(options.ifile):
      self.profile = self.read_profile(options.ifile)
    if(options.ofile):
      self.ofile = options.ofile
    if(options.html):
      self.html_mode = True
    if(options.pdf):
      self.pdf_mode = True
    if(options.plot):
      self.plot_points = True
    if(options.title):
      self.html_title = options.title
    if(options.intsteps):
      self.intsteps = float(options.intsteps)
    if(options.step):
      self.step = float(options.step)
    if(options.mult):
      self.vol_mult = float(options.mult)
    if(options.dec):
      self.decimals = int(options.dec)
    if(options.soffset):
      self.soffset = float(options.soffset)
    if(options.voffset):
      self.voffset = float(options.voffset)
    if(options.conv):
      self.conv = float(options.conv)
    if(options.wall):
      self.wall = float(options.wall)
    if(options.noprog):
      self.show_prog = False
    
    # get analytical full volume
    
    self.fv = self.full_volume() / self.conv
          
    # find y min and max in profile
    
    self.ymax = -1e30
    self.ymin = 1e30
    self.xmin = 1e30
    
    if(self.vmode):
      ff = self.num_vol_vert
      for p in self.profile:
        self.ymax = max(self.ymax,p.x)
        self.ymin = min(self.ymin,p.x)
      # set bottom/left tank coordinate to zero
      for p in self.profile:
        p.x -= self.ymin
      self.ymax -= self.ymin
      self.ymin = 0
    else:
      ff = self.num_vol_horiz
      for p in self.profile:
        self.xmin = min(self.xmin,p.x)
        self.ymax = max(self.ymax,p.y)
        self.ymin = min(self.ymin,-p.y)
      self.soffset += self.ymax
      # reset left boundary to zero
      for p in self.profile:
        p.x -= self.xmin
        
    if(options.version):
      self.version()
      quit()

    if(options.extra):
      self.extra = 'Tank interior volume: %.*f units³\n' % (self.decimals,self.fv * self.vol_mult)
      sa = self.tank_surface_area()
      self.extra += 'Tank surface area: %.*f units²\n' % (self.decimals,sa * self.vol_mult)
      if(self.wall != 0):
        self.extra += 'Tank surface volume: %.*f units³\n' % (self.decimals,sa * self.vol_mult * self.wall)

    if(options.draw):
      self.show_save_graph()
    else:
      if(self.fv < 0):
        if(not self.extra):
          self.extra = ''
        self.extra += 'Error: volume < 0\n'
      else:
        if(self.revmode):
          self.to_array('Volume,Height,%')
        else:
          self.to_array('Height,Volume,%')
        if(self.arg):
          # single argument
          if(self.revmode):
            s = self.volume_to_height(self.arg,ff)
          else:  
            s = self.height_to_volume(self.arg-self.soffset,ff)
          self.to_array(s)
        else:
          # generate table
          if(self.revmode):
            self.volume_table(ff)
          else:  
            self.height_table(ff)
    self.output()
  
TankProfiler()