/* tense.f -- translated by f2c (version 20090411).
   You must link the resulting object file with libf2c:
	on Microsoft Windows system, link with libf2c.lib;
	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
	or, if you install libf2c.a in a standard place, with -lf2c -lm
	-- in that order, at the end of the command line, as in
		cc *.o -lf2c -lm
	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

		http://www.netlib.org/f2c/libf2c.zip
*/

#include "f2c.h"

integer curv1_(integer *n, doublereal *x, doublereal *y, doublereal *slp1, 
	doublereal *slpn, doublereal *yp, doublereal *temp, doublereal *sigma,
	 doublereal *xf, doublereal *yf, integer *nf)
{
    /* System generated locals */
    integer ret_val, i__1, i__2;

    /* Builtin functions */
    double exp(doublereal);

    /* Local variables */
    static integer i__, j;
    static doublereal t, c1, c2, c3;
    static integer i1, nm1;
    static doublereal dx1, dx2;
    static integer np1;
    static doublereal del1, del2;
    static integer ibak;
    static doublereal deln, dels, exps, diag1, diag2, delx1, delx2, slpp1, 
	    exps1, delx12, delnn, sinhs, slppn, delnm1, sinhd1, sinhd2, 
	    diagin, spdiag, sigmap, sinhin;

/* --	generates spline under tension for monotonic x */
/* --	NOTE THAT THIS IS A FUNCTION IN ORDER TO RETURN SUCCESS OR FAILURE */
    /* Parameter adjustments */
    --temp;
    --yp;
    --y;
    --x;
    --yf;
    --xf;

    /* Function Body */
    if (*sigma == 0.f) {
/* 	TYPE *,'tension cannot be 0' */
	ret_val = 0;
	return ret_val;
    }
    nm1 = *n - 1;
    np1 = *n + 1;
    delx1 = x[2] - x[1];
    dx1 = (y[2] - y[1]) / delx1;
    if (*sigma < 0.) {
	goto L50;
    }
    slpp1 = *slp1;
    slppn = *slpn;
L10:
    sigmap = abs(*sigma) * (real) (*n - 1) / (x[*n] - x[1]);
/* 	FOR INTERACTIVE USE, WE SLOW THIS ROUTINE DOWN BY DOING SOME */
/* 	EXTRA CHECKING FOR THE USER, THIS COULD PREVENT MUCH TIME LOST */
/* 	FROM A CRASH */
    dx2 = 0.f;
    i__1 = *n;
    for (i__ = 2; i__ <= i__1; ++i__) {
	delx2 = x[i__] - x[i__ - 1];
	if (delx2 <= 0.) {
	    goto L990;
	}
/* ERROR, EXIT */
	dels = sigmap * delx2;
	if (dels > dx2) {
	    dx2 = dels;
	}
/* FIND MAX */
/* L12: */
    }
    if (dx2 > 50.f) {
/* -- IF IT IS TOO BIG, WE GOT TROUBLE; WE WILL AUTO REDUCE TO FIT */
/* 	TYPE	*,' warning, tension reduced to avoid overflows !' */
	sigmap = sigmap * 50.f / dx2;
    }
    dels = sigmap * delx1;
    exps = exp(dels);
    sinhs = (exps - 1.f / exps) * .5f;
    sinhin = 1.f / (delx1 * sinhs);
    diag1 = sinhin * (dels * .5f * (exps + 1.f / exps) - sinhs);
    diagin = 1.f / diag1;
    yp[1] = diagin * (dx1 - slpp1);
    spdiag = sinhin * (sinhs - dels);
    temp[1] = diagin * spdiag;
    if (*n == 2) {
	goto L30;
    }
    i__1 = nm1;
    for (i__ = 2; i__ <= i__1; ++i__) {
	delx2 = x[i__ + 1] - x[i__];
	dx2 = (y[i__ + 1] - y[i__]) / delx2;
	dels = sigmap * delx2;
	exps = exp(dels);
	sinhs = (exps - 1.f / exps) * .5f;
	sinhin = 1.f / (delx2 * sinhs);
	diag2 = sinhin * (dels * ((exps + 1.f / exps) * .5f) - sinhs);
	diagin = 1.f / (diag1 + diag2 - spdiag * temp[i__ - 1]);
	yp[i__] = diagin * (dx2 - dx1 - spdiag * yp[i__ - 1]);
	spdiag = sinhin * (sinhs - dels);
	temp[i__] = diagin * spdiag;
	dx1 = dx2;
	diag1 = diag2;
/* L20: */
    }
L30:
    diagin = 1.f / (diag1 - spdiag * temp[nm1]);
    yp[*n] = diagin * (slppn - dx2 - spdiag * yp[nm1]);
    i__1 = *n;
    for (i__ = 2; i__ <= i__1; ++i__) {
	ibak = np1 - i__;
	yp[ibak] -= temp[ibak] * yp[ibak + 1];
/* L40: */
    }
    goto L200;
L50:
    if (*n == 2) {
	goto L60;
    }
    delx2 = x[3] - x[2];
    delx12 = x[3] - x[1];
    c1 = -(delx12 + delx1) / delx12 / delx1;
    c2 = delx12 / delx1 / delx2;
    c3 = -delx1 / delx12 / delx2;
    slpp1 = c1 * y[1] + c2 * y[2] + c3 * y[3];
    deln = x[*n] - x[nm1];
    delnm1 = x[nm1] - x[*n - 2];
    delnn = x[*n] - x[*n - 2];
    c1 = (delnn + deln) / delnn / deln;
    c2 = -delnn / deln / delnm1;
    c3 = deln / delnn / delnm1;
    slppn = c3 * y[*n - 2] + c2 * y[nm1] + c1 * y[*n];
    goto L10;
L60:
    yp[1] = 0.f;
    yp[2] = 0.f;
L200:
    i1 = 2;
    i__1 = *nf;
    for (j = 1; j <= i__1; ++j) {
	t = xf[j];
L210:
	i__2 = *n;
	for (i__ = i1; i__ <= i__2; ++i__) {
	    if (x[i__] - t <= 0.) {
		goto L220;
	    } else {
		goto L230;
	    }
L220:
	    ;
	}
	i__ = *n;
L230:
	if (x[i__ - 1] <= t || t <= x[1]) {
	    goto L240;
	}
	i1 = 2;
	goto L210;
L240:
	del1 = t - x[i__ - 1];
	del2 = x[i__] - t;
	dels = x[i__] - x[i__ - 1];
	exps1 = exp(sigmap * del1);
	sinhd1 = (exps1 - 1.f / exps1) * .5f;
	exps = exp(sigmap * del2);
	sinhd2 = (exps - 1.f / exps) * .5f;
	exps = exps1 * exps;
	sinhs = (exps - 1.f / exps) * .5f;
	yf[j] = (yp[i__] * sinhd1 + yp[i__ - 1] * sinhd2) / sinhs + ((y[i__] 
		- yp[i__]) * del1 + (y[i__ - 1] - yp[i__ - 1]) * del2) / dels;
	i1 = i__;
/* L300: */
    }
    ret_val = 1;
    return ret_val;
L990:
/* 	TYPE *,'independent vector not monotonically increasing' */
/* ERROR */
    ret_val = 0;
    return ret_val;
} /* curv1_ */

