← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/python-scipy/oneiric-proposed

« back to all changes in this revision

Viewing changes to scipy/special/specfun/specfun.f

  • Committer: Bazaar Package Importer
  • Author(s): Varun Hiremath
  • Date: 2011-04-06 21:26:25 UTC
  • mfrom: (9.2.1 sid)
  • Revision ID: james.westby@ubuntu.com-20110406212625-3izdplobqe6fzeql
Tags: 0.9.0+dfsg1-1
http://bugs.debian.org/614407
http://bugs.debian.org/579041
http://bugs.debian.org/569008
http://bugs.debian.org/617028
* New upstream release (Closes: #614407, #579041, #569008)
* Convert to dh_python2 (Closes: #617028)

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/special/specfun/specfun.f
1
1
C       COMPUTATION OF SPECIAL FUNCTIONS
2
 
C 
 
2
C
3
3
C          Shanjie Zhang and Jianming Jin
4
4
C
5
 
C       Copyrighted but permission granted to use code in programs. 
 
5
C       Copyrighted but permission granted to use code in programs.
6
6
C       Buy their book "Computation of Special Functions", 1996, John Wiley & Sons, Inc.
7
7
C
8
8
C
9
9
C      Compiled into a single source file and changed REAL To DBLE throughout.
10
10
C
11
11
C      Changed according to ERRATA also.
12
 
C      
 
12
C
13
13
C      Changed GAMMA to GAMMA2 and PSI to PSI_SPEC to avoid potential conflicts.
14
14
C
15
15
 
65
65
        END
66
66
 
67
67
 
68
 
        
 
68
 
69
69
C       **********************************
70
70
 
71
71
        SUBROUTINE CFS(Z,ZF,ZD)
135
135
C       ==========================================================
136
136
C       Purpose: Compute the associated Legendre functions of the
137
137
C                second kind, Qmn(x) and Qmn'(x)
138
 
C       Input :  x  --- Argument of Qmn(x) 
 
138
C       Input :  x  --- Argument of Qmn(x)
139
139
C                m  --- Order of Qmn(x)  ( m = 0,1,2,… )
140
140
C                n  --- Degree of Qmn(x) ( n = 0,1,2,… )
141
141
C                mm --- Physical dimension of QM and QD
224
224
        SUBROUTINE CLPMN(MM,M,N,X,Y,CPM,CPD)
225
225
C
226
226
C       =========================================================
227
 
C       Purpose: Compute the associated Legendre functions Pmn(z)   
228
 
C                and their derivatives Pmn'(z) for a complex 
 
227
C       Purpose: Compute the associated Legendre functions Pmn(z)
 
228
C                and their derivatives Pmn'(z) for a complex
229
229
C                argument
230
230
C       Input :  x  --- Real part of z
231
231
C                y  --- Imaginary part of z
336
336
        END
337
337
 
338
338
 
339
 
        
 
339
 
340
340
C       **********************************
341
341
C       SciPy: Changed P from a character array to an integer array.
342
342
        SUBROUTINE JDZO(NT,N,M,P,ZO)
357
357
C                P(L)  --- 0 (TM) or 1 (TE), a code for designating the
358
358
C                          zeros of Jn(x)  or Jn'(x).
359
359
C                          In the waveguide applications, the zeros
360
 
C                          of Jn(x) correspond to TM modes and 
 
360
C                          of Jn(x) correspond to TM modes and
361
361
C                          those of Jn'(x) correspond to TE modes
362
 
C       Routine called:    BJNDD for computing Jn(x), Jn'(x) and  
 
362
C       Routine called:    BJNDD for computing Jn(x), Jn'(x) and
363
363
C                          Jn''(x)
364
364
C       =============================================================
365
365
C
367
367
        INTEGER P(1400), P1(70)
368
368
        DIMENSION N(1400),M(1400),ZO(0:1400),N1(70),M1(70),
369
369
     &            ZOC(0:70),BJ(101),DJ(101),FJ(101)
 
370
        X = 0
 
371
        ZOC(0) = 0
370
372
        IF (NT.LT.600) THEN
371
373
           XM=-1.0+2.248485*NT**0.5-.0159382*NT+3.208775E-4
372
374
     &        *NT**1.5
449
451
        END
450
452
 
451
453
 
452
 
        
 
454
 
453
455
C       **********************************
454
456
 
455
457
        SUBROUTINE CBK(M,N,C,CV,QT,CK,BK)
519
521
        END
520
522
 
521
523
 
522
 
        
 
524
 
523
525
C       **********************************
524
526
 
525
527
        SUBROUTINE CJY01(Z,CBJ0,CDJ0,CBJ1,CDJ1,CBY0,CDY0,CBY1,CDY1)
526
528
C
527
529
C       =======================================================
528
 
C       Purpose: Compute Bessel functions J0(z), J1(z), Y0(z), 
 
530
C       Purpose: Compute Bessel functions J0(z), J1(z), Y0(z),
529
531
C                Y1(z), and their derivatives for a complex
530
532
C                argument
531
533
C       Input :  z --- Complex argument
670
672
C                of the second kind with a small argument
671
673
C       Routines called:
672
674
C            (1) LPMNS for computing the associated Legendre
673
 
C                functions of the first kind    
 
675
C                functions of the first kind
674
676
C            (2) LQMNS for computing the associated Legendre
675
 
C                functions of the second kind  
 
677
C                functions of the second kind
676
678
C            (3) KMN for computing expansion coefficients
677
679
C                and joining factors
678
680
C       ======================================================
767
769
        END
768
770
 
769
771
 
770
 
        
 
772
 
771
773
C       **********************************
772
774
 
773
775
        SUBROUTINE BERNOB(N,BN)
780
782
C
781
783
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
782
784
        DIMENSION BN(0:N)
783
 
        TPI=6.283185307179586D0 
 
785
        TPI=6.283185307179586D0
784
786
        BN(0)=1.0D0
785
787
        BN(1)=-0.5D0
786
788
        BN(2)=1.0D0/6.0D0
847
849
10               SK=SK+CK(K+1)*CK(L-K+1)
848
850
15            S=S+SK*AP(I-L+1)
849
851
20      AP(I+1)=-R*S
850
 
        QS0=AP(M+1)     
 
852
        QS0=AP(M+1)
851
853
        DO 30 L=1,M
852
854
           R=1.0D0
853
855
           DO 25 K=1,L
859
861
        END
860
862
 
861
863
 
862
 
        
 
864
 
863
865
C       **********************************
864
866
 
865
867
        SUBROUTINE CV0(KD,M,Q,A0)
1026
1028
        END
1027
1029
 
1028
1030
 
1029
 
        
 
1031
 
1030
1032
C       **********************************
1031
1033
 
1032
1034
        SUBROUTINE CVQM(M,Q,A0)
1083
1085
        END
1084
1086
 
1085
1087
 
1086
 
        
 
1088
 
1087
1089
C       **********************************
1088
1090
 
1089
1091
        SUBROUTINE CSPHJY(N,Z,NM,CSJ,CDJ,CSY,CDY)
1169
1171
        INTEGER FUNCTION MSTA1(X,MP)
1170
1172
C
1171
1173
C       ===================================================
1172
 
C       Purpose: Determine the starting point for backward  
1173
 
C                recurrence such that the magnitude of    
 
1174
C       Purpose: Determine the starting point for backward
 
1175
C                recurrence such that the magnitude of
1174
1176
C                Jn(x) at that point is about 10^(-MP)
1175
1177
C       Input :  x     --- Argument of Jn(x)
1176
1178
C                MP    --- Value of magnitude
1177
 
C       Output:  MSTA1 --- Starting point   
 
1179
C       Output:  MSTA1 --- Starting point
1178
1180
C       ===================================================
1179
1181
C
1180
1182
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1183
1185
        F0=ENVJ(N0,A0)-MP
1184
1186
        N1=N0+5
1185
1187
        F1=ENVJ(N1,A0)-MP
1186
 
        DO 10 IT=1,20             
1187
 
           NN=N1-(N1-N0)/(1.0D0-F0/F1)                  
 
1188
        DO 10 IT=1,20
 
1189
           NN=N1-(N1-N0)/(1.0D0-F0/F1)
1188
1190
           F=ENVJ(NN,A0)-MP
1189
1191
           IF(ABS(NN-N1).LT.1) GO TO 20
1190
1192
           N0=N1
1442
1444
        END
1443
1445
 
1444
1446
 
1445
 
        
 
1447
 
1446
1448
C       **********************************
1447
1449
 
1448
1450
        SUBROUTINE RMN2L(M,N,C,X,DF,KD,R2F,R2D,ID)
1453
1455
C                c and a large cx
1454
1456
C       Routine called:
1455
1457
C                SPHY for computing the spherical Bessel
1456
 
C                functions of the second kind      
 
1458
C                functions of the second kind
1457
1459
C       ========================================================
1458
1460
C
1459
1461
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1471
1473
        R0=REG
1472
1474
        DO 10 J=1,2*M+IP
1473
1475
10         R0=R0*J
1474
 
        R=R0    
 
1476
        R=R0
1475
1477
        SUC=R*DF(1)
1476
1478
        SW=0.0D0
1477
1479
        DO 15 K=2,NM
1503
1505
           ID=10
1504
1506
           RETURN
1505
1507
        ENDIF
1506
 
        B0=KD*M/X**3.0D0/(1.0-KD/(X*X))*R2F                
 
1508
        B0=KD*M/X**3.0D0/(1.0-KD/(X*X))*R2F
1507
1509
        SUD=0.0D0
1508
1510
        EPS2=0.0D0
1509
1511
        DO 60 K=1,NM
1520
1522
           EPS2=DABS(SUD-SW)
1521
1523
           IF (K.GT.NM1.AND.EPS2.LT.DABS(SUD)*EPS) GO TO 65
1522
1524
60         SW=SUD
1523
 
65      R2D=B0+A0*C*SUD       
 
1525
65      R2D=B0+A0*C*SUD
1524
1526
        ID2=INT(LOG10(EPS2/DABS(SUD)+EPS))
1525
1527
        ID=MAX(ID1,ID2)
1526
1528
        RETURN
1527
1529
        END
1528
1530
 
1529
1531
 
1530
 
        
 
1532
 
1531
1533
C       **********************************
1532
1534
 
1533
1535
        SUBROUTINE PSI_SPEC(X,PS)
1667
1669
        END
1668
1670
 
1669
1671
 
1670
 
        
 
1672
 
1671
1673
C       **********************************
1672
1674
 
1673
1675
        SUBROUTINE LPMNS(M,N,X,PM,PD)
1719
1721
           PM(K)=PM2
1720
1722
           PMK=PM1
1721
1723
25         PM1=PM2
1722
 
        PD(0)=((1.0D0-M)*PM(1)-X*PM(0))/(X*X-1.0)  
 
1724
        PD(0)=((1.0D0-M)*PM(1)-X*PM(0))/(X*X-1.0)
1723
1725
        DO 30 K=1,N
1724
1726
30         PD(K)=(K*X*PM(K)-(K+M)*PM(K-1))/(X*X-1.0D0)
1725
1727
        DO 35 K=1,N
1822
1824
C                        the second kind
1823
1825
C       Routines called:
1824
1826
C            (1) SDMN for computing expansion coefficients dk
1825
 
C            (2) RMN1 for computing prolate and oblate radial 
 
1827
C            (2) RMN1 for computing prolate and oblate radial
1826
1828
C                functions of the first kind
1827
1829
C            (3) RMN2L for computing prolate and oblate radial
1828
1830
C                functions of the second kind for a large argument
1829
 
C            (4) RMN2SP for computing the prolate radial function 
 
1831
C            (4) RMN2SP for computing the prolate radial function
1830
1832
C                of the second kind for a small argument
1831
1833
C       ==============================================================
1832
1834
C
1847
1849
        END
1848
1850
 
1849
1851
 
1850
 
        
 
1852
 
1851
1853
C       **********************************
1852
1854
 
1853
1855
        SUBROUTINE JYNDD(N,X,BJN,DJN,FJN,BYN,DYN,FYN)
1854
1856
C
1855
1857
C       ===========================================================
1856
1858
C       Purpose: Compute Bessel functions Jn(x) and Yn(x), and
1857
 
C                their first and second derivatives 
 
1859
C                their first and second derivatives
1858
1860
C       Input:   x   ---  Argument of Jn(x) and Yn(x) ( x > 0 )
1859
1861
C                n   ---  Order of Jn(x) and Yn(x)
1860
1862
C       Output:  BJN ---  Jn(x)
2012
2014
        END
2013
2015
 
2014
2016
 
2015
 
        
 
2017
 
2016
2018
C       **********************************
2017
2019
 
2018
2020
        SUBROUTINE CISIA(X,CI,SI)
2314
2316
C                          KF=1 for Ai(x) and Ai'(x)
2315
2317
C                          KF=2 for Bi(x) and Bi'(x)
2316
2318
C       Output:  XA(m) --- a, the m-th zero of Ai(x) or
2317
 
C                          b, the m-th zero of Bi(x) 
 
2319
C                          b, the m-th zero of Bi(x)
2318
2320
C                XB(m) --- a', the m-th zero of Ai'(x) or
2319
2321
C                          b', the m-th zero of Bi'(x)
2320
2322
C                XC(m) --- Ai(a') or Bi(b')
2395
2397
        END
2396
2398
 
2397
2399
 
2398
 
        
 
2400
 
2399
2401
C       **********************************
2400
2402
 
2401
2403
        SUBROUTINE ERROR(X,ERR)
2473
2475
              IF (CDABS(CR/CS).LT.1.0D-15) GO TO 15
2474
2476
10         CONTINUE
2475
2477
15         CER=2.0D0*C0*CS/DSQRT(PI)
2476
 
        ELSE                              
2477
 
           CL=1.0D0/Z1              
 
2478
        ELSE
 
2479
           CL=1.0D0/Z1
2478
2480
           CR=CL
2479
2481
C
2480
2482
C          Asymptotic series; maximum K must be at most ~ R^2.
2535
2537
C
2536
2538
C       ============================================================
2537
2539
C       Purpose: Compute a sequence of characteristic values of
2538
 
C                Mathieu functions 
 
2540
C                Mathieu functions
2539
2541
C       Input :  M  --- Maximum order of Mathieu functions
2540
2542
C                q  --- Parameter of Mathieu functions
2541
2543
C                KD --- Case code
2778
2780
C
2779
2781
C       ========================================================
2780
2782
C       Purpose: Compute the expansion coefficients for the
2781
 
C                asymptotic expansion of Bessel functions 
 
2783
C                asymptotic expansion of Bessel functions
2782
2784
C                with large orders
2783
2785
C       Input :  Km   --- Maximum k
2784
 
C       Output:  A(L) --- Cj(k) where j and k are related to L 
 
2786
C       Output:  A(L) --- Cj(k) where j and k are related to L
2785
2787
C                         by L=j+1+[k*(k+1)]/2; j,k=0,1,...,Km
2786
2788
C       ========================================================
2787
2789
C
2887
2889
C       Purpose: Compute lambda function with arbitrary order v,
2888
2890
C                and their derivative
2889
2891
C       Input :  x --- Argument of lambda function
2890
 
C                v --- Order of lambda function 
 
2892
C                v --- Order of lambda function
2891
2893
C       Output:  VL(n) --- Lambda function of order n+v0
2892
 
C                DL(n) --- Derivative of lambda function 
 
2894
C                DL(n) --- Derivative of lambda function
2893
2895
C                VM --- Highest order computed
2894
2896
C       Routines called:
2895
 
C            (1) MSTA1 and MSTA2 for computing the starting 
 
2897
C            (1) MSTA1 and MSTA2 for computing the starting
2896
2898
C                point for backward recurrence
2897
2899
C            (2) GAM0 for computing gamma function (|x| ≤ 1)
2898
2900
C       =========================================================
3014
3016
        END
3015
3017
 
3016
3018
 
3017
 
        
 
3019
 
3018
3020
C       **********************************
3019
3021
 
3020
3022
        SUBROUTINE CHGUIT(A,B,X,HU,ID)
3116
3118
        END
3117
3119
 
3118
3120
 
3119
 
        
 
3121
 
3120
3122
C       **********************************
3121
3123
 
3122
3124
        SUBROUTINE KMN(M,N,C,CV,KD,DF,DN,CK1,CK2)
3136
3138
        IP=1
3137
3139
        IF (N-M.EQ.2*INT((N-M)/2)) IP=0
3138
3140
        K=0
3139
 
        DO 10 I=1,NN+3     
 
3141
        DO 10 I=1,NN+3
3140
3142
           IF (IP.EQ.0) K=-2*(I-1)
3141
3143
           IF (IP.EQ.1) K=-(2*I-3)
3142
3144
           GK0=2.0D0*M+K
3204
3206
        END
3205
3207
 
3206
3208
 
3207
 
        
 
3209
 
3208
3210
C       **********************************
3209
3211
 
3210
3212
        SUBROUTINE LAGZO(N,X,W)
3298
3300
        END
3299
3301
 
3300
3302
 
3301
 
        
 
3303
 
3302
3304
C       **********************************
3303
3305
 
3304
3306
        SUBROUTINE CJYVA(V,Z,VM,CBJ,CDJ,CBY,CDY)
3305
3307
C
3306
3308
C       ===========================================================
3307
 
C       Purpose: Compute Bessel functions Jv(z), Yv(z) and their 
 
3309
C       Purpose: Compute Bessel functions Jv(z), Yv(z) and their
3308
3310
C                derivatives for a complex argument
3309
3311
C       Input :  z --- Complex argument
3310
3312
C                v --- Order of Jv(z) and Yv(z)
3345
3347
           ELSE
3346
3348
              CDJ(0)=(1.0D+300,0.0D0)
3347
3349
           ENDIF
3348
 
           VM=V                     
 
3350
           VM=V
3349
3351
           RETURN
3350
3352
        ENDIF
3351
3353
        LB0=0.0D0
3554
3556
        END
3555
3557
 
3556
3558
 
3557
 
        
 
3559
 
3558
3560
C       **********************************
3559
3561
 
3560
3562
        SUBROUTINE CJYVB(V,Z,VM,CBJ,CDJ,CBY,CDY)
3561
3563
C
3562
3564
C       ===========================================================
3563
 
C       Purpose: Compute Bessel functions Jv(z), Yv(z) and their 
 
3565
C       Purpose: Compute Bessel functions Jv(z), Yv(z) and their
3564
3566
C                derivatives for a complex argument
3565
3567
C       Input :  z --- Complex argument
3566
3568
C                v --- Order of Jv(z) and Yv(z)
3716
3718
        END
3717
3719
 
3718
3720
 
3719
 
        
 
3721
 
3720
3722
C       **********************************
3721
3723
 
3722
3724
        SUBROUTINE JY01A(X,BJ0,DJ0,BJ1,DJ1,BY0,DY0,BY1,DY1)
3852
3854
C       Purpose: Compute the incomplete gamma function
3853
3855
C                r(a,x), Г(a,x) and P(a,x)
3854
3856
C       Input :  a   --- Parameter ( a ≤ 170 )
3855
 
C                x   --- Argument 
 
3857
C                x   --- Argument
3856
3858
C       Output:  GIN --- r(a,x)
3857
3859
C                GIM --- Г(a,x)
3858
3860
C                GIP --- P(a,x)
3895
3897
        END
3896
3898
 
3897
3899
 
3898
 
        
 
3900
 
3899
3901
C       **********************************
3900
3902
 
3901
3903
        SUBROUTINE ITIKB(X,TI,TK)
4075
4077
           ELSE
4076
4078
              DJ(0)=1.0D+300
4077
4079
           ENDIF
4078
 
           VM=V  
 
4080
           VM=V
4079
4081
           RETURN
4080
4082
        ENDIF
4081
4083
        BJV0=0.0D0
4225
4227
        END
4226
4228
 
4227
4229
 
4228
 
        
 
4230
 
4229
4231
C       **********************************
4230
4232
 
4231
4233
        SUBROUTINE JYNB(N,X,NM,BJ,DJ,BY,DY)
4253
4255
              DJ(K) = 0.0D0
4254
4256
 10           DY(K) = 1.0D+300
4255
4257
           DJ(1)=0.5D0
4256
 
        ELSE   
 
4258
        ELSE
4257
4259
           DJ(0)=-BJ(1)
4258
4260
           DO 40 K=1,NM
4259
4261
 40           DJ(K)=BJ(K-1)-K/X*BJ(K)
4278
4280
C                BY(n-NMIN) --- Yn(x)   ; if indexing starts at 0
4279
4281
C                NM --- Highest order computed
4280
4282
C       Routines called:
4281
 
C                MSTA1 and MSTA2 to calculate the starting 
 
4283
C                MSTA1 and MSTA2 to calculate the starting
4282
4284
C                point for backward recurrence
4283
4285
C       =====================================================
4284
4286
C
4544
4546
        END
4545
4547
 
4546
4548
 
4547
 
        
 
4549
 
4548
4550
C       **********************************
4549
4551
 
4550
4552
        SUBROUTINE JYNA(N,X,NM,BJ,DJ,BY,DY)
4561
4563
C                NM --- Highest order computed
4562
4564
C       Routines called:
4563
4565
C            (1) JY01B to calculate J0(x), J1(x), Y0(x) & Y1(x)
4564
 
C            (2) MSTA1 and MSTA2 to calculate the starting 
 
4566
C            (2) MSTA1 and MSTA2 to calculate the starting
4565
4567
C                point for backward recurrence
4566
4568
C       =========================================================
4567
4569
C
4632
4634
        END
4633
4635
 
4634
4636
 
4635
 
        
 
4637
 
4636
4638
C       **********************************
4637
4639
 
4638
4640
        SUBROUTINE PBDV(V,X,DV,DP,PDF,PDD)
4644
4646
C                v --- Order of Dv(x)
4645
4647
C       Output:  DV(na) --- Dn+v0(x)
4646
4648
C                DP(na) --- Dn+v0'(x)
4647
 
C                ( na = |n|, v0 = v-n, |v0| < 1, 
 
4649
C                ( na = |n|, v0 = v-n, |v0| < 1,
4648
4650
C                  n = 0,±1,±2,… )
4649
4651
C                PDF --- Dv(x)
4650
4652
C                PDD --- Dv'(x)
4748
4750
        END
4749
4751
 
4750
4752
 
4751
 
        
 
4753
 
4752
4754
C       **********************************
4753
4755
 
4754
4756
        SUBROUTINE ITSH0(X,TH0)
4763
4765
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
4764
4766
        DIMENSION A(25)
4765
4767
        PI=3.141592653589793D0
4766
 
        R=1.0D0            
 
4768
        R=1.0D0
4767
4769
        IF (X.LE.30.0) THEN
4768
4770
           S=0.5D0
4769
4771
           DO 10 K=1,100
4857
4859
        END
4858
4860
 
4859
4861
 
4860
 
        
 
4862
 
4861
4863
C       **********************************
4862
4864
 
4863
4865
        SUBROUTINE GAMMA2(X,GA)
4992
4994
        END
4993
4995
 
4994
4996
 
4995
 
        
 
4997
 
4996
4998
C       **********************************
4997
4999
 
4998
5000
        SUBROUTINE LAMN(N,X,NM,BL,DL)
5042
5044
35         DL(N)=-0.5D0*X/(N+1.0D0)*UK
5043
5045
           RETURN
5044
5046
        ENDIF
5045
 
        IF (N.EQ.0) NM=1          
 
5047
        IF (N.EQ.0) NM=1
5046
5048
        M=MSTA1(X,200)
5047
5049
        IF (M.LT.NM) THEN
5048
5050
           NM=M
5150
5152
        END
5151
5153
 
5152
5154
 
5153
 
        
 
5155
 
5154
5156
C       **********************************
5155
5157
 
5156
5158
        SUBROUTINE CVF(KD,M,Q,A,MJ,F)
5199
5201
        END
5200
5202
 
5201
5203
 
5202
 
        
 
5204
 
5203
5205
C       **********************************
5204
5206
 
5205
5207
        SUBROUTINE CLPN(N,X,Y,CPN,CPD)
5412
5414
        END
5413
5415
 
5414
5416
 
5415
 
        
 
5417
 
5416
5418
C       **********************************
5417
5419
 
5418
5420
        SUBROUTINE ELIT(HK,PHI,FE,EE)
5673
5675
              DO 15 J=1,500
5674
5676
                 RG=RG*(A+J-1.0D0)/(J*(B+J-1.0D0))*X
5675
5677
                 HG=HG+RG
5676
 
                 IF (DABS(RG/HG).LT.1.0D-15) GO TO 25
 
5678
                 IF (HG.NE.0D0.AND.DABS(RG/HG).LT.1.0D-15) GO TO 25
5677
5679
15            CONTINUE
5678
5680
           ELSE
5679
5681
              CALL GAMMA2(A,TA)
5710
5712
        END
5711
5713
 
5712
5714
 
5713
 
        
 
5715
 
5714
5716
C       **********************************
5715
5717
 
5716
5718
        SUBROUTINE STVH0(X,SH0)
5952
5954
        END
5953
5955
 
5954
5956
 
5955
 
        
 
5957
 
5956
5958
C       **********************************
5957
5959
 
5958
5960
        SUBROUTINE CCHG(A,B,Z,CHG)
6069
6071
        END
6070
6072
 
6071
6073
 
6072
 
        
 
6074
 
6073
6075
C       **********************************
6074
6076
 
6075
6077
        SUBROUTINE HYGFZ(A,B,C,Z,ZHF)
6076
6078
C
6077
6079
C       ======================================================
6078
 
C       Purpose: Compute the hypergeometric function for a 
 
6080
C       Purpose: Compute the hypergeometric function for a
6079
6081
C                complex argument, F(a,b,c,z)
6080
6082
C       Input :  a --- Parameter
6081
6083
C                b --- Parameter
6146
6148
        ELSE IF (A0.LE.1.0D0) THEN
6147
6149
           IF (X.LT.0.0D0) THEN
6148
6150
              Z1=Z/(Z-1.0D0)
6149
 
              IF (C.GT.A.AND.B.LT.A.AND.B.GT.0.0) THEN  
 
6151
              IF (C.GT.A.AND.B.LT.A.AND.B.GT.0.0) THEN
6150
6152
                 A=BB
6151
6153
                 B=AA
6152
6154
              ENDIF
6360
6362
        END
6361
6363
 
6362
6364
 
6363
 
        
 
6365
 
6364
6366
C       **********************************
6365
6367
 
6366
6368
        SUBROUTINE ITAIRY(X,APT,BPT,ANT,BNT)
6481
6483
C                NM --- Highest order computed
6482
6484
C       Routines called:
6483
6485
C            (1) IK01A for computing I0(x),I1(x),K0(x) & K1(x)
6484
 
C            (2) MSTA1 and MSTA2 for computing the starting 
 
6486
C            (2) MSTA1 and MSTA2 for computing the starting
6485
6487
C                point for backward recurrence
6486
6488
C       ========================================================
6487
6489
C
6549
6551
        END
6550
6552
 
6551
6553
 
6552
 
        
 
6554
 
6553
6555
C       **********************************
6554
6556
 
6555
6557
        SUBROUTINE CJYNA(N,Z,NM,CBJ,CDJ,CBY,CDY)
6566
6568
C                NM --- Highest order computed
6567
6569
C       Rouitines called:
6568
6570
C            (1) CJY01 to calculate J0(z), J1(z), Y0(z), Y1(z)
6569
 
C            (2) MSTA1 and MSTA2 to calculate the starting 
 
6571
C            (2) MSTA1 and MSTA2 to calculate the starting
6570
6572
C                point for backward recurrence
6571
6573
C       =======================================================
6572
6574
C
6635
6637
        CG1=CBY1
6636
6638
        DO 90 K=2,NM
6637
6639
           CYK=2.0D0*(K-1.0D0)/Z*CG1-CG0
6638
 
           IF (CDABS(CYK).GT.1.0D+290) GO TO 90            
 
6640
           IF (CDABS(CYK).GT.1.0D+290) GO TO 90
6639
6641
           YAK=CDABS(CYK)
6640
6642
           YA1=CDABS(CG0)
6641
6643
           IF (YAK.LT.YA0.AND.YAK.LT.YA1) LB=K
6696
6698
        END
6697
6699
 
6698
6700
 
6699
 
        
 
6701
 
6700
6702
C       **********************************
6701
6703
 
6702
6704
        SUBROUTINE CJYNB(N,Z,NM,CBJ,CDJ,CBY,CDY)
6833
6835
        END
6834
6836
 
6835
6837
 
6836
 
        
 
6838
 
6837
6839
C       **********************************
6838
6840
 
6839
6841
        SUBROUTINE IKNB(N,X,NM,BI,DI,BK,DK)
6849
6851
C                DK(n) --- Kn'(x)
6850
6852
C                NM --- Highest order computed
6851
6853
C       Routines called:
6852
 
C                MSTA1 and MSTA2 for computing the starting point 
 
6854
C                MSTA1 and MSTA2 for computing the starting point
6853
6855
C                for backward recurrence
6854
6856
C       ===========================================================
6855
6857
C
6931
6933
        SUBROUTINE LPMN(MM,M,N,X,PM,PD)
6932
6934
C
6933
6935
C       =====================================================
6934
 
C       Purpose: Compute the associated Legendre functions 
 
6936
C       Purpose: Compute the associated Legendre functions
6935
6937
C                Pmn(x) and their derivatives Pmn'(x)
6936
6938
C       Input :  x  --- Argument of Pmn(x)
6937
6939
C                m  --- Order of Pmn(x),  m = 0,1,2,...,n
7056
7058
        END
7057
7059
 
7058
7060
 
7059
 
        
 
7061
 
7060
7062
C       **********************************
7061
7063
 
7062
7064
        SUBROUTINE CY01(KF,Z,ZF,ZD)
7207
7209
        SUBROUTINE FFK(KS,X,FR,FI,FM,FA,GR,GI,GM,GA)
7208
7210
C
7209
7211
C       =======================================================
7210
 
C       Purpose: Compute modified Fresnel integrals F±(x) 
 
7212
C       Purpose: Compute modified Fresnel integrals F±(x)
7211
7213
C                and K±(x)
7212
7214
C       Input :  x   --- Argument of F±(x) and K±(x)
7213
7215
C                KS  --- Sign code
7378
7380
        END
7379
7381
 
7380
7382
 
7381
 
        
 
7383
 
7382
7384
C       **********************************
7383
7385
 
7384
7386
        SUBROUTINE AIRYB(X,AI,BI,AD,BD)
7599
7601
        END
7600
7602
 
7601
7603
 
7602
 
        
 
7604
 
7603
7605
C       **********************************
7604
7606
 
7605
7607
        SUBROUTINE SCKB(M,N,C,DF,CK)
7649
7651
30            R1=R1*I
7650
7652
35         CK(K+1)=FAC*SUM/R1
7651
7653
        RETURN
7652
 
        END 
7653
 
 
7654
 
 
7655
 
        
 
7654
        END
 
7655
 
 
7656
 
 
7657
 
7656
7658
C       **********************************
7657
7659
 
7658
7660
        SUBROUTINE CPDLA(N,Z,CDN)
7680
7682
        END
7681
7683
 
7682
7684
 
7683
 
        
 
7685
 
7684
7686
C       **********************************
7685
7687
 
7686
7688
        SUBROUTINE FCSZO(KF,NT,ZO)
7692
7694
C                        KF=1 for C(z) or KF=2 for S(z)
7693
7695
C                NT  --- Total number of zeros
7694
7696
C       Output:  ZO(L) --- L-th zero of C(z) or S(z)
7695
 
C       Routines called: 
 
7697
C       Routines called:
7696
7698
C            (1) CFC for computing Fresnel integral C(z)
7697
7699
C            (2) CFS for computing Fresnel integral S(z)
7698
7700
C       ==============================================================
7740
7742
        END
7741
7743
 
7742
7744
 
7743
 
        
 
7745
 
7744
7746
C       **********************************
7745
7747
 
7746
7748
        SUBROUTINE E1XA(X,E1)
7747
7749
C
7748
7750
C       ============================================
7749
7751
C       Purpose: Compute exponential integral E1(x)
7750
 
C       Input :  x  --- Argument of E1(x) 
 
7752
C       Input :  x  --- Argument of E1(x)
7751
7753
C       Output:  E1 --- E1(x) ( x > 0 )
7752
7754
C       ============================================
7753
7755
C
7769
7771
 
7770
7772
C       **********************************
7771
7773
 
7772
 
        SUBROUTINE LPMV(V,M,X,PMV)
 
7774
        SUBROUTINE LPMV0(V,M,X,PMV)
7773
7775
C
7774
7776
C       =======================================================
7775
7777
C       Purpose: Compute the associated Legendre function
7776
 
C                Pmv(x) with an integer order and an arbitrary 
 
7778
C                Pmv(x) with an integer order and an arbitrary
7777
7779
C                nonnegative degree v
7778
7780
C       Input :  x   --- Argument of Pm(x)  ( -1 ≤ x ≤ 1 )
7779
7781
C                m   --- Order of Pmv(x)
7864
7866
        RETURN
7865
7867
        END
7866
7868
 
7867
 
 
7868
 
        
 
7869
C       **********************************
 
7870
 
 
7871
        SUBROUTINE LPMV(V,M,X,PMV)
 
7872
C
 
7873
C       =======================================================
 
7874
C       Purpose: Compute the associated Legendre function
 
7875
C                Pmv(x) with an integer order and an arbitrary
 
7876
C                degree v, using down-recursion for large degrees
 
7877
C       Input :  x   --- Argument of Pm(x)  ( -1 ≤ x ≤ 1 )
 
7878
C                m   --- Order of Pmv(x)
 
7879
C                v   --- Degree of Pmv(x)
 
7880
C       Output:  PMV --- Pmv(x)
 
7881
C       Routine called:  LPMV0
 
7882
C       =======================================================
 
7883
C
 
7884
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
7885
        IF (X.EQ.-1.0D0.AND.V.NE.INT(V)) THEN
 
7886
           IF (M.EQ.0) PMV=-1.0D+300
 
7887
           IF (M.NE.0) PMV=1.0D+300
 
7888
           RETURN
 
7889
        ENDIF
 
7890
        VX=V
 
7891
        MX=M
 
7892
        IF (V.LT.0) THEN
 
7893
           VX=-VX-1
 
7894
        ENDIF
 
7895
        NEG_M=0
 
7896
        IF (M.LT.0.AND.(VX+M+1).GT.0D0) THEN
 
7897
C          XXX: does not handle the cases where AMS 8.2.5
 
7898
C               does not help
 
7899
           NEG_M=1
 
7900
           MX=-M
 
7901
        ENDIF
 
7902
        NV=INT(VX)
 
7903
        V0=VX-NV
 
7904
        IF (NV.GT.2.AND.NV.GT.MX) THEN
 
7905
C          Up-recursion on degree, AMS 8.5.3
 
7906
           CALL LPMV0(V0+MX, MX, X, P0)
 
7907
           CALL LPMV0(V0+MX+1, MX, X, P1)
 
7908
           PMV = P1
 
7909
           DO 10 J=MX+2,NV
 
7910
              PMV = ((2*(V0+J)-1)*X*P1 - (V0+J-1+MX)*P0) / (V0+J-MX)
 
7911
              P0 = P1
 
7912
              P1 = PMV
 
7913
10         CONTINUE
 
7914
        ELSE
 
7915
           CALL LPMV0(VX, MX, X, PMV)
 
7916
        ENDIF
 
7917
        IF (NEG_M.NE.0.AND.ABS(PMV).LT.1.0D+300) THEN
 
7918
C          AMS 8.2.5, for integer order
 
7919
           CALL GAMMA2(VX-MX+1, G1)
 
7920
           CALL GAMMA2(VX+MX+1, G2)
 
7921
           PMV = PMV*G1/G2 * (-1)**MX
 
7922
        ENDIF
 
7923
        END
 
7924
 
 
7925
 
7869
7926
C       **********************************
7870
7927
 
7871
7928
        SUBROUTINE CGAMA(X,Y,KF,GR,GI)
7998
8055
        END
7999
8056
 
8000
8057
 
8001
 
        
 
8058
 
8002
8059
C       **********************************
8003
8060
 
8004
8061
        SUBROUTINE CHGUS(A,B,X,HU,ID)
8047
8104
        END
8048
8105
 
8049
8106
 
8050
 
        
 
8107
 
8051
8108
C       **********************************
8052
8109
 
8053
8110
        SUBROUTINE ITTH0(X,TTH)
8185
8242
C       ====================================================
8186
8243
C
8187
8244
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
8188
 
        PI=3.141592653589793D0           
 
8245
        PI=3.141592653589793D0
8189
8246
        EPS=1.0D-12
8190
8247
        EP=DEXP(-.25*X*X)
8191
8248
        A0=DABS(X)**VA*EP
8207
8264
        END
8208
8265
 
8209
8266
 
8210
 
        
 
8267
 
8211
8268
C       **********************************
8212
8269
 
8213
8270
        SUBROUTINE IK01A(X,BI0,DI0,BI1,DI1,BK0,DK0,BK1,DK1)
8322
8379
        SUBROUTINE CPBDN(N,Z,CPB,CPD)
8323
8380
C
8324
8381
C       ==================================================
8325
 
C       Purpose: Compute the parabolic cylinder functions 
 
8382
C       Purpose: Compute the parabolic cylinder functions
8326
8383
C                 Dn(z) and Dn'(z) for a complex argument
8327
8384
C       Input:   z --- Complex argument of Dn(z)
8328
8385
C                n --- Order of Dn(z)  ( n=0,±1,±2,… )
8410
8467
        END
8411
8468
 
8412
8469
 
8413
 
        
 
8470
 
8414
8471
C       **********************************
8415
8472
 
8416
8473
        SUBROUTINE IK01B(X,BI0,DI0,BI1,DI1,BK0,DK0,BK1,DK1)
8508
8565
        END
8509
8566
 
8510
8567
 
8511
 
        
 
8568
 
8512
8569
C       **********************************
8513
8570
 
8514
8571
        SUBROUTINE LPN(N,X,PN,PD)
8730
8787
        ENDIF
8731
8788
85      IF (FC(1).LT.0.0D0) THEN
8732
8789
           DO 90 J=1,KM
8733
 
90            FC(J)=-FC(J)             
 
8790
90            FC(J)=-FC(J)
8734
8791
        ENDIF
8735
8792
        RETURN
8736
8793
        END
8737
8794
 
8738
8795
 
8739
 
        
 
8796
 
8740
8797
C       **********************************
8741
8798
 
8742
8799
        SUBROUTINE SPHI(N,X,NM,SI,DI)
8892
8949
        END
8893
8950
 
8894
8951
 
8895
 
        
 
8952
 
8896
8953
C       **********************************
8897
8954
 
8898
8955
        SUBROUTINE RMN1(M,N,C,X,DF,KD,R1F,R1D)
8904
8961
C       Routines called:
8905
8962
C            (1) SCKB for computing expansion coefficients c2k
8906
8963
C            (2) SPHJ for computing the spherical Bessel
8907
 
C                functions of the first kind     
 
8964
C                functions of the first kind
8908
8965
C       =======================================================
8909
8966
C
8910
8967
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
8919
8976
        R0=REG
8920
8977
        DO 10 J=1,2*M+IP
8921
8978
10         R0=R0*J
8922
 
        R=R0    
 
8979
        R=R0
8923
8980
        SUC=R*DF(1)
8924
8981
        SW=0.0D0
8925
8982
        DO 15 K=2,NM
8958
9015
        CX=C*X
8959
9016
        NM2=2*NM+M
8960
9017
        CALL SPHJ(NM2,CX,NM2,SJ,DJ)
8961
 
        A0=(1.0D0-KD/(X*X))**(0.5D0*M)/SUC  
 
9018
        A0=(1.0D0-KD/(X*X))**(0.5D0*M)/SUC
8962
9019
        R1F=0.0D0
8963
9020
        SW=0.0D0
8964
9021
        LG=0
8976
9033
           IF (K.GT.NM1.AND.DABS(R1F-SW).LT.DABS(R1F)*EPS) GO TO 55
8977
9034
50         SW=R1F
8978
9035
55      R1F=R1F*A0
8979
 
        B0=KD*M/X**3.0D0/(1.0-KD/(X*X))*R1F    
 
9036
        B0=KD*M/X**3.0D0/(1.0-KD/(X*X))*R1F
8980
9037
        SUD=0.0D0
8981
9038
        SW=0.0D0
8982
9039
        DO 60 K=1,NM
8997
9054
        END
8998
9055
 
8999
9056
 
9000
 
        
 
9057
 
9001
9058
C       **********************************
9002
9059
 
9003
9060
        SUBROUTINE DVSA(VA,X,PD)
9049
9106
        END
9050
9107
 
9051
9108
 
9052
 
        
 
9109
 
9053
9110
C       **********************************
9054
9111
 
9055
9112
        SUBROUTINE E1Z(Z,CE1)
9193
9250
        END
9194
9251
 
9195
9252
 
9196
 
        
 
9253
 
9197
9254
C       **********************************
9198
9255
 
9199
9256
        SUBROUTINE GMN(M,N,C,X,BK,GF,GD)
9232
9289
        END
9233
9290
 
9234
9291
 
9235
 
        
 
9292
 
9236
9293
C       **********************************
9237
9294
 
9238
9295
        SUBROUTINE ITJYA(X,TJ,TY)
9388
9445
        END
9389
9446
 
9390
9447
 
9391
 
        
 
9448
 
9392
9449
C       **********************************
9393
9450
 
9394
9451
        SUBROUTINE RCTY(N,X,NM,RY,DY)
9751
9808
        END
9752
9809
 
9753
9810
 
9754
 
        
 
9811
 
9755
9812
C       **********************************
9756
9813
 
9757
9814
        SUBROUTINE CYZO(NT,KF,KC,ZO,ZV)
9758
9815
C
9759
9816
C       ===========================================================
9760
9817
C       Purpose : Compute the complex zeros of Y0(z), Y1(z) and
9761
 
C                 Y1'(z), and their associated values at the zeros 
 
9818
C                 Y1'(z), and their associated values at the zeros
9762
9819
C                 using the modified Newton's iteration method
9763
9820
C       Input:    NT --- Total number of zeros/roots
9764
9821
C                 KF --- Function choice code
9833
9890
        END
9834
9891
 
9835
9892
 
9836
 
        
 
9893
 
9837
9894
C       **********************************
9838
9895
 
9839
9896
        SUBROUTINE KLVNB(X,BER,BEI,GER,GEI,DER,DEI,HER,HEI)
10009
10066
        END
10010
10067
 
10011
10068
 
10012
 
        
 
10069
 
10013
10070
C       **********************************
10014
10071
 
10015
10072
        SUBROUTINE CSPHIK(N,Z,NM,CSI,CDI,CSK,CDK)
10033
10090
        DOUBLE PRECISION A0,PI
10034
10091
        DIMENSION CSI(0:N),CDI(0:N),CSK(0:N),CDK(0:N)
10035
10092
        PI=3.141592653589793D0
10036
 
        A0=CDABS(Z)            
 
10093
        A0=CDABS(Z)
10037
10094
        NM=N
10038
10095
        IF (A0.LT.1.0D-60) THEN
10039
10096
           DO 10 K=0,N
10199
10256
        END
10200
10257
 
10201
10258
 
10202
 
        
 
10259
 
10203
10260
C       **********************************
10204
10261
 
10205
10262
        SUBROUTINE OTHPL(KF,N,X,PL,DPL)
10323
10380
        END
10324
10381
 
10325
10382
 
10326
 
        
 
10383
 
10327
10384
C       **********************************
10328
10385
 
10329
10386
        SUBROUTINE RSWFO(M,N,C,X,CV,KF,R1F,R1D,R2F,R2D)
10376
10433
        END
10377
10434
 
10378
10435
 
10379
 
        
 
10436
 
10380
10437
C       **********************************
10381
10438
 
10382
10439
        SUBROUTINE CH12N(N,Z,NM,CHF1,CHD1,CHF2,CHD2)
10441
10498
        END
10442
10499
 
10443
10500
 
10444
 
        
 
10501
 
10445
10502
C       **********************************
10446
10503
 
10447
10504
        SUBROUTINE JYZO(N,NT,RJ0,RJ1,RY0,RY1)
10520
10577
           X=1.19477+1.08933*N
10521
10578
        ELSE
10522
10579
           X=N+0.93158*N**0.33333+0.26035/N**0.33333
10523
 
        ENDIF           
 
10580
        ENDIF
10524
10581
        L=0
10525
10582
        XGUESS=X
10526
10583
20      X0=X
10544
10601
           X=2.67257+1.16099*N
10545
10602
        ELSE
10546
10603
           X=N+1.8211*N**0.33333+0.94001/N**0.33333
10547
 
        ENDIF  
 
10604
        ENDIF
10548
10605
        L=0
10549
10606
        XGUESS=X
10550
10607
25      X0=X
10565
10622
        END
10566
10623
 
10567
10624
 
10568
 
        
 
10625
 
10569
10626
C       **********************************
10570
10627
 
10571
10628
        SUBROUTINE IKV(V,X,VM,BI,DI,BK,DK)
10714
10771
        END
10715
10772
 
10716
10773
 
10717
 
        
 
10774
 
10718
10775
C       **********************************
10719
10776
 
10720
10777
        SUBROUTINE SDMN(M,N,C,CV,KD,DF)
10742
10799
5             DF(I)=0D0
10743
10800
           DF((N-M)/2+1)=1.0D0
10744
10801
           RETURN
10745
 
        ENDIF   
 
10802
        ENDIF
10746
10803
        CS=C*C*KD
10747
10804
        IP=1
10748
10805
        K=0
10776
10833
12                  DF(K1)=DF(K1)*1.0D-100
10777
10834
                 F1=F1*1.0D-100
10778
10835
                 F0=F0*1.0D-100
10779
 
              ENDIF  
 
10836
              ENDIF
10780
10837
           ELSE
10781
10838
              KB=K
10782
10839
              FL=DF(K+1)
10788
10845
              ELSE IF (KB.EQ.2) THEN
10789
10846
                 DF(2)=F2
10790
10847
                 FS=-((D(2)-CV)*F2+G(2)*F1)/A(2)
10791
 
              ELSE 
 
10848
              ELSE
10792
10849
                 DF(2)=F2
10793
10850
                 DO 20 J=3,KB+1
10794
10851
                    F=-((D(J-1)-CV)*F2+G(J-1)*F1)/A(J-1)
10798
10855
15                        DF(K1)=DF(K1)*1.0D-100
10799
10856
                       F=F*1.0D-100
10800
10857
                       F2=F2*1.0D-100
10801
 
                    ENDIF  
 
10858
                    ENDIF
10802
10859
                    F1=F2
10803
10860
20                  F2=F
10804
10861
                 FS=F
10837
10894
 
10838
10895
 
10839
10896
 
10840
 
        
 
10897
 
10841
10898
C       **********************************
10842
10899
 
10843
10900
        SUBROUTINE AJYIK(X,VJ1,VJ2,VY1,VY2,VI1,VI2,VK1,VK2)
11168
11225
        END
11169
11226
 
11170
11227
 
11171
 
        
 
11228
 
11172
11229
C       **********************************
11173
11230
 
11174
11231
        SUBROUTINE CIKVA(V,Z,VM,CBI,CDI,CBK,CDK)
11187
11244
C                VM --- Highest order computed
11188
11245
C       Routines called:
11189
11246
C            (1) GAMMA2 for computing the gamma function
11190
 
C            (2) MSTA1 and MSTA2 for computing the starting 
 
11247
C            (2) MSTA1 and MSTA2 for computing the starting
11191
11248
C                point for backward recurrence
11192
11249
C       ============================================================
11193
11250
C
11336
11393
        END
11337
11394
 
11338
11395
 
11339
 
        
 
11396
 
11340
11397
C       **********************************
11341
11398
 
11342
11399
        SUBROUTINE CFC(Z,ZF,ZD)
11399
11456
        END
11400
11457
 
11401
11458
 
11402
 
        
 
11459
 
11403
11460
C       **********************************
11404
11461
 
11405
11462
        SUBROUTINE FCS(X,C,S)
11808
11865
        END
11809
11866
 
11810
11867
 
11811
 
        
 
11868
 
11812
11869
C       **********************************
11813
11870
 
11814
11871
        SUBROUTINE GAIH(X,GA)
11958
12015
        END
11959
12016
 
11960
12017
 
11961
 
        
 
12018
 
11962
12019
C       **********************************
11963
12020
 
11964
12021
        SUBROUTINE CLQMN(MM,M,N,X,Y,CQM,CQD)
12081
12138
           DO 5 I=1,N-M+1
12082
12139
5             EG(I)=(I+M)*(I+M-1.0D0)
12083
12140
           GO TO 70
12084
 
        ENDIF                                           
 
12141
        ENDIF
12085
12142
        ICM=(N-M+2)/2
12086
12143
        NM=10+INT(0.5*(N-M)+C)
12087
12144
        CS=C*C*KD
12339
12396
        END
12340
12397
 
12341
12398
 
12342
 
        
 
12399
 
12343
12400
C       **********************************
12344
12401
 
12345
12402
        SUBROUTINE MTU12(KF,KC,M,Q,X,F1R,D1R,F2R,D2R)
12385
12442
        ELSE
12386
12443
           QM=17.0+3.1*SQRT(Q)-.126*Q+.0037*SQRT(Q)*Q
12387
12444
        ENDIF
12388
 
        KM=INT(QM+0.5*M)              
 
12445
        KM=INT(QM+0.5*M)
12389
12446
        CALL FCOEF(KD,M,Q,A,FG)
12390
12447
        IC=INT(M/2)+1
12391
12448
        IF (KD.EQ.4) IC=M/2
12465
12522
        END
12466
12523
 
12467
12524
 
12468
 
        
 
12525
 
12469
12526
C       **********************************
12470
12527
 
12471
12528
        SUBROUTINE CIK01(Z,CBI0,CDI0,CBI1,CDI1,CBK0,CDK0,CBK1,CDK1)
12472
12529
C
12473
12530
C       ==========================================================
12474
 
C       Purpose: Compute modified Bessel functions I0(z), I1(z), 
12475
 
C                K0(z), K1(z), and their derivatives for a 
 
12531
C       Purpose: Compute modified Bessel functions I0(z), I1(z),
 
12532
C                K0(z), K1(z), and their derivatives for a
12476
12533
C                complex argument
12477
12534
C       Input :  z --- Complex argument
12478
12535
C       Output:  CBI0 --- I0(z)
12682
12739
        SY(0)=-DCOS(X)/X
12683
12740
        F0=SY(0)
12684
12741
        DY(0)=(DSIN(X)+DCOS(X)/X)/X
12685
 
        IF (N.LT.1) THEN 
 
12742
        IF (N.LT.1) THEN
12686
12743
           RETURN
12687
12744
        ENDIF
12688
12745
        SY(1)=(SY(0)-DSIN(X))/X
12690
12747
        DO 15 K=2,N
12691
12748
           F=(2.0D0*K-1.0D0)*F1/X-F0
12692
12749
           SY(K)=F
12693
 
           IF (DABS(F).GE.1.0D+300) GO TO 20              
 
12750
           IF (DABS(F).GE.1.0D+300) GO TO 20
12694
12751
           F0=F1
12695
12752
15         F1=F
12696
12753
20      NM=K-1
12700
12757
        END
12701
12758
 
12702
12759
 
12703
 
        
 
12760
 
12704
12761
C       **********************************
12705
12762
 
12706
12763
        SUBROUTINE JELP(U,HK,ESN,ECN,EDN,EPH)
12862
12919
                 BJ0=SR*(PU0*DCOS(T0)-QU0*DSIN(T0))
12863
12920
                 BJ1=SR*(PU1*DCOS(T1)-QU1*DSIN(T1))
12864
12921
C                Forward recurrence for J_|v| (Abm & Stg 9.1.27)
12865
 
C                It's OK for the limited range -8.0 ≤ v ≤ 12.5, 
 
12922
C                It's OK for the limited range -8.0 ≤ v ≤ 12.5,
12866
12923
C                since x >= 20 here; but would be unstable for v <~ -20
12867
12924
                 BF0=BJ0
12868
12925
                 BF1=BJ1