%%	options

copyright owner	=	Dirk Krause
copyright year	=	2012-2014
license		=	bsd

%%	header

#include "dk3all.h"

#ifdef __cplusplus
extern "C" {
#endif

/**	Reset spline segment data.
	All components - including cb - are reset.
	@param	pseg	Spline segment data to reset.
*/
void
dk3xsp_reset(dk3_xspline_segment_t *pseg);

/**	Set compatility bug flag.
	@param	pseg	Spline segment data to modify.
	@param	v	New flag value (1 or 0).
*/
void
dk3xsp_set_cb(dk3_xspline_segment_t *pseg, int v);

/**	Set control points.
	@param	pseg	Spline segment data to modify.
	@param	a	Control point A (may be NULL).
	@param	b	Control point B.
	@param	c	Control point C.
	@param	d	Control point D (may be NULL).
*/
void
dk3xsp_set(
  dk3_xspline_segment_t		*pseg,
  dk3_fig_spline_point_t const	*a,
  dk3_fig_spline_point_t const	*b,
  dk3_fig_spline_point_t const	*c,
  dk3_fig_spline_point_t const	*d
);

/**	Run calculation for a given t.
	@param	pseg	Spline segment.
	@param	t	The t parameter.
	@param	fder	Flag: Calculate derivatives (1=yes, 0=no).
	@return	1 on success, 0 on error.
*/
int
dk3xsp_calculate(
  dk3_xspline_segment_t		*pseg,
  double			 t,
  int				 fder
);

/**	Retrieve calculation result.
	@param	pseg	Spline segment.
	@return	Calculation result: x value.
*/
double
dk3xsp_get_x(dk3_xspline_segment_t const *pseg);

/**	Retrieve calculation result.
	@param	pseg	Spline segment.
	@return	Calculation result: y value.
*/
double
dk3xsp_get_y(dk3_xspline_segment_t const *pseg);

/**	Retrieve calculation result.
	@param	pseg	Spline segment.
	@return	Calculation result: dx/dt value.
*/
double
dk3xsp_get_dxdt(dk3_xspline_segment_t const *pseg);

/**	Retrieve calculation result.
	@param	pseg	Spline segment.
	@return	Calculation result: dy/dt value.
*/
double
dk3xsp_get_dydt(dk3_xspline_segment_t const *pseg);

/**	Calculate X-spline segment interval length.
	@param	pseg	Segment to calculate the length for.
	@param	prec	Precision.
	@param	ec	Pointer to error code variable, may be NULL.
	@return	Length (positive value) on success, negative value on error.
*/
double
dk3xsp_segment_length(
  dk3_xspline_segment_t 	*pseg,
  double			 prec,
  int				*ec
);

/**	Calculate X-spline partial segment length.
	@param	pseg	Segment to calculate the length for.
	@param	tstart	Parameter t for start.
	@param	dir	Direction (0=forward t to 1, 1=backward 0 to t).
	@param	prec	Precision.
	@param	ec	Pointer to error code variable, may be NULL.
	@return	Length (positive value) on success, negative value on error.
*/
double
dk3xsp_segment_partial_length(
  dk3_xspline_segment_t		*pseg,
  double			 tstart,
  int				 dir,
  double			 prec,
  int				*ec
);

/**	Check whether an s value is very close to 0.
	@param	s	Value to check.
	@return	1 if s==0, 0 otherwise.
*/
int
dk3xsp_is_null(double s);

/**	Calculate length of open X-spline.
	@param	points	Array of spline points.
	@param	np	Number of points in array.
	@param	prec	Numeric precision.
	@param	fslc	Full segment length cache (np - 1) elements.
	@param	cb	Flag: Compatibility bug.
	@param	ec	Pointer to error code variable, may be NULL.
	@return	Length (positive value) on success, negative value on error.
*/
double
dk3xsp_open_spline_length(
  dk3_fig_spline_point_t const	*points,
  size_t			 np,
  double			 prec,
  double			*fslc,
  int				 cb,
  int				*ec
);

/**	Calculate length of open X-spline.
	@param	points	Array of spline points.
	@param	np	Number of points in array.
	@param	prec	Numeric precision.
	@param	tstart	Start/end t parameter.
	@param	fslc	Full segment length cache, (np-1) elements.
	@param	dir	Direction (0=forward t to 1, 1=backward 0 to t).
	@param	cb	Flag: Compatibility bug.
	@param	ec	Pointer to error code variable, may be NULL.
	@return	Length (positive value) on success, negative value on error.
*/
double
dk3xsp_open_spline_partial_length(
  dk3_fig_spline_point_t const	*points,
  size_t			 np,
  double			 prec,
  double			 tstart,
  double			*fslc,
  int				 dir,
  int				 cb,
  int				*ec
);

/**	Calculate t value for given distance to end point.
	@param	points	X-spline points.
	@param	np	Number of x-spline points.
	@param	dist	Distance to endpoint.
	@param	prec	Precision.
	@param	fslc	Cache for full segment lengths.
	@param	dir	End-point type: 0=end of spline, 1=start of spline.
	@param	cb	Flag: Compatibility with spline bug.
	@param	ec	Pointer to error code variable.
	@return	Non-negative value on success, negative value on error.
*/
double
dk3xsp_get_t_for_distance(
  dk3_fig_spline_point_t const	*points,
  size_t			 np,
  double			 dist,
  double			 prec,
  double			*fslc,
  int				 dir,
  int				 cb,
  int				*ec
);

/**	Run calculations for a given t value in the range 0<=t<=(np-1).
	@param	pseg	Segment to use for calculation.
	@param	points	Spline points.
	@param	np	Number of points.
	@param	t	Parameter t.
	@param	cb	Flag: Compatibility for spline bug.
	@param	der	Flag: Calculate derivatives dx/dt and dy/dt too.
	@param	ec	Pointer to error code variable.
	@return	1 on success, 0 on error.
*/
int
dk3xsp_calculate_position(
  dk3_xspline_segment_t		*pseg,
  dk3_fig_spline_point_t const	*points,
  size_t			 np,
  double			 t,
  int				 cb,
  int				 der,
  int				*ec
);

/**	Find movement direction for a specified t value.
	@param	points	X-spline points.
	@param	np	Number of points.
	@param	t	The t parameter 0<=t<=(np-1).
	@param	cb	Flag: Compatibility for spline bug.
	@param	ec	Pointer to error code variable.
	@return	Non-negative value 0<=res<=2*pi on success,
	negative value on error.
*/
double
dk3xsp_direction(
  dk3_fig_spline_point_t const	*points,
  size_t			 np,
  double			 t,
  int				 cb,
  int				*ec
);

/**	Find position and movement direction for a specified t value.
	@param	xpos	Pointer to variable for x position.
	@param	ypos	Pointer to variable for y position.
	@param	points	X-spline points.
	@param	np	Number of points.
	@param	t	The t parameter 0<=t<=(np-1).
	@param	cb	Flag: Compatibility for spline bug.
	@param	ec	Pointer to error code variable.
	@return	Non-negative value 0<=res<=2*pi on success,
	negative value on error.
*/
double
dk3xsp_direction_and_position(
  double			*xpos,
  double			*ypos,
  dk3_fig_spline_point_t const	*points,
  size_t			 np,
  double			 t,
  int				 cb,
  int				*ec
);

#ifdef __cplusplus
}
#endif

