// PlotStress plots stresses along y=0 (horizontal line) and x=0 (vertical line)
func (o *PlateHole) PlotStress(t, L float64, npts int) {

	d := utl.LinSpace(o.r, L, npts)
	Sx := make([]float64, npts)
	Sy := make([]float64, npts)
	Sxy := make([]float64, npts)

	plt.Subplot(2, 1, 1)
	for i := 0; i < npts; i++ {
		Sx[i], Sy[i], _, Sxy[i] = o.Stress(t, []float64{d[i], 0}) // y=0
	}
	plt.Plot(d, Sx, "color='r',label='$\\sigma_x$ @ $y=0$'")
	plt.Plot(d, Sy, "color='g',label='$\\sigma_y$ @ $y=0$'")
	plt.Plot(d, Sxy, "color='b',label='$\\sigma_{xy}$ @ $y=0$'")
	plt.Gll("$x$", "stresses", "")

	plt.Subplot(2, 1, 2)
	for i := 0; i < npts; i++ {
		Sx[i], Sy[i], _, Sxy[i] = o.Stress(t, []float64{0, d[i]}) // x=0
	}
	plt.Plot(Sx, d, "color='r',label='$\\sigma_x$ @ $x=0$'")
	plt.Plot(Sy, d, "color='g',label='$\\sigma_y$ @ $x=0$'")
	plt.Plot(Sxy, d, "color='b',label='$\\sigma_{xy}$ @ $x=0$'")
	plt.Gll("stresses", "$y$", "")
}

// PlotTwoNurbs plots two NURBS for comparison
func PlotTwoNurbs(dirout, fn string, b, c *Nurbs, npts int, ids bool, extra func()) {
	plt.Reset()
	if io.FnExt(fn) == ".eps" {
		plt.SetForEps(1.5, 500)
	} else {
		plt.SetForPng(1.5, 500, 150)
	}

	plt.Subplot(3, 1, 1)
	b.DrawCtrl2d(ids, "", "")
	b.DrawElems2d(npts, ids, "", "")
	if extra != nil {
		extra()
	}
	plt.Equal()

	plt.Subplot(3, 1, 2)
	c.DrawCtrl2d(ids, "", "")
	c.DrawElems2d(npts, ids, "", "")
	plt.Equal()

	plt.Subplot(3, 1, 3)
	b.DrawElems2d(npts, ids, ", lw=3", "")
	c.DrawElems2d(npts, ids, ", color='red', marker='+', markevery=10", "color='green', size=7, va='bottom'")
	plt.Equal()

	plt.SaveD(dirout, fn)
}

func Test_functions03(tst *testing.T) {

	//verbose()
	chk.PrintTitle("functions03")

	eps := 1e-2
	f := func(x float64) float64 { return Sabs(x, eps) }
	ff := func(x float64) float64 { return SabsD1(x, eps) }

	np := 401
	//x  := utl.LinSpace(-5e5, 5e5, np)
	//x  := utl.LinSpace(-5e2, 5e2, np)
	x := utl.LinSpace(-5e1, 5e1, np)
	Y := make([]float64, np)
	y := make([]float64, np)
	g := make([]float64, np)
	h := make([]float64, np)
	tolg, tolh := 1e-6, 1e-5
	with_err := false
	for i := 0; i < np; i++ {
		Y[i] = math.Abs(x[i])
		y[i] = Sabs(x[i], eps)
		g[i] = SabsD1(x[i], eps)
		h[i] = SabsD2(x[i], eps)
		gnum := numderiv(f, x[i])
		hnum := numderiv(ff, x[i])
		errg := math.Abs(g[i] - gnum)
		errh := math.Abs(h[i] - hnum)
		clrg, clrh := "[1;32m", "[1;32m"
		if errg > tolg {
			clrg, with_err = "[1;31m", true
		}
		if errh > tolh {
			clrh, with_err = "[1;31m", true
		}
		io.Pf("errg = %s%23.15e   errh = %s%23.15e[0m\n", clrg, errg, clrh, errh)
	}

	if with_err {
		chk.Panic("errors found")
	}

	if false {
		//if true {
		plt.Subplot(3, 1, 1)
		plt.Plot(x, y, "'k--', label='abs'")
		plt.Plot(x, y, "'b-', label='sabs'")
		plt.Gll("x", "y", "")

		plt.Subplot(3, 1, 2)
		plt.Plot(x, g, "'b-', label='sabs'")
		plt.Gll("x", "dy/dx", "")

		plt.Subplot(3, 1, 3)
		plt.Plot(x, h, "'b-', label='sabs'")
		plt.Gll("x", "d2y/dx2", "")

		plt.Show()
	}
}

