?? analytic.f90
字號:
! analytic.f90!! Analytische Dipol Feldberechnung!! Copyright (C) 2007 Paul Panserrieu, < peutetre@cs.tu-berlin.de >!! This program is free software: you can redistribute it and/or modify! it under the terms of the GNU General Public License as published by! the Free Software Foundation, either version 3 of the License.! ! last modified: 25-10-2007 04:47:54 PM CESTMODULE analyticUSE fdtd_gitterIMPLICIT NONECONTAINS! return 1 if dipol is in this direction else 0INTEGER FUNCTION direction(d, komponent) TYPE(dipol), INTENT(IN) :: d INTEGER, INTENT(IN) :: komponent IF (d%E(komponent) .NE. 0.0d0) THEN direction = 1 ELSE direction = 0 ENDIFEND FUNCTION direction! distanz dipol <-> messzelle, |r|DOUBLE PRECISION FUNCTION r(g, d, messzelle) TYPE(gitter), INTENT(IN) :: g TYPE(dipol), INTENT(IN) :: d INTEGER, DIMENSION(1:3), INTENT(IN) :: messzelle r = SQRT(( g%dx * (messzelle(1) - d%px - 0.5d0 * direction(d, 1))) ** 2.0 & +(g%dy * (messzelle(2) - d%py - 0.5d0 * direction(d, 2))) ** 2.0 & +(g%dz * (messzelle(3) - d%pz - 0.5d0 * direction(d, 3))) ** 2.0)END FUNCTION r! components of vector rDOUBLE PRECISION FUNCTION rk(g, d, messzelle, k) TYPE(gitter), INTENT(IN) :: g TYPE(dipol), INTENT(IN) :: d INTEGER, DIMENSION(1:3), INTENT(IN) :: messzelle INTEGER, INTENT(IN) :: k IF (k .EQ. 1) THEN rk = g%dx * (messzelle(1) - d%px - 0.5d0 * direction(d, 1)) ELSEIF (k .EQ. 2) THEN rk = g%dy * (messzelle(2) - d%py - 0.5d0 * direction(d, 2)) ELSE rk = g%dz * (messzelle(3) - d%pz - 0.5d0 * direction(d, 3)) ENDIF END FUNCTION rkSUBROUTINE stimulus(derivs, zelle, g, d, dipol_type, zt_E, zt_H) DOUBLE PRECISION, INTENT(INOUT), DIMENSION(1:2,1:3) :: derivs INTEGER, INTENT(IN), DIMENSION(1:3) :: zelle TYPE(gitter), INTENT(IN) :: g TYPE(dipol), INTENT(IN) :: d INTEGER, INTENT(IN) :: dipol_type DOUBLE PRECISION, INTENT(IN) :: zt_E, zt_H SELECT CASE (dipol_type) CASE(1) CALL cosinusoidal(derivs, zelle, zt_E, zt_H, d, g) CASE(2) CALL sinusoidal(derivs, zelle, zt_E, zt_H, d, g) CASE(3) CALL mix(derivs, zelle, zt_E, zt_H, d, g) CASE DEFAULT WRITE(*,*) "Falscher 'dipol_type' in stimulus()" END SELECTEND SUBROUTINE stimulusDOUBLE PRECISION FUNCTION cos_deriv_0 (d, zeit, distanz) TYPE(dipol), INTENT(IN) :: d DOUBLE PRECISION, INTENT(IN) :: zeit DOUBLE PRECISION, INTENT(IN) :: distanz cos_deriv_0 = cos(d%omega * (zeit - (distanz/C)) + d%phi) END FUNCTION cos_deriv_0DOUBLE PRECISION FUNCTION cos_deriv_1 (d, zeit, distanz) TYPE(dipol), INTENT(IN) :: d DOUBLE PRECISION, INTENT(IN) :: zeit DOUBLE PRECISION, INTENT(IN) :: distanz cos_deriv_1 = - d%omega * sin(d%omega * (zeit - (distanz/C)) + d%phi) END FUNCTION cos_deriv_1DOUBLE PRECISION FUNCTION cos_deriv_2 (d, zeit, distanz) TYPE(dipol), INTENT(IN) :: d DOUBLE PRECISION, INTENT(IN) :: zeit DOUBLE PRECISION, INTENT(IN) :: distanz cos_deriv_2 = - (d%omega ** 2.0) * cos(d%omega * (zeit - (distanz/C)) + d%phi) END FUNCTION cos_deriv_2SUBROUTINE cosinusoidal(derivs, zelle, zt_E, zt_H, d, g) DOUBLE PRECISION, INTENT(INOUT), DIMENSION(1:2,1:3) :: derivs INTEGER, INTENT(IN), DIMENSION(1:3) :: zelle DOUBLE PRECISION, INTENT(IN) :: zt_E, zt_H TYPE(dipol), INTENT(IN) :: d TYPE(gitter), INTENT(IN) :: g DOUBLE PRECISION :: ampl, fac IF(zt_E >= r(g, d, zelle)/C) THEN CALL find_ampl(d%E, ampl) fac = EPS * 3.0d0 *(g%dx * g%dx) * ampl derivs(1,1) = - fac * cos_deriv_0(d, zt_E, r(g, d, zelle)) derivs(1,2) = - fac * cos_deriv_1(d, zt_E, r(g, d, zelle)) derivs(1,3) = - fac * cos_deriv_2(d, zt_E, r(g, d, zelle)) derivs(2,1) = - fac * cos_deriv_0(d, zt_H, r(g, d, zelle)) derivs(2,2) = - fac * cos_deriv_1(d, zt_H, r(g, d, zelle)) derivs(2,3) = - fac * cos_deriv_2(d, zt_H, r(g, d, zelle)) ELSE derivs = 0.0d0 ENDIFEND SUBROUTINE cosinusoidalDOUBLE PRECISION FUNCTION sin_deriv_0 (d, zeit, distanz) TYPE(dipol), INTENT(IN) :: d DOUBLE PRECISION, INTENT(IN) :: zeit DOUBLE PRECISION, INTENT(IN) :: distanz sin_deriv_0 = sin(d%omega * (zeit - (distanz/C)) + d%phi) END FUNCTION sin_deriv_0DOUBLE PRECISION FUNCTION sin_deriv_1 (d, zeit, distanz) TYPE(dipol), INTENT(IN) :: d DOUBLE PRECISION, INTENT(IN) :: zeit DOUBLE PRECISION, INTENT(IN) :: distanz sin_deriv_1 = d%omega * cos(d%omega * (zeit - (distanz/C)) + d%phi) END FUNCTION sin_deriv_1DOUBLE PRECISION FUNCTION sin_deriv_2 (d, zeit, distanz) TYPE(dipol), INTENT(IN) :: d DOUBLE PRECISION, INTENT(IN) :: zeit DOUBLE PRECISION, INTENT(IN) :: distanz sin_deriv_2 = - d%omega * d%omega * sin(d%omega * (zeit - (distanz/C)) + d%phi) END FUNCTION sin_deriv_2SUBROUTINE sinusoidal(derivs, zelle, zt_E, zt_H, d, g) DOUBLE PRECISION, INTENT(INOUT), DIMENSION(1:2,1:3) :: derivs INTEGER, INTENT(IN), DIMENSION(1:3) :: zelle DOUBLE PRECISION, INTENT(IN) :: zt_E, zt_H TYPE(dipol), INTENT(IN) :: d TYPE(gitter), INTENT(IN) :: g DOUBLE PRECISION :: ampl, fac IF(zt_E >= r(g, d, zelle)/C) THEN CALL find_ampl(d%E, ampl) fac = EPS * 3.0d0 * (g%dx * g%dx) * ampl derivs(1,1) = - fac * sin_deriv_0(d, zt_E, r(g, d, zelle)) derivs(1,2) = - fac * sin_deriv_1(d, zt_E, r(g, d, zelle)) derivs(1,3) = - fac * sin_deriv_2(d, zt_E, r(g, d, zelle)) derivs(2,1) = - fac * sin_deriv_0(d, zt_H, r(g, d, zelle)) derivs(2,2) = - fac * sin_deriv_1(d, zt_H, r(g, d, zelle)) derivs(2,3) = - fac * sin_deriv_2(d, zt_H, r(g, d, zelle)) ELSE derivs = 0.0d0 ENDIFEND SUBROUTINE sinusoidalDOUBLE PRECISION FUNCTION mix_deriv_0 (d, zeit, distanz) TYPE(dipol), INTENT(IN) :: d DOUBLE PRECISION, INTENT(IN) :: zeit DOUBLE PRECISION, INTENT(IN) :: distanz mix_deriv_0 = SIN(d%omega * (zeit - (distanz/C)) + d%phi) & * (1.0d0 - COS(d%omega * (zeit - (distanz/C)) + d%phi)) END FUNCTION mix_deriv_0DOUBLE PRECISION FUNCTION mix_deriv_1 (d, zeit, distanz) TYPE(dipol), INTENT(IN) :: d DOUBLE PRECISION, INTENT(IN) :: zeit DOUBLE PRECISION, INTENT(IN) :: distanz mix_deriv_1 = d%omega * (COS(d%omega * (zeit - (distanz/C)) + d%phi) & + SIN(d%omega * (zeit - (distanz/C)) + d%phi) ** 2.0 & - COS(d%omega * (zeit - (distanz/C)) + d%phi) ** 2.0) END FUNCTION mix_deriv_1DOUBLE PRECISION FUNCTION mix_deriv_2 (d, zeit, distanz) TYPE(dipol), INTENT(IN) :: d DOUBLE PRECISION, INTENT(IN) :: zeit DOUBLE PRECISION, INTENT(IN) :: distanz mix_deriv_2 = (d%omega ** 2.0) * SIN(d%omega * (zeit - (distanz/C)) + d%phi) & * (4.0d0 * COS(d%omega * (zeit - (distanz/C)) + d%phi) - 1.0d0) END FUNCTION mix_deriv_2SUBROUTINE mix(derivs, zelle, zt_E, zt_H, d, g) DOUBLE PRECISION, INTENT(INOUT), DIMENSION(1:2,1:3) :: derivs INTEGER, INTENT(IN), DIMENSION(1:3) :: zelle DOUBLE PRECISION, INTENT(IN) :: zt_E, zt_H TYPE(dipol), INTENT(IN) :: d TYPE(gitter), INTENT(IN) :: g DOUBLE PRECISION :: ampl, fac IF(zt_E >= r(g, d, zelle)/C) THEN CALL find_ampl(d%E, ampl) fac = ampl * EPS * 3.0d0 *(g%dx * g%dx) derivs(1,1) = - fac * mix_deriv_0(d, zt_E, r(g, d, zelle)) derivs(1,2) = - fac * mix_deriv_1(d, zt_E, r(g, d, zelle)) derivs(1,3) = - fac * mix_deriv_2(d, zt_E, r(g, d, zelle)) derivs(2,1) = - fac * mix_deriv_0(d, zt_H, r(g, d, zelle)) derivs(2,2) = - fac * mix_deriv_1(d, zt_H, r(g, d, zelle)) derivs(2,3) = - fac * mix_deriv_2(d, zt_H, r(g, d, zelle)) ELSE derivs = 0.0d0 ENDIFEND SUBROUTINE mixSUBROUTINE feldberechnung(f, s, zelle, zt_E, d, g) DOUBLE PRECISION, INTENT(INOUT), DIMENSION(1:6) :: f ! Feld array DOUBLE PRECISION, INTENT(IN), DIMENSION(1:2,1:3) :: s ! Stimulus INTEGER, INTENT(IN), DIMENSION(1:3) :: zelle ! zelle DOUBLE PRECISION, INTENT(IN) :: zt_E ! Zeit_E TYPE(dipol), INTENT(IN) :: d ! Dipolstruktur TYPE(gitter), INTENT(IN) :: g INTEGER :: pos1, pos2, pos3, pos, i pos1 = 1; pos2 = 2; pos3 = 3 IF (zt_E >= r(g, d, zelle)/C) THEN pos = direction(d, 1) + 2 * direction(d, 2) + 3 * direction(d, 3) IF (pos .EQ. pos1) THEN CALL permut(pos1, 1) CALL permut(pos2, 1) CALL permut(pos3, 1) ELSEIF(pos .EQ. pos2) THEN CALL permut(pos1, -1) CALL permut(pos2, -1) CALL permut(pos3, -1) ENDIF ! Ex, Ey, Ez f(pos1) = g%dx * rk(g, d, zelle, 1) * rk(g, d, zelle, 3) & / (4.0d0 * PI * EPS * r(g, d, zelle) ** 5.0) & * ( 3.0d0 * (s(1,1) + (r(g, d, zelle)/C) * s(1,2)) & + (r(g, d, zelle)/C)**2.0 * s(1,3) ) f(pos2) = g%dx * rk(g, d, zelle, 2) * rk(g, d, zelle, 3) & / (4.0d0 * PI * EPS * r(g, d, zelle) ** 5.0) & * ( 3.0d0 * (s(1,1) + (r(g, d, zelle)/C) * s(1,2)) & + (r(g, d, zelle)/C)**2.0 * s(1,3) ) f(pos3) = g%dx / (4.0d0 * PI * EPS * r(g, d, zelle) ** 5) & * ( ( 2.0d0 * rk(g, d, zelle, 3) ** 2.0 & - (rk(g, d, zelle, 1)**2.0 + rk(g, d, zelle, 2) ** 2.0) & ) * (s(1,1) + (r(g, d, zelle)/C) * s(1,2)) & - (rk(g, d, zelle, 1)**2.0 + rk(g, d, zelle, 2) ** 2.0) & * (r(g, d, zelle)/C) ** 2.0 * s(1,3) ) ! Hx, Hy, Hz f(pos1+3) = - g%dx * rk(g, d, zelle, 2) & / (4.0d0 * PI * r(g, d, zelle) ** 3.0) & * (s(2,2) + (r(g, d, zelle)/C) * s(2,3)) f(pos2+3) = g%dx * rk(g, d, zelle, 1) & / (4.0d0 * PI * r(g, d, zelle) ** 3.0) & * (s(2,2) + (r(g, d, zelle)/C) * s(2,3)) f(pos3+3) = 0.0d0 ELSE f = 0.0d0 ENDIFEND SUBROUTINE feldberechnungSUBROUTINE permut(i, sinn) INTEGER, INTENT(INOUT) :: i INTEGER, INTENT(IN) :: sinn INTEGER, DIMENSION(1:3) :: triplet triplet = (/1, 2, 3/); triplet = CSHIFT(triplet, SIGN(1, sinn)) i = triplet(i)END SUBROUTINE permutSUBROUTINE find_ampl(array, ampl) DOUBLE PRECISION, INTENT(IN), DIMENSION(1:3) :: array DOUBLE PRECISION, INTENT(INOUT) :: ampl INTEGER :: i DO i = 1, 3, 1 IF (array(i) .NE. 0.0d0) THEN ampl = array(i) ENDIF ENDDOEND SUBROUTINE find_amplEND MODULE analytic ?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -