API tools faq
paste
Advertisement
Guest User

Untitled

a guest
May 23rd, 2019
91
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 15.81 KB | None | 0 0
raw download clone embed print report
  1. program Helmholtz
  2. implicit none
  3.  
  4. real(kind=8), parameter :: ystar(2) = (0,0), kappa=2, eta=kappa, eulerconst = 0.5772156649, pi = 4.D0*DATAN(1.D0)
  5. integer :: n = 1, k, j
  6. real(kind=8) :: testpoint(2) = (3,3)
  7. real(kind=8) , allocatable :: matrix(:,:), rightpart (:), density(:)
  8.  
  9. do while (n<=200)
  10. ALLOCATE (matrix(4*n,4*n), rightpart(4*n), density(4*n))
  11. do k=0, 2*n-1
  12. do j=0, 2*n-1
  13. Matrix(k + 1, j+1) = Lre(k, j);
  14. Matrix(k + 1, j + 2 * n+1) = -Lim(k, j);
  15. Matrix(k + 2 * n+1, j+1) = Lim(k, j);
  16. Matrix(k + 2 * n+1, j + 2 * n+1) = Lre(k, j);
  17. if (k == j) then
  18. Matrix(k+1, j + 2 * n+1) = Matrix(k+1, j + 2 * n+1) - aim(t(k));
  19. Matrix(k + 2 * n+1, j+1) = Matrix(k + 2 * n+1, j+1) + aim(t(k));
  20. endif
  21. rightPart(k+1) = fre(t(k));
  22. rightPart(k + 2 * n+1) = fim(t(k));
  23. end do
  24. end do
  25. density = solve_wbs(ge_wpp(Matrix,rightPart))
  26. !write(*,*) exactIm(testPoint)
  27. write(*,*) "n = ", n, " ", Abs(resultRe(testPoint, density) - exactRe(testPoint)), " ", &
  28. Abs(resultIm(testPoint, density) - exactIm(testPoint))
  29. DEALLOCATE (matrix, rightpart, density)
  30. n = n*2
  31. end do
  32.  
  33. contains
  34.  
  35. real(kind=8) function x1(t)
  36. real(kind=8) t
  37. x1 = 2.*cos(t)
  38. return
  39. end
  40.  
  41. real(kind=8) function x2(t)
  42. real(kind=8) t
  43. x2 = 2.*sin(t)
  44. return
  45. end
  46.  
  47. real(kind=8) function x1d(t)
  48. real(kind=8) t
  49. x1d = -2.*sin(t)
  50. return
  51. end
  52.  
  53. real(kind=8) function x2d(t)
  54. real(kind=8) t
  55. x2d = 2.*cos(t)
  56. return
  57. end
  58.  
  59. real(kind=8) function x1dd(t)
  60. real(kind=8) t
  61. x1dd = -2.*cos(t)
  62. return
  63. end
  64.  
  65. real(kind=8) function x2dd(t)
  66. real(kind=8) t
  67. x2dd = -2.*sin(t)
  68. return
  69. end
  70.  
  71. real(kind=8) function x1ddd(t)
  72. real(kind=8) t
  73. x1ddd = 2.*sin(t)
  74. return
  75. end
  76.  
  77. real(kind=8) function x2ddd(t)
  78. real(kind=8) t
  79. x2ddd = -2.*cos(t)
  80. return
  81. end
  82.  
  83. real(kind=8) function diffmodule(t, tau)
  84. real(kind=8) t, tau
  85. diffmodule = SQRT((x1(t) - x1(tau)) * (x1(t) - x1(tau)) + (x2(t) - x2(tau)) * (x2(t) - x2(tau)))
  86. return
  87. end
  88.  
  89. real(kind=8) function dermodule(t)
  90. real(kind=8) t
  91. dermodule = SQRT(x1d(t) * x1d(t) + x2d(t) * x2d(t))
  92. return
  93. end
  94.  
  95. real(kind=8) function doubledermodule(t)
  96. real(kind=8) t
  97. doubledermodule = SQRT(x1dd(t) * x1dd(t) + x2dd(t) * x2dd(t))
  98. return
  99. end
  100.  
  101. real(kind=8) function h(t, tau)
  102. real(kind=8) t, tau, numerator, denominator
  103. numerator = ((x1(tau) - x1(t)) * x2d(t) - (x2(tau) - x2(t)) * x1d(t))
  104. denominator = diffModule(t, tau)
  105. h = numerator / denominator
  106. return
  107. end
  108.  
  109. real(kind=8) function him(t, tau)
  110. real(kind=8) t, tau
  111. him = kappa * h(t, tau) * BESSEL_J1(kappa * diffModule(t, tau)) * derModule(tau) / 2.0
  112. return
  113. end
  114.  
  115. real(kind=8) function hre(t, tau)
  116. real(kind=8) t, tau
  117. hre = -kappa * h(t, tau) * BESSEL_Y1(kappa * diffModule(t, tau)) * derModule(tau) / 2.0;
  118. return
  119. end
  120.  
  121. real(kind=8) function h1(t, tau)
  122. real(kind=8) t, tau
  123. if (t /= tau) then
  124. h1 = -kappa * h(t, tau) * BESSEL_J1(kappa * diffModule(t, tau)) * derModule(tau) / (2.0 * pi)
  125. return
  126. else
  127. h1 = 0.0
  128. return
  129. end if
  130. end
  131.  
  132. real(kind=8) function h2re(t, tau)
  133. real(kind=8) t, tau
  134. if (t /= tau) then
  135. h2re = hre(t, tau) - h1(t, tau) * log(4.0 * sin((t - tau) / 2.0) * sin((t - tau) / 2.0))
  136. return
  137. else
  138. h2re = (x2d(t) * x1dd(t) - x1d(t) * x2dd(t)) / (2.0 * pi * derModule(t))
  139. return
  140. end if
  141. end
  142.  
  143. real(kind=8) function h2im(t, tau)
  144. real(kind=8) t, tau
  145. if (t /= tau) then
  146. h2im = him(t, tau)
  147. return
  148. else
  149. h2im = 0.0
  150. return
  151. end if
  152. end
  153.  
  154. real(kind=8) function mim(t, tau)
  155. real(kind=8) t, tau
  156. mim = BESSEL_J0(kappa * diffModule(t, tau)) / 2.0
  157. return
  158. end
  159.  
  160. real(kind=8) function mre(t, tau)
  161. real(kind=8) t, tau
  162. mre = -BESSEL_Y0(kappa * diffModule(t, tau)) / 2.0
  163. return
  164. end
  165.  
  166. real(kind=8) function m1(t, tau)
  167. real(kind=8) t, tau
  168. if (t /= tau) then
  169. m1 = -BESSEL_J0(kappa * diffModule(t, tau)) / (2.0 * pi)
  170. return
  171. else
  172. m1 = -1 / (2.0 * pi)
  173. return
  174. end if
  175. end
  176.  
  177. real(kind=8) function m2re(t, tau)
  178. real(kind=8) t, tau
  179. if (t /= tau) then
  180. m2re = mre(t, tau) - m1(t, tau) * log(4.0 * sin((t - tau) / 2.0) * sin((t - tau) / 2.0))
  181. return
  182. else
  183. m2re = -(eulerconst + log(kappa * derModule(t) / 2.0)) / pi
  184. return
  185. end if
  186. end
  187.  
  188. real(kind=8) function m2im(t, tau)
  189. real(kind=8) t, tau
  190. if (t /= tau) then
  191. m2im = mim(t, tau)
  192. return
  193. else
  194. m2im = 1.0 / 2.0
  195. return
  196. end if
  197. end
  198.  
  199. real(kind=8) function tildan(t, tau)
  200. real(kind=8) t, tau, numerator, denominator
  201. numerator = ((x1(t) - x1(tau)) * x1d(t) + (x2(t) - x2(tau)) * x2d(t)) * &
  202. ((x1(t) - x1(tau)) * x1d(tau) + (x2(t) - x2(tau)) * x2d(tau))
  203. denominator = diffModule(t, tau) * diffModule(t, tau)
  204. tildan = numerator / denominator
  205. return
  206. end
  207.  
  208. real(kind=8) function nre(t, tau)
  209. real(kind=8) t, tau
  210. nre = 1.0 / (4.0 * pi * sin((t - tau) / 2.0) * sin((t - tau) / 2.0)) + tildan(t, tau) * &
  211. (-kappa * kappa * BESSEL_Y0(kappa * diffModule(t, tau)) + &
  212. 2.0 * kappa * BESSEL_Y1(kappa * diffModule(t, tau)) / diffModule(t, tau)) / 2.0 &
  213. - kappa * (x1d(t) * x1d(tau) + x2d(t) * x2d(tau)) * BESSEL_Y1(kappa * diffModule(t, tau)) / (2.0 * diffModule(t, tau))
  214. return
  215. end
  216.  
  217. real(kind=8) function nim(t, tau)
  218. real(kind=8) t, tau
  219. nim = tildan(t, tau) * (kappa * kappa * BESSEL_J0(kappa * diffModule(t, tau)) - &
  220. 2 * kappa * BESSEL_J1(kappa * diffModule(t, tau)) / diffModule(t, tau)) / 2.0 + &
  221. kappa * (x1d(t) * x1d(tau) + x2d(t) * x2d(tau)) * BESSEL_J1(kappa * diffModule(t, tau)) / (2.0 * diffModule(t, tau))
  222. return
  223. end
  224.  
  225. real(kind=8) function n1(t, tau)
  226. real(kind=8) t, tau
  227. if (t /= tau) then
  228. n1 = -tildaN(t, tau) * (kappa * kappa * BESSEL_J0(kappa * diffModule(t, tau)) - 2 * kappa * BESSEL_J1(kappa * &
  229. diffModule(t, tau)) / diffModule(t, tau)) / (2.0 * pi) - &
  230. kappa * BESSEL_J1(kappa * diffModule(t, tau)) * (x1d(t) * x1d(tau) + x2d(t) * x2d(tau)) / &
  231. (2.0 * pi * diffModule(t, tau))
  232. return
  233. else
  234. n1 = -kappa * kappa * derModule(t) * derModule(t) / (4.0 * pi)
  235. return
  236. end if
  237. end
  238.  
  239. real(kind=8) function n2re(t, tau)
  240. real(kind=8) t, tau
  241. if (t /= tau) then
  242. n2re = Nre(t, tau) - N1(t, tau) * log(4.0 * sin((t - tau) / 2.0) * sin((t - tau) / 2.0))
  243. return
  244. else
  245. n2re = -kappa * kappa * derModule(t) * derModule(t) / (4.0 * pi) - eulerconst * kappa * kappa * derModule(t) &
  246. * derModule(t) / (2.0 * pi) - kappa * kappa * derModule(t) * derModule(t) * &
  247. log(kappa * derModule(t) / 2.0) / (2.0 * pi) + 1 / (12.0 * pi) &
  248. + (x1d(t) * x1dd(t) + x2d(t) * x2dd(t)) * (x1d(t) * x1dd(t) + x2d(t) * x2dd(t)) / &
  249. (2.0 * pi * derModule(t)*derModule(t)*derModule(t)*derModule(t)) &
  250. - doubleDerModule(t) * doubleDerModule(t) / (4.0 * pi * derModule(t) * derModule(t)) &
  251. - (x1d(t) * x1ddd(t) + x2d(t) * x2ddd(t)) / (6.0 * pi * derModule(t) * derModule(t))
  252. return
  253. end if
  254. end
  255.  
  256. real(kind=8) function n2im(t, tau)
  257. real(kind=8) t, tau
  258. if (t /= tau) then
  259. n2im = Nim(t, tau)
  260. return
  261. else
  262. n2im = kappa * kappa * derModule(t) * derModule(t) / 4.0
  263. return
  264. end if
  265. end
  266.  
  267. real(kind=8) function k1re(t, tau)
  268. real(kind=8) t, tau
  269. k1re = kappa * kappa * (x1d(t) * x1d(tau) + x2d(t) * x2d(tau)) * M1(t, tau) - N1(t, tau)
  270. return
  271. end
  272.  
  273. real(kind=8) function k1im(t, tau)
  274. real(kind=8) t, tau
  275. k1im = -eta * H1(t, tau)
  276. return
  277. end
  278.  
  279. real(kind=8) function k2re(t, tau)
  280. real(kind=8) t, tau
  281. k2re = kappa * kappa * (x1d(t) * x1d(tau) + x2d(t) * x2d(tau)) * M2re(t, tau) - N2re(t, tau) + eta * H2im(t, tau)
  282. return
  283. end
  284.  
  285. real(kind=8) function k2im(t, tau)
  286. real(kind=8) t, tau
  287. k2im = kappa * kappa * (x1d(t) * x1d(tau) + x2d(t) * x2d(tau)) * M2im(t, tau) - N2im(t, tau) - eta * H2re(t, tau)
  288. return
  289. end
  290.  
  291. real(kind=8) function tcapital(j)
  292. integer j
  293. if (j == 0) then
  294. tcapital = -n / 2.0
  295. return
  296. else if (mod(j,2) == 0) then
  297. tcapital = 0.0
  298. return
  299. else
  300. tcapital = 1 / (2.0 * n * sin(t(j) / 2.0) * sin(t(j) / 2.0))
  301. return
  302. end if
  303. end
  304.  
  305. real(kind=8) function t(j)
  306. integer j
  307. t = j * pi / n
  308. end
  309.  
  310. real(kind=8) function r(j)
  311. integer j, m
  312. real(kind=8) sum
  313. sum = 0.0
  314. do m = 1, n-1
  315. sum = sum + cos(m * j * pi / n) / m
  316. end do
  317. r = -2 * pi * sum / n + (-1)**j * pi / (n * n)
  318. end
  319.  
  320. real(kind=8) function aim(t)
  321. real(kind=8) t
  322. aim = eta * derModule(t)
  323. return
  324. end
  325.  
  326. real(kind=8) function lre(k, j)
  327. integer k, j
  328. lre = tcapital(abs(k - j)) + r(abs(k - j)) * K1re(t(k), t(j)) + pi * K2re(t(k), t(j)) / n
  329. return
  330. end
  331.  
  332. real(kind=8) function lim(k, j)
  333. integer k, j
  334. lim = r(abs(k - j)) * K1im(t(k), t(j)) + pi * K2im(t(k), t(j)) / n
  335. return
  336. end
  337.  
  338. real(kind=8) function uinfre(xhat, density)
  339. real(kind=8), dimension(:) :: xhat, density
  340. real(kind=8) :: sum = 0.0
  341. integer j
  342. do j = 0, 2 * n - 1
  343. sum = sum + (kappa * (xHat(1) * x2d(t(j)) - xHat(2) * x1d(t(j))) + eta) * derModule(t(j)) * &
  344. Ere(xHat, j, density)/(derModule(t(j))* derModule(t(j)))
  345. end do
  346. uinfre = pi * sum / (n * SQRT(8 * pi * kappa))
  347. return
  348. end
  349.  
  350. real(kind=8) function ere(xhat, j, density)
  351. real(kind=8), dimension(:) :: xhat, density
  352. real(kind=8) scalarprod
  353. integer j
  354. scalarprod = (xHat(1) * x1(t(j)) + xHat(2) * x2(t(j))) / SQRT(xHat(1) * xHat(1) + xHat(2) * xHat(2))
  355. ere = cos(kappa * scalarprod + pi / 4.0) * density(j+1) + sin(kappa * scalarprod + pi / 4.0) * density(j + 2 * n+1)
  356. return
  357. end
  358.  
  359. real(kind=8) function uinfim(xhat, density)
  360. real(kind=8), dimension(:) :: xhat, density
  361. real(kind=8) :: sum = 0.0
  362. integer j
  363. do j = 0, 2 * n - 1
  364. sum = sum + (kappa * (xHat(1) * x2d(t(j)) - xHat(2) * x1d(t(j))) + eta) * (derModule(t(j))) * Eim(xHat, j, density)
  365. end do
  366. uinfim = pi * sum / (n * SQRT(8 * pi * kappa))
  367. return
  368. end
  369.  
  370. real(kind=8) function eim(xhat, j, density)
  371. real(kind=8), dimension(:) :: xhat, density
  372. real(kind=8) scalarprod
  373. integer j
  374. scalarprod = (xHat(1) * x1(t(j)) + xHat(2) * x2(t(j))) / SQRT(xHat(1) * xHat(1) + xHat(2) * xHat(2))
  375. eim = -sin(kappa * scalarprod + pi / 4.0) * density(j+1) + cos(kappa * scalarprod + pi / 4.0) * density(j + 2 * n + 1)
  376. return
  377. end
  378.  
  379. real(kind=8) function fre(t)
  380. real(kind=8) t
  381. real(kind=8), dimension(2) :: x
  382. x(1) = x1(t)
  383. x(2) = x2(t)
  384. fre = kappa * BESSEL_Y1(kappa * module(x, yStar)) * ((x1(t) - yStar(1)) * x2d(t) - (x2(t) - yStar(2)) * x1d(t)) / &
  385. (2 * module(x, yStar));
  386. end
  387.  
  388. real(kind=8) function fim(t)
  389. real(kind=8) t
  390. real(kind=8), dimension(2) :: x
  391. x(1) = x1(t)
  392. x(2) = x2(t)
  393. fim = -kappa * BESSEL_J1(kappa * module(x, yStar)) * ((x1(t) - yStar(1)) * x2d(t) - (x2(t) - yStar(2)) * x1d(t)) / &
  394. (2 * module(x, yStar));
  395. end
  396.  
  397. real(kind=8) function module(x, y)
  398. real(kind=8), dimension(:) :: x,y
  399. module = Sqrt((x(1) - y(1)) * (x(1) - y(1)) + (x(2) - y(2)) * (x(2) - y(2)))
  400. return
  401. end
  402.  
  403. function solve_wbs(u) result(x) ! solve with backward substitution
  404. real(kind=8) :: u(:,:)
  405. integer :: i,n
  406. real(kind=8) , allocatable :: x(:)
  407. n = size(u,1)
  408. allocate(x(n))
  409. forall (i=n:1:-1) x(i) = ( u(i,n+1) - sum(u(i,i+1:n)*x(i+1:n)) ) / u(i,i)
  410. end function
  411.  
  412. function ge_wpp(a,b) result(u) ! gaussian eliminate with partial pivoting
  413. real(kind=8) :: a(:,:),b(:),upi
  414. integer :: i,j,n,p
  415. real(kind=8) , allocatable :: u(:,:)
  416. n = size(a,1)
  417. u = reshape( [a,b], [n,n+1] )
  418. do j=1,n
  419. p = maxloc(abs(u(j:n,j)),1) + j-1 ! maxloc returns indices between (1,n-j+1)
  420. if (p /= j) u([p,j],j) = u([j,p],j)
  421. u(j+1:,j) = u(j+1:,j)/u(j,j)
  422. do i=j+1,n+1
  423. upi = u(p,i)
  424. if (p /= j) u([p,j],i) = u([j,p],i)
  425. u(j+1:n,i) = u(j+1:n,i) - upi*u(j+1:n,j)
  426. end do
  427. end do
  428. end function
  429.  
  430. real(kind=8) function resultRe(x, density)
  431. real(kind=8), dimension(:) :: x,density
  432. real(kind=8), dimension(2) :: temp
  433. real(kind=8) sum ,h
  434. integer j
  435. sum = 0
  436. do j=0,2*n-1
  437. temp(1) = x1(t(j))
  438. temp(2) = x2(t(j))
  439. h = ((x1(t(j)) - x(1)) * x2d(t(j)) - (x2(t(j)) - x(2)) * x1d(t(j))) / (module(x, temp) * derModule(t(j)));
  440. sum = sum + (density(j+1) * (kappa * BESSEL_Y1(kappa * module(x, temp)) * h + eta * BESSEL_J0(kappa * &
  441. module(x, temp))) - density(j + 2 * n+1) * (eta * BESSEL_Y0(kappa * module(x, temp)) - &
  442. kappa * BESSEL_J1(kappa * module(x, temp)) * h)) * derModule(t(j))
  443. end do
  444. resultRe = pi * sum / (4 * n)
  445. return
  446. end
  447.  
  448. real(kind=8) function exactRe(x)
  449. real(kind=8), dimension(:) :: x
  450. exactRe =-BESSEL_Y0(kappa*module(x, yStar))/4
  451. return
  452. end
  453.  
  454. real(kind=8) function resultIm(x, density)
  455. real(kind=8), dimension(:) :: x,density
  456. real(kind=8), dimension(2) :: temp
  457. real(kind=8) sum ,h
  458. integer j
  459. sum = 0
  460. do j=0,2*n-1
  461. temp(1) = x1(t(j))
  462. temp(2) = x2(t(j))
  463. h = ((x1(t(j)) - x(1)) * x2d(t(j)) - (x2(t(j)) - x(2)) * x1d(t(j))) / (module(x, temp) * derModule(t(j)));
  464. sum = sum + (density(j+1) * (eta * BESSEL_Y0(kappa * module(x, temp)) - &
  465. kappa * BESSEL_J1(kappa * module(x, temp)) * h) &
  466. + density(j + 2*n + 1) * (kappa * BESSEL_Y1(kappa * module(x, temp)) * h + &
  467. eta * BESSEL_J0(k * module(x, temp)))) * derModule(t(j))
  468. end do
  469. resultIm = pi * sum / (4 * n)
  470. return
  471. end
  472.  
  473. real(kind=8) function exactIm(x)
  474. real(kind=8), dimension(:) :: x
  475. exactIm = BESSEL_J0(kappa*module(x, yStar))/4
  476. return
  477. end
  478.  
  479. end program Helmholtz
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
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!