00001 #include <math.h>
00002 #include <stdio.h>
00003 #include <stdlib.h>
00004 #include <string.h>
00005 #include <assert.h>
00006 #include "f2c.h"
00007
00008 #ifdef _MSC_VER
00009 #pragma warning (disable: 4244)
00010 #endif
00011
00012
00013 extern void
00014 s_wsfe(cilist * f)
00015 {;
00016 }
00017 extern void
00018 e_wsfe(void)
00019 {;
00020 }
00021 extern void
00022 do_fio(integer * c, char *s, ftnlen l)
00023 {;
00024 }
00025
00026
00027
00028
00029 extern int
00030 s_rnge(char *var, int index, char *routine, int lineno)
00031 {
00032 fprintf(stderr,
00033 "array index out-of-bounds for %s[%d] in routine %s:%d\n", var,
00034 index, routine, lineno);
00035 fflush(stderr);
00036 assert(2+2 == 5);
00037 return 0;
00038 }
00039
00040
00041 #ifdef KR_headers
00042 extern double sqrt();
00043 float
00044 f__cabs(real, imag)
00045 float real, imag;
00046 #else
00047 #undef abs
00048
00049 float
00050 f__cabs(float real, float imag)
00051 #endif
00052 {
00053 float temp;
00054
00055 if (real < 0)
00056 real = -real;
00057 if (imag < 0)
00058 imag = -imag;
00059 if (imag > real) {
00060 temp = real;
00061 real = imag;
00062 imag = temp;
00063 }
00064 if ((imag + real) == real)
00065 return ((float) real);
00066
00067 temp = imag / real;
00068 temp = real * sqrt(1.0 + temp * temp);
00069 return (temp);
00070 }
00071
00072
00073 VOID
00074 #ifdef KR_headers
00075 s_cnjg(r, z)
00076 complex *r, *z;
00077 #else
00078 s_cnjg(complex * r, complex * z)
00079 #endif
00080 {
00081 r->r = z->r;
00082 r->i = -z->i;
00083 }
00084
00085
00086 #ifdef KR_headers
00087 float
00088 r_imag(z)
00089 complex *z;
00090 #else
00091 float
00092 r_imag(complex * z)
00093 #endif
00094 {
00095 return (z->i);
00096 }
00097
00098
00099 #define log10e 0.43429448190325182765
00100
00101 #ifdef KR_headers
00102 double log();
00103 float
00104 r_lg10(x)
00105 real *x;
00106 #else
00107 #undef abs
00108
00109 float
00110 r_lg10(real * x)
00111 #endif
00112 {
00113 return (log10e * log(*x));
00114 }
00115
00116
00117 #ifdef KR_headers
00118 float
00119 r_sign(a, b)
00120 real *a, *b;
00121 #else
00122 float
00123 r_sign(real * a, real * b)
00124 #endif
00125 {
00126 float x;
00127 x = (*a >= 0 ? *a : -*a);
00128 return (*b >= 0 ? x : -x);
00129 }
00130
00131
00132 #ifdef KR_headers
00133 double floor();
00134 integer
00135 i_dnnt(x)
00136 real *x;
00137 #else
00138 #undef abs
00139
00140 integer
00141 i_dnnt(real * x)
00142 #endif
00143 {
00144 return ((*x) >= 0 ? floor(*x + .5) : -floor(.5 - *x));
00145 }
00146
00147
00148 #ifdef KR_headers
00149 double pow();
00150 double
00151 pow_dd(ap, bp)
00152 doublereal *ap, *bp;
00153 #else
00154 #undef abs
00155
00156 double
00157 pow_dd(doublereal * ap, doublereal * bp)
00158 #endif
00159 {
00160 return (pow(*ap, *bp));
00161 }
00162
00163
00164 #ifdef KR_headers
00165 float
00166 pow_ri(ap, bp)
00167 real *ap;
00168 integer *bp;
00169 #else
00170 float
00171 pow_ri(real * ap, integer * bp)
00172 #endif
00173 {
00174 float pow, x;
00175 integer n;
00176 unsigned long u;
00177
00178 pow = 1;
00179 x = *ap;
00180 n = *bp;
00181
00182 if (n != 0) {
00183 if (n < 0) {
00184 n = -n;
00185 x = 1 / x;
00186 }
00187 for (u = n;;) {
00188 if (u & 01)
00189 pow *= x;
00190 if (u >>= 1)
00191 x *= x;
00192 else
00193 break;
00194 }
00195 }
00196 return (pow);
00197 }
00198
00199
00200
00201
00202
00203 #define NO_OVERWRITE
00204
00205
00206 #ifndef NO_OVERWRITE
00207
00208 #undef abs
00209 #ifdef KR_headers
00210 extern char *F77_aloc();
00211 extern void free();
00212 extern void exit_();
00213 #else
00214
00215 extern char *F77_aloc(ftnlen, char *);
00216 #endif
00217
00218 #endif
00219
00220 VOID
00221 #ifdef KR_headers
00222 s_cat(lp, rpp, rnp, np, ll)
00223 char *lp, *rpp[];
00224 ftnlen rnp[], *np, ll;
00225 #else
00226 s_cat(char *lp, char *rpp[], ftnlen rnp[], ftnlen * np, ftnlen ll)
00227 #endif
00228 {
00229 ftnlen i, nc;
00230 char *rp;
00231 ftnlen n = *np;
00232 #ifndef NO_OVERWRITE
00233 ftnlen L, m;
00234 char *lp0, *lp1;
00235
00236 lp0 = 0;
00237 lp1 = lp;
00238 L = ll;
00239 i = 0;
00240 while (i < n) {
00241 rp = rpp[i];
00242 m = rnp[i++];
00243 if (rp >= lp1 || rp + m <= lp) {
00244 if ((L -= m) <= 0) {
00245 n = i;
00246 break;
00247 }
00248 lp1 += m;
00249 continue;
00250 }
00251 lp0 = lp;
00252 lp = lp1 = F77_aloc(L = ll, "s_cat");
00253 break;
00254 }
00255 lp1 = lp;
00256 #endif
00257 for (i = 0; i < n; ++i) {
00258 nc = ll;
00259 if (rnp[i] < nc)
00260 nc = rnp[i];
00261 ll -= nc;
00262 rp = rpp[i];
00263 while (--nc >= 0)
00264 *lp++ = *rp++;
00265 }
00266 while (--ll >= 0)
00267 *lp++ = ' ';
00268 #ifndef NO_OVERWRITE
00269 if (lp0) {
00270 memmove(lp0, lp1, L);
00271 free(lp1);
00272 }
00273 #endif
00274 }
00275
00276
00277
00278
00279 #ifdef KR_headers
00280 integer
00281 s_cmp(a0, b0, la, lb)
00282 char *a0, *b0;
00283 ftnlen la, lb;
00284 #else
00285 integer
00286 s_cmp(char *a0, char *b0, ftnlen la, ftnlen lb)
00287 #endif
00288 {
00289 register unsigned char *a, *aend, *b, *bend;
00290 a = (unsigned char *) a0;
00291 b = (unsigned char *) b0;
00292 aend = a + la;
00293 bend = b + lb;
00294
00295 if (la <= lb) {
00296 while (a < aend)
00297 if (*a != *b)
00298 return (*a - *b);
00299 else {
00300 ++a;
00301 ++b;
00302 }
00303
00304 while (b < bend)
00305 if (*b != ' ')
00306 return (' ' - *b);
00307 else
00308 ++b;
00309 }
00310
00311 else {
00312 while (b < bend)
00313 if (*a == *b) {
00314 ++a;
00315 ++b;
00316 }
00317 else
00318 return (*a - *b);
00319 while (a < aend)
00320 if (*a != ' ')
00321 return (*a - ' ');
00322 else
00323 ++a;
00324 }
00325 return (0);
00326 }
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338 #ifdef KR_headers
00339 VOID
00340 s_copy(a, b, la, lb)
00341 register char *a, *b;
00342 ftnlen la, lb;
00343 #else
00344 void
00345 s_copy(register char *a, register char *b, ftnlen la, ftnlen lb)
00346 #endif
00347 {
00348 register char *aend, *bend;
00349
00350 aend = a + la;
00351
00352 if (la <= lb)
00353 #ifndef NO_OVERWRITE
00354 if (a <= b || a >= b + la)
00355 #endif
00356 while (a < aend)
00357 *a++ = *b++;
00358 #ifndef NO_OVERWRITE
00359 else
00360 for (b += la; a < aend;)
00361 *--aend = *--b;
00362 #endif
00363
00364 else {
00365 bend = b + lb;
00366 #ifndef NO_OVERWRITE
00367 if (a <= b || a >= bend)
00368 #endif
00369 while (b < bend)
00370 *a++ = *b++;
00371 #ifndef NO_OVERWRITE
00372 else {
00373 a += lb;
00374 while (b < bend)
00375 *--a = *--bend;
00376 a += lb;
00377 }
00378 #endif
00379 while (a < aend)
00380 *a++ = ' ';
00381 }
00382 }
00383
00384
00385 #ifdef KR_headers
00386 float f__cabs();
00387 float
00388 z_abs(z)
00389 complex *z;
00390 #else
00391 float f__cabs(float, float);
00392 float
00393 z_abs(complex * z)
00394 #endif
00395 {
00396 return (f__cabs(z->r, z->i));
00397 }
00398
00399
00400 #ifdef KR_headers
00401 extern void sig_die();
00402 VOID
00403 z_div(c, a, b)
00404 complex *a, *b, *c;
00405 #else
00406 extern void sig_die(char *, int);
00407 void
00408 z_div(complex * c, complex * a, complex * b)
00409 #endif
00410 {
00411 float ratio, den;
00412 float abr, abi;
00413
00414 if ((abr = b->r) < 0.)
00415 abr = -abr;
00416 if ((abi = b->i) < 0.)
00417 abi = -abi;
00418 if (abr <= abi) {
00419
00420
00421
00422 ratio = b->r / b->i;
00423 den = b->i * (1 + ratio * ratio);
00424 c->r = (a->r * ratio + a->i) / den;
00425 c->i = (a->i * ratio - a->r) / den;
00426 }
00427
00428 else {
00429 ratio = b->i / b->r;
00430 den = b->r * (1 + ratio * ratio);
00431 c->r = (a->r + a->i * ratio) / den;
00432 c->i = (a->i - a->r * ratio) / den;
00433 }
00434
00435 }
00436
00437
00438 #ifdef KR_headers
00439 double sqrt();
00440 double f__cabs();
00441 VOID
00442 z_sqrt(r, z)
00443 complex *r, *z;
00444 #else
00445 #undef abs
00446
00447 extern float f__cabs(float, float);
00448 void
00449 z_sqrt(complex * r, complex * z)
00450 #endif
00451 {
00452 float mag;
00453
00454 if ((mag = f__cabs(z->r, z->i)) == 0.)
00455 r->r = r->i = 0.;
00456 else if (z->r > 0) {
00457 r->r = sqrt(0.5 * (mag + z->r));
00458 r->i = z->i / r->r / 2;
00459 }
00460 else {
00461 r->i = sqrt(0.5 * (mag - z->r));
00462 if (z->i < 0)
00463 r->i = -r->i;
00464 r->r = z->i / r->i / 2;
00465 }
00466 }
00467
00468 #ifdef __cplusplus
00469 extern "C" {
00470 #endif
00471
00472 #ifdef KR_headers
00473 integer pow_ii(ap, bp) integer *ap, *bp;
00474 #else
00475 integer pow_ii(integer * ap, integer * bp)
00476 #endif
00477 {
00478 integer pow, x, n;
00479 unsigned long u;
00480
00481 x = *ap;
00482 n = *bp;
00483
00484 if (n <= 0) {
00485 if (n == 0 || x == 1)
00486 return 1;
00487 if (x != -1)
00488 return x != 0 ? 1 / x : 0;
00489 n = -n;
00490 } u = n;
00491 for (pow = 1;;) {
00492 if (u & 01)
00493 pow *= x;
00494 if (u >>= 1)
00495 x *= x;
00496 else
00497 break;
00498 }
00499 return (pow);
00500 }
00501 #ifdef __cplusplus
00502 }
00503 #endif
00504
00505 #ifdef KR_headers
00506 extern void f_exit();
00507 VOID
00508 s_stop(s, n)
00509 char *s;
00510 ftnlen n;
00511 #else
00512 #undef abs
00513 #undef min
00514 #undef max
00515 #ifdef __cplusplus
00516 extern "C" {
00517 #endif
00518 #ifdef __cplusplus
00519 extern "C" {
00520 #endif
00521 void f_exit(void);
00522
00523 int s_stop(char *s, ftnlen n)
00524 #endif
00525 {
00526 int i;
00527
00528 if (n > 0) {
00529 fprintf(stderr, "STOP ");
00530 for (i = 0; i < n; ++i)
00531 putc(*s++, stderr);
00532 fprintf(stderr, " statement executed\n");
00533 }
00534 #ifdef NO_ONEXIT
00535 f_exit();
00536 #endif
00537 exit(0);
00538
00539
00540
00541
00542
00543 return 0;
00544 }
00545 #ifdef __cplusplus
00546 }
00547 #endif
00548 #ifdef __cplusplus
00549 }
00550 #endif