#ifndef DK3XSP_MAX_ITERATION_SEGMENTS
/**	Maximum number of sub-segments in segment length iteration.
*/
#define	DK3XSP_MAX_ITERATION_SEGMENTS 1073741824L
#endif

%%	module

#include "dk3all.h"
#include "dk3fig.h"
#include "dk3xsp.h"

$!trace-include



$*
Basic X-spline functions.
$*



$!trace-off



/**	Initialize spline point (set all components to 0.0).
	@param	ppt	Point to initialize.
*/
static
void
dk3xsp_init_spline_point(dk3_fig_spline_point_t *ppt)
{
  ppt->x = 0.0; ppt->y = 0.0; ppt->s = 0.0;
}



/**	Copy spline point or initialize a destination point.
	@param	pd	Destination pointer.
	@param	ps	Source pointer, may be NULL.
*/
static
void
dk3xsp_copy_spline_point(
  dk3_fig_spline_point_t *pd, dk3_fig_spline_point_t const *ps
)
{
  if(pd) {
    if(ps) {
      pd->x = ps->x; pd->y = ps->y; pd->s = ps->s;
    } else {
      dk3xsp_init_spline_point(pd);
    }
  }
}



void
dk3xsp_reset(dk3_xspline_segment_t *pseg)
{
  if(pseg) {
    dk3mem_res((void *)pseg, sizeof(dk3_xspline_segment_t));
    dk3xsp_init_spline_point(&(pseg->a));
    dk3xsp_init_spline_point(&(pseg->b));
    dk3xsp_init_spline_point(&(pseg->c));
    dk3xsp_init_spline_point(&(pseg->d));
    pseg->pa = pseg->pb = pseg->pc = pseg->pd = 2.0;
    pseg->qa = pseg->qb = pseg->qc = pseg->qd = 0.0;
    pseg->x = pseg->y = pseg->dxdt = pseg->dydt = 0.0;
    pseg->duadt = pseg->dubdt = pseg->ducdt = pseg->duddt = 0.0;
    pseg->cb = 0;
  }
}



void
dk3xsp_set_cb(dk3_xspline_segment_t *pseg, int v)
{
  if(pseg) {
    pseg->cb = ((v) ? 1 : 0);
  }
}



/**	The f(u) function.
	@param	u	The u argument.
	@param	p	The p parameter.
	@return	Function result.
*/
static
double
dk3xsp_f(double u, double p)
{
  double back;
  back = u * u * u *(10.0 - p + u * (2.0 * p - 15.0 + u * (6.0 - p)));
  $? "= f(%lg, %lg) = %lg", u, p, back
  return back;
}



/**	The f(u) function derived.
	@param	u	The u argument.
	@param	p	The p parameter.
	@return	Function result.
*/
static
double
dk3xsp_dfdu(double u, double p)
{
  double back;
  back = u * u * (30.0 - 3.0*p + u * (8.0*p - 60.0 + u * (30.0 - 5.0*p)));
  $? "= df/dt (%lg, %lg) = %lg", u, p, back
  return back;
}



/**	The g(u) function.
	@param	u	Argument u.
	@param	p	Parameter p.
	@param	q	Parameter q.
	@return	Function result.
*/
static
double
dk3xsp_g(double u, double p, double q)
{
  double	back;
  back = u*(q+u*(2.0*q+u*(10.0-12.0*q-p+u*(2.0*p+14.0*q-15+u*(6.0-5.0*q-p)))));
  $? "= g(%lg, %lg, %lg) = %lg", u, p, q, back
  return back;
}



/**	The g(u) function derived.
	@param	u	Argument u.
	@param	p	Parameter p.
	@param	q	Parameter q.
	@return	Function result.
*/
static
double
dk3xsp_dgdu(double u, double p, double q)
{
  double back;
  back =
  q+u*(4.0*q+u*(30.0-3.0*p-36.0*q+u*(56.0*q+8.0*p-60.0+u*(30.0-5.0*p-25.0*q))));
  $? "= dg/dt (%lg, %lg, %lg) = %lg", u, p, q, back
  return back;
}



/**	The h(u) function.
	@param	u	Argument u.
	@param	q	Parameter q.
	@return	Function result.
*/
static
double
dk3xsp_h(double u, double q)
{
  double back;
  back = u*(q+u*(2.0*q-u*u*(q*(u+2.0))));
  $? "= h(%lg, %lg) = %lg", u, q, back
  return back;
}



/**	The h(u) function derived.
	@param	u	Argument u.
	@param	q	Parameter q.
	@return	Function result.
*/
static
double
dk3xsp_dhdu(double u, double q)
{
  double back;
  back = q+u*(4.0*q-u*u*(q*(8.0+5.0*u)));
  $? "= dh/dt (%lg, %lg) = %lg", u, q, back
  return back;
}



double
dk3xsp_get_x(dk3_xspline_segment_t const *pseg)
{
  double back = 0.0;
  if(pseg) {
    back = pseg->x;
  }
  return back;
}



double
dk3xsp_get_y(dk3_xspline_segment_t const *pseg)
{
  double back = 0.0;
  if(pseg) {
    back = pseg->y;
  }
  return back;
}



double
dk3xsp_get_dxdt(dk3_xspline_segment_t const *pseg)
{
  double back = 0.0;
  if(pseg) {
    back = pseg->dxdt;
  }
  return back;
}



double
dk3xsp_get_dydt(dk3_xspline_segment_t const *pseg)
{
  double back = 0.0;
  if(pseg) {
    back = pseg->dydt;
  }
  return back;
}



