求助-复变量的贝塞尔第一类和第二类函数的Fortran子程序

jlg1206 · 发表于 2021-1-25 12:05:15

求助群里大神们，谁那里有可以求解复变量的贝塞尔第一类和第二类函数的Fortran子程序？肯求提供，感激涕零！

kerb · 发表于 2021-4-17 12:20:37

搜netlib里面的amos，里面有

风平老涡 · 发表于 2021-4-8 20:59:14

jlg1206 发表于 2021-4-8 19:16
好久没登陆了，非常感谢您的帮助。但是我运行了一下您提供的这个程序，编译的时候还是出现很多比如变量未 ...

附件应该有俩个文件。

jlg1206 · 发表于 2021-4-8 19:16:08

风平老涡发表于 2021-2-9 00:05
[mw_shl_code=fortran,true]SUBROUTINE cbsslj(z,cnu,w)
!---------------------------------------------- ...

好久没登陆了，非常感谢您的帮助。但是我运行了一下您提供的这个程序，编译的时候还是出现很多比如变量未声明这样的问题。

风平老涡 · 发表于 2021-2-9 00:06:24

[Fortran] 纯文本查看 复制代码

001FUNCTION anorm(z) RESULT(fn_val)
002! Replaces the statement function anorm in the F77 code.
003 
004COMPLEX (dp), INTENT(IN)  :: z
005REAL (dp)                 :: fn_val
006 
007fn_val = MAX( ABS( REAL(z, KIND=dp)), ABS(AIMAG(z) ) )
008RETURN
009END FUNCTION anorm
010 
011 
012 
013FUNCTION dgam1(x) RESULT(fn_val)
014!-----------------------------------------------------------------------
015!     EVALUATION OF 1/GAMMA(1 + X) - 1  FOR -0.5 <= X <= 1.5
016!-----------------------------------------------------------------------
017 
018!     THE FOLLOWING ARE THE FIRST 49 COEFFICIENTS OF THE MACLAURIN
019!     EXPANSION FOR 1/GAMMA(1 + X) - 1. THE COEFFICIENTS ARE
020!     CORRECT TO 40 DIGITS. THE COEFFICIENTS WERE OBTAINED BY
021!     ALFRED H. MORRIS JR. (NAVAL SURFACE WARFARE CENTER) AND ARE
022!     GIVEN HERE FOR REFERENCE. ONLY THE FIRST 14 COEFFICIENTS ARE
023!     USED IN THIS CODE.
024 
025!                           -----------
026 
027!     DATA A(1)  / .5772156649015328606065120900824024310422D+00/,
028!    *     A(2)  /-.6558780715202538810770195151453904812798D+00/,
029!    *     A(3)  /-.4200263503409523552900393487542981871139D-01/,
030!    *     A(4)  / .1665386113822914895017007951021052357178D+00/,
031!    *     A(5)  /-.4219773455554433674820830128918739130165D-01/,
032!    *     A(6)  /-.9621971527876973562114921672348198975363D-02/,
033!    *     A(7)  / .7218943246663099542395010340446572709905D-02/,
034!    *     A(8)  /-.1165167591859065112113971084018388666809D-02/,
035!    *     A(9)  /-.2152416741149509728157299630536478064782D-03/,
036!    *     A(10) / .1280502823881161861531986263281643233949D-03/
037!     DATA A(11) /-.2013485478078823865568939142102181838229D-04/,
038!    *     A(12) /-.1250493482142670657345359473833092242323D-05/,
039!    *     A(13) / .1133027231981695882374129620330744943324D-05/,
040!    *     A(14) /-.2056338416977607103450154130020572836513D-06/,
041!    *     A(15) / .6116095104481415817862498682855342867276D-08/,
042!    *     A(16) / .5002007644469222930055665048059991303045D-08/,
043!    *     A(17) /-.1181274570487020144588126565436505577739D-08/,
044!    *     A(18) / .1043426711691100510491540332312250191401D-09/,
045!    *     A(19) / .7782263439905071254049937311360777226068D-11/,
046!    *     A(20) /-.3696805618642205708187815878085766236571D-11/
047!     DATA A(21) / .5100370287454475979015481322863231802727D-12/,
048!    *     A(22) /-.2058326053566506783222429544855237419746D-13/,
049!    *     A(23) /-.5348122539423017982370017318727939948990D-14/,
050!    *     A(24) / .1226778628238260790158893846622422428165D-14/,
051!    *     A(25) /-.1181259301697458769513764586842297831212D-15/,
052!    *     A(26) / .1186692254751600332579777242928674071088D-17/,
053!    *     A(27) / .1412380655318031781555803947566709037086D-17/,
054!    *     A(28) /-.2298745684435370206592478580633699260285D-18/,
055!    *     A(29) / .1714406321927337433383963370267257066813D-19/,
056!    *     A(30) / .1337351730493693114864781395122268022875D-21/
057!     DATA A(31) /-.2054233551766672789325025351355733796682D-21/,
058!    *     A(32) / .2736030048607999844831509904330982014865D-22/,
059!    *     A(33) /-.1732356445910516639057428451564779799070D-23/,
060!    *     A(34) /-.2360619024499287287343450735427531007926D-25/,
061!    *     A(35) / .1864982941717294430718413161878666898946D-25/,
062!    *     A(36) /-.2218095624207197204399716913626860379732D-26/,
063!    *     A(37) / .1297781974947993668824414486330594165619D-27/,
064!    *     A(38) / .1180697474966528406222745415509971518560D-29/,
065!    *     A(39) /-.1124584349277088090293654674261439512119D-29/,
066!    *     A(40) / .1277085175140866203990206677751124647749D-30/
067!     DATA A(41) /-.7391451169615140823461289330108552823711D-32/,
068!    *     A(42) / .1134750257554215760954165259469306393009D-34/,
069!    *     A(43) / .4639134641058722029944804907952228463058D-34/,
070!    *     A(44) /-.5347336818439198875077418196709893320905D-35/,
071!    *     A(45) / .3207995923613352622861237279082794391090D-36/,
072!    *     A(46) /-.4445829736550756882101590352124643637401D-38/,
073!    *     A(47) /-.1311174518881988712901058494389922190237D-38/,
074!    *     A(48) / .1647033352543813886818259327906394145400D-39/,
075!    *     A(49) /-.1056233178503581218600561071538285049997D-40/
076 
077!                           -----------
078 
079!     C = A(1) - 1 IS ALSO FREQUENTLY NEEDED. C HAS THE VALUE ...
080 
081!     DATA C /-.4227843350984671393934879099175975689578D+00/
082 
083!-----------------------------------------------------------------------
084REAL (dp), INTENT(IN) :: x
085REAL (dp)             :: fn_val
086 
087! Local variables
088REAL (dp) :: d, t, w, z
089REAL (dp), PARAMETER :: a0 = .611609510448141581788D-08, a1  &
090        = .624730830116465516210D-08, b1 = .203610414066806987300D+00, b2  &
091        = .266205348428949217746D-01, b3 = .493944979382446875238D-03, b4  &
092        = -.851419432440314906588D-05, b5 = -.643045481779353022248D-05, b6  &
093        = .992641840672773722196D-06, b7 = -.607761895722825260739D-07, b8  &
094        = .195755836614639731882D-09
095REAL (dp), PARAMETER :: p0 = .6116095104481415817861D-08, p1  &
096        = .6871674113067198736152D-08, p2 = .6820161668496170657, p3  &
097        = .4686843322948848031080D-10, p4 = .1572833027710446286995D-11, p5  &
098        = -.1249441572276366213222D-12, p6 = .4343529937408594255178D-14, q1  &
099        = .3056961078365221025009D+00, q2 = .5464213086042296536016D-01, q3  &
100        = .4956830093825887312, q4 = .2692369466186361192876D-03
101REAL (dp), PARAMETER :: c = -.422784335098467139393487909917598D+00, c0  &
102        = .577215664901532860606512090082402D+00, c1  &
103        = -.655878071520253881077019515145390D+00, c2  &
104        = -.420026350340952355290039348754298D-01, c3  &
105        = .166538611382291489501700795102105D+00, c4  &
106        = -.421977345555443367482083012891874D-01, c5  &
107        = -.962197152787697356211492167234820D-02, c6  &
108        = .721894324666309954239501034044657D-02, c7  &
109        = -.116516759185906511211397108401839D-02, c8  &
110        = -.215241674114950972815729963053648D-03, c9  &
111        = .128050282388116186153198626328164D-03, c10  &
112        = -.201348547807882386556893914210218D-04, c11  &
113        = -.125049348214267065734535947383309D-05, c12  &
114        = .113302723198169588237412962033074D-05, c13  &
115        = -.205633841697760710345015413002057D-06
116!----------------------------
117t = x
118d = x - 0.5_dp
119IF (d > 0._dp) t = d - 0.5_dp
120IF (t < 0.0_dp) THEN
121  GO TO 30
122ELSE IF (t > 0.0_dp) THEN
123  GO TO 20
124END IF
125 
126fn_val = 0._dp
127RETURN
128!------------
129 
130!                CASE WHEN 0 < T <= 0.5
131 
132!              W IS A MINIMAX APPROXIMATION FOR
133!              THE SERIES A(15) + A(16)*T + ...
134 
135!------------
13620 w = ((((((p6*t + p5)*t + p4)*t + p3)*t + p2)*t + p1)*t + p0) /   &
137       ((((q4*t+q3)*t + q2)*t + q1)*t + 1._dp)
138z = (((((((((((((w*t + c13)*t + c12)*t + c11)*t + c10)*t + c9)*t + c8)*t + c7)*t  &
139    + c6)*t + c5)*t + c4)*t + c3)*t + c2)*t + c1) * t + c0
140 
141IF (d <= 0._dp) THEN
142  fn_val = x * z
143  RETURN
144END IF
145fn_val = (t/x) * ((z-0.5_dp)-0.5_dp)
146RETURN
147!------------
148 
149!                CASE WHEN -0.5 <= T < 0
150 
151!              W IS A MINIMAX APPROXIMATION FOR
152!              THE SERIES A(15) + A(16)*T + ...
153 
154!------------
15530 w = (a1*t + a0) / ((((((((b8*t + b7)*t + b6)*t + b5)*t + b4)*t + b3)*t + b2)*t + b1)*t+1._dp)
156z = (((((((((((((w*t + c13)*t + c12)*t + c11)*t + c10)*t + c9)*t + c8)*t + c7)*t  &
157    + c6)*t + c5)*t + c4)*t + c3)*t + c2)*t + c1) * t + c
158 
159IF (d <= 0._dp) THEN
160  fn_val = x * ((z+0.5_dp)+0.5_dp)
161  RETURN
162END IF
163fn_val = t * z / x
164RETURN
165END FUNCTION dgam1
166 
167 
168 
169FUNCTION dpdel(x) RESULT(fn_val)
170!-----------------------------------------------------------------------
171 
172!     COMPUTATION OF THE FUNCTION DEL(X) FOR  X >= 10  WHERE
173!     LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X)
174 
175!                         --------
176 
177!     THE SERIES FOR DPDEL ON THE INTERVAL 0.0 TO 1.0 DERIVED BY
178!     A.H. MORRIS FROM THE CHEBYSHEV SERIES IN THE SLATEC LIBRARY
179!     OBTAINED BY WAYNE FULLERTON (LOS ALAMOS).
180 
181!-----------------------------------------------------------------------
182REAL (dp), INTENT(IN) :: x
183REAL (dp)             :: fn_val
184 
185! Local variables
186REAL (dp), PARAMETER :: a(15) = (/ .833333333333333333333333333333D-01,  &
187        -.277777777777777777777777752282D-04,  &
188        .793650793650793650791732130419D-07,  &
189        -.595238095238095232389839236182D-09,  &
190        .841750841750832853294451671990D-11,  &
191        -.191752691751854612334149171243D-12,  &
192        .641025640510325475730918472625D-14,  &
193        -.295506514125338232839867823991D-15,  &
194        .179643716359402238723287696452D-16,  &
195        -.139228964661627791231203060395D-17,  &
196        .133802855014020915603275339093D-18,  &
197        -.154246009867966094273710216533D-19,  &
198        .197701992980957427278370133333D-20,  &
199        -.234065664793997056856992426667D-21,  &
200        .171348014966398575409015466667D-22 /)
201REAL (dp) :: t, w
202INTEGER   :: i, k
203!-----------------------------------------------------------------------
204t = (10._dp/x) ** 2
205w = a(15)
206DO i = 1, 14
207  k = 15 - i
208  w = t * w + a(k)
209END DO
210fn_val = w / x
211RETURN
212END FUNCTION dpdel
213 
214END MODULE Complex_Bessel