// PlotFltFlt plots flt-flt contour
// use iFlt==-1 || jFlt==-1 to plot all combinations
func (o *Optimiser) PlotFltFltContour(sols0 []*Solution, iFlt, jFlt, iOva int, pp *PlotParams) {
	best, _ := GetBestFeasible(o, iOva)
	plotAll := iFlt < 0 || jFlt < 0
	plotCommands := func(i, j int) {
		o.PlotContour(i, j, iOva, pp)
		if sols0 != nil {
			o.PlotAddFltFlt(i, j, sols0, &pp.FmtSols0)
		}
		o.PlotAddFltFlt(i, j, o.Solutions, &pp.FmtSols)
		if best != nil {
			plt.PlotOne(best.Flt[i], best.Flt[j], pp.FmtBest.GetArgs(""))
		}
		if pp.Extra != nil {
			pp.Extra()
		}
		if pp.AxEqual {
			plt.Equal()
		}
	}
	if plotAll {
		idx := 1
		ncol := o.Nflt - 1
		for row := 0; row < o.Nflt; row++ {
			idx += row
			for col := row + 1; col < o.Nflt; col++ {
				plt.Subplot(ncol, ncol, idx)
				plt.SplotGap(0.0, 0.0)
				plotCommands(col, row)
				if col > row+1 {
					plt.SetXnticks(0)
					plt.SetYnticks(0)
				} else {
					plt.Gll(io.Sf("$x_{%d}$", col), io.Sf("$x_{%d}$", row), "leg=0")
				}
				idx++
			}
		}
		idx = ncol*(ncol-1) + 1
		plt.Subplot(ncol, ncol, idx)
		plt.AxisOff()
		// TODO: fix formatting of open marker, add star to legend
		plt.DrawLegend([]plt.Fmt{pp.FmtSols0, pp.FmtSols, pp.FmtBest}, 8, "center", false, "")
	} else {
		plotCommands(iFlt, jFlt)
		if pp.Xlabel == "" {
			plt.Gll(io.Sf("$x_{%d}$", iFlt), io.Sf("$x_{%d}$", jFlt), pp.LegPrms)
		} else {
			plt.Gll(pp.Xlabel, pp.Ylabel, pp.LegPrms)
		}
	}
	plt.SaveD(pp.DirOut, pp.FnKey+pp.FnExt)
}

// Plot plot results
func Plot(dirout, fn string, res *Results, yfcn Cb_ycorr, xa, xb float64, argsAna, argsNum string, extra func()) {

	// data
	if res == nil {
		return
	}
	ndim := len(res.Y)
	if ndim < 1 {
		return
	}

	// closed-form solution
	var xc []float64
	var Yc [][]float64
	if yfcn != nil {
		np := 101
		dx := (xb - xa) / float64(np-1)
		xc = make([]float64, np)
		Yc = utl.DblsAlloc(np, ndim)
		for i := 0; i < np; i++ {
			xc[i] = xa + dx*float64(i)
			yfcn(Yc[i], xc[i])
		}
	}

	// plot
	if argsAna == "" {
		argsAna = "'y-', lw=6, label='analytical', clip_on=0"
	}
	if argsNum == "" {
		argsNum = "'b-', marker='.', lw=1, clip_on=0"
	}
	for j := 0; j < ndim; j++ {
		plt.Subplot(ndim+1, 1, j+1)
		if yfcn != nil {
			plt.Plot(xc, Yc[j], argsAna)
		}
		plt.Plot(res.X, res.Y[j], argsNum+","+io.Sf("label='%s'", res.Method))
		plt.Gll("$x$", "$y$", "")
	}
	plt.Subplot(ndim+1, 1, ndim+1)
	plt.Plot(res.X, res.Dx, io.Sf("'b-', marker='.', lw=1, clip_on=0, label='%s'", res.Method))
	plt.SetYlog()
	plt.Gll("$x$", "$\\log(\\delta x)$", "")

	// write file
	if extra != nil {
		extra()
	}
	plt.SaveD(dirout, fn)
}