void
dk3xsp_set(
  dk3_xspline_segment_t		*pseg,
  dk3_fig_spline_point_t const	*a,
  dk3_fig_spline_point_t const	*b,
  dk3_fig_spline_point_t const	*c,
  dk3_fig_spline_point_t const	*d
)
{
  $? "+ dk3xsp_set"
  if(pseg) {
    pseg->ha = ((a) ? 1 : 0);
    pseg->hd = ((d) ? 1 : 0);
    dk3xsp_copy_spline_point(&(pseg->a), a);
    dk3xsp_copy_spline_point(&(pseg->b), b);
    dk3xsp_copy_spline_point(&(pseg->c), c);
    dk3xsp_copy_spline_point(&(pseg->d), d);
    pseg->pa = pseg->pb = pseg->pc = pseg->pd = 0.0;
    pseg->qa = pseg->qb = pseg->qc = pseg->qd = 0.0;
    pseg->duadt = pseg->dubdt = pseg->ducdt = pseg->duddt = 0.0;
    if((pseg->b).s < 0.0) {
      if(pseg->cb) {
        pseg->qa = pseg->qc = 0.0 - ((pseg->b).s);
      } else {
        pseg->qa = pseg->qc = -0.5 * ((pseg->b).s);
      }
      pseg->duadt = -1.0; pseg->ducdt = 1.0;
    } else {
      pseg->pa = pseg->pc = 2.0 * (1.0 + (pseg->b).s) * (1.0 + (pseg->b).s);
      pseg->ducdt = 1.0 / (1.0 + (pseg->b).s);
      pseg->duadt = -1.0 * pseg->ducdt;
    }
    if((pseg->c).s < 0.0) {
      if(pseg->cb) {
        pseg->qb = pseg->qd = 0.0 - ((pseg->c).s);
      } else {
        pseg->qb = pseg->qd = -0.5 * ((pseg->c).s);
      }
      pseg->dubdt = -1.0; pseg->duddt = 1.0;
    } else {
      pseg->pb = pseg->pd = 2.0 * (1.0 + (pseg->c).s) * (1.0 + (pseg->c).s);
      pseg->duddt = 1.0 / (1.0 + (pseg->c).s);
      pseg->dubdt = -1.0 * pseg->duddt;
    }
    $? ". pa = %lg\tqa = %lg\tdua/dt = %lg", pseg->pa, pseg->qa, pseg->duadt
    $? ". pb = %lg\tqb = %lg\tdub/dt = %lg", pseg->pb, pseg->qb, pseg->dubdt
    $? ". pc = %lg\tqc = %lg\tduc/dt = %lg", pseg->pc, pseg->qc, pseg->ducdt
    $? ". pd = %lg\tqd = %lg\tdud/dt = %lg", pseg->pd, pseg->qd, pseg->duddt
  } $? "- dk3xsp_set"
}



int
dk3xsp_calculate(
  dk3_xspline_segment_t		*pseg,
  double			 t,
  int				 fder
)
{
  double	 ua;	/* Argument u for blending function from point A. */
  double	 ub;	/* = B. */
  double	 uc;	/* = C. */
  double	 ud;	/* = D. */
  double	 fa;	/* Result from blending function from point A. */
  double	 fb;	/* = B. */
  double	 fc;	/* = C. */
  double	 fd;	/* = D. */
  double	 zx;	/* Counter in fraction to calculate x. */
  double	 zy;	/* Counter in fraction to calculate y. */
  double	 n;	/* Denominator in both fractions. */
  double	 dfadt;	/* Derivative result for point A blending function. */
  double	 dfbdt;	/* = B. */
  double	 dfcdt;	/* = C. */
  double	 dfddt;	/* = D. */
  double	 nder;	/* Derivative of fraction denominator. */
  double	 zxder;	/* Derivative of x fraction counter. */
  double	 zyder;	/* Derivative of y fraction counter. */
  double	 nsq;	/* Denominator square value. */
  int		 ec = 0;	/* Error code for mathematical operations. */
  int		 back = 0;
  $? "+ dk3xsp_calculate %lg (fder=%d)", t, fder
  if(pseg) {
    back = 1;
    /*
    	Initialize variables.
    */
    ua = ub = uc = ud = fa = fb = fc = fd = zx = zy = n = 0.0;
    /*
    	Set ua, ub, uc, ud.
    */
    if((pseg->b).s < 0.0) {
      ua = -1.0 * t;
      uc = t;
    } else {
      ua = ((pseg->b).s - t) / (1.0 + (pseg->b).s);
      uc = ((pseg->b).s + t) / (1.0 + (pseg->b).s);
    }
    if((pseg->c).s < 0.0) {
      ub = 1.0 - t;
      ud = t - 1.0;
    } else {
      ub = (1.0 + (pseg->c).s - t) / (1.0 + (pseg->c).s);
      ud = (t + (pseg->c).s - 1.0) / (1.0 + (pseg->c).s);
    }
    $? ". ua = %lg", ua
    $? ". ub = %lg", ub
    $? ". uc = %lg", uc
    $? ". ud = %lg", ud
    /*
    	Calculate fa, fb, fc, fd.
    */
    if((pseg->b).s < 0.0) {
      fa = dk3xsp_h(ua, pseg->qa);
      fc = dk3xsp_g(uc, pseg->pc, pseg->qc);
    } else {
      if(t < (pseg->b).s) {
        fa = dk3xsp_f(ua, pseg->pa);
      }
      fc = dk3xsp_f(uc, pseg->pc);
    }
    if((pseg->c).s < 0.0) {
      fb = dk3xsp_g(ub, pseg->pb, pseg->qb);
      fd = dk3xsp_h(ud, pseg->qd);
    } else {
      fb = dk3xsp_f(ub, pseg->pb);
      if(t > (1.0 - (pseg->c).s)) {
        fd =  dk3xsp_f(ud, pseg->pd);
      }
    }
    $? ". fa = %lg", fa
    $? ". fb = %lg", fb
    $? ". fc = %lg", fc
    $? ". fd = %lg", fd
    /*
    	Calculate counters and denominator.
    */
    zx = dk3ma_d_add_ok(
      dk3ma_d_add_ok(
        ((pseg->ha) ? (fa * (pseg->a).x) : 0.0),
	(fb * (pseg->b).x),
	&ec
      ),
      dk3ma_d_add_ok(
        (fc * (pseg->c).x),
	((pseg->hd) ? (fd * (pseg->d).x) : 0.0),
	&ec
      ),
      &ec
    );
    zy = dk3ma_d_add_ok(
      dk3ma_d_add_ok(
        ((pseg->ha) ? (fa * (pseg->a).y) : 0.0),
	(fb * (pseg->b).y),
	&ec
      ),
      dk3ma_d_add_ok(
        (fc * (pseg->c).y),
	((pseg->hd) ? (fd * (pseg->d).y) : 0.0),
	&ec
      ),
      &ec
    );
    n  = dk3ma_d_add_ok(
      dk3ma_d_add_ok(((pseg->ha) ? fa : 0.0), fb, &ec),
      dk3ma_d_add_ok(fc, ((pseg->hd) ? fd : 0.0), &ec),
      &ec
    );
    $? ". zx = %lg", zx
    $? ". zy = %lg", zy
    $? ". n  = %lg", n
    /*
    	Calculate x and y value.
    */
    pseg->x = dk3ma_d_div_ok(zx, n, &ec);
    pseg->y = dk3ma_d_div_ok(zy, n, &ec);
    $? ". x = %lg", pseg->x
    $? ". y = %lg", pseg->y
    /*
    	Calculate derivatives if required.
    */
    if(fder) {
      /*
      	Initialize variables.
      */
      dfadt = dfbdt = dfcdt = dfddt = nder = zxder = zyder = 0.0;
      /*
      	Calculate dfadt, dfbdt, dfcdt, dfddt.
      */
      if((pseg->b).s < 0.0) {
        dfadt = dk3ma_d_mul_ok(dk3xsp_dhdu(ua, pseg->qa), pseg->duadt, &ec);
	dfcdt = dk3ma_d_mul_ok(
	  dk3xsp_dgdu(uc, pseg->pc, pseg->qc),
	  pseg->ducdt,
	  &ec
	);
      } else {
        if(t < (pseg->b).s) {
          dfadt = dk3ma_d_mul_ok(dk3xsp_dfdu(ua, pseg->pa), pseg->duadt, &ec);
	}
	dfcdt = dk3ma_d_mul_ok(dk3xsp_dfdu(uc, pseg->pc), pseg->ducdt, &ec);
      }
      if((pseg->c).s < 0.0) {
        dfbdt = dk3ma_d_mul_ok(
	  dk3xsp_dgdu(ub, pseg->pb, pseg->qb),
	  pseg->dubdt,
	  &ec
	);
	dfddt = dk3ma_d_mul_ok(dk3xsp_dhdu(ud, pseg->qd), pseg->duddt, &ec);
      } else {
        dfbdt = dk3ma_d_mul_ok(dk3xsp_dfdu(ub, pseg->pb), pseg->dubdt, &ec);
	if(t > (1.0 - (pseg->c).s)) {
	  dfddt = dk3ma_d_mul_ok(dk3xsp_dfdu(ud, pseg->pd), pseg->duddt, &ec);
	}
      }
      $? ". dfadt = %lg", dfadt
      $? ". dfbdt = %lg", dfbdt
      $? ". dfcdt = %lg", dfcdt
      $? ". dfddt = %lg", dfddt
      /*
      	Calculate counters and denominator.
      */
      zxder = dk3ma_d_add_ok(
        dk3ma_d_add_ok(
	  ((pseg->ha) ? dk3ma_d_mul_ok(dfadt, (pseg->a).x, &ec) : 0.0),
	  dk3ma_d_mul_ok(dfbdt, (pseg->b).x, &ec),
	  &ec
	),
	dk3ma_d_add_ok(
	  dk3ma_d_mul_ok(dfcdt, (pseg->c).x, &ec),
	  ((pseg->hd) ? dk3ma_d_mul_ok(dfddt, (pseg->d).x, &ec) : 0.0),
	  &ec
	),
	&ec
      );
      zyder = dk3ma_d_add_ok(
        dk3ma_d_add_ok(
	  ((pseg->ha) ? dk3ma_d_mul_ok(dfadt, (pseg->a).y, &ec) : 0.0),
	  dk3ma_d_mul_ok(dfbdt, (pseg->b).y, &ec),
	  &ec
	),
	dk3ma_d_add_ok(
	  dk3ma_d_mul_ok(dfcdt, (pseg->c).y, &ec),
	  ((pseg->hd) ? dk3ma_d_mul_ok(dfddt, (pseg->d).y, &ec) : 0.0),
	  &ec
	),
	&ec
      );
      nder  = dk3ma_d_add_ok(
        dk3ma_d_add_ok(((pseg->ha) ? dfadt : 0.0), dfbdt, &ec),
	dk3ma_d_add_ok(dfcdt, ((pseg->hd) ? dfddt : 0.0), &ec),
	&ec
      );
      nsq   = n * n;
      $? ". zxder = %lg", zxder
      $? ". zyder = %lg", zyder
      $? ". nder  = %lg", nder
      $? ". n^2   = %lg", nsq
      pseg->dxdt = dk3ma_d_div_ok(
        dk3ma_d_sub_ok(
	  dk3ma_d_mul_ok(n, zxder, &ec),
	  dk3ma_d_mul_ok(zx, nder, &ec),
	  &ec
	),
	nsq,
	&ec
      );
      pseg->dydt = dk3ma_d_div_ok(
        dk3ma_d_sub_ok(
	  dk3ma_d_mul_ok(n, zyder, &ec),
	  dk3ma_d_mul_ok(zy, nder, &ec),
	  &ec
	),
	nsq,
	&ec
      );
      $? ". dx/dt = %lg", pseg->dxdt
      $? ". dy/dt = %lg", pseg->dydt
    }
    if(ec) {
      back = 0;
    }
  } $? "- dk3xsp_calculate %d", back
  return back;
}