/* =========================================================================== */
integer kurv1_(integer *n, doublereal *x, doublereal *y, doublereal *slp1, 
	doublereal *slpn, doublereal *xp, doublereal *yp, doublereal *temp, 
	doublereal *sigma, doublereal *t, doublereal *xs, doublereal *ys, 
	integer *nf)
{
    /* Initialized data */

    static doublereal degrad = .017453292;

    /* System generated locals */
    integer ret_val, i__1, i__2;
    doublereal d__1, d__2;

    /* Builtin functions */
    double sqrt(doublereal), atan2(doublereal, doublereal), cos(doublereal), 
	    sin(doublereal), exp(doublereal);

    /* Local variables */
    static integer i__, k;
    static doublereal s, c1, c2, c3;
    static integer i1;
    static doublereal tn, xq, sx, sy;
    static integer nm1;
    static doublereal dx1, dy1, dx2, dy2, sum, del1, del2, deln, dels, delx, 
	    dely, exps, diag1, diag2, dels1, dels2, delx1, dely1, delx2, 
	    dely2, slpp1, exps1, dels12, delnn, sinhs, slppn, delnm1, sinhd1, 
	    sinhd2, diagin, spdiag, sigmap, sinhin;

/* --	generates spline under tension for any open (x,y) curve */
/* --	NOTE THAT THIS IS A FUNCTION IN ORDER TO RETURN SUCCESS OR FAILURE */
    /* Parameter adjustments */
    --temp;
    --yp;
    --xp;
    --y;
    --x;
    --ys;
    --xs;
    --t;

    /* Function Body */
    ret_val = 1;
    if (*sigma == 0.f) {
/* 	TYPE *,'tension cannot be 0' */
	ret_val = 0;
	return ret_val;
    }
    nm1 = *n - 1;
    delx1 = x[2] - x[1];
    dely1 = y[2] - y[1];
    dels1 = sqrt(delx1 * delx1 + dely1 * dely1);
    dx1 = delx1 / dels1;
    dy1 = dely1 / dels1;
    s = dels1;
/* --	check if slopes were input (yes if sigma>0) */
    if (*sigma > 0.) {
/* --	convert input slopes from degress to radians */
	slpp1 = *slp1 * degrad;
	slppn = *slpn * degrad;
    } else {
/* determine slopes */
	if (*n == 2) {
/* --	if only 2 points and no slopes, we have to assume a straight line */
	    xp[1] = 0.f;
	    xp[2] = 0.f;
	    yp[1] = 0.f;
	    yp[2] = 0.f;
	    sigmap = abs(*sigma) * (real) (*n - 1) / s;
	} else {
/* Computing 2nd power */
	    d__1 = x[3] - x[2];
/* Computing 2nd power */
	    d__2 = y[3] - y[2];
	    dels2 = sqrt(d__1 * d__1 + d__2 * d__2);
	    dels12 = dels1 + dels2;
	    c1 = -(dels12 + dels1) / dels12 / dels1;
	    c2 = dels12 / dels1 / dels2;
	    c3 = -dels1 / dels12 / dels2;
	    sx = c1 * x[1] + c2 * x[2] + c3 * x[3];
	    sy = c1 * y[1] + c2 * y[2] + c3 * y[3];
	    slpp1 = atan2(sy, sx);
	    xq = slpp1 / degrad;
/* Computing 2nd power */
	    d__1 = x[*n - 2] - x[nm1];
/* Computing 2nd power */
	    d__2 = y[*n - 2] - y[nm1];
	    delnm1 = sqrt(d__1 * d__1 + d__2 * d__2);
/* Computing 2nd power */
	    d__1 = x[nm1] - x[*n];
/* Computing 2nd power */
	    d__2 = y[nm1] - y[*n];
	    deln = sqrt(d__1 * d__1 + d__2 * d__2);
	    delnn = delnm1 + deln;
	    c1 = (delnn + deln) / delnn / deln;
	    c2 = -delnn / deln / delnm1;
	    c3 = deln / delnn / delnm1;
	    sx = c3 * x[*n - 2] + c2 * x[nm1] + c1 * x[*n];
	    sy = c3 * y[*n - 2] + c2 * y[nm1] + c1 * y[*n];
	    slppn = atan2(sy, sx);
	}
    }
/* --	slopes are now done */
/* 	TYPE *,'slopes are',SLPP1,SLPPN */
    xp[1] = dx1 - cos(slpp1);
    yp[1] = dy1 - sin(slpp1);
    temp[1] = dels1;
    if (*n > 2) {
	i__1 = nm1;
	for (i__ = 2; i__ <= i__1; ++i__) {
	    delx2 = x[i__ + 1] - x[i__];
	    dely2 = y[i__ + 1] - y[i__];
	    dels2 = sqrt(delx2 * delx2 + dely2 * dely2);
	    dx2 = delx2 / dels2;
	    dy2 = dely2 / dels2;
	    xp[i__] = dx2 - dx1;
	    yp[i__] = dy2 - dy1;
	    temp[i__] = dels2;
	    delx1 = delx2;
	    dely1 = dely2;
	    dels1 = dels2;
	    dx1 = dx2;
	    dy1 = dy2;
	    s += dels1;
	}
    }
    xp[*n] = cos(slppn) - dx1;
    yp[*n] = sin(slppn) - dy1;
    sigmap = abs(*sigma) * (real) (*n - 1) / s;
    dels = sigmap * temp[1];
    exps = exp(dels);
/* 	TYPE *,'DELS,EXPS =',DELS,EXPS */
    sinhs = (exps - 1.f / exps) * .5f;
    sinhin = 1.f / (temp[1] * sinhs);
    diag1 = sinhin * (dels * .5f * (exps + 1.f / exps) - sinhs);
    diagin = 1.f / diag1;
    xp[1] = diagin * xp[1];
    yp[1] = diagin * yp[1];
    spdiag = sinhin * (sinhs - dels);
    temp[1] = diagin * spdiag;
    if (*n > 2) {
	i__1 = nm1;
	for (i__ = 2; i__ <= i__1; ++i__) {
	    dels = sigmap * temp[i__];
	    exps = exp(dels);
/* 	TYPE *,'DELS,EXPS =',DELS,EXPS */
	    sinhs = (exps - 1.f / exps) * .5f;
	    sinhin = 1.f / (temp[i__] * sinhs);
	    diag2 = sinhin * (dels * ((exps + 1.f / exps) * .5f) - sinhs);
	    diagin = 1.f / (diag1 + diag2 - spdiag * temp[i__ - 1]);
	    xp[i__] = diagin * (xp[i__] - spdiag * xp[i__ - 1]);
	    yp[i__] = diagin * (yp[i__] - spdiag * yp[i__ - 1]);
	    spdiag = sinhin * (sinhs - dels);
	    temp[i__] = diagin * spdiag;
	    diag1 = diag2;
	}
    }
    diagin = 1.f / (diag1 - spdiag * temp[nm1]);
    xp[*n] = diagin * (xp[*n] - spdiag * xp[nm1]);
    yp[*n] = diagin * (yp[*n] - spdiag * yp[nm1]);
    for (i__ = *n - 1; i__ >= 1; --i__) {
	xp[i__] -= temp[i__] * xp[i__ + 1];
	yp[i__] -= temp[i__] * yp[i__ + 1];
    }
/* --	done with the setup, now do the interpolations */

    i__ = 2;
    sum = 0.f;
    delx = x[2] - x[1];
    dely = y[2] - y[1];
    dels = sqrt(delx * delx + dely * dely);
    i__1 = *nf;
    for (k = 1; k <= i__1; ++k) {
/* --	loop over the nf input values of T */
	tn = (d__1 = t[k] * s, abs(d__1));
L10:
	if (tn < sum) {
/* --	drop back */
	    --i__;
	    if (i__ >= 2) {
		delx = x[i__] - x[i__ - 1];
		dely = y[i__] - y[i__ - 1];
		dels = sqrt(delx * delx + dely * dely);
		sum -= dels;
		goto L10;
	    } else {
		xs[k] = x[1];
		ys[k] = y[1];
		goto L50;
	    }
	}
	if (tn > sum + dels) {
/* --	go forward */
	    i1 = i__;
	    i__2 = *n;
	    for (i__ = i1; i__ <= i__2; ++i__) {
		delx = x[i__] - x[i__ - 1];
		dely = y[i__] - y[i__ - 1];
		dels = sqrt(delx * delx + dely * dely);
		if (sum + dels > tn) {
		    goto L40;
		}
		sum += dels;
	    }
	    xs[k] = x[*n];
	    ys[k] = y[*n];
	    goto L50;
L40:
	    ;
	}
	del1 = tn - sum;
	del2 = dels - del1;
	exps1 = exp(sigmap * del1);
	sinhd1 = (exps1 - 1.f / exps1) * .5f;
	exps = exp(sigmap * del2);
	sinhd2 = (exps - 1.f / exps) * .5f;
	exps = exps1 * exps;
	sinhs = (exps - 1.f / exps) * .5f;
	xs[k] = (xp[i__] * sinhd1 + xp[i__ - 1] * sinhd2) / sinhs + ((x[i__] 
		- xp[i__]) * del1 + (x[i__ - 1] - xp[i__ - 1]) * del2) / dels;
	ys[k] = (yp[i__] * sinhd1 + yp[i__ - 1] * sinhd2) / sinhs + ((y[i__] 
		- yp[i__]) * del1 + (y[i__ - 1] - yp[i__ - 1]) * del2) / dels;
L50:
	;
    }
    return ret_val;
} /* kurv1_ */