// PlotT plots F, G and H for varying t and fixed coordinates x
func PlotT(o Func, dirout, fname string, t0, tf float64, xcte []float64, np int, args string, withG, withH, save, show bool, extra func()) {
	t := utl.LinSpace(t0, tf, np)
	y := make([]float64, np)
	for i := 0; i < np; i++ {
		y[i] = o.F(t[i], xcte)
	}
	var g, h []float64
	nrow := 1
	if withG {
		g = make([]float64, np)
		for i := 0; i < np; i++ {
			g[i] = o.G(t[i], xcte)
		}
		nrow += 1
	}
	if withH {
		h = make([]float64, np)
		for i := 0; i < np; i++ {
			h[i] = o.H(t[i], xcte)
		}
		nrow += 1
	}
	os.MkdirAll(dirout, 0777)
	if withG || withH {
		plt.Subplot(nrow, 1, 1)
	}
	plt.Plot(t, y, args)
	if extra != nil {
		extra()
	}
	plt.Gll("t", "y", "")
	pidx := 2
	if withG {
		plt.Subplot(nrow, 1, pidx)
		plt.Plot(t, g, args)
		plt.Gll("t", "dy/dt", "")
		pidx += 1
	}
	if withH {
		plt.Subplot(nrow, 1, pidx)
		plt.Plot(t, h, args)
		plt.Gll("t", "d2y/dt2", "")
	}
	if save {
		plt.Save(dirout + "/" + fname)
	}
	if show {
		plt.Show()
	}
}

func do_plot_nurbs_refined(b, c *Nurbs) {
	plt.SetForEps(1.5, 400)
	plt.Subplot(3, 1, 1)
	b.DrawCtrl2D(true)
	b.DrawElems2D(21, true, "", "")
	plt.Equal()

	plt.Subplot(3, 1, 2)
	c.DrawCtrl2D(true)
	c.DrawElems2D(21, true, "", "")
	plt.Equal()

	plt.Subplot(3, 1, 3)
	b.DrawElems2D(21, true, ", lw=3", "")
	c.DrawElems2D(21, true, ", color='red', marker='+', markevery=10", "color='magenta', size=8, va='bottom'")
	plt.Equal()
}

// subplot sets subplot
func (o *Plotter) Subplot() {
	if o.Split {
		plt.Clf()
		return
	}
	plt.Subplot(o.Nrow, o.Ncol, o.Pidx)
	o.Pidx += 1
}

func plot_normal(μ, σ float64) {

	var dist DistNormal
	dist.Init(&VarData{M: μ, S: σ})

	n := 101
	x := utl.LinSpace(-2, 2, n)
	y := make([]float64, n)
	Y := make([]float64, n)
	for i := 0; i < n; i++ {
		y[i] = dist.Pdf(x[i])
		Y[i] = dist.Cdf(x[i])
	}
	plt.Subplot(2, 1, 1)
	plt.Plot(x, y, io.Sf("clip_on=0,zorder=10,label=r'$\\mu=%g,\\;\\sigma=%g$'", μ, σ))
	plt.Gll("$x$", "$f(x)$", "leg_out=1, leg_ncol=2")
	plt.Subplot(2, 1, 2)
	plt.Plot(x, Y, io.Sf("clip_on=0,zorder=10,label=r'$\\mu=%g,\\;\\sigma=%g$'", μ, σ))
	plt.Gll("$x$", "$F(x)$", "leg_out=1, leg_ncol=2")
}