$*
Length calculations for single segments.
$*



int
dk3xsp_is_null(double s)
{
  int		 back = 0;
  if(fabs(s) < 1.0e-6) {
    back = 1;
  }
  return back;
}



$!trace-on



/**	Convert long to double.
	@param	l	Long value to convert.
	@return	Conversion result.
*/
static
double
dk3xsp_l_to_d(long l)
{
  return ((double)l);
}

double
dk3xsp_segment_length(
  dk3_xspline_segment_t 	*pseg,
  double			 prec,
  int				*ec
)
{
  double		 back = -1.0;
  double		 lastpassres;		/* Result of previous pass. */
  double		 passres;		/* Result of current pass. */
  double		 deltax;		/* Delta x between points. */
  double		 deltay;		/* Delta y between points. */
  double		 lastx;			/* Previous x. */
  double		 lasty;			/* Previous y. */
  double		 x = 0.0;		/* Current x. */
  double		 y = 0.0;		/* Current y. */
  double		 t;			/* Current t. */
  double		 dsegs;			/* Double value of segs. */
  long			 segs = 64L;		/* Number of segments. */
  long			 i;			/* Current segment. */
  int			 mec = 0;		/* Mathematical error code. */
  int			 mc = 1;		/* Flag: Must continue. */
  int			 cc = 1;		/* Flag: Can continue. */
  $? "+ dk3xsp_segment_length prec=%lg", prec
  if((pseg) && (prec > 0.0)) {
    if(dk3xsp_is_null((pseg->b).s) && dk3xsp_is_null((pseg->c).s)) {
      /*
      	Simple line.
      */
      deltax = dk3ma_d_sub_ok((pseg->c).x, (pseg->b).x, &mec);
      deltay = dk3ma_d_sub_ok((pseg->c).y, (pseg->b).y, &mec);
      back = sqrt(
        dk3ma_d_add_ok(
	  dk3ma_d_mul_ok(deltax, deltax, &mec),
	  dk3ma_d_mul_ok(deltay, deltay, &mec),
	  &mec
	)
      );	$? ". straight-forward line"
    } else {	$? ". real spline"
      /*
      	Spline calculation really necessary.
      */
      lastpassres = 0.0 - 2.0 * prec;
      if(lastpassres > -1.0) { lastpassres = -1.0; }
      while((cc) && (mc) && (!(mec))) {		$? ". segs = %ld", segs
        passres = 0.0;
	if(dk3xsp_calculate(pseg, 0.0, 0)) {
	  lastx = dk3xsp_get_x(pseg);
	  lasty = dk3xsp_get_y(pseg);
	  dsegs = dk3xsp_l_to_d(segs);
	  for(i = 1L; i <= segs; i++) {
	    if(i == segs) {
	      t = 1.0;
	    } else {
	      t = dk3xsp_l_to_d(i) / dsegs;
	    }
	    if(dk3xsp_calculate(pseg, t, 0)) {
	      x = dk3xsp_get_x(pseg);
	      y = dk3xsp_get_y(pseg);
	    } else {
	      cc = 0;
	      mec = DK3_ERROR_MATH_OVERFLOW;
	    }
	    if(cc) {
	      deltax = dk3ma_d_sub_ok(x, lastx, &mec);
	      deltay = dk3ma_d_sub_ok(y, lasty, &mec);
	      passres = dk3ma_d_add_ok(
	        passres,
	        sqrt(
	          dk3ma_d_add_ok(
		    dk3ma_d_mul_ok(deltax, deltax, &mec),
		    dk3ma_d_mul_ok(deltay, deltay, &mec),
		    &mec
		  )
	        ),
	        &mec
	      );
	      lastx = x; lasty = y;
	    }
	  }			$? ". passres = %lg", passres
	  if(fabs(passres - lastpassres) < prec) {	$? ". success"
	    mc = 0;		/* Finished successfully. */
	    back = passres;
	  } else {	$? ". diff = %lg", fabs(passres - lastpassres)
	    lastpassres = passres;
	    segs = segs * 2L;
	    if(segs > DK3XSP_MAX_ITERATION_SEGMENTS) {	$? ". too many passes"
	      cc = 0;
	      if(ec) { *ec = DK3_ERROR_ITERATION; }
	    }
	  }
	} else {
	  cc = 0;
	  mec = DK3_ERROR_MATH_OVERFLOW;
	}
      }
    }
    if(mec) { back = -1.0; if(ec) { *ec = mec; } }
  } else {
    if(ec) { *ec = DK3_ERROR_INVALID_ARGS; }
  } $? "- dk3xsp_segment_length %lg", back
  return back;
}



