1C Output from Public domain Ratfor, version 1.01 2 subroutine qpsedg8xf(tgiyxdw1, dufozmt7, wy1vqfzu) 3 implicit logical (a-z) 4 integer wy1vqfzu, tgiyxdw1(*), dufozmt7(*) 5 integer urohxe6t, bpvaqm5z, ayfnwr1v 6 ayfnwr1v = 1 7 urohxe6t = wy1vqfzu 823000 if(.not.(urohxe6t .ge. 1))goto 23002 9 do23003 bpvaqm5z=1,urohxe6t 10 tgiyxdw1(ayfnwr1v) = bpvaqm5z 11 ayfnwr1v = ayfnwr1v+1 1223003 continue 1323004 continue 1423001 urohxe6t=urohxe6t-1 15 goto 23000 1623002 continue 17 ayfnwr1v = 1 18 do23005 urohxe6t=1,wy1vqfzu 19 do23007 bpvaqm5z=urohxe6t,wy1vqfzu 20 dufozmt7(ayfnwr1v) = bpvaqm5z 21 ayfnwr1v = ayfnwr1v+1 2223007 continue 2323008 continue 2423005 continue 2523006 continue 26 return 27 end 28 integer function viamf(cz8qdfyj, rvy1fpli, wy1vqfzu, tgiyxdw1, duf 29 *ozmt7) 30 integer cz8qdfyj, rvy1fpli, wy1vqfzu, tgiyxdw1(*), dufozmt7(*) 31 integer urohxe6t, imk5wjxg 32 imk5wjxg = wy1vqfzu*(wy1vqfzu+1)/2 33 do23009 urohxe6t=1,imk5wjxg 34 if((tgiyxdw1(urohxe6t).eq.cz8qdfyj .and. dufozmt7(urohxe6t).eq.rvy 35 *1fpli) .or. (tgiyxdw1(urohxe6t).eq.rvy1fpli .and. dufozmt7(urohxe6 36 *t).eq.cz8qdfyj))then 37 viamf = urohxe6t 38 return 39 endif 4023009 continue 4123010 continue 42 viamf = 0 43 return 44 end 45 subroutine vm2af(mat, a, dimm, tgiyxdw1, dufozmt7, kuzxj1lo, wy1vq 46 *fzu, rb1onzwu) 47 implicit logical (a-z) 48 integer dimm, tgiyxdw1(dimm), dufozmt7(dimm), kuzxj1lo, wy1vqfzu, 49 *rb1onzwu 50 double precision mat(dimm,kuzxj1lo), a(wy1vqfzu,wy1vqfzu,kuzxj1lo) 51 integer ayfnwr1v, yq6lorbx, gp1jxzuh, imk5wjxg 52 imk5wjxg = wy1vqfzu * (wy1vqfzu + 1) / 2 53 if(rb1onzwu .eq. 1 .or. dimm .ne. imk5wjxg)then 54 ayfnwr1v = 1 5523015 if(.not.(ayfnwr1v .le. kuzxj1lo))goto 23017 56 yq6lorbx = 1 5723018 if(.not.(yq6lorbx .le. wy1vqfzu))goto 23020 58 gp1jxzuh = 1 5923021 if(.not.(gp1jxzuh .le. wy1vqfzu))goto 23023 60 a(gp1jxzuh,yq6lorbx,ayfnwr1v) = 0.0d0 6123022 gp1jxzuh=gp1jxzuh+1 62 goto 23021 6323023 continue 6423019 yq6lorbx=yq6lorbx+1 65 goto 23018 6623020 continue 6723016 ayfnwr1v=ayfnwr1v+1 68 goto 23015 6923017 continue 70 endif 71 do23024 ayfnwr1v=1,kuzxj1lo 72 do23026 yq6lorbx=1,dimm 73 a(tgiyxdw1(yq6lorbx),dufozmt7(yq6lorbx),ayfnwr1v) = mat(yq6lorbx,a 74 *yfnwr1v) 75 if(rb1onzwu .eq. 0)then 76 a(dufozmt7(yq6lorbx),tgiyxdw1(yq6lorbx),ayfnwr1v) = mat(yq6lorbx,a 77 *yfnwr1v) 78 endif 7923026 continue 8023027 continue 8123024 continue 8223025 continue 83 return 84 end 85 subroutine mux22f(wpuarq2m, tlgduey8, lfu2qhid, dimu, tgiyxdw1, du 86 *fozmt7, kuzxj1lo, wy1vqfzu, wk1200) 87 implicit logical (a-z) 88 integer dimu, tgiyxdw1(*), dufozmt7(*), kuzxj1lo, wy1vqfzu 89 double precision wpuarq2m(dimu,kuzxj1lo), tlgduey8(kuzxj1lo,wy1vqf 90 *zu), lfu2qhid(wy1vqfzu,kuzxj1lo), wk1200(wy1vqfzu,wy1vqfzu) 91 double precision q6zdcwxk 92 integer ayfnwr1v, yq6lorbx, bpvaqm5z, one, rb1onzwu 93 one = 1 94 rb1onzwu = 1 95 ayfnwr1v = 1 9623030 if(.not.(ayfnwr1v .le. kuzxj1lo))goto 23032 97 call vm2af(wpuarq2m(1,ayfnwr1v), wk1200, dimu, tgiyxdw1, dufozmt7, 98 * one, wy1vqfzu, rb1onzwu) 99 yq6lorbx = 1 10023033 if(.not.(yq6lorbx .le. wy1vqfzu))goto 23035 101 q6zdcwxk = 0.0d0 102 bpvaqm5z = yq6lorbx 10323036 if(.not.(bpvaqm5z .le. wy1vqfzu))goto 23038 104 q6zdcwxk = q6zdcwxk + wk1200(yq6lorbx,bpvaqm5z) * tlgduey8(ayfnwr1 105 *v,bpvaqm5z) 10623037 bpvaqm5z=bpvaqm5z+1 107 goto 23036 10823038 continue 109 lfu2qhid(yq6lorbx,ayfnwr1v) = q6zdcwxk 11023034 yq6lorbx=yq6lorbx+1 111 goto 23033 11223035 continue 11323031 ayfnwr1v=ayfnwr1v+1 114 goto 23030 11523032 continue 116 return 117 end 118 subroutine vbksf(wpuarq2m, bvecto, wy1vqfzu, kuzxj1lo, wk1200, tgi 119 *yxdw1, dufozmt7, dimu) 120 implicit logical (a-z) 121 integer wy1vqfzu, kuzxj1lo, tgiyxdw1(*), dufozmt7(*), dimu 122 double precision wpuarq2m(dimu,kuzxj1lo), bvecto(wy1vqfzu,kuzxj1lo 123 *), wk1200(wy1vqfzu,wy1vqfzu) 124 double precision q6zdcwxk 125 integer ayfnwr1v, yq6lorbx, gp1jxzuh, rb1onzwu, one 126 rb1onzwu = 1 127 one = 1 128 ayfnwr1v = 1 12923039 if(.not.(ayfnwr1v .le. kuzxj1lo))goto 23041 130 call vm2af(wpuarq2m(1,ayfnwr1v), wk1200, dimu, tgiyxdw1, dufozmt7, 131 * one, wy1vqfzu, rb1onzwu) 132 yq6lorbx = wy1vqfzu 13323042 if(.not.(yq6lorbx .ge. 1))goto 23044 134 q6zdcwxk = bvecto(yq6lorbx,ayfnwr1v) 135 gp1jxzuh = yq6lorbx+1 13623045 if(.not.(gp1jxzuh .le. wy1vqfzu))goto 23047 137 q6zdcwxk = q6zdcwxk - wk1200(yq6lorbx,gp1jxzuh) * bvecto(gp1jxzuh, 138 *ayfnwr1v) 13923046 gp1jxzuh=gp1jxzuh+1 140 goto 23045 14123047 continue 142 bvecto(yq6lorbx,ayfnwr1v) = q6zdcwxk / wk1200(yq6lorbx,yq6lorbx) 14323043 yq6lorbx=yq6lorbx-1 144 goto 23042 14523044 continue 14623040 ayfnwr1v=ayfnwr1v+1 147 goto 23039 14823041 continue 149 return 150 end 151 subroutine vcholf(wmat, bvecto, wy1vqfzu, dvhw1ulq, isolve) 152 implicit logical (a-z) 153 integer isolve 154 integer wy1vqfzu, dvhw1ulq 155 double precision wmat(wy1vqfzu,wy1vqfzu), bvecto(wy1vqfzu) 156 double precision q6zdcwxk, dsqrt 157 integer ayfnwr1v, yq6lorbx, gp1jxzuh 158 dvhw1ulq=1 159 do23048 ayfnwr1v=1,wy1vqfzu 160 q6zdcwxk = 0d0 161 do23050 gp1jxzuh=1,ayfnwr1v-1 162 q6zdcwxk = q6zdcwxk + wmat(gp1jxzuh,ayfnwr1v) * wmat(gp1jxzuh,ayfn 163 *wr1v) 16423050 continue 16523051 continue 166 wmat(ayfnwr1v,ayfnwr1v) = wmat(ayfnwr1v,ayfnwr1v) - q6zdcwxk 167 if(wmat(ayfnwr1v,ayfnwr1v) .le. 0d0)then 168 dvhw1ulq = 0 169 return 170 endif 171 wmat(ayfnwr1v,ayfnwr1v) = dsqrt(wmat(ayfnwr1v,ayfnwr1v)) 172 do23054 yq6lorbx=ayfnwr1v+1,wy1vqfzu 173 q6zdcwxk = 0d0 174 do23056 gp1jxzuh=1,ayfnwr1v-1 175 q6zdcwxk = q6zdcwxk + wmat(gp1jxzuh,ayfnwr1v) * wmat(gp1jxzuh,yq6l 176 *orbx) 17723056 continue 17823057 continue 179 wmat(ayfnwr1v,yq6lorbx) = (wmat(ayfnwr1v,yq6lorbx) - q6zdcwxk) / w 180 *mat(ayfnwr1v,ayfnwr1v) 18123054 continue 18223055 continue 18323048 continue 18423049 continue 185 if(isolve .eq. 0)then 186 do23060 ayfnwr1v=2,wy1vqfzu 187 do23062 yq6lorbx=1,ayfnwr1v-1 188 wmat(ayfnwr1v,yq6lorbx) = 0.0d0 18923062 continue 19023063 continue 191 return 19223060 continue 19323061 continue 194 endif 195 do23064 yq6lorbx=1,wy1vqfzu 196 q6zdcwxk = bvecto(yq6lorbx) 197 do23066 gp1jxzuh=1,yq6lorbx-1 198 q6zdcwxk = q6zdcwxk - wmat(gp1jxzuh,yq6lorbx) * bvecto(gp1jxzuh) 19923066 continue 20023067 continue 201 bvecto(yq6lorbx) = q6zdcwxk / wmat(yq6lorbx,yq6lorbx) 20223064 continue 20323065 continue 204 yq6lorbx = wy1vqfzu 20523068 if(.not.(yq6lorbx .ge. 1))goto 23070 206 q6zdcwxk = bvecto(yq6lorbx) 207 gp1jxzuh = yq6lorbx+1 20823071 if(.not.(gp1jxzuh .le. wy1vqfzu))goto 23073 209 q6zdcwxk = q6zdcwxk - wmat(yq6lorbx,gp1jxzuh) * bvecto(gp1jxzuh) 21023072 gp1jxzuh=gp1jxzuh+1 211 goto 23071 21223073 continue 213 bvecto(yq6lorbx) = q6zdcwxk / wmat(yq6lorbx,yq6lorbx) 21423069 yq6lorbx=yq6lorbx-1 215 goto 23068 21623070 continue 217 return 218 end 219 subroutine mux17f(wpuarq2m, he7mqnvy, wy1vqfzu, xjc4ywlh, kuzxj1lo 220 *, wk1200, wk3400, tgiyxdw1, dufozmt7, dimu, rutyk8mg) 221 implicit logical (a-z) 222 integer dimu, wy1vqfzu, xjc4ywlh, kuzxj1lo, tgiyxdw1(*), dufozmt7( 223 **), rutyk8mg 224 double precision wpuarq2m(dimu,kuzxj1lo), he7mqnvy(rutyk8mg,xjc4yw 225 *lh), wk1200(wy1vqfzu,wy1vqfzu), wk3400(wy1vqfzu,xjc4ywlh) 226 double precision q6zdcwxk 227 integer ayfnwr1v, yq6lorbx, gp1jxzuh, bpvaqm5z 228 do23074 yq6lorbx=1,wy1vqfzu 229 do23076 ayfnwr1v=1,wy1vqfzu 230 wk1200(ayfnwr1v,yq6lorbx) = 0.0d0 23123076 continue 23223077 continue 23323074 continue 23423075 continue 235 do23078 ayfnwr1v=1,kuzxj1lo 236 do23080 bpvaqm5z=1,dimu 237 wk1200(tgiyxdw1(bpvaqm5z), dufozmt7(bpvaqm5z)) = wpuarq2m(bpvaqm5z 238 *,ayfnwr1v) 23923080 continue 24023081 continue 241 do23082 gp1jxzuh=1,xjc4ywlh 242 do23084 yq6lorbx=1,wy1vqfzu 243 wk3400(yq6lorbx,gp1jxzuh) = he7mqnvy((ayfnwr1v-1)*wy1vqfzu+yq6lorb 244 *x,gp1jxzuh) 24523084 continue 24623085 continue 24723082 continue 24823083 continue 249 do23086 gp1jxzuh=1,xjc4ywlh 250 do23088 yq6lorbx=1,wy1vqfzu 251 q6zdcwxk = 0d0 252 do23090 bpvaqm5z=yq6lorbx,wy1vqfzu 253 q6zdcwxk = q6zdcwxk + wk1200(yq6lorbx,bpvaqm5z) * wk3400(bpvaqm5z, 254 *gp1jxzuh) 25523090 continue 25623091 continue 257 he7mqnvy((ayfnwr1v-1)*wy1vqfzu+yq6lorbx,gp1jxzuh) = q6zdcwxk 25823088 continue 25923089 continue 26023086 continue 26123087 continue 26223078 continue 26323079 continue 264 return 265 end 266 subroutine vrinvf9(wpuarq2m, ldr, wy1vqfzu, dvhw1ulq, ks3wejcv, wo 267 *rk) 268 implicit logical (a-z) 269 integer ldr, wy1vqfzu, dvhw1ulq 270 double precision wpuarq2m(ldr,wy1vqfzu), ks3wejcv(wy1vqfzu,wy1vqfz 271 *u), work(wy1vqfzu,wy1vqfzu) 272 double precision q6zdcwxk 273 integer yq6lorbx, gp1jxzuh, col, uaoynef0 274 dvhw1ulq = 1 275 yq6lorbx = 1 27623092 if(.not.(yq6lorbx .le. wy1vqfzu))goto 23094 277 col = 1 27823095 if(.not.(col .le. wy1vqfzu))goto 23097 279 work(yq6lorbx,col) = 0.0d0 28023096 col=col+1 281 goto 23095 28223097 continue 28323093 yq6lorbx=yq6lorbx+1 284 goto 23092 28523094 continue 286 col = 1 28723098 if(.not.(col .le. wy1vqfzu))goto 23100 288 yq6lorbx = col 28923101 if(.not.(yq6lorbx .ge. 1))goto 23103 290 if(yq6lorbx .eq. col)then 291 q6zdcwxk = 1.0d0 292 else 293 q6zdcwxk = 0.0d0 294 endif 295 gp1jxzuh = yq6lorbx+1 29623106 if(.not.(gp1jxzuh .le. col))goto 23108 297 q6zdcwxk = q6zdcwxk - wpuarq2m(yq6lorbx,gp1jxzuh) * work(gp1jxzuh, 298 *col) 29923107 gp1jxzuh=gp1jxzuh+1 300 goto 23106 30123108 continue 302 if(wpuarq2m(yq6lorbx,yq6lorbx) .eq. 0.0d0)then 303 dvhw1ulq = 0 304 else 305 work(yq6lorbx,col) = q6zdcwxk / wpuarq2m(yq6lorbx,yq6lorbx) 306 endif 30723102 yq6lorbx=yq6lorbx-1 308 goto 23101 30923103 continue 31023099 col=col+1 311 goto 23098 31223100 continue 313 yq6lorbx = 1 31423111 if(.not.(yq6lorbx .le. wy1vqfzu))goto 23113 315 col = yq6lorbx 31623114 if(.not.(col .le. wy1vqfzu))goto 23116 317 if(yq6lorbx .lt. col)then 318 uaoynef0 = col 319 else 320 uaoynef0 = yq6lorbx 321 endif 322 q6zdcwxk = 0.0d0 323 gp1jxzuh = uaoynef0 32423119 if(.not.(gp1jxzuh .le. wy1vqfzu))goto 23121 325 q6zdcwxk = q6zdcwxk + work(yq6lorbx,gp1jxzuh) * work(col,gp1jxzuh) 32623120 gp1jxzuh=gp1jxzuh+1 327 goto 23119 32823121 continue 329 ks3wejcv(yq6lorbx,col) = q6zdcwxk 330 ks3wejcv(col,yq6lorbx) = q6zdcwxk 33123115 col=col+1 332 goto 23114 33323116 continue 33423112 yq6lorbx=yq6lorbx+1 335 goto 23111 33623113 continue 337 return 338 end 339 subroutine tldz5ion(xx, lfu2qhid) 340 implicit logical (a-z) 341 double precision xx, lfu2qhid 342 double precision x, y, hofjnx2e, q6zdcwxk, xd4mybgj(6) 343 integer yq6lorbx 344 xd4mybgj(1)= 76.18009172947146d0 345 xd4mybgj(2)= -86.50532032941677d0 346 xd4mybgj(3)= 24.01409824083091d0 347 xd4mybgj(4)= -1.231739572450155d0 348 xd4mybgj(5)= 0.1208650973866179d-2 349 xd4mybgj(6)= -0.5395239384953d-5 350 x = xx 351 y = xx 352 hofjnx2e = x+5.50d0 353 hofjnx2e = hofjnx2e - (x+0.50d0) * dlog(hofjnx2e) 354 q6zdcwxk=1.000000000190015d0 355 yq6lorbx=1 35623122 if(.not.(yq6lorbx .le. 6))goto 23124 357 y = y + 1.0d0 358 q6zdcwxk = q6zdcwxk + xd4mybgj(yq6lorbx)/y 35923123 yq6lorbx=yq6lorbx+1 360 goto 23122 36123124 continue 362 lfu2qhid = -hofjnx2e + dlog(2.5066282746310005d0 * q6zdcwxk / x) 363 return 364 end 365 subroutine enbin9(bzmd6ftv, hdqsx7bk, nm0eljqk, n2kersmx, n, dvhw1 366 *ulq, zy1mchbf, ux3nadiw, rsynp1go, sguwj9ty) 367 implicit logical (a-z) 368 integer n, dvhw1ulq, zy1mchbf, sguwj9ty 369 double precision bzmd6ftv(n, zy1mchbf), hdqsx7bk(n, zy1mchbf), nm0 370 *eljqk(n, zy1mchbf), n2kersmx, ux3nadiw, rsynp1go 371 integer ayfnwr1v, kij0gwer 372 double precision oxjgzv0e, btiehdm2, ydb, vjz5sxty, esql7umk, pvcj 373 *l2na, mwuvskg1, ft3ijqmy, hmayv1xt, q6zdcwxk, plo6hkdr 374 real csi9ydge 375 if(n2kersmx .le. 0.80d0 .or. n2kersmx .ge. 1.0d0)then 376 dvhw1ulq = 0 377 return 378 endif 379 btiehdm2 = 100.0d0 * rsynp1go 380 oxjgzv0e = 0.001d0 381 dvhw1ulq = 1 382 kij0gwer=1 38323127 if(.not.(kij0gwer.le.zy1mchbf))goto 23129 384 ayfnwr1v=1 38523130 if(.not.(ayfnwr1v.le.n))goto 23132 386 vjz5sxty = nm0eljqk(ayfnwr1v,kij0gwer) / hdqsx7bk(ayfnwr1v,kij0gwe 387 *r) 388 if((vjz5sxty .lt. oxjgzv0e) .or. (nm0eljqk(ayfnwr1v,kij0gwer) .gt. 389 * 1.0d5))then 390 bzmd6ftv(ayfnwr1v,kij0gwer) = -nm0eljqk(ayfnwr1v,kij0gwer) * (1.0d 391 *0 + hdqsx7bk(ayfnwr1v,kij0gwer)/(hdqsx7bk(ayfnwr1v,kij0gwer) + nm0 392 *eljqk(ayfnwr1v,kij0gwer))) / hdqsx7bk(ayfnwr1v,kij0gwer)**2 393 if(bzmd6ftv(ayfnwr1v,kij0gwer) .gt. -btiehdm2)then 394 bzmd6ftv(ayfnwr1v,kij0gwer) = -btiehdm2 395 endif 396 goto 20 397 endif 398 q6zdcwxk = 0.0d0 399 pvcjl2na = hdqsx7bk(ayfnwr1v,kij0gwer) / (hdqsx7bk(ayfnwr1v,kij0gw 400 *er) + nm0eljqk(ayfnwr1v,kij0gwer)) 401 mwuvskg1 = 1.0d0 - pvcjl2na 402 csi9ydge = hdqsx7bk(ayfnwr1v,kij0gwer) 403 if(pvcjl2na .lt. btiehdm2)then 404 pvcjl2na = btiehdm2 405 endif 406 if(mwuvskg1 .lt. btiehdm2)then 407 mwuvskg1 = btiehdm2 408 endif 409 esql7umk = 100.0d0 + 15.0d0 * nm0eljqk(ayfnwr1v,kij0gwer) 410 if(esql7umk .lt. sguwj9ty)then 411 esql7umk = sguwj9ty 412 endif 413 ft3ijqmy = pvcjl2na ** csi9ydge 414 ux3nadiw = ft3ijqmy 415 plo6hkdr = (1.0d0 - ux3nadiw) / hdqsx7bk(ayfnwr1v,kij0gwer)**2 416 q6zdcwxk = q6zdcwxk + plo6hkdr 417 ydb = 1.0d0 418 ft3ijqmy = hdqsx7bk(ayfnwr1v,kij0gwer) * mwuvskg1 * ft3ijqmy 419 ux3nadiw = ux3nadiw + ft3ijqmy 420 plo6hkdr = (1.0d0 - ux3nadiw) / (hdqsx7bk(ayfnwr1v,kij0gwer) + ydb 421 *)**2 422 q6zdcwxk = q6zdcwxk + plo6hkdr 423 ydb = 2.0d0 42423143 if(((ux3nadiw .le. n2kersmx) .or. (plo6hkdr .gt. 1.0d-4)) .and. (y 425 *db .lt. esql7umk))then 426 ft3ijqmy = (hdqsx7bk(ayfnwr1v,kij0gwer) - 1.0d0 + ydb) * mwuvskg1 427 ** ft3ijqmy / ydb 428 ux3nadiw = ux3nadiw + ft3ijqmy 429 plo6hkdr = (1.0d0 - ux3nadiw) / (hdqsx7bk(ayfnwr1v,kij0gwer) + ydb 430 *)**2 431 q6zdcwxk = q6zdcwxk + plo6hkdr 432 ydb = ydb + 1.0d0 433 goto 23143 434 endif 43523144 continue 436 bzmd6ftv(ayfnwr1v,kij0gwer) = -q6zdcwxk 43720 hmayv1xt = 0.0d0 43823131 ayfnwr1v=ayfnwr1v+1 439 goto 23130 44023132 continue 44123128 kij0gwer=kij0gwer+1 442 goto 23127 44323129 continue 444 return 445 end 446 subroutine enbin8(bzmd6ftv, hdqsx7bk, hsj9bzaq, n2kersmx, kuzxj1lo 447 *, dvhw1ulq, zy1mchbf, ux3nadiw, rsynp1go) 448 implicit logical (a-z) 449 integer kuzxj1lo, dvhw1ulq, zy1mchbf 450 double precision bzmd6ftv(kuzxj1lo, zy1mchbf), hdqsx7bk(kuzxj1lo, 451 *zy1mchbf), hsj9bzaq(kuzxj1lo, zy1mchbf), n2kersmx, ux3nadiw, rsynp 452 *1go 453 integer ayfnwr1v, kij0gwer, esql7umk 454 double precision ft3ijqmy, tad5vhsu, o3jyipdf, pq0hfucn, q6zdcwxk, 455 * d1, d2, plo6hkdr, hnu1vjyw 456 logical pok1, pok2, pok12 457 double precision oxjgzv0e, onemse, nm0eljqk, btiehdm2, ydb, kbig 458 d1 = 0.0d0 459 d2 = 0.0d0 460 btiehdm2 = -100.0d0 * rsynp1go 461 esql7umk = 3000 462 if(n2kersmx .le. 0.80d0 .or. n2kersmx .ge. 1.0d0)then 463 dvhw1ulq = 0 464 return 465 endif 466 kbig = 1.0d4 467 oxjgzv0e = 0.001d0 468 hnu1vjyw = 1.0d0 - rsynp1go 469 onemse = 1.0d0 / (1.0d0 + oxjgzv0e) 470 dvhw1ulq = 1 471 kij0gwer=1 47223147 if(.not.(kij0gwer.le.zy1mchbf))goto 23149 473 ayfnwr1v=1 47423150 if(.not.(ayfnwr1v.le.kuzxj1lo))goto 23152 475 if(hdqsx7bk(ayfnwr1v,kij0gwer) .gt. kbig)then 476 hdqsx7bk(ayfnwr1v,kij0gwer) = kbig 477 endif 478 if(hsj9bzaq(ayfnwr1v,kij0gwer) .lt. oxjgzv0e)then 479 hsj9bzaq(ayfnwr1v,kij0gwer) = oxjgzv0e 480 endif 481 if((hsj9bzaq(ayfnwr1v,kij0gwer) .gt. onemse))then 482 nm0eljqk = hdqsx7bk(ayfnwr1v,kij0gwer) * (1.0d0/hsj9bzaq(ayfnwr1v, 483 *kij0gwer) - 1.0d0) 484 bzmd6ftv(ayfnwr1v,kij0gwer) = -nm0eljqk * (1.0d0 + hdqsx7bk(ayfnwr 485 *1v,kij0gwer)/(hdqsx7bk(ayfnwr1v,kij0gwer) + nm0eljqk)) / hdqsx7bk( 486 *ayfnwr1v,kij0gwer)**2 487 if(bzmd6ftv(ayfnwr1v,kij0gwer) .gt. btiehdm2)then 488 bzmd6ftv(ayfnwr1v,kij0gwer) = btiehdm2 489 endif 490 goto 20 491 endif 492 q6zdcwxk = 0.0d0 493 pok1 = .true. 494 pok2 = hsj9bzaq(ayfnwr1v,kij0gwer) .lt. (1.0d0-rsynp1go) 495 pok12 = pok1 .and. pok2 496 if(pok12)then 497 d2 = hdqsx7bk(ayfnwr1v,kij0gwer) * dlog(hsj9bzaq(ayfnwr1v,kij0gwer 498 *)) 499 ux3nadiw = dexp(d2) 500 else 501 ux3nadiw = 0.0d0 502 endif 503 plo6hkdr = (1.0d0 - ux3nadiw) / hdqsx7bk(ayfnwr1v,kij0gwer)**2 504 q6zdcwxk = q6zdcwxk + plo6hkdr 505 call tldz5ion(hdqsx7bk(ayfnwr1v,kij0gwer), o3jyipdf) 506 ydb = 1.0d0 507 call tldz5ion(ydb + hdqsx7bk(ayfnwr1v,kij0gwer), tad5vhsu) 508 pq0hfucn = 0.0d0 509 if(pok12)then 510 d1 = dlog(1.0d0 - hsj9bzaq(ayfnwr1v,kij0gwer)) 511 ft3ijqmy = dexp(ydb * d1 + d2 + tad5vhsu - o3jyipdf - pq0hfucn) 512 else 513 ft3ijqmy = 0.0d0 514 endif 515 ux3nadiw = ux3nadiw + ft3ijqmy 516 plo6hkdr = (1.0d0 - ux3nadiw) / (hdqsx7bk(ayfnwr1v,kij0gwer) + ydb 517 *)**2 518 q6zdcwxk = q6zdcwxk + plo6hkdr 519 ydb = 2.0d0 52023165 if((ux3nadiw .le. n2kersmx) .or. (plo6hkdr .gt. 1.0d-4))then 521 tad5vhsu = tad5vhsu + dlog(ydb + hdqsx7bk(ayfnwr1v,kij0gwer) - 1.0 522 *d0) 523 pq0hfucn = pq0hfucn + dlog(ydb) 524 if(pok12)then 525 ft3ijqmy = dexp(ydb * d1 + d2 + tad5vhsu - o3jyipdf - pq0hfucn) 526 else 527 ft3ijqmy = 0.0d0 528 endif 529 ux3nadiw = ux3nadiw + ft3ijqmy 530 plo6hkdr = (1.0d0 - ux3nadiw) / (hdqsx7bk(ayfnwr1v,kij0gwer) + ydb 531 *)**2 532 q6zdcwxk = q6zdcwxk + plo6hkdr 533 ydb = ydb + 1.0d0 534 if(ydb .gt. 1.0d3)then 535 goto 21 536 endif 537 goto 23165 538 endif 53923166 continue 54021 bzmd6ftv(ayfnwr1v,kij0gwer) = -q6zdcwxk 54120 tad5vhsu = 0.0d0 54223151 ayfnwr1v=ayfnwr1v+1 543 goto 23150 54423152 continue 54523148 kij0gwer=kij0gwer+1 546 goto 23147 54723149 continue 548 return 549 end 550 subroutine mbessi0(bvecto, kuzxj1lo, kpzavbj3, d0, d1, d2, zjkrtol 551 *8, qaltf0nz) 552 implicit logical (a-z) 553 integer kuzxj1lo, kpzavbj3, zjkrtol8, c5aesxkus 554 double precision bvecto(kuzxj1lo), d0(kuzxj1lo), d1(kuzxj1lo), d2( 555 *kuzxj1lo), qaltf0nz 556 integer ayfnwr1v, gp1jxzuh 557 double precision f0, t0, m0, f1, t1, m1, f2, t2, m2 558 double precision toobig 559 toobig = 20.0d0 560 zjkrtol8 = 0 561 if(.not.(kpzavbj3 .eq. 0 .or. kpzavbj3 .eq. 1 .or. kpzavbj3 .eq. 2 562 *))then 563 zjkrtol8 = 1 564 return 565 endif 566 do23173 gp1jxzuh=1,kuzxj1lo 567 if(dabs(bvecto(gp1jxzuh)) .gt. toobig)then 568 zjkrtol8 = 1 569 return 570 endif 571 t1 = bvecto(gp1jxzuh) / 2.0d0 572 f1 = t1 573 t0 = t1 * t1 574 f0 = 1.0d0 + t0 575 t2 = 0.50d0 576 f2 = t2 577 c5aesxkus = 15 578 if(dabs(bvecto(gp1jxzuh)) .gt. 10)then 579 c5aesxkus = 25 580 endif 581 if(dabs(bvecto(gp1jxzuh)) .gt. 15)then 582 c5aesxkus = 35 583 endif 584 if(dabs(bvecto(gp1jxzuh)) .gt. 20)then 585 c5aesxkus = 40 586 endif 587 if(dabs(bvecto(gp1jxzuh)) .gt. 30)then 588 c5aesxkus = 55 589 endif 590 do23185 ayfnwr1v=1,c5aesxkus 591 m0 = (bvecto(gp1jxzuh) / (2.0d0*(ayfnwr1v+1.0d0))) ** 2.0 592 m1 = m0 * (1.0d0 + 1.0d0/ayfnwr1v) 593 m2 = m1 * (2.0d0*ayfnwr1v + 1.0d0) / (2.0d0*ayfnwr1v - 1.0d0) 594 t0 = t0 * m0 595 t1 = t1 * m1 596 t2 = t2 * m2 597 f0 = f0 + t0 598 f1 = f1 + t1 599 f2 = f2 + t2 600 if((dabs(t0) .lt. qaltf0nz) .and. (dabs(t1) .lt. qaltf0nz) .and. ( 601 *dabs(t2) .lt. qaltf0nz))then 602 goto 23186 603 endif 60423185 continue 60523186 continue 606 if(0 .le. kpzavbj3)then 607 d0(gp1jxzuh) = f0 608 endif 609 if(1 .le. kpzavbj3)then 610 d1(gp1jxzuh) = f1 611 endif 612 if(2 .le. kpzavbj3)then 613 d2(gp1jxzuh) = f2 614 endif 61523173 continue 61623174 continue 617 return 618 end 619