// Draw draws or save figure with plot
//  dirout -- directory to save figure
//  fname  -- file name; e.g. myplot.eps or myplot.png. Use "" to skip saving
//  show   -- shows figure
//  extra  -- is called just after Subplot command and before any plotting
//  Note: subplots will be split if using 'eps' files
func Draw(dirout, fname string, show bool, extra ExtraPlt) {
	var fnk string // filename key
	var ext string // extension
	var eps bool   // is eps figure
	if fname != "" {
		fnk = io.FnKey(fname)
		ext = io.FnExt(fname)
		eps = ext == ".eps"
	}
	nplots := len(Splots)
	nr, nc := utl.BestSquare(nplots)
	var k int
	for i := 0; i < nr; i++ {
		for j := 0; j < nc; j++ {
			if !eps {
				plt.Subplot(nr, nc, k+1)
			}
			if extra != nil {
				extra(i+1, j+1, nplots)
			}
			if Splots[k].Title != "" {
				plt.Title(Splots[k].Title, Splots[k].Topts)
			}
			data := Splots[k].Data
			for _, d := range data {
				if d.Style.L == "" {
					d.Style.L = d.Alias
				}
				x, y := d.X, d.Y
				if math.Abs(Splots[k].Xscale) > 0 {
					x = make([]float64, len(d.X))
					la.VecCopy(x, Splots[k].Xscale, d.X)
				}
				if math.Abs(Splots[k].Yscale) > 0 {
					y = make([]float64, len(d.Y))
					la.VecCopy(y, Splots[k].Yscale, d.Y)
				}
				plt.Plot(x, y, d.Style.GetArgs("clip_on=0"))
			}
			plt.Gll(Splots[k].Xlbl, Splots[k].Ylbl, Splots[k].GllArgs)
			if eps {
				savefig(dirout, fnk, ext, k)
				plt.Clf()
			}
			k += 1
		}
	}
	if !eps && fname != "" {
		savefig(dirout, fnk, ext, -1)
	}
	if show {
		plt.Show()
	}
}

func plot_frechet(l, c, a float64, xmin, xmax float64) {

	var dist DistFrechet
	dist.Init(&VarData{L: l, C: c, A: a})

	n := 101
	x := utl.LinSpace(xmin, xmax, n)
	y := make([]float64, n)
	Y := make([]float64, n)
	for i := 0; i < n; i++ {
		y[i] = dist.Pdf(x[i])
		Y[i] = dist.Cdf(x[i])
	}
	plt.Subplot(2, 1, 1)
	plt.Plot(x, y, io.Sf("clip_on=0,zorder=10,label=r'$(%g,%g,%g)$'", l, c, a))
	plt.Gll("$x$", "$f(x)$", "leg_out=1, leg_ncol=4, leg_hlen=1")
	plt.SetYnticks(11)
	plt.Subplot(2, 1, 2)
	plt.Plot(x, Y, io.Sf("clip_on=0,zorder=10,label=r'$(%g,%g,%g)$'", l, c, a))
	plt.Gll("$x$", "$F(x)$", "leg_out=1, leg_ncol=4, leg_hlen=1")
}

func plot_uniform(A, B float64, xmin, xmax float64) {

	var dist DistUniform
	dist.Init(&VarData{Min: A, Max: B})

	n := 101
	x := utl.LinSpace(xmin, xmax, n)
	y := make([]float64, n)
	Y := make([]float64, n)
	for i := 0; i < n; i++ {
		y[i] = dist.Pdf(x[i])
		Y[i] = dist.Cdf(x[i])
	}
	plt.Subplot(2, 1, 1)
	plt.Plot(x, y, io.Sf("clip_on=0,zorder=10,label=r'$(%g,%g)$'", A, B))
	plt.Gll("$x$", "$f(x)$", "leg_out=1, leg_ncol=4, leg_hlen=1")
	plt.SetYnticks(11)
	plt.Subplot(2, 1, 2)
	plt.Plot(x, Y, io.Sf("clip_on=0,zorder=10,label=r'$(%g,%g)$'", A, B))
	plt.Gll("$x$", "$F(x)$", "leg_out=1, leg_ncol=4, leg_hlen=1")
}