double
dk3xsp_segment_partial_length(
  dk3_xspline_segment_t		*pseg,
  double			 tstart,
  int				 dir,
  double			 prec,
  int				*ec
)
{
  double		 back = -1.0;
  double		 lastpassres;		/* Result of previous pass. */
  double		 passres;		/* Result of current pass. */
  double		 deltax;		/* Delta x between points. */
  double		 deltay;		/* Delta y between points. */
  double		 lastx;			/* Previous x. */
  double		 lasty;			/* Previous y. */
  double		 x = 0.0;		/* Current x. */
  double		 y = 0.0;		/* Current y. */
  double		 t;			/* Current t. */
  double		 dsegs;			/* Double value of segs. */
  double		 tdiff;			/* 1.0 - tstart. */
  long			 segs = 64L;		/* Number of segments. */
  long			 i;			/* Current segment. */
  int			 mec = 0;		/* Mathematical error code. */
  int			 mc = 1;		/* Flag: Must continue. */
  int			 cc = 1;		/* Flag: Can continue. */
  $? "+ dk3xsp_segment_partial_length t=%lg prec=%lg dir=%d", tstart, prec,dir
  if((pseg) && (prec > 0.0) && (tstart >= 0.0) && (tstart<= 1.0)) { $? ". args"
    if(dk3xsp_is_null((pseg->b).s) && dk3xsp_is_null((pseg->c).s)) {
      /*
      	Simple line.
      */
      $? ". simple line"
      if(dk3xsp_calculate(pseg, tstart, 0)) {
        x = dk3xsp_get_x(pseg);
	y = dk3xsp_get_y(pseg);
      } else {
        cc = 0;
	mec = DK3_ERROR_MATH_OVERFLOW;
      }
      if(dir) {
        deltax = dk3ma_d_sub_ok((pseg->b).x, x, &mec);
	deltay = dk3ma_d_sub_ok((pseg->b).y, y, &mec);
      } else {
        deltax = dk3ma_d_sub_ok((pseg->c).x, x, &mec);
	deltay = dk3ma_d_sub_ok((pseg->c).y, y, &mec);
      }
      back = sqrt(
        dk3ma_d_add_ok(
	  dk3ma_d_mul_ok(deltax, deltax, &mec),
	  dk3ma_d_mul_ok(deltay, deltay, &mec),
	  &mec
	)
      );
    } else {
      /*
      	Spline calculations really necessary.
      */
      $? ". real spline"
      if(dir) {			$? ". from 0 to t"
        lastpassres = 0.0 - 2.0 * prec;
	if(lastpassres > -1.0) { lastpassres = -1.0; }
	while((cc) && (mc) && (!(mec))) {
	  passres = 0.0;
	  dsegs = dk3xsp_l_to_d(segs);
	  if(dk3xsp_calculate(pseg, 0.0, 0)) {
	    lastx = dk3xsp_get_x(pseg);
	    lasty = dk3xsp_get_y(pseg);
	    for(i = 1L; i <= segs; i++) {
	      if(i == segs) {
	        t = tstart;
	      } else {
	        t = (tstart * dk3xsp_l_to_d(i)) / dsegs;
	      } $? ". t=%lg  (%ld/%ld)", t, i, segs
	      if(dk3xsp_calculate(pseg, t, 0)) {
	        $? ". lastx = %lg", lastx
	        $? ". lasty = %lg", lasty
	        x = dk3xsp_get_x(pseg); $? ". x = %lg", x
	        y = dk3xsp_get_y(pseg); $? ". y = %lg", y
	        deltax = dk3ma_d_sub_ok(x, lastx, &mec); $? ". deltax = %lg", deltax
	        deltay = dk3ma_d_sub_ok(y, lasty, &mec); $? ". deltay = %lg", deltay
	        $? ". passres = %lg", passres
	        passres = dk3ma_d_add_ok(
	          passres,
		  sqrt(
		    dk3ma_d_add_ok(
		      dk3ma_d_mul_ok(deltax, deltax, &mec),
		      dk3ma_d_mul_ok(deltay, deltay, &mec),
		      &mec
		    )
		  ),
		  &mec
	        ); $? ". passres = %lg", passres
	        lastx = x; lasty = y;
	      } else {
	        cc = 0;
	        mec = DK3_ERROR_MATH_OVERFLOW;
	      }
	    }
	    if(fabs(passres - lastpassres) < prec) {
	      mc = 0;
	      back = passres;
	    } else {
	      lastpassres = passres;
	      segs = segs * 2L;
	      if(segs > DK3XSP_MAX_ITERATION_SEGMENTS) {
	        cc = 0;
	        if(ec) { *ec = DK3_ERROR_ITERATION; }
	      }
	    }
	  } else {
	    cc = 0;
	    mec = DK3_ERROR_MATH_OVERFLOW;
	  }
	}
      } else {			$? ". from t to 1"
        tdiff = 1.0 - tstart;
	lastpassres = 0.0 - 2.0 * prec;
	if(lastpassres > -1.0) { lastpassres = -1.0; }
	while((cc) && (mc) && (!(mec))) {
	  passres = 0.0;
	  dsegs = dk3xsp_l_to_d(segs);
	  if(dk3xsp_calculate(pseg, tstart, 0)) {
	    lastx = dk3xsp_get_x(pseg);
	    lasty = dk3xsp_get_y(pseg);
	    for(i = 1L; i <= segs; i++) {
	      if(i == segs) {	$? ". t = 1"
		t = 1.0;
	      } else {	
	        t = tstart + ((tdiff * dk3xsp_l_to_d(i)) / dsegs);
	      }
	      $? ". t = %lg", t
	      if(dk3xsp_calculate(pseg, t, 0)) {
		x = dk3xsp_get_x(pseg);
		y = dk3xsp_get_y(pseg);
	      } else {		$? "! calculation"
		cc = 0;
		mec = DK3_ERROR_MATH_OVERFLOW;
	      }
	      if(cc) {
	        $? ". lastx = %lg", lastx
		$? ". lasty = %lg", lasty
		$? ". x = %lg", x
		$? ". y = %lg", y
	        deltax = dk3ma_d_sub_ok(x, lastx, &mec);
		deltay = dk3ma_d_sub_ok(y, lasty, &mec);
		$? ". deltax = %lg", deltax
		$? ". deltay = %lg", deltay
		$? ". passres = %lg", passres
		passres = dk3ma_d_add_ok(
		  passres,
		  sqrt(
		    dk3ma_d_add_ok(
		      dk3ma_d_mul_ok(deltax, deltax, &mec),
		      dk3ma_d_mul_ok(deltay, deltay, &mec),
		      &mec
		    )
		  ),
		  &mec
		); $? ". passres = %lg", passres
		lastx = x; lasty = y;
	      }
	    }
	  } else {
	    cc = 0;
	    mec = DK3_ERROR_MATH_OVERFLOW;
	  }
	  if(fabs(passres - lastpassres) < prec) {
	    mc = 0;
	    back = passres;
	  } else {
	    lastpassres = passres;
	    segs = segs * 2L;
	    if(segs > DK3XSP_MAX_ITERATION_SEGMENTS) {
	      cc = 0;
	      if(ec) { *ec = DK3_ERROR_ITERATION; }
	    }
	  }
	}
      }
    }
    if(mec) { back = -1.0; if(ec) { *ec = mec; } }
  } else {
    if(ec) { *ec = DK3_ERROR_INVALID_ARGS; }
  }
  $? "- dk3xsp_segment_partial_length %lg", back
  return back;
}