/* ============================================================================= */
integer kurvp1_(integer *n, doublereal *x, doublereal *y, doublereal *xp, 
	doublereal *yp, doublereal *temp, doublereal *sigma, doublereal *t, 
	doublereal *xs, doublereal *ys, integer *nf)
{
    /* System generated locals */
    integer ret_val, i__1, i__2;
    doublereal d__1;

    /* Builtin functions */
    double sqrt(doublereal), exp(doublereal);

    /* Local variables */
    static integer i__, k;
    static doublereal s;
    static integer i1;
    static doublereal tn;
    static integer im1, ip1, nm1;
    static doublereal dx1, dy1;
    static integer np1;
    static doublereal dx2, dy2, sum, del1, del2;
    static integer ibak;
    static doublereal dels, delx, dely, exps, diag1, diag2, dels1, dels2, 
	    delx1, dely1, delx2, dely2, exps1, sinhs, sinhd1, spdig1, sinhd2, 
	    diagin, spdiag, sigmap, sinhin;

/* --	generates spline under tension for any closed loop (x,y) curve */
/* --	NOTE THAT THIS IS A FUNCTION IN ORDER TO RETURN SUCCESS OR FAILURE */
    /* Parameter adjustments */
    --yp;
    --xp;
    --y;
    --x;
    --temp;
    --ys;
    --xs;
    --t;

    /* Function Body */
    ret_val = 1;
    if (*sigma == 0.f) {
/* 	TYPE *,'tension cannot be 0' */
	ret_val = 0;
	return ret_val;
    }
    nm1 = *n - 1;
    np1 = *n + 1;
    delx1 = x[2] - x[1];
    dely1 = y[2] - y[1];
    dels1 = sqrt(delx1 * delx1 + dely1 * dely1);
    dx1 = delx1 / dels1;
    dy1 = dely1 / dels1;
    xp[1] = dx1;
    yp[1] = dy1;
    temp[1] = dels1;
    s = dels1;
    i__1 = *n;
    for (i__ = 2; i__ <= i__1; ++i__) {
	ip1 = i__ + 1;
	if (i__ == *n) {
	    ip1 = 1;
	}
	delx2 = x[ip1] - x[i__];
	dely2 = y[ip1] - y[i__];
	dels2 = sqrt(delx2 * delx2 + dely2 * dely2);
	dx2 = delx2 / dels2;
	dy2 = dely2 / dels2;
	xp[i__] = dx2 - dx1;
	yp[i__] = dy2 - dy1;
	temp[i__] = dels2;
	delx1 = delx2;
	dely1 = dely2;
	dels1 = dels2;
	dx1 = dx2;
	dy1 = dy2;
	s += dels1;
    }
    xp[1] -= dx1;
    yp[1] -= dy1;
    sigmap = abs(*sigma) * (real) (*n) / s;
    dels = sigmap * temp[*n];
    exps = exp(dels);
    sinhs = (exps - 1.f / exps) * .5f;
    sinhin = 1.f / (temp[*n] * sinhs);
    diag1 = sinhin * (dels * .5f * (exps + 1.f / exps) - sinhs);
    diagin = 1.f / diag1;
    spdig1 = sinhin * (sinhs - dels);
    spdiag = 0.f;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	dels = sigmap * temp[i__];
	exps = exp(dels);
	sinhs = (exps - 1.f / exps) * .5f;
	sinhin = 1.f / (temp[i__] * sinhs);
	diag2 = sinhin * (dels * ((exps + 1.f / exps) * .5f) - sinhs);
	if (i__ == *n) {
	    goto L30;
	}
	diagin = 1.f / (diag1 + diag2 - spdiag * temp[i__ - 1]);
	xp[i__] = diagin * (xp[i__] - spdiag * xp[i__ - 1]);
	yp[i__] = diagin * (yp[i__] - spdiag * yp[i__ - 1]);
/* 	   TEMP(N+1) = -DIAGIN*TEMP(NM1+I)*SPDIAG !was a bug, n+1 should be n+i */
	temp[*n + i__] = -diagin * temp[nm1 + i__] * spdiag;
	if (i__ == 1) {
	    temp[np1] = -diagin * spdig1;
	}
	spdiag = sinhin * (sinhs - dels);
	temp[i__] = diagin * spdiag;
	diag1 = diag2;
    }
