i/series
 Home Manual Packages Global Index Keywords Quick Reference ``` /* 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. All rights reserved. */ 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. SEE ALSO: series_r, series_n */ { /* 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 */ list= where(mask); if (numberof(list)) { rrx= rr(list); s= merge((rrx^(1+nn(list))-1.0)/(rrx-1.0), s(where(!mask)), mask); } mask= (nn<0 & rnot); /* n<0 case */ list= where(mask); if (numberof(list)) { rrx= 1.0/rr(list); s= merge((rrx^(1-nn(list))-1.0)/(rrx-1.0), s(where(!mask)), mask); } 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. SEE ALSO: series_s, series_n */ { /* 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 snn); list= where(mask); 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; r= merge(sx^pow - pow + A*(nx/sx)^pow + B/sx, r(where(!mask)), mask); } mask= nbig & (ss<=nn); list= where(mask); if (numberof(list)) { sx= ss(list); nx= nn(list); sn= (sx-1.0)/(nx-1.0); n2r= 1.0/(nx*nx); r= merge(1.0 - 1.0/sx + n2r*sn*sn*(nx+1.0 - sn), r(where(!mask)), mask); } /* Polish the approximation using Newton-Raphson iterations. There are never many of these, so do the entire vector together. */ mask= nbig & (ss!=nn); list= where(mask); if (numberof(list)) { rx= r(list); ss= ss(list); nn= nn(list); for (;;) { rr= rx-1.0; rn= rx^(nn-1); rrss= rr*ss; delta= rrss - (rx*rn-1.0); if (allof(abs(delta)<=1.e-9*abs(rrss))) break; rx+= delta/(nn*rn-ss); } /* try to get it to machine precision */ if (anyof(delta)) rx+= delta/(nn*rn-ss); r= merge(rx, r(where(!mask)), mask); } 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. SEE ALSO: series_s, series_r */ { n= 1.0 + (r-1.0)*s; bad= (n<=0.0 | s<1.0); list= where(bad); if (numberof(list)) n= merge(array(1.0, numberof(list)), n(where(!bad)), bad); mask= r==1.0 & !bad; list= where(mask); if (numberof(list)) n= merge((s+0.0*n)(list), n(where(!mask)), mask); mask= !mask & !bad; list= where(mask); if (numberof(list)) n= merge(log(n(list))/log(r(list)), n(where(!mask)), mask); return long(ceil(n))-1; } ```