API tools faq
paste
Guest User

exit issues

a guest
Nov 7th, 2025
19
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 31.19 KB | None | 0 0
raw download clone embed print report
  1. C Maybe we need a comment
  2. C
  3. PROGRAM STBRHY
  4. C
  5. C PROGRAM FOR CALCULATING THE STARKBROADENING OF HYDROGEN ON THE
  6. C BASIS OF THE UNIFIED THEORY INCLUDING LOWER STATE INTERACTION
  7. C
  8. C ALLOCATE ARRAYS AND MATRICIES
  9. DIMENSION DQQ(30),PFAC(6),STRONG(34),FF(1100),SJUU(136,31),
  10. & SJLL(6,5),RDM(16,3),SLUU(16),SLLL(3),NQUU(34,2),NQLL(34,2),
  11. & MMUU(34,2),MMLL(34,2),NQCU(30),NQCL(30),NXX(30),NST(30),
  12. & DST(30),W(300,5)
  13. C DEFINE VARIABLES WHICH CAN BE ACCESSED ANYWHERE IN THE PROGRAM
  14. COMMON /FDAT/ P1,P2,B1,A2,B2,PPFF,CCREAL,CCIMAG
  15. COMMON /PFW /FIELD(301)
  16. COMMON /11/ SLOP(595,22,2)
  17. COMMON /22/ DIPOL(34,2),NUR(2),IPI,NCQ,BET,FPAR(6,30),ADIFF(34,2),NSTEP
  18. COMMON /33/ AMATR(34,68), BMATR(34,34), CMATR(34,34)
  19. C DEFINE ALIASES FOR CONVENIENCE (Q: WHAT IS FF?)
  20. EQUIVALENCE (SJUU,AMATR),(FF(601),NQUU),(FF(669),NQLL),
  21. & (FF(737),MMUU),(FF(805),MMLL),(FF(873),NQCU),(FF(903),NQCL),
  22. & (FF(933),SJLL),(FF(963),RDM),(FF(1011),SLUU),(FF(1027),STRONG)
  23. FIELD(1) = 0.0
  24. C READ IN W AS A (300,5) MATRIX
  25. READ 100, ((W(I,J), J = 1,5), I = 1,300)
  26. 100 FORMAT (5E12.4)
  27. C
  28. C READ IN THE PHYSICAL PARAMETERS WE'RE INTERESTED IN
  29. 120 READ 150, DEN,TEMP,NNUU,NNLL,GIN,DGG,NTOT,NFAC,(PFAC(I), I=1,6)
  30. 150 FORMAT (2E10.2,I3,I2,2F10.2,2I5/6F10.5)
  31. IF (EOF,60) 577, 170
  32. C
  33. C CALCULATION OF MICROFIELD DISTRIBUTION FUNCTION
  34. 170 SHIELD = 0.0898 * DEN**(1./6.)/SQRT(TEMP)
  35. DO 161 J = 1,5
  36. 161 DST(J) = 0.2 * FLOAT(J-1)
  37. DO 163 I = 1,300
  38. DO 162 J = 1,5
  39. 162 DQQ(J) = W(I,J)
  40. CALL POLY5 (DST,DQQ,5,STRONG,SLUU,SHIELD,0.1,1)
  41. 163 FIELD(I+1) = SLUU(1)
  42. C
  43. C CALCULATION OF AVERAGE ION FIELD
  44. C
  45. DO 111 I = 1,300
  46. AM = I
  47. 111 FF(I) = FIELD(I+1) * AM/10.
  48. CALL WEDDLE (0.1,300,FF,BAV,0.)
  49. DY = 1./300.
  50. Y = 0.
  51. DO 113 I = 1,10
  52. Y = Y + DY
  53. B = 1./Y
  54. 113 FF(I) = WFLD(B) * B**3
  55. CALL WEDDLE (DY,10,FF,DB,0.)
  56. BAV = BAV + DB
  57. READ 175,NPUNCH,(NUR(I),I=1,2),IPI,NNCCQ
  58. 175 FORMAT (5I10)
  59. C
  60. C FOR NPUNCH > 1 K-MATRIX IS READ IN
  61. C FOR NPUNCH < 1 K-MATRIX IS CALCULATED
  62. C FOR NPUNCH = 1 K-MATRIX IS CALCULATED AND PUNCHED OUT
  63. C
  64. IF (NPUNCH.LE.1) GO TO 178
  65. DO 173 MBLOCK = 1,2
  66. NUPP = NUR(MBLOCK)
  67. NRRAB = (NUPP+1)*NUPP/2
  68. DO 174 J = IPI,NNCCQ
  69. READ 172, (FF(I),I=1,NRRAB)
  70. 172 FORMAT (5E16.9)
  71. DO 171 I = 1,NRRAB
  72. 171 SLOP(I,J,MBLOCK) = FF(I)
  73. 174 CONTINUE
  74. 173 CONTINUE
  75. C
  76. C CALCULATION OF WAVELENGTH IN STANDARD AIR
  77. C SOMETIMES I SERIOUSLY ASK MYSELF WHY I DO THESE THINGS
  78. C GENUINELY WHAT THE FUCK???
  79. C
  80. 178 AUU = NNUU
  81. ALL = NNLL
  82. SL1 = 109678.758 * (1./(ALL*ALL) - 1./(AUU*AUU))
  83. SL2 = SL1 * SL1
  84. ALDD = 1.000064328+2949810./(1.46E+10-SL2)+25540./(4.1E+9-SL2)
  85. ALAM = 1.E+8/(ALDD * SL1)
  86. PRINT 180, DEN, TEMP, NNUU, NNLL, ALAM, SHIELD
  87. 180 FORMAT /1H1,* DENSITY =* E10.2* TEMPERATURE =*E10.2,
  88. & * QUANTUMNUMBERS N UPPER =*I3* N LOWER =*I2* WAVELENGTH =*
  89. & F8.2* ANGSTROM*/10X,*SHIELDING PARAMETER =*F6.3,10X,
  90. & * INCLUDING LOWER STATE INTERACTION *//
  91. IF (NNUU.LE.5.AND.NNLL.LE.3) GO TO 200
  92. IF (NNUU.LE.8.AND.NNLL.LE.2) GO TO 200
  93. IF (NNUU.LE.16.AND.NNLL.EQ.1) GO TO 200
  94. PRINT 190
  95. 190 FORMAT (* OVERFLOW OF MATRICIES, PROGRAM NOT EXECUTED*)
  96. 577 CALL EXIT
  97. C
  98. C RADIAL MATRIXELEMENTS AND TOTAL LINESTRENGTH
  99. C
  100. 200 IPI = 2
  101. STOT = 0.
  102. DO 330 KLL = 1,NNLL
  103. LLL = KKL - 1
  104. ALL = LLL
  105. DO 310 KUU = 1,NNUU
  106. LUU = KUU - 1
  107. AUU = LUU
  108. RRMM = 0
  109. IF (IABS(KUU - KLL).NE.1) GO TO 310
  110. RRMM = RADMAT(NNUU,LUU,NNLL,LLL) * S3J(ALL,AUU,1.,0.,0.,0.)
  111. STOT = STOT + (2.*ALL+1.)*(2.*AUU+1.) * RRMM * RRMM
  112. 310 RDM(KUU,KLL) = RRMM
  113. 330 CONTINUE
  114. SQST = SQRT(STOT)
  115. C
  116. C TRANSFORMATION MATRICIES
  117. C
  118. DO 370 NKK = 1,2
  119. NN = NNLL
  120. IF (NKK.EQ.2) NN = NNUU
  121. ANN = NN
  122. AN1 = 0.5 * (ANN - 1.)
  123. KQQ = 2 * NN - 1
  124. NRUN = 0
  125. DO 370 KL = 1,NN
  126. AL = KL - 1
  127. FACC = SQRT(2.*AL+1.)
  128. IF(NKK.EQ.1) SLLL(KL) = FACC
  129. IF(NKK.EQ.2) SLUU(KL) = FACC
  130. DO 370 KM = 1,KL
  131. AM = KM - 1
  132. NRUN = NRUN + 1
  133. QQ = -ANN
  134. ALIM = ANN - 0.9 - AM
  135. DO 370 KQ = 1,KQQ
  136. QQ = QQ + 1.
  137. SIIJ = 0.
  138. ABQQ = ABS(QQ)
  139. IF (ABQQ.GT.ALIM) GO TO 360
  140. NQ = ABQQ + 0.1
  141. IF (MOD(NN+NQ+KM+1,2).NE.1) GO TO 360
  142. AMQ = 0.5 * (AM - QQ)
  143. APQ = 0.5 * (AM + QQ)
  144. SIIJ = FACC * S3J(AN1,AN1,AL,AMQ,APQ,-AM)
  145. 360 IF (NKK.EQ.1) SJLL(NRUN,KQ) = SIIJ
  146. IF (NKK.EQ.2) SJUU(NRUN,KQ) = SIIJ
  147. 370 CONTINUE
  148. C
  149. C DIPOLMATRIXELEMENTS IN PARABOLIC STATES AND MATRIX ORDERING
  150. C
  151. KQUU = 2 * NNUU - 1
  152. KQLL = 2 * NNLL - 1
  153. DO 500 MBLOCK = 1,2
  154. AMLCK = MBLOCK
  155. MBL = MBLOCK - 1
  156. AMBL = MBL
  157. NRUN = 0
  158. NCQ = 0
  159. NAQU = -NNUU
  160. DO 480 KUU = 1,KQUU
  161. NAQU = NAQU + 1
  162. KQU = NNUU + NAQU
  163. NAQL = -NNLL
  164. DO 480 KLL = 1,KQLL
  165. NAQL = NAQL + 1
  166. KQL = NNLL + NAQL
  167. NXDIFF = NNUU * NAQU - NNLL * NAQL
  168. AXDIF = NXDIFF
  169. IF (NXDIF.LT.0) GO TO 400
  170. NCQ = NCQ + 1
  171. IF (MBLOCK.EQ.2) GO TO 400
  172. NQCU(NCQ) = NAQU
  173. NQCL(NCQ) = NAQL
  174. NXX(NCQ) = NXDIF
  175. DQQ(NCQ) = 0.
  176. 400 NPHAS = (NNUU + NNLL + MBL - NAQU - NAQL - 2) / 2
  177. PHASE = (-1.)**MOD(NPHAS, 2)
  178. KMINL = NNLL - IABS(NAQL)
  179. KMINU = NNUU - IABS(NAQU)
  180. KML = -KMINL
  181. MKML = 2 * KMINL - 1
  182. DO 470 KAML = 1, MKML
  183. KML = KML + 1
  184. IF (MOD(NNLL+NAQL+KML,2).NE.1) GO TO 470
  185. KMU = KML + MBL
  186. IF (IABS(KMU).GE.KMINU) GO TO 470
  187. IF (MOD(NNUU+NAQU+KMU,2).NE.1) GO TO 470
  188. NRUN = NRUN + 1
  189. NQUU(NRUN,MBLOCK) = NAQU
  190. NQLL(NRUN,MBLOCK) = NAQL
  191. MMUU(NRUN,MBLOCK) = KMU
  192. MMLL(NRUN,MBLOCK) = KML
  193. ADIFF(NRUN,MBLOCK) = AXDIF
  194. AMU = KMU
  195. AML = KML
  196. KBBU = IABS(KMU) + 1
  197. KBBL = IABS(KML) + 1
  198. RESULT = 0.
  199. DO 450 KKLU = KBBU,NNUU
  200. ALU = KKLU - 1
  201. LLU = (KKLU - 1) * KKLU / 2 + KBBU
  202. FAC1 = SLUU(KKLU) * SJUU(LLU, KQU)
  203. DO 450 KKLL = KBBL, NNLL
  204. IF (IABS(KKLU - KKLL).NE.1) GO TO 450
  205. ALL = KKLL - 1
  206. LLL = (KKLL - 1) * KKLL / 2 + KBBL
  207. FAC2 = SLLL(KKLL) * SJLL(LLL, KQL)
  208. FACC = FAC2 * FAC1 * S3J(ALL, ALU, 1., AML, -AMU, AMBL)
  209. RESULT = RESULT + FACC * RDM(KKLU, KKLL)
  210. 450 CONTINUE
  211. SIIJ = RESULT * PHASE / SQST
  212. DIPOL(NRUN,MBLOCK) = SIIJ
  213. IF (NXDIF.GE.0) DQQ(NCQ) = DQQ(NCQ) + AMLCK * RESULT * RESULT / STOT
  214. PRINT 425, NAQU, KMU, NAQL, KML, SIIJ
  215. 425 FORMAT (16X,*(*,2I3,* I D I*,2I3,* ) =*E16.8)
  216. 470 CONTINUE
  217. 480 CONTINUE
  218. NUR(MBLOCK) = NRUN
  219. PRINT 393
  220. 393 FORMAT (/)
  221. 500 CONTINUE
  222. C
  223. C REORDERING OF THE STARK COMPONENTS
  224. C
  225. NCM1 = NCQ - 1
  226. DO 333 II = 1, NCM1
  227. IPP = II + 1
  228. DO 333 JJ = IPP, NCQ
  229. IF (NXX(II).LE.NXX(JJ)) GO TO 333
  230. TEMPP = NXX(II)
  231. NXX(II) = NXX(JJ)
  232. NXX(JJ) = TEMPP
  233. TEMPP = DQQ(II)
  234. DQQ(II) = DQQ(JJ)
  235. DQQ(JJ) = TEMPP
  236. TEMPP = NQCU(II)
  237. NQCU(II) = NQCU(JJ)
  238. NQCU(JJ) = TEMPP
  239. TEMPP = NQCL(II)
  240. NQCL(II) = NQCL(JJ)
  241. NQCL(JJ) = TEMPP
  242. 333 CONTINUE
  243. PRINT 410, (NQCU(I), NQCL(I), NXX(I), DQQ(I), I = 1, NCQ)
  244. 410 FORMAT(15X,* Q UPPER = *I3,* Q LOWER =*I3,* X =*I5,*
  245. & DIPOL=* E16.8)
  246. PRINT 393
  247. JJ = 1
  248. KMST = 0
  249. 315 KMST = KMST + 1
  250. DST(KMST) = 0.
  251. DO 325 II = JJ, NCQ
  252. IF (NXX(II) .NE. NXX(JJ)) GO TO 325
  253. NST(KMST) = NXX(II)
  254. DST(KMST) = DST(KMST) + DQQ(II)
  255. J = II
  256. 325 CONTINUE
  257. JJ = J + 1
  258. IF (JJ .GT. NCQ) GO TO 345
  259. GO TO 315
  260. 345 PRINT 355, (NST(I), DST(I), I = 1, KMST)
  261. 355 FORMAT (16X, *X =*I3,F20.10)
  262. PRINT 393
  263. C
  264. C K OPERATOR MATRIX
  265. C
  266. IF (NPUNCH - 1) 445, 455, 650
  267. 455 IPP = 2
  268. PUNCH 465, IPP, (NUR(I), I=1,2), IPI, NCQ, NNUU, NNLL
  269. 465 FORMAT (5I10, 15X, 2I5)
  270. 445 NTIME = KLOCK(0)
  271. DO 600 MBLOCK = 1, 2
  272. MBL = 1 - MBLOCK
  273. AMBL = MBL
  274. NUPP = NUR(MBLOCK)
  275. DO 600 NRC = IPI, NCQ
  276. NCQU = NQCU(NRC)
  277. NCQL = NQCL(NRC)
  278. KCQU = NNUU + NCQU
  279. KCQL = NNLL + NCQL
  280. IMULIM = NNUU - IABS(NCQU)
  281. IMLLIM = NNLL - IABS(NCQL)
  282. DO 590 NRA = 1, NUPP
  283. NAQU = NQUU(NRA, MBLOCK)
  284. KAMU = MMUU(NRA, MBLOCK)
  285. AMU = KAMU
  286. IABAMU = IABS(KAMU)
  287. NAQL = NQLL(NRA, MBLOCK)
  288. KAML = MMLL(NRA, MBLOCK)
  289. AML = KAML
  290. IABAML = IABS(KAML)
  291. KAQU = NNUU + NAQU
  292. KAQL = NNLL + NAQL
  293. KAAU = IABS(KAMU) + 1
  294. KAAL = IABS(KAML) + 1
  295. DO 580 NRB = 1, NRA
  296. NRRAB = (NRA - 1) * NRA / 2 + NRB
  297. NBQU = NQUU(NRB, MBLOCK)
  298. KBMU = MMUU(NRB, MBLOCK)
  299. BMU = KBMU
  300. IABBMU = IABS(KBMU)
  301. NBQL = NQLL(NRB, MBLOCK)
  302. KBML = MMLL(NRB, MBLOCK)
  303. BML = KBML
  304. IABBML = IABS(KBML)
  305. KBQU = NNUU + NBQU
  306. KBQL = NNLL + NBQL
  307. KBBU = IABS(KBMU) + 1
  308. KBBL = IABS(KBML) + 1
  309. NPHAS = MBL - NCQU - NCQL - (NAQU + NAQL + NBQU + NBQL) / 2
  310. SIIJ = 0.
  311. DO 540 JLAL = KAAL, NNLL
  312. ALL = JLAL - 1
  313. LALL = (JLAL - 1) * JLAL / 2 + IABAML + 1
  314. FAC1 = SJLL(LALL, KAQL)
  315. DO 540 JLAU = KAAU, NNUU
  316. ALU = JLAU - 1
  317. LALU = (JLAU - 1) * JLAU / 2 + IABAMU + 1
  318. FAC2 = SJUU(LALU, KAQU)
  319. FAC12 = FAC1 * FAC2
  320. DO 540 JLBL = KBBL, NNLL
  321. BLL = JLBL - 1
  322. LBLL = (JLBL - 1) * JLBL / 2 + IABBML + 1
  323. FAC3 = SJLL(LBLL, KBQL)
  324. FAC13 = FAC12 * FAC3
  325. MCLUP = MINO(JLBL, JLAL) - 1
  326. DO 540 JLBU = KBBU, NNUU
  327. BLU = JLBU - 1
  328. LBLU = (JLBU - 1) * JLBU / 2 + IABBMU + 1
  329. FAC4 = SJUU(LBLU, KBQU)
  330. FAC14 = FAC13 * FAC4
  331. MCUUP = MINO(JLBU, JLAU) - 1
  332. IUPP = MINO(IABS(JLAU + JLAL), IABS(JLBU + JLBL)) - 1
  333. MCL = -IMLLIM - 1
  334. DO 540 JMCLL = 1, IMLLIM
  335. MCL = MCL + 2
  336. IABCML = IABS(MCL)
  337. IF (IABCML .GT. MCLUP) GO TO 540
  338. AMCL = MCL
  339. LALL = (JLAL - 1) * JLAL / 2 + IABCML + 1
  340. LBLL = (JLBL - 1) * JLBL / 2 + IABCML + 1
  341. FAC5 = SJLL(LALL, KCQL) * SJLL(LBLL, KCQL)
  342. FAC15 = FAC14 * FAC5
  343. IF (ABS(FAC15) .LT. 1.E-8) GO TO 540
  344. MCU = -IMULIM - 1
  345. DO 535 JMCUU = 1, IMULIM
  346. MCU = MCU + 2
  347. IABCMU = IABS(MCU)
  348. IF (IABCMU .GT. MCUUP) GO TO 535
  349. AMCU = MCU
  350. LALU = (JLAU - 1) * JLAU / 2 + IABCMU + 1
  351. LBLU = (JLBU - 1) * JLBU / 2 + IABCMU + 1
  352. FAC6 = SJUU(LALU, KCQU) * SJUU(LBLU, KCQU)
  353. IF (ABS(FAC6) .LT. 1.E-8) GO TO 535
  354. MGGL = MCL - MCU
  355. AMGL = MGGL
  356. ILOW = MAXO(MBL, IABS(MGGL), IABS(JLAU - JLAL), IABS(JLBU - JLBL)) + 1
  357. IF (ILOW .GT. IUPP) GO TO 535
  358. FAC16 = FAC15 * FAC6 * ((-1.) ** MOD(NPHAS - MGGL, 2))
  359. DO 530 JLDD = ILOW, IUPP
  360. IF (MOD(JLAL + JLAU + JLDD, 2) .NE. 0) GO TO 515
  361. IF (KAMU .EQ. 0 .AND. KAML .EQ. 0) GO TO 530
  362. IF (MCL .EQ. 0 .AND. MCU .EQ. 0) GO TO 530
  363. 515 IF (MOD(JLBL + JLBU + JLDD, 2) .NE. 0) GO TO 520
  364. IF (KBMU .EQ. 0 .AND. KBML .EQ. 0) GO TO 530
  365. IF (MCL .EQ. 0 .AND. MCU .EQ. 0) GO TO 530
  366. 520 ALDD = JLDD - 1
  367. FAC7 = S3J(ALL, ALU, ALDD, -AMCL, AMCU, AMGL) * (2. * ALDD + 1.)
  368. FAC7 = FAC7 * S3J(ALL, ALU, ALDD, -AML, AMU, AMBL)
  369. FAC7 = FAC7 * S3J(BLL, BLU, ALDD, -AMCL, AMCU, AMGL)
  370. FACC = FAC7 * S3J(BLL, BLU, ALDD, -BML, BMU, AMBL) * FAC16
  371. SIIJ = SIIJ + FACC
  372. 530 CONTINUE
  373. 535 CONTINUE
  374. 540 CONTINUE
  375. SLOP(NRRAB, NRC, MBLOCK) = SIIJ
  376. 580 FF(NRRAB) = SIIJ
  377. 590 CONTINUE
  378. IF (NPUNCH .EQ. 1) PUNCH 630, (FF(I), I=1, NRRAB)
  379. 630 FORMAT (5E16.9)
  380. 600 CONTINUE
  381. NTIME = -NTIME + KLOCK(0)
  382. PRINT 617, NTIME
  383. 617 FORMAT (16X, *COMPUTER TIME FOR CALCULATING K-OPERATOR &
  384. & MATRIX =*I8/)
  385. C
  386. C BASIC CONSTANTS AND ARRAY FOR G-FUNCTION CONSTANTS
  387. C
  388. 650 SDEN = SQRT(DEN)
  389. FAC = 2064.936 * TEMP * SQRT(TEMP / DEN)
  390. CFAC = 4.5645E-7 * SDEN / TEMP
  391. DEBROG = 2.1027E-6 / SQRT(TEMP)
  392. ANN = NNUU
  393. RMIN = DEBROG + ANN * ANN * 7.9376E-9
  394. BET = 5.6558E-5 * DEN ** (1. / 6.)
  395. ASY = 0.0
  396. DO 270 NRC = IPI, NCQ
  397. QC = NXX(NRC)
  398. C = CFAC * QC
  399. P1 = -1.671086 * FAC * C * SQRT(C)
  400. BS = 3. * QC * DEBROG / RMIN
  401. STRONG(NRC) = 0.269-2. * (((1. - COSF(BS)) / BS + SINF(BS)) / BS
  402. & - COSINT(BS))
  403. PPFF = -1.128379 * FAC * C * C
  404. P2 = PPFF * (STRONG(NRC) - 2. * LOGF(2. * C))
  405. FPAR(1, NRC) = P1
  406. FPAR(2, NRC) = P2
  407. FPAR(3, NRC) = 0.5 * (P2 / P1) ** 2
  408. FIN = LOGF(QC)
  409. A2 = P2 * ((PFAC(3) * FIN * PFAC(2)) * FIN + PFAC(1))
  410. B2 = (PFAC(6) * FIN + PFAC(5)) * FIN + PFAC(4)
  411. IF (B2 .LT. 0.) A2 = 0.
  412. FPAR(4, NRC) = A2
  413. FPAR(5, NRC) = B2
  414. FPAR(6, NRC) = PPFF
  415. 270 ASY = ASY + 2. * P1 * DQQ(NRC)
  416. PRINT 220
  417. 220 FORMAT (/13X*P1*18X*P2*18X*B1*18X*A2*18X*B2*17X*STRONG*/)
  418. PRINT 240, ((FPAR(K, I), K = 1, 5), STRONG(I), I = IPI, NCQ)
  419. 240 FORMAT (6E20.4)
  420. PRINT 280, FAC, CFAC, BET, AST, DEBROG, BAV, NFAC
  421. 280 FORMAT /* FAC =*E11.4,* CFAC =*E11.4,* BET =*E11.4,* ASY =*E11.4,
  422. & * DEBROG =E11.4,* BAV =*F7.4,* INTEGRATIONFACTOR =*I2//
  423. & 5X*DOM*8X*DLAM*8X*ITOT*6X*IHOLTS*8X*ASY*10X*WING*7X*WHOLTS*7X
  424. & *WWOO*8X*WWBB*8X*WWNG*8X*TIME STEPS*/
  425. ADLFAC = 4.23538E-15 * SDEN * ALAM * ALAM
  426. C
  427. C CALCULATION OF THE IONFIELD INTEGRAL
  428. C
  429. N12 = 12 * NFAC
  430. AN12 = N12
  431. N30 = 30 * NFAC
  432. AN30 = N30
  433. DMCRT = BET * BAV * NXX(NCQ) * 5.
  434. G = GIN - DGG
  435. DO 950 MM = 1, NTOT
  436. NTIME = KLOCK(0)
  437. NSTEP = 0
  438. G = G + DGG
  439. DOM = 10. ** G
  440. DLAM = ADLFAC * DOM
  441. WING52 = -0.2992067103 * ASY / (SQRT(DOM) * DOM * DOM)
  442. FHOLTS = 0.
  443. AWING = 0.
  444. DO 815 NRC = IPI, NCQ
  445. P1 = FPAR(1, NRC)
  446. P2 = FPAR(2, NRC)
  447. B1 = FPAR(3, NRC)
  448. A2 = FPAR(4, NRC)
  449. B2 = FPAR(5, NRC)
  450. PPFF = FPAR(6, NRC)
  451. CALL FOUTR(DOM)
  452. AWING = AWING + DQQ(NRC) * CCIMAG * 2. / (DOM * DOM)
  453. QC = NXX(NRC)
  454. BETFAC = BET * QC
  455. BCRIT = DOM / BETFAC
  456. FHOLTS = FHOLTS + DQQ(NRC) * WFLD(BCRIT) / BETFAC
  457. IF (DOM .GT. DMCRT) GO TO 985
  458. AIRES = 0.
  459. IF (DOM .GT. (-3. * P2)) GO TO 840
  460. ANQ1 = NXX(NCQ)
  461. BCRIT = (DOM - P2) / (ANQ1 * BET)
  462. DB = BCRIT / AN12
  463. B = 0.
  464. DO 820 J = 1, N12
  465. B = B + DB
  466. 820 FF(J) = AIIM(DOM, B) * WFLD(B)
  467. CALL WEDDLE(DB, N12, FF, AIII, 0.)
  468. AIRES = AIII
  469. DY = 1. / (BCRIT * AN30)
  470. Y = 0.
  471. DO 830 J = 1, N30
  472. Y = Y + DY
  473. B = 1. / Y
  474. 830 FF(J) = B * B * AIIM(DOM, B) * WFLD(B)
  475. CALL WEDDLE(DY, N30, FF, AIII, 0.)
  476. AIRES = AIRES + AIII
  477. GO TO 980
  478. 840 BCRCR = DOM / BET
  479. EPSPS = -P2 / BET
  480. DO 957 NQ = IPI, KMST
  481. ANQ = NST(NQ)
  482. BCR = BCRCR / ANQ
  483. EPS = EPSPS / ANQ
  484. IF (NQ .EQ. IPI) GO TO 907
  485. SL1 = 1. / (GAM - BCR)
  486. GO TO 908
  487. 907 SL1 = 0.
  488. 908 SL2 = 1. / EPS
  489. SL3 = 1. / (BCR + EPS)
  490. SL4 = 1. / (BCR - EPS)
  491. IF (NQ .EQ. KMST) GO TO 911
  492. ANQ1 = NST(NQ + 1)
  493. GAM = 0.5 * (BCR - EPS + (BCRCR + EPSPS) / ANQ1)
  494. GO TO 912
  495. 911 GAM = 0.5 * (BCR - EPS)
  496. 912 SL5 = 1. / (BCR - GAM)
  497. CRIT = SL2 - SL5
  498. 904 Y = SL1
  499. IF (NQ .EQ. IPI) GO TO 913
  500. B = BCR + 1. / Y
  501. FA = AIIM(DOM, B) * WFLD(B) / (Y * Y)
  502. GO TO 914
  503. 913 FA = 0.
  504. 914 DY = (SL2 - SL1) / AN12
  505. DO 917 J = 1, N12
  506. Y = Y + DY
  507. Y1 = 1. / Y
  508. B = BCR + Y1
  509. 917 FF(J) = Y1 * Y1 * AIIM(DOM, B) * WFLD(B)
  510. CALL WEDDLE(DY, N12, FF, AIII, FA)
  511. AIRES = AIRES + AIII
  512. Y = SL3
  513. B = 1. / Y
  514. FA = B * B * AIIM(DOM, B) * WFLD(B)
  515. DY = (SL4 - SL3) / AN12
  516. DO 927 J = 1, N12
  517. Y = Y + DY
  518. B = 1. / Y
  519. 927 FF(J) = B * B * AIIM(DOM, B) * WFLD(B)
  520. CALL WEDDLE(DY, N12, FF, AIII, FA)
  521. AIRES = AIRES + AIII
  522. IF (CRIT .LE. 0.) GO TO 977
  523. Y = SL5
  524. B = BCR - 1. / Y
  525. FA = AIIM(DOM, B) * WFLD(B) / (Y * Y)
  526. DY = CRIT/AN12
  527. DO 937 J = 1, N12
  528. Y = Y + DY
  529. Y1 = 1. / Y
  530. B = BCR - Y1
  531. 937 FF(J) = Y1 * Y1 * AIIM(DOM, B) * WFLD(B)
  532. CALL WEDDLE(DY, N12, FF, AIII, FA)
  533. AIRES = AIRES + AIII
  534. 957 CONTINUE
  535. IF (GAM .LT. 5.) GO TO 968
  536. Y = 1. / GAM
  537. DY = (0.2 - Y) / AN12
  538. FA = GAM * GAM * AIIM(DOM, GAM) * WFLD(GAM)
  539. DO 967 J = 1, N12
  540. Y = Y + DY
  541. B = 1. / Y
  542. 967 FF(J) = B * B * AIIM(DOM, B) * WFLD(B)
  543. CALL WEDDLE(DY, N12, FF, AIII, FA)
  544. AIRES = AIRES + AIII
  545. SL4 = 0.2
  546. GO TO 977
  547. 968 SL4 = 1. / GAM
  548. 977 B = 0.
  549. DB = 1. / (SL4 * AN30)
  550. DO 947 J = 1, N30
  551. B = B + DB
  552. 947 FF(J) = AIIM(DOM, B) * WFLD(B)
  553. CALL WEDDLE(DB, N30, FF, AIII, 0.)
  554. AIRES = AIRES + AIII
  555. 980 WWBB = (AIIM(DOM, BAV) + FHOLTS) / WING52
  556. GO TO 990
  557. 985 AIRES = AIIM(DOM, BAV) + FHOLTS
  558. WWBB = 0.
  559. 990 WING = AIRES / WING52
  560. WINHOL = FHOLTS / WING52
  561. WWOO = (AIIM(DOM, 0.) + FHOLTS) / WING52
  562. WWNG = (AWING + FHOLTS) / WING52
  563. NTIME = -NTIME + KLOCK(0)
  564. 950 PRINT 978, /DOM, DLAM, AIRES, FHOLTS, WING52, WING, WINHOL, WWOO, WWBB,
  565. & WWNG, NTIME, NSTEP/
  566. 978 FORMAT (10E12.4, I10, I6)
  567. GO TO 120
  568. END
  569. C
  570. C
  571. FUNCTION AIIM(DOM, B)
  572. C
  573. C CALCULATION OF I(DOM,B) FOR THE GENERAL CASE INCLUDING LOWER
  574. C STATE INTERACTION
  575. DIMENSION DRRF(34), DIIF(34), FMATR(34,34)
  576. COMMON /FDAT/ P1, P2, B1, A2, B2, PPFF, CCREAL, CCIMAG
  577. COMMON /11/ SLOP(595,22,2)
  578. COMMON /22/ DIPOL(34,2), NUR(2), IPI, NCQ, BET, FPAR(6,30), &
  579. & ADIFF(34,2), NSTEP
  580. COMMON /33/ AMATR(34,68), BMATR(34,34), CMATR(34,34)
  581. AIIM = 0.0
  582. NSTEP = NSTEP + 1
  583. DO 800 MBL = 1, 2
  584. AML = MBL
  585. NUPP = NUR(MBL)
  586. DO 750 NRB = 1, NUPP
  587. DOMRB = DOM - BET * B * ADIFF(NRB, MBL)
  588. DO 220 NRC = IPI, NCQ
  589. P1 = FPAR(1, NRC)
  590. P2 = FPAR(2, NRC)
  591. B1 = FPAR(3, NRC)
  592. A2 = FPAR(4, NRC)
  593. B2 = FPAR(5, NRC)
  594. PPFF = FPAR(6, NRC)
  595. CALL FOUTR(DOMRB)
  596. DRRF(NRC) = CCREAL
  597. 220 DIIF(NRC) = CCIMAG
  598. DO 700 NRA = 1, NUPP
  599. IF (NRA .LE. NRB) NRRAB = (NRB - 1) * NRB / 2 + NRA
  600. IF (NRA .GT. NRB) NRRAB = (NRA - 1) * NRA / 2 + NRB
  601. ARRR = 0.0
  602. AIII = 0.0
  603. DO 600 NRC = IPI, NCQ
  604. TEMP = SLOP(NRRAB, NRC, MBL)
  605. ARRR = ARRR + TEMP * DRRF(NRC)
  606. 600 AIII = AIII + TEMP * DIIF(NRC)
  607. AIII = AIII * 6.2831853072
  608. AMATR(NRB, NRA) = AIII
  609. CMATR(NRB, NRA) = AIII
  610. 700 BMATR(NRB, NRA) = ARRR * 6.2831853072
  611. 750 BMATR(NRB, NRB) = BMATR(NRB, NRB) + DOMRB
  612. CALL ENVERS(NNUPP)
  613. DO 320 NRB = 1, NUPP
  614. DO 320 NRA = 1, NUPP
  615. TEMP = 0.
  616. DO 300 NRC = 1, NUPP
  617. 300 TEMP = TEMP + BMATR(NRB, NRC) * AMATR(NRC, NRA + NUPP)
  618. 320 FMATR(NRB, NRA) = TEMP
  619. DO 420 NRB = 1, NUPP
  620. DO 420 NRA = 1, NUPP
  621. TEMP = 0.
  622. DO 400 NRC = 1, NUPP
  623. 400 TEMP = TEMP + FMATR(NRB, NRC) * BMATR(NRC, NRA)
  624. 420 AMATR(NRB, NRA) = TEMP + CMATR(NRB, NRA)
  625. CALL ENVERS(NUPP)
  626. DO 795 NRB = 1, NUPP
  627. DO 795 NRA = 1, NUPP
  628. 795 AIIM = AIIM + AML * DIPOL(NRB, MBL) * DIPOL(NRA, MBL) &
  629. & * AMATR(NRB, NRA + NUPP)
  630. 800 CONTINUE
  631. AIIM = AIIM * 0.3183099
  632. RETURN
  633. END
  634. C
  635. C
  636. SUBROUTINE POLY5 (X,Y,NUMX,XN,YN,XNI,DX,NXNUM)
  637. C A 5 POINT POLYNOMIAL INTERPOLATION ROUTINE
  638. C FOR EACH POINT TO BE INTERPOLATED THE 5 NEAREST KNOWN POINTS
  639. C ARE CHOSEN, AND A 5TH DEGREE POLYNOMIAL IS FITTED TO THEM.
  640. C X AND Y ARE THE ARRAYS OF KNOWN POINTS ON THE CURVE
  641. C NUMX IS THE NUMBER OF KNOWN POINTS
  642. C TWO NEW ARRAYS, XN AND YN, WILL BE GENERATED FOR THE NXNUM
  643. C VALUES OF X AND F(X) STARTING WITH X=XNI IN INCREMENTS OF DX
  644. C
  645. DIMENSION X(5),Y(5),XN(5),YN(5)
  646. J=2
  647. DO 8 I=1,NXNUM
  648. EYE=I-1
  649. XN(I)=XNI+DX*EYE
  650. IF (XN(I)-X(1)).LT.0 GO TO 1
  651. IF (XN(I)-X(1)).EQ.0 GO TO 2
  652. IF (XN(I)-X(1)).GT.0 GO TO 3
  653. 2 YN(I)=Y(1)
  654. GO TO 8
  655. 3 IF (XN(I).LE.X(J)) GO TO 4
  656. IF (J.GE.NUMX) GO TO 7
  657. J = J + 1
  658. GO TO 3
  659. 1 L = 1
  660. GO TO 70
  661. 7 L = NUMX - 4
  662. GO TO 70
  663. 4 L = J - 1
  664. IF (J.GT.2) L = J - 2
  665. IF (J.GT.3) L = J - 3
  666. IF ((J+1).GT.NUMX) L = J - 4
  667. 70 YN(I)=0.0
  668. LLL=L+4
  669. DO 75 K=L,LLL
  670. TERM=1.0
  671. DO 74 M=L,LLL
  672. IF (K.EQ.M) GO TO 74
  673. TDEN=X(K)-X(M)
  674. TNUM=XN(I)-X(M)
  675. TERM=TERM*TNUM/TDEN
  676. 74 CONTINUE
  677. TERM = Y(K)*TERM
  678. 75 YN(I)=YN(I)+TERM
  679. 8 CONTINUE
  680. RETURN
  681. END
  682. C
  683. C
  684. SUBROUTINE FOUTR (DOM)
  685. C FOURIER TRANSFORM OF THERMAL AVERAGE FOR UNIFIED THEORY
  686. C
  687. COMMON /FDAT/ P1, P2, B1, A2, B2, PPFF, CCREAL, CCIMAG
  688. ARG = ABSF(DOM)
  689. Z = B1 * ARG
  690. IF (Z .LE. 0.001) GO TO 600
  691. AR2 = ARG * ARG
  692. IF (Z .LE. 40.) GO TO 300
  693. FAC1 = -0.2992067103 * P1 / (SQRT(ARG) * AR2)
  694. CC = FAC1 * ((1. - 1.3125 / Z) * 0.625 / Z + 1.)
  695. SS = CC - FAC1 * 1.25 / Z
  696. GO TO 500
  697. 300 CALL BSJY01(Z, AJ0, Y0, AJ1, Y1)
  698. FAC1 = Y1 / (2. * Z) + AJ1 - Y0
  699. FAC2 = AJ0 + Y1 - AJ1 / (2. * Z)
  700. CINE = COSF(Z)
  701. SINE = SINF(Z)
  702. FCC = P2 * B1 * B1
  703. CC = FCC * (CINE * FAC1 + SINE * FAC2)
  704. SS = FCC * (CINE * FAC2 - SINE * FAC1)
  705. IF (A2 .EQ. 0.) GO TO 500
  706. Z = B2 * ARG
  707. IF (Z .GT. 10.) GO TO 400
  708. CALL BSJY01(Z, AJ0, Y0, AJ1, Y1)
  709. FAC1 = ((AJ1 - Y0) * 16. * Z - 36. * AJ0 - 28. * Y1) * Z &
  710. & + 15. * Y0 - 3. * AJ1
  711. FAC2 = ((AJ0 + Y1) * 16. * Z - 36. * Y0 + 28. * AJ1) * Z &
  712. & - 15. * AJ0 - 3. * Y1
  713. CINE = COSF(Z)
  714. SINE = SINF(Z)
  715. FCC = A2 * B2 / 6.
  716. CC = CC + FCC * (CINE * FAC1 + SINE * FAC2)
  717. SS = SS + FCC * (CINE * FAC2 - SINE * FAC1)
  718. GO TO 500
  719. 400 FAC1 = 0.1322319336 * A2 * B2 * Z ** (-3.5)
  720. FCC = FAC1 * (1. - (3.9375 / Z + 1.) * 4.375 / Z)
  721. CC = CC + FCC
  722. SS = SS - FCC - FAC1 * 8.75 / Z
  723. 500 IF (DOM .LT. 0.) SS = -SS
  724. CCREAL = -AR2 * SS
  725. CCIMAG = AR2 * CC
  726. RETURN
  727. 600 CCREAL = 0.3183099 * (P2 * B1 - A2) * DOM
  728. CCIMAG = -0.3183099 * P2
  729. RETURN
  730. END
  731. C
  732. C
  733. SUBROUTINE ENVERS (NN)
  734. C
  735. COMMON /33/ C(34,68), BMATR(34,34), CMATR(34,34)
  736. N = NN
  737. IF (N .GT. 1) GO TO 10
  738. C(1,2) = 1.0 / C(1,1)
  739. RETURN
  740. 10 DO 300 K = 1, N
  741. KN = K + N
  742. NP = KN - 1
  743. KP1 = K + 1
  744. TEMP = 1. / C(K,K)
  745. DO 100 I = KP1, NP
  746. 100 C(K,I) = C(K,I) * TEMP
  747. C(K,KN) = TEMP
  748. DO 300 L = 1, N
  749. IF (K .EQ. L) GO TO 300
  750. TEMP = C(L,K)
  751. DO 200 I = KP1, NP
  752. 200 C(L,I) = C(L,I) - C(K,I) * TEMP
  753. C(L,KN) = -C(K,KN) * TEMP
  754. 300 CONTINUE
  755. RETURN
  756. END
  757. C
  758. C
  759. FUNCTION RADMAT(N1,L1,N2,L2)
  760. C
  761. C
  762. DIMENSION F(2)
  763. RADMAT = 0.0
  764. IF(IABS(L1 - L2).NE.1) RETURN
  765. NA = N1
  766. NB = N2
  767. L = L1
  768. IF(L1.GT.L2) GO TO 10
  769. NA = N2
  770. NB = N1
  771. L = L2
  772. 10 ANA = NA
  773. BNB = NB
  774. AL = L
  775. IF(N1.EQ.N2) GO TO 60
  776. RADMAT = FCTRL(ANA + AL) * FCTRL(BNB + AL - 1.)
  777. RADMAT = RADMAT/(FCTRL(ANA - AL - 1.) * FCTRL(BNB - AL))
  778. RADMAT = SQRT(RADMAT)
  779. I = BNB + AL
  780. AL = 2. * AL
  781. RADMAT = RADMAT * ((-1.0)**XMODF(I,2))/(4.*FCTRL(AL - 1.))
  782. ABAB = (ANA - BNB)/(ANA + BNB)
  783. AB2 = 1./(ANA - BNB)**2
  784. RADMAT = RADMAT * ABAB**(NA+NB)*(4.*ANA*BNB*AB2)**(L+1)
  785. Z = -4. * ANA * BNB * AB2
  786. NR = NA - L - 1
  787. NRP = NB - L
  788. ANRP = NRP
  789. DO 30 MM = 1,2
  790. F(MM) = 1.0
  791. IF(MM.EQ.2) NR = NR + 2
  792. K = MINO (NR,NRP)
  793. IF(K.EQ.0) GO TO 30
  794. ANR = NR
  795. PROD = 1.0
  796. DO 40 M = 1,K
  797. AM = M - 1
  798. PROD = PROD * (AM-ANR)*(AM-ANRP)*Z/((AL+AM)*(AM+1.))
  799. 40 F(MM) = F(MM) + PROD
  800. 30 CONTINUE
  801. RADMAT = RADMAT * (F(1) - ABAB * ABAB * F(2))
  802. RETURN
  803. 60 RADMAT = -1.5 * ANA * SQRT((ANA+AL)*(ANA-AL))
  804. RETURN
  805. END
  806. C
  807. C
  808. FUNCTION FCTRL(A)
  809. C
  810. C CALCULATION OF FACTORIALS
  811. DIMENSION FCTI(20)
  812. DATA (FCTI(I),I=1,20) / 1.0,2.0,6.0,24.0,120.0,720.0,5040.0,
  813. & 40320.0,362880.0,3628800.0,39916800.0,479001600.0,
  814. & 6227020800.0,87178291200.0,1307674368000.0,
  815. & 2.0922789888E13, 3.55687428096E14, 6.402373705728E15,
  816. & 1.2164510040883E17, 2.4329020081766E18/
  817. I = A + 0.1
  818. IF(I) 50,60,70
  819. 50 FCTRL = 0.0
  820. RETURN
  821. 60 FCTRL = 1.0
  822. RETURN
  823. 70 IF (I.GT.20) GO TO 130
  824. FCTRL=FCTI(I)
  825. RETURN
  826. 130 F=20.0
  827. FCTRL=FCTI(20)
  828. DO 131 J=21,I
  829. F=F+1.0
  830. 131 FCTRL=FCTRL*F
  831. RETURN
  832. END
  833. C
  834. C
  835. FUNCTION S3J (FJ1, FJ2, FJ3, FM1, FM2, FM3)
  836. C
  837. C CALCULATION OF THE 3J-SYMBOL, FJ1+FJ2+FJ3 HAS TO BE LESS THAN 38
  838. C
  839. DIMENSION GAM(39)
  840. DATA (GAM(I),I=1,39) / 1.,1.,2.,6.,24.,120.,720.,5040.,40320.,
  841. & 362880.,3628800.,39916800.,479001600.,6227020800.,87178291200.,
  842. & 1307674368000.,2.0922789888E+13,3.5568742810E+14,
  843. & 6.4023737057E+15,1.2164510041E+17,2.4329020082E+18,
  844. & 5.1090942172E+19,1.1240007278E+21,2.5852016739E+22,
  845. & 6.2044840173E+23,1.5511210043E+25,4.0329146113E+26,
  846. & 1.0888869450E+28,3.0488834461E+29,8.8417619937E+30,
  847. & 2.6525285981E+32,8.2228386542E+33,2.6313083693E+35,
  848. & 8.6833176188E+36,2.9523279904E+38,1.0333147966E+40,
  849. & 3.7199332679E+41,1.3763753091E+43,5.2302261747E+44/
  850. S3J=0.0
  851. FJA = FJ1+0.001
  852. FJB= FJ2+0.001
  853. FJC = FJ3+0.001
  854. IF (FJA.LT.0.0.OR.FJB.LT.0.0.OR.FJC.LT.0.0) RETURN
  855. IF (AMOD(FJA,0.5).GT.0.002) RETURN
  856. IF (AMOD(FJB,0.5).GT.0.002) RETURN
  857. IF (AMOD(FJC,0.5).GT.0.002) RETURN
  858. FMA = FM1
  859. FMB = FM2
  860. FMC = FM3
  861. IF (ABS(FMA).GT.FJA.OR.ABS(FMB).GT.FJB.OR.ABS(FMC).GT.FJC) RETURN
  862. B = FJA - FMA
  863. P = FJB - FMB
  864. D = FJC - FMC
  865. G = FJA + FJB - FJC
  866. IF ((ABS(FMA+FMB+FMC)+AMOD(B,1.)+AMOD(P,1.)+AMOD(D,1.)+DIM(0.,G)
  867. & +DIM(ABS(FJA-FJB),FJC)).GT.0.01) RETURN
  868. IHH = FJA + FJB + FJC + 2.
  869. IAA = FJB+FJC+FMA+1.
  870. ICC = -FJA+FJB+FJC+1.
  871. IFF = FJA-FJB+FJC+1.
  872. IGG = G+1.
  873. J1M = B+1.
  874. J2M = P+1.
  875. J3M = D+1.
  876. J1P = FJA+FMA+1.
  877. J2P = FJB+FMB+1.
  878. J3P = FJC+FMC+1.
  879. E2 = GAM(J1M)*GAM(J1P)*GAM(J2M)*GAM(J2P)
  880. E1 = (GAM(ICC)*GAM(IFF)/GAM(IHH))*GAM(IGG)*GAM(J3M)*GAM(J3P)
  881. E1 = SQRT(E2/E1)
  882. E = FJA - FJB + FMC
  883. IAEE = ABS(E) + 0.001
  884. IEE = IAEE + 1
  885. I1 = IAEE
  886. 100 I2 = ICC - 1
  887. IF (J3M.LT.ICC) I2 = J3M - 1
  888. DO 150 I=I1,I2
  889. E2 = GAM(I+1)*GAM(ICC-I)*GAM(J3M-I)/GAM(IAA-I)
  890. 150 S3J = S3J+(((-1.)**(I-(I/2)*2))/E2)*GAM(J1M+I)/GAM(IEE+I)
  891. I = ABS(FJA+FMB-FMC) + 0.01
  892. S3J = S3J * ((-1.)**(I-(I/2)*2))/E1
  893. RETURN
  894. END
  895. C
  896. C
  897. SUBROUTINE WEDDLE (DX,N,F,A,F0)
  898. C
  899. C INTEGRATION SUBROUTINE
  900. C
  901. DIMENSION F(N)
  902. A = 0.0
  903. K = N - 1
  904. DO 15 I = 1,6
  905. SM = 0.0
  906. DO 6 J = I,K,6
  907. 6 SM = SM + F(J)
  908. GO TO (8, 10, 12, 10, 8, 14), I
  909. 8 A = A + 5.0 * SM
  910. GO TO 15
  911. 10 A = A + SM
  912. GO TO 15
  913. 12 A = A + 6.0 * SM
  914. GO TO 15
  915. 14 A = A + 2.0 * SM
  916. 15 CONTINUE
  917. A = 0.3 * DX * (A + F0 + F(N))
  918. RETURN
  919. END
  920. C
  921. FUNCTION WFLD(B)
  922. C
  923. C CALCULATION OF THE ION MICROFIELD DISTRIBUTION FUNCTION USING
  924. C 5POINT INTERPOLATION FOR THE DATA READ INTO THE MAIN PROGRAM
  925. C
  926. COMMON /PFW/ FIELD(301)
  927. WFLD = 0.0
  928. IF (B.LE.30.0) GO TO 200
  929. SBS = 1./(B * SQRT(B))
  930. WFLD = ((21.6 * SBS + 7.639) * SBS + 1.496) * SBS/B
  931. RETURN
  932. 200 IF (B.LE.0.0) RETURN
  933. J = (B + 0.2) * 10.0
  934. L = J - 1
  935. IF (J.GT.2) L = J - 2
  936. IF (J.GT.3) L = J - 3
  937. IF (J.GT.300) L = 297
  938. 70 LLL = L + 4
  939. DO 75 K = L,LLL
  940. AK = K - 1
  941. TERM = 1.0
  942. DO 74 M = L,LLL
  943. IF (K.EQ.M) GO TO 74
  944. AM = M - 1
  945. TERM = TERM * (10.*B - AM)/(AK - AM)
  946. 74 CONTINUE
  947. TERM = TERM * FIELD(K)
  948. 75 WFLD = WFLD + TERM
  949. RETURN
  950. END
  951. C
  952. C
  953. FUNCTION COSINT(X)
  954. C
  955. C CALCULATION OF THE COSINE INTEGRAL
  956. C
  957. DOUBLE PRECISION Y2,PROD,SM,PT,DK
  958. IF(X.LE.0.) GO TO 50
  959. X2 = X * X
  960. IF(X.GT.20.) GO TO 30
  961. Y2 = DBLE(X2)
  962. PROD = -Y2 * 0.5
  963. SM = PROD * 0.5
  964. DO 10 K = 2,50
  965. DK = 2 * K
  966. PROD = -PROD * Y2/(DK*(DK - 1.))
  967. SM = SM + PROD/DK
  968. PT = ABS(PROD * 1.D+10)
  969. IF(ABS(SM).GT.PT) GO TO 20
  970. 10 CONTINUE
  971. 20 SS = SNGL(SM)
  972. COSINT = SS + 0.5772156649 + LOGF(X)
  973. RETURN
  974. 30 FA = 1.
  975. FB = 1.
  976. P0 = 1.
  977. X2 = 1./X2
  978. DO 40 K = 1,10
  979. AK = 2 * K
  980. P0 = -P0 * AK * X2
  981. FA = FA + P0
  982. P0 = P0 * (AK + 1.)
  983. FB = FB + P0
  984. PA = ABS(P0 + 1.E+10)
  985. IF(PA.LE.FB) GO TO 45
  986. 40 CONTINUE
  987. 45 FX = FA/X
  988. GX = FB * X2
  989. COSINT = FX * SIN(X) - GX * COS(X)
  990. RETURN
  991. 50 WRITE (61,100) X
  992. 100 FORMAT (* X LESS OR EQUAL TO ZERO, X = *E17.9)
  993. RETURN
  994. END
  995. C
  996. C
  997. SUBROUTINE BSJY01(X, AJ0, Y0, AJ1, Y1)
  998. C
  999. C CALCULATION OF THE BESSEL FUNCTIONS J0, Y0, J1, AND Y1 FOR AN
  1000. C ARGUMENT X
  1001. C
  1002. DIMENSION A(7), B(7), C(7), D(7), E(7), F(7), G(7), H(7)
  1003. DATA (A(I), I = 1,7) / 0.00021, -0.0039444, 0.0444479,
  1004. & -0.3163866, 1.2656208, -2.2499997, 1.0/
  1005. DATA (B(I), I = 1,7) / -0.00024846, 0.00427916, -0.04261214,
  1006. & 0.25300117, -0.74350384, 0.60559366, 0.36746691/
  1007. DATA (C(I), I = 1,7) / 0.00014476, -0.00072805, 0.00137237,
  1008. & -0.00009512, -0.0055274, -0.00000077, 0.79788456/
  1009. DATA (D(I), I = 1,7) / 0.00013558, -0.00029333, -0.00054125,
  1010. & 0.00262573, -0.00003954, -0.04166397, -0.78539816/
  1011. DATA (E(I), I = 1,7) / 0.00001109, -0.00031761, 0.00443319,
  1012. & -0.03954289, 0.21093573, -0.56249985, 0.5/
  1013. DATA (F(I), I = 1,7) / 0.0027873, -0.0400976, 0.3123951,
  1014. & -1.3164827, 2.1682709, 0.2212091, -0.6366198
  1015. DATA (G(I), I = 1,7) / -0.00020033, 0.00113653, -0.00249511,
  1016. & 0.00017105, 0.01659667, 0.00000156, 0.79788456/
  1017. DATA (H(I), I = 1,7) / -0.00029166, 0.00079824, 0.00074348,
  1018. & -0.00637879, 0.0000565, 0.12499612, -2.35619449/
  1019. AX = ABSF(X)
  1020. IF (AX.GT.0.0) GO TO 10
  1021. AJ0 = 1.
  1022. Y0 = -1.E+30
  1023. AJ1 = 0.
  1024. Y1 = -1.E+30
  1025. RETURN
  1026. 10 IF (AX.GT.3.0) GO TO 50
  1027. XX = (AX/3.0) ** 2
  1028. AJ0 = A(1)
  1029. Y0 = B(1)
  1030. AJ1 = E(1)
  1031. Y1 = F(1)
  1032. DO 20 M = 2,7
  1033. AJ0 = AJ0 * XX + A(M)
  1034. Y0 = Y0 * XX + B(M)
  1035. AJ1 = AJ1 * XX + E(M)
  1036. 20 Y1 = Y1 * XX + F(M)
  1037. AJ1 = AJ1 * X
  1038. ALF = 0.6366197724 * LOGF(0.5 * AX)
  1039. Y0 = Y0 + ALF * AJ0
  1040. Y1 = Y1/X + ALF * AJ1
  1041. RETURN
  1042. 50 X3 = 3.0/AX
  1043. F0 = C(1)
  1044. TH0 = D(1)
  1045. F1 = G(1)
  1046. TH1 = H(1)
  1047. DO 60 M = 2,7
  1048. F0 = F0 * X3 + C(M)
  1049. TH0 = TH0 * X3 + D(M)
  1050. F1 = F1 * X3 + G(M)
  1051. 60 TH1 = TH1 * X3 + H(M)
  1052. TH0 = TH0 + AX
  1053. TH1 = TH1 + AX
  1054. XS = 1./SQRT(AX)
  1055. AJ0= XS * F0 * COSF(TH0)
  1056. Y0 = XS * F0 * SINF(TH0)
  1057. AJ1= XS * F1 * COSF(TH1)
  1058. Y1 = XS * F1 * SINF(TH1)
  1059. RETURN
  1060. END
  1061. C
  1062. C
  1063.  
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!