$*
Length calculations for entire splines.
$*



double
dk3xsp_open_spline_length(
  dk3_fig_spline_point_t const	*points,
  size_t			 np,
  double			 prec,
  double			*fslc,
  int				 cb,
  int				*ec
)
{
  dk3_xspline_segment_t		 seg;
  double			 val;
  double			 back = 0.0;
  size_t			 i;
  int				 mec = 0;
  $? "+ dk3xsp_open_spline_length prec=%lg", prec
  if((points) && (1 < np) && (prec > 0.0)) {
    dk3xsp_reset(&seg);
    if(cb) { dk3xsp_set_cb(&seg, 1); }
    for(i = 0; ((i < (np - 1)) && (back > -1.0)); i++) {
      $? ". segment %u (%u points)", (unsigned)i, (unsigned)np
      val = -1.0;
      if(fslc) {
        val = fslc[i];
      }
      if(val < -0.5) {
        dk3xsp_set(
          &seg,
	  ((i > 1) ? (&(points[i-1])) : NULL),
	  &(points[i]),
	  &(points[i+1]),
	  ((i < (np - 2)) ? &(points[i+2]) : NULL)
        );
        val = dk3xsp_segment_length(&seg, (prec / (double)(np - 1)), ec);
	if(val >= 0.0) {
	  if(fslc) {
	    fslc[i] = val;
	  }
	}
      }
      if(val < 0.0) {
          back = -1.0; i = np - 1;
      } else {
        back = dk3ma_d_add_ok(back, val, &mec);
	if(mec) {
	  back = -1.0; i = np - 1;
	  if(ec) { *ec = mec; }
	}
      }
    }
  } else {
    if(ec) { *ec = DK3_ERROR_INVALID_ARGS; }
  } $? "- dk3xsp_open_spline_length %lg", back
  return back;
}



/**	Find Index of segment start point.
	@param	tstart	Parameter t.
	@param	ec	Pointer to error code variable, required.
	@return	Start segment index.
*/
static
size_t
dk3xsp_find_segment_start(double tstart, int *ec)
{
#if VERSION_BEFORE_20140809
  unsigned long		ul;
  size_t		back = 0;
  int			mec = 0;
  $? "+ dk3xsp_find_segment_start %lg", tstart
  ul = dk3ma_d_to_ul_ok(floor(tstart), &mec);
  if(mec) {
    *ec = mec;
  } else {
    back = (size_t)ul;
    if((unsigned long)back != ul) {
      *ec = DK3_ERROR_MATH_OVERFLOW;
    }
  } $? "- dk3xsp_find_segment_start %u", (unsigned)back
  return back;
#else
  size_t		back	= 0;
  int			mec	= 0;
  $? "+ dk3xsp_find_segment_start %lg", tstart
  back = dk3ma_d_to_sz_ok(floor(tstart), &mec);
  if (mec) { if (ec) { *ec = mec; } }
  $? "- dk3xsp_find_segment_start %u (ec=%d)", (unsigned)back, mec
  return back;
#endif
}



int
dk3xsp_calculate_position(
  dk3_xspline_segment_t		*pseg,
  dk3_fig_spline_point_t const	*points,
  size_t			 np,
  double			 t,
  int				 cb,
  int				 der,
  int				*ec
)
{
  size_t	 i = 0;
  int		 back = 0;
  int		 mec = 0;
  if((pseg) && (points) && (np > 1) && ((t >= 0.0) && (t <= (double)(np - 1)))) 
  {
    i = dk3xsp_find_segment_start(t, &mec);
    if(!(mec)) {
      dk3xsp_reset(pseg);
      dk3xsp_set_cb(pseg, cb);
      if(i <= (np - 2)) {
        /* Not final point */
        dk3xsp_set(
	  pseg,
	  ((i > 0) ? &(points[i - 1]) : NULL),
	  &(points[i]),
	  &(points[i+1]),
	  ((i < (np - 2)) ? &(points[i + 2]) : NULL)
	);
	back = dk3xsp_calculate(pseg, (t - (double)i), der);
      } else {
        /* In final point */
	dk3xsp_set(
	  pseg,
	  ((np > 2) ? &(points[np - 3]): NULL),
	  &(points[np - 2]),
	  &(points[np - 1]),
	  NULL
	);
	back = dk3xsp_calculate(pseg, 1.0, der);
      }
    } else {
      if(ec) { *ec = mec; }
    }
  }
  else
  {
    if(ec) { *ec = DK3_ERROR_INVALID_ARGS;}
  }
  return back;
}