// PlotNurbsDerivs plots derivatives of basis functions la and lb
func PlotNurbsDerivs(dirout, fn string, b *Nurbs, la, lb int) {
	npts := 41
	plt.Reset()
	if io.FnExt(fn) == ".eps" {
		plt.SetForEps(1.5, 500)
	} else {
		plt.SetForPng(1.5, 600, 150)
	}

	plt.Subplot(4, 2, 1)
	t0 := time.Now()
	b.PlotDeriv(la, 0, "", npts, 0) // 0 => CalcBasisAndDerivs
	io.Pfcyan("time elapsed (calcbasis) = %v\n", time.Now().Sub(t0))
	plt.Equal()

	plt.Subplot(4, 2, 2)
	t0 = time.Now()
	b.PlotDeriv(la, 0, "", npts, 1) // 1 => NumericalDeriv
	io.Pfcyan("time elapsed (numerical) = %v\n", time.Now().Sub(t0))
	plt.Equal()

	plt.Subplot(4, 2, 3)
	b.PlotDeriv(la, 1, "", npts, 0) // 0 => CalcBasisAndDerivs
	plt.Equal()

	plt.Subplot(4, 2, 4)
	b.PlotDeriv(la, 1, "", npts, 1) // 0 => NumericalDeriv
	plt.Equal()

	plt.Subplot(4, 2, 5)
	b.PlotDeriv(lb, 0, "", npts, 0) // 0 => CalcBasisAndDerivs
	plt.Equal()

	plt.Subplot(4, 2, 6)
	b.PlotDeriv(lb, 0, "", npts, 1) // 0 => NumericalDeriv
	plt.Equal()

	plt.Subplot(4, 2, 7)
	b.PlotDeriv(lb, 1, "", npts, 0) // 0 => CalcBasisAndDerivs
	plt.Equal()

	plt.Subplot(4, 2, 8)
	b.PlotDeriv(lb, 1, "", npts, 1) // 0 => NumericalDeriv
	plt.Equal()

	plt.SaveD(dirout, fn)
}

func Test_bspline01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("bspline01")

	var s1 Bspline
	T1 := []float64{0, 0, 0, 1, 1, 1}
	s1.Init(T1, 2)
	s1.SetControl([][]float64{{0, 0}, {0.5, 1}, {1, 0}})

	var s2 Bspline
	T2 := []float64{0, 0, 0, 0.5, 1, 1, 1}
	s2.Init(T2, 2)
	s2.SetControl([][]float64{{0, 0}, {0.25, 0.5}, {0.75, 0.5}, {1, 0}})

	if chk.Verbose {
		npts := 201
		plt.SetForPng(1.5, 600, 150)
		plt.SplotGap(0.2, 0.4)

		str0 := ",lw=2"
		str1 := ",ls='none',marker='+',color='cyan',markevery=10"
		str2 := ",ls='none',marker='x',markevery=10"
		str3 := ",ls='none',marker='+',markevery=10"
		str4 := ",ls='none',marker='4',markevery=10"

		plt.Subplot(3, 2, 1)
		s1.Draw2d(str0, "", npts, 0) // 0 => CalcBasis
		s1.Draw2d(str1, "", npts, 1) // 1 => RecursiveBasis

		plt.Subplot(3, 2, 2)
		plt.SetAxis(0, 1, 0, 1)
		s2.Draw2d(str0, "", npts, 0) // 0 => CalcBasis
		s2.Draw2d(str1, "", npts, 1) // 1 => RecursiveBasis

		plt.Subplot(3, 2, 3)
		s1.PlotBasis("", npts, 0)   // 0 => CalcBasis
		s1.PlotBasis(str2, npts, 1) // 1 => CalcBasisAndDerivs
		s1.PlotBasis(str3, npts, 2) // 2 => RecursiveBasis

		plt.Subplot(3, 2, 4)
		s2.PlotBasis("", npts, 0)   // 0 => CalcBasis
		s2.PlotBasis(str2, npts, 1) // 1 => CalcBasisAndDerivs
		s2.PlotBasis(str3, npts, 2) // 2 => RecursiveBasis

		plt.Subplot(3, 2, 5)
		s1.PlotDerivs("", npts, 0)   // 0 => CalcBasisAndDerivs
		s1.PlotDerivs(str4, npts, 1) // 1 => NumericalDeriv

		plt.Subplot(3, 2, 6)
		s2.PlotDerivs("", npts, 0)   // 0 => CalcBasisAndDerivs
		s2.PlotDerivs(str4, npts, 1) // 1 => NumericalDeriv

		plt.SaveD("/tmp/gosl/gm", "bspline01.png")
	}
}

