```/*
SERIES.I
Routines for handling geometric series (e.g.- tapered zoning).

\$Id: series.i,v 1.1 1993/08/27 18:50:06 munro Exp \$
*/
/*    Copyright (c) 1994.  The Regents of the University of California.

func series_s (r, n)
/* DOCUMENT series_s(r, n)
returns the sum s of the finite geometric series
1 + r + r^2 + r^3 + ... + r^n
Using n<0 is equivalent to using the reciprocal of r, that is,
series_s(r, -n) == series_s(1./r, n)
R or N or both may be arrays, as long as they are conformable.
*/
{
/* force conformability immediately */
dims= dimsof(r, n);
rr= r + array(0.0, dims);
nn= n + array(0, dims);

/* form result array, initialized to n==0 result */
s= array(1.0, dims);

/* subdivide into n==0, n>0, n<0, and r==1.0 cases */
rnot= rr!=1.0;

mask= (nn>0 & rnot);   /* n>0 case */
if (numberof(list)) {
rrx= rr(list);
}

mask= (nn<0 & rnot);   /* n<0 case */
if (numberof(list)) {
rrx= 1.0/rr(list);
}

list= where(!rnot);   /* r==1.0 case */
if (numberof(list)) s= merge(s(where(rnot)), abs(nn(list)), rnot);

return s;
}

func series_r (s, n)
/* DOCUMENT series_r(s, n)
returns the ratio r of the finite geometric series, given the sum s:
1 + r + r^2 + r^3 + ... + r^n = s
Using n<0 will return the the reciprocal of n>0 result, that is,
series_r(s, -n) == 1.0/series_r(s, n)
If n==0, returns s-1 (the n==1 result).
S or N or both may be arrays, as long as they are conformable.
*/
{
/* force conformability immediately */
dims= dimsof(s, n);
ss= s + array(0.0, dims);
nn= n + array(0, dims);

/* form result array, initialized to abs(n)<2 result */
r= ss-1.0;

/* work only with n>0 -- take reciprocals at end if necessary */
nneg= nn<0;
nn= abs(nn)+1;
nbig= nn>2;

/* compute an approximate result which has exact values and
derivatives at s==1, s==n, and s->infinity --
different approximations apply for s>n and s<n */
if (numberof(list)) {
sx= ss(list);
nx= nn(list);
pow= 1.0/(nx-1.0);
npow= nx^pow - 1.0;
n2r= 1.0/(nx-2.0);
A= (2.0-nx*npow)*n2r;
B= (2.0*npow-nx*pow)*nx*n2r;
}
if (numberof(list)) {
sx= ss(list);
nx= nn(list);
sn= (sx-1.0)/(nx-1.0);
n2r= 1.0/(nx*nx);
}

/* Polish the approximation using Newton-Raphson iterations.
There are never many of these, so do the entire vector together.  */
if (numberof(list)) {
rx= r(list);
ss= ss(list);
nn= nn(list);
for (;;) {
rr= rx-1.0;
rn= rx^(nn-1);
rx+= delta/(nn*rn-ss);
}
/* try to get it to machine precision */
if (anyof(delta)) rx+= delta/(nn*rn-ss);
}

list= where(nneg);
if (numberof(list)) r= merge(1.0/r(list), r(where(!nneg)), nneg);

return r;
}

func series_n (r, s)
/* DOCUMENT series_n(r, s)
returns the minimum number n of terms required for the geometric
series
1 + r + r^2 + r^3 + ... + r^n = s
to reach at least the given value s.  An alternate viewpoint is
that n is the minimum number of terms required to achieve the
sum s, with a ratio no larger than r.
Returns 0 if r<1 and s>1/(1-r), or if s<1.
The routine makes the most sense for r>1 and s substantially
greater than 1.  The intended use is to determine the minimum
number of zones required to span a given thickness t with a given
minimum zone size z, and maximum taper ratio r (assumed >1 here):
n= series_n(r, t/z);
With this n, you have the option of adjusting r or z downwards
(using series_r or series_s, respectively) to achieve the final
desired zoning.
R or S or both may be arrays, as long as they are conformable.
*/
{
n= 1.0 + (r-1.0)*s;
if (numberof(list))