double
dk3xsp_open_spline_partial_length(
  dk3_fig_spline_point_t const	*points,
  size_t			 np,
  double			 prec,
  double			 tstart,
  double			*fslc,
  int				 dir,
  int				 cb,
  int				*ec
)
{
  dk3_xspline_segment_t		 seg;
  double			 back = -1.0;
  double			 val;
  size_t			 startseg;
  size_t			 i;
  int				 mec = 0;
  $? "+ dk3xsp_open_spline_partial_length t=%lg prec=%lg", tstart, prec
  if((points) && (1 < np) && (prec > 0.0)) {
    if((tstart >= 0.0) && (tstart <= ((double)(np - 1)))) {
      dk3xsp_reset(&seg);
      if(cb) { dk3xsp_set_cb(&seg, 1); }
      startseg = dk3xsp_find_segment_start(tstart, &mec);
      if(mec) {				$? "! mec"
        if(ec) { *ec = mec; }
      } else {				$? ". startseg = %u", (unsigned)startseg
        if(startseg == (np - 1)) {	$? ". calculation for final point"
	  /*
	  	tstart points to the final point.
		So we have 0 for forward length or the full
		spline length for backward length.
	  */
	  if(dir) {	$? ". full length for left side of end point"
	    back = dk3xsp_open_spline_length(points, np, prec, fslc, cb, ec);
	  } else {	$? ". zero for end to end"
	    back = 0.0;
	  }
	} else {			$? ". not final point"
	  if(tstart > 0.0) {		$? ". point within the spline"
	    /*
	    	tstart points to somewhere in the spline.
		We have to calculate one segment partially, all other
		segments completely.
	    */
	    dk3xsp_set(
	      &seg,
	      ((startseg > 0) ? &(points[startseg-1]) : NULL),
	      &(points[startseg]),
	      &(points[startseg+1]),
	      ((startseg < (np - 2)) ? &(points[startseg+2]) : NULL)
	    );
	    val = dk3xsp_segment_partial_length(
	      &seg,(tstart - (double)startseg),dir,(prec / (double)(np - 1)),ec
	    );
	    if(val >= 0.0) {	$? ". partial segment length: %lg", val
	      back = val;
	      /*
	      	Partial segment calculated successfully.
		Now calculate other segments.
	      */
	      if(dir) {		$? ". full segments before partial one"
	        /*
			Calculate all segments before the partial one.
		*/
	        for(i = 0; ((i < startseg) && (back > -1.0)); i++) {
		  val = -1.0;	$? ". segment %u", (unsigned)i
		  if(fslc) {
		    val = fslc[i];
		  }
		  if(val < -0.5) {
		    dk3xsp_set(
		      &seg,
		      ((i > 0) ? &(points[i-1]) : NULL),
		      &(points[i]),
		      &(points[i+1]),
		      ((i < (np - 2)) ? &(points[i+2]) : NULL)
		    );
		    val = dk3xsp_segment_length(
		      &seg, (prec / (double)(np - 1)), ec
		    );
		    if(val >= 0.0) {
		      if(fslc) {
		        fslc[i] = val;
		      }
		    }
		  }
		  if(val >= 0.0) { $? ". segment %u length %lg",(unsigned)i,val
		    back = dk3ma_d_add_ok(back, val, &mec);
		    if(mec) {
		      back = -1.0; i = np;
		      if(ec) { *ec = mec; }
		    }
		  } else {
		    back = -1.0; i = np;
		  }
		}
	      } else {		$? ". full segments after partial one"
	        /*
			Calculate all segments after the partial one.
		*/
	        for(i = (startseg + 1); ((i < (np - 1)) && (back > -1)); i++) {
		  val = -1.0;	$? ". segment %u", (unsigned)i
		  if(fslc) {
		    val = fslc[i];
		  }
		  if(val < -0.5) {
		    dk3xsp_set(
		      &seg,
		      ((i > 0) ? &(points[i-1]) : NULL),
		      &(points[i]),
		      &(points[i+1]),
		      ((i < (np - 2)) ? &(points[i+2]) : NULL)
		    );
		    val = dk3xsp_segment_length(
		      &seg, (prec / (double)(np - 1)), ec
		    );
		    if(val >= 0.0) {
		      if(fslc) {
		        fslc[i] = val;
		      }
		    }
		  }
		  if(val >= 0.0) { $? ". segment %u length %lg",(unsigned)i,val
		    back = dk3ma_d_add_ok(back, val, &mec);
		    if(mec) {
		      back = -1.0; i = np;
		      if(ec) { *ec = mec; }
		    }
		  } else {			$? "! val"
		    back = -1.0; i = np;
		  }
		}
	      }	$? ". full segments done"
	    } else {				$? "! val"
	    }
	  } else {				$? ". start point"
	    /*
	    	tstart points to the start point.
		So we have the full spline length for forward length,
		0 for backward length.
	    */
	    if(0 == dir) {	$? ". full length if right from start"
	      back = dk3xsp_open_spline_length(points, np, prec, fslc, cb, ec);
	    } else {		$? ". zero from start to start"
	      back = 0.0;
	    }
	  }
	}
      }
    } else {
      if(ec) { *ec = DK3_ERROR_INVALID_ARGS; }
    }
  } else {
    if(ec) { *ec = DK3_ERROR_INVALID_ARGS; }
  } $? "- dk3xsp_open_spline_partial_length %lg", back
  return back;
}



double
dk3xsp_get_t_for_distance(
  dk3_fig_spline_point_t const	*points,
  size_t			 np,
  double			 dist,
  double			 prec,
  double			*fslc,
  int				 dir,
  int				 cb,
  int				*ec
)
{
  double		 back = -1.0;
  double		 xl;		/* Left x. */
  double		 xr;		/* Right x. */
  double		 xc;		/* Center x. */
  double		 yl;		/* Left y. */
  double		 yr;		/* Right y. */
  double		 yc;		/* Center y. */
  double		 l;		/* Spline length. */
  double		 fprec;		/* factor * precision. */
  double		 ayl = 0.0;	/* abs yl. */
  double		 ayr = 0.0;	/* abs yr. */
  unsigned long		 passno = 0UL;	/* Current pass. */
  int			 mec = 0;	/* Error code variable. */
  int			 cc  = 1;	/* Flag: Can continue. */
  int			 mc  = 1;	/* Flag: Must continue. */
  $? "+ dk3xsp_get_t_for_distance d=%lg prec=%lg", dist, prec
  fprec = 0.2 * prec;
  if((points) && (1 < np) && (prec > 0.0) && (dist >= 0.0)) {
    l = dk3xsp_open_spline_length(points, np, fprec, fslc, cb, &mec);
    $? ". mec = %d", mec
    if(l > dist) {
      xl = 0.0;
      xr = (double)(np - 1);
      if(dir) {
        yl = 0.0 - dist;
	yr = dk3ma_d_sub_ok(l, dist, &mec);	$? ". mec = %d", mec
      } else {
        yl = dk3ma_d_sub_ok(l, dist, &mec);	$? ". mec = %d", mec
	yr = 0.0 - dist;
      }
      while((cc) && (mc) && (!(mec))) {
        xc = 0.5 * dk3ma_d_add_ok(xl, xr, &mec);
	$? ". mec = %d", mec
	yc = dk3ma_d_sub_ok(
	  dk3xsp_open_spline_partial_length(
	    points, np, fprec, xc, fslc, dir, cb, ec
	  ),
	  dist,
	  &mec
	);	$? ". mec = %d", mec
	$? ". xl = %lg   yl = %lg", xl, yl
	$? ". xr = %lg   yr = %lg", xr, yr
	$? ". xc = %lg   yc = %lg", xc, yc
	if(dir) {
	  if(yc > 0.0) {
	    xr = xc; yr = yc;
	  } else {
	    xl = xc; yl = yc;
	  }
	} else {
	  if(yc > 0.0) {
	    xl = xc; yl = yc;
	  } else {
	    xr = xc; yr = yc;
	  }
	}
	ayl = fabs(yl);
	ayr = fabs(yr);
	if((ayl < prec) && (ayr < prec)) {
	  mc = 0;
	} else {
	  passno++;	$? ". pass %lu", passno
	  if(passno > 1000UL) {		$? "! iterations"
	    mec = DK3_ERROR_ITERATION;
	  } else {
	  }
	}
      }
      if((cc) && (!(mec))) {	$? ". SUCCESS %lu", passno
	back = xl + ((xr - xl) * ayl)/(ayl + ayr);
      }
      if(mec) {			$? "! mec = %d", mec
        back = -1.0;
	if(ec) { *ec = mec; }
      }
    } else {
      if(ec) { *ec = DK3_ERROR_INVALID_ARGS; }
    }
  } else {
    if(ec) { *ec = DK3_ERROR_INVALID_ARGS; }
  } $? "- dk3xsp_get_t_for_distance %lg  mec=%d", back, mec
  return back;
}



