home *** CD-ROM | disk | FTP | other *** search
- /* 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_ */
-
-