func Test_dist_lognormal_03(tst *testing.T) {

	//verbose()
	chk.PrintTitle("dist_lognormal_03. random numbers")

	μ := 1.0
	σ := 0.25

	nsamples := 1000
	X := make([]float64, nsamples)
	for i := 0; i < nsamples; i++ {
		X[i] = Lognormal(μ, σ)
	}

	nstations := 41
	xmin := 0.0
	xmax := 3.0
	dx := (xmax - xmin) / float64(nstations-1)

	var hist Histogram
	hist.Stations = utl.LinSpace(xmin, xmax, nstations)
	hist.Count(X, true)

	prob := make([]float64, nstations)
	for i := 0; i < nstations-1; i++ {
		prob[i] = float64(hist.Counts[i]) / (float64(nsamples) * dx)
	}

	io.Pf(TextHist(hist.GenLabels("%.3f"), hist.Counts, 60))
	io.Pforan("dx = %v\n", dx)

	area := 0.0
	for i := 0; i < nstations-1; i++ {
		area += dx * prob[i]
	}
	io.Pforan("area = %v\n", area)
	chk.Scalar(tst, "area", 1e-15, area, 1)

	if chk.Verbose {
		plt.SetForEps(1.5, 300)
		plot_lognormal(μ, σ)
		plt.Subplot(2, 1, 1)
		hist.PlotDensity(nil, "")
		plt.SaveD("/tmp/gosl", "rnd_dist_lognormal_03.eps")
	}
}

// PlotNurbsBasis plots basis functions la and lb
func PlotNurbsBasis(dirout, fn string, b *Nurbs, la, lb int) {
	npts := 41
	plt.Reset()
	if io.FnExt(fn) == ".eps" {
		plt.SetForEps(1.5, 500)
	} else {
		plt.SetForPng(1.5, 600, 150)
	}

	plt.Subplot(3, 2, 1)
	b.DrawCtrl2d(false, "", "")
	b.DrawElems2d(npts, false, "", "")
	t0 := time.Now()
	b.PlotBasis(la, "", 11, 0) // 0 => CalcBasis
	io.Pfcyan("time elapsed (calcbasis) = %v\n", time.Now().Sub(t0))
	plt.Equal()

	plt.Subplot(3, 2, 2)
	b.DrawCtrl2d(false, "", "")
	b.DrawElems2d(npts, false, "", "")
	b.PlotBasis(lb, "", 11, 0) // 0 => CalcBasis
	plt.Equal()

	plt.Subplot(3, 2, 3)
	b.DrawCtrl2d(false, "", "")
	b.DrawElems2d(npts, false, "", "")
	b.PlotBasis(la, "", 11, 1) // 1 => CalcBasisAndDerivs
	plt.Equal()

	plt.Subplot(3, 2, 4)
	b.DrawCtrl2d(false, "", "")
	b.DrawElems2d(npts, false, "", "")
	b.PlotBasis(lb, "", 11, 1) // 1 => CalcBasisAndDerivs
	plt.Equal()

	plt.Subplot(3, 2, 5)
	b.DrawCtrl2d(false, "", "")
	b.DrawElems2d(npts, false, "", "")
	t0 = time.Now()
	b.PlotBasis(la, "", 11, 2) // 2 => RecursiveBasis
	io.Pfcyan("time elapsed (recursive) = %v\n", time.Now().Sub(t0))
	plt.Equal()

	plt.Subplot(3, 2, 6)
	b.DrawCtrl2d(false, "", "")
	b.DrawElems2d(npts, false, "", "")
	b.PlotBasis(lb, "", 11, 2) // 2 => RecursiveBasis
	plt.Equal()

	plt.SaveD(dirout, fn)
}

