diff options
Diffstat (limited to 'macros/luatex/latex/bezierplot')
-rw-r--r-- | macros/luatex/latex/bezierplot/README | 21 | ||||
-rw-r--r-- | macros/luatex/latex/bezierplot/bezierplot-doc.pdf | bin | 0 -> 251276 bytes | |||
-rw-r--r-- | macros/luatex/latex/bezierplot/bezierplot-doc.tex | 200 | ||||
-rwxr-xr-x | macros/luatex/latex/bezierplot/bezierplot.lua | 940 | ||||
-rw-r--r-- | macros/luatex/latex/bezierplot/bezierplot.sty | 17 |
5 files changed, 1178 insertions, 0 deletions
diff --git a/macros/luatex/latex/bezierplot/README b/macros/luatex/latex/bezierplot/README new file mode 100644 index 0000000000..493bfe5c65 --- /dev/null +++ b/macros/luatex/latex/bezierplot/README @@ -0,0 +1,21 @@ +DESCRIPTION: +bezierplot is a Lua program as well as a (Lua)LaTeX package. +Given a smooth function, bezierplot returns a smooth bezier path written +in TikZ notation (which also matches METAPOST) that approximates the +graph of the function. For polynomial functions of degree <= 3 and +inverses of them, the approximation is exact. bezierplot finds special +points such as extreme points and inflection points and reduces the +number of used points. + +VERSION: +1.4 2019-06-18 + +LICENSE: +The package and the program are distributed on CTAN under the terms of +the LaTeX Project Public License (LPPL) version 1.3c. + +Copyright (c) 2018 Linus Romer + +Please write to +linus dot romer at gmx dot ch +to submit bug reports, request new features, etc. diff --git a/macros/luatex/latex/bezierplot/bezierplot-doc.pdf b/macros/luatex/latex/bezierplot/bezierplot-doc.pdf Binary files differnew file mode 100644 index 0000000000..c93fd5f102 --- /dev/null +++ b/macros/luatex/latex/bezierplot/bezierplot-doc.pdf diff --git a/macros/luatex/latex/bezierplot/bezierplot-doc.tex b/macros/luatex/latex/bezierplot/bezierplot-doc.tex new file mode 100644 index 0000000000..a1239b6e57 --- /dev/null +++ b/macros/luatex/latex/bezierplot/bezierplot-doc.tex @@ -0,0 +1,200 @@ +% !TEX program = pdfLaTeX --shell-escape +\documentclass[a4paper]{article} +\usepackage{tikz,multicol,bezierplot,amsmath,cancel} +\usepackage[margin=3.5cm,top=1.75cm]{geometry} +\usepackage{fetamont} +\title{bezierplot}\author{Linus Romer} +\DeclareDocumentCommand{\graphcomparison}{ m m }{ + \begin{center} + \begin{tikzpicture}[scale=.4] + \draw (0,-5) node[below]{\tiny\texttt{\detokenize{#1}\quad | \detokenize{#2}}}; + \draw[step=1,thin] (-5,-5) grid (5,5); + \draw[thick,->] (-5,0) -- (5.5,0) node[below]{$x$}; + \draw[thick,->] (0,-5) -- (0,5.5) node[left]{$y$}; + \foreach \x in {-4,-3,-2,-1,1,2,3,4} {\draw (\x,1pt) -- (\x,-1pt) node[below]{\tiny \x};} + \foreach \y in {-4,-3,-2,-1,1,2,3,4} {\draw (1pt,\y) -- (-1pt,\y) node[left]{\tiny \y};} + \draw[color=red,domain=-5:5,range=-5:5,samples=1000] plot function{#2}; + \draw \bezierplot{#1}; + \end{tikzpicture} + \end{center} +} +\begin{document} +\maketitle\noindent +\section{Introduction} +\texttt{bezierplot} is a Lua program as well as a (Lua)\LaTeX{} package. This document describes both. + +Given a smooth function, \texttt{bezierplot} returns a smooth bezier path written in Ti\emph{k}Z notation (which also matches \MP{}) that approximates the graph of the function. For polynomial functions of degree $\leq 3$ and inverses of them, the approximation is exact (up to numeric precision). \texttt{bezierplot} finds special graph points such as extreme points and inflection points and reduces the number of used points. + +The following example will show a comparison of \textsc{gnuplot} with \verb|bezierplot| for the function $y=\sqrt{x}$ for $0\leq x \leq 5$: +\begin{center} + \begin{tikzpicture}[scale=1.4] + \draw (0,0) .. controls (0,0.7454) and (1.6667,1.4907) .. (5,2.2361); + \draw (0,0) circle(.02) -- (0,0.745) circle( .02); + \draw (1.6667,1.4907) circle(.02) -- (5,2.2361) circle( .02); + \draw (2.5,.5) node[above]{\verb|bezierplot|}; + \begin{scope}[shift={(5.2,0)}] + \draw[domain=0:5,samples=51] plot function{x**0.5}; + \foreach \x in {0,0.1,...,5.05} {\draw (\x,{\x^0.5}) circle (0.02);} + \draw (2.5,.5) node[above]{\textsc{gnuplot}}; + \end{scope} + \end{tikzpicture} +\end{center} +\textsc{gnuplot} used 51 samples (no smoothing) and is still quite inexact at the beginning, whereas \verb|bezierplot| uses 4 points only and is exact (up to numeric precision)! +\section{Installation} +As \texttt{bezierplot} is written in Lua, the installation depends whether you are using Lua\LaTeX{} or another \LaTeX{} engine. +\subsection{Installation For Lua\LaTeX{}} +If you have installed \texttt{bezierplot} by a package manager, the installation is already complete. The manual installation of \texttt{bezierplot} is done in 2 steps: +\begin{itemize} + \item copy the files \texttt{bezierplot.lua} and \texttt{bezierplot.sty} somewhere in your \texttt{texmf} tree (e.g. to \verb|~/texmf/tex/lualatex/bezierplot/bezierplot.sty| and\\ + \verb|~/texmf/scripts/bezierplot/bezierplot.lua|) + \item update the ls-R databases by running \texttt{mktexlsr} +\end{itemize} +\subsection{Additional Installation Steps For Other \LaTeX{} Engines} +You will have to call \texttt{bezierplot} as an external program via the option \texttt{--shell-escape} (\texttt{--write18} for MiK\TeX{}). Therefore, \texttt{bezierplot.lua} has to be copied with the name \texttt{bezierplot} to a place, where your OS can find it. Under Linux this usually means copying to the directory \texttt{/usr/local/bin/}, but for Windows this will probably include more steps (like adding to the \texttt{PATH}). Of course, Lua has to be installed as well. As soon as you can call \texttt{bezierplot} from a command line (e.g. by typing \verb|bezierplot "x^2"|), it should also work with other \LaTeX{} engines. +\section{Loading} +The \texttt{bezierplot} package is loaded with \verb|\usepackage{bezierplot}|. There are no loading options for the package. +\section{Usage} +\begin{multicols}{2} +\noindent A minimal example of Lua\LaTeX{} document could be: +\begin{verbatim} +\documentclass{article} +\usepackage{tikz,bezierplot} +\begin{document} +\tikz \draw \bezierplot{x^2}; +\end{document} +\end{verbatim} +\begin{center} + \tikz \draw[scale=.7] \bezierplot{x^2}; +\end{center} +\end{multicols} +\noindent +The command \verb|\bezierplot| has 6 optional arguments in the sense of +\begin{center} + \verb|\bezierplot[XMIN][XMAX][YMIN][YMAX][SAMPLES]{FUNCTION}| +\end{center} +The defaults are \verb|XMIN| = \verb|YMIN| $= -5$, \verb|XMAX| = \verb|YMAX| $= 5$ and \verb|SAMPLES| $= 0$ (this will set as few samples as possible). +\begin{center} + \begin{tikzpicture}[scale=.7] + \draw \bezierplot[-1][2]{x^2}; + \draw (0,0) node[below]{\verb|\bezierplot[-1][2]{x^2}|}; + \begin{scope}[shift={(10,0)}] + \draw \bezierplot[-1][2][0.5][3]{x^2}; + \draw (0,0) node[below]{\verb|\bezierplot[-1][2][0.5][3]{x^2}|}; + \end{scope} + \end{tikzpicture} +\end{center} +You may reverse the graph by making \verb|XMIN| bigger than \verb|XMAX|. E.g. +\begin{verbatim} + \bezierplot[-5][5]{0.5*x+1} +\end{verbatim} +returns \verb|(-5,-1.5) -- (5,3.5)|, whereas +\begin{verbatim} + \bezierplot[5][-5]{0.5*x+1} +\end{verbatim} +returns the reversed path \verb|(5,3.5) -- (-5,-1.5)|. This is useful, if you want to cycle a path to a closed area: +\begin{multicols}{2} +\begin{verbatim} +\begin{tikzpicture} + \fill[black!30] \bezierplot[-1][1]{2-x^2} + -- \bezierplot[1][-1]{x^3-x} -- cycle; + \draw \bezierplot[-1.1][1.1]{2-x^2}; + \draw \bezierplot[-1.1][1.1]{x^3-x}; +\end{tikzpicture} +\end{verbatim} +\begin{center} + \begin{tikzpicture} + \fill[black!30] \bezierplot[-1][1]{2-x^2} -- \bezierplot[1][-1]{x^3-x} -- cycle; + \draw \bezierplot[-1.1][1.1]{2-x^2}; + \draw \bezierplot[-1.1][1.1]{x^3-x}; + \end{tikzpicture} +\end{center} +\end{multicols} +\subsection{Running Raw \texttt{bezierplot}} +Of course, you can run \verb|bezierplot.lua| in a terminal without using \LaTeX{}, e.g. +\begin{verbatim} +lua bezierplot.lua "3*x^0.8+2" +\end{verbatim} +will return +\begin{verbatim} +(0,2) .. controls (0.03,2.282) and (0.268,3.244) .. (1,5) +\end{verbatim} +You can set the window of the graph and the number of samples as follows: +\begin{verbatim} +lua bezierplot.lua "FUNCTION" XMIN XMAX YMIN YMAX SAMPLES +\end{verbatim} +e.g. +\begin{verbatim} +lua bezierplot.lua "FUNCTION" 0 1 -3 2.5 201 +\end{verbatim} +will set $0\leq x\leq 1$ and $-3\leq y\leq 2.5$ and $201$ equidistant samples. You may also omit the $y$--range, hence +\begin{verbatim} +lua bezierplot.lua "FUNCTION" 0 1 +\end{verbatim} +will set $0\leq x\leq 1$ and leave the default $-5\leq y\leq 5$. The variables \verb|XMIN|, \verb|XMAX|, \verb|YMIN| and \verb|YMAX| may also be computable expressions like \verb|2*pi+6|: +\begin{verbatim} +lua bezierplot.lua "sin(x)" -pi pi +\end{verbatim} +You may use \verb|huge| for $\infty$: +\begin{verbatim} +lua bezierplot "1/x" 0 1 0 huge +\end{verbatim} +As \verb|huge| is very huge and \verb|bezierplot| uses recursive calls for nontrivial functions and non--fixed samples, this can last very long: +\begin{verbatim} +lua bezierplot "1/x" -5 5 -huge huge +\end{verbatim} +But if you set fixed samples, it will be fast again (as this does not use recursive calls): +\begin{verbatim} +lua bezierplot "1/x" -5 5 -huge huge 100 +\end{verbatim} +\subsection{Notation Of Functions} +The function term given to \verb|bezierplot| must contain at most one variable: $x$. E.g. \verb|"2.3*(x-1)^2-3"|. You must not omit \verb|*| operators: +\begin{center} + wrong:\quad $\cancel{\texttt{2x(x+1)}}$ \hfil correct:\quad \texttt{2*x*(x+1)} +\end{center} +You have two possibilities to write powers: \verb|"x^2"| and \verb|"x**2"| both mean $x^2$. + +\medskip + +The following functions and constants are possible: +\begin{center} +\begin{tabular}{ll} + \verb|abs| & absolute value (remember: your function should still be smooth)\\ + \verb|acos| & $\cos^{-1}$ inverse function of cosine in radians\\ + \verb|asin| & $\sin^{-1}$ inverse function of sine in radians\\ + \verb|atan| & $\tan^{-1}$ inverse function of tangent in radians\\ + \verb|cbrt| & cube root $\sqrt[3]{\quad}$ that works for negative numbers, too\\ + \verb|cos| & cosine for angles in radians\\ + \verb|exp| & the exponential function $e^{(\;)}$\\ + \verb|huge| & the numerical $\infty$\\ + \verb|e| & the euler constant $e=\mathrm{exp}(1)$\\ + \verb|log| & the natural logarithm $\mathrm{log}_e(\;)$\\ + \verb|pi| & Archimedes’ constant $\pi\approx 3.14$\\ + \verb|sgn| & sign function\\ + \verb|sin| & sine for angles in radians\\ + \verb|sqrt| & square root $\sqrt{\quad}$\\ + \verb|tan| & tangent for angles in radians +\end{tabular} +\end{center} +% +\newpage +% +\section{Examples of \texttt{bezierplot} in Comparison with \textsc{gnuplot}} +The following graphs are drawn with \texttt{bezierplot} (black) and \textsc{gnuplot} (red). You may not recognize the red behind the black unless you zoom in. \textsc{gnuplot} used 1000 samples per example. The functions are given below the pictures (left: bezierplot, right: \textsc{gnuplot}). +\begin{multicols}{3} +\graphcomparison{0.32*x-0.7}{0.32*x-0.7} +\graphcomparison{-x^2+4}{-x**2+4} +\graphcomparison{(x+1)*x*(x-1)}{(x+1)*x*(x-1)} +\graphcomparison{x^0.5}{x**0.5} +%\graphcomparison{x^(1/3)}{x**(1/3.)} +\graphcomparison{cbrt(x)}{sgn(x)*abs(x)**(1/3.)} +\graphcomparison{x^3*(x-1)}{x**3*(x-1)} +\graphcomparison{2*cos(3*x+4)+3}{2*cos(3*x+4)+3} +\graphcomparison{tan(x)}{tan(x)} +\graphcomparison{x+0.5*sin(x)}{x+0.5*sin(x)} +%\graphcomparison{1/(x-2)+1}{1/(x-2)+1} +\graphcomparison{2*x^2/(3*x-3)}{2*x**2/(3*x-3)} +\graphcomparison{4-exp(x)}{4-exp(x)} +\graphcomparison{log(x+4)}{log(x+4)} +\end{multicols} + +\end{document} diff --git a/macros/luatex/latex/bezierplot/bezierplot.lua b/macros/luatex/latex/bezierplot/bezierplot.lua new file mode 100755 index 0000000000..382da2d74e --- /dev/null +++ b/macros/luatex/latex/bezierplot/bezierplot.lua @@ -0,0 +1,940 @@ +#!/usr/bin/env lua +-- Linus Romer, published 2018 under LPPL Version 1.3c +-- version 1.4 2019-06-18 +abs = math.abs +acos = math.acos +asin = math.asin +atan = math.atan +cos = math.cos +exp = math.exp +e = math.exp(1) +log = math.log +pi = math.pi +sin = math.sin +sqrt = math.sqrt +tan = math.tan +huge = math.huge + +-- just a helper for debugging: +local function printdifftable(t) + for i = 1,#t do + for j = 1, 6 do + if j < 5 then + io.write(t[i][j].." ") + else + if t[i][j] then + io.write("true ") + else + io.write("false ") + end + end + end + io.write("\n") + end +end + +-- cube root defined for all real numbers x +function cbrt(x) + if x < 0 then + return -(-x)^(1/3) + else + return x^(1/3) + end +end + +function sgn(x) + if x<0 then + return -1 + elseif x>0 then + return 1 + else + return 0 + end +end + +local function round(num, decimals) + local result = tonumber(string.format("%." .. (decimals or 0) .. "f", num)) + if abs(result) == 0 then + return 0 + else + return result + end +end + +-- check if a point (x,y) satisfies xmin <= x <= xmax and ymin <= < <= ymax +local function is_in_window(x,y,xmin,xmax,ymin,ymax) + if x >= xmin and x <= xmax and y >= ymin and y <= ymax then + return true + else + return false + end +end + +local function evaluate(s) + local tempfunc = assert(load("return " .. s)) + return tempfunc() +end + + +-- 5-stencil method +-- return from a graph from f in the form {{x,y},...} +-- the derivatives in form {{x,y,dy/dx,ddy/ddx,extrema,inflection},...} +local function diffgraph(func,graph,h) + local dgraph = {} + local l = #graph + if l < 4 then -- this is not worth the pain... + for i = 1, l do + table.insert(dgraph,{graph[i][1],graph[i][2],0,0,false,false}) + end + else + local yh = func(graph[1][1]-h) + local yhh = func(graph[1][1]-2*h) + if yhh > -math.huge and yhh < math.huge -- if defined at all + and yh > -math.huge and yh < math.huge then + dgraph[1] = {graph[1][1],graph[1][2], + (yhh-8*yh+8*graph[2][2]-graph[3][2])/(12*h), + (-yhh+16*yh-30*graph[1][2]+16*graph[2][2]-graph[3][2]) + /(12*h^2)} + dgraph[2] = {graph[2][1],graph[2][2], + (yh-8*graph[1][2]+8*graph[3][2]-graph[4][2])/(12*h), + (-yh+16*graph[1][2]-30*graph[2][2]+16*graph[3][2]-graph[4][2]) + /(12*h^2)} + else -- take neighbour values + dgraph[1] = {graph[1][1],graph[1][2], + (graph[1][2]-8*graph[2][2]+8*graph[4][2]-graph[5][2])/(12*h), + (-graph[1][2]+16*graph[2][2]-30*graph[3][2] + +16*graph[4][2]-graph[5][2])/(12*h^2)} + dgraph[2] = {graph[2][1],graph[2][2], + (graph[1][2]-8*graph[2][2]+8*graph[4][2]-graph[5][2])/(12*h), + (-graph[1][2]+16*graph[2][2]-30*graph[3][2] + +16*graph[4][2]-graph[5][2])/(12*h^2)} + end + for i = 3, l-2 do + table.insert(dgraph,{graph[i][1],graph[i][2], + (graph[i-2][2]-8*graph[i-1][2]+8*graph[i+1][2]-graph[i+2][2]) + /(12*h), + (-graph[i-2][2]+16*graph[i-1][2]-30*graph[i][2] + +16*graph[i+1][2]-graph[i+2][2]) + /(12*h^2)}) + end + yh = func(graph[l][1]+h) + yhh = func(graph[l][1]+2*h) + if yhh > -math.huge and yhh < math.huge -- if defined at all + and yh > -math.huge and yh < math.huge then + dgraph[l-1] = {graph[l-1][1],graph[l-1][2], + (graph[l-3][2]-8*graph[l-2][2]+8*graph[l][2]-yh)/(12*h), + (-graph[l-3][2]+16*graph[l-2][2]-30*graph[l-1][2] + +16*graph[l][2]-yh)/(12*h^2)} + dgraph[l] = {graph[l][1],graph[l][2], + (graph[l-2][2]-8*graph[l-1][2]+8*yh-yhh)/(12*h), + (-graph[l-2][2]+16*graph[l-1][2]-30*graph[l][2] + +16*yh-yhh)/(12*h^2)} + else + -- take neighbour values + dgraph[l] = {graph[l][1],graph[l][2], + (graph[l-4][2]-8*graph[l-3][2]+8*graph[l-1][2]-graph[l][2]) + /(12*h), + (-graph[l-4][2]+16*graph[l-3][2]-30*graph[l-2][2] + +16*graph[l-1][2]-graph[l][2])/(12*h^2)} + dgraph[l-1] = {graph[l-1][1],graph[l-2][2], + (graph[l-4][2]-8*graph[l-3][2]+8*graph[l-1][2]-graph[l][2]) + /(12*h), + (-graph[l-4][2]+16*graph[l-3][2]-30*graph[l-2][2] + +16*graph[l-1][2]-graph[l][2])/(12*h^2)} + end + -- add information about being extremum / inflection point (true/false) + for i = 1, l do + dgraph[i][5] = false -- dy/dx == 0 ? default, may change later + dgraph[i][6] = false -- ddy/ddx == 0 ? default, may change later + end + for i = 1, l-1 do + -- if no gap is inbetween + if not (dgraph[i+1][1] - dgraph[i][1] > 1.5*h) then + -- check for dy/dx == 0 + -- if not already determined as near dy/dx=0 + if dgraph[i][3] == 0 then + if dgraph[i+1][3] == 0 then --take the later + dgraph[i+1][5] = true + dgraph[i][5] = false + else + dgraph[i][5] = true + end + elseif abs(dgraph[i][3]*dgraph[i+1][3]) + ~= dgraph[i][3]*dgraph[i+1][3] then -- this must be near + if abs(dgraph[i][4]) <= abs(dgraph[i+1][4]) then + dgraph[i][5] = true + else + dgraph[i+1][5] = true + end + end + -- check for ddy/ddx == 0 + -- if not already determined as near ddy/ddx=0 + if (not dgraph[i][6]) and (abs(dgraph[i][4]*dgraph[i+1][4]) + ~= dgraph[i][4]*dgraph[i+1][4]) then -- this must be near + if abs(dgraph[i][4]) <= abs(dgraph[i+1][4]) then + dgraph[i][6] = true + else + dgraph[i+1][6] = true + end + end + end + end + end + return dgraph +end + +-- simplified diffgraph function, the function is derived only once +-- return from a graph from f in the form {{x,y},...} +-- the derivatives in form {{x,y,dy/dx},...} +-- we start with index 1 and then always jump indexjump to the next +-- index +local function diffgraphsimple(func,graph,h,indexjump) + local dgraph = {} + local l = math.floor((#graph-1)/indexjump)*indexjump+1 + if l < 2 then -- this is not worth the pain... + dgraph = {graph[1][1],graph[1][2],0} + else + local yh = func(graph[1][1]-h) + local yhh = func(graph[1][1]-2*h) + if yhh > -math.huge and yhh < math.huge -- if defined at all + and yh > -math.huge and yh < math.huge then + dgraph[1] = {graph[1][1],graph[1][2], + (yhh-8*yh+8*graph[2][2]-graph[3][2])/(12*h)} + else -- take neighbour values + dgraph[1] = {graph[1][1],graph[1][2], + (graph[1][2]-8*graph[2][2]+8*graph[4][2]-graph[5][2]) + /(12*h)} + end + for i = 1+indexjump, l-1, indexjump do + table.insert(dgraph,{graph[i][1],graph[i][2], + (graph[i-2][2]-8*graph[i-1][2]+8*graph[i+1][2]-graph[i+2][2]) + /(12*h)}) + end + yh = func(graph[l][1]+h) + yhh = func(graph[l][1]+2*h) + if yhh > -math.huge and yhh < math.huge -- if defined at all + and yh > -math.huge and yh < math.huge then + + table.insert(dgraph,{graph[l][1],graph[l][2], + (graph[l-2][2]-8*graph[l-1][2]+8*yh-yhh)/(12*h)}) + else + -- take neighbour values + table.insert(dgraph,{graph[l][1],graph[l][2], + (graph[l-4][2]-8*graph[l-3][2]+8*graph[l-1][2]-graph[l][2]) + /(12*h)}) + end + end + return dgraph +end + +-- diffgraph for cubic function ax^3+bx^2+cx+d +-- returns the derivatives in form {{x,y,dy/dx,ddy/ddx},...} +-- if isinverse = true then the coordinates will be inversed +local function diffgraphcubic(graph,a,b,c,d,isinverse) + local dgraph = {} + local l = #graph + for i = 1, l do + if isinverse then + dgraph[#dgraph+1] = {graph[i][2],graph[i][1],c + +graph[i][2]*(2*b+3*a*graph[i][2]),6*a*graph[i][2]+2*b} + else + dgraph[#dgraph+1] = {graph[i][1],graph[i][2],c + +graph[i][1]*(2*b+3*a*graph[i][1]),6*a*graph[i][1]+2*b} + end + end + return dgraph +end + +-- checks for 100 x, if the function given by funcstring +-- fits the graph g (up to maxerror) after filling in +-- the parameters a, b, c, d +-- if the graph is inverted, then isinverse has to be set true +local function do_parameters_fit(a,b,c,d,funcstring,funcgraph,maxerror,isinverse) + if not (a > -math.huge and a < math.huge and b > -math.huge and b < math.huge and + c > -math.huge and c < math.huge and d > -math.huge and d < math.huge) then + return false + end + local funcx = string.gsub(funcstring, "a", a) + local funcx = string.gsub(funcx, "b", b) + local funcx = string.gsub(funcx, "c", c) + local funcx = string.gsub(funcx, "d", d) + local func = assert(load("local x = ...; return "..funcx)) + for i = 1, #funcgraph, math.max(1,math.floor(0.01*#funcgraph)) do + if isinverse then + if abs(func(funcgraph[i][2])-funcgraph[i][1]) + > maxerror then + return false + end + else + if abs(func(funcgraph[i][1])-funcgraph[i][2]) + > maxerror then + return false + end + end + end + return true +end + +-- f(x)=a*x^3+b*x+c +local function parameters_cubic(xp,yp,xq,yq,xr,yr,xs,ys) + local a = (((xp^2 * xq) * yr) - ((xp^2 * xq) * ys) + - ((xp^2 * xr) * yq) + ((xp^2 * xr) * ys) + ((xp^2 * xs) * yq) + - ((xp^2 * xs) * yr) - ((xp * xq^2) * yr) + ((xp * xq^2) * ys) + + ((xp * xr^2) * yq) - ((xp * xr^2) * ys) - ((xp * xs^2) * yq) + + ((xp * xs^2) * yr) + ((xq^2 * xr) * yp) - ((xq^2 * xr) * ys) + - ((xq^2 * xs) * yp) + ((xq^2 * xs) * yr) - ((xq * xr^2) * yp) + + ((xq * xr^2) * ys) + ((xq * xs^2) * yp) - ((xq * xs^2) * yr) + + ((xr^2 * xs) * yp) - ((xr^2 * xs) * yq) - ((xr * xs^2) * yp) + + ((xr * xs^2) * yq)) / + (((xp^3 * xq^2) * xr) - ((xp^3 * xq^2) * xs) + - ((xp^3 * xq) * xr^2) + ((xp^3 * xq) * xs^2) + + ((xp^3 * xr^2) * xs) - ((xp^3 * xr) * xs^2) + - ((xp^2 * xq^3) * xr) + ((xp^2 * xq^3) * xs) + + ((xp^2 * xq) * xr^3) - ((xp^2 * xq) * xs^3) + - ((xp^2 * xr^3) * xs) + ((xp^2 * xr) * xs^3) + + ((xp * xq^3) * xr^2) - ((xp * xq^3) * xs^2) + - ((xp * xq^2) * xr^3) + ((xp * xq^2) * xs^3) + + ((xp * xr^3) * xs^2) - ((xp * xr^2) * xs^3) + - ((xq^3 * xr^2) * xs) + ((xq^3 * xr) * xs^2) + + ((xq^2 * xr^3) * xs) - ((xq^2 * xr) * xs^3) + - ((xq * xr^3) * xs^2) + ((xq * xr^2) * xs^3)) + local b = ((((-xp^3) * xq) * yr) + ((xp^3 * xq) * ys) + + ((xp^3 * xr) * yq) - ((xp^3 * xr) * ys) - ((xp^3 * xs) * yq) + + ((xp^3 * xs) * yr) + ((xp * xq^3) * yr) - ((xp * xq^3) * ys) + - ((xp * xr^3) * yq) + ((xp * xr^3) * ys) + ((xp * xs^3) * yq) + - ((xp * xs^3) * yr) - ((xq^3 * xr) * yp) + ((xq^3 * xr) * ys) + + ((xq^3 * xs) * yp) - ((xq^3 * xs) * yr) + ((xq * xr^3) * yp) + - ((xq * xr^3) * ys) - ((xq * xs^3) * yp) + ((xq * xs^3) * yr) + - ((xr^3 * xs) * yp) + ((xr^3 * xs) * yq) + ((xr * xs^3) * yp) + - ((xr * xs^3) * yq)) / + (((xp^3 * xq^2) * xr) - ((xp^3 * xq^2) * xs) + - ((xp^3 * xq) * xr^2) + ((xp^3 * xq) * xs^2) + + ((xp^3 * xr^2) * xs) - ((xp^3 * xr) * xs^2) + - ((xp^2 * xq^3) * xr) + ((xp^2 * xq^3) * xs) + + ((xp^2 * xq) * xr^3) - ((xp^2 * xq) * xs^3) + - ((xp^2 * xr^3) * xs) + ((xp^2 * xr) * xs^3) + + ((xp * xq^3) * xr^2) - ((xp * xq^3) * xs^2) + - ((xp * xq^2) * xr^3) + ((xp * xq^2) * xs^3) + + ((xp * xr^3) * xs^2) - ((xp * xr^2) * xs^3) + - ((xq^3 * xr^2) * xs) + ((xq^3 * xr) * xs^2) + + ((xq^2 * xr^3) * xs) - ((xq^2 * xr) * xs^3) + - ((xq * xr^3) * xs^2) + ((xq * xr^2) * xs^3)) + local c = (((xp^3 * xq^2) * yr) - ((xp^3 * xq^2) * ys) + - ((xp^3 * xr^2) * yq) + ((xp^3 * xr^2) * ys) + + ((xp^3 * xs^2) * yq) - ((xp^3 * xs^2) * yr) + - ((xp^2 * xq^3) * yr) + ((xp^2 * xq^3) * ys) + + ((xp^2 * xr^3) * yq) - ((xp^2 * xr^3) * ys) + - ((xp^2 * xs^3) * yq) + ((xp^2 * xs^3) * yr) + + ((xq^3 * xr^2) * yp) - ((xq^3 * xr^2) * ys) + - ((xq^3 * xs^2) * yp) + ((xq^3 * xs^2) * yr) + - ((xq^2 * xr^3) * yp) + ((xq^2 * xr^3) * ys) + + ((xq^2 * xs^3) * yp) - ((xq^2 * xs^3) * yr) + + ((xr^3 * xs^2) * yp) - ((xr^3 * xs^2) * yq) + - ((xr^2 * xs^3) * yp) + ((xr^2 * xs^3) * yq)) / + (((xp^3 * xq^2) * xr) - ((xp^3 * xq^2) * xs) + - ((xp^3 * xq) * xr^2) + ((xp^3 * xq) * xs^2) + + ((xp^3 * xr^2) * xs) - ((xp^3 * xr) * xs^2) + - ((xp^2 * xq^3) * xr) + ((xp^2 * xq^3) * xs) + + ((xp^2 * xq) * xr^3) - ((xp^2 * xq) * xs^3) + - ((xp^2 * xr^3) * xs) + ((xp^2 * xr) * xs^3) + + ((xp * xq^3) * xr^2) - ((xp * xq^3) * xs^2) + - ((xp * xq^2) * xr^3) + ((xp * xq^2) * xs^3) + + ((xp * xr^3) * xs^2) - ((xp * xr^2) * xs^3) + - ((xq^3 * xr^2) * xs) + ((xq^3 * xr) * xs^2) + + ((xq^2 * xr^3) * xs) - ((xq^2 * xr) * xs^3) + - ((xq * xr^3) * xs^2) + ((xq * xr^2) * xs^3)) + local d = ((((xp^(3) * xq^(2)) * xr) * ys) + - (((xp^(3) * xq^(2)) * xs) * yr) - (((xp^(3) * xq) * xr^(2)) * ys) + + (((xp^(3) * xq) * xs^(2)) * yr) + (((xp^(3) * xr^(2)) * xs) * yq) + - (((xp^(3) * xr) * xs^(2)) * yq) - (((xp^(2) * xq^(3)) * xr) * ys) + + (((xp^(2) * xq^(3)) * xs) * yr) + (((xp^(2) * xq) * xr^(3)) * ys) + - (((xp^(2) * xq) * xs^(3)) * yr) - (((xp^(2) * xr^(3)) * xs) * yq) + + (((xp^(2) * xr) * xs^(3)) * yq) + (((xp * xq^(3)) * xr^(2)) * ys) + - (((xp * xq^(3)) * xs^(2)) * yr) - (((xp * xq^(2)) * xr^(3)) * ys) + + (((xp * xq^(2)) * xs^(3)) * yr) + (((xp * xr^(3)) * xs^(2)) * yq) + - (((xp * xr^(2)) * xs^(3)) * yq) - (((xq^(3) * xr^(2)) * xs) * yp) + + (((xq^(3) * xr) * xs^(2)) * yp) + (((xq^(2) * xr^(3)) * xs) * yp) + - (((xq^(2) * xr) * xs^(3)) * yp) - (((xq * xr^(3)) * xs^(2)) * yp) + + (((xq * xr^(2)) * xs^(3)) * yp)) / + (((xp^(3) * xq^(2)) * xr) - + ((xp^(3) * xq^(2)) * xs) - ((xp^(3) * xq) * xr^(2)) + + ((xp^(3) * xq) * xs^(2)) + ((xp^(3) * xr^(2)) * xs) + - ((xp^(3) * xr) * xs^(2)) - ((xp^(2) * xq^(3)) * xr) + + ((xp^(2) * xq^(3)) * xs) + ((xp^(2) * xq) * xr^(3)) + - ((xp^(2) * xq) * xs^(3)) - ((xp^(2) * xr^(3)) * xs) + + ((xp^(2) * xr) * xs^(3)) + ((xp * xq^(3)) * xr^(2)) + - ((xp * xq^(3)) * xs^(2)) - ((xp * xq^(2)) * xr^(3)) + + ((xp * xq^(2)) * xs^(3)) + ((xp * xr^(3)) * xs^(2)) + - ((xp * xr^(2)) * xs^(3)) - ((xq^(3) * xr^(2)) * xs) + + ((xq^(3) * xr) * xs^(2)) + ((xq^(2) * xr^(3)) * xs) + - ((xq^(2) * xr) * xs^(3)) - ((xq * xr^(3)) * xs^(2)) + + ((xq * xr^(2)) * xs^(3))) + return a, b, c, d +end + +-- f(x)=a*x+b +local function parameters_affine(xp,yp,xq,yq) + local a = (yp - yq) / (xp - xq) + local b = ((xp * yq) - (xq * yp)) / (xp - xq) + return a, b +end + +-- what is the sum of the squared error +-- when comparing the bezier path +-- p.. control q and r .. s +-- with the graph g from index starti to endi +-- (looking at the points at roughly t=.33 and t=.67) +local function squareerror(f,g,starti,endi,qx,qy,rx,ry) + local result = 0 + for t = .1, .9, .1 do + x = (1-t)^3*g[starti][1]+3*t*(1-t)^2*qx+3*t^2*(1-t)*rx+t^3*g[endi][1] + y = (1-t)^3*g[starti][2]+3*t*(1-t)^2*qy+3*t^2*(1-t)*ry+t^3*g[endi][2] + result = result + (y-f(x))^2 + end + return result +end + +-- converts a table with bezier point information +-- to a string with rounded values +-- the path is reversed, if rev is true +-- e.g. if bezierpoints = {{0,1},{2,3,4,5,6,7},{8,9,10,11,12,13}} +-- then +-- (0,1) .. controls (2,3) and (4,5) .. (6,7) .. controls +-- (8,9) and (10,11) .. (12,13) +-- will be returned +-- the notation "pgfplots" will change the notation to +-- YES: \addplot coordinates {(0,1) (6,7) (2,3) (4,5) (6,7) (12,13) (8,9) (10,11)} +-- NO: 0 1 \\ 6 7 \\ 2 3 \\ 4 5 \\ \\ 6 7 \\ 12 13 \\ 8 9 \\ 10 11 \\ +-- As pgfplots does not connect the bezier segments +-- reverse paths are not implemented +local function beziertabletostring(bezierpoints,rndx,rndy,rev,notation) + local bezierstring = "" + local b = {{round(bezierpoints[1][1],rndx),round(bezierpoints[1][2],rndy)}} -- rounded and then + -- reduced points (if identical after rounding) + -- rounding + for i = 2, #bezierpoints do + -- check if x--coordinates are identical + if round(bezierpoints[i][#bezierpoints[i]-1],rndx) ~= b[#b][#b[#b]-1] then + b[#b+1] = {} + for j = 1, #bezierpoints[i] do + if j % 2 == 0 then -- x coordinate + b[#b][j] = round(bezierpoints[i][j],rndx) + else + b[#b][j] = round(bezierpoints[i][j],rndy) + end + end + end + end + if #b > 1 then -- if not empty or single point + -- check if bezierstring contains only straight lines + local onlystraightlines = true + for i = 1, #b do + if #b[i] > 2 then + onlystraightlines = false + break + end + end + if onlystraightlines then + if rev then + bezierstring = "(" .. b[#b][1] .. "," .. b[#b][2] ..")" + for i = #b-1, 1, -1 do + bezierstring = bezierstring .. + " -- (" .. b[i][1] .. "," .. b[i][2] ..")" + end + else + if notation == "pgfplots" then + bezierstring = "\\addplot coordinates {(" + .. b[1][1] .. "," .. b[1][2] .. ") (" + .. b[2][1] .. "," .. b[2][2] .. ") (" + .. b[1][1] .. "," .. b[1][2] .. ") (" + .. b[2][1] .. "," .. b[2][2] .. ") }" + else -- notation = tikz + bezierstring = "(" .. b[1][1] .. "," .. b[1][2] ..")" + for i = 2, #b do + bezierstring = bezierstring .. + " -- (" .. b[i][1] .. "," .. b[i][2] ..")" + end + end + end + else + if rev then + bezierstring = "(" .. b[#b][#b[#b]-1] .. "," + .. b[#b][#b[#b]] ..")" -- initial point + for i = #b, 2, -1 do + if #b[i] >= 6 then -- cubic bezier spline + bezierstring = bezierstring .. " .. controls (" + .. b[i][3] .. "," .. b[i][4] ..") and (" + .. b[i][1] .. "," .. b[i][2] .. ") .. (" + .. b[i-1][#b[i-1]-1] .. "," .. b[i-1][#b[i-1]]..")" + else + bezierstring = bezierstring .. " (" + .. b[i-1][#b[i-1]-1] .. "," .. b[i-1][#b[i-1]] ..")" + end + end + else + if notation == "pgfplots" then + bezierstring = "\\addplot coordinates {" + for i = 1, #b-1 do + if #b[i+1] >= 6 then -- cubic bezier spline + bezierstring = bezierstring .. "(" + .. b[i][#b[i]-1] .. "," .. b[i][#b[i]] .. ") (" + .. b[i+1][5] .. "," .. b[i+1][6] .. ") (" + .. b[i+1][1] .. "," .. b[i+1][2] .. ") (" + .. b[i+1][3] .. "," .. b[i+1][4] .. ") " + end + end + bezierstring = bezierstring .. "}" + else -- notation = tikz + bezierstring = "(" .. b[1][1] .. "," + .. b[1][2] ..")" -- initial point + for i = 2, #b do + if #b[i] >= 6 then -- cubic bezier spline + bezierstring = bezierstring .. " .. controls (" + .. b[i][1] .. "," .. b[i][2] ..") and (" + .. b[i][3] .. "," .. b[i][4] .. ") .. (" + .. b[i][5] .. "," .. b[i][6]..")" + else + bezierstring = bezierstring .. " (" + .. b[i][1] .. "," .. b[i][2] ..")" + end + end + end + end + end + end + return bezierstring +end + +-- take end points of a graph g of the function f +-- (from indices starti to endi) +-- without extrema or inflection points inbetween +-- and try to approximate it with a cubic bezier curve +-- (round to rndx and rndy when printing) +-- if maxerror <= 0, the function will not be recursive anymore +local function graphtobezierapprox(f,g,starti,endi,maxerror) + local px = g[starti][1] + local py = g[starti][2] + local dp = g[starti][3] + local sx = g[endi][1] + local sy = g[endi][2] + local ds = g[endi][3] + -- we compute the corner point c, where the controls would meet + local cx = ((dp * px) - (ds * sx) - py + sy) / (dp - ds) + local cy = (dp * ((ds * px) - (ds * sx) - py + sy) / (dp - ds)) + py + -- now we slide q between p and c & r between s and c + -- and search for the best qx and best rx + local qx = px+.01*(cx-px) + local qy = py+.01*(cy-py) + local rx = sx+.01*(cx-sx) + local ry = sy+.01*(cy-sy) + local err = squareerror(f,g,starti,endi,qx,qy,rx,ry) + for i = 2, 99 do + for j = 2, 99 do + xa = px+i*.01*(cx-px) + ya = py+i*.01*(cy-py) + xb = sx+j*.01*(cx-sx) + yb = sy+j*.01*(cy-sy) + -- now check, if xa and xb fit better + -- than the last qx and rx did + -- (sum of squares must be smaller) + local newerror = squareerror(f,g,starti,endi,xa,ya,xb,yb) + if newerror < err then + qx = xa + qy = ya + rx = xb + ry = yb + err = newerror + end + end + end + if maxerror > 0 then + -- check if it is close enough: (recycling err, xa, ya) + err = 0 + for t = .1, .9, .1 do + xa = (1-t)^3*g[starti][1]+3*t*(1-t)^2*qx+3*t^2*(1-t)*rx+t^3*g[endi][1] + ya = (1-t)^3*g[starti][2]+3*t*(1-t)^2*qy+3*t^2*(1-t)*ry+t^3*g[endi][2] + if abs(ya-f(xa)) > err then + err = abs(ya-f(xa)) + err = abs(ya-f(xa)) + end + end + if (err <= maxerror) + and qx > -math.huge and qx < math.huge + and qy > -math.huge and qy < math.huge + and rx > -math.huge and ry < math.huge + and sx > -math.huge and sy < math.huge + then + return {qx,qy,rx,ry,sx,sy} + else + -- search for an intermediate point where the graph has the same + -- slope as the line from the start point to the end point: + local interindex = math.floor(.5*starti+.5*endi) -- will change + for i = starti + 1, endi - 1 do + if abs(g[i][3]-(g[endi][2]-g[starti][2]) + /(g[endi][1]-g[starti][1])) + < abs(g[interindex][3]-(g[endi][2]-g[starti][2]) + /(g[endi][1]-g[starti][1])) then + interindex = i + end + end + local left = graphtobezierapprox(f,g,starti,interindex,maxerror) + local right = graphtobezierapprox(f,g,interindex,endi,maxerror) + for i=1, #right do --now append the right to the left: + left[#left+1] = right[i] + end + return left + end + elseif qx > -math.huge and qx < math.huge + and qy > -math.huge and qy < math.huge + and rx > -math.huge and ry < math.huge + and sx > -math.huge and sy < math.huge + then + return {qx,qy,rx,ry,sx,sy} + end +end + +-- like above but exact for quadratic and cubic (if not inverse) +-- resp. exact for squareroot and cuberoot (if inverse) +local function graphtobezier(g,starti,endi,isinverse) + local px = g[starti][1] + local py = g[starti][2] + local dp = g[starti][3] + local sx = g[endi][1] + local sy = g[endi][2] + local ds = g[endi][3] + local qx = px+(sx-px)/3 + local rx = px+2*(sx-px)/3 + local qy = py+(qx-px)*dp + local ry = sy+(rx-sx)*ds + if isinverse then + return {qy,qx,ry,rx,sy,sx} + else + return {qx,qy,rx,ry,sx,sy} + end +end + +-- main function +function bezierplot(functionstring,xminstring,xmaxstring,yminstring,ymaxstring,samplesstring,notation) + local fstringreplaced = string.gsub(functionstring, "%*%*", "^") + local f = assert(load("local x = ...; return " .. fstringreplaced)) + local xmin = evaluate(xminstring) + local xmax = evaluate(xmaxstring) + local ymin = evaluate(yminstring) + local ymax = evaluate(ymaxstring) + local samples = evaluate(samplesstring) + local isreverse = false + if xmin > xmax then + isreverse = true + elseif xmin == xmax then + xmax = xmin + 10 + end + xmin, xmax = math.min(xmin,xmax), math.max(xmin,xmax) + if ymin == ymax then + ymax = ymin + 10 + end + ymin, ymax = math.min(ymin,ymax), math.max(ymin,ymax) + local xsteps = 50000 + -- if samples < 2 the samples will be chosen as wisely as possible + local arbitrary_samples = true + if samples >= 2 then + arbitrary_samples = false + xsteps = (samples-1)*math.max(2,math.floor(xsteps/(samples-1))) + end + local xstep = (xmax-xmin)/xsteps + -- the output of the x coordinates will be rounded to rndx digits + local rndx = math.max(0,math.floor(5.5-log(xmax-xmin)/log(10))) + local xerror = abs(xmax-xmin)/(10^rndx) + -- the output of the y coordinates will be rounded to rndy digits + local rndy = math.max(0,math.floor(5.5-log(ymax-ymin)/log(10))) + local yerror = (ymax-ymin)/(10^rndy) + -- determine parts of the graph that are inside window + local graphs = {} -- graph split to the connected parts + local graph = {} -- graphs concatenated (needed for function type) + local outside = true -- value is outside window + local i = 0 + local j = 0 + local yminreal -- determine the real minimimum of the y coord. + local ymaxreal -- just decring + local yminrealfound = false + local ymaxrealfound = false + for n = 0, xsteps do + local x = xmin + n/xsteps*(xmax-xmin) + if n == xsteps then + x = xmax + end + local y = f(x) + if (y >= ymin-.1*yerror and ymin ~= -huge or y > ymin and ymin == -huge) + and (y <= ymax+.1*yerror and ymax ~= huge or y < ymax and ymax == huge) + then -- inside + if outside then -- if it was outside before + outside = false + j = 0 + i = i + 1 + graphs[i] = {} + end + j = j + 1 + graphs[i][j] = {x,y} + graph[#graph+1] = {x,y} + if not yminrealfound or yminrealfound and y < yminreal then + yminreal = y + yminrealfound = true + end + if not ymaxrealfound or ymaxrealfound and y > ymaxreal then + ymaxreal = y + ymaxrealfound = true + end + else + outside = true + end + end + + -- some redefinitions + if #graph ~= 0 and yminreal ~= ymaxreal then + ymin = yminreal + ymax = ymaxreal + rndy = math.max(0,math.floor(5.5-log(ymax-ymin)/log(10))) + yerror = (ymax-ymin)/(10^rndy) + end + + -- check for the function type (for this, we need the concatenated + -- parts of the graph) + -- go through the connected parts + local functiontype = "unknown" + local a, b, c, d -- possible function parameter + -- check for affine functions: + local l = #graph + a, b = parameters_affine(graph[1][1],graph[1][2], + graph[l][1],graph[l][2]) + if do_parameters_fit(a,b,0,0,"a*x+b",graph,yerror,false) then + functiontype = "affine" + else -- check for cubic functions (includes quadratic functions) + a, b, c, d = parameters_cubic(graph[1][1],graph[1][2], + graph[math.floor(l/3)][1],graph[math.floor(l/3)][2], + graph[math.floor(2*l/3)][1],graph[math.floor(2*l/3)][2], + graph[l][1],graph[l][2]) + if do_parameters_fit(a,b,c,d,"a*x^3+b*x^2+c*x+d",graph, + yerror,false) then + functiontype = "cubic" + else -- check for cuberoot functions (includes squareroots) + a, b, c, d = parameters_cubic(graph[1][2],graph[1][1], + graph[math.floor(l/3)][2],graph[math.floor(l/3)][1], + graph[math.floor(2*l/3)][2],graph[math.floor(2*l/3)][1], + graph[l][2],graph[l][1]) + if do_parameters_fit(a,b,c,d,"a*x^3+b*x^2+c*x+d",graph, + xerror,true) then + functiontype = "cuberoot" + end + end + end + + local bezierpoints = {} + -- the bezier path (0,1) .. controls + -- (2,3) and (4,5) .. (6,7) .. controls + -- (8,9) and (10,11) .. (12,13) + -- will be stored as + -- bezierpoints={{0,1},{2,3,4,5,6,7},{8,9,10,11,12,13}} + + if functiontype == "affine" then + if arbitrary_samples then + bezierpoints = {{graph[1][1],graph[1][2]},{graph[#graph][1], + graph[#graph][2]}} + else -- we can here savely assume that graphs has only one part, + -- therefore graphs[1]=graph + for i = 1, #graph, math.floor(xsteps/(samples-1)) do + bezierpoints[#bezierpoints+1] = {graph[i][1],graph[i][2]} + end + end + elseif functiontype == "cubic" then + local extrema_inflections = {} -- store the extrema and + -- inflection points for arbitrary samples + if arbitrary_samples then + if math.abs(a) < yerror*1e-10 then -- quadratic case (one extremum) + if is_in_window(-c/(2*b),(-c^2+4*b*d)/(4*b),xmin,xmax, + ymin,ymax) then + extrema_inflections = {{-c/(2*b),(-c^2+4*b*d)/(4*b)}} + end + else -- cubic case (two extrema and one inflection point) + -- we order the points with the help of sgn + -- check for first extrema + if is_in_window((-sgn(a)*sqrt(-3*a*c+b^2)-b)/(3*a), + (2*b^3+27*a^2*d-9*a*b*c+sqrt(b^2-3*a*c)*sgn(a)* + (2*b^2-6*a*c))/(27*a^2),xmin,xmax,ymin,ymax) then + extrema_inflections[#extrema_inflections+1] = + {(-sgn(a)*sqrt(-3*a*c+b^2)-b)/(3*a),(2*b^3+27*a^2*d- + 9*a*b*c+sqrt(b^2-3*a*c)*sgn(a)*(2*b^2-6*a*c))/(27*a^2)} + end + -- check for inflection point (has to be inbetween) + if is_in_window(-b/(3*a),(2*b^3+27*a^2*d-9*a*b*c) + /(27*a^2),xmin,xmax,ymin,ymax) then + extrema_inflections[#extrema_inflections+1]={-b/(3*a), + (2*b^3+27*a^2*d-9*a*b*c)/(27*a^2)} + end + -- check for second extrema + if is_in_window((sgn(a)*sqrt(-3*a*c+b^2)-b)/(3*a), + (2*b^3+27*a^2*d-9*a*b*c+sqrt(b^2-3*a*c)*sgn(a)* + (-2*b^2+6*a*c))/(27*a^2),xmin,xmax,ymin,ymax) then + extrema_inflections[#extrema_inflections+1] = + {(sgn(a)*sqrt(-3*a*c+b^2)-b)/(3*a),(2*b^3+27*a^2*d- + 9*a*b*c+sqrt(b^2-3*a*c)*sgn(a)*(-2*b^2+6*a*c))/(27*a^2)} + end + end + end + for part = 1, #graphs do + bezierpoints[#bezierpoints+1] = {graphs[part][1][1], + graphs[part][1][2]} -- initial points + local graphsamples = {}-- will be the graph reduced to the + -- samples (or the most important points) + local dg -- will be the differentiated graph + if arbitrary_samples then -- add extrema and inflection + -- points to the border points + graphsamples = {{graphs[part][1][1], + graphs[part][1][2]}} + for j = 1, #extrema_inflections do + if extrema_inflections[j][1] > math.min( + graphs[part][1][1] ,graphs[part][#graphs[part]][1]) + and extrema_inflections[j][1] < math.max( + graphs[part][1][1] ,graphs[part][#graphs[part]][1]) + then + graphsamples[#graphsamples+1] = + {extrema_inflections[j][1], + extrema_inflections[j][2]} + end + end + graphsamples[#graphsamples+1] = + {graphs[part][#graphs[part]][1], + graphs[part][#graphs[part]][2]} + else + for i = 1, #graphs[part], xsteps/(samples-1) do + graphsamples[#graphsamples+1] = + {graphs[part][i][1],graphs[part][i][2]} + end + end + dg = diffgraphcubic(graphsamples,a,b,c,d,false) + for i = 2, #dg do + bezierpoints[#bezierpoints+1] = graphtobezier(dg,i-1,i,false) + end + end + elseif functiontype == "cuberoot" then + local inflection = {} -- store the inflection point + if arbitrary_samples and math.abs(a) ~= 0 + and is_in_window((2*b^3+27*a^2*d-9*a*b*c)/(27*a^2),-b/(3*a), + xmin,xmax,ymin,ymax) then + inflection = {(2*b^3+27*a^2*d-9*a*b*c)/(27*a^2),-b/(3*a)} + end + -- (there cannot be more than one part) + bezierpoints[#bezierpoints+1] = {graphs[1][1][1], + graphs[1][1][2]} -- initial points + local graphsamples = {}-- will be the graph reduced to the + -- samples (or the most important points) + local dg -- will be the differentiated graph + if arbitrary_samples then -- add inflection point (if exis.) + graphsamples = {{graphs[1][1][1], + graphs[1][1][2]}} + if #inflection > 0 and inflection[1] > math.min( + graphs[1][1][1],graphs[1][#graphs[1]][1]) + and inflection[1] < math.max( + graphs[1][1][1],graphs[1][#graphs[1]][1]) + then + graphsamples[#graphsamples+1] = + {inflection[1],inflection[2]} + end + graphsamples[#graphsamples+1] = + {graphs[1][#graphs[1]][1], + graphs[1][#graphs[1]][2]} + else + for i = 1, #graphs[1], xsteps/(samples-1) do + graphsamples[#graphsamples+1] = + {graphs[1][i][1],graphs[1][i][2]} + end + end + dg = diffgraphcubic(graphsamples,a,b,c,d,true) + for i = 2, #dg do + bezierpoints[#bezierpoints+1] = graphtobezier(dg,i-1,i,true) + end + else + ---------- generic case (no special function) ---------------- + if arbitrary_samples then + -- go through the connected parts + for part = 1, #graphs do + local dg = diffgraph(f,graphs[part],xstep) + --printdifftable(dg) -- for debugging + bezierpoints[#bezierpoints+1] = {dg[1][1],dg[1][2]} + local startindex = 1 + for k = 2, #dg do + if dg[k][5] or dg[k][6] then -- extrema and inflection points + local tobeadded = graphtobezierapprox( + f,dg,startindex,k,10*yerror) + -- tobeadded may contain a multiple of 6 entries + -- e.g. {1,2,3,4,5,6,7,8,9,10,11,12} + for i = 1, math.floor(#tobeadded/6) do + bezierpoints[#bezierpoints+1] = {} + for j = 1, 6 do + bezierpoints[#bezierpoints][j] = tobeadded[(i-1)*6+j] + end + end + startindex = k + end + end + if startindex ~= #dg then -- if no special points inbetween + local tobeadded = graphtobezierapprox(f,dg, + startindex,#dg,10*yerror) + -- tobeadded may contain a multiple of 6 entries + -- e.g. {1,2,3,4,5,6,7,8,9,10,11,12} + for i = 1, math.floor(#tobeadded/6) do + bezierpoints[#bezierpoints+1] = {} + for j = 1, 6 do + bezierpoints[#bezierpoints][j] = tobeadded[(i-1)*6+j] + end + end + end + end + else -- fixed samples in the generic case + -- go through the connected parts + for part = 1, #graphs do + local dg = diffgraphsimple(f,graphs[part],xstep, + math.floor(0.5+xsteps/(samples-1))) + bezierpoints[#bezierpoints+1] = {dg[1][1],dg[1][2]} -- initial points + for i = 2, #dg do + bezierpoints[#bezierpoints+1] = graphtobezier(dg,i-1,i,false) + end + end + end + end + return beziertabletostring(bezierpoints,rndx,rndy,isreverse,notation) +end + +-- main program -- + +if not pcall(debug.getlocal, 4, 1) then +--if debug.getinfo(3) == nil then + if #arg >= 1 then + local xmin = -5 + local xmax = 5 + if #arg >= 2 then + xmin = arg[2] + end + if #arg >= 3 then + xmax = arg[3] + end + local ymin = -5 + local ymax = 5 + if #arg >= 4 then + ymin = arg[4] + end + if #arg >= 5 then + ymax = arg[5] + end + local samples = 0 + if #arg >= 6 then + samples = arg[6] + end + local notation = "tikz" + if #arg >= 7 then + notation = arg[7] + end + print(bezierplot(arg[1],xmin,xmax,ymin,ymax,samples,notation)) + end +end + + + diff --git a/macros/luatex/latex/bezierplot/bezierplot.sty b/macros/luatex/latex/bezierplot/bezierplot.sty new file mode 100644 index 0000000000..94be30b4dc --- /dev/null +++ b/macros/luatex/latex/bezierplot/bezierplot.sty @@ -0,0 +1,17 @@ +\NeedsTeXFormat{LaTeX2e} +\ProvidesPackage{bezierplot}[2019/06/18 bezierplot] +\RequirePackage{xparse} +\RequirePackage{iftex} +\ifLuaTeX + \directlua{require("bezierplot")} + \DeclareExpandableDocumentCommand{\xbezierplot}{O{-5} O{5} O{-5} O{5} O{0} O{tikz} m}{% + \directlua{tex.sprint(bezierplot("#7","#1","#2","#3","#4","#5","#6"))} + } +\else + \let\xpandblinpt\@@input + \DeclareExpandableDocumentCommand{\xbezierplot}{O{-5} O{5} O{-5} O{5} O{0} O{tikz} m}{% + \xpandblinpt|"bezierplot '#7' '#1' '#2' '#3' '#4' '#5' '#6'" + } +\fi +\providecommand\bezierplot{\romannumeral`\^^@\xbezierplot} +\endinput
\ No newline at end of file |