double
dk3xsp_direction(
  dk3_fig_spline_point_t const	*points,
  size_t			 np,
  double			 t,
  int				 cb,
  int				*ec
)
{
  dk3_xspline_segment_t	 seg;
  double		 back = -10.0 * M_PI;
  double		 dxdt;
  double		 dydt;
  double		 delta;
  double		 tleft;
  double		 tright;
  double		 xleft;
  double		 xright;
  double		 yleft;
  double		 yright;
  int			 mec = 0;
  $? "+ dk3xsp_direction %lg (%u)", t, (unsigned)(np - 1)
  if((points) && (1 < np)) {
    if((t >= 0.0) && (t <= ((double)(np - 1)))) {
      dk3xsp_reset(&seg);
      dk3xsp_set_cb(&seg, ((cb) ? 1 : 0));
      if(dk3xsp_calculate_position(&seg, points, np, t, cb, 1, &mec)) {
        dxdt = dk3xsp_get_dxdt(&seg);
	dydt = dk3xsp_get_dydt(&seg);
	if((fabs(dxdt) > 1.0e-6) || (fabs(dydt) > 1.0e-6)) {
	  back = dk3ma_d_atan2(dydt, dxdt);
	} else {			$? ". must iterate"
	  delta = 0.001;
	  tleft = tright = t;
	  if(t > delta) { tleft = t - delta; }
	  if(t < (((double)(np - 1)) - delta)) { tright = t + delta; }
	  if(dk3xsp_calculate_position(&seg, points, np, tleft, cb, 0, &mec)) {
	  xleft = dk3xsp_get_x(&seg);
	  yleft = dk3xsp_get_y(&seg);
	  if(dk3xsp_calculate_position(&seg, points, np, tright, cb, 0, &mec)) {
	    xright = dk3xsp_get_x(&seg);
	    yright = dk3xsp_get_y(&seg);
	    dxdt = dk3ma_d_sub_ok(xright, xleft, &mec);
	    dydt = dk3ma_d_sub_ok(yright, yleft, &mec);
	    back = dk3ma_d_atan2(dydt, dxdt);
	  } else {		$? "! calculation"
	  }
	  } else {		$? "! calculation"
	  }
	}
      } else {				$? "! calculation failed"
      }
      if(mec) { if(ec) { *ec = mec; } }
    } else {
      if(ec) { *ec = DK3_ERROR_INVALID_ARGS; }
    }
  } else {
    if(ec) { *ec = DK3_ERROR_INVALID_ARGS; }
  } $? "- dk3xsp_direction %lg", back
  return back;
}



double
dk3xsp_direction_and_position(
  double			*xpos,
  double			*ypos,
  dk3_fig_spline_point_t const	*points,
  size_t			 np,
  double			 t,
  int				 cb,
  int				*ec
)
{
  dk3_xspline_segment_t	 seg;
  double		 back = -10.0 * M_PI;
  double		 dxdt;
  double		 dydt;
  double		 delta;
  double		 tleft;
  double		 tright;
  double		 xleft;
  double		 xright;
  double		 yleft;
  double		 yright;
  int			 mec = 0;
  $? "+ dk3xsp_direction %lg (%u)", t, (unsigned)(np - 1)
  if((points) && (1 < np) && (xpos) && (ypos)) {
    if((t >= 0.0) && (t <= ((double)(np - 1)))) {
      dk3xsp_reset(&seg);
      dk3xsp_set_cb(&seg, ((cb) ? 1 : 0));
      if(dk3xsp_calculate_position(&seg, points, np, t, cb, 1, &mec)) {
	*xpos = dk3xsp_get_x(&seg);
	*ypos = dk3xsp_get_y(&seg);
        dxdt = dk3xsp_get_dxdt(&seg);
	dydt = dk3xsp_get_dydt(&seg);	$? ". dxdt=%lg dydt=%lg", dxdt, dydt
	if((fabs(dxdt) > 1.0e-6) || (fabs(dydt) > 1.0e-6)) {
	  back = dk3ma_d_atan2(dydt, dxdt);	$? ". back = %lg", back
	} else {			$? ". must iterate"
	  delta = 0.001;
	  tleft = tright = t;
	  if(t > delta) { tleft = t - delta; }
	  if(t < (((double)(np - 1)) - delta)) { tright = t + delta; }
	  if(dk3xsp_calculate_position(&seg, points, np, tleft, cb, 0, &mec)) {
	  xleft = dk3xsp_get_x(&seg);
	  yleft = dk3xsp_get_y(&seg);
	  if(dk3xsp_calculate_position(&seg, points, np, tright, cb, 0, &mec)) {
	    xright = dk3xsp_get_x(&seg);
	    yright = dk3xsp_get_y(&seg);
	    dxdt = dk3ma_d_sub_ok(xright, xleft, &mec);
	    dydt = dk3ma_d_sub_ok(yright, yleft, &mec);
	    back = dk3ma_d_atan2(dydt, dxdt);
	  } else {		$? "! calculation"
	  }
	  } else {		$? "! calculation"
	  }
	}
      } else {				$? "! calculation failed"
      }
      if(mec) { if(ec) { *ec = mec; } }
    } else {
      if(ec) { *ec = DK3_ERROR_INVALID_ARGS; }
    }
  } else {
    if(ec) { *ec = DK3_ERROR_INVALID_ARGS; }
  } $? "- dk3xsp_direction %lg", ((180.0 * back)/M_PI)
  return back;
}