风平老涡 · 发表于 2021-2-9 00:05:53

[Fortran] 纯文本查看 复制代码

001SUBROUTINE cbdb(cz,cnu,fn,w)
002!-----------------------------------------------------------------------
003 
004!         CALCULATION OF J   (CZ) BY THE DEBYE APPROXIMATION
005!                         CNU
006!                         ------------------
007 
008!     IT IS ASSUMED THAT REAL(CZ) .GE. 0 AND THAT REAL(CNU) = FN + K
009!     WHERE K IS AN INTEGER.
010 
011!-----------------------------------------------------------------------
012 
013COMPLEX (dp), INTENT(IN)   :: cz, cnu
014REAL (dp), INTENT(IN)      :: fn
015COMPLEX (dp), INTENT(OUT)  :: w
016 
017! Local variables
018REAL (dp)     :: is, inu, izn
019COMPLEX (dp)  :: c1, c2, eta, nu, p, p1, q, r, s, s1, s2, sm, t, z, zn
020!----------------------
021!     C = 1/SQRT(2*PI)
022!     BND = PI/3
023!----------------------
024REAL (dp), PARAMETER  :: c = .398942280401433_dp, pi = 3.14159265358979_dp,  &
025                         pi2 = 6.28318530717959_dp, bnd = 1.04719755119660_dp
026COMPLEX (dp), PARAMETER :: j = (0.0, 1.0)
027REAL (dp)  :: alpha, am, aq, ar
028REAL (dp)  :: phi, sgn, theta
029REAL (dp)  :: u, v, x, y
030INTEGER    :: ind, k, l, m
031 
032!----------------------
033!             COEFFICIENTS OF THE FIRST 16 POLYNOMIALS
034!                   IN THE DEBYE APPROXIMATION
035!----------------------
036 
037REAL (dp)  :: a(136) = (/ 1.0_dp, -.208333333333333_dp, .125000000000000_dp, .334201388888889_dp, &
038  -.401041666666667_dp, .703125000000000D-01,-.102581259645062D+01, .184646267361111D+01, &
039  -.891210937500000_dp, .732421875000000D-01, .466958442342625D+01,-.112070026162230D+02, &
040   .878912353515625D+01,-.236408691406250D+01, .112152099609375_dp,-.282120725582002D+02, &
041   .846362176746007D+02,-.918182415432400D+02, .425349987453885D+02,-.736879435947963D+01, &
042   .227108001708984_dp, .212570130039217D+03,-.765252468141182D+03, .105999045252800D+04, &
043  -.699579627376133D+03, .218190511744212D+03,-.264914304869516D+02, .572501420974731_dp, &
044  -.191945766231841D+04, .806172218173731D+04,-.135865500064341D+05, .116553933368645D+05, &
045  -.530564697861340D+04, .120090291321635D+04,-.108090919788395D+03, .172772750258446D+01, &
046   .202042913309661D+05,-.969805983886375D+05, .192547001232532D+06,-.203400177280416D+06, &
047   .122200464983017D+06,-.411926549688976D+05, .710951430248936D+04,-.493915304773088D+03, &
048   .607404200127348D+01,-.242919187900551D+06, .131176361466298D+07,-.299801591853811D+07, &
049   .376327129765640D+07,-.281356322658653D+07, .126836527332162D+07,-.331645172484564D+06, &
050   .452187689813627D+05,-.249983048181121D+04, .243805296995561D+02, .328446985307204D+07, &
051  -.197068191184322D+08, .509526024926646D+08,-.741051482115327D+08, .663445122747290D+08, &
052  -.375671766607634D+08, .132887671664218D+08,-.278561812808645D+07, .308186404612662D+06, &
053  -.138860897537170D+05, .110017140269247D+03,-.493292536645100D+08, .325573074185766D+09, &
054  -.939462359681578D+09, .155359689957058D+10,-.162108055210834D+10, .110684281682301D+10, &
055  -.495889784275030D+09, .142062907797533D+09,-.244740627257387D+08, .224376817792245D+07, &
056  -.840054336030241D+05, .551335896122021D+03, .814789096118312D+09,-.586648149205185D+10, &
057   .186882075092958D+11,-.346320433881588D+11, .412801855797540D+11,-.330265997498007D+11, &
058   .179542137311556D+11,-.656329379261928D+10, .155927986487926D+10,-.225105661889415D+09, &
059   .173951075539782D+08,-.549842327572289D+06, .303809051092238D+04,-.146792612476956D+11, &
060   .114498237732026D+12,-.399096175224466D+12, .819218669548577D+12,-.109837515608122D+13, &
061   .100815810686538D+13,-.645364869245377D+12, .287900649906151D+12,-.878670721780233D+11, &
062   .176347306068350D+11,-.216716498322380D+10, .143157876718889D+09,-.387183344257261D+07, &
063   .182577554742932D+05, .286464035717679D+12,-.240629790002850D+13, .910934118523990D+13, &
064  -.205168994109344D+14, .305651255199353D+14,-.316670885847852D+14, .233483640445818D+14, &
065  -.123204913055983D+14, .461272578084913D+13,-.119655288019618D+13, .205914503232410D+12, &
066  -.218229277575292D+11, .124700929351271D+10,-.291883881222208D+08, .118838426256783D+06, &
067  -.601972341723401D+13, .541775107551060D+14,-.221349638702525D+15, .542739664987660D+15, &
068  -.889496939881026D+15, .102695519608276D+16,-.857461032982895D+15, .523054882578445D+15, &
069  -.232604831188940D+15, .743731229086791D+14,-.166348247248925D+14, .248500092803409D+13, &
070  -.229619372968246D+12, .114657548994482D+11,-.234557963522252D+09, .832859304016289D+06 /)
071 
072z = cz
073nu = cnu
074inu = AIMAG(cnu)
075IF (inu < 0.0D0) THEN
076  z = CONJG(z)
077  nu = CONJG(nu)
078END IF
079x = REAL(z, KIND=dp)
080y = AIMAG(z)
081 
082!          TANH(GAMMA) = SQRT(1 - (Z/NU)**2) = W/NU
083!          T = EXP(NU*(TANH(GAMMA) - GAMMA))
084 
085zn = z / nu
086izn = AIMAG(zn)
087IF (ABS(izn) <= 0.1D0*ABS(REAL(zn, KIND=dp))) THEN
088   
089  s = (1.0D0-zn) * (1.0D0+zn)
090  eta = 1.0D0 / s
091  q = SQRT(s)
092  s = 1.0D0 / (nu*q)
093  t = zn / (1.0D0 + q)
094  t = EXP(nu*(q + LOG(t)))
095ELSE
096   
097  s = (nu-z) * (nu+z)
098  eta = (nu*nu) / s
099  w = SQRT(s)
100  q = w / nu
101  IF (REAL(q, KIND=dp) < 0.0D0) w = -w
102  s = 1.0D0 / w
103  t = z / (nu+w)
104  t = EXP(w + nu*LOG(t))
105END IF
106 
107is = AIMAG(s)
108r = SQRT(s)
109c1 = r * t
110ar = REAL(r, KIND=dp) * REAL(r, KIND=dp) + AIMAG(r) * AIMAG(r)
111aq = -1.0D0 / (REAL(q, KIND=dp)*REAL(q, KIND=dp) + AIMAG(q)*AIMAG(q))
112 
113phi = ATAN2(y,x) / 3.0D0
114q = nu - z
115theta = ATAN2(AIMAG(q),REAL(q, KIND=dp)) - phi
116ind = 0
117IF (ABS(theta) >= 2.0D0*bnd) THEN
118   
119  ind = 1
120  CALL dcrec(REAL(t, KIND=dp), AIMAG(t),u,v)
121  c2 = -j * r * CMPLX(u,v, KIND=dp)
122  IF (is >= 0.0D0) THEN
123    IF (is > 0.0D0) GO TO 10
124    IF (REAL(s, KIND=dp) <= 0.0D0) GO TO 10
125  END IF
126  c2 = -c2
127END IF
128 
129!          SUMMATION OF THE SERIES S1 AND S2
130 
13110 sm = s * s
132p = (a(2)*eta + a(3)) * s
133p1 = ((a(4)*eta + a(5))*eta + a(6)) * sm
134s1 = (1.0D0 + p) + p1
135IF (ind /= 0) s2 = (1.0D0-p) + p1
136sgn = 1.0D0
137am = ar * ar
138m = 4
139l = 6
140 
141!          P = VALUE OF THE M-TH POLYNOMIAL
142 
14320 l = l + 1
144alpha = a(l)
145p = CMPLX(a(l),0.0D0, KIND=dp)
146DO  k = 2, m
147  l = l + 1
148  alpha = a(l) + aq * alpha
149  p = a(l) + eta * p
150END DO
151 
152!          ONLY THE S1 SUM IS FORMED WHEN IND = 0
153 
154sm = s * sm
155p = p * sm
156s1 = s1 + p
157IF (ind /= 0) THEN
158  sgn = -sgn
159  s2 = s2 + sgn * p
160END IF
161am = ar * am
162IF (1.0D0 + alpha*am /= 1.0D0) THEN
163  m = m + 1
164  IF (m <= 16) GO TO 20
165END IF
166 
167!          FINAL ASSEMBLY
168 
169s1 = c * c1 * s1
170IF (ind == 0) THEN
171  w = s1
172ELSE
173   
174  s2 = c * c2 * s2
175  q = nu + z
176  theta = ATAN2(AIMAG(q),REAL(q, KIND=dp)) - phi
177  IF (ABS(theta) <= bnd) THEN
178    w = s1 + s2
179  ELSE
180     
181    alpha = pi2
182    IF (izn < 0.0D0) alpha = -alpha
183    t = alpha * CMPLX(ABS(inu),-fn, KIND=dp)
184    alpha = EXP(REAL(t))
185    u = AIMAG(t)
186    r = CMPLX(COS(u),SIN(u), KIND=dp)
187    t = s1 - (alpha*r) * s1
188    IF (x == 0.0D0 .AND. inu == 0.0D0) t = -t
189     
190    IF (y < 0.0D0) THEN
191      IF (izn >= 0.0D0 .AND. theta <= SIGN(pi,theta)) s2 =  &
192                                                      s2 * ( CONJG(r)/alpha)
193      IF (x == 0.0D0) GO TO 40
194      IF (izn >= 0.0D0) THEN
195        IF (is < 0.0D0) GO TO 40
196      END IF
197    END IF
198     
199    w = s2 + t
200    GO TO 50
201    40 w = s2 - t
202  END IF
203END IF
204 
20550 IF (inu < 0.0D0) w = CONJG(w)
206RETURN
207END SUBROUTINE cbdb
208 
209 
210 
211SUBROUTINE cbja(cz,cnu,w)
212!-----------------------------------------------------------------------
213!        COMPUTATION OF J(NU,Z) BY THE ASYMPTOTIC EXPANSION
214!-----------------------------------------------------------------------
215 
216COMPLEX (dp), INTENT(IN)   :: cz
217COMPLEX (dp), INTENT(IN)   :: cnu
218COMPLEX (dp), INTENT(OUT)  :: w
219 
220! Local variables
221REAL (dp)     :: eps, inu, m
222COMPLEX (dp)  :: a, a1, arg, e, eta, nu, p, q, t, z, zr, zz
223!--------------------------
224REAL (dp) :: r, rnu, tol, u, v
225REAL (dp) :: x, y
226INTEGER   :: i, ind
227 
228!--------------------------
229!     PIHALF = PI/2
230!     C = 2*PI**(-1/2)
231!--------------------------
232REAL (dp), PARAMETER    :: pihalf = 1.5707963267949_dp, c = 1.12837916709551_dp
233COMPLEX (dp), PARAMETER :: j = (0.0_dp, 1.0_dp)
234!--------------------------
235 
236!     ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE
237!            SMALLEST NUMBER SUCH THAT 1.0 + EPS .GT. 1.0 .
238 
239eps = EPSILON(0.0_dp)
240 
241!--------------------------
242z = cz
243x = REAL(z, KIND=dp)
244y = AIMAG(z)
245nu = cnu
246ind = 0
247IF (ABS(x) <= 1.d-2*ABS(y)) THEN
248  IF (AIMAG(nu) < 0.0D0 .AND. ABS(REAL(nu)) < 1.d-2*ABS(AIMAG(nu))) THEN
249    ind = 1
250    nu = CONJG(nu)
251    z = CONJG(z)
252    y = -y
253  END IF
254END IF
255 
256IF (x < -1.d-2*y) z = -z
257zz = z + z
258CALL dcrec(REAL(zz, KIND=dp),AIMAG(zz),u,v)
259zr = CMPLX(u,v, KIND=dp)
260eta = -zr * zr
261 
262p = (0.0D0,0.0D0)
263q = (0.0D0,0.0D0)
264a1 = nu * nu - 0.25D0
265a = a1
266t = a1
267m = 1.0D0
268tol = eps * anorm(a1)
269DO  i = 1, 16
270  a = a - 2.0D0 * m
271  m = m + 1.0D0
272  t = t * a * eta / m
273  p = p + t
274  a = a - 2.0D0 * m
275  m = m + 1.0D0
276  t = t * a / m
277  q = q + t
278  IF (anorm(t) <= tol) GO TO 20
279END DO
280 
28120 p = p + 1.0D0
282q = (q+a1) * zr
283w = z - pihalf * nu
284IF (ABS(AIMAG(w)) <= 1.0D0) THEN
285  arg = w - 0.5D0 * pihalf
286  w = c * SQRT(zr) * (p*COS(arg) - q*SIN(arg))
287ELSE
288  e = EXP(-j*w)
289  t = q - j * p
290  IF (AIMAG(z) > 0.0D0 .AND. REAL(z, KIND=dp) <= 1.d-2*AIMAG(z).AND.  &
291      ABS(REAL(nu, KIND=dp)) < 1.d-2*AIMAG(nu)) t = 0.5D0 * t
292  CALL dcrec(REAL(e, KIND=dp),AIMAG(e),u,v)
293  w = 0.5D0 * c * SQRT(j*zr) * ((p-j*q)*e + t*CMPLX(u,v, KIND=dp))
294END IF
295 
296IF (x < -1.d-2*y) THEN
297  IF (y < 0.0D0) nu = -nu
298   
299!     COMPUTATION OF EXP(I*PI*NU)
300   
301  rnu = REAL(nu, KIND=dp)
302  inu = AIMAG(nu)
303  r = EXP(-2.0D0*pihalf*inu)
304  u = r * dcos1(rnu)
305  v = r * dsin1(rnu)
306  w = w * CMPLX(u,v, KIND=dp)
307END IF
308 
309IF (ind /= 0) w = CONJG(w)
310RETURN
311END SUBROUTINE cbja

