/* scale.f -- translated by f2c (version of 3 February 1990  3:36:42).
   You must link the resulting object file with the libraries:
	-lF77 -lI77 -lm -lc   (in that order)
*/

#include "f2c.h"

/* Common Block Declarations */

struct {
    doublereal twopi, xlog2, xlog10, root2, rad, boltz, charge, ctok, gmin, 
	    reltol, abstol, vntol, trtol, chgtol, eps0, epssil, epsox, pivtol,
	     pivrel;
} knstnt_;

#define knstnt_1 knstnt_

/*<       subroutine scale(xmin,xmax,n,xminp,xmaxp,del) >*/
/* Subroutine */ int scale_(xmin, xmax, n, xminp, xmaxp, del)
doublereal *xmin, *xmax;
integer *n;
doublereal *xminp, *xmaxp, *del;
{
    /* Initialized data */

    static doublereal vint[5] = { 1.,2.,5.,10.,20. };
    static doublereal eps = 1e-12;

    /* System generated locals */
    doublereal d_1, d_2, d_3;

    /* Builtin functions */
    double d_lg10(), exp();

    /* Local variables */
    static doublereal a, b;
    static integer i;
    static doublereal xfact;
    static integer m1, m2, np, nx;
    static doublereal fm1, fm2;
    static integer nal;

/*<       implicit double precision (a-h,o-z) >*/

/*     this routine determines the 'optimal' scale to use for the plot of 
*/
/* some output variable. */


/*  adapted from algorithm 463 of 'collected algorithms of the cacm' */

/* spice version 2g.6  sccsid=knstnt 3/15/83 */
/*<       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok, >*/
/*<      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox, >*/
/*<      2   pivtol,pivrel >*/
/*<       integer xxor >*/
/*<       dimension vint(5) >*/
/*<       data vint / 1.0d0,2.0d0,5.0d0,10.0d0,20.0d0 / >*/
/*<       data eps / 1.0d-12 / >*/


/* ...  trap too-small data spread */
/* *********************************************************** */
/*  temporily check 'equality' this way */
/*<       if(xmin.eq.0.0d0.and.xmax.eq.0.0d0) go to 4 >*/
    if (*xmin == 0. && *xmax == 0.) {
	goto L4;
    }
/*<       if(dabs((xmax-xmin)/dmax1(dabs(xmin),dabs(xmax))).ge.1.0d-4) >*/
/*<      1  go to 10 >*/
/* Computing MAX */
    d_2 = abs(*xmin), d_3 = abs(*xmax);
    if ((d_1 = (*xmax - *xmin) / max(d_3,d_2), abs(d_1)) >= 1e-4) {
	goto L10;
    }
/*<     4 continue >*/
L4:
/*<       if (xmin.ge.0.0d0) go to 5 >*/
    if (*xmin >= 0.) {
	goto L5;
    }
/*<       xmax=0.5d0*xmin+eps >*/
    *xmax = *xmin * .5 + eps;
/*<       xmin=1.5d0*xmin-eps >*/
    *xmin = *xmin * 1.5 - eps;
/*<       go to 10 >*/
    goto L10;
/*<     5 xmax=1.5d0*xmin+eps >*/
L5:
    *xmax = *xmin * 1.5 + eps;
/*<       xmin=0.5d0*xmin-eps >*/
    *xmin = *xmin * .5 - eps;
/* ...  find approximate interval size, normalized to [1,10] */
/*<    10 a=(xmax-xmin)/dble(n) >*/
L10:
    a = (*xmax - *xmin) / (doublereal) (*n);
/*<       nal=idint(dlog10(a)) >*/
    nal = (integer) d_lg10(&a);
/*<       if (a.lt.1.0d0) nal=nal-1 >*/
    if (a < 1.) {
	--nal;
    }
/*<       xfact=dexp(xlog10*dble(nal)) >*/
    xfact = exp(knstnt_1.xlog10 * (doublereal) nal);
/*<       b=a/xfact >*/
    b = a / xfact;
/* ...  find closest permissible interval size */
/*<       do 20 i=1,3 >*/
    for (i = 1; i <= 3; ++i) {
/*<       if (b.lt.(vint(i)+eps)) go to 30 >*/
	if (b < vint[i - 1] + eps) {
	    goto L30;
	}
/*<    20 continue >*/
/* L20: */
    }
/*<       i=4 >*/
    i = 4;
/* ...  compute interval size */
/*<    30 del=vint(i)*xfact >*/
L30:
    *del = vint[i - 1] * xfact;
/*<       fm1=xmin/del >*/
    fm1 = *xmin / *del;
/*<       m1=fm1 >*/
    m1 = (integer) fm1;
/*<       if (fm1.lt.0.0d0) m1=m1-1 >*/
    if (fm1 < 0.) {
	--m1;
    }
/*<       if (dabs(dble(m1)+1.0d0-fm1).lt.eps) m1=m1+1 >*/
    if ((d_1 = (doublereal) m1 + 1. - fm1, abs(d_1)) < eps) {
	++m1;
    }
/* ...  compute new maximum and minimum limits */
/*<       xminp=del*dble(m1) >*/
    *xminp = *del * (doublereal) m1;
/*<       fm2=xmax/del >*/
    fm2 = *xmax / *del;
/*<       m2=fm2+1.0d0 >*/
    m2 = (integer) (fm2 + 1.);
/*<       if (fm2.lt.(-1.0d0)) m2=m2-1 >*/
    if (fm2 < -1.) {
	--m2;
    }
/*<       if (dabs(fm2+1.0d0-dble(m2)).lt.eps) m2=m2-1 >*/
    if ((d_1 = fm2 + 1. - (doublereal) m2, abs(d_1)) < eps) {
	--m2;
    }
/*<       xmaxp=del*dble(m2) >*/
    *xmaxp = *del * (doublereal) m2;
/*<       np=m2-m1 >*/
    np = m2 - m1;
/* ...  check whether another loop required */
/*<       if (np.le.n) go to 40 >*/
    if (np <= *n) {
	goto L40;
    }
/*<       i=i+1 >*/
    ++i;
/*<       go to 30 >*/
    goto L30;
/* ...  do final adjustments and correct for roundoff error(s) */
/*<    40 nx=(n-np)/2 >*/
L40:
    nx = (*n - np) / 2;
/*<       xminp=dmin1(xmin,xminp-dble(nx)*del) >*/
/* Computing MAX */
    d_1 = *xmin, d_2 = *xminp - (doublereal) nx * *del;
    *xminp = min(d_2,d_1);
/*<       xmaxp=dmax1(xmax,xminp+dble(n)*del) >*/
/* Computing MAX */
    d_1 = *xmax, d_2 = *xminp + (doublereal) (*n) * *del;
    *xmaxp = max(d_2,d_1);
/*<       return >*/
    return 0;
/*<       end >*/
} /* scale_ */

