#!/usr/bin/env texlua

-- $Id: luagraphlib.lua 27 2024-07-26 18:01:43Z reinhard $

-- Copyright (C) 2020 Reinhard Kotucha <reinhard.kotucha@web.de>
-- 
-- 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., 51 Franklin Street, Fifth Floor, Boston,
-- MA 02110-1301, USA.

module(..., package.seeall)

local pdfcomments = true
local latelua = false

local fmt  = string.format
local pi   = math.pi
local sin  = math.sin
local cos  = math.cos
local tan  = math.tan
local deg  = math.deg
local rad  = math.rad
local sqrt = math.sqrt
local abs  = math.abs
local arc_tan = math.atan2


-- debugging

function wlog (...) -- write to log file.  Similar to \wlog in TeX
  local args = {...}
  local arglist = {}
  for _,arg in ipairs(args) do
    if type(arg) == 'nil' then
      arglist[#arglist+1] = '<nil>'
    else
      arglist[#arglist+1] = arg
    end
  end
  texio.write_nl('=> '..table.concat(arglist, ' ')..'\n')
end


-- utility functions

local function round (n, prec)
  -- prec is the number of fractional digits
  local rounded = string.format('%.'..prec..'f', n)
  -- strip trailing zeroes
  rounded = (tostring(rounded):gsub('0+$', ''))
  -- strip decimal point if number is integer
  rounded = (rounded:gsub('%.$', ''))
  return rounded
end


-- Polar/Cartesian transformations

-- In PostScript angles are in degrees.  For further processing angles
-- are sometimes needed in radians.  We also provide functions which
-- return angles in radians in order to avoid back and forth
-- conversions.

function cartesian_rad (r, phi) -- angle in radians
  return r * cos(phi), r * sin(phi)
end


function cartesian (r, phi) -- angle in degrees
  return r * cos(rad(phi)), r * sin(rad(phi))
end


function polar_rad (x, y) -- angle in radians
  local r = sqrt(x^2 + y^2)
  local phi = arc_tan(y,x)
  return r, phi
end


function polar (x, y)  -- angle in degrees
  local r, phi = polar_rad (x, y)
  return r, deg(phi)
end


local function is_odd (n)
  if n%2 == 1 then return true end
end


function deepcopy (orig) -- http://lua-users.org/wiki/CopyTable
  local orig_type = type(orig)
  local copy
  if orig_type == 'table' then
    copy = {}
    for orig_key, orig_value in next, orig, nil do
      copy[deepcopy(orig_key)] = deepcopy(orig_value)
    end
    setmetatable(copy, deepcopy(getmetatable(orig)))
  else -- number, string, boolean, etc
    copy = orig
  end
  return copy
end


function readdatafile (file, lambda)
  local retval = {}
  
  if not lfs.isfile(file) then
    error('\nERROR: File "'..file..'" not found.')
  end

  local fh = assert(io.open(file))

  for line in fh:lines() do
    -- Replace possible non-space delimiters by spaces.  This allows
    -- to process various formats like GNU Octave ASCII data files,
    -- "CSV", etc. without interaction.
    line = line:gsub('[\t:;,%%!]', ' ')

    -- Remove all spaces at the beginning of a line.
    line = line:gsub('^%s+', '')

    -- Ignore all lines which don't begin with a number.  All comments
    -- and empty lines in data files are ignored.
    if line:match('^%-?%.?[0-9]') then

      -- The result must be a table with two columns.  If no lambda
      -- expression is provided, the first two colums of the data file
      -- are returned.  With a lambda expression as second argument
      -- arbitrary columns can be selected from a data file and can be
      -- pre-processed in any way possible.
      local a, b
      local cols = line:explode (' +')

      -- Convert strings to numbers.  Because lambda expressions can
      -- access any column we convert all table entries from strings
      -- to numbers if possible instead of only the return values.
      for i, col in ipairs(cols) do
        cols[i] = tonumber(cols[i])
      end

      if lambda then
        a, b = lambda(cols)
        if a and b then
          retval[#retval + 1] = {a, b}
        end
      else
        retval[#retval + 1] = cols
      end
    end
  end
  fh:close()
  return retval
end


-- Convert other units to millimeters
-- We cannot use 'in' for 'inch' because 'in' is a keyword in Lua.

function pt (v)   return v * 25.4/72.27 end
function bp (v)   return v * 25.4/72 end
function pc (v)   return 12 * pt(v) end
function inch (v) return v * 25.4 end
function mm (v)   return v end
function cm (v)   return 10 * v end
function sp (v)   return pt(v) / 65536 end
function dd (v)   return pt(v) * 1238/1175 end
function cc (v)   return 12 * dd(v) end
-- For definitions of 'nd' and 'nc' see 'The LaTeX3 Interfaces'.
function nd (v)   return 0.375 * mm(v) end
function nc (v)   return 12 * nd(v) end


-- Write s to PDF file.  s is either a string or a list of strings
function to_pdf (s)
  if type(s) == 'string' then
    if latelua then
      pdf.print('page', s..'\n')
    else
      tex.print([[\pdfextension literal {]]..s..'}')
    end
  elseif type(s) == 'table' then
    if latelua then
      pdf.print('page', table.concat(s, ' ')..'\n')
    else
      tex.print([[\pdfextension literal {]]..table.concat(s, ' ')..'}')
    end
  end
end

-- Here a path is a list of strings where each string has the form
--   { x y 'm' } |
--   { x y 'l' } |
--   { x1 y1 x1 y2 x3 y3 'c' } 
--
-- closepath is not part of the path and can be added later when required.

local function path2pdf (p)
  for  _,line in ipairs(p) do
    for i,token in ipairs(line) do
      if tonumber(token) then
        -- round to avoid fp numbers
        line[i] = round(line[i], 3)
      end
    end
    to_pdf(table.concat(line, ' '))
  end
end


-- add string s as PDF comment 
function pdfcomment (s)
    if pdfcomments then
    to_pdf{'%% ', s}
  end
end


-- matrix operations

local function setmatrix (matrix)
  for i,v in ipairs(matrix) do
    matrix[i] = round(matrix[i], 4) -- we need 4 fractional digits here
  end
  matrix[#matrix+1] = 'cm'
  to_pdf(matrix)
end


function scale (x, y)
  setmatrix{x, 0, 0, y, 0, 0}
end

function translate (x, y)
  setmatrix{1, 0, 0, 1, x, y}
end


-- rotate, phi in radians.
function rotate_rad (phi)
  setmatrix{ cos(phi), sin(phi),
            -sin(phi), cos(phi), 0, 0}
end


-- PostScript assumes angles in degrees
function rotate (phi)
  rotate_rad (rad(phi))
end


-- graphics state

--[==============================================================[
  Unlike PostScript, PDF's doesn't maintain paths in its graphics
  state and there is no concept of currentpoint.
  
  We maintain an extension of the graphics stack here in order to
  provide rmoveto, rlineto, and rcurveto.  Furthermore, when paths are
  maintained by the graphics state stack we can change the order of
  fill and stroke, etc.

  A path is an object which is usually constructed by procedures like
  [r]moveto, [r]lineto, [r]curveto.  closepath is not part of the
  path.

  A path can also be constructed by a custom function, for instance
  from files generated by external programs.

  The value of path at the top of the graphics state stack is used by
  painting operators like fill, stroke, and paint.

  We also maintain currentpoint in the graphics state stack.

    200 200 moveto
    gsave
      0 100 rlineto stroke
    grestore
    100 0 rlineto stroke


--]==============================================================]



-- initialize graphics state stack

gstate = {} -- graphics state stack

local initial_gstate = {
  currentpoint = {x = 0, y = 0},
  path = {},
}

gstate[0] = initial_gstate

-- gsave and grestore

function gsave ()
  to_pdf{'q'}
  gstate[#gstate+1] = deepcopy (gstate[#gstate])
end


function grestore ()
  to_pdf{'Q'}
  gstate[#gstate] = nil
end


-- maintain currentpoint

function currentpoint ()
  return gstate[#gstate].currentpoint.x, gstate[#gstate].currentpoint.y
end

function setcurrentpoint (currentpoint)
  gstate[#gstate].currentpoint = currentpoint
end


-- maintain path

function current_path ()
  return gstate[#gstate].path
end

local function pop_current_path ()
  gstate[#gstate].path = {}
end

-- add path_obj to top level path

local function push_path (path_obj)
  local path = current_path ()
  path[#path+1] = path_obj
end


-- drawing attributes

function setlinewidth (lw)
  to_pdf{lw, 'w'}
end


function setlinecap (cap)
  local linecap = {
    butt = 0,
    round = 1,
    square = 2
  }
  if type(cap) == 'string' then
    to_pdf{linecap[cap], 'J'}
  else
    to_pdf{cap, 'J'}
  end
end


function setlinejoin (join)
  local linejoin = {
    miter = 0,
    round = 1,
    bevel = 2
  }
  if type(join) == 'string' then
    to_pdf{linejoin[join], 'j'}
  else
    to_pdf{join, 'j'}
  end
end


function setmiterlimit (limit)
  if limit < 1 then
    error ('miterlimit < 1')
  end
  to_pdf{limit, 'M'}
end


function setdash (pattern, phase)
  -- pattern: table of numbers,
  -- phase:   number
  if not phase then phase = 0 end
  to_pdf{'['..table.concat(pattern, ' ')..']', phase, 'd'}
end


local function setcolor (color_arg, stroke)
  local color = deepcopy(color_arg)
  local colorspace = {'g', 'nil', 'rg', 'k'}
  if stroke then
    colorspace = {'G', 'nil', 'RG', 'K'}
  end
  color[#color+1] = colorspace[#color]
  to_pdf(color)
end      

function setstrokecolor (color) setcolor(color, true)  end
function setfillcolor   (color) setcolor(color, false) end

-- basic drawing commands

function moveto (x, y)
  push_path {x, y, 'm'}
  setcurrentpoint {x = x, y = y}
end

function lineto (x, y)
  push_path {x, y, 'l'}
  setcurrentpoint {x = x, y = y}
end

function curveto (x1, y1, x2, y2, x3, y3)
  push_path {x1, y1, x2, y2, x3, y3, 'c'}
  setcurrentpoint {x = x3, y = y3}
end

function rmoveto (dx, dy)
  local x, y = currentpoint()
  push_path {x + dx, y + dy, 'm'}
  setcurrentpoint {x = x + dx, y = y + dy}
end

function rlineto (dx, dy)
  local x, y = currentpoint()
  push_path {x + dx, y + dy, 'l'}
  setcurrentpoint {x = x + dx, y = y + dy}
end

function rcurveto (x1, y1, x2, y2, x3, y3)
  local x, y = currentpoint()
  push_path {x + x1, y + y1, x + x2, y + y2, x + x3, y + y3, 'c'}
  setcurrentpoint {x = x + x3, y = y + y3}
end

local closepath_p  = false

function closepath ()
  closepath_p = true
end

-- higher level drawing commands

function rect (x, y, w, h)
  push_path {x, y, w, h, 're'}
end

function rectn (x, y, w, h)
  moveto (x, y)
  rlineto (0, h)
  rlineto (w, 0)
  rlineto (0, -h)
  closepath()
end  

function circle (r)
  -- bezier control points for 90° segment at 4r/3 * (sqrt(2)-1) 
  local k = r * 4/3 * (sqrt(2) - 1)

  rmoveto(r, 0)
  rcurveto(0, k, -r+k, r, -r, r)
  rcurveto(-k, 0, -r, -r+k, -r, -r)
  rcurveto(0, -k, r-k, -r, r, -r)
  rcurveto(k, 0, r, r-k, r, r)
  closepath()
end


--[[
function circle (r)
  -- this function is less efficient due to many polar-cartesian transformations
  local x, y = currentpoint()
  rmoveto(r, 0)
  arc(x, y, r, 0, 360)
  closepath()
end
--]]

function circlen (r) -- negative circle (drawn clockwise)
  local x, y = currentpoint()
  rmoveto(r, 0)
  arcn(x, y, r, 360, 0)
  closepath()
end


-- Draw arc counterclockwise.  Make sure that
--
--  *  ang2 > ang1 
--  *  ang2 - ang1 < 360°
--
-- Maybe it's best to add multiples of 360° first until both angles
-- are positive.

function arc (x, y, r, ang1, ang2, arrowspec)

  -- make sure that ang2 > ang1
  while ang2 < ang1 do ang2 = ang2 + 360 end
  -- make sure that ang2 - ang1 <= 360
  if ang2 - ang1 > 360 then ang2 = ang2 - 360 end
  
  --  print('ang1 = '..ang1..', ang2 = '..ang2..', ang2 - ang1 = '..ang2 - ang1)

  -- local arrowtarget = cartesian (r, ang2)

  -- if arrowspec and arrowspec.len then
  --   -- subract 2/3 of arrowspec.len from circumference at end of the arc.
  --   --
  --   local delta = .99 * arrowspec.len
  --   if arrowspec.inset > 0 then
  --     delta = .99 * (arrowspec.len - arrowspec.inset)
  --   end
  --   local delta_phi = 360 * delta / (2 * pi * r)
    
  --   ang2 = ang2 - delta_phi    
  -- end
    
  local phi1 = rad(ang1)
  local phi2 = rad(ang2)

  local L = (ang2-ang1)/360                -- relative length of arc
  local k = 4*(sqrt(2) - 1)/3              -- position of control point
  local theta = arc_tan(k)                 -- angle of control point
  local R = sqrt(r^2 + (r*tan(L*theta))^2) -- magnitude of control point
  
  local function cartesian (a, phi)  -- cartesian() with radians
    return a * cos(phi), a * sin(phi)
  end
  
  local x1, y1 = cartesian(R, phi1 + L * theta)
  local x2, y2 = cartesian(R, phi1 + L * (pi/2-theta))
  local x3, y3 = cartesian(r, phi1 + L * (pi/2))

  curveto(x1+x, y1+y, x2+x, y2+y, x3+x, y3+y)

  local x4, y4 = cartesian(R, phi1 + L * (pi/2+theta))
  local x5, y5 = cartesian(R, phi1 + L * (pi-theta))
  local x6, y6 = cartesian(r, phi1 + L * pi)
  
  curveto(x4+x, y4+y, x5+x, y5+y, x6+x, y6+y)

  local x7, y7 = cartesian(R, phi1 + L * (pi+theta))
  local x8, y8 = cartesian(R, phi1 + L * (1.5*pi-theta))
  local x9, y9 = cartesian(r, phi1 + L * (1.5*pi))

  curveto(x7+x, y7+y, x8+x, y8+y, x9+x, y9+y)

  local x10, y10 = cartesian(R, phi1 + L * (1.5*pi+theta))
  local x11, y11 = cartesian(R, phi1 + L * (2*pi-theta))
  local x12, y12 = cartesian(r, phi1 + L * (2*pi))

  curveto(x10+x, y10+y, x11+x, y11+y, x12+x, y12+y)

--  if arrowspec then
--    local arrowhead = arrowshape (arrowspec)
--    moveto(arrowtarget)
--    arrowhead(0,0,9)
--  end
end

local function mark_xy (x, y)
  gsave()
    scale(.5,.5)
    tex.print([[\MTxcross]])
  grestore()
end

function arcn (x, y, r, ang1, ang2, arrowspec)

  -- make sure that ang2 < ang1
  while ang2 > ang1 do ang1 = ang1 + 360 end
  -- make sure that ang2 - ang1 <= 360
  if ang1 - ang2 > 360 then ang1 = ang1 - 360 end

  -- local cpx, cpy = currentpoint()
  
  -- print('\ncurrentpoint:'..cpx+x..', '..cpy+y)

  -- local arrow_x, arrow_y = cartesian (r, ang2)
  -- arrow_x = arrow_x - x

  -- print('\narrowpos:'..arrow_x..', '..arrow_y)
  
  
  if arrowspec and arrowspec.len then
    -- subract 2/3 of arrowspec.len from circumference at end of the arc.
    --
    local delta = .99 * arrowspec.len
    if arrowspec.inset > 0 then
      delta = .99 * (arrowspec.len - arrowspec.inset)
    end
    local delta_phi = 360 * delta / (2 * pi * r)
    ang2 = ang2 + delta_phi    
  end

  local phi1 = rad(ang1)
  local phi2 = rad(ang2)

  local L = (ang1-ang2)/360                -- relative length of arc
  local k = 4*(sqrt(2) - 1)/3              -- position of control point
  local theta = arc_tan(k)                 -- angle of control point
  local R = sqrt(r^2 + (r*tan(L*theta))^2) -- magnitude of control point
  
  local function cartesian (a, phi)  -- cartesian() with radians
    return a * cos(phi), a * sin(phi)    
  end
  
  local x1, y1 = cartesian(R, phi1 - L * theta)
  local x2, y2 = cartesian(R, phi1 - L * (pi/2-theta))
  local x3, y3 = cartesian(r, phi1 - L * (pi/2))

  curveto(x1+x, y1+y, x2+x, y2+y, x3+x, y3+y)

  local x4, y4 = cartesian(R, phi1 - L * (pi/2+theta))
  local x5, y5 = cartesian(R, phi1 - L * (pi-theta))
  local x6, y6 = cartesian(r, phi1 - L * pi)
  
  curveto(x4+x, y4+y, x5+x, y5+y, x6+x, y6+y)

  local x7, y7 = cartesian(R, phi1 - L * (pi+theta))
  local x8, y8 = cartesian(R, phi1 - L * (1.5*pi-theta))
  local x9, y9 = cartesian(r, phi1 - L * (1.5*pi))

  curveto(x7+x, y7+y, x8+x, y8+y, x9+x, y9+y)

  local x10, y10 = cartesian(R, phi1 - L * (1.5*pi+theta))
  local x11, y11 = cartesian(R, phi1 - L * (2*pi-theta))
  local x12, y12 = cartesian(r, phi1 - L * (2*pi))

  curveto(x10+x, y10+y, x11+x, y11+y, x12+x, y12+y)

  if arrowspec then
    local arrowhead = arrowshape (arrowspec)

    -- gsave()
    -- setlinewidth(pt(0.4))
    -- stroke()
    -- moveto(arrow_x, arrow_y)
    -- print('\narrowpos:'..arrow_x..', '..arrow_y)
    -- mark_xy (arrow_x, arrow_y)

    -- arrowhead(0,0,9)
    -- grestore()
  end
end


-- path painting operators

function stroke ()
  path2pdf(current_path())
  if closepath_p then
    to_pdf {'s'}
    closepath_p = false
  else
    to_pdf {'S'}
  end
  pop_current_path()
end 


function fill ()
  path2pdf(current_path())
  to_pdf {'h f'}
  closepath_p = false
  pop_current_path()
end


function eofill ()
  path2pdf(current_path())
  to_pdf {'h f*'}
  closepath_p = false
  pop_current_path()
end


function paint (s)
  local eofill = false
  if s and s:match('^eo') then eofill = true end
  
  path2pdf(current_path())
  if closepath_p then
    if eofill then to_pdf {'b*'} else to_pdf {'b'} end
    closepath_p = false
  else
    if eofill then to_pdf {'B*'} else to_pdf {'B'} end
  end
  pop_current_path()
end

-- clipping operators

function clip ()
  to_pdf {'W n'}
end

function eoclip ()
  to_pdf {'W* n'}
end

function rectclip (x, y, w, h)
  push_path {x, y, w, h, 're W n'}
end

-- text

function framedbox (wd, ht, dp, attrstring)
  wd = wd:gsub('pt$','')
  ht = ht:gsub('pt$','')
  dp = dp:gsub('pt$','')

  local attributes = attrstring:explode(';')
  local attr = {}

  if #attrstring > 0 then
    for _,attribute in pairs(attributes) do
      local row = attribute:explode('=')

      if row[2]:match('%[') then -- value is a table (color)
        row[2]=row[2]:gsub('%[', '')
        row[2]=row[2]:gsub('%]', '')
        attr[row[1]] = row[2]:explode(',')
      else
        if tonumber(row[2]) then
          attr[row[1]] = tonumber(row[2])
        else
          attr[row[1]] = row[2]
        end
      end
    end
  end

  local shift = 0 -- left and right alignments.  Text is centered by
                  -- default.
  if attr.align:match('^l') then
    shift = wd/2
  elseif attr.align:match('^r') then
    shift = -wd/2
  end
  
  local sep = attr.sep
    
  gsave()
    if attr.strokecolor then setstrokecolor(attr.strokecolor) end
    if attr.fillcolor then setfillcolor(attr.fillcolor) end
    if attr.linewidth then setlinewidth(attr.linewidth) end

    translate(shift,0)

    -- Below, sw, se, ne, and nw denote southwest, southeast,
    -- northeast, and northwest.
    
    local corner = {
      sw = {
        x = -wd/2 - sep,
        y = -dp - sep,
      },
      se = {
        x = wd/2 + sep,
        y = -dp - sep,
      },
      ne = {
        x = wd/2 + sep,
        y = sep + ht,
      },
      nw = {
        x = -wd/2 - sep,
        y = ht + sep,
      },
    }

    -- If radius == 0 we draw a normal rectangle which
    -- obeys setlinecap(), setlinejoin(), and setmiterlimit().

    if attr.radius > 0 then
      local r = attr.radius

      if r > wd/2 then
        r = wd/2
      elseif r > (ht+dp)/2 then
        r = (ht+dp)/2
      end

      moveto(corner.sw.x - r, corner.sw.y + r)
      arc(corner.sw.x, corner.sw.y, r, 180, 270)
      lineto(corner.se.x, corner.se.y - r) 
      arc(corner.se.x, corner.se.y, r, 270, 360)
      lineto(corner.ne.x + r, corner.ne.y) 
      arc(corner.ne.x, corner.ne.y, r, 0, 90)
      lineto(corner.nw.x + r, corner.nw.y + r) 
      arc(corner.nw.x, corner.nw.y, r, 90, 180)
    else
      moveto(corner.sw.x, corner.sw.y)
      lineto(corner.se.x, corner.se.y)
      lineto(corner.ne.x, corner.ne.y)
      lineto(corner.nw.x, corner.nw.y)
    end
    closepath()
    paint()
  grestore()
end


local function serialize (t)
  local strings = {}
  for k, v in pairs(t) do
    if type(v) == 'table' then
      strings[#strings+1] = k..'=['..table.concat(v, ',')..']'
    elseif type(v) == 'number' or type(v) == 'string' then
      strings[#strings+1] = k..'='..v
    end
  end
  return table.concat(strings, ';')
end
  

function t2s (t) -- convert table to string.
  return '['..table.concat(t, ':')..']'
end


function text (s, attr)
  local bp = 25.4/72.27
  if not attr then attr = {} end

  local a = {}
  a.scale       = attr.scale or 1
  a.slant       = attr.slant or 0
  a.extend      = attr.extend or 1
  a.align       = attr.align or 'c'
  a.strut       = attr.strut or ''
  a.textcolor   = attr.textcolor or {0}
  a.strokecolor = attr.strokecolor or {0}
  a.fillcolor   = attr.fillcolor or {.5}
  a.linewidth   = attr.linewidth or .4
  a.radius      = attr.radius or 0
  a.sep         = attr.sep or 0
  a.framed      = attr.framed or false
  
  local textbox
  
  if a.align == 'c' or a.align == 'center' then
    textbox = '\\clap{'..s..'}'
  elseif a.align == 'l' or a.align == 'left' then
    textbox = '\\rlap{'..s..'}'
  elseif a.align == 'r' or attr.align == 'right' then
    textbox = '\\llap{'..s..'}'
  end

  gsave()
    setmatrix{bp*a.scale*a.extend, 0, a.slant, bp*a.scale, 0, 0}
    
    if a.framed then
      tex.print('\\luaframedbox{'..s..'}{'..a.strut..'}{'..serialize(a)..'}')
    end
    setfillcolor(a.textcolor)
    tex.print(textbox)
  grestore()
end


-- path operations

function scalepath (x, y, x_orig, y_orig)
  x_orig = x_orig or 0
  y_orig = y_orig or 0
  local cp = current_path()
  for _, line in ipairs (cp) do
    for i,v in ipairs(line) do
      if tonumber (v) then
        if is_odd (i) then
          line[i] = (v - x_orig) * x + x_orig
        else
          line[i] = (v - y_orig) * y + y_orig
        end
      end
    end
  end
end


function translatepath (x, y)
  local cp = current_path()
  for _,line in ipairs (cp) do
    for i,v in ipairs(line) do
      if tonumber (v) then
        if is_odd (i) then
          line[i] = v + x
        else
          line[i] = v + y
        end
      end
    end
  end
end

function reversepath ()
  local path = current_path()

  local revpath={}
  local procedures={} 
  local coordinates={}

  local function push (t,...)
    local args = {...}
    for _,arg in ipairs(args) do
      t[#t+1] = arg
    end
  end

  local function pop (coordinates)
    local list={}
    for _,v in ipairs(coordinates[#coordinates]) do -- copy t[able]
      list[#list+1] = v
    end
    coordinates[#coordinates]= nil
    return list
  end
  
  for i,line in ipairs(path) do
    local proc = line[#line]
    procedures[#procedures+1] = proc

    if proc == 'm' then
      push(coordinates, {line[1], line[2]})
    elseif 
      proc == 'l' then
        push(coordinates, {line[1], line[2]})
    elseif   proc == 'c' then
      push(coordinates, {line[1], line[2]})
      push(coordinates, {line[3], line[4]})
      push(coordinates, {line[5], line[6]})
    end
  end

  for i,proc in ipairs(procedures) do
    if proc == 'm' then
      local x, y = table.unpack(pop(coordinates))
      push(revpath, {x, y, proc})
    elseif proc == 'l' then
      local x, y = table.unpack(pop(coordinates))
      push(revpath, {x, y, proc})
    elseif proc == 'c' then
      local x1, y1 = table.unpack(pop(coordinates))
      local x2, y2 = table.unpack(pop(coordinates))
      local x3, y3 = table.unpack(pop(coordinates))
      push(revpath, {x1, y1, x2, y2, x3, y3, proc})
    end
  end
--  for i,v in ipairs(revpath) do print (table.concat(v, ' ')) end
  gstate[#gstate].path = revpath
end


function datafile (filename, x, y)
  local data = assert(io.open(filename, 'r'))
  local coords = {}

  for line in data:lines() do
    line = line:gsub('[\t,;]', ' ')    -- delimiters
    line = line:gsub('^ *[%%#].*', '') -- strip comments
    line = line:gsub('^ *$', '')       -- strip whitespaces
    
    if #line > 0 then
      local cols = line:explode(' +')
      coords[#coords+1] = {cols[x], cols[y]}
    end
  end      
  data:close()
  return coords
end


function to_path (t)
  local firstpoint = true

  for _,cols in ipairs(t) do
    if firstpoint then
      push_path {cols[1], cols[2], 'm'}
      firstpoint = false
    else
      push_path {cols[1], cols[2], 'l'}
    end
  end
end


function polar_to_cartesian (t)
  local cartesian = {}
  
  for _,polar in ipairs(t) do
    local r = polar[1]
    local phi = polar[2] * pi/180
    cartesian[#cartesian+1] = {r * cos(phi), r * sin(phi)}
  end
  return cartesian
end

  
function smooth ()
  local path = current_path()
  local newpath = {}

  local function push (t)
    newpath[#newpath+1] = t
  end
  
  local pts={} -- put path coordinates into a table
  for _,p in ipairs(path) do
    pts[#pts+1] = { x = p[1], y = p[2]}
  end
  
  for i in ipairs(pts) do
    if i == 1 then
      local x = pts[i].x
      local y = pts[i].y 
      
      local x_r = pts[i+1].x - 2 * ((pts[i+1].x - pts[i].x)/3)
      
      -- dx/dy
      local d_r = (pts[i+1].y - pts[i].y) / (pts[i+1].x - pts[i].x)
    
      local y_r = y + (d_r * (pts[i+1].x - pts[i].x)/3)
    
      push {x, y, 'm'}
      push {x_r, y_r}
 
    elseif i > 1 and i < #pts then
      local x = pts[i].x
      local x_l = pts[i].x - (pts[i].x - pts[i-1].x)/3
      local x_r = pts[i].x + (pts[i+1].x - pts[i].x)/3
      
      -- dx/dy
      local d = (pts[i+1].y - pts[i-1].y) / (pts[i+1].x - pts[i-1].x)
      
      local y = pts[i].y 
      local y_l = y - (d * (pts[i].x - pts[i-1].x)/3)
      local y_r = y + (d * (pts[i+1].x - pts[i].x)/3)
      
      push {x_l, y_l}
      push {x, y, 'c'}
      push {x_r, y_r}
      
    elseif i == #pts then
      local x = pts[i].x
      local y = pts[i].y 
      
      local x_l = pts[i].x - (pts[i].x - pts[i-1].x)/3
      
      -- dx/dy
      local d_l  = (pts[i].y - pts[i-1].y) / (pts[i].x - pts[i-1].x)
      
      local y_l = y - (d_l * (pts[i].x - pts[i-1].x)/3)
      
      push {x_l, y_l}
      push {x, y, 'c'}
    end
  end
  gstate[#gstate].path = newpath
end


-- higher level functions


function polylineto (...) -- not tested
  local coords = {...}
  if is_odd(coords) then
    error('polylineto: odd number of arguments not allowed.')
  end
  for i=0, #coords/2 - 1 do
    lineto(coords[coords*2*i], coords[coords*2*i+1])
  end
end

-- TODO: ... either x1, y1, x2, y2, or {x1, y1}, {x2, y2},


function arrowhead (x, y, len, wd, inset, offset)
  local wd = wd or len/3
  gsave()
    translate(offset, 0)
    moveto(0,0)
    lineto(-len, wd/2)
    lineto(-len + inset, 0)
    lineto(-len, -wd/2)
    closepath()
    paint()
  grestore()
end


function arrowshape (len, wd, inset, offset)
  -- len: length of arrowhead
  -- width: width of arrowhead
  -- inset: retraction
  -- offset: distance between arrow tip and center

  -- retval: table of coordinate pairs
  --         offset 
  --         angle 

  return function (x, y, angle)
    local wd = wd or len/3
    gsave()
      translate(x, y)
      rotate(angle)
      moveto(offset,0)
      rlineto(-len, wd/2)
      rlineto(inset, -wd/2)
      rlineto(-inset, -wd/2)
      closepath()
      fill()
    grestore()
    return x, y, angle, ofset  
  end
end


function arrowshape (shape)
  local len    = shape.len
  local wd     = shape.wd or shape.len * 2/3
  local inset  = shape.inset or 0
  local offset = shape.offset or 0

  -- len: length of arrowhead
  -- width: width of arrowhead
  -- inset: retraction
  -- offset: distance between arrow tip and center

  -- retval: table of coordinate pairs
  --         offset 
  --         angle 

  return function (x, y, angle)
    gsave()
      translate(x, y)
      rotate(angle)
      moveto(offset,0)
      rlineto(-len, wd/2)
      rlineto(inset, -wd/2)
      rlineto(-inset, -wd/2)
      closepath()
      fill()
    grestore()
    return x, y, angle, offset  
  end
end



------------------------------------------------------------
-- Import
------------------------------------------------------------

-- import symbols to global namespace
function import()
  local t =_G[_NAME]
  for k,v in pairs(t) do 
    if type(v) ~= 'table' then _G[k] = v end
  end
end


-- Local Variables:
--  mode: Lua
--  TeX-engine: luatex
--  lua-indent-level: 2
--  indent-tabs-mode: nil
--  coding: utf-8-unix
-- End:
-- vim:set tabstop=2 expandtab:



--[[
function rectstroke (x, y, w, h) -- not tested
  push_path {x, y, w, h, 're S'}
end

function rectfill (x, y, w, h) -- doesnt't work
  push_path {x, y, w, h, 're h f'}
end
--]]


--[[
local function printf (...)
  print(fmt(...))
end
--]]


--[[

function datafile_to_path (filename, x, y)
  local data = assert(io.open(filename, 'r'))
  local firstpoint = true
  for line in data:lines() do
    line = line:gsub('\t', '')
    if line:match('^ *[%-%.0-9]') then -- ignore any line which starts with a number
                                       -- this needs improvement!
      local cols = line:explode(' +')
      if firstpoint then
        push_path {cols[x], cols[y], 'm'}
        firstpoint = false
      else
        push_path {cols[x], cols[y], 'l'}
      end
    end
  end
  data:close()
end
--]]


--[[
function table_to_path (t)
  local firstpoint = true
  for _, row in ipairs(t) do
    if firstpoint then
      push_path {row[1], row[2], 'm'}
      firstpoint = false
    else
      push_path {row[1], row[2], 'l'}
    end
  end
end
--]]





--[[
  if token == nil then  -- plain Lua, no TeX
    io.write('=> '..fmt(...)..'\n')
  else
    texio.write_nl('=> '..fmt(...)..'\n')
  end
--]]