风平老涡 · 发表于 2021-2-9 00:05:17

[Fortran] 纯文本查看 复制代码

001SUBROUTINE cbsslj(z,cnu,w)
002!-----------------------------------------------------------------------
003 
004!         EVALUATION OF THE COMPLEX BESSEL FUNCTION J   (Z)
005!                                                    CNU
006!-----------------------------------------------------------------------
007 
008!     WRITTEN BY
009!         ANDREW H. VAN TUYL AND ALFRED H. MORRIS, JR.
010!         NAVAL SURFACE WARFARE CENTER
011!         OCTOBER, 1991
012 
013!     A MODIFICATION OF THE PROCEDURE DEVELOPED BY ALLEN V. HERSHEY
014!     (NAVAL SURFACE WARFARE CENTER) IN 1978 FOR HANDLING THE DEBYE
015!     APPROXIMATION IS EMPLOYED.
016 
017!-----------------------------------------------------------------------
018 
019COMPLEX (dp), INTENT(IN)   :: z
020COMPLEX (dp), INTENT(IN)   :: cnu
021COMPLEX (dp), INTENT(OUT)  :: w
022 
023COMPLEX (dp)  :: c, nu, s, sm1, sm2, t, tsc, w0, w1, zn, zz
024!-----------------------
025REAL (dp) :: a, cn1, cn2, e, fn
026REAL (dp) :: pn, qm, qn, qnp1
027REAL (dp) :: r, rn2, r2, sn, t1, t2
028REAL (dp) :: u, v, x, y
029INTEGER   :: i, k, m, n
030REAL (dp), PARAMETER  :: pi = 3.141592653589793238462643383279502884197_dp
031!-----------------------
032x = REAL(z, KIND=dp)
033y = AIMAG(z)
034r = cpabs(x,y)
035cn1 = REAL(cnu, KIND=dp)
036cn2 = AIMAG(cnu)
037rn2 = cn1 * cn1 + cn2 * cn2
038pn = INT(cn1)
039fn = cn1 - pn
040sn = 1.0_dp
041 
042!          CALCULATION WHEN ORDER IS AN INTEGER
043 
044IF (fn == 0.0_dp .AND. cn2 == 0.0_dp) THEN
045  n = pn
046  pn = ABS(pn)
047  cn1 = pn
048  IF (n < 0 .AND. n /= (n/2)*2) sn = -1.0_dp
049END IF
050 
051!          SELECTION OF METHOD
052 
053IF (r > 17.5D0) THEN
054  IF (r > 17.5D0 + 0.5D0*rn2) GO TO 10
055  GO TO 20
056END IF
057 
058!          USE MACLAURIN EXPANSION AND RECURSION
059 
060IF (cn1 < 0.0D0) THEN
061  qn = -1.25D0 * (r + 0.5D0*ABS(cn2) - ABS(y-0.5D0*cn2))
062  IF (cn1 < qn) THEN
063    qn = 1.25D0 * (r - MAX(1.2D0*r,ABS(y-cn2)))
064    IF (cn1 < qn) THEN
065      qn = MIN(pn, REAL(-INT(1.25D0*(r-ABS(cn2))), KIND=dp))
066      GO TO 60
067    END IF
068  END IF
069END IF
070 
071r2 = r * r
072qm = 0.0625D0 * r2 * r2 - cn2 * cn2
073qn = MAX(pn, REAL(INT(SQRT(MAX(0.0D0,qm))), KIND=dp))
074GO TO 60
075 
076!          USE ASYMPTOTIC EXPANSION
077 
07810 CALL cbja(z,cnu,w)
079RETURN
080 
081!          CALCULATION FOR 17.5 < ABS(Z) <= 17.5 + 0.5*ABS(CNU)**2
082 
08320 n = 0
084IF (ABS(cn2) < 0.8D0*ABS(y)) THEN
085  qm = -1.25D0 * (r + 0.5D0*ABS(cn2) - ABS(y-0.5D0*cn2))
086  IF (cn1 < qm) THEN
087    qm = 1.25D0 * (r - MAX(1.2D0*r, ABS(y-cn2)))
088    IF (cn1 < qm) n = 1
089  END IF
090END IF
091 
092qn = pn
093a = 4.d-3 * r * r
094zz = z
095IF (x < 0.0D0) zz = -z
096 
097!          CALCULATION OF ZONE OF EXCLUSION OF DEBYE APPROXIMATION
098 
09930 nu = CMPLX(qn+fn,cn2, KIND=dp)
100zn = nu / z
101t2 = AIMAG(zn) * AIMAG(zn)
102u = 1.0D0 - REAL(zn, KIND=dp)
103t1 = u * u + t2
104u = 1.0D0 + DBLE(zn)
105t2 = u * u + t2
106u = t1 * t2
107v = a * u / (t1*t1 + t2*t2)
108IF (u*v*v <= 1.0D0) THEN
109   
110!          THE ARGUMENT LIES INSIDE THE ZONE OF EXCLUSION
111   
112  qn = qn + 1.0D0
113  IF (n == 0) GO TO 30
114   
115!          USE MACLAURIN EXPANSION WITH FORWARD RECURRENCE
116   
117  qn = MIN(pn, REAL(-INT(1.25D0*(r-ABS(cn2))), KIND=dp))
118ELSE
119   
120!          USE BACKWARD RECURRENCE STARTING FROM THE ASYMPTOTIC EXPANSION
121   
122  qnp1 = qn + 1.0D0
123  IF (ABS(qn) < ABS(pn)) THEN
124    IF (r >= 17.5D0 + 0.5D0*(qnp1*qnp1 + cn2*cn2)) THEN
125       
126      nu = CMPLX(qn+fn,cn2, KIND=dp)
127      CALL cbja(zz,nu,sm1)
128      nu = CMPLX(qnp1+fn,cn2, KIND=dp)
129      CALL cbja(zz,nu,sm2)
130      GO TO 40
131    END IF
132  END IF
133   
134!          USE BACKWARD RECURRENCE STARTING FROM THE DEBYE APPROXIMATION
135   
136  nu = CMPLX(qn+fn,cn2, KIND=dp)
137  CALL cbdb(zz,nu,fn,sm1)
138  IF (qn == pn) GO TO 50
139  nu = CMPLX(qnp1+fn,cn2, KIND=dp)
140  CALL cbdb(zz,nu,fn,sm2)
141   
142  40 nu = CMPLX(qn+fn,cn2, KIND=dp)
143  tsc = 2.0D0 * nu * sm1 / zz - sm2
144  sm2 = sm1
145  sm1 = tsc
146  qn = qn - 1.0D0
147  IF (qn /= pn) GO TO 40
148   
149  50 w = sm1
150  IF (sn < 0.0D0) w = -w
151  IF (x >= 0.0D0) RETURN
152   
153  nu = pi * CMPLX(-cn2,cn1, KIND=dp)
154  IF (y < 0.0D0) nu = -nu
155  w = EXP(nu) * w
156  RETURN
157END IF
158 
159!          USE MACLAURIN EXPANSION WITH FORWARD OR BACKWARD RECURRENCE.
160 
16160 m = qn - pn
162IF (ABS(m) <= 1) THEN
163  nu = CMPLX(cn1,cn2, KIND=dp)
164  CALL cbjm(z,nu,w)
165ELSE
166  nu = CMPLX(qn+fn,cn2, KIND=dp)
167  CALL cbjm(z,nu,w1)
168  w0 = 0.25D0 * z * z
169  IF (m <= 0) THEN
170     
171!          FORWARD RECURRENCE
172     
173    m = ABS(m)
174    nu = nu + 1.0D0
175    CALL cbjm(z,nu,w)
176    DO  i = 2, m
177      c = nu * (nu+1.0D0)
178      t = (c/w0) * (w-w1)
179      w1 = w
180      w = t
181      nu = nu + 1.0D0
182    END DO
183  ELSE
184     
185!          BACKWARD RECURRENCE
186     
187    nu = nu - 1.0D0
188    CALL cbjm(z,nu,w)
189    DO  i = 2, m
190      c = nu * (nu+1.0D0)
191      t = (w0/c) * w1
192      w1 = w
193      w = w - t
194      nu = nu - 1.0D0
195    END DO
196  END IF
197END IF
198 
199!          FINAL ASSEMBLY
200 
201IF (fn == 0.0D0 .AND. cn2 == 0.0D0) THEN
202  k = pn
203  IF (k == 0) RETURN
204  e = sn / dgamma(pn+1.0D0)
205  w = e * w * (0.5D0*z) ** k
206  RETURN
207END IF
208 
209s = cnu * LOG(0.5D0*z)
210w = EXP(s) * w
211IF (rn2 <= 0.81D0) THEN
212  w = w * cgam0(cnu)
213  RETURN
214END IF
215CALL dcgama(0,cnu,t)
216w = cdiv(w,cnu*t)
217RETURN
218END SUBROUTINE cbsslj
219 
220 
221 
222SUBROUTINE cbjm(z,cnu,w)
223!-----------------------------------------------------------------------
224 
225!       COMPUTATION OF  (Z/2)**(-CNU) * GAMMA(CNU + 1) * J(CNU,Z)
226 
227!                           -----------------
228 
229!     THE MACLAURIN EXPANSION IS USED.  IT IS ASSUMED THAT CNU IS NOT
230!     A NEGATIVE INTEGER.
231 
232!-----------------------------------------------------------------------
233 
234COMPLEX (dp), INTENT(IN)   :: z
235COMPLEX (dp), INTENT(IN)   :: cnu
236COMPLEX (dp), INTENT(OUT)  :: w
237 
238COMPLEX (dp)  :: nu, nup1, p, s, sn, t, ti
239!--------------------------
240REAL (dp)  :: a, a0, eps, inu, m, rnu
241INTEGER    :: i, imin, k, km1, km2
242 
243!--------------------------
244 
245!     ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE
246!            SMALLEST NUMBER SUCH THAT 1.0 + EPS .GT. 1.0 .
247 
248eps = EPSILON(0.0_dp)
249 
250!--------------------------
251s = -0.25D0 * (z*z)
252nu = cnu
253rnu = REAL(nu, KIND=dp)
254inu = AIMAG(nu)
255a = 0.5D0 + (0.5D0+rnu)
256nup1 = CMPLX(a,inu, KIND=dp)
257 
258IF (a > 0.0D0) THEN
259  m = 1.0D0
260  t = s / nup1
261  w = 1.0D0 + t
262ELSE
263   
264!     ADD 1.0 AND THE FIRST K-1 TERMS
265   
266  k = INT(-a) + 2
267  km1 = k - 1
268  w = (1.0D0,0.0D0)
269  t = w
270  DO  i = 1, km1
271    m = i
272    t = t * (s/(m*(nu+m)))
273    w = w + t
274    IF (anorm(t) <= eps*anorm(w)) GO TO 20
275  END DO
276  GO TO 50
277   
278!     CHECK IF THE (K-1)-ST AND K-TH TERMS CAN BE IGNORED.
279!     IF SO THEN THE SUMMATION IS COMPLETE.
280   
281  20 IF (i /= km1) THEN
282    imin = i + 1
283    IF (imin < k-5) THEN
284      ti = t
285       
286      m = km1
287      t = s / (nu+m)
288      a0 = anorm(t) / m
289      t = t * (s/(nu+(m+1.0D0)))
290      a = anorm(t) / (m*(m+1.0D0))
291      a = MAX(a,a0)
292       
293      t = (1.0D0,0.0D0)
294      km2 = k - 2
295      DO  i = imin, km2
296        m = i
297        t = t * (s/(m*(nu+m)))
298        IF (a*anorm(t) < 0.5D0) RETURN
299      END DO
300      t = t * ti
301      imin = km2
302    END IF
303     
304!     ADD THE (K-1)-ST TERM
305     
306    a = 1.0D0
307    p = (1.0D0,0.0D0)
308    sn = p
309    DO  i = imin, km1
310      m = i
311      a = a * m
312      p = p * (nu+m)
313      sn = s * sn
314    END DO
315    t = t * (cdiv(sn,p)/a)
316    w = w + t
317  END IF
318END IF
319 
320!     ADD THE REMAINING TERMS
321 
32250 m = m + 1.0D0
323t = t * (s/(m*(nu+m)))
324w = w + t
325IF (anorm(t) > eps*anorm(w)) GO TO 50
326 
327RETURN
328END SUBROUTINE cbjm