func do_plot_nurbs_derivs(b *Nurbs, la, lb int) {
	np := 11
	plt.SetForEps(1.5, 500)

	plt.Subplot(4, 2, 1)
	t0 := time.Now()
	b.PlotDeriv(la, 0, "", np, 0) // 0 => CalcBasisAndDerivs
	io.Pfcyan("time elapsed (calcbasis) = %v\n", time.Now().Sub(t0))
	plt.Equal()

	plt.Subplot(4, 2, 2)
	t0 = time.Now()
	b.PlotDeriv(la, 0, "", np, 1) // 1 => NumericalDeriv
	io.Pfcyan("time elapsed (numerical) = %v\n", time.Now().Sub(t0))
	plt.Equal()

	plt.Subplot(4, 2, 3)
	b.PlotDeriv(la, 1, "", np, 0) // 0 => CalcBasisAndDerivs
	plt.Equal()

	plt.Subplot(4, 2, 4)
	b.PlotDeriv(la, 1, "", np, 1) // 0 => NumericalDeriv
	plt.Equal()

	plt.Subplot(4, 2, 5)
	b.PlotDeriv(lb, 0, "", np, 0) // 0 => CalcBasisAndDerivs
	plt.Equal()

	plt.Subplot(4, 2, 6)
	b.PlotDeriv(lb, 0, "", np, 1) // 0 => NumericalDeriv
	plt.Equal()

	plt.Subplot(4, 2, 7)
	b.PlotDeriv(lb, 1, "", np, 0) // 0 => CalcBasisAndDerivs
	plt.Equal()

	plt.Subplot(4, 2, 8)
	b.PlotDeriv(lb, 1, "", np, 1) // 0 => NumericalDeriv
	plt.Equal()
}

func Test_fun01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("fun01. Decreasing Reference Model")

	ya := 1.0
	yb := -0.5
	λ1 := 1.0

	o, err := New("ref-dec-gen", []*Prm{
		&Prm{N: "bet", V: 5.0},
		&Prm{N: "a", V: -λ1},
		&Prm{N: "b", V: -1.0},
		&Prm{N: "c", V: ya},
		&Prm{N: "A", V: 0.0},
		&Prm{N: "B", V: λ1},
		&Prm{N: "xini", V: 0.0},
		&Prm{N: "yini", V: yb},
	})
	if err != nil {
		tst.Errorf("test failed: %v\n", err)
		return
	}

	tmax := 3.0
	xcte := []float64{0, 0, 0}
	if chk.Verbose {
		plt.SetForPng(1.2, 400, 150)
		PlotT(o, "", "", 0.0, tmax, xcte, 201, "", "", "", "", "label='f'", "label='g'", "label='h'")
		plt.Subplot(3, 1, 1)
		plt.Plot([]float64{0, tmax}, []float64{ya, ya - λ1*tmax}, "'k-'")
		plt.Equal()
		plt.SaveD("/tmp/gosl/fun", "ref-dec-gen.png")
	}
	//
	sktol := 1e-10
	dtol := 1e-10
	dtol2 := 1e-10
	ver := chk.Verbose
	CheckDerivT(tst, o, 0.0, tmax, xcte, 11, nil, sktol, dtol, dtol2, ver)
}

func do_plot_nurbs_basis(b *Nurbs, la, lb int) {
	npts := 21
	plt.SetForEps(1.2, 500)

	plt.Subplot(3, 2, 1)
	b.DrawCtrl2D(false)
	b.DrawElems2D(npts, false, "", "")
	t0 := time.Now()
	b.PlotBasis(la, "", 11, 0) // 0 => CalcBasis
	io.Pfcyan("time elapsed (calcbasis) = %v\n", time.Now().Sub(t0))
	plt.Equal()

	plt.Subplot(3, 2, 2)
	b.DrawCtrl2D(false)
	b.DrawElems2D(npts, false, "", "")
	b.PlotBasis(lb, "", 11, 0) // 0 => CalcBasis
	plt.Equal()

	plt.Subplot(3, 2, 3)
	b.DrawCtrl2D(false)
	b.DrawElems2D(npts, false, "", "")
	b.PlotBasis(la, "", 11, 1) // 1 => CalcBasisAndDerivs
	plt.Equal()

	plt.Subplot(3, 2, 4)
	b.DrawCtrl2D(false)
	b.DrawElems2D(npts, false, "", "")
	b.PlotBasis(lb, "", 11, 1) // 1 => CalcBasisAndDerivs
	plt.Equal()

	plt.Subplot(3, 2, 5)
	b.DrawCtrl2D(false)
	b.DrawElems2D(npts, false, "", "")
	t0 = time.Now()
	b.PlotBasis(la, "", 11, 2) // 2 => RecursiveBasis
	io.Pfcyan("time elapsed (recursive) = %v\n", time.Now().Sub(t0))
	plt.Equal()

	plt.Subplot(3, 2, 6)
	b.DrawCtrl2D(false)
	b.DrawElems2D(npts, false, "", "")
	b.PlotBasis(lb, "", 11, 2) // 2 => RecursiveBasis
	plt.Equal()
}

