00001
00004
00005
00006
00007
00008
00009
00010
00011 #include "prototypes.h"
00012 #include "prototypes2.h"
00013
00014
00015
00016 static integer c__4 = 4;
00017 static integer c__15 = 15;
00018 static integer c__16 = 16;
00019 static integer c__5 = 5;
00020 static integer c__14 = 14;
00021 static integer c__9 = 9;
00022 static integer c__1 = 1;
00023 static integer c__2 = 2;
00024 static integer c__25 = 25;
00025 static doublereal c_b147 = .5;
00026 static doublereal c_b148 = 0.;
00027 static integer c__0 = 0;
00028
00029 int zuchk_(doublereal *yr, doublereal *yi, integer *nz, doublereal *ascle, doublereal *tol);
00030 int zunhj_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *ipmtr, doublereal *tol, doublereal *phir, doublereal *phii, doublereal *argr, doublereal *argi, doublereal *zeta1r, doublereal *zeta1i, doublereal *zeta2r, doublereal *zeta2i, doublereal *asumr, doublereal *asumi, doublereal *bsumr, doublereal *bsumi);
00031 int zuoik_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *ikflg, integer *n, doublereal *yr, doublereal *yi, integer *nuf, doublereal *tol, doublereal *elim, doublereal *alim);
00032
00033 int zacon_(double*, double*, double*, integer*,
00034 integer*, integer*, double*, double*,
00035 integer*, double*, double*, double*,
00036 double*, double*);
00037 int zbknu_(double*, double*, double*, integer*,
00038 integer*, double*, double*,
00039 integer*, double*, double*, double*);
00040 int zbunk_(double*, double*, double*, integer*,
00041 integer*, integer*, double*, double*,
00042 integer*, double*, double*, double*);
00043 int zuoik_(double*, double*, double*, integer*,
00044 integer*, integer*, double*, double*,
00045 integer*, double*, double*, double*);
00046
00047 doublereal d1mach_(const integer*);
00048 integer i1mach_(const integer*);
00049 doublereal myzabs_(const doublereal*, const doublereal*);
00050
00051 static doublereal ascle, csgni, csgnr, cspni, cspnr;
00052 extern int zs1s2_(doublereal *zrr, doublereal *zri, doublereal *s1r, doublereal *s1i, doublereal *s2r, doublereal *s2i, integer *nz, doublereal *ascle, doublereal *alim, integer *iuf);
00053 extern int zbknu_(double*, double*, double*, integer*,
00054 integer*, double*, double*,
00055 integer*, double*, double*, double*),
00056 zseri_(doublereal *zr, doublereal *zi, doublereal *fnu,
00057 integer *kode, integer *n, doublereal *yr, doublereal *yi,
00058 integer *nz,
00059 doublereal *tol, doublereal *elim, doublereal *alim);
00060 extern int
00061 zmlri_(doublereal *zr, doublereal *zi, doublereal *fnu,
00062 integer *kode, integer *n,
00063 doublereal *yr, doublereal *yi, integer *nz, doublereal *tol),
00064 zasyi_(doublereal *zr, doublereal *zi, doublereal *fnu,
00065 integer *kode, integer *n,
00066 doublereal *yr, doublereal *yi, integer *nz,
00067 doublereal *rl, doublereal *tol,
00068 doublereal *elim, doublereal *alim);
00069
00070 extern int zbinu_(doublereal *zr, doublereal *zi,
00071 doublereal *fnu, integer *kode,
00072 integer *n, doublereal *cyr,
00073 doublereal *cyi, integer *nz,
00074 doublereal *rl, doublereal *fnul,
00075 doublereal *tol, doublereal *elim,
00076 doublereal *alim),
00077 zbknu_(double*, double*, double*, integer*,
00078 integer*, double*, double*,
00079 integer*, double*, double*, double*);
00080
00081 ;
00082 extern int zshch_(doublereal *zr, doublereal *zi, doublereal *cshr, doublereal *cshi, doublereal *cchr, doublereal *cchi);
00083 extern int zkscl_(doublereal *zrr, doublereal *zri, doublereal *fnu, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *rzr, doublereal *rzi, doublereal *ascle, doublereal *tol, doublereal *elim);
00084 extern int zsqrt_(const doublereal *ar, const doublereal *ai, doublereal *br, doublereal *bi);
00085
00086 extern int zlog_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, integer *ierr);
00087 extern int zdiv_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, doublereal *cr, doublereal *ci);
00088 extern int zexp_(const doublereal *ar, const doublereal *ai, doublereal *br, doublereal *bi),
00089 zmlt_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, doublereal *cr, doublereal *ci);
00090 extern doublereal dgamln_(doublereal *z__, integer *ierr);
00091
00092 extern int zunhj_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *ipmtr, doublereal *tol, doublereal *phir, doublereal *phii, doublereal *argr, doublereal *argi, doublereal *zeta1r, doublereal *zeta1i, doublereal *zeta2r, doublereal *zeta2i, doublereal *asumr, doublereal *asumi, doublereal *bsumr, doublereal *bsumi);
00093
00094 extern int zbuni_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *n, doublereal *yr, doublereal *yi, integer *nz, integer *nui, integer *nlast, doublereal *fnul, doublereal *tol, doublereal *elim, doublereal *alim),
00095 zseri_(doublereal *zr, doublereal *zi, doublereal *fnu,
00096 integer *kode, integer *n, doublereal *yr, doublereal *yi,
00097 integer *nz,
00098 doublereal *tol, doublereal *elim, doublereal *alim),
00099 zmlri_(doublereal *zr, doublereal *zi, doublereal *fnu,
00100 integer *kode, integer *n,
00101 doublereal *yr, doublereal *yi, integer *nz, doublereal *tol),
00102 zasyi_(doublereal *zr, doublereal *zi, doublereal *fnu,
00103 integer *kode, integer *n,
00104 doublereal *yr, doublereal *yi, integer *nz,
00105 doublereal *rl, doublereal *tol,
00106 doublereal *elim, doublereal *alim),
00107 zuoik_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *ikflg, integer *n, doublereal *yr, doublereal *yi, integer *nuf, doublereal *tol, doublereal *elim, doublereal *alim),
00108 zwrsk_(doublereal *zrr, doublereal *zri, doublereal *fnu, integer *kode, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *cwr, doublereal *cwi, doublereal *tol, doublereal *elim, doublereal *alim);
00109
00110 extern int zuni1_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *n, doublereal *yr, doublereal *yi, integer *nz, integer *nlast, doublereal *fnul, doublereal *tol, doublereal *elim, doublereal *alim), zuni2_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *n, doublereal *yr, doublereal *yi, integer *nz, integer *nlast, doublereal *fnul, doublereal *tol, doublereal *elim, doublereal *alim);
00111
00112 extern int xerror_(char* mess, integer* nmess,
00113 integer* l1, integer* l2, ftnlen mess_len);
00114 extern int zunk1_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *mr, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *tol, doublereal *elim, doublereal *alim), zunk2_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *mr, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *tol, doublereal *elim, doublereal *alim);
00115
00116
00117 extern int zlog_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, integer *ierr),
00118 zdiv_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, doublereal *cr, doublereal *ci),
00119 zmlt_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, doublereal *cr, doublereal *ci);
00120
00121 extern int zunik_(doublereal *zrr, doublereal *zri, doublereal *fnu, integer *ikflg, integer *ipmtr, doublereal *tol, integer *init, doublereal *phir, doublereal *phii, doublereal *zeta1r, doublereal *zeta1i, doublereal *zeta2r, doublereal *zeta2i, doublereal *sumr, doublereal *sumi, doublereal *cwrkr, doublereal *cwrki), zuoik_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *ikflg, integer *n, doublereal *yr, doublereal *yi, integer *nuf, doublereal *tol, doublereal *elim, doublereal *alim);
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132 int zbesh_ (doublereal* zr, doublereal* zi, doublereal* fnu,
00133 integer* kode, integer* m, integer* n,
00134 doublereal* cyr, doublereal* cyi,
00135 integer* nz, integer* ierr)
00136 {
00137
00138
00139 static doublereal hpi = 1.57079632679489662;
00140
00141
00142 integer i__1, i__2;
00143 doublereal d__1, d__2;
00144
00145
00146
00147
00148 static doublereal alim, elim, atol, rhpi;
00149 static integer inuh;
00150 static doublereal fnul, rtol;
00151 static integer i__, k;
00152 static doublereal ascle, csgni;
00153 static doublereal csgnr;
00154 static integer k1;
00155 static integer k2;
00156 static doublereal aa, bb, fn;
00157 static integer mm;
00158 static doublereal az;
00159 static integer ir, nn;
00160 static doublereal rl;
00161 static integer mr, nw;
00162 static doublereal dig, arg, aln, fmm, r1m5, ufl, sgn;
00163 static integer nuf, inu;
00164 static doublereal tol, sti, zni, zti, str, znr;
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325 --cyi;
00326 --cyr;
00327
00328
00329
00330
00331 *ierr = 0;
00332 *nz = 0;
00333 if (*zr == 0. && *zi == 0.) {
00334 *ierr = 1;
00335 }
00336 if (*fnu < 0.) {
00337 *ierr = 1;
00338 }
00339 if (*m < 1 || *m > 2) {
00340 *ierr = 1;
00341 }
00342 if (*kode < 1 || *kode > 2) {
00343 *ierr = 1;
00344 }
00345 if (*n < 1) {
00346 *ierr = 1;
00347 }
00348 if (*ierr != 0) {
00349 return 0;
00350 }
00351 nn = *n;
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364 d__1 = d1mach_(&c__4);
00365 tol = max(d__1,1e-18);
00366 k1 = i1mach_(&c__15);
00367 k2 = i1mach_(&c__16);
00368 r1m5 = d1mach_(&c__5);
00369
00370 i__1 = abs(k1), i__2 = abs(k2);
00371 k = min(i__1,i__2);
00372 elim = ((doublereal) ((real) k) * r1m5 - 3.) * 2.303;
00373 k1 = i1mach_(&c__14) - 1;
00374 aa = r1m5 * (doublereal) ((real) k1);
00375 dig = min(aa,18.);
00376 aa *= 2.303;
00377
00378 d__1 = -aa;
00379 alim = elim + max(d__1,-41.45);
00380 fnul = (dig - 3.) * 6. + 10.;
00381 rl = dig * 1.2 + 3.;
00382 fn = *fnu + (doublereal) ((real) (nn - 1));
00383 mm = 3 - *m - *m;
00384 fmm = (doublereal) ((real) mm);
00385 znr = fmm * *zi;
00386 zni = -fmm * *zr;
00387
00388
00389
00390 az = myzabs_(zr, zi);
00391 aa = .5 / tol;
00392 bb = (doublereal) ((real) i1mach_(&c__9)) * .5;
00393 aa = min(aa,bb);
00394 if (az > aa) {
00395 goto L260;
00396 }
00397 if (fn > aa) {
00398 goto L260;
00399 }
00400 aa = sqrt(aa);
00401 if (az > aa) {
00402 *ierr = 3;
00403 }
00404 if (fn > aa) {
00405 *ierr = 3;
00406 }
00407
00408
00409
00410 ufl = d1mach_(&c__1) * 1e3;
00411 if (az < ufl) {
00412 goto L230;
00413 }
00414 if (*fnu > fnul) {
00415 goto L90;
00416 }
00417 if (fn <= 1.) {
00418 goto L70;
00419 }
00420 if (fn > 2.) {
00421 goto L60;
00422 }
00423 if (az > tol) {
00424 goto L70;
00425 }
00426 arg = az * .5;
00427 aln = -fn * log(arg);
00428 if (aln > elim) {
00429 goto L230;
00430 }
00431 goto L70;
00432 L60:
00433 zuoik_(&znr, &zni, fnu, kode, &c__2, &nn, &cyr[1], &cyi[1], &nuf, &tol, &
00434 elim, &alim);
00435 if (nuf < 0) {
00436 goto L230;
00437 }
00438 *nz += nuf;
00439 nn -= nuf;
00440
00441
00442
00443
00444 if (nn == 0) {
00445 goto L140;
00446 }
00447 L70:
00448 if (znr < 0. || (znr == 0. && zni < 0. && *m == 2)) {
00449 goto L80;
00450 }
00451
00452
00453
00454
00455 zbknu_(&znr, &zni, fnu, kode, &nn, &cyr[1], &cyi[1], nz, &tol, &elim, &
00456 alim);
00457 goto L110;
00458
00459
00460
00461 L80:
00462 mr = -mm;
00463 zacon_(&znr, &zni, fnu, kode, &mr, &nn, &cyr[1], &cyi[1], &nw, &rl, &fnul,
00464 &tol, &elim, &alim);
00465 if (nw < 0) {
00466 goto L240;
00467 }
00468 *nz = nw;
00469 goto L110;
00470 L90:
00471
00472
00473
00474 mr = 0;
00475 if (znr >= 0. && (znr != 0. || zni >= 0. || *m != 2)) {
00476 goto L100;
00477 }
00478 mr = -mm;
00479 if (znr != 0. || zni >= 0.) {
00480 goto L100;
00481 }
00482 znr = -znr;
00483 zni = -zni;
00484 L100:
00485 zbunk_(&znr, &zni, fnu, kode, &mr, &nn, &cyr[1], &cyi[1], &nw, &tol, &
00486 elim, &alim);
00487 if (nw < 0) {
00488 goto L240;
00489 }
00490 *nz += nw;
00491 L110:
00492
00493
00494
00495
00496
00497 d__1 = -fmm;
00498 sgn = d_sign(&hpi, &d__1);
00499
00500
00501
00502
00503 inu = (integer) ((real) (*fnu));
00504 inuh = inu / 2;
00505 ir = inu - (inuh << 1);
00506 arg = (*fnu - (doublereal) ((real) (inu - ir))) * sgn;
00507 rhpi = 1. / sgn;
00508
00509
00510 csgni = rhpi * cos(arg);
00511 csgnr = -rhpi * sin(arg);
00512 if (inuh % 2 == 0) {
00513 goto L120;
00514 }
00515
00516
00517 csgnr = -csgnr;
00518 csgni = -csgni;
00519 L120:
00520 zti = -fmm;
00521 rtol = 1. / tol;
00522 ascle = ufl * rtol;
00523 i__1 = nn;
00524 for (i__ = 1; i__ <= i__1; ++i__) {
00525
00526
00527
00528
00529
00530
00531 aa = cyr[i__];
00532 bb = cyi[i__];
00533 atol = 1.;
00534
00535 d__1 = abs(aa), d__2 = abs(bb);
00536 if (max(d__1,d__2) > ascle) {
00537 goto L135;
00538 }
00539 aa *= rtol;
00540 bb *= rtol;
00541 atol = tol;
00542 L135:
00543 str = aa * csgnr - bb * csgni;
00544 sti = aa * csgni + bb * csgnr;
00545 cyr[i__] = str * atol;
00546 cyi[i__] = sti * atol;
00547 str = -csgni * zti;
00548 csgni = csgnr * zti;
00549 csgnr = str;
00550
00551 }
00552 return 0;
00553 L140:
00554 if (znr < 0.) {
00555 goto L230;
00556 }
00557 return 0;
00558 L230:
00559 *nz = 0;
00560 *ierr = 2;
00561 return 0;
00562 L240:
00563 if (nw == -1) {
00564 goto L230;
00565 }
00566 *nz = 0;
00567 *ierr = 5;
00568 return 0;
00569 L260:
00570 *nz = 0;
00571 *ierr = 4;
00572 return 0;
00573 }
00574
00575
00576 doublereal d1mach_(const integer* i__)
00577 {
00578
00579
00580 static struct {
00581 integer e_1[10];
00582 doublereal fill_2[1];
00583 doublereal e_3;
00584 } equiv_4 = { {2002288515, 1050897, 1487780761, 2146426097,
00585 -1209488034, 1017118298, -1209488034, 1018166874, 1352628735,
00586 1070810131}, {0}, 0. };
00587
00588
00589
00590 doublereal ret_val;
00591
00592
00593 #define log10 ((integer *)&equiv_4 + 8)
00594 #define dmach ((doublereal *)&equiv_4)
00595 #define large ((integer *)&equiv_4 + 2)
00596 #define small ((integer *)&equiv_4)
00597 #define diver ((integer *)&equiv_4 + 6)
00598 #define right ((integer *)&equiv_4 + 4)
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001 if (*i__ < 1 || *i__ > 5) {
01002 xerror_("D1MACH -- I OUT OF BOUNDS", &c__25, &c__1, &c__2, (ftnlen)25)
01003 ;
01004 }
01005
01006 ret_val = dmach[*i__ - 1];
01007 return ret_val;
01008
01009 }
01010
01011 #undef right
01012 #undef diver
01013 #undef small
01014 #undef large
01015 #undef dmach
01016 #undef log10
01017
01018
01019 doublereal dgamln_(doublereal* z__, integer* ierr)
01020 {
01021
01022
01023 static doublereal gln[100] = { 0.,0.,.693147180559945309,
01024 1.791759469228055,3.17805383034794562,4.78749174278204599,
01025 6.579251212010101,8.5251613610654143,10.6046029027452502,
01026 12.8018274800814696,15.1044125730755153,17.5023078458738858,
01027 19.9872144956618861,22.5521638531234229,25.1912211827386815,
01028 27.8992713838408916,30.6718601060806728,33.5050734501368889,
01029 36.3954452080330536,39.339884187199494,42.335616460753485,
01030 45.380138898476908,48.4711813518352239,51.6066755677643736,
01031 54.7847293981123192,58.0036052229805199,61.261701761002002,
01032 64.5575386270063311,67.889743137181535,71.257038967168009,
01033 74.6582363488301644,78.0922235533153106,81.5579594561150372,
01034 85.0544670175815174,88.5808275421976788,92.1361756036870925,
01035 95.7196945421432025,99.3306124547874269,102.968198614513813,
01036 106.631760260643459,110.320639714757395,114.034211781461703,
01037 117.771881399745072,121.533081515438634,125.317271149356895,
01038 129.123933639127215,132.95257503561631,136.802722637326368,
01039 140.673923648234259,144.565743946344886,148.477766951773032,
01040 152.409592584497358,156.360836303078785,160.331128216630907,
01041 164.320112263195181,168.327445448427652,172.352797139162802,
01042 176.395848406997352,180.456291417543771,184.533828861449491,
01043 188.628173423671591,192.739047287844902,196.866181672889994,
01044 201.009316399281527,205.168199482641199,209.342586752536836,
01045 213.532241494563261,217.736934113954227,221.956441819130334,
01046 226.190548323727593,230.439043565776952,234.701723442818268,
01047 238.978389561834323,243.268849002982714,247.572914096186884,
01048 251.890402209723194,256.221135550009525,260.564940971863209,
01049 264.921649798552801,269.291097651019823,273.673124285693704,
01050 278.067573440366143,282.474292687630396,286.893133295426994,
01051 291.323950094270308,295.766601350760624,300.220948647014132,
01052 304.686856765668715,309.164193580146922,313.652829949879062,
01053 318.152639620209327,322.663499126726177,327.185287703775217,
01054 331.717887196928473,336.261181979198477,340.815058870799018,
01055 345.379407062266854,349.954118040770237,354.539085519440809,
01056 359.134205369575399 };
01057 static doublereal cf[22] = { .0833333333333333333,-.00277777777777777778,
01058 7.93650793650793651e-4,-5.95238095238095238e-4,
01059 8.41750841750841751e-4,-.00191752691752691753,
01060 .00641025641025641026,-.0295506535947712418,.179644372368830573,
01061 -1.39243221690590112,13.402864044168392,-156.848284626002017,
01062 2193.10333333333333,-36108.7712537249894,691472.268851313067,
01063 -15238221.5394074162,382900751.391414141,-10882266035.7843911,
01064 347320283765.002252,-12369602142269.2745,488788064793079.335,
01065 -21320333960919373.9 };
01066 static doublereal con = 1.83787706640934548;
01067
01068
01069 integer i__1;
01070 doublereal ret_val = 0;
01071
01072
01073
01074
01075 static doublereal zinc, zmin, zdmy;
01076 static integer i__, k;
01077 static doublereal s, wdtol;
01078 static doublereal t1, fz, zm;
01079 static integer mz, nz;
01080 static doublereal zp;
01081 static integer i1m;
01082 static doublereal fln, tlg, rln, trm, tst, zsq;
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127 *ierr = 0;
01128 if (*z__ <= 0.) {
01129 goto L70;
01130 }
01131 if (*z__ > 101.) {
01132 goto L10;
01133 }
01134 nz = (integer) ((real) (*z__));
01135 fz = *z__ - (real) nz;
01136 if (fz > 0.) {
01137 goto L10;
01138 }
01139 if (nz > 100) {
01140 goto L10;
01141 }
01142 ret_val = gln[nz - 1];
01143 return ret_val;
01144 L10:
01145 wdtol = d1mach_(&c__4);
01146 wdtol = max(wdtol,5e-19);
01147 i1m = i1mach_(&c__14);
01148 rln = d1mach_(&c__5) * (real) i1m;
01149 fln = min(rln,20.);
01150 fln = max(fln,3.);
01151 fln += -3.;
01152 zm = fln * .3875 + 1.8;
01153 mz = (integer) ((real) zm) + 1;
01154 zmin = (real) mz;
01155 zdmy = *z__;
01156 zinc = 0.;
01157 if (*z__ >= zmin) {
01158 goto L20;
01159 }
01160 zinc = zmin - (real) nz;
01161 zdmy = *z__ + zinc;
01162 L20:
01163 zp = 1. / zdmy;
01164 t1 = cf[0] * zp;
01165 s = t1;
01166 if (zp < wdtol) {
01167 goto L40;
01168 }
01169 zsq = zp * zp;
01170 tst = t1 * wdtol;
01171 for (k = 2; k <= 22; ++k) {
01172 zp *= zsq;
01173 trm = cf[k - 1] * zp;
01174 if (abs(trm) < tst) {
01175 goto L40;
01176 }
01177 s += trm;
01178
01179 }
01180 L40:
01181 if (zinc != 0.) {
01182 goto L50;
01183 }
01184 tlg = log(*z__);
01185 ret_val = *z__ * (tlg - 1.) + (con - tlg) * .5 + s;
01186 return ret_val;
01187 L50:
01188 zp = 1.;
01189 nz = (integer) ((real) zinc);
01190 i__1 = nz;
01191 for (i__ = 1; i__ <= i__1; ++i__) {
01192 zp *= *z__ + (real) (i__ - 1);
01193
01194 }
01195 tlg = log(zdmy);
01196 ret_val = zdmy * (tlg - 1.) - log(zp) + (con - tlg) * .5 + s;
01197 return ret_val;
01198
01199
01200 L70:
01201 *ierr = 1;
01202 return ret_val;
01203 }
01204
01205
01206 integer i1mach_(const integer* i__)
01207 {
01208
01209
01210 static struct {
01211 integer e_1[16];
01212 } equiv_0 = { {5, 6, 0, 0, 32, 4, 2, 31, 2147483647, 2, 24, -125, 127,
01213 53, -1021, 1023} };
01214
01215
01216
01217 static char fmt_9000[] = "(\0021ERROR 1 IN I1MACH - I OUT OF BOUND\
01218 S\002)";
01219
01220
01221 integer ret_val = 0;
01222
01223
01224
01225
01226 #define imach ((integer *)&equiv_0)
01227 #define output ((integer *)&equiv_0 + 3)
01228
01229
01230 static cilist io___72 = { 0, 0, 0, fmt_9000, 0 };
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539
01540
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574
01575
01576
01577
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600
01601
01602
01603
01604
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734
01735
01736
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
01748
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871 if (*i__ < 1 || *i__ > 16) {
01872 goto L10;
01873 }
01874
01875 ret_val = imach[*i__ - 1];
01876 return ret_val;
01877
01878 L10:
01879 io___72.ciunit = *output;
01880 s_wsfe(&io___72);
01881 e_wsfe();
01882
01883
01884
01885
01886 s_stop("", (ftnlen)0);
01887 return ret_val;
01888 }
01889
01890 #undef output
01891 #undef imach
01892
01893
01894
01895 int xerror_(char* mess, integer* nmess,
01896 integer* l1, integer* l2, ftnlen mess_len)
01897 {
01898
01899 static char fmt_900[] = "(/)";
01900
01901
01902 integer i__1, i__2;
01903
01904
01905
01906
01907 static integer kmin, i__, k, nn, nr;
01908
01909
01910 static cilist io___76 = { 0, 6, 0, fmt_900, 0 };
01911 static cilist io___79 = { 0, 6, 0, 0, 0 };
01912 static cilist io___80 = { 0, 6, 0, fmt_900, 0 };
01913
01914
01915
01916
01917
01918
01919
01920
01921 nn = *nmess / 70;
01922 nr = *nmess - nn * 70;
01923 if (nr != 0) {
01924 ++nn;
01925 }
01926 k = 1;
01927 s_wsfe(&io___76);
01928 e_wsfe();
01929 i__1 = nn;
01930 for (i__ = 1; i__ <= i__1; ++i__) {
01931
01932 i__2 = k + 69;
01933 kmin = min(i__2,*nmess);
01934 s_wsle(&io___79);
01935 do_lio(&c__9, &c__1, mess + (k - 1), kmin - (k - 1));
01936 e_wsle();
01937 k += 70;
01938
01939 }
01940 s_wsfe(&io___80);
01941 e_wsfe();
01942 return 0;
01943 }
01944
01945 doublereal myzabs_(const doublereal* zr, const doublereal* zi)
01946 {
01947
01948 doublereal ret_val;
01949
01950
01951
01952
01953 static doublereal q, s, u, v;
01954
01955
01956
01957
01958
01959
01960
01961
01962
01963 u = abs(*zr);
01964 v = abs(*zi);
01965 s = u + v;
01966
01967
01968
01969
01970 s *= 1.;
01971 if (s == 0.) {
01972 goto L20;
01973 }
01974 if (u > v) {
01975 goto L10;
01976 }
01977 q = u / v;
01978 ret_val = v * sqrt(q * q + 1.);
01979 return ret_val;
01980 L10:
01981 q = v / u;
01982 ret_val = u * sqrt(q * q + 1.);
01983 return ret_val;
01984 L20:
01985 ret_val = 0.;
01986 return ret_val;
01987 }
01988
01989
01990 int zacai_(doublereal* zr, doublereal* zi, doublereal* fnu,
01991 integer* kode, integer* mr, integer* n,
01992 doublereal* yr, doublereal* yi,
01993 integer* nz,
01994 doublereal* rl, doublereal*tol,
01995 doublereal* elim, doublereal* alim)
01996 {
01997
01998
01999 static doublereal pi = 3.14159265358979324;
02000
02001
02002
02003
02004 static doublereal dfnu;
02005 static doublereal az;
02006 static integer nn, nw;
02007 static doublereal yy, c1i, c2i;
02008 static doublereal c1r, c2r, arg;
02009 static integer iuf;
02010 static doublereal cyi[2], fmr, sgn;
02011 static integer inu;
02012 static doublereal cyr[2], zni, znr;
02013
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029
02030
02031
02032 --yi;
02033 --yr;
02034
02035
02036 *nz = 0;
02037 znr = -(*zr);
02038 zni = -(*zi);
02039 az = myzabs_(zr, zi);
02040 nn = *n;
02041 dfnu = *fnu + (doublereal) ((real) (*n - 1));
02042 if (az <= 2.) {
02043 goto L10;
02044 }
02045 if (az * az * .25 > dfnu + 1.) {
02046 goto L20;
02047 }
02048 L10:
02049
02050
02051
02052 zseri_(&znr, &zni, fnu, kode, &nn, &yr[1], &yi[1], &nw, tol, elim, alim);
02053 goto L40;
02054 L20:
02055 if (az < *rl) {
02056 goto L30;
02057 }
02058
02059
02060
02061 zasyi_(&znr, &zni, fnu, kode, &nn, &yr[1], &yi[1], &nw, rl, tol, elim,
02062 alim);
02063 if (nw < 0) {
02064 goto L80;
02065 }
02066 goto L40;
02067 L30:
02068
02069
02070
02071 zmlri_(&znr, &zni, fnu, kode, &nn, &yr[1], &yi[1], &nw, tol);
02072 if (nw < 0) {
02073 goto L80;
02074 }
02075 L40:
02076
02077
02078
02079 zbknu_(&znr, &zni, fnu, kode, &c__1, cyr, cyi, &nw, tol, elim, alim);
02080 if (nw != 0) {
02081 goto L80;
02082 }
02083 fmr = (doublereal) ((real) (*mr));
02084 sgn = -d_sign(&pi, &fmr);
02085 csgnr = 0.;
02086 csgni = sgn;
02087 if (*kode == 1) {
02088 goto L50;
02089 }
02090 yy = -zni;
02091 csgnr = -csgni * sin(yy);
02092 csgni *= cos(yy);
02093 L50:
02094
02095
02096
02097
02098 inu = (integer) ((real) (*fnu));
02099 arg = (*fnu - (doublereal) ((real) inu)) * sgn;
02100 cspnr = cos(arg);
02101 cspni = sin(arg);
02102 if (inu % 2 == 0) {
02103 goto L60;
02104 }
02105 cspnr = -cspnr;
02106 cspni = -cspni;
02107 L60:
02108 c1r = cyr[0];
02109 c1i = cyi[0];
02110 c2r = yr[1];
02111 c2i = yi[1];
02112 if (*kode == 1) {
02113 goto L70;
02114 }
02115 iuf = 0;
02116 ascle = d1mach_(&c__1) * 1e3 / *tol;
02117 zs1s2_(&znr, &zni, &c1r, &c1i, &c2r, &c2i, &nw, &ascle, alim, &iuf);
02118 *nz += nw;
02119 L70:
02120 yr[1] = cspnr * c1r - cspni * c1i + csgnr * c2r - csgni * c2i;
02121 yi[1] = cspnr * c1i + cspni * c1r + csgnr * c2i + csgni * c2r;
02122 return 0;
02123 L80:
02124 *nz = -1;
02125 if (nw == -2) {
02126 *nz = -2;
02127 }
02128 return 0;
02129 }
02130
02131 int zacon_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *mr, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *rl, doublereal *fnul, doublereal *tol, doublereal *elim, doublereal *alim)
02132 {
02133
02134
02135 static doublereal pi = 3.14159265358979324;
02136 static doublereal zeror = 0.;
02137 static doublereal coner = 1.;
02138
02139
02140 integer i__1;
02141
02142
02143
02144
02145 static doublereal cscl, cscr, csrr[3], cssr[3], razn;
02146 static integer i__, kflag;
02147 static doublereal ascle, bscle, csgni, csgnr, cspni, cspnr;
02148 static doublereal fn;
02149 static integer nn, nw;
02150 static doublereal yy, c1i, c2i, c1m;
02151 static doublereal as2, c1r, c2r, s1i, s2i, s1r, s2r, cki, arg, ckr, cpn;
02152 static integer iuf;
02153 static doublereal cyi[2], fmr, csr, azn, sgn;
02154 static integer inu;
02155 static doublereal bry[3], cyr[2], pti, spn, sti, zni, rzi, ptr, str, znr,
02156 rzr, sc1i, sc2i, sc1r, sc2r;
02157
02158
02159
02160
02161
02162
02163
02164
02165
02166
02167
02168
02169
02170
02171
02172
02173
02174 --yi;
02175 --yr;
02176
02177
02178 *nz = 0;
02179 znr = -(*zr);
02180 zni = -(*zi);
02181 nn = *n;
02182 zbinu_(&znr, &zni, fnu, kode, &nn, &yr[1], &yi[1], &nw, rl, fnul, tol,
02183 elim, alim);
02184 if (nw < 0) {
02185 goto L90;
02186 }
02187
02188
02189
02190 nn = min(2,*n);
02191 zbknu_(&znr, &zni, fnu, kode, &nn, cyr, cyi, &nw, tol, elim, alim);
02192 if (nw != 0) {
02193 goto L90;
02194 }
02195 s1r = cyr[0];
02196 s1i = cyi[0];
02197 fmr = (doublereal) ((real) (*mr));
02198 sgn = -d_sign(&pi, &fmr);
02199 csgnr = zeror;
02200 csgni = sgn;
02201 if (*kode == 1) {
02202 goto L10;
02203 }
02204 yy = -zni;
02205 cpn = cos(yy);
02206 spn = sin(yy);
02207 zmlt_(&csgnr, &csgni, &cpn, &spn, &csgnr, &csgni);
02208 L10:
02209
02210
02211
02212
02213 inu = (integer) ((real) (*fnu));
02214 arg = (*fnu - (doublereal) ((real) inu)) * sgn;
02215 cpn = cos(arg);
02216 spn = sin(arg);
02217 cspnr = cpn;
02218 cspni = spn;
02219 if (inu % 2 == 0) {
02220 goto L20;
02221 }
02222 cspnr = -cspnr;
02223 cspni = -cspni;
02224 L20:
02225 iuf = 0;
02226 c1r = s1r;
02227 c1i = s1i;
02228 c2r = yr[1];
02229 c2i = yi[1];
02230 ascle = d1mach_(&c__1) * 1e3 / *tol;
02231 if (*kode == 1) {
02232 goto L30;
02233 }
02234 zs1s2_(&znr, &zni, &c1r, &c1i, &c2r, &c2i, &nw, &ascle, alim, &iuf);
02235 *nz += nw;
02236 sc1r = c1r;
02237 sc1i = c1i;
02238 L30:
02239 zmlt_(&cspnr, &cspni, &c1r, &c1i, &str, &sti);
02240 zmlt_(&csgnr, &csgni, &c2r, &c2i, &ptr, &pti);
02241 yr[1] = str + ptr;
02242 yi[1] = sti + pti;
02243 if (*n == 1) {
02244 return 0;
02245 }
02246 cspnr = -cspnr;
02247 cspni = -cspni;
02248 s2r = cyr[1];
02249 s2i = cyi[1];
02250 c1r = s2r;
02251 c1i = s2i;
02252 c2r = yr[2];
02253 c2i = yi[2];
02254 if (*kode == 1) {
02255 goto L40;
02256 }
02257 zs1s2_(&znr, &zni, &c1r, &c1i, &c2r, &c2i, &nw, &ascle, alim, &iuf);
02258 *nz += nw;
02259 sc2r = c1r;
02260 sc2i = c1i;
02261 L40:
02262 zmlt_(&cspnr, &cspni, &c1r, &c1i, &str, &sti);
02263 zmlt_(&csgnr, &csgni, &c2r, &c2i, &ptr, &pti);
02264 yr[2] = str + ptr;
02265 yi[2] = sti + pti;
02266 if (*n == 2) {
02267 return 0;
02268 }
02269 cspnr = -cspnr;
02270 cspni = -cspni;
02271 azn = myzabs_(&znr, &zni);
02272 razn = 1. / azn;
02273 str = znr * razn;
02274 sti = -zni * razn;
02275 rzr = (str + str) * razn;
02276 rzi = (sti + sti) * razn;
02277 fn = *fnu + 1.;
02278 ckr = fn * rzr;
02279 cki = fn * rzi;
02280
02281
02282
02283 cscl = 1. / *tol;
02284 cscr = *tol;
02285 cssr[0] = cscl;
02286 cssr[1] = coner;
02287 cssr[2] = cscr;
02288 csrr[0] = cscr;
02289 csrr[1] = coner;
02290 csrr[2] = cscl;
02291 bry[0] = ascle;
02292 bry[1] = 1. / ascle;
02293 bry[2] = d1mach_(&c__2);
02294 as2 = myzabs_(&s2r, &s2i);
02295 kflag = 2;
02296 if (as2 > bry[0]) {
02297 goto L50;
02298 }
02299 kflag = 1;
02300 goto L60;
02301 L50:
02302 if (as2 < bry[1]) {
02303 goto L60;
02304 }
02305 kflag = 3;
02306 L60:
02307 bscle = bry[kflag - 1];
02308 s1r *= cssr[kflag - 1];
02309 s1i *= cssr[kflag - 1];
02310 s2r *= cssr[kflag - 1];
02311 s2i *= cssr[kflag - 1];
02312 csr = csrr[kflag - 1];
02313 i__1 = *n;
02314 for (i__ = 3; i__ <= i__1; ++i__) {
02315 str = s2r;
02316 sti = s2i;
02317 s2r = ckr * str - cki * sti + s1r;
02318 s2i = ckr * sti + cki * str + s1i;
02319 s1r = str;
02320 s1i = sti;
02321 c1r = s2r * csr;
02322 c1i = s2i * csr;
02323 str = c1r;
02324 sti = c1i;
02325 c2r = yr[i__];
02326 c2i = yi[i__];
02327 if (*kode == 1) {
02328 goto L70;
02329 }
02330 if (iuf < 0) {
02331 goto L70;
02332 }
02333 zs1s2_(&znr, &zni, &c1r, &c1i, &c2r, &c2i, &nw, &ascle, alim, &iuf);
02334 *nz += nw;
02335 sc1r = sc2r;
02336 sc1i = sc2i;
02337 sc2r = c1r;
02338 sc2i = c1i;
02339 if (iuf != 3) {
02340 goto L70;
02341 }
02342 iuf = -4;
02343 s1r = sc1r * cssr[kflag - 1];
02344 s1i = sc1i * cssr[kflag - 1];
02345 s2r = sc2r * cssr[kflag - 1];
02346 s2i = sc2i * cssr[kflag - 1];
02347 str = sc2r;
02348 sti = sc2i;
02349 L70:
02350 ptr = cspnr * c1r - cspni * c1i;
02351 pti = cspnr * c1i + cspni * c1r;
02352 yr[i__] = ptr + csgnr * c2r - csgni * c2i;
02353 yi[i__] = pti + csgnr * c2i + csgni * c2r;
02354 ckr += rzr;
02355 cki += rzi;
02356 cspnr = -cspnr;
02357 cspni = -cspni;
02358 if (kflag >= 3) {
02359 goto L80;
02360 }
02361 ptr = abs(c1r);
02362 pti = abs(c1i);
02363 c1m = max(ptr,pti);
02364 if (c1m <= bscle) {
02365 goto L80;
02366 }
02367 ++kflag;
02368 bscle = bry[kflag - 1];
02369 s1r *= csr;
02370 s1i *= csr;
02371 s2r = str;
02372 s2i = sti;
02373 s1r *= cssr[kflag - 1];
02374 s1i *= cssr[kflag - 1];
02375 s2r *= cssr[kflag - 1];
02376 s2i *= cssr[kflag - 1];
02377 csr = csrr[kflag - 1];
02378 L80:
02379 ;
02380 }
02381 return 0;
02382 L90:
02383 *nz = -1;
02384 if (nw == -2) {
02385 *nz = -2;
02386 }
02387 return 0;
02388 }
02389
02390
02391 int zasyi_(doublereal *zr, doublereal *zi, doublereal *fnu,
02392 integer *kode, integer *n,
02393 doublereal *yr, doublereal *yi, integer *nz,
02394 doublereal *rl, doublereal *tol,
02395 doublereal *elim, doublereal *alim)
02396 {
02397
02398
02399 static doublereal pi = 3.14159265358979324;
02400 static doublereal rtpi = .159154943091895336;
02401 static doublereal zeror = 0.;
02402 static doublereal zeroi = 0.;
02403 static doublereal coner = 1.;
02404 static doublereal conei = 0.;
02405
02406
02407 integer i__1, i__2;
02408 doublereal d__1, d__2;
02409
02410
02411
02412
02413 static doublereal dfnu, atol;
02414 static integer i__, j, k, m;
02415 static doublereal s;
02416 static integer koded;
02417 static doublereal aa, bb;
02418 static integer ib;
02419 static doublereal ak, bk;
02420 static integer il, jl;
02421 static doublereal az;
02422 static integer nn;
02423 static doublereal p1i, s2i, p1r, s2r, cki, dki, fdn, arg, aez, arm, ckr,
02424 dkr, czi, ezi, sgn;
02425 static integer inu;
02426 static doublereal raz, czr, ezr, sqk, sti, rzi, tzi, str, rzr, tzr, ak1i,
02427 ak1r, cs1i, cs2i, cs1r, cs2r, dnu2, rtr1;
02428
02429
02430
02431
02432
02433
02434
02435
02436
02437
02438
02439
02440
02441 --yi;
02442 --yr;
02443
02444
02445
02446 *nz = 0;
02447 az = myzabs_(zr, zi);
02448 arm = d1mach_(&c__1) * 1e3;
02449 rtr1 = sqrt(arm);
02450 il = min(2,*n);
02451 dfnu = *fnu + (doublereal) ((real) (*n - il));
02452
02453
02454
02455 raz = 1. / az;
02456 str = *zr * raz;
02457 sti = -(*zi) * raz;
02458 ak1r = rtpi * str * raz;
02459 ak1i = rtpi * sti * raz;
02460 zsqrt_(&ak1r, &ak1i, &ak1r, &ak1i);
02461 czr = *zr;
02462 czi = *zi;
02463 if (*kode != 2) {
02464 goto L10;
02465 }
02466 czr = zeror;
02467 czi = *zi;
02468 L10:
02469 if (abs(czr) > *elim) {
02470 goto L100;
02471 }
02472 dnu2 = dfnu + dfnu;
02473 koded = 1;
02474 if (abs(czr) > *alim && *n > 2) {
02475 goto L20;
02476 }
02477 koded = 0;
02478 zexp_(&czr, &czi, &str, &sti);
02479 zmlt_(&ak1r, &ak1i, &str, &sti, &ak1r, &ak1i);
02480 L20:
02481 fdn = 0.;
02482 if (dnu2 > rtr1) {
02483 fdn = dnu2 * dnu2;
02484 }
02485 ezr = *zr * 8.;
02486 ezi = *zi * 8.;
02487
02488
02489
02490
02491
02492 aez = az * 8.;
02493 s = *tol / aez;
02494 jl = (integer) ((real) (*rl + *rl)) + 2;
02495 p1r = zeror;
02496 p1i = zeroi;
02497 if (*zi == 0.) {
02498 goto L30;
02499 }
02500
02501
02502
02503
02504 inu = (integer) ((real) (*fnu));
02505 arg = (*fnu - (doublereal) ((real) inu)) * pi;
02506 inu = inu + *n - il;
02507 ak = -sin(arg);
02508 bk = cos(arg);
02509 if (*zi < 0.) {
02510 bk = -bk;
02511 }
02512 p1r = ak;
02513 p1i = bk;
02514 if (inu % 2 == 0) {
02515 goto L30;
02516 }
02517 p1r = -p1r;
02518 p1i = -p1i;
02519 L30:
02520 i__1 = il;
02521 for (k = 1; k <= i__1; ++k) {
02522 sqk = fdn - 1.;
02523 atol = s * abs(sqk);
02524 sgn = 1.;
02525 cs1r = coner;
02526 cs1i = conei;
02527 cs2r = coner;
02528 cs2i = conei;
02529 ckr = coner;
02530 cki = conei;
02531 ak = 0.;
02532 aa = 1.;
02533 bb = aez;
02534 dkr = ezr;
02535 dki = ezi;
02536 i__2 = jl;
02537 for (j = 1; j <= i__2; ++j) {
02538 zdiv_(&ckr, &cki, &dkr, &dki, &str, &sti);
02539 ckr = str * sqk;
02540 cki = sti * sqk;
02541 cs2r += ckr;
02542 cs2i += cki;
02543 sgn = -sgn;
02544 cs1r += ckr * sgn;
02545 cs1i += cki * sgn;
02546 dkr += ezr;
02547 dki += ezi;
02548 aa = aa * abs(sqk) / bb;
02549 bb += aez;
02550 ak += 8.;
02551 sqk -= ak;
02552 if (aa <= atol) {
02553 goto L50;
02554 }
02555
02556 }
02557 goto L110;
02558 L50:
02559 s2r = cs1r;
02560 s2i = cs1i;
02561 if (*zr + *zr >= *elim) {
02562 goto L60;
02563 }
02564 tzr = *zr + *zr;
02565 tzi = *zi + *zi;
02566 d__1 = -tzr;
02567 d__2 = -tzi;
02568 zexp_(&d__1, &d__2, &str, &sti);
02569 zmlt_(&str, &sti, &p1r, &p1i, &str, &sti);
02570 zmlt_(&str, &sti, &cs2r, &cs2i, &str, &sti);
02571 s2r += str;
02572 s2i += sti;
02573 L60:
02574 fdn = fdn + dfnu * 8. + 4.;
02575 p1r = -p1r;
02576 p1i = -p1i;
02577 m = *n - il + k;
02578 yr[m] = s2r * ak1r - s2i * ak1i;
02579 yi[m] = s2r * ak1i + s2i * ak1r;
02580
02581 }
02582 if (*n <= 2) {
02583 return 0;
02584 }
02585 nn = *n;
02586 k = nn - 2;
02587 ak = (doublereal) ((real) k);
02588 str = *zr * raz;
02589 sti = -(*zi) * raz;
02590 rzr = (str + str) * raz;
02591 rzi = (sti + sti) * raz;
02592 ib = 3;
02593 i__1 = nn;
02594 for (i__ = ib; i__ <= i__1; ++i__) {
02595 yr[k] = (ak + *fnu) * (rzr * yr[k + 1] - rzi * yi[k + 1]) + yr[k + 2];
02596 yi[k] = (ak + *fnu) * (rzr * yi[k + 1] + rzi * yr[k + 1]) + yi[k + 2];
02597 ak += -1.;
02598 --k;
02599
02600 }
02601 if (koded == 0) {
02602 return 0;
02603 }
02604 zexp_(&czr, &czi, &ckr, &cki);
02605 i__1 = nn;
02606 for (i__ = 1; i__ <= i__1; ++i__) {
02607 str = yr[i__] * ckr - yi[i__] * cki;
02608 yi[i__] = yr[i__] * cki + yi[i__] * ckr;
02609 yr[i__] = str;
02610
02611 }
02612 return 0;
02613 L100:
02614 *nz = -1;
02615 return 0;
02616 L110:
02617 *nz = -2;
02618 return 0;
02619 }
02620
02621 int zbinu_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *n, doublereal *cyr, doublereal *cyi, integer *nz, doublereal *rl, doublereal *fnul, doublereal *tol, doublereal *elim, doublereal *alim)
02622 {
02623
02624
02625 static doublereal zeror = 0.;
02626 static doublereal zeroi = 0.;
02627
02628
02629 integer i__1;
02630
02631
02632 static doublereal dfnu;
02633 static integer i__, nlast;
02634 static doublereal az;
02635 static integer nn, nw;
02636 static doublereal cwi[2], cwr[2];
02637 static integer nui, inw;
02638
02639
02640
02641
02642
02643
02644
02645
02646
02647 --cyi;
02648 --cyr;
02649
02650
02651
02652 *nz = 0;
02653 az = myzabs_(zr, zi);
02654 nn = *n;
02655 dfnu = *fnu + (doublereal) ((real) (*n - 1));
02656 if (az <= 2.) {
02657 goto L10;
02658 }
02659 if (az * az * .25 > dfnu + 1.) {
02660 goto L20;
02661 }
02662 L10:
02663
02664
02665
02666 zseri_(zr, zi, fnu, kode, &nn, &cyr[1], &cyi[1], &nw, tol, elim, alim);
02667 inw = abs(nw);
02668 *nz += inw;
02669 nn -= inw;
02670 if (nn == 0) {
02671 return 0;
02672 }
02673 if (nw >= 0) {
02674 goto L120;
02675 }
02676 dfnu = *fnu + (doublereal) ((real) (nn - 1));
02677 L20:
02678 if (az < *rl) {
02679 goto L40;
02680 }
02681 if (dfnu <= 1.) {
02682 goto L30;
02683 }
02684 if (az + az < dfnu * dfnu) {
02685 goto L50;
02686 }
02687
02688
02689
02690 L30:
02691 zasyi_(zr, zi, fnu, kode, &nn, &cyr[1], &cyi[1], &nw, rl, tol, elim, alim)
02692 ;
02693 if (nw < 0) {
02694 goto L130;
02695 }
02696 goto L120;
02697 L40:
02698 if (dfnu <= 1.) {
02699 goto L70;
02700 }
02701 L50:
02702
02703
02704
02705 zuoik_(zr, zi, fnu, kode, &c__1, &nn, &cyr[1], &cyi[1], &nw, tol, elim,
02706 alim);
02707 if (nw < 0) {
02708 goto L130;
02709 }
02710 *nz += nw;
02711 nn -= nw;
02712 if (nn == 0) {
02713 return 0;
02714 }
02715 dfnu = *fnu + (doublereal) ((real) (nn - 1));
02716 if (dfnu > *fnul) {
02717 goto L110;
02718 }
02719 if (az > *fnul) {
02720 goto L110;
02721 }
02722 L60:
02723 if (az > *rl) {
02724 goto L80;
02725 }
02726 L70:
02727
02728
02729
02730 zmlri_(zr, zi, fnu, kode, &nn, &cyr[1], &cyi[1], &nw, tol);
02731 if (nw < 0) {
02732 goto L130;
02733 }
02734 goto L120;
02735 L80:
02736
02737
02738
02739
02740
02741
02742 zuoik_(zr, zi, fnu, kode, &c__2, &c__2, cwr, cwi, &nw, tol, elim, alim);
02743 if (nw >= 0) {
02744 goto L100;
02745 }
02746 *nz = nn;
02747 i__1 = nn;
02748 for (i__ = 1; i__ <= i__1; ++i__) {
02749 cyr[i__] = zeror;
02750 cyi[i__] = zeroi;
02751
02752 }
02753 return 0;
02754 L100:
02755 if (nw > 0) {
02756 goto L130;
02757 }
02758 zwrsk_(zr, zi, fnu, kode, &nn, &cyr[1], &cyi[1], &nw, cwr, cwi, tol, elim,
02759 alim);
02760 if (nw < 0) {
02761 goto L130;
02762 }
02763 goto L120;
02764 L110:
02765
02766
02767
02768 nui = (integer) ((real) (*fnul - dfnu)) + 1;
02769 nui = max(nui,0);
02770 zbuni_(zr, zi, fnu, kode, &nn, &cyr[1], &cyi[1], &nw, &nui, &nlast, fnul,
02771 tol, elim, alim);
02772 if (nw < 0) {
02773 goto L130;
02774 }
02775 *nz += nw;
02776 if (nlast == 0) {
02777 goto L120;
02778 }
02779 nn = nlast;
02780 goto L60;
02781 L120:
02782 return 0;
02783 L130:
02784 *nz = -1;
02785 if (nw == -2) {
02786 *nz = -2;
02787 }
02788 return 0;
02789 }
02790
02791
02792 int zbknu_(doublereal *zr, doublereal *zi, doublereal *fnu,
02793 integer *kode, integer *n,
02794 doublereal *yr, doublereal *yi,
02795 integer *nz, doublereal *tol,
02796 doublereal *elim, doublereal *alim)
02797 {
02798
02799
02800 static integer kmax = 30;
02801 static doublereal czeror = 0.;
02802 static doublereal czeroi = 0.;
02803 static doublereal coner = 1.;
02804 static doublereal conei = 0.;
02805 static doublereal ctwor = 2.;
02806 static doublereal r1 = 2.;
02807 static doublereal dpi = 3.14159265358979324;
02808 static doublereal rthpi = 1.25331413731550025;
02809 static doublereal spi = 1.90985931710274403;
02810 static doublereal hpi = 1.57079632679489662;
02811 static doublereal fpi = 1.89769999331517738;
02812 static doublereal tth = .666666666666666666;
02813 static doublereal cc[8] = { .577215664901532861,-.0420026350340952355,
02814 -.0421977345555443367,.00721894324666309954,
02815 -2.15241674114950973e-4,-2.01348547807882387e-5,
02816 1.13302723198169588e-6,6.11609510448141582e-9 };
02817
02818
02819 integer i__1;
02820 doublereal d__1;
02821
02822
02823
02824
02825 static doublereal cchi, cchr, alas, cshi;
02826 static integer inub, idum;
02827 static doublereal cshr, fmui, rcaz, csrr[3], cssr[3], fmur;
02828 static doublereal smui;
02829 static doublereal smur;
02830 static integer i__, j, k, iflag;
02831 static doublereal s;
02832 static integer kflag;
02833 static doublereal coefi;
02834 static integer koded;
02835 static doublereal ascle, coefr, helim, celmr, csclr, crscr;
02836 static doublereal a1, a2, etest;
02837 static doublereal g1, g2;
02838 static doublereal t1, t2;
02839 static doublereal aa, bb, fc, ak, bk;
02840 static integer ic;
02841 static doublereal fi, fk, as;
02842 static integer kk;
02843 static doublereal fr, pi, qi, tm, pr, qr;
02844 static integer nw;
02845 static doublereal p1i, p2i, s1i, s2i, p2m, p1r, p2r, s1r, s2r, cbi, cbr,
02846 cki, caz, csi, ckr, fhs, fks, rak, czi, dnu, csr, elm, zdi, bry[3]
02847 , pti, czr, sti, zdr, cyr[2], rzi, ptr, cyi[2];
02848 static integer inu;
02849 static doublereal str, rzr, dnu2;
02850
02851
02852
02853
02854
02855
02856
02857
02858
02859
02860
02861
02862
02863
02864 --yi;
02865 --yr;
02866
02867
02868
02869 caz = myzabs_(zr, zi);
02870 csclr = 1. / *tol;
02871 crscr = *tol;
02872 cssr[0] = csclr;
02873 cssr[1] = 1.;
02874 cssr[2] = crscr;
02875 csrr[0] = crscr;
02876 csrr[1] = 1.;
02877 csrr[2] = csclr;
02878 bry[0] = d1mach_(&c__1) * 1e3 / *tol;
02879 bry[1] = 1. / bry[0];
02880 bry[2] = d1mach_(&c__2);
02881 *nz = 0;
02882 iflag = 0;
02883 koded = *kode;
02884 rcaz = 1. / caz;
02885 str = *zr * rcaz;
02886 sti = -(*zi) * rcaz;
02887 rzr = (str + str) * rcaz;
02888 rzi = (sti + sti) * rcaz;
02889 inu = (integer) ((real) (*fnu + .5));
02890 dnu = *fnu - (doublereal) ((real) inu);
02891 if (abs(dnu) == .5) {
02892 goto L110;
02893 }
02894 dnu2 = 0.;
02895 if (abs(dnu) > *tol) {
02896 dnu2 = dnu * dnu;
02897 }
02898 if (caz > r1) {
02899 goto L110;
02900 }
02901
02902
02903
02904 fc = 1.;
02905 zlog_(&rzr, &rzi, &smur, &smui, &idum);
02906 fmur = smur * dnu;
02907 fmui = smui * dnu;
02908 zshch_(&fmur, &fmui, &cshr, &cshi, &cchr, &cchi);
02909 if (dnu == 0.) {
02910 goto L10;
02911 }
02912 fc = dnu * dpi;
02913 fc /= sin(fc);
02914 smur = cshr / dnu;
02915 smui = cshi / dnu;
02916 L10:
02917 a2 = dnu + 1.;
02918
02919
02920
02921 t2 = exp(-dgamln_(&a2, &idum));
02922 t1 = 1. / (t2 * fc);
02923 if (abs(dnu) > .1) {
02924 goto L40;
02925 }
02926
02927
02928
02929 ak = 1.;
02930 s = cc[0];
02931 for (k = 2; k <= 8; ++k) {
02932 ak *= dnu2;
02933 tm = cc[k - 1] * ak;
02934 s += tm;
02935 if (abs(tm) < *tol) {
02936 goto L30;
02937 }
02938
02939 }
02940 L30:
02941 g1 = -s;
02942 goto L50;
02943 L40:
02944 g1 = (t1 - t2) / (dnu + dnu);
02945 L50:
02946 g2 = (t1 + t2) * .5;
02947 fr = fc * (cchr * g1 + smur * g2);
02948 fi = fc * (cchi * g1 + smui * g2);
02949 zexp_(&fmur, &fmui, &str, &sti);
02950 pr = str * .5 / t2;
02951 pi = sti * .5 / t2;
02952 zdiv_(&c_b147, &c_b148, &str, &sti, &ptr, &pti);
02953 qr = ptr / t1;
02954 qi = pti / t1;
02955 s1r = fr;
02956 s1i = fi;
02957 s2r = pr;
02958 s2i = pi;
02959 ak = 1.;
02960 a1 = 1.;
02961 ckr = coner;
02962 cki = conei;
02963 bk = 1. - dnu2;
02964 if (inu > 0 || *n > 1) {
02965 goto L80;
02966 }
02967
02968
02969
02970 if (caz < *tol) {
02971 goto L70;
02972 }
02973 zmlt_(zr, zi, zr, zi, &czr, &czi);
02974 czr *= .25;
02975 czi *= .25;
02976 t1 = caz * .25 * caz;
02977 L60:
02978 fr = (fr * ak + pr + qr) / bk;
02979 fi = (fi * ak + pi + qi) / bk;
02980 str = 1. / (ak - dnu);
02981 pr *= str;
02982 pi *= str;
02983 str = 1. / (ak + dnu);
02984 qr *= str;
02985 qi *= str;
02986 str = ckr * czr - cki * czi;
02987 rak = 1. / ak;
02988 cki = (ckr * czi + cki * czr) * rak;
02989 ckr = str * rak;
02990 s1r = ckr * fr - cki * fi + s1r;
02991 s1i = ckr * fi + cki * fr + s1i;
02992 a1 = a1 * t1 * rak;
02993 bk = bk + ak + ak + 1.;
02994 ak += 1.;
02995 if (a1 > *tol) {
02996 goto L60;
02997 }
02998 L70:
02999 yr[1] = s1r;
03000 yi[1] = s1i;
03001 if (koded == 1) {
03002 return 0;
03003 }
03004 zexp_(zr, zi, &str, &sti);
03005 zmlt_(&s1r, &s1i, &str, &sti, &yr[1], &yi[1]);
03006 return 0;
03007
03008
03009
03010 L80:
03011 if (caz < *tol) {
03012 goto L100;
03013 }
03014 zmlt_(zr, zi, zr, zi, &czr, &czi);
03015 czr *= .25;
03016 czi *= .25;
03017 t1 = caz * .25 * caz;
03018 L90:
03019 fr = (fr * ak + pr + qr) / bk;
03020 fi = (fi * ak + pi + qi) / bk;
03021 str = 1. / (ak - dnu);
03022 pr *= str;
03023 pi *= str;
03024 str = 1. / (ak + dnu);
03025 qr *= str;
03026 qi *= str;
03027 str = ckr * czr - cki * czi;
03028 rak = 1. / ak;
03029 cki = (ckr * czi + cki * czr) * rak;
03030 ckr = str * rak;
03031 s1r = ckr * fr - cki * fi + s1r;
03032 s1i = ckr * fi + cki * fr + s1i;
03033 str = pr - fr * ak;
03034 sti = pi - fi * ak;
03035 s2r = ckr * str - cki * sti + s2r;
03036 s2i = ckr * sti + cki * str + s2i;
03037 a1 = a1 * t1 * rak;
03038 bk = bk + ak + ak + 1.;
03039 ak += 1.;
03040 if (a1 > *tol) {
03041 goto L90;
03042 }
03043 L100:
03044 kflag = 2;
03045 a1 = *fnu + 1.;
03046 ak = a1 * abs(smur);
03047 if (ak > *alim) {
03048 kflag = 3;
03049 }
03050 str = cssr[kflag - 1];
03051 p2r = s2r * str;
03052 p2i = s2i * str;
03053 zmlt_(&p2r, &p2i, &rzr, &rzi, &s2r, &s2i);
03054 s1r *= str;
03055 s1i *= str;
03056 if (koded == 1) {
03057 goto L210;
03058 }
03059 zexp_(zr, zi, &fr, &fi);
03060 zmlt_(&s1r, &s1i, &fr, &fi, &s1r, &s1i);
03061 zmlt_(&s2r, &s2i, &fr, &fi, &s2r, &s2i);
03062 goto L210;
03063
03064
03065
03066
03067
03068
03069 L110:
03070 zsqrt_(zr, zi, &str, &sti);
03071 zdiv_(&rthpi, &czeroi, &str, &sti, &coefr, &coefi);
03072 kflag = 2;
03073 if (koded == 2) {
03074 goto L120;
03075 }
03076 if (*zr > *alim) {
03077 goto L290;
03078 }
03079
03080 str = exp(-(*zr)) * cssr[kflag - 1];
03081 sti = -str * sin(*zi);
03082 str *= cos(*zi);
03083 zmlt_(&coefr, &coefi, &str, &sti, &coefr, &coefi);
03084 L120:
03085 if (abs(dnu) == .5) {
03086 goto L300;
03087 }
03088
03089
03090
03091 ak = cos(dpi * dnu);
03092 ak = abs(ak);
03093 if (ak == czeror) {
03094 goto L300;
03095 }
03096 fhs = (d__1 = .25 - dnu2, abs(d__1));
03097 if (fhs == czeror) {
03098 goto L300;
03099 }
03100
03101
03102
03103
03104
03105
03106 t1 = (doublereal) ((real) (i1mach_(&c__14) - 1));
03107 t1 = t1 * d1mach_(&c__5) * 3.321928094;
03108 t1 = max(t1,12.);
03109 t1 = min(t1,60.);
03110 t2 = tth * t1 - 6.;
03111 if (*zr != 0.) {
03112 goto L130;
03113 }
03114 t1 = hpi;
03115 goto L140;
03116 L130:
03117 t1 = atan(*zi / *zr);
03118 t1 = abs(t1);
03119 L140:
03120 if (t2 > caz) {
03121 goto L170;
03122 }
03123
03124
03125
03126 etest = ak / (dpi * caz * *tol);
03127 fk = coner;
03128 if (etest < coner) {
03129 goto L180;
03130 }
03131 fks = ctwor;
03132 ckr = caz + caz + ctwor;
03133 p1r = czeror;
03134 p2r = coner;
03135 i__1 = kmax;
03136 for (i__ = 1; i__ <= i__1; ++i__) {
03137 ak = fhs / fks;
03138 cbr = ckr / (fk + coner);
03139 ptr = p2r;
03140 p2r = cbr * p2r - p1r * ak;
03141 p1r = ptr;
03142 ckr += ctwor;
03143 fks = fks + fk + fk + ctwor;
03144 fhs = fhs + fk + fk;
03145 fk += coner;
03146 str = abs(p2r) * fk;
03147 if (etest < str) {
03148 goto L160;
03149 }
03150
03151 }
03152 goto L310;
03153 L160:
03154 fk += spi * t1 * sqrt(t2 / caz);
03155 fhs = (d__1 = .25 - dnu2, abs(d__1));
03156 goto L180;
03157 L170:
03158
03159
03160
03161 a2 = sqrt(caz);
03162 ak = fpi * ak / (*tol * sqrt(a2));
03163 aa = t1 * 3. / (caz + 1.);
03164 bb = t1 * 14.7 / (caz + 28.);
03165 ak = (log(ak) + caz * cos(aa) / (caz * .008 + 1.)) / cos(bb);
03166 fk = ak * .12125 * ak / caz + 1.5;
03167 L180:
03168
03169
03170
03171 k = (integer) ((real) fk);
03172 fk = (doublereal) ((real) k);
03173 fks = fk * fk;
03174 p1r = czeror;
03175 p1i = czeroi;
03176 p2r = *tol;
03177 p2i = czeroi;
03178 csr = p2r;
03179 csi = p2i;
03180 i__1 = k;
03181 for (i__ = 1; i__ <= i__1; ++i__) {
03182 a1 = fks - fk;
03183 ak = (fks + fk) / (a1 + fhs);
03184 rak = 2. / (fk + coner);
03185 cbr = (fk + *zr) * rak;
03186 cbi = *zi * rak;
03187 ptr = p2r;
03188 pti = p2i;
03189 p2r = (ptr * cbr - pti * cbi - p1r) * ak;
03190 p2i = (pti * cbr + ptr * cbi - p1i) * ak;
03191 p1r = ptr;
03192 p1i = pti;
03193 csr += p2r;
03194 csi += p2i;
03195 fks = a1 - fk + coner;
03196 fk -= coner;
03197
03198 }
03199
03200
03201
03202
03203 tm = myzabs_(&csr, &csi);
03204 ptr = 1. / tm;
03205 s1r = p2r * ptr;
03206 s1i = p2i * ptr;
03207 csr *= ptr;
03208 csi = -csi * ptr;
03209 zmlt_(&coefr, &coefi, &s1r, &s1i, &str, &sti);
03210 zmlt_(&str, &sti, &csr, &csi, &s1r, &s1i);
03211 if (inu > 0 || *n > 1) {
03212 goto L200;
03213 }
03214 zdr = *zr;
03215 zdi = *zi;
03216 if (iflag == 1) {
03217 goto L270;
03218 }
03219 goto L240;
03220 L200:
03221
03222
03223
03224 tm = myzabs_(&p2r, &p2i);
03225 ptr = 1. / tm;
03226 p1r *= ptr;
03227 p1i *= ptr;
03228 p2r *= ptr;
03229 p2i = -p2i * ptr;
03230 zmlt_(&p1r, &p1i, &p2r, &p2i, &ptr, &pti);
03231 str = dnu + .5 - ptr;
03232 sti = -pti;
03233 zdiv_(&str, &sti, zr, zi, &str, &sti);
03234 str += 1.;
03235 zmlt_(&str, &sti, &s1r, &s1i, &s2r, &s2i);
03236
03237
03238
03239
03240 L210:
03241 str = dnu + 1.;
03242 ckr = str * rzr;
03243 cki = str * rzi;
03244 if (*n == 1) {
03245 --inu;
03246 }
03247 if (inu > 0) {
03248 goto L220;
03249 }
03250 if (*n > 1) {
03251 goto L215;
03252 }
03253 s1r = s2r;
03254 s1i = s2i;
03255 L215:
03256 zdr = *zr;
03257 zdi = *zi;
03258 if (iflag == 1) {
03259 goto L270;
03260 }
03261 goto L240;
03262 L220:
03263 inub = 1;
03264 if (iflag == 1) {
03265 goto L261;
03266 }
03267 L225:
03268 p1r = csrr[kflag - 1];
03269 ascle = bry[kflag - 1];
03270 i__1 = inu;
03271 for (i__ = inub; i__ <= i__1; ++i__) {
03272 str = s2r;
03273 sti = s2i;
03274 s2r = ckr * str - cki * sti + s1r;
03275 s2i = ckr * sti + cki * str + s1i;
03276 s1r = str;
03277 s1i = sti;
03278 ckr += rzr;
03279 cki += rzi;
03280 if (kflag >= 3) {
03281 goto L230;
03282 }
03283 p2r = s2r * p1r;
03284 p2i = s2i * p1r;
03285 str = abs(p2r);
03286 sti = abs(p2i);
03287 p2m = max(str,sti);
03288 if (p2m <= ascle) {
03289 goto L230;
03290 }
03291 ++kflag;
03292 ascle = bry[kflag - 1];
03293 s1r *= p1r;
03294 s1i *= p1r;
03295 s2r = p2r;
03296 s2i = p2i;
03297 str = cssr[kflag - 1];
03298 s1r *= str;
03299 s1i *= str;
03300 s2r *= str;
03301 s2i *= str;
03302 p1r = csrr[kflag - 1];
03303 L230:
03304 ;
03305 }
03306 if (*n != 1) {
03307 goto L240;
03308 }
03309 s1r = s2r;
03310 s1i = s2i;
03311 L240:
03312 str = csrr[kflag - 1];
03313 yr[1] = s1r * str;
03314 yi[1] = s1i * str;
03315 if (*n == 1) {
03316 return 0;
03317 }
03318 yr[2] = s2r * str;
03319 yi[2] = s2i * str;
03320 if (*n == 2) {
03321 return 0;
03322 }
03323 kk = 2;
03324 L250:
03325 ++kk;
03326 if (kk > *n) {
03327 return 0;
03328 }
03329 p1r = csrr[kflag - 1];
03330 ascle = bry[kflag - 1];
03331 i__1 = *n;
03332 for (i__ = kk; i__ <= i__1; ++i__) {
03333 p2r = s2r;
03334 p2i = s2i;
03335 s2r = ckr * p2r - cki * p2i + s1r;
03336 s2i = cki * p2r + ckr * p2i + s1i;
03337 s1r = p2r;
03338 s1i = p2i;
03339 ckr += rzr;
03340 cki += rzi;
03341 p2r = s2r * p1r;
03342 p2i = s2i * p1r;
03343 yr[i__] = p2r;
03344 yi[i__] = p2i;
03345 if (kflag >= 3) {
03346 goto L260;
03347 }
03348 str = abs(p2r);
03349 sti = abs(p2i);
03350 p2m = max(str,sti);
03351 if (p2m <= ascle) {
03352 goto L260;
03353 }
03354 ++kflag;
03355 ascle = bry[kflag - 1];
03356 s1r *= p1r;
03357 s1i *= p1r;
03358 s2r = p2r;
03359 s2i = p2i;
03360 str = cssr[kflag - 1];
03361 s1r *= str;
03362 s1i *= str;
03363 s2r *= str;
03364 s2i *= str;
03365 p1r = csrr[kflag - 1];
03366 L260:
03367 ;
03368 }
03369 return 0;
03370
03371
03372
03373 L261:
03374 helim = *elim * .5;
03375 elm = exp(-(*elim));
03376 celmr = elm;
03377 ascle = bry[0];
03378 zdr = *zr;
03379 zdi = *zi;
03380 ic = -1;
03381 j = 2;
03382 i__1 = inu;
03383 for (i__ = 1; i__ <= i__1; ++i__) {
03384 str = s2r;
03385 sti = s2i;
03386 s2r = str * ckr - sti * cki + s1r;
03387 s2i = sti * ckr + str * cki + s1i;
03388 s1r = str;
03389 s1i = sti;
03390 ckr += rzr;
03391 cki += rzi;
03392 as = myzabs_(&s2r, &s2i);
03393 alas = log(as);
03394 p2r = -zdr + alas;
03395 if (p2r < -(*elim)) {
03396 goto L263;
03397 }
03398 zlog_(&s2r, &s2i, &str, &sti, &idum);
03399 p2r = -zdr + str;
03400 p2i = -zdi + sti;
03401 p2m = exp(p2r) / *tol;
03402 p1r = p2m * cos(p2i);
03403 p1i = p2m * sin(p2i);
03404 zuchk_(&p1r, &p1i, &nw, &ascle, tol);
03405 if (nw != 0) {
03406 goto L263;
03407 }
03408 j = 3 - j;
03409 cyr[j - 1] = p1r;
03410 cyi[j - 1] = p1i;
03411 if (ic == i__ - 1) {
03412 goto L264;
03413 }
03414 ic = i__;
03415 goto L262;
03416 L263:
03417 if (alas < helim) {
03418 goto L262;
03419 }
03420 zdr -= *elim;
03421 s1r *= celmr;
03422 s1i *= celmr;
03423 s2r *= celmr;
03424 s2i *= celmr;
03425 L262:
03426 ;
03427 }
03428 if (*n != 1) {
03429 goto L270;
03430 }
03431 s1r = s2r;
03432 s1i = s2i;
03433 goto L270;
03434 L264:
03435 kflag = 1;
03436 inub = i__ + 1;
03437 s2r = cyr[j - 1];
03438 s2i = cyi[j - 1];
03439 j = 3 - j;
03440 s1r = cyr[j - 1];
03441 s1i = cyi[j - 1];
03442 if (inub <= inu) {
03443 goto L225;
03444 }
03445 if (*n != 1) {
03446 goto L240;
03447 }
03448 s1r = s2r;
03449 s1i = s2i;
03450 goto L240;
03451 L270:
03452 yr[1] = s1r;
03453 yi[1] = s1i;
03454 if (*n == 1) {
03455 goto L280;
03456 }
03457 yr[2] = s2r;
03458 yi[2] = s2i;
03459 L280:
03460 ascle = bry[0];
03461 zkscl_(&zdr, &zdi, fnu, n, &yr[1], &yi[1], nz, &rzr, &rzi, &ascle, tol,
03462 elim);
03463 inu = *n - *nz;
03464 if (inu <= 0) {
03465 return 0;
03466 }
03467 kk = *nz + 1;
03468 s1r = yr[kk];
03469 s1i = yi[kk];
03470 yr[kk] = s1r * csrr[0];
03471 yi[kk] = s1i * csrr[0];
03472 if (inu == 1) {
03473 return 0;
03474 }
03475 kk = *nz + 2;
03476 s2r = yr[kk];
03477 s2i = yi[kk];
03478 yr[kk] = s2r * csrr[0];
03479 yi[kk] = s2i * csrr[0];
03480 if (inu == 2) {
03481 return 0;
03482 }
03483 t2 = *fnu + (doublereal) ((real) (kk - 1));
03484 ckr = t2 * rzr;
03485 cki = t2 * rzi;
03486 kflag = 1;
03487 goto L250;
03488 L290:
03489
03490
03491
03492 koded = 2;
03493 iflag = 1;
03494 kflag = 2;
03495 goto L120;
03496
03497
03498
03499 L300:
03500 s1r = coefr;
03501 s1i = coefi;
03502 s2r = coefr;
03503 s2i = coefi;
03504 goto L210;
03505
03506
03507 L310:
03508 *nz = -2;
03509 return 0;
03510 }
03511
03512 int zbuni_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *n, doublereal *yr, doublereal *yi, integer *nz, integer *nui, integer *nlast, doublereal *fnul, doublereal *tol, doublereal *elim, doublereal *alim)
03513 {
03514
03515 integer i__1;
03516
03517
03518 static doublereal dfnu, fnui;
03519 static integer i__, k, iflag;
03520 static doublereal ascle, csclr, cscrr;
03521 static integer iform;
03522 static doublereal ax, ay;
03523 static integer nl, nw;
03524 static doublereal c1i, c1m;
03525 static doublereal c1r, s1i, s2i, s1r, s2r, cyi[2], gnu, raz, cyr[2], sti,
03526 bry[3], rzi, str, rzr;
03527
03528
03529
03530
03531
03532
03533
03534
03535
03536
03537
03538
03539
03540
03541 --yi;
03542 --yr;
03543
03544
03545 *nz = 0;
03546 ax = abs(*zr) * 1.7321;
03547 ay = abs(*zi);
03548 iform = 1;
03549 if (ay > ax) {
03550 iform = 2;
03551 }
03552 if (*nui == 0) {
03553 goto L60;
03554 }
03555 fnui = (doublereal) ((real) (*nui));
03556 dfnu = *fnu + (doublereal) ((real) (*n - 1));
03557 gnu = dfnu + fnui;
03558 if (iform == 2) {
03559 goto L10;
03560 }
03561
03562
03563
03564
03565 zuni1_(zr, zi, &gnu, kode, &c__2, cyr, cyi, &nw, nlast, fnul, tol, elim,
03566 alim);
03567 goto L20;
03568 L10:
03569
03570
03571
03572
03573
03574 zuni2_(zr, zi, &gnu, kode, &c__2, cyr, cyi, &nw, nlast, fnul, tol, elim,
03575 alim);
03576 L20:
03577 if (nw < 0) {
03578 goto L50;
03579 }
03580 if (nw != 0) {
03581 goto L90;
03582 }
03583 str = myzabs_(cyr, cyi);
03584
03585
03586
03587 bry[0] = d1mach_(&c__1) * 1e3 / *tol;
03588 bry[1] = 1. / bry[0];
03589 bry[2] = bry[1];
03590 iflag = 2;
03591 ascle = bry[1];
03592 csclr = 1.;
03593 if (str > bry[0]) {
03594 goto L21;
03595 }
03596 iflag = 1;
03597 ascle = bry[0];
03598 csclr = 1. / *tol;
03599 goto L25;
03600 L21:
03601 if (str < bry[1]) {
03602 goto L25;
03603 }
03604 iflag = 3;
03605 ascle = bry[2];
03606 csclr = *tol;
03607 L25:
03608 cscrr = 1. / csclr;
03609 s1r = cyr[1] * csclr;
03610 s1i = cyi[1] * csclr;
03611 s2r = cyr[0] * csclr;
03612 s2i = cyi[0] * csclr;
03613 raz = 1. / myzabs_(zr, zi);
03614 str = *zr * raz;
03615 sti = -(*zi) * raz;
03616 rzr = (str + str) * raz;
03617 rzi = (sti + sti) * raz;
03618 i__1 = *nui;
03619 for (i__ = 1; i__ <= i__1; ++i__) {
03620 str = s2r;
03621 sti = s2i;
03622 s2r = (dfnu + fnui) * (rzr * str - rzi * sti) + s1r;
03623 s2i = (dfnu + fnui) * (rzr * sti + rzi * str) + s1i;
03624 s1r = str;
03625 s1i = sti;
03626 fnui += -1.;
03627 if (iflag >= 3) {
03628 goto L30;
03629 }
03630 str = s2r * cscrr;
03631 sti = s2i * cscrr;
03632 c1r = abs(str);
03633 c1i = abs(sti);
03634 c1m = max(c1r,c1i);
03635 if (c1m <= ascle) {
03636 goto L30;
03637 }
03638 ++iflag;
03639 ascle = bry[iflag - 1];
03640 s1r *= cscrr;
03641 s1i *= cscrr;
03642 s2r = str;
03643 s2i = sti;
03644 csclr *= *tol;
03645 cscrr = 1. / csclr;
03646 s1r *= csclr;
03647 s1i *= csclr;
03648 s2r *= csclr;
03649 s2i *= csclr;
03650 L30:
03651 ;
03652 }
03653 yr[*n] = s2r * cscrr;
03654 yi[*n] = s2i * cscrr;
03655 if (*n == 1) {
03656 return 0;
03657 }
03658 nl = *n - 1;
03659 fnui = (doublereal) ((real) nl);
03660 k = nl;
03661 i__1 = nl;
03662 for (i__ = 1; i__ <= i__1; ++i__) {
03663 str = s2r;
03664 sti = s2i;
03665 s2r = (*fnu + fnui) * (rzr * str - rzi * sti) + s1r;
03666 s2i = (*fnu + fnui) * (rzr * sti + rzi * str) + s1i;
03667 s1r = str;
03668 s1i = sti;
03669 str = s2r * cscrr;
03670 sti = s2i * cscrr;
03671 yr[k] = str;
03672 yi[k] = sti;
03673 fnui += -1.;
03674 --k;
03675 if (iflag >= 3) {
03676 goto L40;
03677 }
03678 c1r = abs(str);
03679 c1i = abs(sti);
03680 c1m = max(c1r,c1i);
03681 if (c1m <= ascle) {
03682 goto L40;
03683 }
03684 ++iflag;
03685 ascle = bry[iflag - 1];
03686 s1r *= cscrr;
03687 s1i *= cscrr;
03688 s2r = str;
03689 s2i = sti;
03690 csclr *= *tol;
03691 cscrr = 1. / csclr;
03692 s1r *= csclr;
03693 s1i *= csclr;
03694 s2r *= csclr;
03695 s2i *= csclr;
03696 L40:
03697 ;
03698 }
03699 return 0;
03700 L50:
03701 *nz = -1;
03702 if (nw == -2) {
03703 *nz = -2;
03704 }
03705 return 0;
03706 L60:
03707 if (iform == 2) {
03708 goto L70;
03709 }
03710
03711
03712
03713
03714 zuni1_(zr, zi, fnu, kode, n, &yr[1], &yi[1], &nw, nlast, fnul, tol, elim,
03715 alim);
03716 goto L80;
03717 L70:
03718
03719
03720
03721
03722
03723 zuni2_(zr, zi, fnu, kode, n, &yr[1], &yi[1], &nw, nlast, fnul, tol, elim,
03724 alim);
03725 L80:
03726 if (nw < 0) {
03727 goto L50;
03728 }
03729 *nz = nw;
03730 return 0;
03731 L90:
03732 *nlast = *n;
03733 return 0;
03734 }
03735
03736 int zbunk_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *mr, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *tol, doublereal *elim, doublereal *alim)
03737 {
03738 static doublereal ax, ay;
03739
03740
03741
03742
03743
03744
03745
03746
03747
03748
03749
03750
03751 --yi;
03752 --yr;
03753
03754
03755 *nz = 0;
03756 ax = abs(*zr) * 1.7321;
03757 ay = abs(*zi);
03758 if (ay > ax) {
03759 goto L10;
03760 }
03761
03762
03763
03764
03765 zunk1_(zr, zi, fnu, kode, mr, n, &yr[1], &yi[1], nz, tol, elim, alim);
03766 goto L20;
03767 L10:
03768
03769
03770
03771
03772
03773 zunk2_(zr, zi, fnu, kode, mr, n, &yr[1], &yi[1], nz, tol, elim, alim);
03774 L20:
03775 return 0;
03776 }
03777
03778 int zdiv_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, doublereal *cr, doublereal *ci)
03779 {
03780 static doublereal ca, cb, cc, cd, bm;
03781
03782
03783
03784
03785
03786
03787
03788
03789 bm = 1. / myzabs_(br, bi);
03790 cc = *br * bm;
03791 cd = *bi * bm;
03792 ca = (*ar * cc + *ai * cd) * bm;
03793 cb = (*ai * cc - *ar * cd) * bm;
03794 *cr = ca;
03795 *ci = cb;
03796 return 0;
03797 }
03798
03799 int zexp_(const doublereal *ar, const doublereal *ai, doublereal *br, doublereal *bi)
03800 {
03801
03802
03803
03804 static doublereal ca, cb, zm;
03805
03806
03807
03808
03809
03810
03811
03812
03813 zm = exp(*ar);
03814 ca = zm * cos(*ai);
03815 cb = zm * sin(*ai);
03816 *br = ca;
03817 *bi = cb;
03818 return 0;
03819 }
03820
03821 int zkscl_(doublereal *zrr, doublereal *zri, doublereal *fnu, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *rzr, doublereal *rzi, doublereal *ascle, doublereal *tol, doublereal *elim)
03822 {
03823
03824
03825 static doublereal zeror = 0.;
03826 static doublereal zeroi = 0.;
03827
03828
03829 integer i__1;
03830
03831
03832
03833
03834 static doublereal alas;
03835 static integer idum;
03836 static integer i__;
03837 static doublereal helim, celmr;
03838 static integer ic;
03839 static doublereal as, fn;
03840 static integer kk, nn, nw;
03841 static doublereal s1i, s2i, s1r, s2r, acs, cki, elm, csi, ckr, cyi[2],
03842 zdi, csr, cyr[2], zdr, str;
03843
03844
03845
03846
03847
03848
03849
03850
03851
03852
03853
03854
03855 --yi;
03856 --yr;
03857
03858
03859
03860 *nz = 0;
03861 ic = 0;
03862 nn = min(2,*n);
03863 i__1 = nn;
03864 for (i__ = 1; i__ <= i__1; ++i__) {
03865 s1r = yr[i__];
03866 s1i = yi[i__];
03867 cyr[i__ - 1] = s1r;
03868 cyi[i__ - 1] = s1i;
03869 as = myzabs_(&s1r, &s1i);
03870 acs = -(*zrr) + log(as);
03871 ++(*nz);
03872 yr[i__] = zeror;
03873 yi[i__] = zeroi;
03874 if (acs < -(*elim)) {
03875 goto L10;
03876 }
03877 zlog_(&s1r, &s1i, &csr, &csi, &idum);
03878 csr -= *zrr;
03879 csi -= *zri;
03880 str = exp(csr) / *tol;
03881 csr = str * cos(csi);
03882 csi = str * sin(csi);
03883 zuchk_(&csr, &csi, &nw, ascle, tol);
03884 if (nw != 0) {
03885 goto L10;
03886 }
03887 yr[i__] = csr;
03888 yi[i__] = csi;
03889 ic = i__;
03890 --(*nz);
03891 L10:
03892 ;
03893 }
03894 if (*n == 1) {
03895 return 0;
03896 }
03897 if (ic > 1) {
03898 goto L20;
03899 }
03900 yr[1] = zeror;
03901 yi[1] = zeroi;
03902 *nz = 2;
03903 L20:
03904 if (*n == 2) {
03905 return 0;
03906 }
03907 if (*nz == 0) {
03908 return 0;
03909 }
03910 fn = *fnu + 1.;
03911 ckr = fn * *rzr;
03912 cki = fn * *rzi;
03913 s1r = cyr[0];
03914 s1i = cyi[0];
03915 s2r = cyr[1];
03916 s2i = cyi[1];
03917 helim = *elim * .5;
03918 elm = exp(-(*elim));
03919 celmr = elm;
03920 zdr = *zrr;
03921 zdi = *zri;
03922
03923
03924
03925
03926 i__1 = *n;
03927 for (i__ = 3; i__ <= i__1; ++i__) {
03928 kk = i__;
03929 csr = s2r;
03930 csi = s2i;
03931 s2r = ckr * csr - cki * csi + s1r;
03932 s2i = cki * csr + ckr * csi + s1i;
03933 s1r = csr;
03934 s1i = csi;
03935 ckr += *rzr;
03936 cki += *rzi;
03937 as = myzabs_(&s2r, &s2i);
03938 alas = log(as);
03939 acs = -zdr + alas;
03940 ++(*nz);
03941 yr[i__] = zeror;
03942 yi[i__] = zeroi;
03943 if (acs < -(*elim)) {
03944 goto L25;
03945 }
03946 zlog_(&s2r, &s2i, &csr, &csi, &idum);
03947 csr -= zdr;
03948 csi -= zdi;
03949 str = exp(csr) / *tol;
03950 csr = str * cos(csi);
03951 csi = str * sin(csi);
03952 zuchk_(&csr, &csi, &nw, ascle, tol);
03953 if (nw != 0) {
03954 goto L25;
03955 }
03956 yr[i__] = csr;
03957 yi[i__] = csi;
03958 --(*nz);
03959 if (ic == kk - 1) {
03960 goto L40;
03961 }
03962 ic = kk;
03963 goto L30;
03964 L25:
03965 if (alas < helim) {
03966 goto L30;
03967 }
03968 zdr -= *elim;
03969 s1r *= celmr;
03970 s1i *= celmr;
03971 s2r *= celmr;
03972 s2i *= celmr;
03973 L30:
03974 ;
03975 }
03976 *nz = *n;
03977 if (ic == *n) {
03978 *nz = *n - 1;
03979 }
03980 goto L45;
03981 L40:
03982 *nz = kk - 2;
03983 L45:
03984 i__1 = *nz;
03985 for (i__ = 1; i__ <= i__1; ++i__) {
03986 yr[i__] = zeror;
03987 yi[i__] = zeroi;
03988
03989 }
03990 return 0;
03991 }
03992
03993 int zlog_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, integer *ierr)
03994 {
03995
03996
03997 static doublereal dpi = 3.141592653589793238462643383;
03998 static doublereal dhpi = 1.570796326794896619231321696;
03999
04000
04001
04002
04003 static doublereal zm, dtheta;
04004
04005
04006
04007
04008
04009
04010
04011
04012
04013 *ierr = 0;
04014 if (*ar == 0.) {
04015 goto L10;
04016 }
04017 if (*ai == 0.) {
04018 goto L20;
04019 }
04020 dtheta = atan(*ai / *ar);
04021 if (dtheta <= 0.) {
04022 goto L40;
04023 }
04024 if (*ar < 0.) {
04025 dtheta -= dpi;
04026 }
04027 goto L50;
04028 L10:
04029 if (*ai == 0.) {
04030 goto L60;
04031 }
04032 *bi = dhpi;
04033 *br = log((abs(*ai)));
04034 if (*ai < 0.) {
04035 *bi = -(*bi);
04036 }
04037 return 0;
04038 L20:
04039 if (*ar > 0.) {
04040 goto L30;
04041 }
04042 *br = log((abs(*ar)));
04043 *bi = dpi;
04044 return 0;
04045 L30:
04046 *br = log(*ar);
04047 *bi = 0.;
04048 return 0;
04049 L40:
04050 if (*ar < 0.) {
04051 dtheta += dpi;
04052 }
04053 L50:
04054 zm = myzabs_(ar, ai);
04055 *br = log(zm);
04056 *bi = dtheta;
04057 return 0;
04058 L60:
04059 *ierr = 1;
04060 return 0;
04061 }
04062
04063
04064 int zmlri_(doublereal *zr, doublereal *zi, doublereal *fnu,
04065 integer *kode, integer *n,
04066 doublereal *yr, doublereal *yi, integer *nz, doublereal *tol)
04067 {
04068
04069
04070 static doublereal zeror = 0.;
04071 static doublereal zeroi = 0.;
04072 static doublereal coner = 1.;
04073 static doublereal conei = 0.;
04074
04075
04076 integer i__1, i__2;
04077 doublereal d__1, d__2, d__3;
04078
04079
04080
04081
04082 static doublereal flam, fkap, scle, tfnf;
04083 static integer idum, ifnu;
04084 static doublereal sumi, sumr;
04085 static integer i__, k, m, itime;
04086 static doublereal ak, bk, ap, at;
04087 static integer kk, km;
04088 static doublereal az;
04089 static doublereal cnormi, cnormr;
04090 static doublereal p1i, p2i, p1r, p2r, ack, cki, fnf, fkk, ckr;
04091 static integer iaz;
04092 static doublereal rho;
04093 static integer inu;
04094 static doublereal pti, raz, sti, rzi, ptr, str, tst, rzr, rho2;
04095
04096
04097
04098
04099
04100
04101
04102
04103
04104
04105
04106 --yi;
04107 --yr;
04108
04109
04110 scle = d1mach_(&c__1) / *tol;
04111 *nz = 0;
04112 az = myzabs_(zr, zi);
04113 iaz = (integer) ((real) az);
04114 ifnu = (integer) ((real) (*fnu));
04115 inu = ifnu + *n - 1;
04116 at = (doublereal) ((real) iaz) + 1.;
04117 raz = 1. / az;
04118 str = *zr * raz;
04119 sti = -(*zi) * raz;
04120 ckr = str * at * raz;
04121 cki = sti * at * raz;
04122 rzr = (str + str) * raz;
04123 rzi = (sti + sti) * raz;
04124 p1r = zeror;
04125 p1i = zeroi;
04126 p2r = coner;
04127 p2i = conei;
04128 ack = (at + 1.) * raz;
04129 rho = ack + sqrt(ack * ack - 1.);
04130 rho2 = rho * rho;
04131 tst = (rho2 + rho2) / ((rho2 - 1.) * (rho - 1.));
04132 tst /= *tol;
04133
04134
04135
04136 ak = at;
04137 for (i__ = 1; i__ <= 80; ++i__) {
04138 ptr = p2r;
04139 pti = p2i;
04140 p2r = p1r - (ckr * ptr - cki * pti);
04141 p2i = p1i - (cki * ptr + ckr * pti);
04142 p1r = ptr;
04143 p1i = pti;
04144 ckr += rzr;
04145 cki += rzi;
04146 ap = myzabs_(&p2r, &p2i);
04147 if (ap > tst * ak * ak) {
04148 goto L20;
04149 }
04150 ak += 1.;
04151
04152 }
04153 goto L110;
04154 L20:
04155 ++i__;
04156 k = 0;
04157 if (inu < iaz) {
04158 goto L40;
04159 }
04160
04161
04162
04163 p1r = zeror;
04164 p1i = zeroi;
04165 p2r = coner;
04166 p2i = conei;
04167 at = (doublereal) ((real) inu) + 1.;
04168 str = *zr * raz;
04169 sti = -(*zi) * raz;
04170 ckr = str * at * raz;
04171 cki = sti * at * raz;
04172 ack = at * raz;
04173 tst = sqrt(ack / *tol);
04174 itime = 1;
04175 for (k = 1; k <= 80; ++k) {
04176 ptr = p2r;
04177 pti = p2i;
04178 p2r = p1r - (ckr * ptr - cki * pti);
04179 p2i = p1i - (ckr * pti + cki * ptr);
04180 p1r = ptr;
04181 p1i = pti;
04182 ckr += rzr;
04183 cki += rzi;
04184 ap = myzabs_(&p2r, &p2i);
04185 if (ap < tst) {
04186 goto L30;
04187 }
04188 if (itime == 2) {
04189 goto L40;
04190 }
04191 ack = myzabs_(&ckr, &cki);
04192 flam = ack + sqrt(ack * ack - 1.);
04193 fkap = ap / myzabs_(&p1r, &p1i);
04194 rho = min(flam,fkap);
04195 tst *= sqrt(rho / (rho * rho - 1.));
04196 itime = 2;
04197 L30:
04198 ;
04199 }
04200 goto L110;
04201 L40:
04202
04203
04204
04205 ++k;
04206
04207 i__1 = i__ + iaz, i__2 = k + inu;
04208 kk = max(i__1,i__2);
04209 fkk = (doublereal) ((real) kk);
04210 p1r = zeror;
04211 p1i = zeroi;
04212
04213
04214
04215 p2r = scle;
04216 p2i = zeroi;
04217 fnf = *fnu - (doublereal) ((real) ifnu);
04218 tfnf = fnf + fnf;
04219 d__1 = fkk + tfnf + 1.;
04220 d__2 = fkk + 1.;
04221 d__3 = tfnf + 1.;
04222 bk = dgamln_(&d__1, &idum) - dgamln_(&d__2, &idum) - dgamln_(&d__3, &idum)
04223 ;
04224 bk = exp(bk);
04225 sumr = zeror;
04226 sumi = zeroi;
04227 km = kk - inu;
04228 i__1 = km;
04229 for (i__ = 1; i__ <= i__1; ++i__) {
04230 ptr = p2r;
04231 pti = p2i;
04232 p2r = p1r + (fkk + fnf) * (rzr * ptr - rzi * pti);
04233 p2i = p1i + (fkk + fnf) * (rzi * ptr + rzr * pti);
04234 p1r = ptr;
04235 p1i = pti;
04236 ak = 1. - tfnf / (fkk + tfnf);
04237 ack = bk * ak;
04238 sumr += (ack + bk) * p1r;
04239 sumi += (ack + bk) * p1i;
04240 bk = ack;
04241 fkk += -1.;
04242
04243 }
04244 yr[*n] = p2r;
04245 yi[*n] = p2i;
04246 if (*n == 1) {
04247 goto L70;
04248 }
04249 i__1 = *n;
04250 for (i__ = 2; i__ <= i__1; ++i__) {
04251 ptr = p2r;
04252 pti = p2i;
04253 p2r = p1r + (fkk + fnf) * (rzr * ptr - rzi * pti);
04254 p2i = p1i + (fkk + fnf) * (rzi * ptr + rzr * pti);
04255 p1r = ptr;
04256 p1i = pti;
04257 ak = 1. - tfnf / (fkk + tfnf);
04258 ack = bk * ak;
04259 sumr += (ack + bk) * p1r;
04260 sumi += (ack + bk) * p1i;
04261 bk = ack;
04262 fkk += -1.;
04263 m = *n - i__ + 1;
04264 yr[m] = p2r;
04265 yi[m] = p2i;
04266
04267 }
04268 L70:
04269 if (ifnu <= 0) {
04270 goto L90;
04271 }
04272 i__1 = ifnu;
04273 for (i__ = 1; i__ <= i__1; ++i__) {
04274 ptr = p2r;
04275 pti = p2i;
04276 p2r = p1r + (fkk + fnf) * (rzr * ptr - rzi * pti);
04277 p2i = p1i + (fkk + fnf) * (rzr * pti + rzi * ptr);
04278 p1r = ptr;
04279 p1i = pti;
04280 ak = 1. - tfnf / (fkk + tfnf);
04281 ack = bk * ak;
04282 sumr += (ack + bk) * p1r;
04283 sumi += (ack + bk) * p1i;
04284 bk = ack;
04285 fkk += -1.;
04286
04287 }
04288 L90:
04289 ptr = *zr;
04290 pti = *zi;
04291 if (*kode == 2) {
04292 ptr = zeror;
04293 }
04294 zlog_(&rzr, &rzi, &str, &sti, &idum);
04295 p1r = -fnf * str + ptr;
04296 p1i = -fnf * sti + pti;
04297 d__1 = fnf + 1.;
04298 ap = dgamln_(&d__1, &idum);
04299 ptr = p1r - ap;
04300 pti = p1i;
04301
04302
04303
04304
04305 p2r += sumr;
04306 p2i += sumi;
04307 ap = myzabs_(&p2r, &p2i);
04308 p1r = 1. / ap;
04309 zexp_(&ptr, &pti, &str, &sti);
04310 ckr = str * p1r;
04311 cki = sti * p1r;
04312 ptr = p2r * p1r;
04313 pti = -p2i * p1r;
04314 zmlt_(&ckr, &cki, &ptr, &pti, &cnormr, &cnormi);
04315 i__1 = *n;
04316 for (i__ = 1; i__ <= i__1; ++i__) {
04317 str = yr[i__] * cnormr - yi[i__] * cnormi;
04318 yi[i__] = yr[i__] * cnormi + yi[i__] * cnormr;
04319 yr[i__] = str;
04320
04321 }
04322 return 0;
04323 L110:
04324 *nz = -2;
04325 return 0;
04326 }
04327
04328 int zmlt_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, doublereal *cr, doublereal *ci)
04329 {
04330 static doublereal ca, cb;
04331
04332
04333
04334
04335
04336
04337
04338
04339 ca = *ar * *br - *ai * *bi;
04340 cb = *ar * *bi + *ai * *br;
04341 *cr = ca;
04342 *ci = cb;
04343 return 0;
04344 }
04345
04346 int zrati_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *n, doublereal *cyr, doublereal *cyi, doublereal *tol)
04347 {
04348
04349
04350 static doublereal czeror = 0.;
04351 static doublereal czeroi = 0.;
04352 static doublereal coner = 1.;
04353 static doublereal conei = 0.;
04354 static doublereal rt2 = 1.41421356237309505;
04355
04356
04357 integer i__1;
04358 doublereal d__1;
04359
04360
04361
04362
04363 static doublereal flam, dfnu, fdnu;
04364 static integer magz, idnu;
04365 static doublereal fnup;
04366 static doublereal test, test1;
04367 static integer i__, k;
04368 static doublereal amagz;
04369 static integer itime;
04370 static doublereal ak;
04371 static integer id, kk;
04372 static doublereal az, cdfnui, cdfnur, ap1, ap2;
04373 static doublereal p1i, p2i, t1i, p1r, p2r, t1r, arg, rak, rho;
04374 static integer inu;
04375 static doublereal pti, tti, rzi, ptr, ttr, rzr, rap1;
04376
04377
04378
04379
04380
04381
04382
04383
04384
04385
04386
04387
04388
04389
04390
04391 --cyi;
04392 --cyr;
04393
04394
04395 az = myzabs_(zr, zi);
04396 inu = (integer) ((real) (*fnu));
04397 idnu = inu + *n - 1;
04398 magz = (integer) ((real) az);
04399 amagz = (doublereal) ((real) (magz + 1));
04400 fdnu = (doublereal) ((real) idnu);
04401 fnup = max(amagz,fdnu);
04402 id = idnu - magz - 1;
04403 itime = 1;
04404 k = 1;
04405 ptr = 1. / az;
04406 rzr = ptr * (*zr + *zr) * ptr;
04407 rzi = -ptr * (*zi + *zi) * ptr;
04408 t1r = rzr * fnup;
04409 t1i = rzi * fnup;
04410 p2r = -t1r;
04411 p2i = -t1i;
04412 p1r = coner;
04413 p1i = conei;
04414 t1r += rzr;
04415 t1i += rzi;
04416 if (id > 0) {
04417 id = 0;
04418 }
04419 ap2 = myzabs_(&p2r, &p2i);
04420 ap1 = myzabs_(&p1r, &p1i);
04421
04422
04423
04424
04425
04426
04427 arg = (ap2 + ap2) / (ap1 * *tol);
04428 test1 = sqrt(arg);
04429 test = test1;
04430 rap1 = 1. / ap1;
04431 p1r *= rap1;
04432 p1i *= rap1;
04433 p2r *= rap1;
04434 p2i *= rap1;
04435 ap2 *= rap1;
04436 L10:
04437 ++k;
04438 ap1 = ap2;
04439 ptr = p2r;
04440 pti = p2i;
04441 p2r = p1r - (t1r * ptr - t1i * pti);
04442 p2i = p1i - (t1r * pti + t1i * ptr);
04443 p1r = ptr;
04444 p1i = pti;
04445 t1r += rzr;
04446 t1i += rzi;
04447 ap2 = myzabs_(&p2r, &p2i);
04448 if (ap1 <= test) {
04449 goto L10;
04450 }
04451 if (itime == 2) {
04452 goto L20;
04453 }
04454 ak = myzabs_(&t1r, &t1i) * .5;
04455 flam = ak + sqrt(ak * ak - 1.);
04456
04457 d__1 = ap2 / ap1;
04458 rho = min(d__1,flam);
04459 test = test1 * sqrt(rho / (rho * rho - 1.));
04460 itime = 2;
04461 goto L10;
04462 L20:
04463 kk = k + 1 - id;
04464 ak = (doublereal) ((real) kk);
04465 t1r = ak;
04466 t1i = czeroi;
04467 dfnu = *fnu + (doublereal) ((real) (*n - 1));
04468 p1r = 1. / ap2;
04469 p1i = czeroi;
04470 p2r = czeror;
04471 p2i = czeroi;
04472 i__1 = kk;
04473 for (i__ = 1; i__ <= i__1; ++i__) {
04474 ptr = p1r;
04475 pti = p1i;
04476 rap1 = dfnu + t1r;
04477 ttr = rzr * rap1;
04478 tti = rzi * rap1;
04479 p1r = ptr * ttr - pti * tti + p2r;
04480 p1i = ptr * tti + pti * ttr + p2i;
04481 p2r = ptr;
04482 p2i = pti;
04483 t1r -= coner;
04484
04485 }
04486 if (p1r != czeror || p1i != czeroi) {
04487 goto L40;
04488 }
04489 p1r = *tol;
04490 p1i = *tol;
04491 L40:
04492 zdiv_(&p2r, &p2i, &p1r, &p1i, &cyr[*n], &cyi[*n]);
04493 if (*n == 1) {
04494 return 0;
04495 }
04496 k = *n - 1;
04497 ak = (doublereal) ((real) k);
04498 t1r = ak;
04499 t1i = czeroi;
04500 cdfnur = *fnu * rzr;
04501 cdfnui = *fnu * rzi;
04502 i__1 = *n;
04503 for (i__ = 2; i__ <= i__1; ++i__) {
04504 ptr = cdfnur + (t1r * rzr - t1i * rzi) + cyr[k + 1];
04505 pti = cdfnui + (t1r * rzi + t1i * rzr) + cyi[k + 1];
04506 ak = myzabs_(&ptr, &pti);
04507 if (ak != czeror) {
04508 goto L50;
04509 }
04510 ptr = *tol;
04511 pti = *tol;
04512 ak = *tol * rt2;
04513 L50:
04514 rak = coner / ak;
04515 cyr[k] = rak * ptr * rak;
04516 cyi[k] = -rak * pti * rak;
04517 t1r -= coner;
04518 --k;
04519
04520 }
04521 return 0;
04522 }
04523
04524 int zs1s2_(doublereal *zrr, doublereal *zri, doublereal *s1r, doublereal *s1i, doublereal *s2r, doublereal *s2i, integer *nz, doublereal *ascle, doublereal *alim, integer *iuf)
04525 {
04526
04527
04528 static doublereal zeror = 0.;
04529 static doublereal zeroi = 0.;
04530
04531
04532
04533
04534 static integer idum;
04535 static doublereal aa, c1i, as1, as2, c1r;
04536 static doublereal aln, s1di, s1dr;
04537
04538
04539
04540
04541
04542
04543
04544
04545
04546
04547
04548
04549
04550
04551
04552 *nz = 0;
04553 as1 = myzabs_(s1r, s1i);
04554 as2 = myzabs_(s2r, s2i);
04555 if (*s1r == 0. && *s1i == 0.) {
04556 goto L10;
04557 }
04558 if (as1 == 0.) {
04559 goto L10;
04560 }
04561 aln = -(*zrr) - *zrr + log(as1);
04562 s1dr = *s1r;
04563 s1di = *s1i;
04564 *s1r = zeror;
04565 *s1i = zeroi;
04566 as1 = zeror;
04567 if (aln < -(*alim)) {
04568 goto L10;
04569 }
04570 zlog_(&s1dr, &s1di, &c1r, &c1i, &idum);
04571 c1r = c1r - *zrr - *zrr;
04572 c1i = c1i - *zri - *zri;
04573 zexp_(&c1r, &c1i, s1r, s1i);
04574 as1 = myzabs_(s1r, s1i);
04575 ++(*iuf);
04576 L10:
04577 aa = max(as1,as2);
04578 if (aa > *ascle) {
04579 return 0;
04580 }
04581 *s1r = zeror;
04582 *s1i = zeroi;
04583 *s2r = zeror;
04584 *s2i = zeroi;
04585 *nz = 1;
04586 *iuf = 0;
04587 return 0;
04588 }
04589
04590
04591 int zseri_(doublereal *zr, doublereal *zi, doublereal *fnu,
04592 integer *kode, integer *n, doublereal *yr, doublereal *yi,
04593 integer *nz,
04594 doublereal *tol, doublereal *elim, doublereal *alim)
04595 {
04596
04597
04598 static doublereal zeror = 0.;
04599 static doublereal zeroi = 0.;
04600 static doublereal coner = 1.;
04601 static doublereal conei = 0.;
04602
04603
04604 integer i__1;
04605
04606
04607
04608
04609 static doublereal dfnu;
04610 static integer idum;
04611 static doublereal atol, fnup;
04612 static integer i__, k, l, m, iflag;
04613 static doublereal s, coefi, ascle, coefr, crscr;
04614 static doublereal aa;
04615 static integer ib;
04616 static doublereal ak;
04617 static integer il;
04618 static doublereal az;
04619 static integer nn;
04620 static doublereal wi[2];
04621 static doublereal rs, ss;
04622 static integer nw;
04623 static doublereal wr[2];
04624 static doublereal s1i, s2i, s1r, s2r, cki, acz, arm, ckr, czi, hzi, raz,
04625 czr, sti, hzr, rzi, str, rzr, ak1i, ak1r, rtr1;
04626
04627
04628
04629
04630
04631
04632
04633
04634
04635
04636
04637
04638
04639
04640
04641
04642 --yi;
04643 --yr;
04644
04645
04646
04647 *nz = 0;
04648 az = myzabs_(zr, zi);
04649 if (az == 0.) {
04650 goto L160;
04651 }
04652 arm = d1mach_(&c__1) * 1e3;
04653 rtr1 = sqrt(arm);
04654 crscr = 1.;
04655 iflag = 0;
04656 if (az < arm) {
04657 goto L150;
04658 }
04659 hzr = *zr * .5;
04660 hzi = *zi * .5;
04661 czr = zeror;
04662 czi = zeroi;
04663 if (az <= rtr1) {
04664 goto L10;
04665 }
04666 zmlt_(&hzr, &hzi, &hzr, &hzi, &czr, &czi);
04667 L10:
04668 acz = myzabs_(&czr, &czi);
04669 nn = *n;
04670 zlog_(&hzr, &hzi, &ckr, &cki, &idum);
04671 L20:
04672 dfnu = *fnu + (doublereal) ((real) (nn - 1));
04673 fnup = dfnu + 1.;
04674
04675
04676
04677 ak1r = ckr * dfnu;
04678 ak1i = cki * dfnu;
04679 ak = dgamln_(&fnup, &idum);
04680 ak1r -= ak;
04681 if (*kode == 2) {
04682 ak1r -= *zr;
04683 }
04684 if (ak1r > -(*elim)) {
04685 goto L40;
04686 }
04687 L30:
04688 ++(*nz);
04689 yr[nn] = zeror;
04690 yi[nn] = zeroi;
04691 if (acz > dfnu) {
04692 goto L190;
04693 }
04694 --nn;
04695 if (nn == 0) {
04696 return 0;
04697 }
04698 goto L20;
04699 L40:
04700 if (ak1r > -(*alim)) {
04701 goto L50;
04702 }
04703 iflag = 1;
04704 ss = 1. / *tol;
04705 crscr = *tol;
04706 ascle = arm * ss;
04707 L50:
04708 aa = exp(ak1r);
04709 if (iflag == 1) {
04710 aa *= ss;
04711 }
04712 coefr = aa * cos(ak1i);
04713 coefi = aa * sin(ak1i);
04714 atol = *tol * acz / fnup;
04715 il = min(2,nn);
04716 i__1 = il;
04717 for (i__ = 1; i__ <= i__1; ++i__) {
04718 dfnu = *fnu + (doublereal) ((real) (nn - i__));
04719 fnup = dfnu + 1.;
04720 s1r = coner;
04721 s1i = conei;
04722 if (acz < *tol * fnup) {
04723 goto L70;
04724 }
04725 ak1r = coner;
04726 ak1i = conei;
04727 ak = fnup + 2.;
04728 s = fnup;
04729 aa = 2.;
04730 L60:
04731 rs = 1. / s;
04732 str = ak1r * czr - ak1i * czi;
04733 sti = ak1r * czi + ak1i * czr;
04734 ak1r = str * rs;
04735 ak1i = sti * rs;
04736 s1r += ak1r;
04737 s1i += ak1i;
04738 s += ak;
04739 ak += 2.;
04740 aa = aa * acz * rs;
04741 if (aa > atol) {
04742 goto L60;
04743 }
04744 L70:
04745 s2r = s1r * coefr - s1i * coefi;
04746 s2i = s1r * coefi + s1i * coefr;
04747 wr[i__ - 1] = s2r;
04748 wi[i__ - 1] = s2i;
04749 if (iflag == 0) {
04750 goto L80;
04751 }
04752 zuchk_(&s2r, &s2i, &nw, &ascle, tol);
04753 if (nw != 0) {
04754 goto L30;
04755 }
04756 L80:
04757 m = nn - i__ + 1;
04758 yr[m] = s2r * crscr;
04759 yi[m] = s2i * crscr;
04760 if (i__ == il) {
04761 goto L90;
04762 }
04763 zdiv_(&coefr, &coefi, &hzr, &hzi, &str, &sti);
04764 coefr = str * dfnu;
04765 coefi = sti * dfnu;
04766 L90:
04767 ;
04768 }
04769 if (nn <= 2) {
04770 return 0;
04771 }
04772 k = nn - 2;
04773 ak = (doublereal) ((real) k);
04774 raz = 1. / az;
04775 str = *zr * raz;
04776 sti = -(*zi) * raz;
04777 rzr = (str + str) * raz;
04778 rzi = (sti + sti) * raz;
04779 if (iflag == 1) {
04780 goto L120;
04781 }
04782 ib = 3;
04783 L100:
04784 i__1 = nn;
04785 for (i__ = ib; i__ <= i__1; ++i__) {
04786 yr[k] = (ak + *fnu) * (rzr * yr[k + 1] - rzi * yi[k + 1]) + yr[k + 2];
04787 yi[k] = (ak + *fnu) * (rzr * yi[k + 1] + rzi * yr[k + 1]) + yi[k + 2];
04788 ak += -1.;
04789 --k;
04790
04791 }
04792 return 0;
04793
04794
04795
04796 L120:
04797
04798
04799
04800
04801 s1r = wr[0];
04802 s1i = wi[0];
04803 s2r = wr[1];
04804 s2i = wi[1];
04805 i__1 = nn;
04806 for (l = 3; l <= i__1; ++l) {
04807 ckr = s2r;
04808 cki = s2i;
04809 s2r = s1r + (ak + *fnu) * (rzr * ckr - rzi * cki);
04810 s2i = s1i + (ak + *fnu) * (rzr * cki + rzi * ckr);
04811 s1r = ckr;
04812 s1i = cki;
04813 ckr = s2r * crscr;
04814 cki = s2i * crscr;
04815 yr[k] = ckr;
04816 yi[k] = cki;
04817 ak += -1.;
04818 --k;
04819 if (myzabs_(&ckr, &cki) > ascle) {
04820 goto L140;
04821 }
04822
04823 }
04824 return 0;
04825 L140:
04826 ib = l + 1;
04827 if (ib > nn) {
04828 return 0;
04829 }
04830 goto L100;
04831 L150:
04832 *nz = *n;
04833 if (*fnu == 0.) {
04834 --(*nz);
04835 }
04836 L160:
04837 yr[1] = zeror;
04838 yi[1] = zeroi;
04839 if (*fnu != 0.) {
04840 goto L170;
04841 }
04842 yr[1] = coner;
04843 yi[1] = conei;
04844 L170:
04845 if (*n == 1) {
04846 return 0;
04847 }
04848 i__1 = *n;
04849 for (i__ = 2; i__ <= i__1; ++i__) {
04850 yr[i__] = zeror;
04851 yi[i__] = zeroi;
04852
04853 }
04854 return 0;
04855
04856
04857
04858
04859 L190:
04860 *nz = -(*nz);
04861 return 0;
04862 }
04863
04864 int zshch_(doublereal *zr, doublereal *zi, doublereal *cshr, doublereal *cshi, doublereal *cchr, doublereal *cchi)
04865 {
04866
04867
04868
04869 static doublereal ch, cn, sh, sn;
04870
04871
04872
04873
04874
04875
04876
04877
04878
04879
04880 sh = sinh(*zr);
04881 ch = cosh(*zr);
04882 sn = sin(*zi);
04883 cn = cos(*zi);
04884 *cshr = sh * cn;
04885 *cshi = ch * sn;
04886 *cchr = ch * cn;
04887 *cchi = sh * sn;
04888 return 0;
04889 }
04890
04891 int zsqrt_(const doublereal *ar, const doublereal *ai, doublereal *br, doublereal *bi)
04892 {
04893
04894
04895 static doublereal drt = .7071067811865475244008443621;
04896 static doublereal dpi = 3.141592653589793238462643383;
04897
04898
04899
04900 static doublereal zm, dtheta;
04901
04902
04903
04904
04905
04906
04907
04908
04909 zm = myzabs_(ar, ai);
04910 zm = sqrt(zm);
04911 if (*ar == 0.) {
04912 goto L10;
04913 }
04914 if (*ai == 0.) {
04915 goto L20;
04916 }
04917 dtheta = atan(*ai / *ar);
04918 if (dtheta <= 0.) {
04919 goto L40;
04920 }
04921 if (*ar < 0.) {
04922 dtheta -= dpi;
04923 }
04924 goto L50;
04925 L10:
04926 if (*ai > 0.) {
04927 goto L60;
04928 }
04929 if (*ai < 0.) {
04930 goto L70;
04931 }
04932 *br = 0.;
04933 *bi = 0.;
04934 return 0;
04935 L20:
04936 if (*ar > 0.) {
04937 goto L30;
04938 }
04939 *br = 0.;
04940 *bi = sqrt((abs(*ar)));
04941 return 0;
04942 L30:
04943 *br = sqrt(*ar);
04944 *bi = 0.;
04945 return 0;
04946 L40:
04947 if (*ar < 0.) {
04948 dtheta += dpi;
04949 }
04950 L50:
04951 dtheta *= .5;
04952 *br = zm * cos(dtheta);
04953 *bi = zm * sin(dtheta);
04954 return 0;
04955 L60:
04956 *br = zm * drt;
04957 *bi = zm * drt;
04958 return 0;
04959 L70:
04960 *br = zm * drt;
04961 *bi = -zm * drt;
04962 return 0;
04963 }
04964
04965 int zuchk_(doublereal *yr, doublereal *yi, integer *nz, doublereal *ascle, doublereal *tol)
04966 {
04967 static doublereal wi, ss, st, wr;
04968
04969
04970
04971
04972
04973
04974
04975
04976
04977
04978
04979
04980
04981
04982
04983
04984 *nz = 0;
04985 wr = abs(*yr);
04986 wi = abs(*yi);
04987 st = min(wr,wi);
04988 if (st > *ascle) {
04989 return 0;
04990 }
04991 ss = max(wr,wi);
04992 st /= *tol;
04993 if (ss < st) {
04994 *nz = 1;
04995 }
04996 return 0;
04997 }
04998
04999 int zunhj_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *ipmtr, doublereal *tol, doublereal *phir, doublereal *phii, doublereal *argr, doublereal *argi, doublereal *zeta1r, doublereal *zeta1i, doublereal *zeta2r, doublereal *zeta2i, doublereal *asumr, doublereal *asumi, doublereal *bsumr, doublereal *bsumi)
05000 {
05001
05002
05003 static doublereal ar[14] = { 1.,.104166666666666667,.0835503472222222222,
05004 .12822657455632716,.291849026464140464,.881627267443757652,
05005 3.32140828186276754,14.9957629868625547,78.9230130115865181,
05006 474.451538868264323,3207.49009089066193,24086.5496408740049,
05007 198923.119169509794,1791902.00777534383 };
05008 static doublereal br[14] = { 1.,-.145833333333333333,
05009 -.0987413194444444444,-.143312053915895062,-.317227202678413548,
05010 -.942429147957120249,-3.51120304082635426,-15.7272636203680451,
05011 -82.2814390971859444,-492.355370523670524,-3316.21856854797251,
05012 -24827.6742452085896,-204526.587315129788,-1838444.9170682099 };
05013 static doublereal c__[105] = { 1.,-.208333333333333333,.125,
05014 .334201388888888889,-.401041666666666667,.0703125,
05015 -1.02581259645061728,1.84646267361111111,-.8912109375,.0732421875,
05016 4.66958442342624743,-11.2070026162229938,8.78912353515625,
05017 -2.3640869140625,.112152099609375,-28.2120725582002449,
05018 84.6362176746007346,-91.8182415432400174,42.5349987453884549,
05019 -7.3687943594796317,.227108001708984375,212.570130039217123,
05020 -765.252468141181642,1059.99045252799988,-699.579627376132541,
05021 218.19051174421159,-26.4914304869515555,.572501420974731445,
05022 -1919.457662318407,8061.72218173730938,-13586.5500064341374,
05023 11655.3933368645332,-5305.64697861340311,1200.90291321635246,
05024 -108.090919788394656,1.7277275025844574,20204.2913309661486,
05025 -96980.5983886375135,192547.001232531532,-203400.177280415534,
05026 122200.46498301746,-41192.6549688975513,7109.51430248936372,
05027 -493.915304773088012,6.07404200127348304,-242919.187900551333,
05028 1311763.6146629772,-2998015.91853810675,3763271.297656404,
05029 -2813563.22658653411,1268365.27332162478,-331645.172484563578,
05030 45218.7689813627263,-2499.83048181120962,24.3805296995560639,
05031 3284469.85307203782,-19706819.1184322269,50952602.4926646422,
05032 -74105148.2115326577,66344512.2747290267,-37567176.6607633513,
05033 13288767.1664218183,-2785618.12808645469,308186.404612662398,
05034 -13886.0897537170405,110.017140269246738,-49329253.664509962,
05035 325573074.185765749,-939462359.681578403,1553596899.57058006,
05036 -1621080552.10833708,1106842816.82301447,-495889784.275030309,
05037 142062907.797533095,-24474062.7257387285,2243768.17792244943,
05038 -84005.4336030240853,551.335896122020586,814789096.118312115,
05039 -5866481492.05184723,18688207509.2958249,-34632043388.1587779,
05040 41280185579.753974,-33026599749.8007231,17954213731.1556001,
05041 -6563293792.61928433,1559279864.87925751,-225105661.889415278,
05042 17395107.5539781645,-549842.327572288687,3038.09051092238427,
05043 -14679261247.6956167,114498237732.02581,-399096175224.466498,
05044 819218669548.577329,-1098375156081.22331,1008158106865.38209,
05045 -645364869245.376503,287900649906.150589,-87867072178.0232657,
05046 17634730606.8349694,-2167164983.22379509,143157876.718888981,
05047 -3871833.44257261262,18257.7554742931747 };
05048 static doublereal alfa[180] = { -.00444444444444444444,
05049 -9.22077922077922078e-4,-8.84892884892884893e-5,
05050 1.65927687832449737e-4,2.4669137274179291e-4,
05051 2.6599558934625478e-4,2.61824297061500945e-4,
05052 2.48730437344655609e-4,2.32721040083232098e-4,
05053 2.16362485712365082e-4,2.00738858762752355e-4,
05054 1.86267636637545172e-4,1.73060775917876493e-4,
05055 1.61091705929015752e-4,1.50274774160908134e-4,
05056 1.40503497391269794e-4,1.31668816545922806e-4,
05057 1.23667445598253261e-4,1.16405271474737902e-4,
05058 1.09798298372713369e-4,1.03772410422992823e-4,
05059 9.82626078369363448e-5,9.32120517249503256e-5,
05060 8.85710852478711718e-5,8.42963105715700223e-5,
05061 8.03497548407791151e-5,7.66981345359207388e-5,
05062 7.33122157481777809e-5,7.01662625163141333e-5,
05063 6.72375633790160292e-5,6.93735541354588974e-4,
05064 2.32241745182921654e-4,-1.41986273556691197e-5,
05065 -1.1644493167204864e-4,-1.50803558053048762e-4,
05066 -1.55121924918096223e-4,-1.46809756646465549e-4,
05067 -1.33815503867491367e-4,-1.19744975684254051e-4,
05068 -1.0618431920797402e-4,-9.37699549891194492e-5,
05069 -8.26923045588193274e-5,-7.29374348155221211e-5,
05070 -6.44042357721016283e-5,-5.69611566009369048e-5,
05071 -5.04731044303561628e-5,-4.48134868008882786e-5,
05072 -3.98688727717598864e-5,-3.55400532972042498e-5,
05073 -3.1741425660902248e-5,-2.83996793904174811e-5,
05074 -2.54522720634870566e-5,-2.28459297164724555e-5,
05075 -2.05352753106480604e-5,-1.84816217627666085e-5,
05076 -1.66519330021393806e-5,-1.50179412980119482e-5,
05077 -1.35554031379040526e-5,-1.22434746473858131e-5,
05078 -1.10641884811308169e-5,-3.54211971457743841e-4,
05079 -1.56161263945159416e-4,3.0446550359493641e-5,
05080 1.30198655773242693e-4,1.67471106699712269e-4,
05081 1.70222587683592569e-4,1.56501427608594704e-4,
05082 1.3633917097744512e-4,1.14886692029825128e-4,
05083 9.45869093034688111e-5,7.64498419250898258e-5,
05084 6.07570334965197354e-5,4.74394299290508799e-5,
05085 3.62757512005344297e-5,2.69939714979224901e-5,
05086 1.93210938247939253e-5,1.30056674793963203e-5,
05087 7.82620866744496661e-6,3.59257485819351583e-6,
05088 1.44040049814251817e-7,-2.65396769697939116e-6,
05089 -4.9134686709848591e-6,-6.72739296091248287e-6,
05090 -8.17269379678657923e-6,-9.31304715093561232e-6,
05091 -1.02011418798016441e-5,-1.0880596251059288e-5,
05092 -1.13875481509603555e-5,-1.17519675674556414e-5,
05093 -1.19987364870944141e-5,3.78194199201772914e-4,
05094 2.02471952761816167e-4,-6.37938506318862408e-5,
05095 -2.38598230603005903e-4,-3.10916256027361568e-4,
05096 -3.13680115247576316e-4,-2.78950273791323387e-4,
05097 -2.28564082619141374e-4,-1.75245280340846749e-4,
05098 -1.25544063060690348e-4,-8.22982872820208365e-5,
05099 -4.62860730588116458e-5,-1.72334302366962267e-5,
05100 5.60690482304602267e-6,2.313954431482868e-5,
05101 3.62642745856793957e-5,4.58006124490188752e-5,
05102 5.2459529495911405e-5,5.68396208545815266e-5,
05103 5.94349820393104052e-5,6.06478527578421742e-5,
05104 6.08023907788436497e-5,6.01577894539460388e-5,
05105 5.891996573446985e-5,5.72515823777593053e-5,
05106 5.52804375585852577e-5,5.3106377380288017e-5,
05107 5.08069302012325706e-5,4.84418647620094842e-5,
05108 4.6056858160747537e-5,-6.91141397288294174e-4,
05109 -4.29976633058871912e-4,1.83067735980039018e-4,
05110 6.60088147542014144e-4,8.75964969951185931e-4,
05111 8.77335235958235514e-4,7.49369585378990637e-4,
05112 5.63832329756980918e-4,3.68059319971443156e-4,
05113 1.88464535514455599e-4,3.70663057664904149e-5,
05114 -8.28520220232137023e-5,-1.72751952869172998e-4,
05115 -2.36314873605872983e-4,-2.77966150694906658e-4,
05116 -3.02079514155456919e-4,-3.12594712643820127e-4,
05117 -3.12872558758067163e-4,-3.05678038466324377e-4,
05118 -2.93226470614557331e-4,-2.77255655582934777e-4,
05119 -2.59103928467031709e-4,-2.39784014396480342e-4,
05120 -2.20048260045422848e-4,-2.00443911094971498e-4,
05121 -1.81358692210970687e-4,-1.63057674478657464e-4,
05122 -1.45712672175205844e-4,-1.29425421983924587e-4,
05123 -1.14245691942445952e-4,.00192821964248775885,
05124 .00135592576302022234,-7.17858090421302995e-4,
05125 -.00258084802575270346,-.00349271130826168475,
05126 -.00346986299340960628,-.00282285233351310182,
05127 -.00188103076404891354,-8.895317183839476e-4,
05128 3.87912102631035228e-6,7.28688540119691412e-4,
05129 .00126566373053457758,.00162518158372674427,.00183203153216373172,
05130 .00191588388990527909,.00190588846755546138,.00182798982421825727,
05131 .0017038950642112153,.00155097127171097686,.00138261421852276159,
05132 .00120881424230064774,.00103676532638344962,
05133 8.71437918068619115e-4,7.16080155297701002e-4,
05134 5.72637002558129372e-4,4.42089819465802277e-4,
05135 3.24724948503090564e-4,2.20342042730246599e-4,
05136 1.28412898401353882e-4,4.82005924552095464e-5 };
05137 static doublereal beta[210] = { .0179988721413553309,
05138 .00559964911064388073,.00288501402231132779,.00180096606761053941,
05139 .00124753110589199202,9.22878876572938311e-4,
05140 7.14430421727287357e-4,5.71787281789704872e-4,
05141 4.69431007606481533e-4,3.93232835462916638e-4,
05142 3.34818889318297664e-4,2.88952148495751517e-4,
05143 2.52211615549573284e-4,2.22280580798883327e-4,
05144 1.97541838033062524e-4,1.76836855019718004e-4,
05145 1.59316899661821081e-4,1.44347930197333986e-4,
05146 1.31448068119965379e-4,1.20245444949302884e-4,
05147 1.10449144504599392e-4,1.01828770740567258e-4,
05148 9.41998224204237509e-5,8.74130545753834437e-5,
05149 8.13466262162801467e-5,7.59002269646219339e-5,
05150 7.09906300634153481e-5,6.65482874842468183e-5,
05151 6.25146958969275078e-5,5.88403394426251749e-5,
05152 -.00149282953213429172,-8.78204709546389328e-4,
05153 -5.02916549572034614e-4,-2.94822138512746025e-4,
05154 -1.75463996970782828e-4,-1.04008550460816434e-4,
05155 -5.96141953046457895e-5,-3.1203892907609834e-5,
05156 -1.26089735980230047e-5,-2.42892608575730389e-7,
05157 8.05996165414273571e-6,1.36507009262147391e-5,
05158 1.73964125472926261e-5,1.9867297884213378e-5,
05159 2.14463263790822639e-5,2.23954659232456514e-5,
05160 2.28967783814712629e-5,2.30785389811177817e-5,
05161 2.30321976080909144e-5,2.28236073720348722e-5,
05162 2.25005881105292418e-5,2.20981015361991429e-5,
05163 2.16418427448103905e-5,2.11507649256220843e-5,
05164 2.06388749782170737e-5,2.01165241997081666e-5,
05165 1.95913450141179244e-5,1.9068936791043674e-5,
05166 1.85533719641636667e-5,1.80475722259674218e-5,
05167 5.5221307672129279e-4,4.47932581552384646e-4,
05168 2.79520653992020589e-4,1.52468156198446602e-4,
05169 6.93271105657043598e-5,1.76258683069991397e-5,
05170 -1.35744996343269136e-5,-3.17972413350427135e-5,
05171 -4.18861861696693365e-5,-4.69004889379141029e-5,
05172 -4.87665447413787352e-5,-4.87010031186735069e-5,
05173 -4.74755620890086638e-5,-4.55813058138628452e-5,
05174 -4.33309644511266036e-5,-4.09230193157750364e-5,
05175 -3.84822638603221274e-5,-3.60857167535410501e-5,
05176 -3.37793306123367417e-5,-3.15888560772109621e-5,
05177 -2.95269561750807315e-5,-2.75978914828335759e-5,
05178 -2.58006174666883713e-5,-2.413083567612802e-5,
05179 -2.25823509518346033e-5,-2.11479656768912971e-5,
05180 -1.98200638885294927e-5,-1.85909870801065077e-5,
05181 -1.74532699844210224e-5,-1.63997823854497997e-5,
05182 -4.74617796559959808e-4,-4.77864567147321487e-4,
05183 -3.20390228067037603e-4,-1.61105016119962282e-4,
05184 -4.25778101285435204e-5,3.44571294294967503e-5,
05185 7.97092684075674924e-5,1.031382367082722e-4,
05186 1.12466775262204158e-4,1.13103642108481389e-4,
05187 1.08651634848774268e-4,1.01437951597661973e-4,
05188 9.29298396593363896e-5,8.40293133016089978e-5,
05189 7.52727991349134062e-5,6.69632521975730872e-5,
05190 5.92564547323194704e-5,5.22169308826975567e-5,
05191 4.58539485165360646e-5,4.01445513891486808e-5,
05192 3.50481730031328081e-5,3.05157995034346659e-5,
05193 2.64956119950516039e-5,2.29363633690998152e-5,
05194 1.97893056664021636e-5,1.70091984636412623e-5,
05195 1.45547428261524004e-5,1.23886640995878413e-5,
05196 1.04775876076583236e-5,8.79179954978479373e-6,
05197 7.36465810572578444e-4,8.72790805146193976e-4,
05198 6.22614862573135066e-4,2.85998154194304147e-4,
05199 3.84737672879366102e-6,-1.87906003636971558e-4,
05200 -2.97603646594554535e-4,-3.45998126832656348e-4,
05201 -3.53382470916037712e-4,-3.35715635775048757e-4,
05202 -3.04321124789039809e-4,-2.66722723047612821e-4,
05203 -2.27654214122819527e-4,-1.89922611854562356e-4,
05204 -1.5505891859909387e-4,-1.2377824076187363e-4,
05205 -9.62926147717644187e-5,-7.25178327714425337e-5,
05206 -5.22070028895633801e-5,-3.50347750511900522e-5,
05207 -2.06489761035551757e-5,-8.70106096849767054e-6,
05208 1.1369868667510029e-6,9.16426474122778849e-6,
05209 1.5647778542887262e-5,2.08223629482466847e-5,
05210 2.48923381004595156e-5,2.80340509574146325e-5,
05211 3.03987774629861915e-5,3.21156731406700616e-5,
05212 -.00180182191963885708,-.00243402962938042533,
05213 -.00183422663549856802,-7.62204596354009765e-4,
05214 2.39079475256927218e-4,9.49266117176881141e-4,
05215 .00134467449701540359,.00148457495259449178,.00144732339830617591,
05216 .00130268261285657186,.00110351597375642682,
05217 8.86047440419791759e-4,6.73073208165665473e-4,
05218 4.77603872856582378e-4,3.05991926358789362e-4,
05219 1.6031569459472163e-4,4.00749555270613286e-5,
05220 -5.66607461635251611e-5,-1.32506186772982638e-4,
05221 -1.90296187989614057e-4,-2.32811450376937408e-4,
05222 -2.62628811464668841e-4,-2.82050469867598672e-4,
05223 -2.93081563192861167e-4,-2.97435962176316616e-4,
05224 -2.96557334239348078e-4,-2.91647363312090861e-4,
05225 -2.83696203837734166e-4,-2.73512317095673346e-4,
05226 -2.6175015580676858e-4,.00638585891212050914,
05227 .00962374215806377941,.00761878061207001043,.00283219055545628054,
05228 -.0020984135201272009,-.00573826764216626498,
05229 -.0077080424449541462,-.00821011692264844401,
05230 -.00765824520346905413,-.00647209729391045177,
05231 -.00499132412004966473,-.0034561228971313328,
05232 -.00201785580014170775,-7.59430686781961401e-4,
05233 2.84173631523859138e-4,.00110891667586337403,
05234 .00172901493872728771,.00216812590802684701,.00245357710494539735,
05235 .00261281821058334862,.00267141039656276912,.0026520307339598043,
05236 .00257411652877287315,.00245389126236094427,.00230460058071795494,
05237 .00213684837686712662,.00195896528478870911,.00177737008679454412,
05238 .00159690280765839059,.00142111975664438546 };
05239 static doublereal gama[30] = { .629960524947436582,.251984209978974633,
05240 .154790300415655846,.110713062416159013,.0857309395527394825,
05241 .0697161316958684292,.0586085671893713576,.0504698873536310685,
05242 .0442600580689154809,.0393720661543509966,.0354283195924455368,
05243 .0321818857502098231,.0294646240791157679,.0271581677112934479,
05244 .0251768272973861779,.0234570755306078891,.0219508390134907203,
05245 .020621082823564624,.0194388240897880846,.0183810633800683158,
05246 .0174293213231963172,.0165685837786612353,.0157865285987918445,
05247 .0150729501494095594,.0144193250839954639,.0138184805735341786,
05248 .0132643378994276568,.0127517121970498651,.0122761545318762767,
05249 .0118338262398482403 };
05250 static doublereal ex1 = .333333333333333333;
05251 static doublereal ex2 = .666666666666666667;
05252 static doublereal hpi = 1.57079632679489662;
05253 static doublereal gpi = 3.14159265358979324;
05254 static doublereal thpi = 4.71238898038468986;
05255 static doublereal zeror = 0.;
05256 static doublereal zeroi = 0.;
05257 static doublereal coner = 1.;
05258 static doublereal conei = 0.;
05259
05260
05261 integer i__1, i__2;
05262 doublereal d__1;
05263
05264
05265
05266
05267 static doublereal rfn13;
05268 static integer idum;
05269 static doublereal atol, btol, tfni;
05270 static integer kmax;
05271 static doublereal azth, tzai, tfnr, rfnu;
05272 static doublereal zthi, test, tzar, zthr, rfnu2;
05273 static integer j, k, l, m;
05274 static doublereal zetai, ptfni, sumai, sumbi, zetar, ptfnr, razth, sumar,
05275 sumbr, rzthi;
05276 static integer l1, l2;
05277 static doublereal rzthr, rtzti;
05278 static doublereal rtztr, ac, ap[30], pi[30];
05279 static integer is, jr, ks, ju;
05280 static doublereal pp, wi, pr[30];
05281 static integer lr;
05282 static doublereal wr;
05283 static doublereal aw2;
05284 static integer kp1;
05285 static doublereal przthi, t2i, w2i, t2r, przthr, w2r, ang, fn13, fn23;
05286 static integer ias;
05287 static doublereal cri[14], dri[14];
05288 static integer ibs;
05289 static doublereal zai, zbi, zci, crr[14], drr[14], raw, zar, upi[14], sti,
05290 zbr, zcr, upr[14], str, raw2;
05291 static integer lrp1;
05292
05293
05294
05295
05296
05297
05298
05299
05300
05301
05302
05303
05304
05305
05306
05307
05308
05309
05310
05311
05312
05313
05314
05315
05316
05317
05318
05319
05320
05321
05322
05323
05324
05325
05326
05327
05328 rfnu = 1. / *fnu;
05329
05330
05331
05332 test = d1mach_(&c__1) * 1e3;
05333 ac = *fnu * test;
05334 if (abs(*zr) > ac || abs(*zi) > ac) {
05335 goto L15;
05336 }
05337 *zeta1r = (d__1 = log(test), abs(d__1)) * 2. + *fnu;
05338 *zeta1i = 0.;
05339 *zeta2r = *fnu;
05340 *zeta2i = 0.;
05341 *phir = 1.;
05342 *phii = 0.;
05343 *argr = 1.;
05344 *argi = 0.;
05345 return 0;
05346 L15:
05347 zbr = *zr * rfnu;
05348 zbi = *zi * rfnu;
05349 rfnu2 = rfnu * rfnu;
05350
05351
05352
05353 fn13 = pow_dd(fnu, &ex1);
05354 fn23 = fn13 * fn13;
05355 rfn13 = 1. / fn13;
05356 w2r = coner - zbr * zbr + zbi * zbi;
05357 w2i = conei - zbr * zbi - zbr * zbi;
05358 aw2 = myzabs_(&w2r, &w2i);
05359 if (aw2 > .25) {
05360 goto L130;
05361 }
05362
05363
05364
05365 k = 1;
05366 pr[0] = coner;
05367 pi[0] = conei;
05368 sumar = gama[0];
05369 sumai = zeroi;
05370 ap[0] = 1.;
05371 if (aw2 < *tol) {
05372 goto L20;
05373 }
05374 for (k = 2; k <= 30; ++k) {
05375 pr[k - 1] = pr[k - 2] * w2r - pi[k - 2] * w2i;
05376 pi[k - 1] = pr[k - 2] * w2i + pi[k - 2] * w2r;
05377 sumar += pr[k - 1] * gama[k - 1];
05378 sumai += pi[k - 1] * gama[k - 1];
05379 ap[k - 1] = ap[k - 2] * aw2;
05380 if (ap[k - 1] < *tol) {
05381 goto L20;
05382 }
05383
05384 }
05385 k = 30;
05386 L20:
05387 kmax = k;
05388 zetar = w2r * sumar - w2i * sumai;
05389 zetai = w2r * sumai + w2i * sumar;
05390 *argr = zetar * fn23;
05391 *argi = zetai * fn23;
05392 zsqrt_(&sumar, &sumai, &zar, &zai);
05393 zsqrt_(&w2r, &w2i, &str, &sti);
05394 *zeta2r = str * *fnu;
05395 *zeta2i = sti * *fnu;
05396 str = coner + ex2 * (zetar * zar - zetai * zai);
05397 sti = conei + ex2 * (zetar * zai + zetai * zar);
05398 *zeta1r = str * *zeta2r - sti * *zeta2i;
05399 *zeta1i = str * *zeta2i + sti * *zeta2r;
05400 zar += zar;
05401 zai += zai;
05402 zsqrt_(&zar, &zai, &str, &sti);
05403 *phir = str * rfn13;
05404 *phii = sti * rfn13;
05405 if (*ipmtr == 1) {
05406 goto L120;
05407 }
05408
05409
05410
05411 sumbr = zeror;
05412 sumbi = zeroi;
05413 i__1 = kmax;
05414 for (k = 1; k <= i__1; ++k) {
05415 sumbr += pr[k - 1] * beta[k - 1];
05416 sumbi += pi[k - 1] * beta[k - 1];
05417
05418 }
05419 *asumr = zeror;
05420 *asumi = zeroi;
05421 *bsumr = sumbr;
05422 *bsumi = sumbi;
05423 l1 = 0;
05424 l2 = 30;
05425 btol = *tol * (abs(*bsumr) + abs(*bsumi));
05426 atol = *tol;
05427 pp = 1.;
05428 ias = 0;
05429 ibs = 0;
05430 if (rfnu2 < *tol) {
05431 goto L110;
05432 }
05433 for (is = 2; is <= 7; ++is) {
05434 atol /= rfnu2;
05435 pp *= rfnu2;
05436 if (ias == 1) {
05437 goto L60;
05438 }
05439 sumar = zeror;
05440 sumai = zeroi;
05441 i__1 = kmax;
05442 for (k = 1; k <= i__1; ++k) {
05443 m = l1 + k;
05444 sumar += pr[k - 1] * alfa[m - 1];
05445 sumai += pi[k - 1] * alfa[m - 1];
05446 if (ap[k - 1] < atol) {
05447 goto L50;
05448 }
05449
05450 }
05451 L50:
05452 *asumr += sumar * pp;
05453 *asumi += sumai * pp;
05454 if (pp < *tol) {
05455 ias = 1;
05456 }
05457 L60:
05458 if (ibs == 1) {
05459 goto L90;
05460 }
05461 sumbr = zeror;
05462 sumbi = zeroi;
05463 i__1 = kmax;
05464 for (k = 1; k <= i__1; ++k) {
05465 m = l2 + k;
05466 sumbr += pr[k - 1] * beta[m - 1];
05467 sumbi += pi[k - 1] * beta[m - 1];
05468 if (ap[k - 1] < atol) {
05469 goto L80;
05470 }
05471
05472 }
05473 L80:
05474 *bsumr += sumbr * pp;
05475 *bsumi += sumbi * pp;
05476 if (pp < btol) {
05477 ibs = 1;
05478 }
05479 L90:
05480 if (ias == 1 && ibs == 1) {
05481 goto L110;
05482 }
05483 l1 += 30;
05484 l2 += 30;
05485
05486 }
05487 L110:
05488 *asumr += coner;
05489 pp = rfnu * rfn13;
05490 *bsumr *= pp;
05491 *bsumi *= pp;
05492 L120:
05493 return 0;
05494
05495
05496
05497 L130:
05498 zsqrt_(&w2r, &w2i, &wr, &wi);
05499 if (wr < 0.) {
05500 wr = 0.;
05501 }
05502 if (wi < 0.) {
05503 wi = 0.;
05504 }
05505 str = coner + wr;
05506 sti = wi;
05507 zdiv_(&str, &sti, &zbr, &zbi, &zar, &zai);
05508 zlog_(&zar, &zai, &zcr, &zci, &idum);
05509 if (zci < 0.) {
05510 zci = 0.;
05511 }
05512 if (zci > hpi) {
05513 zci = hpi;
05514 }
05515 if (zcr < 0.) {
05516 zcr = 0.;
05517 }
05518 zthr = (zcr - wr) * 1.5;
05519 zthi = (zci - wi) * 1.5;
05520 *zeta1r = zcr * *fnu;
05521 *zeta1i = zci * *fnu;
05522 *zeta2r = wr * *fnu;
05523 *zeta2i = wi * *fnu;
05524 azth = myzabs_(&zthr, &zthi);
05525 ang = thpi;
05526 if (zthr >= 0. && zthi < 0.) {
05527 goto L140;
05528 }
05529 ang = hpi;
05530 if (zthr == 0.) {
05531 goto L140;
05532 }
05533 ang = atan(zthi / zthr);
05534 if (zthr < 0.) {
05535 ang += gpi;
05536 }
05537 L140:
05538 pp = pow_dd(&azth, &ex2);
05539 ang *= ex2;
05540 zetar = pp * cos(ang);
05541 zetai = pp * sin(ang);
05542 if (zetai < 0.) {
05543 zetai = 0.;
05544 }
05545 *argr = zetar * fn23;
05546 *argi = zetai * fn23;
05547 zdiv_(&zthr, &zthi, &zetar, &zetai, &rtztr, &rtzti);
05548 zdiv_(&rtztr, &rtzti, &wr, &wi, &zar, &zai);
05549 tzar = zar + zar;
05550 tzai = zai + zai;
05551 zsqrt_(&tzar, &tzai, &str, &sti);
05552 *phir = str * rfn13;
05553 *phii = sti * rfn13;
05554 if (*ipmtr == 1) {
05555 goto L120;
05556 }
05557 raw = 1. / sqrt(aw2);
05558 str = wr * raw;
05559 sti = -wi * raw;
05560 tfnr = str * rfnu * raw;
05561 tfni = sti * rfnu * raw;
05562 razth = 1. / azth;
05563 str = zthr * razth;
05564 sti = -zthi * razth;
05565 rzthr = str * razth * rfnu;
05566 rzthi = sti * razth * rfnu;
05567 zcr = rzthr * ar[1];
05568 zci = rzthi * ar[1];
05569 raw2 = 1. / aw2;
05570 str = w2r * raw2;
05571 sti = -w2i * raw2;
05572 t2r = str * raw2;
05573 t2i = sti * raw2;
05574 str = t2r * c__[1] + c__[2];
05575 sti = t2i * c__[1];
05576 upr[1] = str * tfnr - sti * tfni;
05577 upi[1] = str * tfni + sti * tfnr;
05578 *bsumr = upr[1] + zcr;
05579 *bsumi = upi[1] + zci;
05580 *asumr = zeror;
05581 *asumi = zeroi;
05582 if (rfnu < *tol) {
05583 goto L220;
05584 }
05585 przthr = rzthr;
05586 przthi = rzthi;
05587 ptfnr = tfnr;
05588 ptfni = tfni;
05589 upr[0] = coner;
05590 upi[0] = conei;
05591 pp = 1.;
05592 btol = *tol * (abs(*bsumr) + abs(*bsumi));
05593 ks = 0;
05594 kp1 = 2;
05595 l = 3;
05596 ias = 0;
05597 ibs = 0;
05598 for (lr = 2; lr <= 12; lr += 2) {
05599 lrp1 = lr + 1;
05600
05601
05602
05603
05604 i__1 = lrp1;
05605 for (k = lr; k <= i__1; ++k) {
05606 ++ks;
05607 ++kp1;
05608 ++l;
05609 zar = c__[l - 1];
05610 zai = zeroi;
05611 i__2 = kp1;
05612 for (j = 2; j <= i__2; ++j) {
05613 ++l;
05614 str = zar * t2r - t2i * zai + c__[l - 1];
05615 zai = zar * t2i + zai * t2r;
05616 zar = str;
05617
05618 }
05619 str = ptfnr * tfnr - ptfni * tfni;
05620 ptfni = ptfnr * tfni + ptfni * tfnr;
05621 ptfnr = str;
05622 upr[kp1 - 1] = ptfnr * zar - ptfni * zai;
05623 upi[kp1 - 1] = ptfni * zar + ptfnr * zai;
05624 crr[ks - 1] = przthr * br[ks];
05625 cri[ks - 1] = przthi * br[ks];
05626 str = przthr * rzthr - przthi * rzthi;
05627 przthi = przthr * rzthi + przthi * rzthr;
05628 przthr = str;
05629 drr[ks - 1] = przthr * ar[ks + 1];
05630 dri[ks - 1] = przthi * ar[ks + 1];
05631
05632 }
05633 pp *= rfnu2;
05634 if (ias == 1) {
05635 goto L180;
05636 }
05637 sumar = upr[lrp1 - 1];
05638 sumai = upi[lrp1 - 1];
05639 ju = lrp1;
05640 i__1 = lr;
05641 for (jr = 1; jr <= i__1; ++jr) {
05642 --ju;
05643 sumar = sumar + crr[jr - 1] * upr[ju - 1] - cri[jr - 1] * upi[ju
05644 - 1];
05645 sumai = sumai + crr[jr - 1] * upi[ju - 1] + cri[jr - 1] * upr[ju
05646 - 1];
05647
05648 }
05649 *asumr += sumar;
05650 *asumi += sumai;
05651 test = abs(sumar) + abs(sumai);
05652 if (pp < *tol && test < *tol) {
05653 ias = 1;
05654 }
05655 L180:
05656 if (ibs == 1) {
05657 goto L200;
05658 }
05659 sumbr = upr[lr + 1] + upr[lrp1 - 1] * zcr - upi[lrp1 - 1] * zci;
05660 sumbi = upi[lr + 1] + upr[lrp1 - 1] * zci + upi[lrp1 - 1] * zcr;
05661 ju = lrp1;
05662 i__1 = lr;
05663 for (jr = 1; jr <= i__1; ++jr) {
05664 --ju;
05665 sumbr = sumbr + drr[jr - 1] * upr[ju - 1] - dri[jr - 1] * upi[ju
05666 - 1];
05667 sumbi = sumbi + drr[jr - 1] * upi[ju - 1] + dri[jr - 1] * upr[ju
05668 - 1];
05669
05670 }
05671 *bsumr += sumbr;
05672 *bsumi += sumbi;
05673 test = abs(sumbr) + abs(sumbi);
05674 if (pp < btol && test < btol) {
05675 ibs = 1;
05676 }
05677 L200:
05678 if (ias == 1 && ibs == 1) {
05679 goto L220;
05680 }
05681
05682 }
05683 L220:
05684 *asumr += coner;
05685 str = -(*bsumr) * rfn13;
05686 sti = -(*bsumi) * rfn13;
05687 zdiv_(&str, &sti, &rtztr, &rtzti, bsumr, bsumi);
05688 goto L120;
05689 }
05690
05691 int zuni1_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *n, doublereal *yr, doublereal *yi, integer *nz, integer *nlast, doublereal *fnul, doublereal *tol, doublereal *elim, doublereal *alim)
05692 {
05693
05694
05695 static doublereal zeror = 0.;
05696 static doublereal zeroi = 0.;
05697 static doublereal coner = 1.;
05698
05699
05700 integer i__1;
05701
05702
05703
05704
05705 static doublereal aphi, cscl, phii, crsc, phir;
05706 static integer init;
05707 static doublereal csrr[3], cssr[3], rast, sumi, sumr;
05708 static integer i__, k, m, iflag;
05709 static doublereal ascle, cwrki[16];
05710 static doublereal cwrkr[16];
05711 static doublereal zeta1i, zeta2i, zeta1r, zeta2r;
05712 static integer nd;
05713 static doublereal fn;
05714 static integer nn, nw;
05715 static doublereal c2i, c2m;
05716 static doublereal c1r, c2r, s1i, s2i, rs1, s1r, s2r, cyi[2];
05717 static integer nuf;
05718 static doublereal bry[3], cyr[2], sti, rzi, str, rzr;
05719
05720
05721
05722
05723
05724
05725
05726
05727
05728
05729
05730
05731
05732
05733
05734
05735
05736
05737 --yi;
05738 --yr;
05739
05740
05741
05742 *nz = 0;
05743 nd = *n;
05744 *nlast = 0;
05745
05746
05747
05748
05749
05750 cscl = 1. / *tol;
05751 crsc = *tol;
05752 cssr[0] = cscl;
05753 cssr[1] = coner;
05754 cssr[2] = crsc;
05755 csrr[0] = crsc;
05756 csrr[1] = coner;
05757 csrr[2] = cscl;
05758 bry[0] = d1mach_(&c__1) * 1e3 / *tol;
05759
05760
05761
05762 fn = max(*fnu,1.);
05763 init = 0;
05764 zunik_(zr, zi, &fn, &c__1, &c__1, tol, &init, &phir, &phii, &zeta1r, &
05765 zeta1i, &zeta2r, &zeta2i, &sumr, &sumi, cwrkr, cwrki);
05766 if (*kode == 1) {
05767 goto L10;
05768 }
05769 str = *zr + zeta2r;
05770 sti = *zi + zeta2i;
05771 rast = fn / myzabs_(&str, &sti);
05772 str = str * rast * rast;
05773 sti = -sti * rast * rast;
05774 s1r = -zeta1r + str;
05775 s1i = -zeta1i + sti;
05776 goto L20;
05777 L10:
05778 s1r = -zeta1r + zeta2r;
05779 s1i = -zeta1i + zeta2i;
05780 L20:
05781 rs1 = s1r;
05782 if (abs(rs1) > *elim) {
05783 goto L130;
05784 }
05785 L30:
05786 nn = min(2,nd);
05787 i__1 = nn;
05788 for (i__ = 1; i__ <= i__1; ++i__) {
05789 fn = *fnu + (doublereal) ((real) (nd - i__));
05790 init = 0;
05791 zunik_(zr, zi, &fn, &c__1, &c__0, tol, &init, &phir, &phii, &zeta1r, &
05792 zeta1i, &zeta2r, &zeta2i, &sumr, &sumi, cwrkr, cwrki);
05793 if (*kode == 1) {
05794 goto L40;
05795 }
05796 str = *zr + zeta2r;
05797 sti = *zi + zeta2i;
05798 rast = fn / myzabs_(&str, &sti);
05799 str = str * rast * rast;
05800 sti = -sti * rast * rast;
05801 s1r = -zeta1r + str;
05802 s1i = -zeta1i + sti + *zi;
05803 goto L50;
05804 L40:
05805 s1r = -zeta1r + zeta2r;
05806 s1i = -zeta1i + zeta2i;
05807 L50:
05808
05809
05810
05811 rs1 = s1r;
05812 if (abs(rs1) > *elim) {
05813 goto L110;
05814 }
05815 if (i__ == 1) {
05816 iflag = 2;
05817 }
05818 if (abs(rs1) < *alim) {
05819 goto L60;
05820 }
05821
05822
05823
05824 aphi = myzabs_(&phir, &phii);
05825 rs1 += log(aphi);
05826 if (abs(rs1) > *elim) {
05827 goto L110;
05828 }
05829 if (i__ == 1) {
05830 iflag = 1;
05831 }
05832 if (rs1 < 0.) {
05833 goto L60;
05834 }
05835 if (i__ == 1) {
05836 iflag = 3;
05837 }
05838 L60:
05839
05840
05841
05842 s2r = phir * sumr - phii * sumi;
05843 s2i = phir * sumi + phii * sumr;
05844 str = exp(s1r) * cssr[iflag - 1];
05845 s1r = str * cos(s1i);
05846 s1i = str * sin(s1i);
05847 str = s2r * s1r - s2i * s1i;
05848 s2i = s2r * s1i + s2i * s1r;
05849 s2r = str;
05850 if (iflag != 1) {
05851 goto L70;
05852 }
05853 zuchk_(&s2r, &s2i, &nw, bry, tol);
05854 if (nw != 0) {
05855 goto L110;
05856 }
05857 L70:
05858 cyr[i__ - 1] = s2r;
05859 cyi[i__ - 1] = s2i;
05860 m = nd - i__ + 1;
05861 yr[m] = s2r * csrr[iflag - 1];
05862 yi[m] = s2i * csrr[iflag - 1];
05863
05864 }
05865 if (nd <= 2) {
05866 goto L100;
05867 }
05868 rast = 1. / myzabs_(zr, zi);
05869 str = *zr * rast;
05870 sti = -(*zi) * rast;
05871 rzr = (str + str) * rast;
05872 rzi = (sti + sti) * rast;
05873 bry[1] = 1. / bry[0];
05874 bry[2] = d1mach_(&c__2);
05875 s1r = cyr[0];
05876 s1i = cyi[0];
05877 s2r = cyr[1];
05878 s2i = cyi[1];
05879 c1r = csrr[iflag - 1];
05880 ascle = bry[iflag - 1];
05881 k = nd - 2;
05882 fn = (doublereal) ((real) k);
05883 i__1 = nd;
05884 for (i__ = 3; i__ <= i__1; ++i__) {
05885 c2r = s2r;
05886 c2i = s2i;
05887 s2r = s1r + (*fnu + fn) * (rzr * c2r - rzi * c2i);
05888 s2i = s1i + (*fnu + fn) * (rzr * c2i + rzi * c2r);
05889 s1r = c2r;
05890 s1i = c2i;
05891 c2r = s2r * c1r;
05892 c2i = s2i * c1r;
05893 yr[k] = c2r;
05894 yi[k] = c2i;
05895 --k;
05896 fn += -1.;
05897 if (iflag >= 3) {
05898 goto L90;
05899 }
05900 str = abs(c2r);
05901 sti = abs(c2i);
05902 c2m = max(str,sti);
05903 if (c2m <= ascle) {
05904 goto L90;
05905 }
05906 ++iflag;
05907 ascle = bry[iflag - 1];
05908 s1r *= c1r;
05909 s1i *= c1r;
05910 s2r = c2r;
05911 s2i = c2i;
05912 s1r *= cssr[iflag - 1];
05913 s1i *= cssr[iflag - 1];
05914 s2r *= cssr[iflag - 1];
05915 s2i *= cssr[iflag - 1];
05916 c1r = csrr[iflag - 1];
05917 L90:
05918 ;
05919 }
05920 L100:
05921 return 0;
05922
05923
05924
05925 L110:
05926 if (rs1 > 0.) {
05927 goto L120;
05928 }
05929 yr[nd] = zeror;
05930 yi[nd] = zeroi;
05931 ++(*nz);
05932 --nd;
05933 if (nd == 0) {
05934 goto L100;
05935 }
05936 zuoik_(zr, zi, fnu, kode, &c__1, &nd, &yr[1], &yi[1], &nuf, tol, elim,
05937 alim);
05938 if (nuf < 0) {
05939 goto L120;
05940 }
05941 nd -= nuf;
05942 *nz += nuf;
05943 if (nd == 0) {
05944 goto L100;
05945 }
05946 fn = *fnu + (doublereal) ((real) (nd - 1));
05947 if (fn >= *fnul) {
05948 goto L30;
05949 }
05950 *nlast = nd;
05951 return 0;
05952 L120:
05953 *nz = -1;
05954 return 0;
05955 L130:
05956 if (rs1 > 0.) {
05957 goto L120;
05958 }
05959 *nz = *n;
05960 i__1 = *n;
05961 for (i__ = 1; i__ <= i__1; ++i__) {
05962 yr[i__] = zeror;
05963 yi[i__] = zeroi;
05964
05965 }
05966 return 0;
05967 }
05968
05969 int zuni2_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *n, doublereal *yr, doublereal *yi, integer *nz, integer *nlast, doublereal *fnul, doublereal *tol, doublereal *elim, doublereal *alim)
05970 {
05971
05972
05973 static doublereal zeror = 0.;
05974 static doublereal zeroi = 0.;
05975 static doublereal coner = 1.;
05976 static doublereal cipr[4] = { 1.,0.,-1.,0. };
05977 static doublereal cipi[4] = { 0.,1.,0.,-1. };
05978 static doublereal hpi = 1.57079632679489662;
05979 static doublereal aic = 1.265512123484645396;
05980
05981
05982 integer i__1;
05983
05984
05985
05986
05987 static doublereal daii, cidi, aarg;
05988 static integer ndai;
05989 static doublereal dair, aphi, argi, cscl, phii, crsc, argr;
05990 static integer idum;
05991 static doublereal phir, csrr[3], cssr[3], rast;
05992 static integer i__, j, k, iflag;
05993 static doublereal ascle, asumi, bsumi;
05994 static doublereal asumr, bsumr;
05995 static doublereal zeta1i, zeta2i, zeta1r, zeta2r;
05996 static integer nd;
05997 static doublereal fn;
05998 static integer in, nn, nw;
05999 static doublereal c2i, c2m;
06000 static doublereal c1r, c2r, s1i, s2i, rs1, s1r, s2r, aii, ang, car;
06001 static integer nai;
06002 static doublereal air, zbi, cyi[2], sar;
06003 static integer nuf, inu;
06004 static doublereal bry[3], raz, sti, zbr, zni, cyr[2], rzi, str, znr, rzr;
06005
06006
06007
06008
06009
06010
06011
06012
06013
06014
06015
06016
06017
06018
06019
06020
06021
06022
06023
06024 --yi;
06025 --yr;
06026
06027
06028
06029 *nz = 0;
06030 nd = *n;
06031 *nlast = 0;
06032
06033
06034
06035
06036
06037 cscl = 1. / *tol;
06038 crsc = *tol;
06039 cssr[0] = cscl;
06040 cssr[1] = coner;
06041 cssr[2] = crsc;
06042 csrr[0] = crsc;
06043 csrr[1] = coner;
06044 csrr[2] = cscl;
06045 bry[0] = d1mach_(&c__1) * 1e3 / *tol;
06046
06047
06048
06049 znr = *zi;
06050 zni = -(*zr);
06051 zbr = *zr;
06052 zbi = *zi;
06053 cidi = -coner;
06054 inu = (integer) ((real) (*fnu));
06055 ang = hpi * (*fnu - (doublereal) ((real) inu));
06056 c2r = cos(ang);
06057 c2i = sin(ang);
06058 car = c2r;
06059 sar = c2i;
06060 in = inu + *n - 1;
06061 in = in % 4 + 1;
06062 str = c2r * cipr[in - 1] - c2i * cipi[in - 1];
06063 c2i = c2r * cipi[in - 1] + c2i * cipr[in - 1];
06064 c2r = str;
06065 if (*zi > 0.) {
06066 goto L10;
06067 }
06068 znr = -znr;
06069 zbi = -zbi;
06070 cidi = -cidi;
06071 c2i = -c2i;
06072 L10:
06073
06074
06075
06076 fn = max(*fnu,1.);
06077 zunhj_(&znr, &zni, &fn, &c__1, tol, &phir, &phii, &argr, &argi, &zeta1r, &
06078 zeta1i, &zeta2r, &zeta2i, &asumr, &asumi, &bsumr, &bsumi);
06079 if (*kode == 1) {
06080 goto L20;
06081 }
06082 str = zbr + zeta2r;
06083 sti = zbi + zeta2i;
06084 rast = fn / myzabs_(&str, &sti);
06085 str = str * rast * rast;
06086 sti = -sti * rast * rast;
06087 s1r = -zeta1r + str;
06088 s1i = -zeta1i + sti;
06089 goto L30;
06090 L20:
06091 s1r = -zeta1r + zeta2r;
06092 s1i = -zeta1i + zeta2i;
06093 L30:
06094 rs1 = s1r;
06095 if (abs(rs1) > *elim) {
06096 goto L150;
06097 }
06098 L40:
06099 nn = min(2,nd);
06100 i__1 = nn;
06101 for (i__ = 1; i__ <= i__1; ++i__) {
06102 fn = *fnu + (doublereal) ((real) (nd - i__));
06103 zunhj_(&znr, &zni, &fn, &c__0, tol, &phir, &phii, &argr, &argi, &
06104 zeta1r, &zeta1i, &zeta2r, &zeta2i, &asumr, &asumi, &bsumr, &
06105 bsumi);
06106 if (*kode == 1) {
06107 goto L50;
06108 }
06109 str = zbr + zeta2r;
06110 sti = zbi + zeta2i;
06111 rast = fn / myzabs_(&str, &sti);
06112 str = str * rast * rast;
06113 sti = -sti * rast * rast;
06114 s1r = -zeta1r + str;
06115 s1i = -zeta1i + sti + abs(*zi);
06116 goto L60;
06117 L50:
06118 s1r = -zeta1r + zeta2r;
06119 s1i = -zeta1i + zeta2i;
06120 L60:
06121
06122
06123
06124 rs1 = s1r;
06125 if (abs(rs1) > *elim) {
06126 goto L120;
06127 }
06128 if (i__ == 1) {
06129 iflag = 2;
06130 }
06131 if (abs(rs1) < *alim) {
06132 goto L70;
06133 }
06134
06135
06136
06137
06138 aphi = myzabs_(&phir, &phii);
06139 aarg = myzabs_(&argr, &argi);
06140 rs1 = rs1 + log(aphi) - log(aarg) * .25 - aic;
06141 if (abs(rs1) > *elim) {
06142 goto L120;
06143 }
06144 if (i__ == 1) {
06145 iflag = 1;
06146 }
06147 if (rs1 < 0.) {
06148 goto L70;
06149 }
06150 if (i__ == 1) {
06151 iflag = 3;
06152 }
06153 L70:
06154
06155
06156
06157
06158 zairy_(&argr, &argi, &c__0, &c__2, &air, &aii, &nai, &idum);
06159 zairy_(&argr, &argi, &c__1, &c__2, &dair, &daii, &ndai, &idum);
06160 str = dair * bsumr - daii * bsumi;
06161 sti = dair * bsumi + daii * bsumr;
06162 str += air * asumr - aii * asumi;
06163 sti += air * asumi + aii * asumr;
06164 s2r = phir * str - phii * sti;
06165 s2i = phir * sti + phii * str;
06166 str = exp(s1r) * cssr[iflag - 1];
06167 s1r = str * cos(s1i);
06168 s1i = str * sin(s1i);
06169 str = s2r * s1r - s2i * s1i;
06170 s2i = s2r * s1i + s2i * s1r;
06171 s2r = str;
06172 if (iflag != 1) {
06173 goto L80;
06174 }
06175 zuchk_(&s2r, &s2i, &nw, bry, tol);
06176 if (nw != 0) {
06177 goto L120;
06178 }
06179 L80:
06180 if (*zi <= 0.) {
06181 s2i = -s2i;
06182 }
06183 str = s2r * c2r - s2i * c2i;
06184 s2i = s2r * c2i + s2i * c2r;
06185 s2r = str;
06186 cyr[i__ - 1] = s2r;
06187 cyi[i__ - 1] = s2i;
06188 j = nd - i__ + 1;
06189 yr[j] = s2r * csrr[iflag - 1];
06190 yi[j] = s2i * csrr[iflag - 1];
06191 str = -c2i * cidi;
06192 c2i = c2r * cidi;
06193 c2r = str;
06194
06195 }
06196 if (nd <= 2) {
06197 goto L110;
06198 }
06199 raz = 1. / myzabs_(zr, zi);
06200 str = *zr * raz;
06201 sti = -(*zi) * raz;
06202 rzr = (str + str) * raz;
06203 rzi = (sti + sti) * raz;
06204 bry[1] = 1. / bry[0];
06205 bry[2] = d1mach_(&c__2);
06206 s1r = cyr[0];
06207 s1i = cyi[0];
06208 s2r = cyr[1];
06209 s2i = cyi[1];
06210 c1r = csrr[iflag - 1];
06211 ascle = bry[iflag - 1];
06212 k = nd - 2;
06213 fn = (doublereal) ((real) k);
06214 i__1 = nd;
06215 for (i__ = 3; i__ <= i__1; ++i__) {
06216 c2r = s2r;
06217 c2i = s2i;
06218 s2r = s1r + (*fnu + fn) * (rzr * c2r - rzi * c2i);
06219 s2i = s1i + (*fnu + fn) * (rzr * c2i + rzi * c2r);
06220 s1r = c2r;
06221 s1i = c2i;
06222 c2r = s2r * c1r;
06223 c2i = s2i * c1r;
06224 yr[k] = c2r;
06225 yi[k] = c2i;
06226 --k;
06227 fn += -1.;
06228 if (iflag >= 3) {
06229 goto L100;
06230 }
06231 str = abs(c2r);
06232 sti = abs(c2i);
06233 c2m = max(str,sti);
06234 if (c2m <= ascle) {
06235 goto L100;
06236 }
06237 ++iflag;
06238 ascle = bry[iflag - 1];
06239 s1r *= c1r;
06240 s1i *= c1r;
06241 s2r = c2r;
06242 s2i = c2i;
06243 s1r *= cssr[iflag - 1];
06244 s1i *= cssr[iflag - 1];
06245 s2r *= cssr[iflag - 1];
06246 s2i *= cssr[iflag - 1];
06247 c1r = csrr[iflag - 1];
06248 L100:
06249 ;
06250 }
06251 L110:
06252 return 0;
06253 L120:
06254 if (rs1 > 0.) {
06255 goto L140;
06256 }
06257
06258
06259
06260 yr[nd] = zeror;
06261 yi[nd] = zeroi;
06262 ++(*nz);
06263 --nd;
06264 if (nd == 0) {
06265 goto L110;
06266 }
06267 zuoik_(zr, zi, fnu, kode, &c__1, &nd, &yr[1], &yi[1], &nuf, tol, elim,
06268 alim);
06269 if (nuf < 0) {
06270 goto L140;
06271 }
06272 nd -= nuf;
06273 *nz += nuf;
06274 if (nd == 0) {
06275 goto L110;
06276 }
06277 fn = *fnu + (doublereal) ((real) (nd - 1));
06278 if (fn < *fnul) {
06279 goto L130;
06280 }
06281
06282
06283
06284
06285
06286
06287
06288
06289
06290 in = inu + nd - 1;
06291 in = in % 4 + 1;
06292 c2r = car * cipr[in - 1] - sar * cipi[in - 1];
06293 c2i = car * cipi[in - 1] + sar * cipr[in - 1];
06294 if (*zi <= 0.) {
06295 c2i = -c2i;
06296 }
06297 goto L40;
06298 L130:
06299 *nlast = nd;
06300 return 0;
06301 L140:
06302 *nz = -1;
06303 return 0;
06304 L150:
06305 if (rs1 > 0.) {
06306 goto L140;
06307 }
06308 *nz = *n;
06309 i__1 = *n;
06310 for (i__ = 1; i__ <= i__1; ++i__) {
06311 yr[i__] = zeror;
06312 yi[i__] = zeroi;
06313
06314 }
06315 return 0;
06316 }
06317
06318 int zunik_(doublereal *zrr, doublereal *zri, doublereal *fnu, integer *ikflg, integer *ipmtr, doublereal *tol, integer *init, doublereal *phir, doublereal *phii, doublereal *zeta1r, doublereal *zeta1i, doublereal *zeta2r, doublereal *zeta2i, doublereal *sumr, doublereal *sumi, doublereal *cwrkr, doublereal *cwrki)
06319 {
06320
06321
06322 static doublereal zeror = 0.;
06323 static doublereal zeroi = 0.;
06324 static doublereal coner = 1.;
06325 static doublereal conei = 0.;
06326 static doublereal con[2] = { .398942280401432678,1.25331413731550025 };
06327 static doublereal c__[120] = { 1.,-.208333333333333333,.125,
06328 .334201388888888889,-.401041666666666667,.0703125,
06329 -1.02581259645061728,1.84646267361111111,-.8912109375,.0732421875,
06330 4.66958442342624743,-11.2070026162229938,8.78912353515625,
06331 -2.3640869140625,.112152099609375,-28.2120725582002449,
06332 84.6362176746007346,-91.8182415432400174,42.5349987453884549,
06333 -7.3687943594796317,.227108001708984375,212.570130039217123,
06334 -765.252468141181642,1059.99045252799988,-699.579627376132541,
06335 218.19051174421159,-26.4914304869515555,.572501420974731445,
06336 -1919.457662318407,8061.72218173730938,-13586.5500064341374,
06337 11655.3933368645332,-5305.64697861340311,1200.90291321635246,
06338 -108.090919788394656,1.7277275025844574,20204.2913309661486,
06339 -96980.5983886375135,192547.001232531532,-203400.177280415534,
06340 122200.46498301746,-41192.6549688975513,7109.51430248936372,
06341 -493.915304773088012,6.07404200127348304,-242919.187900551333,
06342 1311763.6146629772,-2998015.91853810675,3763271.297656404,
06343 -2813563.22658653411,1268365.27332162478,-331645.172484563578,
06344 45218.7689813627263,-2499.83048181120962,24.3805296995560639,
06345 3284469.85307203782,-19706819.1184322269,50952602.4926646422,
06346 -74105148.2115326577,66344512.2747290267,-37567176.6607633513,
06347 13288767.1664218183,-2785618.12808645469,308186.404612662398,
06348 -13886.0897537170405,110.017140269246738,-49329253.664509962,
06349 325573074.185765749,-939462359.681578403,1553596899.57058006,
06350 -1621080552.10833708,1106842816.82301447,-495889784.275030309,
06351 142062907.797533095,-24474062.7257387285,2243768.17792244943,
06352 -84005.4336030240853,551.335896122020586,814789096.118312115,
06353 -5866481492.05184723,18688207509.2958249,-34632043388.1587779,
06354 41280185579.753974,-33026599749.8007231,17954213731.1556001,
06355 -6563293792.61928433,1559279864.87925751,-225105661.889415278,
06356 17395107.5539781645,-549842.327572288687,3038.09051092238427,
06357 -14679261247.6956167,114498237732.02581,-399096175224.466498,
06358 819218669548.577329,-1098375156081.22331,1008158106865.38209,
06359 -645364869245.376503,287900649906.150589,-87867072178.0232657,
06360 17634730606.8349694,-2167164983.22379509,143157876.718888981,
06361 -3871833.44257261262,18257.7554742931747,286464035717.679043,
06362 -2406297900028.50396,9109341185239.89896,-20516899410934.4374,
06363 30565125519935.3206,-31667088584785.1584,23348364044581.8409,
06364 -12320491305598.2872,4612725780849.13197,-1196552880196.1816,
06365 205914503232.410016,-21822927757.5292237,1247009293.51271032,
06366 -29188388.1222208134,118838.426256783253 };
06367
06368
06369 integer i__1;
06370 doublereal d__1, d__2;
06371
06372
06373
06374
06375 static integer idum;
06376 static doublereal test;
06377 static integer i__, j, k, l;
06378 static doublereal crfni, crfnr;
06379 static doublereal ac, si, ti, sr, tr, t2i, t2r, rfn, sri, sti, zni, srr,
06380 str, znr;
06381
06382
06383
06384
06385
06386
06387
06388
06389
06390
06391
06392
06393
06394
06395
06396
06397
06398
06399
06400
06401
06402
06403
06404
06405 --cwrki;
06406 --cwrkr;
06407
06408
06409
06410 if (*init != 0) {
06411 goto L40;
06412 }
06413
06414
06415
06416 rfn = 1. / *fnu;
06417
06418
06419
06420 test = d1mach_(&c__1) * 1e3;
06421 ac = *fnu * test;
06422 if (abs(*zrr) > ac || abs(*zri) > ac) {
06423 goto L15;
06424 }
06425 *zeta1r = (d__1 = log(test), abs(d__1)) * 2. + *fnu;
06426 *zeta1i = 0.;
06427 *zeta2r = *fnu;
06428 *zeta2i = 0.;
06429 *phir = 1.;
06430 *phii = 0.;
06431 return 0;
06432 L15:
06433 tr = *zrr * rfn;
06434 ti = *zri * rfn;
06435 sr = coner + (tr * tr - ti * ti);
06436 si = conei + (tr * ti + ti * tr);
06437 zsqrt_(&sr, &si, &srr, &sri);
06438 str = coner + srr;
06439 sti = conei + sri;
06440 zdiv_(&str, &sti, &tr, &ti, &znr, &zni);
06441 zlog_(&znr, &zni, &str, &sti, &idum);
06442 *zeta1r = *fnu * str;
06443 *zeta1i = *fnu * sti;
06444 *zeta2r = *fnu * srr;
06445 *zeta2i = *fnu * sri;
06446 zdiv_(&coner, &conei, &srr, &sri, &tr, &ti);
06447 srr = tr * rfn;
06448 sri = ti * rfn;
06449 zsqrt_(&srr, &sri, &cwrkr[16], &cwrki[16]);
06450 *phir = cwrkr[16] * con[*ikflg - 1];
06451 *phii = cwrki[16] * con[*ikflg - 1];
06452 if (*ipmtr != 0) {
06453 return 0;
06454 }
06455 zdiv_(&coner, &conei, &sr, &si, &t2r, &t2i);
06456 cwrkr[1] = coner;
06457 cwrki[1] = conei;
06458 crfnr = coner;
06459 crfni = conei;
06460 ac = 1.;
06461 l = 1;
06462 for (k = 2; k <= 15; ++k) {
06463 sr = zeror;
06464 si = zeroi;
06465 i__1 = k;
06466 for (j = 1; j <= i__1; ++j) {
06467 ++l;
06468 str = sr * t2r - si * t2i + c__[l - 1];
06469 si = sr * t2i + si * t2r;
06470 sr = str;
06471
06472 }
06473 str = crfnr * srr - crfni * sri;
06474 crfni = crfnr * sri + crfni * srr;
06475 crfnr = str;
06476 cwrkr[k] = crfnr * sr - crfni * si;
06477 cwrki[k] = crfnr * si + crfni * sr;
06478 ac *= rfn;
06479 test = (d__1 = cwrkr[k], abs(d__1)) + (d__2 = cwrki[k], abs(d__2));
06480 if (ac < *tol && test < *tol) {
06481 goto L30;
06482 }
06483
06484 }
06485 k = 15;
06486 L30:
06487 *init = k;
06488 L40:
06489 if (*ikflg == 2) {
06490 goto L60;
06491 }
06492
06493
06494
06495 sr = zeror;
06496 si = zeroi;
06497 i__1 = *init;
06498 for (i__ = 1; i__ <= i__1; ++i__) {
06499 sr += cwrkr[i__];
06500 si += cwrki[i__];
06501
06502 }
06503 *sumr = sr;
06504 *sumi = si;
06505 *phir = cwrkr[16] * con[0];
06506 *phii = cwrki[16] * con[0];
06507 return 0;
06508 L60:
06509
06510
06511
06512 sr = zeror;
06513 si = zeroi;
06514 tr = coner;
06515 i__1 = *init;
06516 for (i__ = 1; i__ <= i__1; ++i__) {
06517 sr += tr * cwrkr[i__];
06518 si += tr * cwrki[i__];
06519 tr = -tr;
06520
06521 }
06522 *sumr = sr;
06523 *sumi = si;
06524 *phir = cwrkr[16] * con[1];
06525 *phii = cwrki[16] * con[1];
06526 return 0;
06527 }
06528
06529 int zunk1_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *mr, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *tol, doublereal *elim, doublereal *alim)
06530 {
06531
06532
06533 static doublereal zeror = 0.;
06534 static doublereal zeroi = 0.;
06535 static doublereal coner = 1.;
06536 static doublereal pi = 3.14159265358979324;
06537
06538
06539 integer i__1;
06540
06541
06542
06543
06544 static doublereal aphi, cscl, phii[2], crsc, phir[2];
06545 static integer init[2];
06546 static doublereal csrr[3], cssr[3], rast, sumi[2], razr;
06547 static doublereal sumr[2];
06548 static integer i__, j, k, m, iflag, kflag;
06549 static doublereal ascle;
06550 static integer kdflg;
06551 static doublereal phidi;
06552 static integer ipard;
06553 static doublereal csgni, phidr;
06554 static integer initd;
06555 static doublereal cspni, cwrki[48] , sumdi;
06556 static doublereal cspnr, cwrkr[48] , sumdr;
06557 static doublereal zeta1i[2], zeta2i[2], zet1di, zet2di, zeta1r[2], zeta2r[
06558 2], zet1dr, zet2dr;
06559 static integer ib, ic;
06560 static doublereal fn;
06561 static integer il, kk, nw;
06562 static doublereal c1i, c2i, c2m;
06563 static doublereal c1r, c2r, s1i, s2i, rs1, s1r, s2r, ang, asc, cki, fnf;
06564 static integer ifn;
06565 static doublereal ckr;
06566 static integer iuf;
06567 static doublereal cyi[2], fmr, csr, sgn;
06568 static integer inu;
06569 static doublereal bry[3], cyr[2], sti, rzi, zri, str, rzr, zrr;
06570
06571
06572
06573
06574
06575
06576
06577
06578
06579
06580
06581
06582
06583
06584
06585 --yi;
06586 --yr;
06587
06588
06589
06590 kdflg = 1;
06591 *nz = 0;
06592
06593
06594
06595
06596 cscl = 1. / *tol;
06597 crsc = *tol;
06598 cssr[0] = cscl;
06599 cssr[1] = coner;
06600 cssr[2] = crsc;
06601 csrr[0] = crsc;
06602 csrr[1] = coner;
06603 csrr[2] = cscl;
06604 bry[0] = d1mach_(&c__1) * 1e3 / *tol;
06605 bry[1] = 1. / bry[0];
06606 bry[2] = d1mach_(&c__2);
06607 zrr = *zr;
06608 zri = *zi;
06609 if (*zr >= 0.) {
06610 goto L10;
06611 }
06612 zrr = -(*zr);
06613 zri = -(*zi);
06614 L10:
06615 j = 2;
06616 i__1 = *n;
06617 for (i__ = 1; i__ <= i__1; ++i__) {
06618
06619
06620
06621 j = 3 - j;
06622 fn = *fnu + (doublereal) ((real) (i__ - 1));
06623 init[j - 1] = 0;
06624 zunik_(&zrr, &zri, &fn, &c__2, &c__0, tol, &init[j - 1], &phir[j - 1],
06625 &phii[j - 1], &zeta1r[j - 1], &zeta1i[j - 1], &zeta2r[j - 1],
06626 &zeta2i[j - 1], &sumr[j - 1], &sumi[j - 1], &cwrkr[(j << 4)
06627 - 16], &cwrki[(j << 4) - 16]);
06628 if (*kode == 1) {
06629 goto L20;
06630 }
06631 str = zrr + zeta2r[j - 1];
06632 sti = zri + zeta2i[j - 1];
06633 rast = fn / myzabs_(&str, &sti);
06634 str = str * rast * rast;
06635 sti = -sti * rast * rast;
06636 s1r = zeta1r[j - 1] - str;
06637 s1i = zeta1i[j - 1] - sti;
06638 goto L30;
06639 L20:
06640 s1r = zeta1r[j - 1] - zeta2r[j - 1];
06641 s1i = zeta1i[j - 1] - zeta2i[j - 1];
06642 L30:
06643 rs1 = s1r;
06644
06645
06646
06647 if (abs(rs1) > *elim) {
06648 goto L60;
06649 }
06650 if (kdflg == 1) {
06651 kflag = 2;
06652 }
06653 if (abs(rs1) < *alim) {
06654 goto L40;
06655 }
06656
06657
06658
06659 aphi = myzabs_(&phir[j - 1], &phii[j - 1]);
06660 rs1 += log(aphi);
06661 if (abs(rs1) > *elim) {
06662 goto L60;
06663 }
06664 if (kdflg == 1) {
06665 kflag = 1;
06666 }
06667 if (rs1 < 0.) {
06668 goto L40;
06669 }
06670 if (kdflg == 1) {
06671 kflag = 3;
06672 }
06673 L40:
06674
06675
06676
06677
06678 s2r = phir[j - 1] * sumr[j - 1] - phii[j - 1] * sumi[j - 1];
06679 s2i = phir[j - 1] * sumi[j - 1] + phii[j - 1] * sumr[j - 1];
06680 str = exp(s1r) * cssr[kflag - 1];
06681 s1r = str * cos(s1i);
06682 s1i = str * sin(s1i);
06683 str = s2r * s1r - s2i * s1i;
06684 s2i = s1r * s2i + s2r * s1i;
06685 s2r = str;
06686 if (kflag != 1) {
06687 goto L50;
06688 }
06689 zuchk_(&s2r, &s2i, &nw, bry, tol);
06690 if (nw != 0) {
06691 goto L60;
06692 }
06693 L50:
06694 cyr[kdflg - 1] = s2r;
06695 cyi[kdflg - 1] = s2i;
06696 yr[i__] = s2r * csrr[kflag - 1];
06697 yi[i__] = s2i * csrr[kflag - 1];
06698 if (kdflg == 2) {
06699 goto L75;
06700 }
06701 kdflg = 2;
06702 goto L70;
06703 L60:
06704 if (rs1 > 0.) {
06705 goto L300;
06706 }
06707
06708
06709
06710 if (*zr < 0.) {
06711 goto L300;
06712 }
06713 kdflg = 1;
06714 yr[i__] = zeror;
06715 yi[i__] = zeroi;
06716 ++(*nz);
06717 if (i__ == 1) {
06718 goto L70;
06719 }
06720 if (yr[i__ - 1] == zeror && yi[i__ - 1] == zeroi) {
06721 goto L70;
06722 }
06723 yr[i__ - 1] = zeror;
06724 yi[i__ - 1] = zeroi;
06725 ++(*nz);
06726 L70:
06727 ;
06728 }
06729 i__ = *n;
06730 L75:
06731 razr = 1. / myzabs_(&zrr, &zri);
06732 str = zrr * razr;
06733 sti = -zri * razr;
06734 rzr = (str + str) * razr;
06735 rzi = (sti + sti) * razr;
06736 ckr = fn * rzr;
06737 cki = fn * rzi;
06738 ib = i__ + 1;
06739 if (*n < ib) {
06740 goto L160;
06741 }
06742
06743
06744
06745
06746 fn = *fnu + (doublereal) ((real) (*n - 1));
06747 ipard = 1;
06748 if (*mr != 0) {
06749 ipard = 0;
06750 }
06751 initd = 0;
06752 zunik_(&zrr, &zri, &fn, &c__2, &ipard, tol, &initd, &phidr, &phidi, &
06753 zet1dr, &zet1di, &zet2dr, &zet2di, &sumdr, &sumdi, &cwrkr[32], &
06754 cwrki[32]);
06755 if (*kode == 1) {
06756 goto L80;
06757 }
06758 str = zrr + zet2dr;
06759 sti = zri + zet2di;
06760 rast = fn / myzabs_(&str, &sti);
06761 str = str * rast * rast;
06762 sti = -sti * rast * rast;
06763 s1r = zet1dr - str;
06764 s1i = zet1di - sti;
06765 goto L90;
06766 L80:
06767 s1r = zet1dr - zet2dr;
06768 s1i = zet1di - zet2di;
06769 L90:
06770 rs1 = s1r;
06771 if (abs(rs1) > *elim) {
06772 goto L95;
06773 }
06774 if (abs(rs1) < *alim) {
06775 goto L100;
06776 }
06777
06778
06779
06780 aphi = myzabs_(&phidr, &phidi);
06781 rs1 += log(aphi);
06782 if (abs(rs1) < *elim) {
06783 goto L100;
06784 }
06785 L95:
06786 if (abs(rs1) > 0.) {
06787 goto L300;
06788 }
06789
06790
06791
06792 if (*zr < 0.) {
06793 goto L300;
06794 }
06795 *nz = *n;
06796 i__1 = *n;
06797 for (i__ = 1; i__ <= i__1; ++i__) {
06798 yr[i__] = zeror;
06799 yi[i__] = zeroi;
06800
06801 }
06802 return 0;
06803
06804
06805
06806 L100:
06807 s1r = cyr[0];
06808 s1i = cyi[0];
06809 s2r = cyr[1];
06810 s2i = cyi[1];
06811 c1r = csrr[kflag - 1];
06812 ascle = bry[kflag - 1];
06813 i__1 = *n;
06814 for (i__ = ib; i__ <= i__1; ++i__) {
06815 c2r = s2r;
06816 c2i = s2i;
06817 s2r = ckr * c2r - cki * c2i + s1r;
06818 s2i = ckr * c2i + cki * c2r + s1i;
06819 s1r = c2r;
06820 s1i = c2i;
06821 ckr += rzr;
06822 cki += rzi;
06823 c2r = s2r * c1r;
06824 c2i = s2i * c1r;
06825 yr[i__] = c2r;
06826 yi[i__] = c2i;
06827 if (kflag >= 3) {
06828 goto L120;
06829 }
06830 str = abs(c2r);
06831 sti = abs(c2i);
06832 c2m = max(str,sti);
06833 if (c2m <= ascle) {
06834 goto L120;
06835 }
06836 ++kflag;
06837 ascle = bry[kflag - 1];
06838 s1r *= c1r;
06839 s1i *= c1r;
06840 s2r = c2r;
06841 s2i = c2i;
06842 s1r *= cssr[kflag - 1];
06843 s1i *= cssr[kflag - 1];
06844 s2r *= cssr[kflag - 1];
06845 s2i *= cssr[kflag - 1];
06846 c1r = csrr[kflag - 1];
06847 L120:
06848 ;
06849 }
06850 L160:
06851 if (*mr == 0) {
06852 return 0;
06853 }
06854
06855
06856
06857 *nz = 0;
06858 fmr = (doublereal) ((real) (*mr));
06859 sgn = -d_sign(&pi, &fmr);
06860
06861
06862
06863 csgni = sgn;
06864 inu = (integer) ((real) (*fnu));
06865 fnf = *fnu - (doublereal) ((real) inu);
06866 ifn = inu + *n - 1;
06867 ang = fnf * sgn;
06868 cspnr = cos(ang);
06869 cspni = sin(ang);
06870 if (ifn % 2 == 0) {
06871 goto L170;
06872 }
06873 cspnr = -cspnr;
06874 cspni = -cspni;
06875 L170:
06876 asc = bry[0];
06877 iuf = 0;
06878 kk = *n;
06879 kdflg = 1;
06880 --ib;
06881 ic = ib - 1;
06882 i__1 = *n;
06883 for (k = 1; k <= i__1; ++k) {
06884 fn = *fnu + (doublereal) ((real) (kk - 1));
06885
06886
06887
06888
06889 m = 3;
06890 if (*n > 2) {
06891 goto L175;
06892 }
06893 L172:
06894 initd = init[j - 1];
06895 phidr = phir[j - 1];
06896 phidi = phii[j - 1];
06897 zet1dr = zeta1r[j - 1];
06898 zet1di = zeta1i[j - 1];
06899 zet2dr = zeta2r[j - 1];
06900 zet2di = zeta2i[j - 1];
06901 sumdr = sumr[j - 1];
06902 sumdi = sumi[j - 1];
06903 m = j;
06904 j = 3 - j;
06905 goto L180;
06906 L175:
06907 if (kk == *n && ib < *n) {
06908 goto L180;
06909 }
06910 if (kk == ib || kk == ic) {
06911 goto L172;
06912 }
06913 initd = 0;
06914 L180:
06915 zunik_(&zrr, &zri, &fn, &c__1, &c__0, tol, &initd, &phidr, &phidi, &
06916 zet1dr, &zet1di, &zet2dr, &zet2di, &sumdr, &sumdi, &cwrkr[(m
06917 << 4) - 16], &cwrki[(m << 4) - 16]);
06918 if (*kode == 1) {
06919 goto L200;
06920 }
06921 str = zrr + zet2dr;
06922 sti = zri + zet2di;
06923 rast = fn / myzabs_(&str, &sti);
06924 str = str * rast * rast;
06925 sti = -sti * rast * rast;
06926 s1r = -zet1dr + str;
06927 s1i = -zet1di + sti;
06928 goto L210;
06929 L200:
06930 s1r = -zet1dr + zet2dr;
06931 s1i = -zet1di + zet2di;
06932 L210:
06933
06934
06935
06936 rs1 = s1r;
06937 if (abs(rs1) > *elim) {
06938 goto L260;
06939 }
06940 if (kdflg == 1) {
06941 iflag = 2;
06942 }
06943 if (abs(rs1) < *alim) {
06944 goto L220;
06945 }
06946
06947
06948
06949 aphi = myzabs_(&phidr, &phidi);
06950 rs1 += log(aphi);
06951 if (abs(rs1) > *elim) {
06952 goto L260;
06953 }
06954 if (kdflg == 1) {
06955 iflag = 1;
06956 }
06957 if (rs1 < 0.) {
06958 goto L220;
06959 }
06960 if (kdflg == 1) {
06961 iflag = 3;
06962 }
06963 L220:
06964 str = phidr * sumdr - phidi * sumdi;
06965 sti = phidr * sumdi + phidi * sumdr;
06966 s2r = -csgni * sti;
06967 s2i = csgni * str;
06968 str = exp(s1r) * cssr[iflag - 1];
06969 s1r = str * cos(s1i);
06970 s1i = str * sin(s1i);
06971 str = s2r * s1r - s2i * s1i;
06972 s2i = s2r * s1i + s2i * s1r;
06973 s2r = str;
06974 if (iflag != 1) {
06975 goto L230;
06976 }
06977 zuchk_(&s2r, &s2i, &nw, bry, tol);
06978 if (nw == 0) {
06979 goto L230;
06980 }
06981 s2r = zeror;
06982 s2i = zeroi;
06983 L230:
06984 cyr[kdflg - 1] = s2r;
06985 cyi[kdflg - 1] = s2i;
06986 c2r = s2r;
06987 c2i = s2i;
06988 s2r *= csrr[iflag - 1];
06989 s2i *= csrr[iflag - 1];
06990
06991
06992
06993 s1r = yr[kk];
06994 s1i = yi[kk];
06995 if (*kode == 1) {
06996 goto L250;
06997 }
06998 zs1s2_(&zrr, &zri, &s1r, &s1i, &s2r, &s2i, &nw, &asc, alim, &iuf);
06999 *nz += nw;
07000 L250:
07001 yr[kk] = s1r * cspnr - s1i * cspni + s2r;
07002 yi[kk] = cspnr * s1i + cspni * s1r + s2i;
07003 --kk;
07004 cspnr = -cspnr;
07005 cspni = -cspni;
07006 if (c2r != 0. || c2i != 0.) {
07007 goto L255;
07008 }
07009 kdflg = 1;
07010 goto L270;
07011 L255:
07012 if (kdflg == 2) {
07013 goto L275;
07014 }
07015 kdflg = 2;
07016 goto L270;
07017 L260:
07018 if (rs1 > 0.) {
07019 goto L300;
07020 }
07021 s2r = zeror;
07022 s2i = zeroi;
07023 goto L230;
07024 L270:
07025 ;
07026 }
07027 k = *n;
07028 L275:
07029 il = *n - k;
07030 if (il == 0) {
07031 return 0;
07032 }
07033
07034
07035
07036
07037
07038 s1r = cyr[0];
07039 s1i = cyi[0];
07040 s2r = cyr[1];
07041 s2i = cyi[1];
07042 csr = csrr[iflag - 1];
07043 ascle = bry[iflag - 1];
07044 fn = (doublereal) ((real) (inu + il));
07045 i__1 = il;
07046 for (i__ = 1; i__ <= i__1; ++i__) {
07047 c2r = s2r;
07048 c2i = s2i;
07049 s2r = s1r + (fn + fnf) * (rzr * c2r - rzi * c2i);
07050 s2i = s1i + (fn + fnf) * (rzr * c2i + rzi * c2r);
07051 s1r = c2r;
07052 s1i = c2i;
07053 fn += -1.;
07054 c2r = s2r * csr;
07055 c2i = s2i * csr;
07056 ckr = c2r;
07057 cki = c2i;
07058 c1r = yr[kk];
07059 c1i = yi[kk];
07060 if (*kode == 1) {
07061 goto L280;
07062 }
07063 zs1s2_(&zrr, &zri, &c1r, &c1i, &c2r, &c2i, &nw, &asc, alim, &iuf);
07064 *nz += nw;
07065 L280:
07066 yr[kk] = c1r * cspnr - c1i * cspni + c2r;
07067 yi[kk] = c1r * cspni + c1i * cspnr + c2i;
07068 --kk;
07069 cspnr = -cspnr;
07070 cspni = -cspni;
07071 if (iflag >= 3) {
07072 goto L290;
07073 }
07074 c2r = abs(ckr);
07075 c2i = abs(cki);
07076 c2m = max(c2r,c2i);
07077 if (c2m <= ascle) {
07078 goto L290;
07079 }
07080 ++iflag;
07081 ascle = bry[iflag - 1];
07082 s1r *= csr;
07083 s1i *= csr;
07084 s2r = ckr;
07085 s2i = cki;
07086 s1r *= cssr[iflag - 1];
07087 s1i *= cssr[iflag - 1];
07088 s2r *= cssr[iflag - 1];
07089 s2i *= cssr[iflag - 1];
07090 csr = csrr[iflag - 1];
07091 L290:
07092 ;
07093 }
07094 return 0;
07095 L300:
07096 *nz = -1;
07097 return 0;
07098 }
07099
07100 int zunk2_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *mr, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *tol, doublereal *elim, doublereal *alim)
07101 {
07102
07103
07104 static doublereal zeror = 0.;
07105 static doublereal zeroi = 0.;
07106 static doublereal coner = 1.;
07107 static doublereal cr1r = 1.;
07108 static doublereal cr1i = 1.73205080756887729;
07109 static doublereal cr2r = -.5;
07110 static doublereal cr2i = -.866025403784438647;
07111 static doublereal hpi = 1.57079632679489662;
07112 static doublereal pi = 3.14159265358979324;
07113 static doublereal aic = 1.26551212348464539;
07114 static doublereal cipr[4] = { 1.,0.,-1.,0. };
07115 static doublereal cipi[4] = { 0.,-1.,0.,1. };
07116
07117
07118 integer i__1;
07119
07120
07121
07122
07123 static doublereal daii, aarg;
07124 static integer ndai;
07125 static doublereal dair, aphi, argi[2], cscl, phii[2], crsc, argr[2];
07126 static integer idum;
07127 static doublereal phir[2], csrr[3], cssr[3], rast, razr;
07128 static integer i__, k, j, iflag, kflag;
07129 static doublereal argdi, ascle;
07130 static integer kdflg;
07131 static doublereal phidi, argdr;
07132 static integer ipard;
07133 static doublereal csgni, phidr, cspni, asumi[2], bsumi[2];
07134 static doublereal cspnr, asumr[2], bsumr[2];
07135 static doublereal zeta1i[2], zeta2i[2], zet1di, zet2di, zeta1r[2], zeta2r[
07136 2], zet1dr, zet2dr;
07137 static integer ib, ic;
07138 static doublereal fn;
07139 static integer il, kk, in, nw;
07140 static doublereal asumdi, bsumdi, yy, asumdr, bsumdr, c1i, c2i, c2m;
07141 static doublereal c1r, c2r, s1i, s2i, rs1, s1r, s2r, aii, ang, asc, car,
07142 cki, fnf;
07143 static integer nai;
07144 static doublereal air;
07145 static integer ifn;
07146 static doublereal csi, ckr;
07147 static integer iuf;
07148 static doublereal cyi[2], fmr, sar, csr, sgn, zbi;
07149 static integer inu;
07150 static doublereal bry[3], cyr[2], pti, sti, zbr, zni, rzi, ptr, zri, str,
07151 znr, rzr, zrr;
07152
07153
07154
07155
07156
07157
07158
07159
07160
07161
07162
07163
07164
07165
07166
07167
07168
07169
07170
07171 --yi;
07172 --yr;
07173
07174
07175
07176 kdflg = 1;
07177 *nz = 0;
07178
07179
07180
07181
07182 cscl = 1. / *tol;
07183 crsc = *tol;
07184 cssr[0] = cscl;
07185 cssr[1] = coner;
07186 cssr[2] = crsc;
07187 csrr[0] = crsc;
07188 csrr[1] = coner;
07189 csrr[2] = cscl;
07190 bry[0] = d1mach_(&c__1) * 1e3 / *tol;
07191 bry[1] = 1. / bry[0];
07192 bry[2] = d1mach_(&c__2);
07193 zrr = *zr;
07194 zri = *zi;
07195 if (*zr >= 0.) {
07196 goto L10;
07197 }
07198 zrr = -(*zr);
07199 zri = -(*zi);
07200 L10:
07201 yy = zri;
07202 znr = zri;
07203 zni = -zrr;
07204 zbr = zrr;
07205 zbi = zri;
07206 inu = (integer) ((real) (*fnu));
07207 fnf = *fnu - (doublereal) ((real) inu);
07208 ang = -hpi * fnf;
07209 car = cos(ang);
07210 sar = sin(ang);
07211 c2r = hpi * sar;
07212 c2i = -hpi * car;
07213 kk = inu % 4 + 1;
07214 str = c2r * cipr[kk - 1] - c2i * cipi[kk - 1];
07215 sti = c2r * cipi[kk - 1] + c2i * cipr[kk - 1];
07216 csr = cr1r * str - cr1i * sti;
07217 csi = cr1r * sti + cr1i * str;
07218 if (yy > 0.) {
07219 goto L20;
07220 }
07221 znr = -znr;
07222 zbi = -zbi;
07223 L20:
07224
07225
07226
07227
07228
07229 j = 2;
07230 i__1 = *n;
07231 for (i__ = 1; i__ <= i__1; ++i__) {
07232
07233
07234
07235 j = 3 - j;
07236 fn = *fnu + (doublereal) ((real) (i__ - 1));
07237 zunhj_(&znr, &zni, &fn, &c__0, tol, &phir[j - 1], &phii[j - 1], &argr[
07238 j - 1], &argi[j - 1], &zeta1r[j - 1], &zeta1i[j - 1], &zeta2r[
07239 j - 1], &zeta2i[j - 1], &asumr[j - 1], &asumi[j - 1], &bsumr[
07240 j - 1], &bsumi[j - 1]);
07241 if (*kode == 1) {
07242 goto L30;
07243 }
07244 str = zbr + zeta2r[j - 1];
07245 sti = zbi + zeta2i[j - 1];
07246 rast = fn / myzabs_(&str, &sti);
07247 str = str * rast * rast;
07248 sti = -sti * rast * rast;
07249 s1r = zeta1r[j - 1] - str;
07250 s1i = zeta1i[j - 1] - sti;
07251 goto L40;
07252 L30:
07253 s1r = zeta1r[j - 1] - zeta2r[j - 1];
07254 s1i = zeta1i[j - 1] - zeta2i[j - 1];
07255 L40:
07256
07257
07258
07259 rs1 = s1r;
07260 if (abs(rs1) > *elim) {
07261 goto L70;
07262 }
07263 if (kdflg == 1) {
07264 kflag = 2;
07265 }
07266 if (abs(rs1) < *alim) {
07267 goto L50;
07268 }
07269
07270
07271
07272 aphi = myzabs_(&phir[j - 1], &phii[j - 1]);
07273 aarg = myzabs_(&argr[j - 1], &argi[j - 1]);
07274 rs1 = rs1 + log(aphi) - log(aarg) * .25 - aic;
07275 if (abs(rs1) > *elim) {
07276 goto L70;
07277 }
07278 if (kdflg == 1) {
07279 kflag = 1;
07280 }
07281 if (rs1 < 0.) {
07282 goto L50;
07283 }
07284 if (kdflg == 1) {
07285 kflag = 3;
07286 }
07287 L50:
07288
07289
07290
07291
07292 c2r = argr[j - 1] * cr2r - argi[j - 1] * cr2i;
07293 c2i = argr[j - 1] * cr2i + argi[j - 1] * cr2r;
07294 zairy_(&c2r, &c2i, &c__0, &c__2, &air, &aii, &nai, &idum);
07295 zairy_(&c2r, &c2i, &c__1, &c__2, &dair, &daii, &ndai, &idum);
07296 str = dair * bsumr[j - 1] - daii * bsumi[j - 1];
07297 sti = dair * bsumi[j - 1] + daii * bsumr[j - 1];
07298 ptr = str * cr2r - sti * cr2i;
07299 pti = str * cr2i + sti * cr2r;
07300 str = ptr + (air * asumr[j - 1] - aii * asumi[j - 1]);
07301 sti = pti + (air * asumi[j - 1] + aii * asumr[j - 1]);
07302 ptr = str * phir[j - 1] - sti * phii[j - 1];
07303 pti = str * phii[j - 1] + sti * phir[j - 1];
07304 s2r = ptr * csr - pti * csi;
07305 s2i = ptr * csi + pti * csr;
07306 str = exp(s1r) * cssr[kflag - 1];
07307 s1r = str * cos(s1i);
07308 s1i = str * sin(s1i);
07309 str = s2r * s1r - s2i * s1i;
07310 s2i = s1r * s2i + s2r * s1i;
07311 s2r = str;
07312 if (kflag != 1) {
07313 goto L60;
07314 }
07315 zuchk_(&s2r, &s2i, &nw, bry, tol);
07316 if (nw != 0) {
07317 goto L70;
07318 }
07319 L60:
07320 if (yy <= 0.) {
07321 s2i = -s2i;
07322 }
07323 cyr[kdflg - 1] = s2r;
07324 cyi[kdflg - 1] = s2i;
07325 yr[i__] = s2r * csrr[kflag - 1];
07326 yi[i__] = s2i * csrr[kflag - 1];
07327 str = csi;
07328 csi = -csr;
07329 csr = str;
07330 if (kdflg == 2) {
07331 goto L85;
07332 }
07333 kdflg = 2;
07334 goto L80;
07335 L70:
07336 if (rs1 > 0.) {
07337 goto L320;
07338 }
07339
07340
07341
07342 if (*zr < 0.) {
07343 goto L320;
07344 }
07345 kdflg = 1;
07346 yr[i__] = zeror;
07347 yi[i__] = zeroi;
07348 ++(*nz);
07349 str = csi;
07350 csi = -csr;
07351 csr = str;
07352 if (i__ == 1) {
07353 goto L80;
07354 }
07355 if (yr[i__ - 1] == zeror && yi[i__ - 1] == zeroi) {
07356 goto L80;
07357 }
07358 yr[i__ - 1] = zeror;
07359 yi[i__ - 1] = zeroi;
07360 ++(*nz);
07361 L80:
07362 ;
07363 }
07364 i__ = *n;
07365 L85:
07366 razr = 1. / myzabs_(&zrr, &zri);
07367 str = zrr * razr;
07368 sti = -zri * razr;
07369 rzr = (str + str) * razr;
07370 rzi = (sti + sti) * razr;
07371 ckr = fn * rzr;
07372 cki = fn * rzi;
07373 ib = i__ + 1;
07374 if (*n < ib) {
07375 goto L180;
07376 }
07377
07378
07379
07380
07381 fn = *fnu + (doublereal) ((real) (*n - 1));
07382 ipard = 1;
07383 if (*mr != 0) {
07384 ipard = 0;
07385 }
07386 zunhj_(&znr, &zni, &fn, &ipard, tol, &phidr, &phidi, &argdr, &argdi, &
07387 zet1dr, &zet1di, &zet2dr, &zet2di, &asumdr, &asumdi, &bsumdr, &
07388 bsumdi);
07389 if (*kode == 1) {
07390 goto L90;
07391 }
07392 str = zbr + zet2dr;
07393 sti = zbi + zet2di;
07394 rast = fn / myzabs_(&str, &sti);
07395 str = str * rast * rast;
07396 sti = -sti * rast * rast;
07397 s1r = zet1dr - str;
07398 s1i = zet1di - sti;
07399 goto L100;
07400 L90:
07401 s1r = zet1dr - zet2dr;
07402 s1i = zet1di - zet2di;
07403 L100:
07404 rs1 = s1r;
07405 if (abs(rs1) > *elim) {
07406 goto L105;
07407 }
07408 if (abs(rs1) < *alim) {
07409 goto L120;
07410 }
07411
07412
07413
07414 aphi = myzabs_(&phidr, &phidi);
07415 rs1 += log(aphi);
07416 if (abs(rs1) < *elim) {
07417 goto L120;
07418 }
07419 L105:
07420 if (rs1 > 0.) {
07421 goto L320;
07422 }
07423
07424
07425
07426 if (*zr < 0.) {
07427 goto L320;
07428 }
07429 *nz = *n;
07430 i__1 = *n;
07431 for (i__ = 1; i__ <= i__1; ++i__) {
07432 yr[i__] = zeror;
07433 yi[i__] = zeroi;
07434
07435 }
07436 return 0;
07437 L120:
07438 s1r = cyr[0];
07439 s1i = cyi[0];
07440 s2r = cyr[1];
07441 s2i = cyi[1];
07442 c1r = csrr[kflag - 1];
07443 ascle = bry[kflag - 1];
07444 i__1 = *n;
07445 for (i__ = ib; i__ <= i__1; ++i__) {
07446 c2r = s2r;
07447 c2i = s2i;
07448 s2r = ckr * c2r - cki * c2i + s1r;
07449 s2i = ckr * c2i + cki * c2r + s1i;
07450 s1r = c2r;
07451 s1i = c2i;
07452 ckr += rzr;
07453 cki += rzi;
07454 c2r = s2r * c1r;
07455 c2i = s2i * c1r;
07456 yr[i__] = c2r;
07457 yi[i__] = c2i;
07458 if (kflag >= 3) {
07459 goto L130;
07460 }
07461 str = abs(c2r);
07462 sti = abs(c2i);
07463 c2m = max(str,sti);
07464 if (c2m <= ascle) {
07465 goto L130;
07466 }
07467 ++kflag;
07468 ascle = bry[kflag - 1];
07469 s1r *= c1r;
07470 s1i *= c1r;
07471 s2r = c2r;
07472 s2i = c2i;
07473 s1r *= cssr[kflag - 1];
07474 s1i *= cssr[kflag - 1];
07475 s2r *= cssr[kflag - 1];
07476 s2i *= cssr[kflag - 1];
07477 c1r = csrr[kflag - 1];
07478 L130:
07479 ;
07480 }
07481 L180:
07482 if (*mr == 0) {
07483 return 0;
07484 }
07485
07486
07487
07488 *nz = 0;
07489 fmr = (doublereal) ((real) (*mr));
07490 sgn = -d_sign(&pi, &fmr);
07491
07492
07493
07494 csgni = sgn;
07495 if (yy <= 0.) {
07496 csgni = -csgni;
07497 }
07498 ifn = inu + *n - 1;
07499 ang = fnf * sgn;
07500 cspnr = cos(ang);
07501 cspni = sin(ang);
07502 if (ifn % 2 == 0) {
07503 goto L190;
07504 }
07505 cspnr = -cspnr;
07506 cspni = -cspni;
07507 L190:
07508
07509
07510
07511
07512
07513
07514 csr = sar * csgni;
07515 csi = car * csgni;
07516 in = ifn % 4 + 1;
07517 c2r = cipr[in - 1];
07518 c2i = cipi[in - 1];
07519 str = csr * c2r + csi * c2i;
07520 csi = -csr * c2i + csi * c2r;
07521 csr = str;
07522 asc = bry[0];
07523 iuf = 0;
07524 kk = *n;
07525 kdflg = 1;
07526 --ib;
07527 ic = ib - 1;
07528 i__1 = *n;
07529 for (k = 1; k <= i__1; ++k) {
07530 fn = *fnu + (doublereal) ((real) (kk - 1));
07531
07532
07533
07534
07535 if (*n > 2) {
07536 goto L175;
07537 }
07538 L172:
07539 phidr = phir[j - 1];
07540 phidi = phii[j - 1];
07541 argdr = argr[j - 1];
07542 argdi = argi[j - 1];
07543 zet1dr = zeta1r[j - 1];
07544 zet1di = zeta1i[j - 1];
07545 zet2dr = zeta2r[j - 1];
07546 zet2di = zeta2i[j - 1];
07547 asumdr = asumr[j - 1];
07548 asumdi = asumi[j - 1];
07549 bsumdr = bsumr[j - 1];
07550 bsumdi = bsumi[j - 1];
07551 j = 3 - j;
07552 goto L210;
07553 L175:
07554 if (kk == *n && ib < *n) {
07555 goto L210;
07556 }
07557 if (kk == ib || kk == ic) {
07558 goto L172;
07559 }
07560 zunhj_(&znr, &zni, &fn, &c__0, tol, &phidr, &phidi, &argdr, &argdi, &
07561 zet1dr, &zet1di, &zet2dr, &zet2di, &asumdr, &asumdi, &bsumdr,
07562 &bsumdi);
07563 L210:
07564 if (*kode == 1) {
07565 goto L220;
07566 }
07567 str = zbr + zet2dr;
07568 sti = zbi + zet2di;
07569 rast = fn / myzabs_(&str, &sti);
07570 str = str * rast * rast;
07571 sti = -sti * rast * rast;
07572 s1r = -zet1dr + str;
07573 s1i = -zet1di + sti;
07574 goto L230;
07575 L220:
07576 s1r = -zet1dr + zet2dr;
07577 s1i = -zet1di + zet2di;
07578 L230:
07579
07580
07581
07582 rs1 = s1r;
07583 if (abs(rs1) > *elim) {
07584 goto L280;
07585 }
07586 if (kdflg == 1) {
07587 iflag = 2;
07588 }
07589 if (abs(rs1) < *alim) {
07590 goto L240;
07591 }
07592
07593
07594
07595 aphi = myzabs_(&phidr, &phidi);
07596 aarg = myzabs_(&argdr, &argdi);
07597 rs1 = rs1 + log(aphi) - log(aarg) * .25 - aic;
07598 if (abs(rs1) > *elim) {
07599 goto L280;
07600 }
07601 if (kdflg == 1) {
07602 iflag = 1;
07603 }
07604 if (rs1 < 0.) {
07605 goto L240;
07606 }
07607 if (kdflg == 1) {
07608 iflag = 3;
07609 }
07610 L240:
07611 zairy_(&argdr, &argdi, &c__0, &c__2, &air, &aii, &nai, &idum);
07612 zairy_(&argdr, &argdi, &c__1, &c__2, &dair, &daii, &ndai, &idum);
07613 str = dair * bsumdr - daii * bsumdi;
07614 sti = dair * bsumdi + daii * bsumdr;
07615 str += air * asumdr - aii * asumdi;
07616 sti += air * asumdi + aii * asumdr;
07617 ptr = str * phidr - sti * phidi;
07618 pti = str * phidi + sti * phidr;
07619 s2r = ptr * csr - pti * csi;
07620 s2i = ptr * csi + pti * csr;
07621 str = exp(s1r) * cssr[iflag - 1];
07622 s1r = str * cos(s1i);
07623 s1i = str * sin(s1i);
07624 str = s2r * s1r - s2i * s1i;
07625 s2i = s2r * s1i + s2i * s1r;
07626 s2r = str;
07627 if (iflag != 1) {
07628 goto L250;
07629 }
07630 zuchk_(&s2r, &s2i, &nw, bry, tol);
07631 if (nw == 0) {
07632 goto L250;
07633 }
07634 s2r = zeror;
07635 s2i = zeroi;
07636 L250:
07637 if (yy <= 0.) {
07638 s2i = -s2i;
07639 }
07640 cyr[kdflg - 1] = s2r;
07641 cyi[kdflg - 1] = s2i;
07642 c2r = s2r;
07643 c2i = s2i;
07644 s2r *= csrr[iflag - 1];
07645 s2i *= csrr[iflag - 1];
07646
07647
07648
07649 s1r = yr[kk];
07650 s1i = yi[kk];
07651 if (*kode == 1) {
07652 goto L270;
07653 }
07654 zs1s2_(&zrr, &zri, &s1r, &s1i, &s2r, &s2i, &nw, &asc, alim, &iuf);
07655 *nz += nw;
07656 L270:
07657 yr[kk] = s1r * cspnr - s1i * cspni + s2r;
07658 yi[kk] = s1r * cspni + s1i * cspnr + s2i;
07659 --kk;
07660 cspnr = -cspnr;
07661 cspni = -cspni;
07662 str = csi;
07663 csi = -csr;
07664 csr = str;
07665 if (c2r != 0. || c2i != 0.) {
07666 goto L255;
07667 }
07668 kdflg = 1;
07669 goto L290;
07670 L255:
07671 if (kdflg == 2) {
07672 goto L295;
07673 }
07674 kdflg = 2;
07675 goto L290;
07676 L280:
07677 if (rs1 > 0.) {
07678 goto L320;
07679 }
07680 s2r = zeror;
07681 s2i = zeroi;
07682 goto L250;
07683 L290:
07684 ;
07685 }
07686 k = *n;
07687 L295:
07688 il = *n - k;
07689 if (il == 0) {
07690 return 0;
07691 }
07692
07693
07694
07695
07696
07697 s1r = cyr[0];
07698 s1i = cyi[0];
07699 s2r = cyr[1];
07700 s2i = cyi[1];
07701 csr = csrr[iflag - 1];
07702 ascle = bry[iflag - 1];
07703 fn = (doublereal) ((real) (inu + il));
07704 i__1 = il;
07705 for (i__ = 1; i__ <= i__1; ++i__) {
07706 c2r = s2r;
07707 c2i = s2i;
07708 s2r = s1r + (fn + fnf) * (rzr * c2r - rzi * c2i);
07709 s2i = s1i + (fn + fnf) * (rzr * c2i + rzi * c2r);
07710 s1r = c2r;
07711 s1i = c2i;
07712 fn += -1.;
07713 c2r = s2r * csr;
07714 c2i = s2i * csr;
07715 ckr = c2r;
07716 cki = c2i;
07717 c1r = yr[kk];
07718 c1i = yi[kk];
07719 if (*kode == 1) {
07720 goto L300;
07721 }
07722 zs1s2_(&zrr, &zri, &c1r, &c1i, &c2r, &c2i, &nw, &asc, alim, &iuf);
07723 *nz += nw;
07724 L300:
07725 yr[kk] = c1r * cspnr - c1i * cspni + c2r;
07726 yi[kk] = c1r * cspni + c1i * cspnr + c2i;
07727 --kk;
07728 cspnr = -cspnr;
07729 cspni = -cspni;
07730 if (iflag >= 3) {
07731 goto L310;
07732 }
07733 c2r = abs(ckr);
07734 c2i = abs(cki);
07735 c2m = max(c2r,c2i);
07736 if (c2m <= ascle) {
07737 goto L310;
07738 }
07739 ++iflag;
07740 ascle = bry[iflag - 1];
07741 s1r *= csr;
07742 s1i *= csr;
07743 s2r = ckr;
07744 s2i = cki;
07745 s1r *= cssr[iflag - 1];
07746 s1i *= cssr[iflag - 1];
07747 s2r *= cssr[iflag - 1];
07748 s2i *= cssr[iflag - 1];
07749 csr = csrr[iflag - 1];
07750 L310:
07751 ;
07752 }
07753 return 0;
07754 L320:
07755 *nz = -1;
07756 return 0;
07757 }
07758
07759 int zuoik_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *ikflg, integer *n, doublereal *yr, doublereal *yi, integer *nuf, doublereal *tol, doublereal *elim, doublereal *alim)
07760 {
07761
07762
07763 static doublereal zeror = 0.;
07764 static doublereal zeroi = 0.;
07765 static doublereal aic = 1.265512123484645396;
07766
07767
07768 integer i__1;
07769
07770
07771
07772
07773 static doublereal aarg, aphi, argi, phii, argr;
07774 static integer idum;
07775 static doublereal phir;
07776 static integer init;
07777 static doublereal sumi, sumr;
07778 static integer i__;
07779 static doublereal ascle;
07780 static integer iform;
07781 static doublereal asumi, bsumi, cwrki[16];
07782 static doublereal asumr, bsumr, cwrkr[16];
07783 static doublereal zeta1i, zeta2i, zeta1r, zeta2r, ax, ay;
07784 static integer nn, nw;
07785 static doublereal fnn, gnn, zbi, czi, gnu, zbr, czr, rcz, sti, zni, zri,
07786 str, znr, zrr;
07787
07788
07789
07790
07791
07792
07793
07794
07795
07796
07797
07798
07799
07800
07801
07802
07803
07804
07805
07806
07807
07808
07809
07810
07811
07812
07813
07814
07815
07816 --yi;
07817 --yr;
07818
07819
07820 *nuf = 0;
07821 nn = *n;
07822 zrr = *zr;
07823 zri = *zi;
07824 if (*zr >= 0.) {
07825 goto L10;
07826 }
07827 zrr = -(*zr);
07828 zri = -(*zi);
07829 L10:
07830 zbr = zrr;
07831 zbi = zri;
07832 ax = abs(*zr) * 1.7321;
07833 ay = abs(*zi);
07834 iform = 1;
07835 if (ay > ax) {
07836 iform = 2;
07837 }
07838 gnu = max(*fnu,1.);
07839 if (*ikflg == 1) {
07840 goto L20;
07841 }
07842 fnn = (doublereal) ((real) nn);
07843 gnn = *fnu + fnn - 1.;
07844 gnu = max(gnn,fnn);
07845 L20:
07846
07847
07848
07849
07850
07851 if (iform == 2) {
07852 goto L30;
07853 }
07854 init = 0;
07855 zunik_(&zrr, &zri, &gnu, ikflg, &c__1, tol, &init, &phir, &phii, &zeta1r,
07856 &zeta1i, &zeta2r, &zeta2i, &sumr, &sumi, cwrkr, cwrki);
07857 czr = -zeta1r + zeta2r;
07858 czi = -zeta1i + zeta2i;
07859 goto L50;
07860 L30:
07861 znr = zri;
07862 zni = -zrr;
07863 if (*zi > 0.) {
07864 goto L40;
07865 }
07866 znr = -znr;
07867 L40:
07868 zunhj_(&znr, &zni, &gnu, &c__1, tol, &phir, &phii, &argr, &argi, &zeta1r,
07869 &zeta1i, &zeta2r, &zeta2i, &asumr, &asumi, &bsumr, &bsumi);
07870 czr = -zeta1r + zeta2r;
07871 czi = -zeta1i + zeta2i;
07872 aarg = myzabs_(&argr, &argi);
07873 L50:
07874 if (*kode == 1) {
07875 goto L60;
07876 }
07877 czr -= zbr;
07878 czi -= zbi;
07879 L60:
07880 if (*ikflg == 1) {
07881 goto L70;
07882 }
07883 czr = -czr;
07884 czi = -czi;
07885 L70:
07886 aphi = myzabs_(&phir, &phii);
07887 rcz = czr;
07888
07889
07890
07891 if (rcz > *elim) {
07892 goto L210;
07893 }
07894 if (rcz < *alim) {
07895 goto L80;
07896 }
07897 rcz += log(aphi);
07898 if (iform == 2) {
07899 rcz = rcz - log(aarg) * .25 - aic;
07900 }
07901 if (rcz > *elim) {
07902 goto L210;
07903 }
07904 goto L130;
07905 L80:
07906
07907
07908
07909 if (rcz < -(*elim)) {
07910 goto L90;
07911 }
07912 if (rcz > -(*alim)) {
07913 goto L130;
07914 }
07915 rcz += log(aphi);
07916 if (iform == 2) {
07917 rcz = rcz - log(aarg) * .25 - aic;
07918 }
07919 if (rcz > -(*elim)) {
07920 goto L110;
07921 }
07922 L90:
07923 i__1 = nn;
07924 for (i__ = 1; i__ <= i__1; ++i__) {
07925 yr[i__] = zeror;
07926 yi[i__] = zeroi;
07927
07928 }
07929 *nuf = nn;
07930 return 0;
07931 L110:
07932 ascle = d1mach_(&c__1) * 1e3 / *tol;
07933 zlog_(&phir, &phii, &str, &sti, &idum);
07934 czr += str;
07935 czi += sti;
07936 if (iform == 1) {
07937 goto L120;
07938 }
07939 zlog_(&argr, &argi, &str, &sti, &idum);
07940 czr = czr - str * .25 - aic;
07941 czi -= sti * .25;
07942 L120:
07943 ax = exp(rcz) / *tol;
07944 ay = czi;
07945 czr = ax * cos(ay);
07946 czi = ax * sin(ay);
07947 zuchk_(&czr, &czi, &nw, &ascle, tol);
07948 if (nw != 0) {
07949 goto L90;
07950 }
07951 L130:
07952 if (*ikflg == 2) {
07953 return 0;
07954 }
07955 if (*n == 1) {
07956 return 0;
07957 }
07958
07959
07960
07961 L140:
07962 gnu = *fnu + (doublereal) ((real) (nn - 1));
07963 if (iform == 2) {
07964 goto L150;
07965 }
07966 init = 0;
07967 zunik_(&zrr, &zri, &gnu, ikflg, &c__1, tol, &init, &phir, &phii, &zeta1r,
07968 &zeta1i, &zeta2r, &zeta2i, &sumr, &sumi, cwrkr, cwrki);
07969 czr = -zeta1r + zeta2r;
07970 czi = -zeta1i + zeta2i;
07971 goto L160;
07972 L150:
07973 zunhj_(&znr, &zni, &gnu, &c__1, tol, &phir, &phii, &argr, &argi, &zeta1r,
07974 &zeta1i, &zeta2r, &zeta2i, &asumr, &asumi, &bsumr, &bsumi);
07975 czr = -zeta1r + zeta2r;
07976 czi = -zeta1i + zeta2i;
07977 aarg = myzabs_(&argr, &argi);
07978 L160:
07979 if (*kode == 1) {
07980 goto L170;
07981 }
07982 czr -= zbr;
07983 czi -= zbi;
07984 L170:
07985 aphi = myzabs_(&phir, &phii);
07986 rcz = czr;
07987 if (rcz < -(*elim)) {
07988 goto L180;
07989 }
07990 if (rcz > -(*alim)) {
07991 return 0;
07992 }
07993 rcz += log(aphi);
07994 if (iform == 2) {
07995 rcz = rcz - log(aarg) * .25 - aic;
07996 }
07997 if (rcz > -(*elim)) {
07998 goto L190;
07999 }
08000 L180:
08001 yr[nn] = zeror;
08002 yi[nn] = zeroi;
08003 --nn;
08004 ++(*nuf);
08005 if (nn == 0) {
08006 return 0;
08007 }
08008 goto L140;
08009 L190:
08010 ascle = d1mach_(&c__1) * 1e3 / *tol;
08011 zlog_(&phir, &phii, &str, &sti, &idum);
08012 czr += str;
08013 czi += sti;
08014 if (iform == 1) {
08015 goto L200;
08016 }
08017 zlog_(&argr, &argi, &str, &sti, &idum);
08018 czr = czr - str * .25 - aic;
08019 czi -= sti * .25;
08020 L200:
08021 ax = exp(rcz) / *tol;
08022 ay = czi;
08023 czr = ax * cos(ay);
08024 czi = ax * sin(ay);
08025 zuchk_(&czr, &czi, &nw, &ascle, tol);
08026 if (nw != 0) {
08027 goto L180;
08028 }
08029 return 0;
08030 L210:
08031 *nuf = -1;
08032 return 0;
08033 }
08034
08035 int zwrsk_(doublereal *zrr, doublereal *zri, doublereal *fnu, integer *kode, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *cwr, doublereal *cwi, doublereal *tol, doublereal *elim, doublereal *alim)
08036 {
08037
08038 integer i__1;
08039
08040
08041
08042
08043 static doublereal ract;
08044 static integer i__;
08045 static doublereal ascle, csclr, cinui, cinur;
08046 static integer nw;
08047 static doublereal c1i, c2i;
08048 static doublereal c1r, c2r, act, acw, cti, ctr, pti, sti, ptr, str;
08049
08050
08051
08052
08053
08054
08055
08056
08057
08058
08059
08060
08061
08062
08063
08064
08065 --yi;
08066 --yr;
08067 --cwr;
08068 --cwi;
08069
08070
08071 *nz = 0;
08072 zbknu_(zrr, zri, fnu, kode, &c__2, &cwr[1], &cwi[1], &nw, tol, elim, alim)
08073 ;
08074 if (nw != 0) {
08075 goto L50;
08076 }
08077 zrati_(zrr, zri, fnu, n, &yr[1], &yi[1], tol);
08078
08079
08080
08081
08082 cinur = 1.;
08083 cinui = 0.;
08084 if (*kode == 1) {
08085 goto L10;
08086 }
08087 cinur = cos(*zri);
08088 cinui = sin(*zri);
08089 L10:
08090
08091
08092
08093
08094
08095
08096 acw = myzabs_(&cwr[2], &cwi[2]);
08097 ascle = d1mach_(&c__1) * 1e3 / *tol;
08098 csclr = 1.;
08099 if (acw > ascle) {
08100 goto L20;
08101 }
08102 csclr = 1. / *tol;
08103 goto L30;
08104 L20:
08105 ascle = 1. / ascle;
08106 if (acw < ascle) {
08107 goto L30;
08108 }
08109 csclr = *tol;
08110 L30:
08111 c1r = cwr[1] * csclr;
08112 c1i = cwi[1] * csclr;
08113 c2r = cwr[2] * csclr;
08114 c2i = cwi[2] * csclr;
08115 str = yr[1];
08116 sti = yi[1];
08117
08118
08119
08120
08121 ptr = str * c1r - sti * c1i;
08122 pti = str * c1i + sti * c1r;
08123 ptr += c2r;
08124 pti += c2i;
08125 ctr = *zrr * ptr - *zri * pti;
08126 cti = *zrr * pti + *zri * ptr;
08127 act = myzabs_(&ctr, &cti);
08128 ract = 1. / act;
08129 ctr *= ract;
08130 cti = -cti * ract;
08131 ptr = cinur * ract;
08132 pti = cinui * ract;
08133 cinur = ptr * ctr - pti * cti;
08134 cinui = ptr * cti + pti * ctr;
08135 yr[1] = cinur * csclr;
08136 yi[1] = cinui * csclr;
08137 if (*n == 1) {
08138 return 0;
08139 }
08140 i__1 = *n;
08141 for (i__ = 2; i__ <= i__1; ++i__) {
08142 ptr = str * cinur - sti * cinui;
08143 cinui = str * cinui + sti * cinur;
08144 cinur = ptr;
08145 str = yr[i__];
08146 sti = yi[i__];
08147 yr[i__] = cinur * csclr;
08148 yi[i__] = cinui * csclr;
08149
08150 }
08151 return 0;
08152 L50:
08153 *nz = -1;
08154 if (nw == -2) {
08155 *nz = -2;
08156 }
08157 return 0;
08158 }
08159