风平老涡 · 发表于 2021-2-9 00:04:31

[Fortran] 纯文本查看 复制代码

001FUNCTION cgam0(z) RESULT(fn_val)
002!-----------------------------------------------------------------------
003!          EVALUATION OF 1/GAMMA(1 + Z)  FOR ABS(Z) < 1.0
004!-----------------------------------------------------------------------
005 
006COMPLEX (dp), INTENT(IN)  :: z
007COMPLEX (dp)              :: fn_val
008 
009COMPLEX (dp)  :: w
010INTEGER       :: i, k, n
011!-----------------------
012REAL (dp)  :: x, y
013REAL, PARAMETER :: a(25) = (/ .577215664901533_dp, -.655878071520254_dp,  &
014        -.420026350340952D-01, .166538611382291_dp, -.421977345555443D-01,  &
015        -.962197152787697D-02, .721894324666310D-02, -.116516759185907D-02,  &
016        -.215241674114951D-03, .128050282388116D-03, -.201348547807882D-04,  &
017        -.125049348214267D-0, .113302723198170D-05, -.205633841697761D-0,  &
018        .611609510448142D-08, .500200764446922D-08, -.118127457048702D-08,  &
019        .104342671169110D-09, .778226343990507D-11, -.369680561864221D-11,  &
020        .510037028745448D-12, -.205832605356651D-13, -.534812253942302D-14,  &
021        .122677862823826D-14, -.118125930169746D-15 /)
022!-----------------------
023n = 25
024x = REAL(z, KIND=dp)
025y = AIMAG(z)
026IF (x*x+y*y <= 0.04D0) n = 14
027 
028k = n
029w = a(n)
030DO  i = 2, n
031  k = k - 1
032  w = a(k) + z * w
033END DO
034fn_val = 1.0D0 + z * w
035RETURN
036END FUNCTION cgam0
037 
038 
039 
040FUNCTION dgamma(a) RESULT(fn_val)
041!-----------------------------------------------------------------------
042 
043!                EVALUATION OF THE GAMMA FUNCTION FOR
044!                     REAL (dp) ARGUMENTS
045 
046!                           -----------
047 
048!     DGAMMA(A) IS ASSIGNED THE VALUE 0 WHEN THE GAMMA FUNCTION CANNOT
049!     BE COMPUTED.
050 
051!-----------------------------------------------------------------------
052!     WRITTEN BY ALFRED H. MORRIS, JR.
053!          NAVAL SURFACE WEAPONS CENTER
054!          DAHLGREN, VIRGINIA
055!-----------------------------------------------------------------------
056REAL (dp), INTENT(IN) :: a
057REAL (dp)             :: fn_val
058 
059! Local variables
060REAL (dp), PARAMETER :: d = 0.41893853320467274178032973640562_dp,  &
061                        pi = 3.14159265358979323846264338327950_dp
062REAL (dp) :: s, t, x, w
063INTEGER   :: j, n
064!-----------------------------------------------------------------------
065!     D = 0.5*(LN(2*PI) - 1)
066!-----------------------------------------------------------------------
067fn_val = 0.0_dp
068x = a
069IF (ABS(a) <= 20._dp) THEN
070!-----------------------------------------------------------------------
071!             EVALUATION OF DGAMMA(A) FOR ABS(A) <= 20
072!-----------------------------------------------------------------------
073  t = 1.0_dp
074  n = x
075  n = n - 1
076 
077!     LET T BE THE PRODUCT OF A-J WHEN A >= 2
078 
079  IF (n < 0) THEN
080    GO TO 40
081  ELSE IF (n == 0) THEN
082    GO TO 30
083  END IF
084 
085  DO j = 1, n
086    x = x - 1._dp
087    t = x * t
088  END DO
089  30 x = x - 1._dp
090  GO TO 60
091 
092!     LET T BE THE PRODUCT OF A+J WHEN A < 1
093 
094  40 t = a
095  IF (a <= 0._dp) THEN
096    n = -n - 1
097    IF (n /= 0) THEN
098      DO j = 1, n
099        x = x + 1._dp
100        t = x * t
101      END DO
102    END IF
103    x = (x+0.5_dp) + 0.5_dp
104    t = x * t
105    IF (t == 0._dp) RETURN
106  END IF
107 
108!     THE FOLLOWING CODE CHECKS IF 1/T CAN OVERFLOW. THIS
109!     CODE MAY BE OMITTED IF DESIRED.
110 
111  IF (ABS(t) < 1.d-33) THEN
112    IF (ABS(t)*HUGE(1.0_dp) <= 1.000000001_dp) RETURN
113    fn_val = 1._dp / t
114    RETURN
115  END IF
116 
117!     COMPUTE DGAMMA(1 + X) FOR 0 <= X < 1
118 
119  60 fn_val = 1._dp / (1._dp + dgam1(x))
120 
121!     TERMINATION
122 
123  IF (a >= 1._dp) THEN
124    fn_val = fn_val * t
125    RETURN
126  END IF
127  fn_val = fn_val / t
128  RETURN
129END IF
130!-----------------------------------------------------------------------
131!           EVALUATION OF DGAMMA(A) FOR ABS(A) > 20
132!-----------------------------------------------------------------------
133IF (ABS(a) >= 1.d3) RETURN
134IF (a <= 0._dp) THEN
135  s = dsin1(a) / pi
136  IF (s == 0._dp) RETURN
137  x = -a
138END IF
139 
140!     COMPUTE THE MODIFIED ASYMPTOTIC SUM
141 
142w = dpdel(x)
143 
144!     FINAL ASSEMBLY
145 
146w = (d+w) + (x-0.5_dp) * (LOG(x)-1._dp)
147IF (w > dxparg(0)) RETURN
148fn_val = EXP(w)
149IF (a < 0._dp) fn_val = (1._dp/(fn_val*s)) / x
150 
151RETURN
152END FUNCTION dgamma
153 
154 
155 
156FUNCTION glog(x) RESULT(fn_val)
157!     -------------------
158!     EVALUATION OF LN(X) FOR X >= 15
159!     -------------------
160REAL (dp), INTENT(IN) :: x
161REAL (dp)             :: fn_val
162 
163! Local variables
164REAL (dp), PARAMETER :: c1 = .286228750476730_dp, c2 = .399999628131494+dp,  &
165                        c3 = .666666666752663_dp
166REAL (dp), PARAMETER :: w(163) = (/ .270805020110221007D+01, .277258872223978124D+01,  &
167  .283321334405621608D+01, .289037175789616469D+01, .294443897916644046D+01, .299573227355399099D+01, .304452243772342300D+01, &
168  .309104245335831585D+01, .313549421592914969D+01, .317805383034794562D+01, .321887582486820075D+01, .325809653802148205D+01, &
169  .329583686600432907D+01, .333220451017520392D+01, .336729582998647403D+01, .340119738166215538D+01, .343398720448514625D+01, &
170  .346573590279972655D+01, .349650756146648024D+01, .352636052461616139D+01, .355534806148941368D+01, .358351893845611000D+01, &
171  .361091791264422444D+01, .363758615972638577D+01, .366356164612964643D+01, .368887945411393630D+01, .371357206670430780D+01, &
172  .373766961828336831D+01, .376120011569356242D+01, .378418963391826116D+01, .380666248977031976D+01, .382864139648909500D+01, &
173  .385014760171005859D+01, .387120101090789093D+01, .389182029811062661D+01, .391202300542814606D+01, .393182563272432577D+01, &
174  .395124371858142735D+01, .397029191355212183D+01, .398898404656427438D+01, .400733318523247092D+01, .402535169073514923D+01, &
175  .404305126783455015D+01, .406044301054641934D+01, .407753744390571945D+01, .409434456222210068D+01, .411087386417331125D+01, &
176  .412713438504509156D+01, .414313472639153269D+01, .415888308335967186D+01, .417438726989563711D+01, .418965474202642554D+01, &
177  .420469261939096606D+01, .421950770517610670D+01, .423410650459725938D+01, .424849524204935899D+01, .426267987704131542D+01, &
178  .427666611901605531D+01, .429045944114839113D+01, .430406509320416975D+01, .431748811353631044D+01, .433073334028633108D+01, &
179  .434380542185368385D+01, .435670882668959174D+01, .436944785246702149D+01, .438202663467388161D+01, .439444915467243877D+01, &
180  .440671924726425311D+01, .441884060779659792D+01, .443081679884331362D+01, .444265125649031645D+01, .445434729625350773D+01, &
181  .446590811865458372D+01, .447733681447820647D+01, .448863636973213984D+01, .449980967033026507D+01, .451085950651685004D+01, &
182  .452178857704904031D+01, .453259949315325594D+01, .454329478227000390D+01, .455387689160054083D+01, .456434819146783624D+01, &
183  .457471097850338282D+01, .458496747867057192D+01, .459511985013458993D+01, .460517018598809137D+01, .461512051684125945D+01, &
184  .462497281328427108D+01, .463472898822963577D+01, .464439089914137266D+01, .465396035015752337D+01, .466343909411206714D+01, &
185  .467282883446190617D+01, .468213122712421969D+01, .469134788222914370D+01, .470048036579241623D+01, .470953020131233414D+01, &
186  .471849887129509454D+01, .472738781871234057D+01, .473619844839449546D+01, .474493212836325007D+01, .475359019110636465D+01, &
187  .476217393479775612D+01, .477068462446566476D+01, .477912349311152939D+01, .478749174278204599D+01, .479579054559674109D+01, &
188  .480402104473325656D+01, .481218435537241750D+01, .482028156560503686D+01, .482831373730230112D+01, .483628190695147800D+01, &
189  .484418708645859127D+01, .485203026391961717D+01, .485981240436167211D+01, .486753445045558242D+01, .487519732320115154D+01, &
190  .488280192258637085D+01, .489034912822175377D+01, .489783979995091137D+01, .490527477843842945D+01, .491265488573605201D+01, &
191  .491998092582812492D+01, .492725368515720469D+01, .493447393313069176D+01, .494164242260930430D+01, .494875989037816828D+01, &
192  .495582705760126073D+01, .496284463025990728D+01, .496981329957600062D+01, .497673374242057440D+01, .498360662170833644D+01, &
193  .499043258677873630D+01, .499721227376411506D+01, .500394630594545914D+01, .501063529409625575D+01, .501727983681492433D+01, &
194  .502388052084627639D+01, .503043792139243546D+01, .503695260241362916D+01, .504342511691924662D+01, .504985600724953705D+01, &
195  .505624580534830806D+01, .506259503302696680D+01, .506890420222023153D+01, .507517381523382692D+01, .508140436498446300D+01, &
196  .508759633523238407D+01, .509375020080676233D+01, .509986642782419842D+01, .510594547390058061D+01, .511198778835654323D+01, &
197  .511799381241675511D+01, .512396397940325892D+01, .512989871492307347D+01, .513579843705026176D+01, .514166355650265984D+01, &
198  .514749447681345304D+01, .515329159449777895D+01, .515905529921452903D+01, .516478597392351405D+01, .517048399503815178D+01, &
199  .517614973257382914D+01 /)
200REAL (dp) :: t, t2, z
201INTEGER   :: n
202!     -------------------
203 
204IF (x < 178.0_dp) THEN
205  n = x
206  t = (x-n) / (x+n)
207  t2 = t * t
208  z = (((c1*t2 + c2)*t2 + c3)*t2 + 2.0) * t
209  fn_val = w(n-14) + z
210  RETURN
211END IF
212 
213fn_val = LOG(x)
214RETURN
215END FUNCTION glog