func Test_bspline03(tst *testing.T) {

	//verbose()
	chk.PrintTitle("bspline03")

	//             0 1 2 3 4 5 6 7 8 9 10
	T := []float64{0, 0, 0, 1, 2, 3, 4, 4, 5, 5, 5}
	var s Bspline
	s.Init(T, 2)
	s.SetControl([][]float64{{0, 0}, {0.5, 1}, {1, 0}, {1.5, 0}, {2, 1}, {2.5, 1}, {3, 0.5}, {3.5, 0}})

	// analytical derivatives
	s.CalcBasisAndDerivs(3.99)
	io.Pfpink("ana: dNdt(t=3.99, i=5) = %v\n", s.GetDeriv(5))
	io.Pfpink("ana: dNdt(t=3.99, i=6) = %v\n", s.GetDeriv(6))
	io.Pfpink("ana: dNdt(t=3.99, i=7) = %v\n", s.GetDeriv(7))
	s.CalcBasisAndDerivs(4.0)
	io.Pforan("ana: dNdt(t=4.00, i=5) = %v\n", s.GetDeriv(5))
	io.Pforan("ana: dNdt(t=4.00, i=6) = %v\n", s.GetDeriv(6))
	io.Pforan("ana: dNdt(t=4.00, i=7) = %v\n", s.GetDeriv(7))

	// numerical derivatives
	io.Pfcyan("num: dNdt(t=3.99, i=5) = %v\n", s.NumericalDeriv(3.99, 5))
	io.Pfcyan("num: dNdt(t=3.99, i=6) = %v\n", s.NumericalDeriv(3.99, 6))
	io.Pfcyan("num: dNdt(t=3.99, i=7) = %v\n", s.NumericalDeriv(3.99, 7))
	io.Pfblue2("num: dNdt(t=4.00, i=5) = %v\n", s.NumericalDeriv(4.00, 5))
	io.Pfblue2("num: dNdt(t=4.00, i=6) = %v\n", s.NumericalDeriv(4.00, 6))
	io.Pfblue2("num: dNdt(t=4.00, i=7) = %v\n", s.NumericalDeriv(4.00, 7))

	ver := false
	tol := 1e-5
	tt := utl.LinSpace(0, 5, 11)
	numd := make([]float64, s.NumBasis())
	anad := make([]float64, s.NumBasis())
	for _, t := range tt {
		for i := 0; i < s.NumBasis(); i++ {
			s.CalcBasisAndDerivs(t)
			anad[i] = s.GetDeriv(i)
			numd[i] = s.NumericalDeriv(t, i)
			// numerical fails @ 4 [4,5,6]
			if t == 4 {
				numd[4] = anad[4]
				numd[5] = anad[5]
				numd[6] = anad[6]
			}
			chk.PrintAnaNum(io.Sf("i=%d t=%v", i, t), tol, anad[i], numd[i], ver)
		}
		chk.Vector(tst, io.Sf("derivs @ %v", t), tol, numd, anad)
	}

	if chk.Verbose {

		npts := 201
		plt.SetForPng(1.5, 600, 150)
		plt.SplotGap(0, 0.3)

		str0 := ",lw=2"
		str1 := ",ls='none',marker='+',color='cyan',markevery=10"
		str2 := ",ls='none',marker='x',markevery=10"
		str3 := ",ls='none',marker='+',markevery=10"
		str4 := ",ls='none',marker='4',markevery=10"

		plt.Subplot(3, 1, 1)
		s.Draw2d(str0, "", npts, 0) // 0 => CalcBasis
		s.Draw2d(str1, "", npts, 1) // 1 => RecursiveBasis

		plt.Subplot(3, 1, 2)
		s.PlotBasis("", npts, 0)   // 0 => CalcBasis
		s.PlotBasis(str2, npts, 1) // 1 => CalcBasisAndDerivs
		s.PlotBasis(str3, npts, 2) // 2 => RecursiveBasis

		plt.Subplot(3, 1, 3)
		s.PlotDerivs("", npts, 0)   // 0 => CalcBasisAndDerivs
		s.PlotDerivs(str4, npts, 1) // 1 => NumericalDeriv

		plt.SaveD("/tmp/gosl/gm", "bspline03.png")
	}
}