This is lpc.stable.c in view mode; [Download] [Up]
/* lpc.stable.f -- translated by f2c (version of 26 January 1990 18:57:16).
You must link the resulting object file with the libraries:
-lF77 -lI77 -lm -lc (in that order)
*/
#include <f2c.h>
/* Table of constant values */
static integer c__1 = 1;
/* Subroutine */ int stabl_(anal, unstab, npoles, anal_len, unstab_len)
char *anal, *unstab;
integer *npoles;
ftnlen anal_len;
ftnlen unstab_len;
{
/* System generated locals */
integer i_1, i_2, i_3, i_4, i_5;
doublereal d_1, d_2;
doublecomplex z_1, z_2;
olist o_1;
/* Builtin functions */
integer f_open(), s_rdue(), do_uio(), e_rdue();
double sqrt(), atan2(), cos(), sin();
integer s_wdue(), e_wdue();
/* Subroutine */ int s_stop();
/* Local variables */
static integer nall;
static doublecomplex zero;
static real a[40];
static integer j, k;
static doublereal r[37];
static doublecomplex w[37];
static integer ndata;
static doublereal y[37];
static integer nlocs;
static doublereal rooti[36];
static integer l1, k4;
static doublereal rootr[36];
static integer ii, ll;
static doublereal th[37];
static doublecomplex ww;
extern /* Subroutine */ int factor_();
static integer kinsid;
static doublereal zz;
static integer kprint, k4m, jjj, kkk;
static doublecomplex one;
static doublereal eps;
/* Fortran I/O blocks */
static cilist io__9 = { 0, 19, 1, 0, 0 };
static cilist io__11 = { 0, 18, 1, 0, 0 };
static cilist io__29 = { 0, 18, 0, 0, 0 };
zero.r = 0., zero.i = 0.;
one.r = 1., one.i = 0.;
ndata = 4;
nlocs = *npoles + ndata << 2;
o_1.oerr = 0;
o_1.ounit = 18;
o_1.ofnmlen = 64;
o_1.ofnm = anal;
o_1.orl = nlocs;
o_1.osta = "old";
o_1.oacc = "direct";
o_1.ofm = "unformatted";
o_1.oblnk = 0;
f_open(&o_1);
o_1.oerr = 0;
o_1.ounit = 19;
o_1.ofnmlen = 64;
o_1.ofnm = unstab;
o_1.orl = 4;
o_1.osta = "old";
o_1.oacc = "direct";
o_1.ofm = "unformatted";
o_1.oblnk = 0;
f_open(&o_1);
k4 = *npoles + 1;
k4m = k4 - 1;
nall = k4m + ndata;
for (jjj = 1; jjj <= 10000; ++jjj) {
io__9.cirec = jjj;
i_1 = s_rdue(&io__9);
if (i_1 != 0) {
goto L999;
}
i_1 = do_uio(&c__1, (char *)&kkk, (ftnlen)sizeof(integer));
if (i_1 != 0) {
goto L999;
}
i_1 = e_rdue();
io__11.cirec = kkk;
i_1 = s_rdue(&io__11);
if (i_1 != 0) {
goto L999;
}
i_2 = nall;
for (ll = 1; ll <= i_2; ++ll) {
i_1 = do_uio(&c__1, (char *)&a[ll - 1], (ftnlen)sizeof(real));
if (i_1 != 0) {
goto L999;
}
}
i_1 = e_rdue();
i_1 = k4m;
for (ii = 1; ii <= i_1; ++ii) {
/* L1601: */
y[ii - 1] = -(doublereal)a[ii + 3];
}
y[k4 - 1] = (float)1.;
kprint = 0;
eps = 1.0000000000000008e-8;
/* L303: */
factor_(y, &k4, rootr, rooti, &kinsid, &kprint, &eps);
i_1 = k4m;
for (j = 1; j <= i_1; ++j) {
/* Computing 2nd power */
d_1 = rootr[j - 1];
/* Computing 2nd power */
d_2 = rooti[j - 1];
r[j - 1] = sqrt(d_1 * d_1 + d_2 * d_2);
th[j - 1] = atan2(rooti[j - 1], rootr[j - 1]);
if (r[j - 1] >= 1.) {
r[j - 1] = (float)1. / r[j - 1];
}
/* L100: */
}
i_1 = k4m;
for (k = 1; k <= i_1; ++k) {
/* L10: */
i_2 = k - 1;
w[i_2].r = zero.r, w[i_2].i = zero.i;
}
i_2 = k4 - 1;
w[i_2].r = one.r, w[i_2].i = one.i;
i_2 = k4m;
for (k = 1; k <= i_2; ++k) {
/* ww=dcmplx(rootr(k),rooti(k)) */
d_1 = r[k - 1] * cos(th[k - 1]);
d_2 = r[k - 1] * sin(th[k - 1]);
z_1.r = d_1, z_1.i = d_2;
ww.r = z_1.r, ww.i = z_1.i;
l1 = k4 - k;
i_1 = k4m;
for (j = l1; j <= i_1; ++j) {
/* L12: */
i_3 = j - 1;
i_4 = j;
i_5 = j - 1;
z_2.r = ww.r * w[i_5].r - ww.i * w[i_5].i, z_2.i = ww.r * w[
i_5].i + ww.i * w[i_5].r;
z_1.r = w[i_4].r - z_2.r, z_1.i = w[i_4].i - z_2.i;
w[i_3].r = z_1.r, w[i_3].i = z_1.i;
}
/* L20: */
i_3 = k4 - 1;
z_2.r = -ww.r, z_2.i = -ww.i;
i_4 = k4 - 1;
z_1.r = z_2.r * w[i_4].r - z_2.i * w[i_4].i, z_1.i = z_2.r * w[
i_4].i + z_2.i * w[i_4].r;
w[i_3].r = z_1.r, w[i_3].i = z_1.i;
}
i_3 = k4;
for (j = 2; j <= i_3; ++j) {
i_4 = j - 1;
zz = w[i_4].r;
a[k4 + ndata + 1 - j - 1] = -zz;
/* L30: */
}
io__29.cirec = kkk;
s_wdue(&io__29);
i_3 = nall;
for (ll = 1; ll <= i_3; ++ll) {
do_uio(&c__1, (char *)&a[ll - 1], (ftnlen)sizeof(real));
}
e_wdue();
/* L199: */
}
L999:
s_stop("", 0L);
} /* stabl_ */
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.