L30:
    temp[nm1] = temp[*n + nm1] - temp[nm1];
    if (*n == 2) {
	goto L50;
    }
    i__1 = *n;
    for (i__ = 3; i__ <= i__1; ++i__) {
	ibak = np1 - i__;
	xp[ibak] -= temp[ibak] * xp[ibak + 1];
	yp[ibak] -= temp[ibak] * yp[ibak + 1];
	temp[ibak] = temp[*n + ibak] - temp[ibak] * temp[ibak + 1];
    }
L50:
    xp[*n] = (xp[*n] - spdig1 * xp[1] - spdiag * xp[nm1]) / (diag1 + diag2 + 
	    spdig1 * temp[1] + spdiag * temp[nm1]);
    yp[*n] = (yp[*n] - spdig1 * yp[1] - spdiag * yp[nm1]) / (diag1 + diag2 + 
	    spdig1 * temp[1] + spdiag * temp[nm1]);
    i__1 = nm1;
    for (i__ = 1; i__ <= i__1; ++i__) {
	xp[i__] += temp[i__] * xp[*n];
	yp[i__] += temp[i__] * yp[*n];
    }
/* --	done with setup, now do the interpolations */

    i__ = 2;
    sum = 0.f;
    delx = x[2] - x[1];
    dely = y[2] - y[1];
    dels = sqrt(delx * delx + dely * dely);
    i__1 = *nf;
    for (k = 1; k <= i__1; ++k) {
/* --	loop over the nf input values of T */
	tn = (d__1 = t[k] * s, abs(d__1));
	if (tn > s) {
	    tn = s;
	}
/* 	type *,'start, k,tn,sum,dels',k,tn,sum,dels */
L10:
	if (tn < sum) {
/* --	drop back */
/* 	type *,'drop back, I=',I */
	    --i__;
	    if (i__ >= 2) {
		delx = x[i__] - x[i__ - 1];
		dely = y[i__] - y[i__ - 1];
		dels = sqrt(delx * delx + dely * dely);
		sum -= dels;
		goto L10;
	    } else {
/* --	I was 1, we must have dels since we never start with 1 */
		i__ = *n;
		sum -= dels;
/* 	type *,'drop back, I was 1' */
		goto L10;
	    }
	}
	if (tn > sum + dels) {
/* --	go forward */
/* 	type *,'go forward, I=',I */
	    i1 = i__;
	    i__2 = *n;
	    for (i__ = i1; i__ <= i__2; ++i__) {
		delx = x[i__] - x[i__ - 1];
		dely = y[i__] - y[i__ - 1];
		dels = sqrt(delx * delx + dely * dely);
		if (sum + dels > tn) {
		    goto L40;
		}
		sum += dels;
	    }
/* --	must be between N and 1 */
	    i__ = 1;
	    dels = s - sum;
	}
L40:
	im1 = i__ - 1;
	if (im1 == 0) {
	    im1 = *n;
	}
/* 	type *,'resultants, i,im1,tn,sum,dels',i,im1,tn,sum,dels */
/* 	type *,'xp,yp',XP(I),YP(I) */
	del1 = tn - sum;
	del2 = dels - del1;
	exps1 = exp(sigmap * del1);
	sinhd1 = (exps1 - 1.f / exps1) * .5f;
	exps = exp(sigmap * del2);
	sinhd2 = (exps - 1.f / exps) * .5f;
	exps = exps1 * exps;
	sinhs = (exps - 1.f / exps) * .5f;
	xs[k] = (xp[i__] * sinhd1 + xp[im1] * sinhd2) / sinhs + ((x[i__] - xp[
		i__]) * del1 + (x[im1] - xp[im1]) * del2) / dels;
	ys[k] = (yp[i__] * sinhd1 + yp[im1] * sinhd2) / sinhs + ((y[i__] - yp[
		i__]) * del1 + (y[im1] - yp[im1]) * del2) / dels;
    }
    return ret_val;
} /* kurvp1_ */