REPÚBLICA BOLIVARIANA DE VENEZUELA UNIVERSIDAD DEL ZULIA
FACULTAD DE INGENIERÍA DIVISION DE POSTGRADO
PROGRAMA DE POSTGRADO EN MATEMÁTICA APLICADA
APLICACIÓN DEL MÉTODO DE DIFERENCIAS FINITAS EN EL DOMINIO DEL TIEMPO PARA LA SIMULACIÓN DE LA PROPAGACIÓN
DE ONDAS ELECTROMAGNÉTICAS EN AMBIENTES CERRADOS
Trabajo de Grado presentado ante la Ilustre Universidad del Zulia
para optar al Grado Académico de
MAGISTER SCIENTIARIUM EN MATEMÁTICA APLICADA
Autor: Ing. Maryory Urdaneta Tutor: Prof. Ángel Ochoa
Co-Tutor: Prof. José Canelón
Maracaibo, Octubre de 2007
APROBACION Este Jurado aprueba el Trabajo de Grado titulado: APLICACIÓN DEL MÉTODO DE DIFERENCIAS FINITAS EN EL DOMINIO DEL TIEMPO PARA LA SIMULACIÓN DE LA PROPAGACIÓN DE ONDAS ELECTROMAGNÉTICAS EN AMBIENTES CERRADOS, que Maryory Urdaneta, C.I: 14.374.608 presenta ante el Consejo técnico de la División de Postgrado de la Facultad de Ingeniería en cumplimiento del Artículo 51, parágrafo 51.6 de la Sección Segunda del Reglamento de Estudios para Graduados de la Universidad del Zulia, como requisito para optar al Grado Académico de:
MAGISTER SCIENTIARUM EN MATEMÁTICA APLICADA
Coordinador del Jurado Nombre y apellidos
C.I. :_________________ Nombre y apellidos Nombre y apellidos C.I. :________________ C.I. :_________________
Director de la División de Postgrado
Nombre y apellido
Maracaibo, Octubre de 2007
Urdaneta, Maryory. Aplicación del método de diferencias finitas en el dominio del tiempo para la simulación de la propagación de ondas electromagnéticas en ambientes cerrados. (2007). Trabajo de Grado. Universidad del Zulia. Facultad de Ingeniería. División de Postgrado. Maracaibo, Tutor: Prof. Ángel Ochoa; Co-Tutor: Prof. José Canelón.
RESUMEN En este trabajo se utiliza el método de diferencias finitas en el dominio del tiempo (FDTD), para simular la propagación bidimensional de ondas electromagnéticas en un espacio cerrado. La aplicación de este método inicia con la formulación discreta de las ecuaciones de Maxwell. Posteriormente se diseñó un algoritmo para la simulación, con condiciones de frontera establecidas. Para la verificación del algoritmo desarrollado se realizó la simulación de un caso de estudio previamente establecido en trabajos anteriores. En este caso de estudio se simuló un piso completo de un edificio con la fuente generadora de ondas electromagnéticas colocada en su interior. Se obtuvo resultados para los campos eléctrico y magnético en todos los puntos del área de estudio y se comparó estos resultados, obtenidos del algoritmo FDTD, con otras técnicas utilizadas para la propagación en ambientes interiores y con valores reales medidos en las áreas de estudio. El software desarrollado en lenguaje Matlab es adaptable a cualquier espacio que se desee estudiar. Palabras clave: Diferencias finitas, propagación en interiores, Matlab E-mail del autor: [email protected] .
Urdaneta, Maryory. Aplicación del método de diferencias finitas en el dominio del tiempo para la simulación de la propagación de ondas electromagnéticas en ambientes cerrados. (2007). Trabajo Especial de Grado. Universidad del Zulia. Facultad de Ingeniería. División de Postgrado. Maracaibo, Tutor: Prof. Ángel Ochoa; Co-Tutor: Prof. José Canelón.
ABSTRACT In this work the finite differences time domain method (FDTD) is used in order to simulate the two dimensions electromagnetic waves propagation in a closed space. The application of this method starts with the discrete formulation of Maxwell’s equations. Lately, it was designed an algorithm for the simulation, with established boundary conditions. The developed algorithm was verified trough the simulation of a case of study established in previous works. In this case of study it was simulated a complete floor of a building with the source placed inside of it. In was obtained results for the electric and magnetic fields, in all the points of the study area, were obtained and compared to those resulting from the application of other techniques used for the inner atmosphere propagation and with measured real values in the same study areas. The software developed in Matlab language is adaptable to any space that is desired to study. Key Words: Finite differences, indoor propagation, Matlab. Author’s e-mail: [email protected] .
AGRADECIMIENTOS
A Dios por darme la vida.
A mi familia por apoyarme siempre.
A la Universidad del Zulia por permitirme ser parte de su comunidad.
A los Profesores Ángel Ochoa, Manuel Briceño, y José Canelón, por su
dedicación y apoyo incondicional.
A mis compañeros de maestría por ser un apoyo durante mis estudios.
A Dios por darme la vida.
A mi familia por apoyarme siempre.
A la Universidad Rafael Belloso Chacín por permitirme ser parte de su
comunidad.
A mis compañeros de maestría por ser un apoyo durante mis estudios.
A mi tutor Pedro González, por su comprensión y dedicación
TABLA DE CONTENIDO
Página RESUMEN...………………….……………………………..………………………...… 3
ABSTRACT …………….……………………………………………………………….. 4
DEDICATORIA………….………………………………………………………………. 5
AGRADECIMIENTO….………………………………………………………………… 6
TABLA DE CONTENIDO….….……………………………..…………...................... 7
LISTA DE FIGURAS.…………………………………………..………………………. 9
INTRODUCCIÓN…….……………………………………………………….………… 10
CAPITULO I. EL PROBLEMA………………………………………………………… 11
1.1. Planteamiento del Problema……………………………...…............................. 11
1.2. Objetivos de la Investigación……………………………...……………………… 12
1.2.1. Objetivo General…………………………………….…………………………... 12
1.2.2. Objetivos Específicos………………………………….................................... 13
1.3. Justificación de la Investigación…..……………………….…………………….. 13
CAPITULO II. MARCO TEÓRICO……………………………..……...……………… 14
2.1. Método de diferencias finitas……………………………………….…………….. 14
2.1.1. Esquema de las diferencias finitas…………………………………………….. 14
2.2. Ecuaciones de Maxwell………………………………………………………....... 17
2.2.1. Modos de propagación………………………………………………………….. 28
2.3. Formulación de las ecuaciones de Maxwell en FDTD……………….………... 30
2.4. Estabilidad Numérica…………………………………………………………....... 37
2.5. Técnica Ray Tracing………………………………………………………………. 38
CAPITULO III. MARCO METODOLÓGICO……….………………......................... 40
3.1. Metodología Aplicada……………………………..………….............................. 40
3.1.1. Definición de constantes………………………………………………………... 41
3.1.2. Particiones en espacio y tiempo……………………………………………….. 41
3.1.3. Definición del dieléctrico………………………………………………………... 41
3.1.4. Formulación de las ecuaciones de campo eléctrico y magnético en dos
dimensiones……………………………………………………………………………...
41
3.1.5. Condiciones de frontera………………………………………………………… 41
CAPITULO IV. DESARROLLO E IMPLEMENTACIÓN DEL SOFTWARE……... 42
4.1.Desarrollo del software.……………………………………………………………. 42
4.2. Caso de estudio……………………………………………………………………. 46
CONCLUSIONES…………………………………………………………….………… 53
RECOMENDACIONES………………………………………………………………... 55
REFERENCIAS BIBLIOGRÁFICAS………….…………………………………....… 56
ANEXOS……………………………………………………………………………….… 58
LISTA DE FIGURAS
Pág. 2.1. Estimaciones para la derivada de ( )f x en P utilizando diferencias hacia
delante, hacia atrás y central……………………………………………………………
14
2.2. Corte transversal de un medio dieléctrico polarizado………………………….. 18
2.3. Corte transversal de un material magnetizado………………………………….. 22
2.4. Componentes E y H en el modo TE bidimensional .……………….………….. 29
2.5. Componentes E y H en el modo TM bidimensional ………….……………….. 30
2.6. Interpretación grafica en espacio-tiempo de la componente en una
dimensión de las ecuaciones de Maxwell y su Discretización……………….….…..
31
2.7. Malla de Yee……………………………..…………………………………………. 32
2.8. Formulación bidimensional de la componente Hz …………………………...…. 33
2.9. Formulación bidimensional de la componente Ex ……………..………...……. 34
2.10. Formulación bidimensional de la componente Ey………………………..………….………… 34
2.11. Formulación bidimensional de la componente Ez ………….…………………. 35
2.12. Formulación bidimensional de la componente Hx……………………………………………... 36
2.13. Formulación bidimensional de la componente Hy……………………………………….…….. 37
3.1. Diagrama de flujo para la implementación del método de las diferencias
finitas en el dominio del tiempo (FDTD)……………………...………………………..
40
4.1. Algoritmo necesario para llevar a cabo la ejecución del programa
desarrollado……………………………………………………………………………….
47
4.2. Representación del caso de estudio……………………………………………… 48
4.3. Ilustración de un pulso gaussiano………………………………………………… 49
4.4. Representación de la malla de Yee para el caso de estudio………………….. 50
4.5. Valor de las pérdidas en los diferentes puntos del caso de estudio…………. 50
4.6. Error para las técnicas FDTD y Ray Tracing……………………………………. 52
INTRODUCCION
Durante los últimos años se ha dedicado un gran esfuerzo al desarrollo de
software de apoyo, para la construcción de redes de radio. Los mayores esfuerzos se
han dedicado a las aplicaciones al aire libre, dejando a un lado la investigación
relacionada con los ambientes interiores.
En la actualidad, la mayoría de los trabajos existentes relacionados con
propagación en interiores se dedica al estudio de la técnica del trazo de rayos, que
presta poca atención a la estructura de las paredes.
En este sentido, el método de las diferencias finitas en el dominio del tiempo
(FDTD), desarrollado por Kane Yee en 1966 para resolver las ecuaciones de Maxwell,
permite obtener un modelo de propagación adecuado para ambientes cerrados, a partir
de la ubicación de la fuente y una completa descripción del ambiente en términos de
sus parámetros dieléctricos.
El propósito de este trabajo es utilizar el método de las diferencias finitas en el
dominio del tiempo para simular la propagación de ondas electromagnéticas en
ambientes cerrados, considerando la estructura de las paredes. Asimismo, se desarrolla
un programa con fines didácticos, para evaluar el efecto de algunos parámetros clave
como estructura de las paredes y geometría del espacio, en la caracterización de los
canales de propagación en ambientes cerrados.
El trabajo está estructurado como se describe a continuación.
El capítulo I presenta el planteamiento del problema, los objetivos generales y
específicos, la justificación de la investigación y la delimitación del trabajo.
En el capítulo II se incluyen fundamentos teóricos del trabajo.
El capítulo III es el marco metodológico, donde se describe la metodología
aplicada.
Por último, el capítulo IV presenta el desarrollo del software y la implementación
del mismo a través de un caso de estudio.
Posteriormente, se incluyen las conclusiones y recomendaciones pertinentes al
trabajo.
CAPITULO I
EL PROBLEMA 1.1. PLANTEAMIENTO Y FORMULACIÓN DEL PROBLEMA
El fenómeno de propagación de ondas electromagnéticas puede ocurrir tanto al
aire libre como en ambientes cerrados. La propagación al aire libre se ve influenciada
por condiciones atmosféricas tales como nubes, lluvia, nieve, entre otros. En ambientes
cerrados, la propagación se ve afectada principalmente por los materiales de
construcción y por la configuración geométrica del espacio.
En el caso de propagación de ondas en ambientes cerrados debe considerarse la
reflexión, refracción y difracción de la onda transmitida, como consecuencia de
obstáculos presentes en la trayectoria de la señal. Esto hace que la señal llegue al
receptor por más de un camino. Este fenómeno que se conoce como múltiples
trayectorias (multipath en inglés [1]) se reduce al utilizar en las áreas de construcción
materiales absorbentes.
Durante las últimas dos décadas se ha desarrollado gran cantidad de software de
apoyo en la construcción de redes de radio. Sin embargo, con el crecimiento acelerado
de los sistemas de telefonía celular, los mayores esfuerzos se han dedicado a las
aplicaciones al aire libre. Durante la última década se ha prestado mayor atención a los
ambientes interiores. K. Chwung, J. Sau y R.D.Murch [2] presentan un nuevo modelo
empírico para la predicción de la propagación en ambientes interiores. R. Valenzuela [3]
introduce la técnica Ray Tracing para predecir la transmisión inalámbrica en ambientes
interiores y D. Lu y D. Rutledge [4] combinan las técnicas Ray Tracing y Plano E/H para
el modelado de canales inalámbricos en ambientes interiores. Estos trabajos se han
dedicado al desarrollo de la técnica del trazo de rayos que según W. Honcharenko y
H.L. Bertoni [5] y C. L.Holloway y P.L.Perini [6], prestan poca atención a la estructura
de las paredes.
12
En la actualidad existen varias herramientas de software diseñadas para simular la
propagación electromagnética, tanto en ambientes abiertos como cerrados. Dentro de
las más comunes se encuentran: EDX Signal Pro V5.0 (http://www.edx.com) y WinProp
(http://www.awe-communications.com), Gíreles Valley (http://www.wirelessvalley.com),
Cindoor (www.gsr.unican.es) [7]. Sin embargo, estas herramientas son costosas (miles
de dólares), y no están diseñadas para fines didácticos. Además, utilizan técnicas que
poco tienen en cuenta la estructura de las paredes.
En este sentido, el método de las diferencias finitas en el dominio del tiempo,
desarrollado por Kane Yee en 1966 [8], es una técnica simple y efectiva, que encuentra
aplicación en el modelado de la distribución del campo electromagnético dentro de un
espacio cerrado, considerando la estructura interna de las paredes.
Una solución por diferencias finitas involucra básicamente tres pasos: i) discretizar
una región de interés en grillas de nodos, ii) aproximar la ecuación diferencial dada por
una ecuación en diferencias finitas equivalente y iii) resolver la ecuación en diferencias
finitas, sujeta a las condiciones de contorno y/o condiciones iniciales.
El método de las diferencias finitas en el dominio del tiempo permitirá predecir el
comportamiento de la onda electromagnética al hacer contacto con los obstáculos
presentes en un ambiente cerrado, tales como paredes, mobiliario, etc.
De acuerdo a la problemática planteada, se tiene la siguiente formulación del
problema, dado un espacio cerrado, su arquitectura, los materiales de construcción y la
ubicación de una fuente, determine los niveles de potencia en los diferentes puntos del
espacio. 1.2. OBJETIVOS DE LA INVESTIGACIÓN
1.2.1. OBJETIVO GENERAL
Desarrollar un software que implemente el método de diferencias finitas en el
dominio del tiempo (FDTD), para simular la propagación de ondas electromagnéticas en
ambientes cerrados.
13
1.2.2. OBJETIVOS ESPECIFICOS
• Formular las ecuaciones de Maxwell de manera discreta, en dos
dimensiones, utilizando el método de diferencias finitas en el dominio del tiempo
(FDTD).
• Desarrollar un software para implementar el método de diferencias finitas en
el dominio del tiempo para simular la propagación bidimensional de ondas
electromagnéticas en ambientes cerrados, con condiciones de frontera establecidas.
• Comparar el desempeño del método de diferencias finitas en el dominio del
tiempo con el Método Ray Tracing utilizando un caso de estudio. Esta evaluación debe
tomar en cuenta tipos de fuentes, ubicación de las fuentes, y condiciones de frontera
establecidas según el material de construcción y la geometría del espacio.
1.3. JUSTIFICACIÓN DE LA INVESTIGACIÓN
En un sistema de comunicaciones es importante establecer la distribución de
potencia, en cada punto del espacio donde se esté llevando a cabo la transmisión de la
señal. Esto permite determinar el número óptimo de estaciones bases requeridas para
lograr que la transmisión de la señal alcance todos los puntos del espacio.
En este sentido, el método de las diferencias finitas en el dominio del tiempo
permite desarrollar un modelo de propagación adecuado para ambientes cerrados, a
partir de la ubicación de la fuente, y una descripción completa del ambiente a estudiar
[9].
Para los estudiantes de Ingeniería Eléctrica es importante contar con un software
que permita predecir el comportamiento de las ondas electromagnéticas y el nivel de
potencia en cualquier punto dentro de un entorno, para establecer la ubicación de las
radio bases necesarias para lograr un determinado rango de cobertura.
CAPITULO II
MARCO TEÓRICO
En este capítulo se establecen las bases teóricas sobre las cuales se desarrollará
la investigación.
2.1. MÉTODO DE DIFERENCIAS FINITAS
El método de Diferencias Finitas (Finite Difference Method, FDM, por sus siglas en
inglés) fue desarrollado por A.Thom en 1920 [10], bajo el nombre “El método de los
cuadrados”, para resolver ecuaciones hidrodinámicas no lineales. Desde entonces, el
método se ha utilizado para resolver diferentes tipos de problemas. Las técnicas de las
diferencias finitas están basadas en una serie de aproximaciones, que permiten
aproximar ecuaciones diferenciales por ecuaciones en diferencias finitas. Esas
aproximaciones en diferencias finitas se expresan en forma algebraica y relacionan el
valor de la variable dependiente en un punto en la región de solución, con el valor en
algunos puntos vecinos.
2.1.1. Esquema de las diferencias finitas
Considere una función ( )f x como la que se muestra en la Figura 2.1.
Figura 2.1. Estimaciones para la derivada de ( )f x en P utilizando diferencias hacia
delante, hacia atrás y central
15
La derivada de la función ( )f x en el punto P mide la rapidez de cambio de la
función en ese punto. En otras palabras, la derivada expresa que tan rápido crece (o
decrece) la función en el punto P. La derivada es equivalente a la pendiente de la recta
tangente a la función en ese punto.
Existen varias alternativas para aproximar la derivada de ( )f x en P; utilizando la
formula de la diferencia hacia delante
' 0 0
0( ) ( )( ) ,f x x f xf x
x+ Δ −Δ
(2.1)
utilizando la formula de la diferencia hacia atrás
' 0 0
0( ) ( )( ) ,f x f x xf x
x− + ΔΔ
(2.2)
o utilizando la formula de la diferencia central
' 0 0
0( ) ( )( )
2f x x f x xf x
x+ Δ − −Δ
⋅Δ
(2.3)
Asimismo, se puede aproximar la segunda derivada de ( )f x en P como:
( )
' ''' 0 0
0
0 0 0 0
0 0 02
( / 2) ( / 2)( )
( ) ( ) ( ) ( )1
( ) 2 ( ) ( )
f x x f x xf xx
f x x f x f x f x xx x x
f x x f x f x xx
+ Δ − −ΔΔ
+ Δ − − −Δ⎡ ⎤= −⎢ ⎥Δ Δ Δ⎣ ⎦+ Δ − + −Δ
= ⋅Δ
(2.4)
Cualquier aproximación de una derivada en términos de valores de la función en
un conjunto finito de puntos se conoce como Aproximación por Diferencias Finitas.
Una manera alternativa de obtener estas aproximaciones es utilizando la serie
Taylor. Así
( ) ( )2 3' '' '''0 0 0 0 0
1 1( ) ( ) ( ) ( ) ( ) ...2! 3!
f x x f x xf x x f x x f x+ Δ = + Δ + Δ + Δ + (2.5)
16
( ) ( )2 3' '' '''0 0 0 0 0
1 1( ) ( ) ( ) ( ) ( ) ...2! 3!
f x x f x xf x x f x x f x−Δ = −Δ + Δ − Δ + (2.6)
Sumando (2.5) y (2.6) se obtiene
( ) ( )2 4''
0 0 0 0( ) ( ) 2 ( ) ( )f x x f x x f x x f x x+ Δ + −Δ = + Δ + Δ (2.7)
donde ( )4xΔ es el error introducido por los truncamientos en las series. Este error es
del orden ( )4xΔ o simplemente ( )4xΔ . Además, ( )4xΔ representa los términos que
no son mayores que ( )4xΔ . Si xΔ es muy pequeño, estos términos pueden despreciarse
para obtener
( )'' 0 0 0
0 2
( ) 2 ( ) ( )( ) f x x f x f x xf xx
+ Δ − + −Δ
Δ (2.8)
que es igual a la ecuación (2.4). Restando la ecuación (2.6) de la ecuación (2.5), y
despreciando los términos de orden ( )3xΔ , se obtiene:
' 0 0
0( ) ( )( )
2f x x f x xf x
x+ Δ − −Δ
Δ (2.9)
que corresponde a la ecuación (2.3). Esto nos muestra que los principales errores en
las ecuaciones (2.3) y (2.4) están en el orden de ( )2xΔ . De manera similar, la fórmula
de diferencias en las ecuaciones (2.1) y (2.2) tienen errores por truncamiento del orden
( )xΔ . Pueden obtenerse aproximaciones por diferencias finitas de mayor orden,
considerando más términos en la serie de Taylor.
Una solución por diferencias finitas se realiza en tres pasos:
1. Diferenciar la variable independiente
2. Aproximar la ecuación diferencial por su equivalente en diferencias finitas
3. Resolver la ecuación en diferencias, sujeta a las condiciones de contorno y/o
condiciones iniciales.
17
2.2. ECUACIONES DE MAXWELL
La electrostática es el estudio de los efectos de las cargas eléctricas en reposo, y
de los campos eléctricos que no cambian con el tiempo.
La intensidad de campo eléctrico Er
se define como la fuerza por unidad de carga
que experimenta una carga de prueba estacionaria muy pequeña, al colocarse en una
región donde existe un campo eléctrico, es decir,
0q
FE Limq→
=r
r (V/m)
(2.10)
Los dos postulados fundamentales de la electrostática en el espacio libre son:
v
o
E ρε
∇ =r
(2.11)
0xE∇ =r
(2.12)
donde E∇r
y xE∇r
son la divergencia y el rotacional de Er
, respectivamente, vρ es la
densidad volumétrica de carga libre (C/m3) y oε es la permitividad del espacio libre.
En un problema electromagnético intervienen cuatro campos vectoriales
fundamentales: intensidad de campo eléctrico Er
(V/m), densidad de flujo eléctrico o
desplazamiento eléctrico Dr
(C/m2), densidad de flujo magnético Br
(T) e intensidad de
campo magnético Hr
(A/m).
En el espacio libre, los únicos vectores necesarios para analizar la electrostática
(los efectos de las cargas eléctricas estacionarias) y la magnetostática (los efectos de
las corrientes eléctricas estacionarias), son la intensidad de campo Er
y la densidad de
flujo de campo magnético Br
, respectivamente. Cuando hay presentes medios
materiales, los vectores desplazamiento Dr
e intensidad de campo magnético Hr
son
los más interesantes para el estudio de problemas eléctricos y magnéticos.
Al aplicar un campo eléctrico externo a un medio material aparece un vector de
polarización Pr
, que se ilustra en la Figura 2.3
18
1
0lim
n r
kk
v
pP
v
Δ
=
Δ →=
Δ
∑ rr
(C/m2) (2.13)
donde n es el número de moléculas por unidad de volumen, y el numerador representa
la suma vectorial de los momentos bipolares inducidos que están contenidos en un
volumen muy pequeño vΔ .
Figura 2.2. Corte transversal de un medio dieléctrico polarizado
En la Figura 2.2 puede verse que las moléculas contribuyen de forma efectiva a la
distribución de cargas superficiales positivas en la frontera a la derecha, y a la
distribución de cargas superficiales negativas en la frontera a la izquierda. Dado que la
densidad superficial de carga depende de la densidad de dipolos eléctricos que
sobresalen más allá de las líneas punteadas en una superficie, la densidad superficial
de carga de polarización equivalente es
ˆ.ps nP aρ =
r (C/m2). (2.14)
Para una superficie S que encierra un volumen V, la carga neta total que sale
fuera de V como resultado de la polarización se obtiene integrando la ecuación (2.14).
19
La carga neta que permanece dentro del volumen es el negativo de esta integral, es
decir
( ). .n pv
S V V
Q P a dS P d dν ρ ν= − = −∇ =∫ ∫ ∫r r r r (2.15)
donde se ha aplicado el teorema de la divergencia para convertir la integral de
superficie cerrada en una integral de volumen. Ahora, se puede definir la densidad
volumétrica de carga de polarización equivalente como
.pv Pρ = −∇r
(C/m3). (2.16)
Considerando que un dieléctrico polarizado da lugar a una densidad volumétrica
de carga equivalente pvρ , es de esperar que la intensidad de campo eléctrico en un
dieléctrico, debido a una distribución de fuentes dada, sea diferente de la intensidad de
campo en el espacio libre. Específicamente, hay que modificar la divergencia en la
ecuación (2.12) para incluir el efecto de pvρ , es decir
( )1v pv
o
E ρ ρε
∇ = + ⋅r
(2.17)
Utilizando la ecuación (2.16) se tiene
( ). o vE Pε ρ∇ + =
r r (2.18)
Se define ahora una nueva cantidad fundamental de campo, la densidad de flujo
eléctrico o desplazamiento eléctrico Dr
, como
oD E Pε= +
r r r (C/m2). (2.19)
El uso del vector D
r permite escribir una relación de divergencia entre el campo
eléctrico y la distribución de cargas libres en cualquier medio, sin tener que tratar de
manera explícita con el vector de polarización Pr
, ni con la densidad de carga de
20
polarización pvρ . Al combinar las ecuaciones (2.18) y (2.19), se obtiene la nueva
ecuación
. vD ρ∇ =r r (C/m3) (2.20)
Cuando las propiedades dieléctricas del medio son lineales e isótropas, la
polarización es directamente proporcional a la intensidad de campo eléctrico, y la
constante de proporcionalidad es independiente de la dirección del campo, esto es
0 eP Eε χ=
r r (2.21)
donde eχ es una cantidad sin dimensiones llamada susceptibilidad eléctrica. Un medio
dieléctrico es lineal si eχ es independiente de Er
, y homogéneo si eχ es independiente
de las coordenadas espaciales. Si se sustituye la ecuación (2.21) en la ecuación (2.19)
se obtiene
0 0(1 )e rD E E Eε χ ε ε ε= + = =
r r r r (C/m2) (2.22)
donde
0
1r eεε χε
= + = (C/m2) (2.23)
es una cantidad sin dimensiones conocida como permitividad relativa o constante
dieléctrica del medio. El coeficiente 0 rε ε ε= es la permitividad absoluta (con frecuencia
llamada simplemente permitividad) del medio, y se mide en farad por metro (F/m).
Las cargas en movimiento producen una corriente que a su vez crea un campo
magnético. En este caso hay que considerar el vector de densidad de flujo magnético,
B. Los dos postulados fundamentales de la magnetostática (campos magnéticos
estáticos) son
0B∇ =
r (2.24)
21
oxB Jμ∇ =r r
(2.25)
donde B∇ y xB∇ son la divergencia y el rotacional de B
r, respectivamente, oμ es la
permeabilidad del espacio libre, 74 10xπ − (H/m), y Jr
es la densidad de corriente (A/m2).
En ausencia de un campo magnético externo, los dipolos magnéticos de los
átomos de la mayoría de los materiales (con excepción de los imanes permanentes)
tienen orientaciones aleatorias, de manera que no hay momento magnético neto. La
aplicación de un campo magnético externo ocasiona, tanto la alineación de los
momentos magnéticos de los electrones giratorios, como un momento magnético
inducido, que se debe a un cambio en el movimiento orbital de los electrones. Para
obtener una fórmula que nos permita determinar el cambio cuantitativo en la densidad
de flujo magnético ocasionado por la presencia de un material magnético, definimos
kmr como el momento bipolar magnético de un átomo. Si hay n átomos por unidad de
volumen, se define un vector de magnetización Mr
como
1
0lim
n v
kk
v
mM
v
Δ
=
Δ →=
Δ
∑ rr
(A/m) (2.26)
que es la densidad de volumen del momento bipolar magnético.
De forma similar a la equivalencia de Pr
de los dipolos eléctricos inducidos, a una
densidad superficial de carga polarizada ˆ.ps nP aρ =rr y a una densidad volumétrica de
carga polarizada .pv Pρ = −∇rr , se puede demostrar de manera analítica la equivalencia
de Mr
de dipolos magnéticos a una densidad superficial de corriente de magnetización
ˆms nJ M x a=
r r (A/m) (2.27)
donde ˆna es la normal unitaria hacia fuera de la frontera, y una densidad de corriente de
volumen de magnetización
mvJ xM= ∇r r
(A/m2). (2.28)
22
En la Figura 2.3 puede observarse que los dipolos magnéticos en la superficie,
contribuyen de manera efectiva a una corriente superficial más allá de las líneas
punteadas. La magnitud de la corriente superficial es directamente proporcional a la
densidad de volumen del momento bipolar magnético, y la dirección de la corriente en
ambas fronteras está expresada correctamente por ˆnMxar
en la figura.
Figura 2.3. Corte transversal de un material magnetizado
Considerando que la aplicación de un campo magnético externo ocasiona tanto
una alineación de los momentos bipolares magnéticos como un momento magnético
inducido en un material magnético, se espera que la densidad de flujo magnético
resultante en presencia de un material magnético, sea diferente de su valor en el
espacio libre. El efecto macroscópico de la magnetización puede estudiarse
sustituyendo la densidad de corriente equivalente de volumen de magnetización mvJr
,
dada por la ecuación (2.28), en la ecuación (2.25). Se obtiene entonces
1
mvo o
BxB J J J xM x M Jμ μ
⎛ ⎞∇ = + = +∇ ⇒∇ − =⎜ ⎟
⎝ ⎠
rr r r r r r
(2.29)
23
Ahora, se define una nueva cantidad de campo fundamental, la intensidad de
campo magnético Hr
, como
o
BH Mμ
= −r
r r (A/m)
(2.30)
Al combinar las ecuaciones (2.29) y (2.30) se obtiene
xH J∇ =r r
(A/m2) (2.31)
donde J
r(A/m2) es la densidad de volumen de la corriente libre.
Cuando las propiedades magnéticas del medio son lineales e isótropas, la
magnetización es directamente proporcional a la intensidad de campo magnético, es
decir
mM Hχ=
r r (2.32)
donde mχ es una cantidad adimensional llamada susceptibilidad magnética. Al sustituir
la ecuación (2.32) en la ecuación (2.30), se obtiene la siguiente relación constitutiva
( )0 01 m rB H H Hμ χ μ μ μ= + = =
r r r r (Wb/m2) (2.33)
Que puede reescribirse como
1H Bμ
=r r
(A/m) (2.34)
donde
0
1r mμμ χμ
= + = (2.35)
es otra cantidad adimensional conocida como permeabilidad relativa del medio. El
parámetro 0 rμ μ μ= es la permeabilidad absoluta (en ocasiones simplemente
24
permeabilidad) del medio y se mide en H/m; mχ y, por consiguiente, rμ pueden ser
funciones de las coordenadas espaciales. En el caso de un medio simple (lineal,
isótropo y homogéneo), mχ y rμ son constantes.
Se observa que Er
(intensidad de campo eléctrico) y Dr
(desplazamiento eléctrico)
en el modelo electrostático, no están relacionados con Br
y Hr
en el modelo
magnetostático. En un medio conductor pueden existir campos eléctricos y magnéticos
estáticos, para formar un campo electromagnetostático. El campo eléctrico estático en
un medio conductor hace que fluya una corriente estacionaria, que a su vez produce un
campo magnético estático.
Los modelos estáticos son sencillos, pero inadecuados para explicar los
fenómenos electromagnéticos variables con el tiempo. Los campos eléctricos y
magnéticos estáticos no producen ondas que se propagan, y transportan energía e
información. Las ondas son la esencia de la acción electromagnética a distancia. En
condiciones variables con el tiempo, es necesario construir un modelo electromagnético
donde los vectores de campo eléctrico E y D estén relacionados con los vectores de
campo magnético B y H.
En este sentido, Michael Faraday consiguió uno de los mayores avances en la
teoría electromagnética, cuando en 1831 descubrió experimentalmente que se inducía
una corriente en una espira conductora cuando cambiaba el flujo magnético que
atravesaba la espira. La relación cuantitativa entre la fuerza electromotriz inducida y la
razón de cambio del flujo asociado, se conoce como ley de Faraday. Es una ley
experimental, y puede considerarse como un postulado.
El postulado fundamental de la inducción electromagnética es
BxEt
∂∇ = − ⋅
∂
rr
(2.36)
Esta ecuación (2.36) se aplica a todos los puntos en el espacio, ya sea el espacio
libre o un medio material. La intensidad de campo eléctrico en una región de densidad
de flujo magnético variable con el tiempo es por consiguiente no conservativa, y no
puede expresarse como el gradiente negativo de un potencial escalar. El postulado
fundamental de la inducción electromagnética nos asegura que un campo magnético
25
variable con el tiempo origina un campo eléctrico. Por lo tanto, se tiene ahora el
siguiente conjunto de ecuaciones, dos de rotacional
BxEt
∂∇ = −
∂
rr
(2.37)
xH J∇ =r r
(2.38)
y dos de divergencia
vD ρ∇ =r r (2.39)
0.B∇ =r
(2.40)
Adicionalmente, siempre debe satisfacerse el principio de conservación de la
carga, expresado por la ecuación de continuidad
vJtρ∂
∇ = − ⋅∂
rr (2.41)
Este principio debe satisfacerse en las ecuaciones (2.37) a (2.40), en una
situación variable con el tiempo. Si se toma la divergencia de la ecuación (2.39)
( ) 0xH J∇ ∇ = = ∇
r r (2.42)
Se observa que este principio no se satisface, por lo que hay que añadir un
término /v tρ∂ ∂ en el lado derecho de la ecuación (2.42) para obtener
( ) 0 vxH Jtρ∂
∇ ∇ = = ∇ +∂
rr r (2.43)
Al utilizar la ecuación (2.39) en la ecuación (2.43) se obtiene
( ) 0 . DxH Jt
⎛ ⎞∂∇ ∇ = = ∇ +⎜ ⎟∂⎝ ⎠
rr r
(2.44)
lo cual implica que
26
DxH Jt
∂∇ = +
∂
rr r
(2.45)
La ecuación (2.45) indica que un campo eléctrico variable con el tiempo producirá
un campo magnético, aunque no exista un flujo de corriente libre (es decir, incluso si
0J =r
). El término adicional /D t∂ ∂ es necesario para que la ecuación (2.45) sea
consistente con el principio de conservación de la carga.
Este término se denomina densidad de corriente de desplazamiento, y su
introducción en la ecuación xH∇r
fue una de las contribuciones principales de James
Clerk Maxwell (1831 – 1879). Para ser consistentes con la ecuación de continuidad en
una situación variable con el tiempo hay que generalizar las ecuaciones
BxEt
∂∇ = −
∂
rr
(2.46)
DxH Jt
∂∇ = +
∂
rr r
(2.47)
vD ρ∇ =r r (2.48)
0B∇ =r
(2.49)
En este caso vρ es la densidad volumétrica de cargas libres y J
r es la densidad de
corrientes libres, que pueden comprender tanto corriente de convección ( )vuρ como
corriente de conducción ( )Eσ .
En medios no-dispersivos e isótropos, ε y μ son escalares que no dependen del
tiempo, por lo que las ecuaciones de Maxwell se reducen a:
HxEt
μ ∂∇ = −
∂
rr
(2.50)
ExH Jt
ε ∂∇ = +∂
rr r
(2.51)
Eε ρ∇ =r r (2.52)
27
0Hμ∇ =r
(2.53)
donde J
r representa la densidad de corriente eléctrica equivalente J Eσ=
r r, donde σ es
la conductividad eléctrica. Por lo tanto, la ecuación (2.51) resulta en
ExH Et
σ ε ∂∇ = + ⋅∂
rr r
(2.54)
Los campos magnéticos tienen asociada una densidad de corriente magnética
equivalente mJr
, dada por *mJ Hσ=r r
, donde *σ representa la resistividad magnética. Al
considerar mJr
, la ley de Faraday puede reescribirse como
*
mH HxE J Ht t
μ μ σ∂ ∂∇ = − − = − − ⋅
∂ ∂
r rr r r
(2.55)
Si se expresa los campos magnético y eléctrico en coordenadas rectangulares
ˆ ˆ ˆ, ,x x y y z zH H e H e H e⎡ ⎤= ⎣ ⎦
r (2.56)
ˆ ˆ ˆ, ,x x y y z zE E e E e E e⎡ ⎤= ⎣ ⎦r
(2.57)
el rotacional de estos campos queda definido como
ˆ ˆ ˆx y z
x y z
e e e
xEx y z
E E E
∂ ∂ ∂∇ =
∂ ∂ ∂
r r
(2.58)
ˆ ˆ ˆx y z
x y z
e e e
xHx y z
H H H
∂ ∂ ∂∇ =
∂ ∂ ∂
r r
(2.59)
28
Sustituyendo las expresiones (2.58) y (2.59) en las ecuaciones (2.54) y (2.55)
tenemos
*
ˆ ˆ ˆ
1ˆ ˆ ˆ ˆ ˆ ˆ, , , ,
x y z
x x y y z z x x y y z z
x y z
e e e
H e H e H e H e H e H et x y z
E E E
σμ μ
∂ ∂ ∂ ∂⎡ ⎤ ⎡ ⎤= − −⎣ ⎦ ⎣ ⎦∂ ∂ ∂ ∂
(2.60)
ˆ ˆ ˆ
1ˆ ˆ ˆ ˆ ˆ ˆ, , , ,
x y z
x x y y z z x x y y z z
x y z
e e e
E e E e E e E e E e E et x y z
H H H
σε ε
∂ ∂ ∂ ∂⎡ ⎤ ⎡ ⎤= − + ⋅⎣ ⎦ ⎣ ⎦∂ ∂ ∂ ∂
(2.61)
Reescribiendo las ecuaciones vectoriales como un sistema de 6 ecuaciones
escalares se obtiene
*1 yx z
x
EH E Ht z y
σμ
∂⎛ ⎞∂ ∂= − −⎜ ⎟∂ ∂ ∂⎝ ⎠
(2.62)
*1y xzy
H EE Ht x z
σμ
∂ ∂∂⎛ ⎞= − −⎜ ⎟∂ ∂ ∂⎝ ⎠
(2.63)
*1 yxzz
EEH Ht y x
σμ
∂⎛ ⎞∂∂= − −⎜ ⎟∂ ∂ ∂⎝ ⎠
(2.64)
1 yx zx
HE H Et y z
σε
∂⎛ ⎞∂ ∂= − −⎜ ⎟∂ ∂ ∂⎝ ⎠
(2.65)
1y x zy
E H H Et z x
σε
∂ ∂ ∂⎛ ⎞= − −⎜ ⎟∂ ∂ ∂⎝ ⎠
(2.66)
1 y xzz
H HE Et x y
σε
∂⎛ ⎞∂∂= − − ⋅⎜ ⎟∂ ∂ ∂⎝ ⎠
(2.67)
2.2.1. Modos de propagación
Para que tenga lugar la propagación en la guía de ondas, la configuración de
campos eléctrico y magnético de la onda electromagnética debe satisfacer ciertas
condiciones. Hay muchas posibles configuraciones, llamadas modos, que se designan
29
según las direcciones de estos campos con respecto de la dirección de propagación.
Así tenemos los modos
• Transversal eléctrico (TE): en este modo el campo eléctrico es perpendicular a la
dirección de propagación, tal como se muestra en la Figura.2.4. En este caso, las
soluciones se derivan de la componente del campo magnético zHr
, con las condiciones
0zE =r
y 0x yH H= =r r
. Por lo tanto, se tiene el siguiente sistema de tres ecuaciones
*1 yxzz
EEH Ht y x
σμ
⎛ ⎞∂∂∂= − −⎜ ⎟⎜ ⎟∂ ∂ ∂⎝ ⎠
rrrr
(2.68)
1x zx
E H Et y
σε⎛ ⎞∂ ∂
= −⎜ ⎟∂ ∂⎝ ⎠
r rr
(2.69)
1y zy
E H Et x
σε
∂ ⎛ ⎞∂= − −⎜ ⎟∂ ∂⎝ ⎠
r rr
(2.70)
Figura 2.4. Componentes E y H en el modo TE
• Transversal magnético (TM): en este modo sólo el campo magnético es
perpendicular a la dirección de propagación, tal como se muestra en la Figura.2.5. En
este caso, las soluciones se derivan a través de la componente del campo eléctrico zEr
,
con la condición de que 0zH =r
y 0x yE E= =r r
, esto es, la componente axial del campo
magnético es cero. Con esto se asegura la transmisión de la potencia en la dirección z,
que es la que se ha seleccionado como la dirección de propagación de la onda.
En consecuencia, se tiene el siguiente sistema de ecuaciones:
30
*1x zx
H E Ht y
σμ⎛ ⎞∂ ∂
= − − −⎜ ⎟∂ ∂⎝ ⎠
r rr
(2.71)
*1y zy
H E Ht x
σμ
∂ ⎛ ⎞∂= −⎜ ⎟∂ ∂⎝ ⎠
r rr
(2.72)
1 y xzz
H HE Et x y
σε
⎛ ⎞∂ ∂∂= − −⎜ ⎟⎜ ⎟∂ ∂ ∂⎝ ⎠
r rrr
(2.73)
Figura 2.5. Componentes E y H en el modo TM
2.3. FORMULACION DE LAS ECUACIONES DE MAXWELL EN FDTD
En este trabajo se resuelven numéricamente las ecuaciones (2.68), (2.69), (2.70),
(2.71), (2.72), (2.73) aplicando el método de diferencias finitas y el algoritmo de Yee.
Considerando el caso sin pérdidas en la ecuación (2.72) se tiene
1y z
H Et xμ
∂ ⎛ ⎞∂= ⎜ ⎟∂ ∂⎝ ⎠
r r
(2.74)
Según la definición clásica de una derivada, la ecuación (2.74) puede escribirse
como:
0 0
1lim limy z
t x
H Et xμΔ → Δ →
Δ Δ=
Δ Δ
r r
(2.75)
31
La Figura 2.6 es una interpretación grafica en espacio-tiempo de la componente
en una dimensión de las ecuaciones de Maxwell. En el límite se obtiene una solución
continua y exacta a la ecuación (2.75) en el punto ( ),x t . Es importante señalar que en
este punto, las derivadas en tiempo y espacio se igualan, más no se iguala el valor
actual de los campos
La técnica FDTD impone una malla rectangular sobre la región de interés, y
resuelve una versión discreta de las ecuaciones de campo en los nodos de la malla. Sin
embargo, como los problemas de dinámica implican campos eléctricos y magnéticos
variantes en el tiempo, el rotacional en las ecuaciones de Maxwell de las leyes de
Ampere y Faraday debe resolverse en cada punto de la malla.
En 1966, K. S. Yee, introdujo un método para resolver estas ecuaciones, utilizando
una malla como la que se muestra en la Figura 2.7. Esta malla se conoce como malla
de Yee
Figura 2.6. Componente en una dimensión de las ecuaciones de Maxwell
32
Figura 2.7. Malla de Yee
Un examen cuidadoso de esta malla muestra que los puntos de solución para el
campo eléctrico, están espacialmente desfasados de los puntos de solución para el
campo magnético. Esto permite calcular el campo magnético o eléctrico en cada nodo,
usando los cuatro valores de campo magnético o eléctrico circundantes.
Además de hacer E y H discretos en el espacio, en el método FDTD los cambios
temporales de los campos se calculan a intervalos de tiempo discretos. Para
separaciones espaciales de malla de xΔ , yΔ , zΔ e incrementos de tiempo tΔ , los
campos se pueden escribir como:
( , , , ) ( , , , ) ( , , )nF x y z t F i x j y k z n t F x y z= Δ Δ Δ Δ = (2.76)
En primer lugar se formularán las ecuaciones para determinar los campos eléctrico
y magnético del modo transversal eléctrico bidimensional (TE).
Para determinar el valor del campo magnético en el eje z, se aplica diferencias
centrales a la ecuación (2.68) alrededor del punto ( 1/ 2, 1/ 2)i j+ + y nt , como se indica
en la Figura 2.8 para obtener
33
Figura 2.8. Formulación bidimensional de la componente Hz
1/ 2 1/ 2
*
( 1/ 2, 1) ( 1/ 2, )( 1/ 2, 1/ 2) ( 1/ 2, 1/ 2) 1
( 1, 1/ 2) ( , 1/ 2)
( 1/ 2, 1/ 2)
x xz z
y
z
n nn n
n ny
n
E i j E i jH i j H i jt y
E i j E i j
xH i j
μ
σ
+ − ⎡ + + − ++ + − + += −⎢
Δ Δ⎢⎣
+ + − +−
Δ⎤+ + ⎦
r rr r
r r
r
(2.77)
De esta ecuación se despeja el término 1/ 2 ( 1/ 2, 1/ 2)
z
nH i j+ + +r
y se obtiene
( 1 / 2 , 1) ( 1 / 2 , )1 / 2 1 / 2( 1 / 2 , 1 / 2 ) ( 1 / 2 , 1 / 2 )
( 1, 1 / 2 ) ( , 1 / 2 )
* ( 1 / 2 , 1 / 2 )
n nE i j E i jtn n x xH i j H i jz z y
n nE i j E i jy yx
nH i jz
μ
σ
⎡ + + − +Δ ⎢+ −+ + = + + + −⎢ Δ⎢⎣
+ + − +−
Δ
⎤+ + ⋅⎥⎦
r rr r
r r
r
(2.78)
Para determinar el valor del campo eléctrico en el eje x, se aplica diferencias
centrales a la ecuación (2.69) alrededor del punto ( 1/ 2, )i j+ y nt como se indica en la
Figura 2.9, para obtener
34
Figura 2.9. Formulación bidimensional de la componente Ex
1/ 2 1/ 2( 1/ 2, ) ( 1/ 2, ) ( 1/ 2, 1/ 2) ( 1/ 2, 1/ 2)1
( 1/ 2, )
x z z
n n n nx
nx
E i j E i j H i j H i jt y
E i j
ε
σ
+ − ⎡+ − + + + − + −= −⎢
Δ Δ⎢⎣⎤+ ⎦
r r r r
r
(2.79)
Despejando el término 1/ 2 ( 1/ 2, )nxE i j+ +r
de la ecuación (2.79) se obtiene
1/ 2 1/ 2 ( 1/ 2, 1/ 2) ( 1/ 2, 1/ 2)( 1/ 2, ) ( 1/ 2, )
( 1/ 2, )
z z
x
n nn nx
nx
H i j H i jtE i j E i jy
E i j
ε
σ
+ −⎡ + + − + −Δ
+ = + + −⎢Δ⎢⎣
⎤+ ⋅⎦
r rr r
r
(2.80)
Para determinar el valor del campo eléctrico en el eje y, se aplica diferencias
centrales a la ecuación (2.70) alrededor del punto ( , 1/ 2)i j + y nt como se indica en la
Figura 2.10, para obtener
Figura 2.10. Formulación bidimensional de la componente Ey
35
1/ 2 1/ 2 ( 1/ 2, 1/ 2) ( 1/ 2, 1/ 2)( , 1/ 2) ( , 1/ 2) 1
( , 1/ 2)
z z
n nn ny y
ny
H i j H i jE i j E i jt x
E i j
ε
σ
+ − ⎡ + + − − ++ − += − −⎢
Δ Δ⎢⎣⎤+ ⎦
r rr r
r
(2.81)
Despejando el término 1/ 2 ( , 1/ 2)n
yE i j+ +r
de la ecuación (2.81) se tiene
1/ 2 1/ 2 ( 1/ 2, 1/ 2) ( 1/ 2, 1/ 2)( , 1/ 2) ( , 1/ 2)
( , 1/ 2)
z z
n nn ny y
ny
H i j H i jtE i j E i jx
E i j
ε
σ
+ −⎡ + + − − +Δ
+ = + + − −⎢Δ⎢⎣
⎤+ ⎦
r rr r
r
(2.82)
Ahora se formularán las ecuaciones para determinar los campos eléctrico y
magnético del modo transversal magnético bidimensional (TM).
Para determinar el valor del campo eléctrico en el eje z, se aplica diferencias
centrales a la ecuación (2.71) alrededor del punto ( , )i j y 1/ 2nt t+ Δ , como se indica en
la Figura 2.11, para obtener
1/ 2 1/ 2( 1/ 2, ) ( 1/ 2, )1 1/ 2 1/ 2( , ) ( , ) 1 ( , 1/ 2) ( , 1 / 2)
1/ 2 ( , )
n nH i j H i jn n n nyE i j E i j H i j H i jyz z x xt yx
nE i jz
ε
σ
+ +⎡ + − −+ + +⎢− + − −= −⎢Δ ΔΔ⎢
⎣
+ ⎤− ⋅⎥⎦
r rr r r r
r
(2.83)
Figura 2.11. Formulación bidimensional de la componente Ez
36
Despejando el término 1( , )nzE i j+r
de la ecuación (2.83) se tiene
1/ 2 1/ 2 1/ 2 1/ 21
1/ 2
( 1/ 2, ) ( 1/ 2, ) ( , 1/ 2) ( , 1/ 2)( , ) ( , )
( , )
y
n n n nyn n x x
z z
nz
H i j H i j H i j H i jtE i j E i jx y
E i j
ε
σ
+ + + ++
+
⎡ + − − + − −Δ= + ⎢ −
Δ Δ⎢⎣⎤− ⎦
r r r rr r
r
(2.84)
Para determinar el valor del campo magnético en el eje x, se aplican diferencias
centrales a la ecuación (2.72) alrededor del punto ( , 1/ 2)i j + , y 12nt t+ Δ , como se ilustra
en la Figura 2.12, para obtener
1 1/ 2 1/ 2* 1/ 2( , 1/ 2) ( , 1/ 2) ( , 1) ( , )1 ( , 1/ 2)x
n n n nx nz z
x
H i j H i j E i j E i j H i jt y
σμ
+ + +++ − + ⎡ ⎤+ −
= − − +⎢ ⎥Δ Δ⎣ ⎦
r r r rr (2.85)
Despejando el término 1( , 1/ 2)
x
nH i j+ +r
de la ecuación (2.85) se tiene
1/ 2 1/ 2
1 * 1/ 2( , 1) ( , )( , 1/ 2) ( , 1/ 2) ( , 1/ 2)x
n nn n nz z
x xE i j E i jtH i j H i j H i j
yσ
μ
+ ++ +⎡ ⎤+ −Δ
+ = + + − − + ⋅⎢ ⎥Δ⎣ ⎦
r rr r r (2.86)
Para determinar el valor del campo magnético en el eje y se aplica diferencias
centrales a la ecuación (2.73) alrededor del punto ( 1/ 2, )i j+ y 12nt t+ Δ , como se ve en
la Figura 2.13 para obtener
Figura 2.12. Formulación bidimensional de la componente Hx
37
Figura 2.13. Formulación bidimensional de la componente Hy
1 1/ 2 1/ 2* 1/ 2( 1/ 2, ) ( 1/ 2, ) ( 1, ) ( , )1 ( 1/ 2, )
n n n ny y nz z
y
H i j H i j E i j E i j H i jt x
σμ
+ + +++ − + ⎡ ⎤+ −
= − +⎢ ⎥Δ Δ⎣ ⎦
r r r rr
(2.87)
Despejando el término 1( 1/ 2, )n
yH i j+ +r
de la ecuación (2.87) tenemos
1/ 2 1/ 2
1 * 1/ 2( 1, ) ( , )( 1/ 2, ) ( 1/ 2, ) ( 1/ 2, )n n
n n nz zy y y
E i j E i jtH i j H i j H i jx
σμ
+ ++ +⎡ ⎤+ −Δ
+ = + + − +⎢ ⎥Δ⎣ ⎦
r rr r r
(2.88)
2.4. ESTABILIDAD NUMÉRICA
El tamaño de la grilla debe elegirse de manera que los campos electromagnéticos
no cambien sustancialmente de un nodo a otro de la grilla. Para esto, la dimensión de la
grilla deberá ser una fracción de la longitud de onda λ. En general se recomienda
utilizar /10λ y considerar x yΔ = Δ .
La estabilidad de la solución se obtiene aplicando el criterio de Courant [8], el cual
establece que
/ 1c t LΔ > (2.89)
donde c es la velocidad de la luz, L una medida lineal del elemento como el largo o
ancho de la celda y tΔ es el intervalo de tiempo.
El criterio de estabilidad de Courant para 2 dimensiones es
38
( ) ( )2 2
11
tc
x y
Δ ≤ ⋅
Δ + Δ
(2.90)
Una vez elegida la grilla, la condición anterior impone una restricción para el
intervalo de tiempo tΔ .
En el caso que x yΔ = Δ = Δ , el criterio se convierte en
2t
cΔ
Δ ≤ (2.91)
2.5. TECNICA RAY TRACING La esencia del método Ray tracing es la generación y la descripción de los rayos.
En este método se considera un paquete que contiene todos los rayos
transmitidos que pueden o no llegar al receptor. El número de rayos considerados y la
distancia del receptor al transmisor determinan la resolución espacial posible y así
mismo la exactitud de este modelo.
Se selecciona un número finito de posibles direcciones de propagación. Un rayo
es enviado por cada una de las direcciones, si uno de los rayos choca contra algún
objeto, entonces se genera un rayo reflejado y uno refractado. La esfera de recepción
con un radio correcto puede describir la región en la cual se recibe uno solo. Si el radio
es muy grande se pueden recibir dos. Si el radio es muy pequeño es probable que
ninguno de los dos rayos llegue a la esfera de recepción.
En dos dimensiones todos los rayos o tubos de rayos son sectores de rayos.
Desde la fuente los rayos son emitidos a lo largo de diferentes direcciones con el mismo
ángulo de sector en el plano. Si el ángulo se selecciona pequeño, esto nos dará una
gran exactitud e implicará un tiempo largo de cálculo.
Cada rayo es emitido desde la fuente y puede ser remontado en un árbol binario.
Las intersecciones con las superficies de objetos se representan como nodos en el
árbol, del rayo incidente se descomponen un rayo reflejado en el objeto y otro
penetrado en el objeto, el punto de difracción es procesado como una fuente.
39
El proceso de descomposición se repite de manera recursiva, y este
procedimiento se repite hasta que los rayos se hacen más débiles que el valor de
umbral o hasta que el rayo salga fuera del área de propagación de interés o que dicho
rayo sea recibido.
CAPITULO III
MARCO METODOLÓGICO
Este capítulo incluye la metodología utilizada para llevar resolver el problema
planteado.
3.1. METODOLOGÍA APLICADA
La Figura 3.1 muestra un diagrama de flujo para la implementación del método
de las diferencias finitas en el dominio del tiempo (FDTD), para la simulación de la
propagación bidimensional de ondas electromagnéticas.
Figura 3.1. Diagrama de flujo para la implementación del método de las diferencias finitas en el dominio del tiempo (FDTD)
Particiones en espacio y tiempo
Calculo del campo eléctrico en todos los puntos espaciales
Inserción de la fuente en el campo eléctrico
Condiciones de frontera
Calculo del campo magnético en todos los puntos espaciales
Gráficas
Paredes de la habitación
Arquitectura de la habitación
Algoritm
o numérico
Constantes del espacio libre y frecuencia de operación
41
A continuación se presenta una explicación detallada de cada uno de los pasos
del diagrama de flujo de la Figura 3.1.
3.1.1. Definición de constantes
Inicialmente es necesario definir las constantes a utilizar. Estas constantes
incluyen, en primer lugar, los parámetros asociados al espacio libre tales como
velocidad de la luz , permeabilidad del espacio libre, y permitividad del espacio libre, y
en segundo lugar, los parámetros asociados al tipo de fuente a utilizar tales como
longitud de onda, período y ancho del pulso.
3.1.2. Particiones en espacio y tiempo
En este punto debe indicarse que criterio será utilizado para la selección del
tamaño de la celda.
3.1.3. Definición del medio dieléctrico En este paso se define la geometría del espacio cerrado a estudiar. Esta
geometría incluye dimensiones de las paredes y puertas, así como también las
propiedades dieléctricas de las mismas (permitividad relativa ( rε ) y conductividad (σ )).
3.1.4. Formulación de las ecuaciones de campo eléctrico y magnético en dos dimensiones Las ecuaciones para el cálculo de los campos eléctrico y magnético deben
definirse para la propagación en dos dimensiones, en este caso debe seleccionarse con
cual modo de propagación se calcularán los campos: trasversal magnético (TE) o
trasversal eléctrico (TE).
3.1.5. Condiciones de frontera
Las condiciones de frontera absorbentes son necesarias para simular
adecuadamente la propagación de los campos electromagnéticos, una vez que éstos
han llegado a la frontera de la malla de simulación.
CAPITULO IV
DESARROLLO DEL SOFTWARE Y ANÁLISIS DE RESULTADOS
Este capítulo presenta el desarrollo del software y el análisis de resultados, para
un caso de estudio previamente seleccionado.
4.1. DESARROLLO DEL SOFTWARE Como punto de partida en el software se definen las constantes a utilizar (Anexo
A). Aquí se indican los parámetros del espacio libre tales como velocidad de la
luz, 83x10 m/sc = , permeabilidad del espacio libre, -70 4 x10 H/mμ = y permitividad del
espacio libre, ,( -120 8,85419x10 F/mε = , los parámetros asociados al tipo de fuente a
utilizar tales como longitud de onda (λ ), período (τ ) y ancho del pulso (δ ).
Luego debe establecerse el criterio para seleccionar el tamaño de cada celda de la
malla y cuantas particiones en espacio y tiempo serán necesarias para construir la
malla de Yee. Para esto deberá tenerse en cuenta el criterio de Courant, mencionado
en el punto 2.4, el cual establece que el ancho y largo recomendable para la malla es
/10λ y la partición temporal deberá ser / 2t cΔ ≤ Δ . En el software desarrollado se
utilizó /10λ y / 2t cΔ = Δ .
Una vez definidas las particiones espacial y temporal deberá definirse la geometría
del espacio cerrado a estudiar. Es decir, se deberá establecer en cuales puntos de la
malla se encuentran ubicadas las paredes que componen el espacio a estudiar.
Además deberá indicarse las propiedades dieléctricas de las mismas como lo son
permitividad relativa ( rε ) y conductividad (σ ).
El siguiente paso en el desarrollo del software será la formulación de las
ecuaciones para el campo eléctrico y magnético. En este caso se formularán para el
modo transversal magnético bidimensional (TM). Recordando las ecuaciones se tiene,
43
1x zH Et yμ
⎛ ⎞∂ ∂= − ⎜ ⎟∂ ∂⎝ ⎠
r r
(2.71)
1y zH Et xμ
∂ ⎛ ⎞∂= ⎜ ⎟∂ ∂⎝ ⎠
r r
(2.72)
1 y xzz
H HE Et x y
σε
⎛ ⎞∂ ∂∂= − −⎜ ⎟⎜ ⎟∂ ∂ ∂⎝ ⎠
r rrr
(2.73)
Las ecuaciones para campo eléctrico y magnético deben tener la mismas
unidades. Considerando que ε y μ tienen unidades de (F/m) y (H/m), respectivamente,
zE y yH resultarán unidades diferentes. Entonces se hace necesario el siguiente
cambio de variable [11]:
0
0
E Eεμ
= ⋅r
% (4.1)
Sustituyendo la ecuación (4.1) en las ecuaciones (2.71), (2.72) y (2.73) se obtiene
0 0
1x zH Et yε μ
⎛ ⎞∂ ∂= − ⎜ ⎟∂ ∂⎝ ⎠
r%
(4.2)
0 0
1y zH Et xε μ
∂ ⎛ ⎞∂= ⎜ ⎟∂ ∂⎝ ⎠
r%
(4.3)
0 0
1 y xzz
r or
H HE Et x y
σε εε ε μ
⎛ ⎞∂ ∂∂= − −⎜ ⎟⎜ ⎟∂ ∂ ∂⎝ ⎠
r r%
% (4.4)
La aproximación por diferencias finitas de las ecuaciones (4.2), (4.3) y (4.4) está
dada por
( ) ( ) ( ) ( )1 1/ 2 1/ 2
0 0
, 1/ 2 , 1/ 2 , 1 ,1x x z z
n n n nH i j H i j E i j E i jt yε μ
+ + +⎛ ⎞+ − + + −= − ⎜ ⎟⎜ ⎟Δ Δ⎝ ⎠
r r% %
(4.5)
44
( ) ( ) ( ) ( )1 1/ 2 1/ 2
0 0
1/ 2, 1/ 2, 1, ,1y x z z
n n n nH i j H i j E i j E i jt xε μ
+ + ++ − + ⎛ ⎞+ −= ⎜ ⎟⎜ ⎟Δ Δ⎝ ⎠
r r% %
(4.6)
( ) ( ) ( ) ( )
( ) ( ) ( )
1/ 2 1/ 2
0 0
1/ 2, 1/ 2,, , 1
, 1/ 2 , 1/ 2,
y yz z
x
z
n nn n
r
n nx n
r o
H i j H i jE i j E i jt x
H i j H i jE i j
y
ε ε μ
σε ε
+ − ⎡ + − −−= ⎢
Δ Δ⎢⎣⎤+ − −
− −⎥Δ ⎥⎦
r r% %
r r
%
(4.7)
respectivamente.
En la ecuación (4.7) de campo eléctrico, se observa puede ver que todos los campos
del lado derecho de la ecuación están evaluados en el paso temporal n. Considerando
que el valor de zE en ese paso, no está almacenado en la memoria del computador (se
hace necesario estimar ese valor, se supone que sólo el valor previo de E en el paso
temporal 1/ 2n − está almacenado en la memoria del computador), se hace necesario
estimar este valor. Un buen estimado es el promedio a través de dos pasos temporales,
es decir
( ) ( ) ( )1/ 2 1/ 2, ,,
2
n nz zn
z
E i j E i jE i j
+ −+=% %
% (4.8)
Suponiendo que el tamaño de las celdas para x e y son iguales, es decir, que
x yΔ = Δ , se tiene
( ) ( )( )
( ) ( ) ( )
01/ 2 1/ 2
0 00 0
1 1/ 2,2, ,
1 12 2
1/ 2, , 1/ 2 , 1/ 2
y
z z
y x
nrn n
rr r
n n nx
tH i jtE i j E i j
xt tx
H i j H i j H i jx y
σε εσ σε ε μ
ε ε ε ε
+ −
⎛ ⎞Δ−⎜ ⎟ ⎡ +Δ⎝ ⎠= + ⎢Δ⎛ ⎞ ⎛ ⎞Δ Δ ⎢⎣+ Δ +⎜ ⎟ ⎜ ⎟
⎝ ⎠ ⎝ ⎠⎤− + − −
− − ⎥Δ Δ ⎥⎦
r
% %
r r r
(4.9)
Aplicando las condiciones de estabilidad dadas en la ecuación (2.91)
2t
cΔ
Δ ≤ (2.91)
45
Recordando la ecuación que define la velocidad de la luz
00 0
1cε μ
= (4.10)
Sustituyendo la ecuación (4.10) en la ecuación (2.91) se tiene
0 0
1 1/ 2txε μ
Δ=
Δ
(4.11)
Aplicando esta condición de estabilidad dada por (4.11) en las ecuaciones (4.5),
(4.6) y (4.9) se tiene
( ) ( ) ( ) ( )1 1/ 2 1/ 21, 1/ 2 , 1/ 2 , 1 ,2x x z z
n n n nH i j H i j E i j E i j+ + +⎡ ⎤+ = + − + −⎣ ⎦r r
% % (4.12)
( ) ( ) ( ) ( )1 1/ 2 1/ 211/ 2, 1/ 2, 1, ,2y x z z
n n n nH i j H i j E i j E i j+ + +⎡ ⎤+ = + + + −⎣ ⎦r r
% % (4.13)
( ) ( ) ( ) ( )
( ) ( )
01/ 2 1/ 2
0 0
12 1/ 2, , 1/ 2, 1/ 2,
1 12 2
, 1/ 2 , 1/ 2
z z y y
x
rn n n n
rr r
n nx
t
E i j E i j H i j H i jt t
H i j H i j
σε εσ σε
ε ε ε ε
+ −
⎛ ⎞Δ−⎜ ⎟⎝ ⎠ ⎡= + + − −⎣⎛ ⎞ ⎛ ⎞Δ Δ+ +⎜ ⎟ ⎜ ⎟
⎝ ⎠ ⎝ ⎠⎤− + − − ⎦
r r% %
r r
(4.14)
que serán las ecuaciones que se implementarán en el software para el calculo de los
campos.
De igual forma deben establecerse las condiciones de frontera absorbentes. Es
decir, deberá establecerse una condición toda vez que los campos han llegado a la
frontera de la malla de simulación. En el algoritmo de las diferencias finitas, los valores
de los campos se determinan mediante un promedio de los campos en los puntos
vecinos. El problema consiste en que en la frontera de la malla de simulación, este
promedio no se puede dar porque no se conoce el valor del campo fuera de la malla.
46
Físicamente, se espera que los campos se propaguen hacia fuera de la malla de
simulación, ya que se considera que fuera de la malla no existen fuentes [18].
La distancia que la onda viaja en un paso temporal está determinada por la
relación:
( )0 0 0tan / 2 / 2,dis cia c t c z c z= Δ = Δ = Δ (4.15)
lo que implica que el campo toma un par de pasos temporales, para viajar un paso
espacial. En este sentido, Sullivan [11] determinó las condiciones de borde en el
extremo izquierdo de la malla como
( ) ( )21,1 2, 2 .n n
z zE E −= (4.16)
De manera similar, para el lado derecho, la condición de borde será
( ) ( )2, 1, 1 .n n
z zE Nx Ny E Nx Ny−= − − (4.17)
Con esta expresión se asigna el valor del campo en la frontera en lugar de
calcularlo. Esta condición es fácil de implementar, ya que solo es necesario guardar el
valor de Ex(2,2) un par de pasos temporales y luego asignarlo a Ex(1,1). Para el lado
derecho de la malla se procede de manera similar.
La Figura 4.1 ilustra el diagrama de flujo necesario para llevar a cabo la corrida del
software. En esta figura se observa que al inicio del software deberán introducirse todos
los valores necesarios tanto de la malla de estudio como de la geometría del espacio a
estudiar. Posteriormente existe un lazo para el cálculo de los campos eléctrico y
magnético a partir de una fuente que generará la onda electromagnética. Este lazo
dependerá de los pasos espaciales que se coloquen al inicio del software, es decir, se
repetirá tantas veces como sea el numero de pasos temporales colocados al inicio del
programa, en este caso, se colocará el doble para los pasos temporales, que será el
tiempo que le tomará a la onda recorrer toda la malla.
48
4.2. CASO DE ESTUDIO: PROPAGACION BIDIMENSIONAL EN UN PISO DE UN EDIFICIO, COMPUESTO POR PAREDES DE CONCRETO Y MADERA Este caso de estudio fue tomado del trabajo realizado por Z. Ji, B. Li, H. Wang, H.
Chen y Y. Zhau, titulado “A New Indoor Ray-Tracing Propagation Prediction Model” [1],
en el cual se emplea la técnica Ray Tracing para estudiar la propagación de ondas en
un piso de un edificio compuesto por 20 salones. En este piso se toman 39 puntos,
distribuidos a lo largo del mismo, para hacer las mediciones respectivas (Figura 4.2). El
ancho del piso es 18.26 m, mientras que el largo es 76.39m. Las paredes que dividen
los salones están hechas de madera. Las paredes laterales de los pasillos están
hechas principalmente de concreto. Todas las puertas de los salones son de madera.
Figura 4.2. Representación del caso de estudio
El campo electromagnético que se genera de la fuente (T ), se encuentra ubicada
dentro del piso y será un pulso gaussiano que tiene la forma
20
0
0, 0( )
, 0t t
tE t
E tδ−⎛ ⎞−⎜ ⎟
⎝ ⎠
⎧ <⎪= ⎨⎪ >⎩ l
que asemeja un pulso de luz (Figura 4.3). El pulso tiene una amplitud máxima de 0E ,
está centrado en ot , sólo toma valores diferentes de cero en la vecindad de 0t , y tiene
una duración media δ . En este caso se seleccionaron 0 0,5 * ºt N pasos temporales= y
25 segδ = .
49
La implementación de las ecuaciones (4.12), (4.13) y (4.14) en un algoritmo
computacional requiere tres vectores bidimensionales para describir la coordenada
espacial, uno para zE , otro para xH y otro para yH .
Figura 4.3. Ilustración de un pulso gaussiano
Se construyó una malla de Yee con 4400 particiones en x y 1100 particiones en y.
Se utilizó una frecuencia de 1,7 Ghz, con lo cual la longitud de onda es 0,1765o mλ = ,
con particiones espaciales /10od λ= y particiones temporales delta_t = d/2c (Figura 4.4).
Las paredes de concreto poseen un espesor de 12d , con εr= 3.25. Las puertas y
paredes de madera tienen un espesor 3d y εr= 1.94.
El algoritmo se ejecuta para 7000 pasos de tiempo, y el tiempo de ejecución es
aproximadamente 3 horas.
La Figura 4.5 muestra los valores de las pérdidas en los puntos indicados en la
Figura 4.2. Se muestran los valores de las pérdidas medidos y los calculados utilizando
el algoritmo FDTD y la técnica Ray Tracing en el área en estudio.
50
Figura 4.4. Representación de la malla de Yee para el caso de estudio
0 5 10 15 20 25 30 35 40-60
-50
-40
-30
-20
-10
0
Puntos de medición
dB
MedicionRay-TracingFDTD
Figura 4.5. Valor de las pérdidas en los diferentes puntos del caso de estudio
x
y
Δx
Malla de Yee
1 2 Nx ………… 0.0176 m
Ny
………
Δy 0.0176 m
1
2
t = x/2c=2,9417e-11Δ
51
De la Figura 4.5 se observa que aquellos puntos que tienen línea de vista directa
con la fuente (18, 19, 20, 21 y 39), tienen valores cercanos de pérdidas en los tres
casos, medido, técnica FDTD y técnica ray tracing, lo cual se puede verificar al
comparar, por ejemplo, los valores arrojados por el punto 20, indicando que el valor de
las pérdidas medidas es de 15,02 dB, las pérdidas arrojadas por la técnica FDTD es
14,93 dB y las pérdidas arrojadas por la técnica FDTD es 14,03 dB, lo cual expresa un
error porcentual de 0,6 % para la técnica FDTD y 2,06 % para la técnica ray tracing.
Asimismo, en la Figura 4.5, se observa que los puntos que no están ubicados
dentro de ningún salón, sino ubicados en el área del pasillo, y se encuentran a una
mayor distancia del punto fuente y que por la forma en la que la onda se propaga, con
influencia de la arquitectura de las paredes, se tienen valores de pérdidas que varían
entre la técnica FDTD y la técnica Ray Tracing, tal y como se observa en los valores de
pérdidas en el punto 26, el cual tiene 26,52 dB como valor medido, 26,07 dB como valor
estimado por la técnica FDTD y 24,01 dB como valor estimado por la técnica ray
tracing, expresando un error porcentual de 1,7 % para la técnica FDTD y 9,46 % para la
técnica ray tracing.
En los puntos ubicados dentro de los salones en los cuales la onda
electromagnética debe atravesar paredes completas, se ve que aquellos puntos que
están ubicados a una mayor distancia del punto fuente son aquellos que tienen
mayores diferencias entre los valores estimados por ambas técnicas, mientras que
aquellos que se encuentran en los salones cercanos al punto fuente tienen menores
diferencias entre los valores estimados por ambas técnicas. Si se observa el punto 3, el
cual es el punto más alejado del punto fuente, los valores de pérdidas son 48,95 dB,
como valor medido, 47,23 dB como valor estimado para la técnica FDTD y 41 dB para
la técnica ray tracing, indicando un error porcentual de 3,51 % para FDTD y 16,24 %
para ray Tracing. Este hecho se debe a que la onda electromagnética deberá atravesar
un mayor número de paredes para llegar a aquellos puntos que se encuentran en los
salones más lejanos a la fuente, mientras que para llegar a los salones más cercanos a
la fuente, la onda deberá atravesar un menor número de paredes.
En la Figura 4.6 se muestran las curvas de error tanto para la técnica FDTD como
para la técnica Ray Tracing. En ella se observa que la técnica Ray Tracing ofrece un
mayor porcentaje de error en la estimación de las pérdidas.
52
0 5 10 15 20 25 30 35 400
5
10
15
20
25
30
35
40
45
50
Puntos de medición
% E
rror
FDTD
Ray-Tracing
Figura 4.6. Error para las técnicas FDTD y Ray Tracing
La Tabla 4.1 muestra los valores estadísticos para el error que ofrecen ambas
técnicas estudiadas. En ella se observa que la técnica FDTD ofrece una mejor
estimación de las pérdidas. Esto se demuestra al observar el valor promedio del error,
el cual en la técnica ray tracing es aproximadamente el 60 % mayor que el valor
promedio del error en la técnica FDTD. De igual forma se observa que el valor máximo
del error en la técnica Ray Tracing es aproximadamente 3 veces el valor máximo del
error en la técnica FDTD. Asimismo se observa que la desviación estándar es mucho
mayor, aproximadamente 3 veces, en la técnica Ray Tracing que en la técnica FDTD, lo
cual nos indica que hay mayor variabilidad en el error dado por la técnica Ray Tracing.
Tabla 4.1.Valores estadísticos para el error en las técnicas FDTD y Ray Tracing
FDTD Ray Tracing
Promedio del error 5,87 13,67
Valor máximo de error 16,34 46,77
Valor mínimo de error 0,56 1,13
Desviación estándar del error 3,66 9,69
CONCLUSIONES
En este trabajo se desarrolló un software para implementar el método de las
diferencias finitas en el dominio del tiempo, para simular la propagación de ondas
electromagnéticas en un espacio cerrado, considerando la arquitectura del espacio y
sus propiedades dieléctricas, y la ubicación de la fuente.
El método de diferencias finitas se aplicó a un caso de estudio, arrojando
resultados que indican que aquellos puntos que se encuentran mas cercanos a la
fuente y tienen línea de vista con la misma, es donde existen menores valores de
pérdidas. En este caso ambas técnicas, FDTD y Ray Tracing, ofrecen estimaciones con
valores muy cercanos a los valores medidos.
Asimismo se demostró que aquellos puntos en los cuales la onda
electromagnética debe atravesar mayor número de paredes, es donde existen mayores
valores de pérdidas, ofreciendo la técnica FDTD una mayor estimación, en dichos
puntos, que la técnica Ray Tracing, al compararse estos resultados con valores
medidos.
El desempeño del método FDTD se comparó con la técnica Ray Tracing,
graficando para ambos casos la curva del error ofrecido para ambas técnicas, la cual
indica que la técnica FDTD ofrece una mejor estimación de las pérdidas, principalmente
en aquellos casos en los que la onda debe atravesar mayor número de obstáculos.
Se demostró que una de las ventajas más significativas de la técnica FDTD es que
toma en cuenta las propiedades dieléctricas de cada una de las paredes que componen
el espacio cerrado, estimando con mayor precisión los valores de las pérdidas. Además,
es posible predecir el comportamiento de las ondas electromagnéticas en cada uno de
54
los puntos del área de estudio.
Una limitación del método es que la estabilidad del mismo depende del tamaño de
la grilla utilizada en la discretización. Es decir, mientras mas grande sea el tamaño de la
grilla, mayor será el tiempo utilizado para realizar los cálculos y mayor será el
requerimiento computacional.
RECOMENDACIONES
Al culminar el trabajo, se plantean las siguientes recomendaciones
Extender el trabajo a la propagación de la onda electromagnética en tres
dimensiones, considerando además el modelado en diferencias finitas de la fuente.
Analizar la posibilidad de simular el algoritmo desarrollado en este trabajo en otros
lenguajes de programación diferentes al Matlab, tales como Fortran, con la finalidad de
comparar los resultados arrojados por distintos software, tanto en tiempo de ejecución
como en exactitud de resultados.
Crear en la Escuela de Ingeniería Eléctrica de LUZ un laboratorio de investigación
en el área de propagación, el cual contenga los equipos de medición necesarios para
validar los trabajos de investigación desarrollados en el área.
REFERENCIAS BIBLIOGRAFICAS
[1] Z. Ji, B. Li, H. Wang, H. Chen y Y. Zhau. “A New Indoor Ray-Tracing Propagation Prediction Model”. Computational Electromagnetics and Its Applications. 540-542, 1999.
[2] K. Chwung, J. Sau y R.D.Murch. “A New Empirical Model Indoor Propagation Prediction”, IEEE Trans. Vehic. Tech.,1997.
[3] R. Valenzuela. “Ray Tracing approach to predicting indoor wireless transmission”, Proceedings of the 43rd IEEE Vehicular Technology Conference, p. 214–218, May 1993.
[4] D. Lu y D. Rutledge. “Indoor Wireless Channel Modeling from 2.4 to 24GHz Using a Combined E/H-Plane 2D Ray Tracing Method”, Int. Symp. on Ant. and Prop., Monterey, CA., 2004.
[5] W. Honcharenko y H.L. Bertoni. “Transmission and reflection characteristics at concrete block walls in the UHF bands proposal for future PCs”, IEEE Trans. On Antennas and Propagation Vol. 42, 232-239, Febrero 1994.
[6] C. L.Holloway, P.L.Perini. “Analysis of Composite Walls and Their Effects on Short-Path Propagation Modeling”, IEEE Trans. On Vehicular Technology Vol.46 NO. 3, 730-738, Agosto 1997.
[7] Castellanos E., Talero J. B., Rugeles J. y Ortega H., “Análisis De Propagación Electromagnética En Espacios Cerrados: Herramienta Software En Matlab Para Predicción Y Simulación”, Revista Colombiana de Tecnologías de Avanzada, Volumen 2 Número 6 año 2005.
[8] Yee K.S. “Numerical Solution of Initial Boundary Value Problems Involving Maxwell’s Equations in Isotropic Media. IEEE Transactions on Antennas and Propagation, Vol 14 Nº3, 302-307, Mayo 1966. [9] L. Talbi y G. Y. Delisle. “Finite Difference Time Domain Characterization of Indoor Radio Propagation”, Progress In Electromagnetics Research, PIER 12, 251-275, 1996.
[10] A. Thom and C.J. Apelt, Field Computations in Engineering and Physics. London: D. Van Nostrand, 1961.
[11] Sullivan, Dennis M, “Electromagnetic Simulation using the FDTD Method”, IEEE Press series on RF and Microwave Technology.
57
[12] Zhong J., Tapan K. S. y Binhong L. “Analysis of the Effects of Walls on Indoor Wave Propagation Using FDTD”.
ANEXOS
59
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % SIMULACION DEL CASO DE ESTUDIO % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; clc; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Definicion de constantes y frecuencia de operación %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c = 3.0E+8; % velocidad de la luz lo = 0.1765; % longitud de onda en metros F=1.7GHz spread = 2.5E-9; %ancho del pulso t0 = 5E-9; t = 0.0; % valor inicial del tiempo eta = 1/2; mu=4*pi*1E-7; %permeabilidad del espacio libre epszero = 8.85419e-12; % permitividad del espacio libre %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Particiones en tiempo y espacio %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% N_x = 4400; % particiones en x N_y = 1100; % particiones en y N_t = 1300; % pasos temporales delta_z = lo/10; % particion espacial delta_t = delta_z/(2*c); % particion temporal %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Definición de las paredes de la habitacion %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% param_diel_a = ones(N_x,N_y); param_diel_b = ones(N_x,N_y); % Variables del medio con pérdidas sigma = 10e-5; epsilon = 3.25; %Constante dieléctrica de la pared de concreto epsilon_vidrio = 3.95; %Constante dieléctrica de la ventana epsilon_madera = 1.94; %Constante dieléctrica de la puerta y pared de madera etaeps = eta/epsilon; etaeps_vidrio = eta/epsilon_vidrio; etaeps_madera = eta/epsilon_madera; param_perd = (delta_t*sigma)/((2*epszero*epsilon)); param_perd_vidrio = (delta_t*sigma)/((2*epszero*epsilon_vidrio)); param_perd_madera = (delta_t*sigma)/((2*epszero*epsilon_madera)); % Medio Dieléctrico for i = 1:N_x for j = 1:N_y param_diel_a(i,j) = 1; param_diel_b(i,j) = eta; end end
60
%Paredes Concreto %%%Parametro a param_diel_a(50:4390,50)= (1-param_perd)/(1+param_perd); %1 param_diel_a(62:4378,62)= (1-param_perd)/(1+param_perd); %2 param_diel_a(50,50:1087)= (1-param_perd)/(1+param_perd); %3 param_diel_a(62,62:1075)= (1-param_perd)/(1+param_perd); %4 param_diel_a(50:4390,1087)= (1-param_perd)/(1+param_perd); %5 param_diel_a(62:4378,1075)= (1-param_perd)/(1+param_perd); %6 param_diel_a(4390,50:1087)= (1-param_perd)/(1+param_perd); %7 param_diel_a(4378,62:1075)= (1-param_perd)/(1+param_perd); %8 param_diel_a(62:385,486)= (1-param_perd)/(1+param_perd); %9 param_diel_a(62:385,476)= (1-param_perd)/(1+param_perd); %10 param_diel_a(385,476:486)= (1-param_perd)/(1+param_perd); %11 param_diel_a(385,389:431)= (1-param_perd)/(1+param_perd); %12 param_diel_a(373,377:431)= (1-param_perd)/(1+param_perd); %13 param_diel_a(373:385,431)= (1-param_perd)/(1+param_perd); %14 param_diel_a(373:543,377)= (1-param_perd)/(1+param_perd); %15 param_diel_a(385:531,389)= (1-param_perd)/(1+param_perd); %16 param_diel_a(543,377:431)= (1-param_perd)/(1+param_perd); %17 param_diel_a(531,389:431)= (1-param_perd)/(1+param_perd); %18 param_diel_a(531:543,431)= (1-param_perd)/(1+param_perd); %19 param_diel_a(531:654,486)= (1-param_perd)/(1+param_perd); %20 param_diel_a(531:654,476)= (1-param_perd)/(1+param_perd); %21 param_diel_a(531,476:486)= (1-param_perd)/(1+param_perd); %22 param_diel_a(654,476:486)= (1-param_perd)/(1+param_perd); %23 param_diel_a(707:1109,486)= (1-param_perd)/(1+param_perd); %24
61
param_diel_a(707:1109,476)= (1-param_perd)/(1+param_perd); %25 param_diel_a(707,476:486)= (1-param_perd)/(1+param_perd); %26 param_diel_a(1109,476:486)= (1-param_perd)/(1+param_perd); %27 param_diel_a(1162:1187,486)= (1-param_perd)/(1+param_perd); %28 param_diel_a(1162:1187,476)= (1-param_perd)/(1+param_perd); %29 param_diel_a(1162,476:486)= (1-param_perd)/(1+param_perd); %30 param_diel_a(1187,476:486)= (1-param_perd)/(1+param_perd); %31 param_diel_a(1240:1642,486)= (1-param_perd)/(1+param_perd); %32 param_diel_a(1240:1642,476)= (1-param_perd)/(1+param_perd); %33 param_diel_a(1240,476:486)= (1-param_perd)/(1+param_perd); %34 param_diel_a(1642,476:486)= (1-param_perd)/(1+param_perd); %35 param_diel_a(1695:1836,486)= (1-param_perd)/(1+param_perd); %36 param_diel_a(1695:1824,476)= (1-param_perd)/(1+param_perd); %37 param_diel_a(1695,476:486)= (1-param_perd)/(1+param_perd); %38 param_diel_a(1836,476:486)= (1-param_perd)/(1+param_perd); %39 param_diel_a(1836,62:486)= (1-param_perd)/(1+param_perd); %40 param_diel_a(1824,62:476)= (1-param_perd)/(1+param_perd); %41 param_diel_a(2504,62:486)= (1-param_perd)/(1+param_perd); %42 param_diel_a(2516,62:476)= (1-param_perd)/(1+param_perd); %43 param_diel_a(2504:2657,486)= (1-param_perd)/(1+param_perd); %44 param_diel_a(2516:2657,476)= (1-param_perd)/(1+param_perd); %45 param_diel_a(2657,476:486)= (1-param_perd)/(1+param_perd); %46 param_diel_a(2710:3112,486)= (1-param_perd)/(1+param_perd); %47 param_diel_a(2710:3112,476)= (1-param_perd)/(1+param_perd); %48 param_diel_a(2710,476:486)= (1-param_perd)/(1+param_perd); %49
62
param_diel_a(3112,476:486)= (1-param_perd)/(1+param_perd); %50 param_diel_a(3165:3190,486)= (1-param_perd)/(1+param_perd); %51 param_diel_a(3165:3190,476)= (1-param_perd)/(1+param_perd); %52 param_diel_a(3165,476:486)= (1-param_perd)/(1+param_perd); %53 param_diel_a(3190,476:486)= (1-param_perd)/(1+param_perd); %54 param_diel_a(3243:3645,486)= (1-param_perd)/(1+param_perd); %55 param_diel_a(3243:3645,476)= (1-param_perd)/(1+param_perd); %56 param_diel_a(3243,476:486)= (1-param_perd)/(1+param_perd); %57 param_diel_a(3645,476:486)= (1-param_perd)/(1+param_perd); %58 param_diel_a(3698:3873,486)= (1-param_perd)/(1+param_perd); %59 param_diel_a(3698:3873,476)= (1-param_perd)/(1+param_perd); %60 param_diel_a(3698,476:486)= (1-param_perd)/(1+param_perd); %61 param_diel_a(3873,476:486)= (1-param_perd)/(1+param_perd); %62 param_diel_a(3873,383:431)= (1-param_perd)/(1+param_perd); %63 param_diel_a(3861,371:431)= (1-param_perd)/(1+param_perd); %64 param_diel_a(3861:3873,431)= (1-param_perd)/(1+param_perd); %65 param_diel_a(3861:3961,371)= (1-param_perd)/(1+param_perd); %66 param_diel_a(3873:3961,383)= (1-param_perd)/(1+param_perd); %67 param_diel_a(3961,371:383)= (1-param_perd)/(1+param_perd); %68 param_diel_a(4014,371:476)= (1-param_perd)/(1+param_perd); %69 param_diel_a(4026,371:476)= (1-param_perd)/(1+param_perd); %70 param_diel_a(4014:4026,371)= (1-param_perd)/(1+param_perd); %71 param_diel_a(4014:4207,486)= (1-param_perd)/(1+param_perd); %72 param_diel_a(4026:4195,476)= (1-param_perd)/(1+param_perd); %73 param_diel_a(4014,476:486)= (1-param_perd)/(1+param_perd); %74
63
param_diel_a(4207,62:486)= (1-param_perd)/(1+param_perd); %75 param_diel_a(4195,62:476)= (1-param_perd)/(1+param_perd); %76 param_diel_a(62:97,651)= (1-param_perd)/(1+param_perd); %77 param_diel_a(62:97,663)= (1-param_perd)/(1+param_perd); %78 param_diel_a(97,651:663)= (1-param_perd)/(1+param_perd); %79 param_diel_a(150:623,651)= (1-param_perd)/(1+param_perd); %80 param_diel_a(150:623,663)= (1-param_perd)/(1+param_perd); %81 param_diel_a(150,651:663)= (1-param_perd)/(1+param_perd); %82 param_diel_a(623,651:663)= (1-param_perd)/(1+param_perd); %83 param_diel_a(676:729,651)= (1-param_perd)/(1+param_perd); %84 param_diel_a(676:729,663)= (1-param_perd)/(1+param_perd); %85 param_diel_a(676,651:663)= (1-param_perd)/(1+param_perd); %86 param_diel_a(729,651:663)= (1-param_perd)/(1+param_perd); %87 param_diel_a(782:1202,651)= (1-param_perd)/(1+param_perd); %88 param_diel_a(782:1202,663)= (1-param_perd)/(1+param_perd); %89 param_diel_a(782,651:663)= (1-param_perd)/(1+param_perd); %90 param_diel_a(1202,651:663)= (1-param_perd)/(1+param_perd); %91 param_diel_a(1255:1308,651)= (1-param_perd)/(1+param_perd); %92 param_diel_a(1255:1308,663)= (1-param_perd)/(1+param_perd); %93 param_diel_a(1255,651:663)= (1-param_perd)/(1+param_perd); %94 param_diel_a(1308,651:663)= (1-param_perd)/(1+param_perd); %95 param_diel_a(1361:1764,651)= (1-param_perd)/(1+param_perd); %96 param_diel_a(1361:1764,663)= (1-param_perd)/(1+param_perd); %97 param_diel_a(1361,651:663)= (1-param_perd)/(1+param_perd); %98 param_diel_a(1764,651:663)= (1-param_perd)/(1+param_perd); %99
64
param_diel_a(1817:1922,651)= (1-param_perd)/(1+param_perd); %100 param_diel_a(1817:1922,663)= (1-param_perd)/(1+param_perd); %101 param_diel_a(1817,651:663)= (1-param_perd)/(1+param_perd); %102 param_diel_a(1922,651:663)= (1-param_perd)/(1+param_perd); %103 param_diel_a(1975:2373,651)= (1-param_perd)/(1+param_perd); %104 param_diel_a(1975:2373,663)= (1-param_perd)/(1+param_perd); %105 param_diel_a(1975,651:663)= (1-param_perd)/(1+param_perd); %106 param_diel_a(2373,651:663)= (1-param_perd)/(1+param_perd); %107 param_diel_a(2426:2531,651)= (1-param_perd)/(1+param_perd); %108 param_diel_a(2426:2531,663)= (1-param_perd)/(1+param_perd); %109 param_diel_a(2426,651:663)= (1-param_perd)/(1+param_perd); %110 param_diel_a(2531,651:663)= (1-param_perd)/(1+param_perd); %111 param_diel_a(2584:2987,651)= (1-param_perd)/(1+param_perd); %112 param_diel_a(2584:2987,663)= (1-param_perd)/(1+param_perd); %113 param_diel_a(2584,651:663)= (1-param_perd)/(1+param_perd); %114 param_diel_a(2987,651:663)= (1-param_perd)/(1+param_perd); %115 param_diel_a(3040:3093,651)= (1-param_perd)/(1+param_perd); %116 param_diel_a(3040:3093,663)= (1-param_perd)/(1+param_perd); %117 param_diel_a(3040,651:663)= (1-param_perd)/(1+param_perd); %118 param_diel_a(3093,651:663)= (1-param_perd)/(1+param_perd); %119 param_diel_a(3146:3649,651)= (1-param_perd)/(1+param_perd); %120 param_diel_a(3146:3649,663)= (1-param_perd)/(1+param_perd); %121 param_diel_a(3146,651:663)= (1-param_perd)/(1+param_perd); %122 param_diel_a(3649,651:663)= (1-param_perd)/(1+param_perd); %123 param_diel_a(3702:3755,651)= (1-param_perd)/(1+param_perd); %124
65
param_diel_a(3702:3755,663)= (1-param_perd)/(1+param_perd); %125 param_diel_a(3702,651:663)= (1-param_perd)/(1+param_perd); %126 param_diel_a(3755,651:663)= (1-param_perd)/(1+param_perd); %127 param_diel_a(3808:4290,651)= (1-param_perd)/(1+param_perd); %128 param_diel_a(3808:4290,663)= (1-param_perd)/(1+param_perd); %129 param_diel_a(3808,651:663)= (1-param_perd)/(1+param_perd); %130 param_diel_a(4290,651:663)= (1-param_perd)/(1+param_perd); %131 param_diel_a(4343:4378,651)= (1-param_perd)/(1+param_perd); %132 param_diel_a(4343:4378,663)= (1-param_perd)/(1+param_perd); %133 param_diel_a(4343,651:663)= (1-param_perd)/(1+param_perd); %134 param_diel_a(452,62:376)= (1-param_perd)/(1+param_perd); %135 param_diel_a(464,62:376)= (1-param_perd)/(1+param_perd); %136 param_diel_a(895,62:476)= (1-param_perd)/(1+param_perd); %137 param_diel_a(907,62:476)= (1-param_perd)/(1+param_perd); %138 param_diel_a(895,663:1075)= (1-param_perd)/(1+param_perd); %139 param_diel_a(907,663:1075)= (1-param_perd)/(1+param_perd); %140 param_diel_a(1817,663:1075)= (1-param_perd)/(1+param_perd); %141 param_diel_a(1829,663:1075)= (1-param_perd)/(1+param_perd); %142 param_diel_a(2519,663:1075)= (1-param_perd)/(1+param_perd); %143 param_diel_a(2531,663:1075)= (1-param_perd)/(1+param_perd); %144 param_diel_a(3450,62:476)= (1-param_perd)/(1+param_perd); %145 param_diel_a(3462,62:476)= (1-param_perd)/(1+param_perd); %146 param_diel_a(3450,663:1075)= (1-param_perd)/(1+param_perd); %147 param_diel_a(3462,663:1075)= (1-param_perd)/(1+param_perd); %148 param_diel_a(3724,62:476)= (1-param_perd)/(1+param_perd); %149
66
param_diel_a(3736,62:476)= (1-param_perd)/(1+param_perd); %150 param_diel_a(3724,663:1075)= (1-param_perd)/(1+param_perd); %151 param_diel_a(3736,663:1075)= (1-param_perd)/(1+param_perd); %152 param_diel_a(3949,62:371)= (1-param_perd)/(1+param_perd); %153 param_diel_a(3961,62:371)= (1-param_perd)/(1+param_perd); %154 %%%Parametro b param_diel_b(50:4390,50)= 0.5/(epsilon*(1+param_perd)); %1 param_diel_b(62:4378,62)= 0.5/(epsilon*(1+param_perd)); %2 param_diel_b(50,50:1087)= 0.5/(epsilon*(1+param_perd)); %3 param_diel_b(62,62:1075)= 0.5/(epsilon*(1+param_perd)); %4 param_diel_b(50:4390,1087)= 0.5/(epsilon*(1+param_perd)); %5 param_diel_b(62:4378,1075)= 0.5/(epsilon*(1+param_perd)); %6 param_diel_b(4390,50:1087)= 0.5/(epsilon*(1+param_perd)); %7 param_diel_b(4378,62:1075)= 0.5/(epsilon*(1+param_perd)); %8 param_diel_b(62:385,486)= 0.5/(epsilon*(1+param_perd)); %9 param_diel_b(62:385,476)= 0.5/(epsilon*(1+param_perd)); %10 param_diel_b(385,476:486)= 0.5/(epsilon*(1+param_perd)); %11 param_diel_b(385,389:431)= 0.5/(epsilon*(1+param_perd)); %12 param_diel_b(373,377:431)= 0.5/(epsilon*(1+param_perd)); %13 param_diel_b(373:385,431)= 0.5/(epsilon*(1+param_perd)); %14 param_diel_b(373:543,377)= 0.5/(epsilon*(1+param_perd)); %15 param_diel_b(385:531,389)= 0.5/(epsilon*(1+param_perd)); %16 param_diel_b(543,377:431)= 0.5/(epsilon*(1+param_perd)); %17 param_diel_b(531,389:431)= 0.5/(epsilon*(1+param_perd)); %18 param_diel_b(531:543,431)= 0.5/(epsilon*(1+param_perd)); %19 param_diel_b(531:654,486)= 0.5/(epsilon*(1+param_perd)); %20
67
param_diel_b(531:654,476)= 0.5/(epsilon*(1+param_perd)); %21 param_diel_b(531,476:486)= 0.5/(epsilon*(1+param_perd)); %22 param_diel_b(654,476:486)= 0.5/(epsilon*(1+param_perd)); %23 param_diel_b(707:1109,486)= 0.5/(epsilon*(1+param_perd)); %24 param_diel_b(707:1109,476)= 0.5/(epsilon*(1+param_perd)); %25 param_diel_b(707,476:486)= 0.5/(epsilon*(1+param_perd)); %26 param_diel_b(1109,476:486)= 0.5/(epsilon*(1+param_perd)); %27 param_diel_b(1162:1187,486)= 0.5/(epsilon*(1+param_perd)); %28 param_diel_b(1162:1187,476)= 0.5/(epsilon*(1+param_perd)); %29 param_diel_b(1162,476:486)= 0.5/(epsilon*(1+param_perd)); %30 param_diel_b(1187,476:486)= 0.5/(epsilon*(1+param_perd)); %31 param_diel_b(1240:1642,486)= 0.5/(epsilon*(1+param_perd)); %32 param_diel_b(1240:1642,476)= 0.5/(epsilon*(1+param_perd)); %33 param_diel_b(1240,476:486)= 0.5/(epsilon*(1+param_perd)); %34 param_diel_b(1642,476:486)= 0.5/(epsilon*(1+param_perd)); %35 param_diel_b(1695:1836,486)= 0.5/(epsilon*(1+param_perd)); %36 param_diel_b(1695:1824,476)= 0.5/(epsilon*(1+param_perd)); %37 param_diel_b(1695,476:486)= 0.5/(epsilon*(1+param_perd)); %38 param_diel_b(1836,476:486)= 0.5/(epsilon*(1+param_perd)); %39 param_diel_b(1836,62:486)= 0.5/(epsilon*(1+param_perd)); %40 param_diel_b(1824,62:476)= 0.5/(epsilon*(1+param_perd)); %41 param_diel_b(2504,62:486)= 0.5/(epsilon*(1+param_perd)); %42 param_diel_b(2516,62:476)= 0.5/(epsilon*(1+param_perd)); %43 param_diel_b(2504:2657,486)= 0.5/(epsilon*(1+param_perd)); %44 param_diel_b(2516:2657,476)= 0.5/(epsilon*(1+param_perd)); %45
68
param_diel_b(2657,476:486)= 0.5/(epsilon*(1+param_perd)); %46 param_diel_b(2710:3112,486)= 0.5/(epsilon*(1+param_perd)); %47 param_diel_b(2710:3112,476)= 0.5/(epsilon*(1+param_perd)); %48 param_diel_b(2710,476:486)= 0.5/(epsilon*(1+param_perd)); %49 param_diel_b(3112,476:486)= 0.5/(epsilon*(1+param_perd)); %50 param_diel_b(3165:3190,486)= 0.5/(epsilon*(1+param_perd)); %51 param_diel_b(3165:3190,476)= 0.5/(epsilon*(1+param_perd)); %52 param_diel_b(3165,476:486)= 0.5/(epsilon*(1+param_perd)); %53 param_diel_b(3190,476:486)= 0.5/(epsilon*(1+param_perd)); %54 param_diel_b(3243:3645,486)= 0.5/(epsilon*(1+param_perd)); %55 param_diel_b(3243:3645,476)= 0.5/(epsilon*(1+param_perd)); %56 param_diel_b(3243,476:486)= 0.5/(epsilon*(1+param_perd)); %57 param_diel_b(3645,476:486)= 0.5/(epsilon*(1+param_perd)); %58 param_diel_b(3698:3873,486)= 0.5/(epsilon*(1+param_perd)); %59 param_diel_b(3698:3873,476)= 0.5/(epsilon*(1+param_perd)); %60 param_diel_b(3698,476:486)= 0.5/(epsilon*(1+param_perd)); %61 param_diel_b(3873,476:486)= 0.5/(epsilon*(1+param_perd)); %62 param_diel_b(3873,383:431)= 0.5/(epsilon*(1+param_perd)); %63 param_diel_b(3861,371:431)= 0.5/(epsilon*(1+param_perd));%64 param_diel_b(3861:3873,431)= 0.5/(epsilon*(1+param_perd)); %65 param_diel_b(3861:3961,371)= 0.5/(epsilon*(1+param_perd)); %66 param_diel_b(3873:3961,383)= 0.5/(epsilon*(1+param_perd)); %67 param_diel_b(3961,371:383)= 0.5/(epsilon*(1+param_perd)); %68 param_diel_b(4014,371:476)= 0.5/(epsilon*(1+param_perd)); %69 param_diel_b(4026,371:476)= 0.5/(epsilon*(1+param_perd)); %70
69
param_diel_b(4014:4026,371)= 0.5/(epsilon*(1+param_perd)); %71 param_diel_b(4014:4207,486)= 0.5/(epsilon*(1+param_perd)); %72 param_diel_b(4026:4195,476)= 0.5/(epsilon*(1+param_perd)); %73 param_diel_b(4014,476:486)= 0.5/(epsilon*(1+param_perd)); %74 param_diel_b(4207,62:486)= 0.5/(epsilon*(1+param_perd)); %75 param_diel_b(4195,62:476)= 0.5/(epsilon*(1+param_perd)); %76 param_diel_b(62:97,651)= 0.5/(epsilon*(1+param_perd)); %77 param_diel_b(62:97,663)= 0.5/(epsilon*(1+param_perd)); %78 param_diel_b(97,651:663)= 0.5/(epsilon*(1+param_perd)); %79 param_diel_b(150:623,651)= 0.5/(epsilon*(1+param_perd)); %80 param_diel_b(150:623,663)= 0.5/(epsilon*(1+param_perd)); %81 param_diel_b(150,651:663)= 0.5/(epsilon*(1+param_perd)); %82 param_diel_b(623,651:663)= 0.5/(epsilon*(1+param_perd)); %83 param_diel_b(676:729,651)= 0.5/(epsilon*(1+param_perd)); %84 param_diel_b(676:729,663)= 0.5/(epsilon*(1+param_perd)); %85 param_diel_b(676,651:663)= 0.5/(epsilon*(1+param_perd)); %86 param_diel_b(729,651:663)= 0.5/(epsilon*(1+param_perd)); %87 param_diel_b(782:1202,651)= 0.5/(epsilon*(1+param_perd)); %88 param_diel_b(782:1202,663)= 0.5/(epsilon*(1+param_perd)); %89 param_diel_b(782,651:663)= 0.5/(epsilon*(1+param_perd)); %90 param_diel_b(1202,651:663)= 0.5/(epsilon*(1+param_perd)); %91 param_diel_b(1255:1308,651)= 0.5/(epsilon*(1+param_perd)); %92 param_diel_b(1255:1308,663)= 0.5/(epsilon*(1+param_perd)); %93 param_diel_b(1255,651:663)= 0.5/(epsilon*(1+param_perd)); %94 param_diel_b(1308,651:663)= 0.5/(epsilon*(1+param_perd)); %95
70
param_diel_b(1361:1764,651)= 0.5/(epsilon*(1+param_perd)); %96 param_diel_b(1361:1764,663)= 0.5/(epsilon*(1+param_perd)); %97 param_diel_b(1361,651:663)= 0.5/(epsilon*(1+param_perd)); %98 param_diel_b(1764,651:663)= 0.5/(epsilon*(1+param_perd)); %99 param_diel_b(1817:1922,651)= 0.5/(epsilon*(1+param_perd)); %100 param_diel_b(1817:1922,663)= 0.5/(epsilon*(1+param_perd)); %101 param_diel_b(1817,651:663)= 0.5/(epsilon*(1+param_perd)); %102 param_diel_b(1922,651:663)= 0.5/(epsilon*(1+param_perd)); %103 param_diel_b(1975:2373,651)= 0.5/(epsilon*(1+param_perd)); %104 param_diel_b(1975:2373,663)= 0.5/(epsilon*(1+param_perd)); %105 param_diel_b(1975,651:663)= 0.5/(epsilon*(1+param_perd)); %106 param_diel_b(2373,651:663)= 0.5/(epsilon*(1+param_perd)); %107 param_diel_b(2426:2531,651)= 0.5/(epsilon*(1+param_perd)); %108 param_diel_b(2426:2531,663)= 0.5/(epsilon*(1+param_perd)); %109 param_diel_b(2426,651:663)= 0.5/(epsilon*(1+param_perd)); %110 param_diel_b(2531,651:663)= 0.5/(epsilon*(1+param_perd)); %111 param_diel_b(2584:2987,651)= 0.5/(epsilon*(1+param_perd)); %112 param_diel_b(2584:2987,663)= 0.5/(epsilon*(1+param_perd)); %113 param_diel_b(2584,651:663)= 0.5/(epsilon*(1+param_perd)); %114 param_diel_b(2987,651:663)= 0.5/(epsilon*(1+param_perd)); %115 param_diel_b(3040:3093,651)= 0.5/(epsilon*(1+param_perd)); %116 param_diel_b(3040:3093,663)= 0.5/(epsilon*(1+param_perd)); %117 param_diel_b(3040,651:663)= 0.5/(epsilon*(1+param_perd)); %118 param_diel_b(3093,651:663)= 0.5/(epsilon*(1+param_perd)); %119 param_diel_b(3146:3649,651)= 0.5/(epsilon*(1+param_perd)); %120
71
param_diel_b(3146:3649,663)= 0.5/(epsilon*(1+param_perd)); %121 param_diel_b(3146,651:663)= 0.5/(epsilon*(1+param_perd)); %122 param_diel_b(3649,651:663)= 0.5/(epsilon*(1+param_perd)); %123 param_diel_b(3702:3755,651)= 0.5/(epsilon*(1+param_perd)); %124 param_diel_b(3702:3755,663)= 0.5/(epsilon*(1+param_perd)); %125 param_diel_b(3702,651:663)= 0.5/(epsilon*(1+param_perd)); %126 param_diel_b(3755,651:663)= 0.5/(epsilon*(1+param_perd)); %127 param_diel_b(3808:4290,651)= 0.5/(epsilon*(1+param_perd)); %128 param_diel_b(3808:4290,663)= 0.5/(epsilon*(1+param_perd)); %129 param_diel_b(3808,651:663)= 0.5/(epsilon*(1+param_perd)); %130 param_diel_b(4290,651:663)= 0.5/(epsilon*(1+param_perd)); %131 param_diel_b(4343:4378,651)= 0.5/(epsilon*(1+param_perd)); %132 param_diel_b(4343:4378,663)= 0.5/(epsilon*(1+param_perd)); %133 param_diel_b(4343,651:663)= 0.5/(epsilon*(1+param_perd)); %134 param_diel_b(452,62:376)= 0.5/(epsilon*(1+param_perd)); %135 param_diel_b(464,62:376)= 0.5/(epsilon*(1+param_perd)); %136 param_diel_b(895,62:476)= 0.5/(epsilon*(1+param_perd)); %137 param_diel_b(907,62:476)= 0.5/(epsilon*(1+param_perd)); %138 param_diel_b(895,663:1075)= 0.5/(epsilon*(1+param_perd)); %139 param_diel_b(907,663:1075)= 0.5/(epsilon*(1+param_perd)); %140 param_diel_b(1817,663:1075)= 0.5/(epsilon*(1+param_perd)); %141 param_diel_b(1829,663:1075)= 0.5/(epsilon*(1+param_perd)); %142 param_diel_b(2519,663:1075)= 0.5/(epsilon*(1+param_perd)); %143 param_diel_b(2531,663:1075)= 0.5/(epsilon*(1+param_perd)); %144 param_diel_b(3450,62:476)= 0.5/(epsilon*(1+param_perd)); %145
72
param_diel_b(3462,62:476)= 0.5/(epsilon*(1+param_perd)); %146 param_diel_b(3450,663:1075)= 0.5/(epsilon*(1+param_perd)); %147 param_diel_b(3462,663:1075)= 0.5/(epsilon*(1+param_perd)); %148 param_diel_b(3724,62:476)= 0.5/(epsilon*(1+param_perd)); %149 param_diel_b(3736,62:476)= 0.5/(epsilon*(1+param_perd)); %150 param_diel_b(3724,663:1075)= 0.5/(epsilon*(1+param_perd)); %151 param_diel_b(3736,663:1075)= 0.5/(epsilon*(1+param_perd)); %152 param_diel_b(3949,62:371)= 0.5/(epsilon*(1+param_perd)); %153 param_diel_b(3961,62:371)= 0.5/(epsilon*(1+param_perd)); %154 %Paredes Madera %%%Parametro a param_diel_a(378,431:476)=(1-param_perd_madera)/(1+param_perd_madera); %155 param_diel_a(381,431:476)=(1-param_perd_madera)/(1+param_perd_madera); %156 param_diel_a(535,431:476)=(1-param_perd_madera)/(1+param_perd_madera); %157 param_diel_a(538,431:476)=(1-param_perd_madera)/(1+param_perd_madera); %158 param_diel_a(654:707,479)=(1-param_perd_madera)/(1+param_perd_madera); %159 param_diel_a(654:707,482)=(1-param_perd_madera)/(1+param_perd_madera); %160 param_diel_a(1109:1162,479)=(1-param_perd_madera)/(1+param_perd_madera); %161 param_diel_a(1109:1162,482)=(1-param_perd_madera)/(1+param_perd_madera); %162 param_diel_a(1187:1240,479)=(1-param_perd_madera)/(1+param_perd_madera); %163 param_diel_a(1187:1240,482)=(1-param_perd_madera)/(1+param_perd_madera); %164 param_diel_a(1642:1695,479)=(1-param_perd_madera)/(1+param_perd_madera); %165 param_diel_a(1642:1695,482)=(1-param_perd_madera)/(1+param_perd_madera); %166 param_diel_a(2657:2710,479)=(1-param_perd_madera)/(1+param_perd_madera); %167 param_diel_a(2657:2710,482)=(1-param_perd_madera)/(1+param_perd_madera); %168 param_diel_a(3112:3165,479)=(1-param_perd_madera)/(1+param_perd_madera); %169
73
param_diel_a(3112:3165,482)=(1-param_perd_madera)/(1+param_perd_madera); %170 param_diel_a(3190:3243,479)=(1-param_perd_madera)/(1+param_perd_madera); %171 param_diel_a(3190:3243,482)=(1-param_perd_madera)/(1+param_perd_madera); %172 param_diel_a(3645:3698,479)=(1-param_perd_madera)/(1+param_perd_madera); %173 param_diel_a(3645:3698,482)=(1-param_perd_madera)/(1+param_perd_madera); %174 param_diel_a(3961:4014,375)=(1-param_perd_madera)/(1+param_perd_madera); %175 param_diel_a(3961:4014,378)=(1-param_perd_madera)/(1+param_perd_madera); %176 param_diel_a(3865,431:476)=(1-param_perd_madera)/(1+param_perd_madera); %177 param_diel_a(3868,431:476)=(1-param_perd_madera)/(1+param_perd_madera); %178 param_diel_a(97:150,655)=(1-param_perd_madera)/(1+param_perd_madera); %179 param_diel_a(97:150,658)=(1-param_perd_madera)/(1+param_perd_madera); %180 param_diel_a(623:676,655)=(1-param_perd_madera)/(1+param_perd_madera); %181 param_diel_a(623:676,658)=(1-param_perd_madera)/(1+param_perd_madera); %182 param_diel_a(729:782,655)=(1-param_perd_madera)/(1+param_perd_madera); %183 param_diel_a(729:782,658)=(1-param_perd_madera)/(1+param_perd_madera); %184 param_diel_a(1202:1255,655)=(1-param_perd_madera)/(1+param_perd_madera); %185 param_diel_a(1202:1255,658)=(1-param_perd_madera)/(1+param_perd_madera); %186 param_diel_a(1308:1361,655)=(1-param_perd_madera)/(1+param_perd_madera); %187 param_diel_a(1308:1361,658)=(1-param_perd_madera)/(1+param_perd_madera); %188 param_diel_a(1764:1817,655)=(1-param_perd_madera)/(1+param_perd_madera); %189 param_diel_a(1764:1817,658)=(1-param_perd_madera)/(1+param_perd_madera); %190 param_diel_a(1922:1975,655)=(1-param_perd_madera)/(1+param_perd_madera); %191 param_diel_a(1922:1975,658)=(1-param_perd_madera)/(1+param_perd_madera); %192 param_diel_a(2373:2426,655)=(1-param_perd_madera)/(1+param_perd_madera); %193 param_diel_a(2373:2426,658)=(1-param_perd_madera)/(1+param_perd_madera); %194
74
param_diel_a(2531:2584,655)=(1-param_perd_madera)/(1+param_perd_madera); %195 param_diel_a(2531:2584,658)=(1-param_perd_madera)/(1+param_perd_madera); %196 param_diel_a(2987:3040,655)=(1-param_perd_madera)/(1+param_perd_madera); %197 param_diel_a(2987:3040,658)=(1-param_perd_madera)/(1+param_perd_madera); %198 param_diel_a(3093:3146,655)=(1-param_perd_madera)/(1+param_perd_madera); %199 param_diel_a(3093:3146,658)=(1-param_perd_madera)/(1+param_perd_madera); %200 param_diel_a(3649:3702,655)=(1-param_perd_madera)/(1+param_perd_madera); %201 param_diel_a(3649:3702,658)=(1-param_perd_madera)/(1+param_perd_madera); %202 param_diel_a(3755:3808,655)=(1-param_perd_madera)/(1+param_perd_madera); %203 param_diel_a(3755:3808,658)=(1-param_perd_madera)/(1+param_perd_madera); %204 param_diel_a(4290:4343,655)=(1-param_perd_madera)/(1+param_perd_madera); %205 param_diel_a(4290:4343,658)=(1-param_perd_madera)/(1+param_perd_madera); %206 param_diel_a(211,60:476)=(1-param_perd_madera)/(1+param_perd_madera); %207 param_diel_a(214,60:476)=(1-param_perd_madera)/(1+param_perd_madera); %208 param_diel_a(620,60:476)=(1-param_perd_madera)/(1+param_perd_madera); %209 param_diel_a(623,60:476)=(1-param_perd_madera)/(1+param_perd_madera); %210 param_diel_a(1173,60:476)=(1-param_perd_madera)/(1+param_perd_madera); %211 param_diel_a(1176,60:476)=(1-param_perd_madera)/(1+param_perd_madera); %212 param_diel_a(1434,60:476)=(1-param_perd_madera)/(1+param_perd_madera); %213 param_diel_a(1437,60:476)=(1-param_perd_madera)/(1+param_perd_madera); %214 param_diel_a(2892,60:476)=(1-param_perd_madera)/(1+param_perd_madera); %215 param_diel_a(2895,60:476)=(1-param_perd_madera)/(1+param_perd_madera); %216 param_diel_a(3176,60:476)=(1-param_perd_madera)/(1+param_perd_madera); %217 param_diel_a(3179,60:476)=(1-param_perd_madera)/(1+param_perd_madera); %218 param_diel_a(322,663:1075)=(1-param_perd_madera)/(1+param_perd_madera); %219
75
param_diel_a(325,663:1075)=(1-param_perd_madera)/(1+param_perd_madera); %220 %%%Parametro b param_diel_b(378,431:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %155 param_diel_b(381,431:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %156 param_diel_b(535,431:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %157 param_diel_b(538,431:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %158 param_diel_b(654:707,479)= 0.5/(epsilon_madera*(1+param_perd_madera)); %159 param_diel_b(654:707,482)= 0.5/(epsilon_madera*(1+param_perd_madera)); %160 param_diel_b(1109:1162,479)= 0.5/(epsilon_madera*(1+param_perd_madera)); %161 param_diel_b(1109:1162,482)= 0.5/(epsilon_madera*(1+param_perd_madera)); %162 param_diel_b(1187:1240,479)= 0.5/(epsilon_madera*(1+param_perd_madera)); %163 param_diel_b(1187:1240,482)= 0.5/(epsilon_madera*(1+param_perd_madera)); %164 param_diel_b(1642:1695,479)= 0.5/(epsilon_madera*(1+param_perd_madera)); %165 param_diel_b(1642:1695,482)= 0.5/(epsilon_madera*(1+param_perd_madera)); %166 param_diel_b(2657:2710,479)= 0.5/(epsilon_madera*(1+param_perd_madera)); %167 param_diel_b(2657:2710,482)= 0.5/(epsilon_madera*(1+param_perd_madera)); %168 param_diel_b(3112:3165,479)= 0.5/(epsilon_madera*(1+param_perd_madera)); %169 param_diel_b(3112:3165,482)= 0.5/(epsilon_madera*(1+param_perd_madera)); %170 param_diel_b(3190:3243,479)= 0.5/(epsilon_madera*(1+param_perd_madera)); %171 param_diel_b(3190:3243,482)= 0.5/(epsilon_madera*(1+param_perd_madera)); %172 param_diel_b(3645:3698,479)= 0.5/(epsilon_madera*(1+param_perd_madera)); %173 param_diel_b(3645:3698,482)= 0.5/(epsilon_madera*(1+param_perd_madera)); %174 param_diel_b(3961:4014,375)= 0.5/(epsilon_madera*(1+param_perd_madera)); %175 param_diel_b(3961:4014,378)= 0.5/(epsilon_madera*(1+param_perd_madera)); %176 param_diel_b(3865,431:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %177 param_diel_b(3868,431:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %178
76
param_diel_b(97:150,655)= 0.5/(epsilon_madera*(1+param_perd_madera)); %179 param_diel_b(97:150,658)= 0.5/(epsilon_madera*(1+param_perd_madera)); %180 param_diel_b(623:676,655)= 0.5/(epsilon_madera*(1+param_perd_madera)); %181 param_diel_b(623:676,658)= 0.5/(epsilon_madera*(1+param_perd_madera)); %182 param_diel_b(729:782,655)= 0.5/(epsilon_madera*(1+param_perd_madera)); %183 param_diel_b(729:782,658)= 0.5/(epsilon_madera*(1+param_perd_madera)); %184 param_diel_b(1202:1255,655)= 0.5/(epsilon_madera*(1+param_perd_madera)); %185 param_diel_b(1202:1255,658)= 0.5/(epsilon_madera*(1+param_perd_madera)); %186 param_diel_b(1308:1361,655)= 0.5/(epsilon_madera*(1+param_perd_madera)); %187 param_diel_b(1308:1361,658)= 0.5/(epsilon_madera*(1+param_perd_madera)); %188 param_diel_b(1764:1817,655)= 0.5/(epsilon_madera*(1+param_perd_madera)); %189 param_diel_b(1764:1817,658)= 0.5/(epsilon_madera*(1+param_perd_madera)); %190 param_diel_b(1922:1975,655)= 0.5/(epsilon_madera*(1+param_perd_madera)); %191 param_diel_b(1922:1975,658)= 0.5/(epsilon_madera*(1+param_perd_madera)); %192 param_diel_b(2373:2426,655)= 0.5/(epsilon_madera*(1+param_perd_madera)); %193 param_diel_b(2373:2426,658)= 0.5/(epsilon_madera*(1+param_perd_madera)); %194 param_diel_b(2531:2584,655)= 0.5/(epsilon_madera*(1+param_perd_madera)); %195 param_diel_b(2531:2584,658)= 0.5/(epsilon_madera*(1+param_perd_madera)); %196 param_diel_b(2987:3040,655)= 0.5/(epsilon_madera*(1+param_perd_madera)); %197 param_diel_b(2987:3040,658)= 0.5/(epsilon_madera*(1+param_perd_madera)); %198 param_diel_b(3093:3146,655)= 0.5/(epsilon_madera*(1+param_perd_madera)); %199 param_diel_b(3093:3146,658)= 0.5/(epsilon_madera*(1+param_perd_madera)); %200 param_diel_b(3649:3702,655)= 0.5/(epsilon_madera*(1+param_perd_madera)); %201 param_diel_b(3649:3702,658)= 0.5/(epsilon_madera*(1+param_perd_madera)); %202 param_diel_b(3755:3808,655)= 0.5/(epsilon_madera*(1+param_perd_madera)); %203
77
param_diel_b(3755:3808,658)= 0.5/(epsilon_madera*(1+param_perd_madera)); %204 param_diel_b(4290:4343,655)= 0.5/(epsilon_madera*(1+param_perd_madera)); %205 param_diel_b(4290:4343,658)= 0.5/(epsilon_madera*(1+param_perd_madera)); %206 param_diel_b(211,60:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %207 param_diel_b(214,60:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %208 param_diel_b(620,60:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %209 param_diel_b(623,60:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %210 param_diel_b(1173,60:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %211 param_diel_b(1176,60:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %212 param_diel_b(1434,60:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %213 param_diel_b(1437,60:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %214 param_diel_b(2892,60:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %215 param_diel_b(2895,60:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %216 param_diel_b(3176,60:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %217 param_diel_b(3179,60:476)= 0.5/(epsilon_madera*(1+param_perd_madera)); %218 param_diel_b(322,663:1075)= 0.5/(epsilon_madera*(1+param_perd_madera)); %219 param_diel_b(325,663:1075)= 0.5/(epsilon_madera*(1+param_perd_madera)); %220 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Condiciones Iniciales de Frontera Absorbente %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Ez_low1 = 0; Ez_low2 = 0; Ez_high2 = 0; Ez_high1 = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Inicialización de los campos %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Ez = zeros(N_x,N_y); %Inicializa en cero el Campo Electrico Hx = Ez; %Inicializa en cero el Campo Magnetico Hy = Ez; %Inicializa en cero el Campo Magnetico %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Lazo para los calculos de E y H
78
for n = 1:N_t %Incrementa la onda electrica en el espacio for i = 2:N_x-1 for j = 2:N_y-1 Ez(i,j) = param_diel_a(i,j)*Ez(i,j) + param_diel_b(i,j)*( Hy(i,j)-Hy(i-1,j)-Hx(i,j)+Hx(i,j-1) ); end end %Fuente Ez(9,7)= Ez(9,7) + 50*exp (-0.5*((t-t0).*(t-t0)/spread^2)); %Condiciones de frontera Ez(1,1) = Ez_low2; Ez_low2 = Ez_low1; Ez_low1 = Ez(2,2); Ez(N_x,N_y) = Ez_high2; Ez_high2 = Ez_high1; Ez_high1 = Ez(N_x-1,N_y-1); %Incrementa la onda magnética en el espacio for i = 1:N_x-1 for j = 1:N_y-1 Hx(i,j) = Hx(i,j) + eta*( Ez(i,j)-Ez(i,j+1) ); Hy(i,j) = Hy(i,j) + eta*( Ez(i+1,j)-Ez(i,j) ); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Gráficas %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% zf = N_x*delta_z; % valor final de z z = delta_z*(1:N_x); % inicializa valores de z end