风平老涡 · 发表于 2021-2-9 00:02:57

[Fortran] 纯文本查看 复制代码

001FUNCTION drexp(x) RESULT(fn_val)
002!-----------------------------------------------------------------------
003!            EVALUATION OF THE FUNCTION EXP(X) - 1
004!-----------------------------------------------------------------------
005REAL (dp), INTENT(IN)  :: x
006REAL (dp)              :: fn_val
007 
008! Local variables
009REAL (dp) :: e, w, z
010REAL (dp) :: a0 = .248015873015873015873016D-04,   &
011    a1 = -.344452080605731005808147D-05, a2 = .206664230430046597475413D-06,  &
012    a3 = -.447300111094328162971036D-08, a4 = .114734027080634968083920D-11,  &
013    b1 = -.249994190011341852652396D+00, b2 = .249987228833107957725728D-01,  &
014    b3 = -.119037506846942249362528D-02, b4 = .228908693387350391768682D-04
015REAL (dp) :: c1 = .1666666666666666666666666666666667D+00,   &
016             c2 = .4166666666666666666666666666666667D-01,   &
017             c3 = .8333333333333333333333333333333333D-02,   &
018             c4 = .1388888888888888888888888888888889D-02,   &
019             c5 = .1984126984126984126984126984126984D-03
020!---------------------------
021IF (ABS(x) <= 0.15_dp) THEN
022 
023!        Z IS A MINIMAX APPROXIMATION OF THE SERIES
024 
025!                C6 + C7*X + C8*X**2 + ....
026 
027!        THIS APPROXIMATION IS ACCURATE TO WITHIN
028!        1 UNIT OF THE 23-RD SIGNIFICANT DIGIT.
029!        THE RESULTING VALUE FOR W IS ACCURATE TO
030!        WITHIN 1 UNIT OF THE 33-RD SIGNIFICANT DIGIT.
031 
032  z = ((((a4*x + a3)*x + a2)*x + a1)*x + a0) /  &
033      ((((b4*x + b3)*x + b2)*x + b1)*x + 1._dp)
034  w = ((((((z*x + c5)*x + c4)*x + c3)*x + c2)*x + c1)*x + 0.5_dp)*x + 1._dp
035  fn_val = x * w
036  RETURN
037END IF
038 
039IF (x >= 0.0_dp) THEN
040  e = EXP(x)
041  fn_val = e * (0.5_dp + (0.5_dp - 1.0_dp/e))
042  RETURN
043END IF
044IF (x >= -77._dp) THEN
045  fn_val = (EXP(x)-0.5_dp) - 0.5_dp
046  RETURN
047END IF
048fn_val = -1._dp
049RETURN
050END FUNCTION drexp
051 
052 
053 
054SUBROUTINE dcgama(mo, z, w)
055!-----------------------------------------------------------------------
056 
057!        EVALUATION OF THE COMPLEX GAMMA AND LOGGAMMA FUNCTIONS
058 
059!                        ---------------
060 
061!     MO IS AN INTEGER.  Z AND W ARE INTERPRETED AS REAL (dp)
062!     COMPLEX NUMBERS.  IT IS ASSUMED THAT Z(1) AND Z(2) ARE THE REAL
063!     AND IMAGINARY PARTS OF THE COMPLEX NUMBER Z, AND THAT W(1) AND
064!     W(2) ARE THE REAL AND IMAGINARY PARTS OF W.
065 
066!                 W = GAMMA(Z)       IF MO = 0
067!                 W = LN(GAMMA(Z))   OTHERWISE
068 
069!-----------------------------------------------------------------------
070INTEGER, INTENT(IN)        :: mo
071COMPLEX (dp), INTENT(IN)   :: z
072COMPLEX (dp), INTENT(OUT)  :: w
073 
074! Local variables
075REAL (dp), PARAMETER :: c0(30)  &
076        = (/ .8333333333333333333333333333333333333333D-01,  &
077        -.2777777777777777777777777777777777777778D-02,  &
078         .7936507936507936507936507936507936507937D-03,  &
079        -.5952380952380952380952380952380952380952D-03,  &
080         .8417508417508417508417508417508417508418D-03,  &
081        -.1917526917526917526917526917526917526918D-02,  &
082         .6410256410256410256410256410256410256410D-02,  &
083        -.2955065359477124183006535947712418300654D-01,  &
084         .1796443723688305731649384900158893966944D+00,  &
085        -.1392432216905901116427432216905901116427D+01,  &
086         .1340286404416839199447895100069013112491D+02,  &
087        -.1568482846260020173063651324520889738281D+03,  &
088         .2193103333333333333333333333333333333333D+04,  &
089        -.3610877125372498935717326521924223073648D+05,  &
090         .6914722688513130671083952507756734675533D+06,  &
091        -.1523822153940741619228336495888678051866D+08,  &
092         .3829007513914141414141414141414141414141D+09,  &
093        -.1088226603578439108901514916552510537473D+11,  &
094         .3473202837650022522522522522522522522523D+12,  &
095        -.1236960214226927445425171034927132488108D+14,  &
096         .4887880647930793350758151625180229021085D+15,  &
097        -.2132033396091937389697505898213683855747D+17,  &
098         .1021775296525700077565287628053585500394D+19,  &
099        -.5357547217330020361082770919196920448485D+20,  &
100         .3061578263704883415043151051329622758194D+22,  &
101        -.1899991742639920405029371429306942902947D+24,  &
102         .1276337403382883414923495137769782597654D+26,  &
103        -.9252847176120416307230242348347622779519D+27,  &
104         .7218822595185610297836050187301637922490D+29,  &
105        -.6045183405995856967743148238754547286066D+31 /),  &
106        dlpi = 1.144729885849400174143427351353058711647_dp,  &
107        hl2p =  .9189385332046727417803297364056176398614_dp,  &
108        pi = 3.141592653589793238462643383279502884197_dp,  &
109        pi2 = 6.283185307179586476925286766559005768394_dp
110REAL (dp) :: a, a1, a2, c, cn, cut, d, eps, et, e2t, h1, h2, q1, q2, s, sn,  &
111             s1, s2, t, t1, t2, u, u1, u2, v1, v2, w1, w2, x, y, y2
112INTEGER   :: j, k, max, n, nm1
113!---------------------------
114!     DLPI = LOG(PI)
115!     HL2P = 0.5 * LOG(2*PI)
116!---------------------------
117 
118!     ****** MAX AND EPS ARE MACHINE DEPENDENT CONSTANTS.
119!            MAX IS THE LARGEST POSITIVE INTEGER THAT MAY
120!            BE USED, AND EPS IS THE SMALLEST NUMBER SUCH
121!            THAT  1._dp + EPS > 1._dp.
122 
123!                      MAX = IPMPAR(3)
124max = HUGE(3)
125eps = EPSILON(1.0_dp)
126 
127!---------------------------
128x = REAL(z, KIND=dp)
129y = AIMAG(z)
130IF (x < 0._dp) THEN
131!-----------------------------------------------------------------------
132!            CASE WHEN THE REAL PART OF Z IS NEGATIVE
133!-----------------------------------------------------------------------
134  y = ABS(y)
135  t = -pi * y
136  et = EXP(t)
137  e2t = et * et
138 
139!     SET  A1 = (1 + E2T)/2  AND  A2 = (1 - E2T)/2
140 
141  a1 = 0.5_dp * (1._dp+e2t)
142  t2 = t + t
143  IF (t2 >= -0.15_dp) THEN
144    a2 = -0.5_dp * drexp(t2)
145  ELSE
146    a2 = 0.5_dp * (0.5_dp+(0.5_dp-e2t))
147  END IF
148 
149!     COMPUTE SIN(PI*X) AND COS(PI*X)
150 
151  u = MAX
152  IF (ABS(x) >= MIN(u,1._dp/eps)) GO TO 80
153  k = ABS(x)
154  u = x + k
155  k = MOD(k,2)
156  IF (u <= -0.5_dp) THEN
157    u = 0.5_dp + (0.5_dp+u)
158    k = k + 1
159  END IF
160  u = pi * u
161  sn = SIN(u)
162  cn = COS(u)
163  IF (k == 1) THEN
164    sn = -sn
165    cn = -cn
166  END IF
167 
168!     SET  H1 + H2*I  TO  PI/SIN(PI*Z)  OR  LOG(PI/SIN(PI*Z))
169 
170  a1 = sn * a1
171  a2 = cn * a2
172  a = a1 * a1 + a2 * a2
173  IF (a == 0._dp) GO TO 80
174  IF (mo == 0) THEN
175 
176    h1 = a1 / a
177    h2 = -a2 / a
178    c = pi * et
179    h1 = c * h1
180    h2 = c * h2
181  ELSE
182 
183    h1 = (dlpi+t) - 0.5_dp * LOG(a)
184    h2 = -ATAN2(a2,a1)
185  END IF
186  IF (AIMAG(z) >= 0._dp) THEN
187    x = 1.0 - x
188    y = -y
189  ELSE
190    h2 = -h2
191    x = 1.0 - x
192  END IF
193END IF
194!-----------------------------------------------------------------------
195!           CASE WHEN THE REAL PART OF Z IS NONNEGATIVE
196!-----------------------------------------------------------------------
197w1 = 0._dp
198w2 = 0._dp
199n = 0
200t = x
201y2 = y * y
202a = t * t + y2
203cut = 225._dp
204IF (eps > 1.d-30) cut = 144._dp
205IF (eps > 1.d-20) cut = 64._dp
206IF (a < cut) THEN
207  IF (a == 0._dp) GO TO 80
208  10 n = n + 1
209  t = t + 1._dp
210  a = t * t + y2
211  IF (a < cut) GO TO 10
212 
213!     LET S1 + S2*I BE THE PRODUCT OF THE TERMS (Z+J)/(Z+N)
214 
215  u1 = (x*t+y2) / a
216  u2 = y / a
217  s1 = u1
218  s2 = n * u2
219  IF (n >= 2) THEN
220    u = t / a
221    nm1 = n - 1
222    DO j = 1, nm1
223      v1 = u1 + j * u
224      v2 = (n-j) * u2
225      c = s1 * v1 - s2 * v2
226      d = s1 * v2 + s2 * v1
227      s1 = c
228      s2 = d
229    END DO
230  END IF
231 
232!     SET  W1 + W2*I = LOG(S1 + S2*I)  WHEN MO IS NONZERO
233 
234  s = s1 * s1 + s2 * s2
235  IF (mo /= 0) THEN
236    w1 = 0.5_dp * LOG(s)
237    w2 = ATAN2(s2,s1)
238  END IF
239END IF
240 
241!     SET  V1 + V2*I = (Z - 0.5) * LOG(Z + N) - Z
242 
243t1 = 0.5_dp * LOG(a) - 1._dp
244t2 = ATAN2(y,t)
245u = x - 0.5_dp
246v1 = (u*t1-0.5_dp) - y * t2
247v2 = u * t2 + y * t1
248 
249!     LET A1 + A2*I BE THE ASYMPTOTIC SUM
250 
251u1 = t / a
252u2 = -y / a
253q1 = u1 * u1 - u2 * u2
254q2 = 2._dp * u1 * u2
255a1 = 0._dp
256a2 = 0._dp
257DO j = 1, 30
258  t1 = a1
259  t2 = a2
260  a1 = a1 + c0(j) * u1
261  a2 = a2 + c0(j) * u2
262  IF (a1 == t1) THEN
263    IF (a2 == t2) GO TO 40
264  END IF
265  t1 = u1 * q1 - u2 * q2
266  t2 = u1 * q2 + u2 * q1
267  u1 = t1
268  u2 = t2
269END DO
270!-----------------------------------------------------------------------
271!                 GATHERING TOGETHER THE RESULTS
272!-----------------------------------------------------------------------
27340 w1 = (((a1+hl2p)-w1)+v1) - n
274w2 = (a2-w2) + v2
275IF (REAL(z, KIND=dp) < 0.0_dp) GO TO 60
276IF (mo == 0) THEN
277 
278!     CASE WHEN THE REAL PART OF Z IS NONNEGATIVE AND MO = 0
279 
280  a = EXP(w1)
281  w1 = a * COS(w2)
282  w2 = a * SIN(w2)
283  IF (n == 0) GO TO 70
284  c = (s1*w1+s2*w2) / s
285  d = (s1*w2-s2*w1) / s
286  w1 = c
287  w2 = d
288  GO TO 70
289END IF
290 
291!     CASE WHEN THE REAL PART OF Z IS NONNEGATIVE AND MO IS NONZERO.
292!     THE ANGLE W2 IS REDUCED TO THE INTERVAL -PI < W2 <= PI.
293 
29450 IF (w2 <= pi) THEN
295  k = 0.5_dp - w2 / pi2
296  w2 = w2 + pi2 * k
297  GO TO 70
298END IF
299k = w2 / pi2 - 0.5_dp
300u = k + 1
301w2 = w2 - pi2 * u
302IF (w2 <= -pi) w2 = pi
303GO TO 70
304 
305!     CASE WHEN THE REAL PART OF Z IS NEGATIVE AND MO IS NONZERO
306 
30760 IF (mo /= 0) THEN
308  w1 = h1 - w1
309  w2 = h2 - w2
310  GO TO 50
311END IF
312 
313!     CASE WHEN THE REAL PART OF Z IS NEGATIVE AND MO = 0
314 
315a = EXP(-w1)
316t1 = a * COS(-w2)
317t2 = a * SIN(-w2)
318w1 = h1 * t1 - h2 * t2
319w2 = h1 * t2 + h2 * t1
320IF (n /= 0) THEN
321  c = w1 * s1 - w2 * s2
322  d = w1 * s2 + w2 * s1
323  w1 = c
324  w2 = d
325END IF
326 
327!     TERMINATION
328 
32970 w = CMPLX(w1, w2, KIND=dp)
330RETURN
331!-----------------------------------------------------------------------
332!             THE REQUESTED VALUE CANNOT BE COMPUTED
333!-----------------------------------------------------------------------
33480 w = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
335RETURN
336END SUBROUTINE dcgama

风平老涡 · 发表于 2021-2-9 00:01:57

[Fortran] 纯文本查看 复制代码

001FUNCTION dsin1(x) RESULT(fn_val)
002!-----------------------------------------------------------------------
003 
004!                REAL (dp) EVALUATION OF SIN(PI*X)
005 
006!                             --------------
007 
008!     THE EXPANSION FOR SIN(PI*A) (ABS(A) <= PI/4) USING A1,...,A13
009!     IS ACCURATE TO WITHIN 2 UNITS OF THE 40-TH SIGNIFICANT DIGIT, AND
010!     THE EXPANSION FOR COS(PI*A) (ABS(A) <= PI/4) USING B1,...,B13
011!     IS ACCURATE TO WITHIN 4 UNITS OF THE 40-TH SIGNIFICANT DIGIT.
012 
013!-----------------------------------------------------------------------
014REAL (dp), INTENT(IN)  :: x
015REAL (dp)              :: fn_val
016 
017! Local variables
018REAL (dp)            :: a, t, w
019REAL (dp), PARAMETER :: pi = 3.141592653589793238462643383279502884197D+00
020REAL (dp), PARAMETER :: a1 = -.1028083791780141522795259479153765743002D+00,   &
021      a2  = .3170868848763100170457042079710451905600D-02,   &
022      a3  = -.4657026956105571623449026167864697920000D-04,  &
023      a4  = .3989844942879455643410226655783424000000D-06,   &
024      a5  = -.2237397227721999776371894030796800000000D-08,  &
025      a6  = .8847045483056962709715066675200000000000D-11,   &
026      a7  = -.2598715447506450292885585920000000000000D-13,  &
027      a8  = .5893449774331011070033920000000000000000D-16 ,  &
028      a9  = -.1062975472045522550784000000000000000000D-18,   &
029      a10 = .1561182648301780992000000000000000000000D-21,    &
030      a11 = -.1903193516670976000000000000000000000000D-24,   &
031      a12 = .1956617650176000000000000000000000000000D-27,    &
032      a13 = -.1711276032000000000000000000000000000000D-30
033REAL (dp), PARAMETER :: b1 = -.3084251375340424568385778437461297229882D+00, &
034      b2  = .1585434424381550085228521039855226435920D-01,   &
035      b3  = -.3259918869273900136414318317506279360000D-03,  &
036      b4  = .3590860448591510079069203991239232000000D-05,   &
037      b5  = -.2461136950494199754009084061808640000000D-07,  &
038      b6  = .1150115912797405152263195572224000000000D-09,   &
039      b7  = -.3898073171259675439899172864000000000000D-12,  &
040      b8  = .1001886461636271969091584000000000000000D-14,   &
041      b9  = -.2019653396886572027084800000000000000000D-17,  &
042      b10 = .3278483561466560512000000000000000000000D-20,   &
043      b11 = -.4377345082051788800000000000000000000000D-23,  &
044      b12 = .4891532381388800000000000000000000000000D-26,   &
045      b13 = -.4617089843200000000000000000000000000000D-29
046INTEGER  :: max, n
047!------------------------
048 
049!     ****** MAX IS A MACHINE DEPENDENT CONSTANT. MAX IS THE
050!            LARGEST POSITIVE INTEGER THAT MAY BE USED.
051 
052!                       MAX = IPMPAR(3)
053max = HUGE(3)
054 
055!------------------------
056a = ABS(x)
057t = MAX
058IF (a >= t) THEN
059  fn_val = 0.0_dp
060  RETURN
061END IF
062 
063n = a
064t = n
065a = a - t
066IF (a <= 0.75_dp) THEN
067  IF (a < 0.25_dp) GO TO 10
068 
069!                    0.25 <= A <= 0.75
070 
071  a = 0.25_dp + (0.25_dp-a)
072  t = 16._dp * a * a
073  fn_val = (((((((((((((b13*t + b12)*t + b11)*t + b10)*t + b9)*t + b8)*t  &
074           + b7)*t + b6)*t + b5)*t + b4)*t + b3)*t + b2)*t + b1)*t +  &
075           0.5_dp) + 0.5_dp
076  GO TO 20
077END IF
078 
079!                 A < 0.25  OR  A > 0.75
080 
081a = 0.25_dp + (0.75_dp-a)
08210 t = 16._dp * a * a
083w = (((((((((((((a13*t + a12)*t + a11)*t + a10)*t + a9)*t + a8)*t + a7)*t  &
084    + a6)*t + a5)*t + a4)*t + a3)*t + a2)*t + a1)*t + 0.5_dp) + 0.5_dp
085fn_val = pi * a * w
086 
087!                        TERMINATION
088 
08920 IF (x < 0.0) fn_val = -fn_val
090IF (MOD(n,2) /= 0) fn_val = -fn_val
091RETURN
092END FUNCTION dsin1
093 
094 
095 
096FUNCTION dcos1 (x) RESULT(fn_val)
097  
098!-----------------------------------------------------------------------
099 
100!                REAL (dp) EVALUATION OF COS(PI*X)
101 
102!                             --------------
103 
104!     THE EXPANSION FOR SIN(PI*A) (ABS(A) .LE. PI/4) USING A1,...,A13
105!     IS ACCURATE TO WITHIN 2 UNITS OF THE 40-TH SIGNIFICANT DIGIT, AND
106!     THE EXPANSION FOR COS(PI*A) (ABS(A) .LE. PI/4) USING B1,...,B13
107!     IS ACCURATE TO WITHIN 4 UNITS OF THE 40-TH SIGNIFICANT DIGIT.
108 
109!-----------------------------------------------------------------------
110 
111REAL (dp), INTENT(IN)  :: x
112REAL (dp)              :: fn_val
113 
114REAL (dp)  :: a, t, w
115INTEGER    :: MAX, n
116!------------------------
117REAL (dp), PARAMETER  :: pi = 3.141592653589793238462643383279502884197_dp
118!------------------------
119REAL (dp), PARAMETER  :: &
120    a1  = -.1028083791780141522795259479153765743002D+00,  &
121    a2  =  .3170868848763100170457042079710451905600D-02,  &
122    a3  = -.4657026956105571623449026167864697920000D-04,  &
123    a4  =  .3989844942879455643410226655783424000000D-06,  &
124    a5  = -.2237397227721999776371894030796800000000D-08,  &
125    a6  =  .8847045483056962709715066675200000000000D-11,  &
126    a7  = -.2598715447506450292885585920000000000000D-13,  &
127    a8  =  .5893449774331011070033920000000000000000D-16,  &
128    a9  = -.1062975472045522550784000000000000000000D-18,  &
129    a10 =  .1561182648301780992000000000000000000000D-21,  &
130    a11 = -.1903193516670976000000000000000000000000D-24,  &
131    a12 =  .1956617650176000000000000000000000000000D-27,  &
132    a13 = -.1711276032000000000000000000000000000000D-30
133!------------------------
134REAL (dp), PARAMETER  :: &
135    b1  = -.3084251375340424568385778437461297229882D+00,  &
136    b2  =  .1585434424381550085228521039855226435920D-01,  &
137    b3  = -.3259918869273900136414318317506279360000D-03,  &
138    b4  =  .3590860448591510079069203991239232000000D-05,  &
139    b5  = -.2461136950494199754009084061808640000000D-07,  &
140    b6  =  .1150115912797405152263195572224000000000D-09,  &
141    b7  = -.3898073171259675439899172864000000000000D-12,  &
142    b8  =  .1001886461636271969091584000000000000000D-14,  &
143    b9  = -.2019653396886572027084800000000000000000D-17,  &
144    b10 =  .3278483561466560512000000000000000000000D-20,  &
145    b11 = -.4377345082051788800000000000000000000000D-23,  &
146    b12 =  .4891532381388800000000000000000000000000D-26,  &
147    b13 = -.4617089843200000000000000000000000000000D-29
148!------------------------
149 
150!     ****** MAX IS A MACHINE DEPENDENT CONSTANT. MAX IS THE
151!            LARGEST POSITIVE INTEGER THAT MAY BE USED.
152 
153MAX = HUGE(0)
154 
155!------------------------
156a = ABS(x)
157t = MAX
158IF (a < t) GO TO 10
159fn_val = 1.d0
160RETURN
161 
16210 n = a
163t = n
164a = a - t
165IF (a > 0.75D0) GO TO 20
166IF (a < 0.25D0) GO TO 21
167 
168!                    0.25 .LE. A .LE. 0.75
169 
170a = 0.25D0 + (0.25D0 - a)
171t = 16.d0*a*a
172w = (((((((((((((a13*t + a12)*t + a11)*t + a10)*t + a9)*t +  &
173    a8)*t + a7)*t + a6)*t + a5)*t + a4)*t + a3)*t + a2)*t + a1)*t + 0.5D0) + 0.5D0
174fn_val = pi*a*w
175GO TO 30
176 
177!                 A .LT. 0.25  OR  A .GT. 0.75
178 
17920 a = 0.25D0 + (0.75D0 - a)
180n = n - 1
18121 t = 16.d0*a*a
182fn_val = (((((((((((((b13*t + b12)*t + b11)*t + b10)*t + b9)*t + b8)*t + &
183         b7)*t + b6)*t + b5)*t + b4)*t + b3)*t + b2)*t + b1)*t + 0.5D0) + 0.5D0
184 
185!                        TERMINATION
186 
18730 IF (MOD(n,2) /= 0) fn_val = -fn_val
188RETURN
189END FUNCTION dcos1

[通用算法] 求助-复变量的贝塞尔第一类和第二类函数的Fortran子程序

规矩勋章