ON THE APPLICATION OF DIGITAL SIGNAL FOR STUDYING OPTICAL EMISSION SPECTROSCOPY

ELISABETH RESTREPO PARRA
M.Sc. Física, Universidad Nacional de Colombia – Sede Manizales, erestrepopa@unal.edu.co

JUAN CARLOS RIAÑO ROJAS
M.Sc. Matemáticas, Universidad Nacional de Colombia – Sede Manizales, jcrianoro@unal.edu.co

ANDRÉS FELIPE SERNA MORALES
Universidad Nacional de Colombia – Sede Manizales, afsernam@unal.edu.co

DAVID CÁRDENAS PEÑA ]]> Universidad Nacional de Colombia – Sede Manizales, dcardenasp@unal.edu.co

FLAVIO AUGUSTO PRIETO ÓRTIZ
Universidad Nacional de Colombia-Bogotá, faprietoo@unal.edu.co

Recibido para revisar junio 6 de 2008, aceptado septiembre 23 de 2008, versión final octubre 21 de 2008

RESUMEN: En este trabajo el primer paso consistió en llevar a cabo un suavizado y un filtrado de la señal, con el fin de eliminar el ruido y la información no relevante. La segunda parte consiste en la determinación automática de las longitudes de onda de los picos o máximos. Luego, se lleva a cabo la identificación elemental química de los picos empleando una base de datos especializada. Después de esto, el software permite la determinación de características necesarias para los cálculos como intensidad, FWHM (por su sigla en inglés Full Width at half Maximum, en español Ancho total a la Mitad del Máximo), los valores del continuo y el perfil de las líneas. Estos parámetros fueron empleados junto con métodos espectroscópicos tales como la gráfica de Boltzmann, relación entre intensidades de líneas de diferente grado de ionización y relación línea-continuo para determinar T_exc= 4778.32 K, n_e=2.54 10¹⁷ cm^-3 y T_e= 2668 K respectivamente. Estos valores son similares a los reportes, para este tipo de plasmas. El beneficio principal del software consintió en llevar a cabo todo el estudio del plasma en tiempo muy corto (del orden de segundos), mientras que su análisis manual puede tomar varios días.

PALABRAS CLAVE: Espectroscopía óptica de emisión, películas delgadas, plasma, extracción automática de caracteres.

ABSTRACT:In this work, an automatic extraction of characteristics from spectra taken during coatings production by plasma assisted techniques is proposed. The optical spectra also depend on these external parameters. The first step consists on carried out a filtering and smoothing of the signal, in order to eliminate the noise and non important information. The second step refers to determine automatically the wavelength of peaks. Then, the identification of chemical elements for each peak was done from a specialized data bases. After that, the software allows to determine peak characteristics such as intensity, FWHM (Full With at Half Maximum), continuous and the line profile. These parameters were employed together with spectroscopic methods such as Boltzmann plot, different ionization degree line to line ratio and line to continuous ratio in order to determine T_exc= 4778.32 K, n_e=2.54 10¹⁷ cm^-3 y T_e= 2668 K respectively. These values are similar to those reported in the literature, for this kind of plasmas. The principal advantage of this software was the possibility to analyze the plasma in very short time (in order of seconds), while the manual process can take several days.

KEYWORDS: Optical emission spectroscopy, thin films, plasma, automatic extraction of characteristics.

1. INTRODUCCIÓN

La espectroscopía óptica de emisión- OES (por su sigla en inglés Optical Emission Spectroscopy) ha emergido como una técnica poderosa de diagnóstico para el monitoreo en tiempo real de procesos de plasmas [1-3]. Este método provee un análisis in-situ no invasivo de plasmas sin interferir con el proceso [4-6]. La caracterización por OES ha sido ampliamente empleada en procesamiento de materiales por medio de diferentes técnicas asistidas por plasma, especialmente en procesos de crecimientos de recubrimientos. Yong M. Kim et al. [7] realizaron un estudio espacial de un plasma de Ar-N₂ producido por la técnica de magnetron sputtering desbalanceado. En este caso se consideraron las líneas de emisión producidas por el cátodo de titanio y por la mezcla de gases. G. Zambrano et al. [8] llevaron a cabo una investigación por medio de OES a una descarga de magnetron sputtering r.f. empleada para producir multicapas nanométricas de carburo de tungsteno (WC) y Diamod-like-carbon (DLC). Por medio de estos espectros se obtuvo información de las especies presentes en el plasma y las densidades de estas especies que se pudieron relacionar con la estructura y la composición de las capas depositadas. S. Peter et al [9] estudiaron procesos CVD empleados para la deposición de TiN por medio de una descarga DC con mezcla de TiCl₄-H₂-N₂-Ar. El análisis de los espectros permitió determinar las especies reactivas en la vecindad del crecimiento de la película. La evolución temporal de la descarga se atribuyó al movimiento lento del frente de ionización desde el cátodo al ánodo. A. Raveh et al [10] monitorearon el proceso de deposición de películas graduales de TiAlN por la técnica de sputtering reactivo asistido por resonancia ciclotrónica electrónica (ECR). Los resultados de OES indicaron que cuando se adiciona micro-ondas a un plasma r.f., este produce cambios en la concentración de Ti y Al en la fase gaseosa. Existen muchos otros trabajos en donde se puede observar la gran utilidad de la técnica OES como método de caracterización de plasmas. Sin embargo, aunque el montaje experimental no es muy complicado, un problema frecuente en esta técnica es la identificación de las líneas debido a diferentes ambigüedades, haciéndose dispendiosa la determinación de las características. Además, el procesamiento manual de los espectros es largo y algunas veces no muy confiable. Es por esto que se ha buscado la manera de optimizar este proceso. Algunos autores han aplicado metodologías y técnicas de modelamiento y automatización del estudio de espectros de plasmas. Roawen Chen [11] empleó varias técnicas de modelamiento tales como análisis de componentes principales multivariantes (por sus siglas en inglés MPCA) y mínimos cuadrados parciales (por su sigla en inglés PMC), para relacionar características de los espectros con el desempeño de un plasma etching (plasma de ataque), empleado en la fabricación de dispositivos de circuitos integrados submicrométricos. M. Gatto et al [12] desarrollaron un algoritmo para estudiar los espectros de emisión generados por los gases emitidos por un transistor bipolar de hetero-union (por sus siglas en inglés heterojunction bipolar transistor HBT). Las líneas de Galio permitieron determinar cuándo se alcanza la capa de GaInAs.

Debido a la necesidad que se tiene de optimizar el estudio de los espectros ópticos de emisión, en este trabajo se presenta la implementación de un sistema automático empleado en su procesamiento. El software fue desarrollado empleando herramientas del lenguaje builder c++.

2. MONTAJE EXPERIMENTAL

Los espectros fueron tomados empleando un espectrómetro de alta resolución HR4000 el cual incluye una rejilla de difracción HC1 de 1200 líneas/mm y una ranura de entrada de 20 mm. Los espectros son capturados en un rango entre 100 nm y 1200 nm, con una resolución de 1 nm. El plasma fue producido por una descarga r.f. en una mezcla de argón y nitrógeno, con una presión de 10^-3 mbar y un flujo de 2.8 cm³-min^-1.

La potencia de la descarga fue de 200 W y el material del blanco es TiAl (Titanio-Aluminio).

3. FUNDAMENTOS TEÓRICOS

3.1 Ecuación de Boltzmann para determinar T_{exc ]]>} La densidad de partículas que se encuentran en el estado i dentro del plasma se denomina . Se puede determinar la densidad total de partículas como [13].

(1)

Por ejemplo, si se considera i=2, se tiene que

(2)

g₂ y E₂ son el peso estadístico y la energía del nivel 2 respectivamente, k es la constante de Boltzmann, T_exc es la temperatura promedio de los electrones en el estado de energía 2, también denominada temperatura de excitación y es la función de partición dada por [14]

(3)

Cuya sumatoria se realiza sobre todos los niveles de energía i. Esta función de partición actúa como factor de normalización que asegura el cumplimiento de (2).

Por otro lado, la intensidad de una línea espectral está dada por [15]:

(4)

Siendo A₂₁ el coeficiente de Einstein, h la constante de Planck y la frecuencia de la línea , con c la velocidad de la luz y la longitud de onda. Reemplazando (2) en la ecuación (4), en función de la longitud de onda se tiene, la intensidad de una línea espectral se puede escribir de la forma:

]]> (5)

N, h, c y U son los mismos para todas las líneas atómicas del espectro, entonces:

Si se grafica

(6)

]]> Se obtiene una curva cuya pendiente es

(7)

Este método es ampliamente empleado en la literatura [16-18].

3.2 Relación línea-continuo para determinar T_e
Partiendo de la ecuación (4) para la intensidad de una línea y de la distribución de Saha [19]:

(8)

Siendo E_i la energía de ionización, el potencial mínimo de ionización, n, n_ion y n_e las densidades de las partículas neutras, los iones y los electrones respectivamente, m_e la masa del electrón y T_ion, la temperatura de los iones del plasma. Tomando N=n, despejando n/U de (8) y reemplazando en (5), se tiene que la intensidad de la línea se puede escribir como

(9)

El continuo de un espectro está definido por la expresión [20].

(10)

En condiciones de equilibrio termodinámico local (LTE), se puede considerar que [21]. Dividiendo (9) entre (10) se tiene

(11)

Con el fin de simplificar

(12)

(13)

(14)

(15)

(16)

(17)

Se grafican las funciones y de (16) y (17), luego se determina T_e a partir de su intercepto. Este método se ha reportado en la literatura por Sola et al. [22]

3.3 Relación entre intensidades de líneas de diferente grado de ionización para el cálculo de n_e
Para obtener la densidad electrónica se emplea la relación entre dos líneas que pertenecen al mismo elemento pero que tienen diferente grado de ionización. Sean dos líneas , de un elemento neutro, debida a la transición entre los niveles de energía q y p e del mismo elemento ionizado una vez, debida a la transición entre los niveles de energía r y s. La relación entre sus intensidades, de acuerdo con (5), está dada por:

(18)

Siendo la densidad de partículas neutras y la densidad de partículas ionizadas. La ecuación en función de la longitud de onda está dada por:

(19)

Lo que es lo mismo

(20)

De la ecuación de Saha se obtiene [19]:

(21)

Sustituyendo (21) en (20)

(22)

Despejando n_e se tiene que:

(23)

Esta es la expresión final empleada para determinar la densidad electrónica n_e. Esta técnica ha sido empleada por diversos autores [23].

4. IMPLEMENTACIÓN DEL ALGORITMO

En la figura 1 se presenta el diagrama de bloques del algoritmo implementado. El sistema automático tiene las siguientes etapas:

Figura 1. Diagrama de bloques empleado para la implementación del software automatizado para el procesamiento de espectros
Figure 1. Blocks diagram employed in the implementation of automated software for the spectra processing

Etapa 1: Se filtra la señal y al mismo tiempo se suaviza. Esta etapa es necesaria para reducir el ruido de la señal, que en gran parte se debe a los componentes electrónicos del sistema de medición [24]. Además, se pretende reducir el impacto de los datos irrelevantes o de los picos más pequeños, los cuales no son muy útiles a la hora de realizar el cálculo de los parámetros, pues la relación señal/ruido para estos casos no es pequeña. Al filtrar la señal, se hace menos engorrosa la etapa de identificación de picos, pues solo se tienen en cuenta aquellas líneas con una intensidad lo suficientemente alta como para que sean útiles a la hora de procesar la señal. La escogencia de los filtros se realiza de tal manera que no afecte los picos altos. Con este objetivo se emplean dos posibles filtros. Estos fueron escogidos debido a que son fáciles de implementar y además, este tipo de espectros no requieren de filtros más complejos. El primero se basa en un filtro Gaussiano, puesto que en condiciones ideales las señales gaussianas poseen propiedades importantes como la concentración de los datos alrededor de la media [25]. En este caso se lleva a cabo una convolución espectral, con una curva de media , una varianza dadas por [26].

(24)

(25)

Siendo cada punto del eje vertical y N el número de puntos. El tamaño de ventaneo es definido por el usuario. El segundo es un filtro mediana, que es un filtro no lineal [27], el cual puede tener la ventaja sobre el filtro gaussiano, de no correr los picos máximos con respecto a su posición inicial. La mediana se calcula a partir de una ventana del espectro de tamaño impar que puede ser escogida por el usuario. Este valor de la mediana es asignado como el nuevo valor del punto central de la ventana.

Etapa 2: Con el fin de determinar los máximos del espectro, se establece un umbral usando la hipótesis de normalidad gausiana del espectro, como se muestra en la figura 2. Aquí se representa una distribución gaussiana normal. Los datos que se encuentren por encima de , pertenecen a la cola derecha de la distribución y son los que se ubican más lejos de la media, por lo tanto tendrán las intensidades más altas. Para determinar el umbral se implementaron dos procedimientos:

El primero consiste en realizar un proceso de umbral medio [28] (PMT – por sus siglas en inglés process mean threshold), que se basa en calcular la magnitud de la media de todo el espectro y su desviación estándar , como se muestra en la figura 3(a). La segunda posibilidad consiste en un umbral medio local, que usa una ventana media móvil con un tamaño de 101 píxeles y su desviación estándar. Este genera un umbral que puede adaptarse al fondo irregular de los espectros (ver figura 3(b)). El usuario tiene la posibilidad de adicionar o eliminar números decimales de la desviación estándar para modificar el umbral.

Figura 2. Distribución normal de datos
Figure 2. Normal data distribution

Figura 3. (a) Umbral medio del proceso (b) Umbral medio local
Figure 3. (a) Process mean threshold (b) Local mean threshold

Etapa 3: Después del filtrado y de la ubicación del umbral, se procede a la identificación de los máximos del espectro que están por encima del umbral, es decir por encima del 84 por ciento de la altura de los datos (Ver figura 1), que es la probabilidad de estar en la cola derecha .

Además, se identifican los mínimos para cada máximo, uno izquierdo y otro derecho. Con el fin de separar los máximos de los mínimos, se obtienen los puntos cuyos vecinos tiene valores menores, empleando la expresión .

Después de esto, se crean los vectores que contienen los puntos máximos y mínimos, el izquierdo y el derecho más próximos, adjuntándole a cada máximo sus dos mininos correspondientes, como se muestra en la figura 4.

Figura 4. Identificación de máximos del espectro
Figure 4. Identification of spectrum maxima

Etapa 4: Se determina la intensidad de cada pico calculando el área bajo la curva de la línea espectral, empleando la expresión

(26)

Etapa 5: Para realizar la identificación elemental, se emplea la base de datos de NIST (Nacional Institute of Standard System Technology), NIST Atomic Spectra Database [29] y herramientas SQL (Structured Query Language). Se comparan los máximos obtenidos en el espectro con las longitudes de onda de la base de datos, con una desviación de 1 nm, ya que esta es la resolución del equipo empleado para la toma del espectro. Luego se seleccionan las líneas que posean la mayor intensidad, que está directamente relacionada con la probabilidad de que ocurra la transición entre los dos niveles de energía correspondientes, y que en la base de datos se encuentra en la columna de intensidades relativas.

Etapa 6: Se calcula la temperatura de excitación ( ). Con este fin, se adiciona una desviación estándar de 1.5 a la mediana, esto es, datos por encima de un 93 porciento, escogiéndose los picos más altos. Las líneas se separan dependiendo de su grado de ionización. Luego, se escoge el elemento con mayor cantidad de líneas, que en este caso es el Ar I (Argón neutro). Después de esto, se separan las líneas generadas por transiciones de energía cuyos niveles inferiores sean similares. Esta es una importante condición física, para que se puedan tener valores muy cercanos a los reales. Con este objetivo se define una posible diferencia entre niveles de energía inferiores 2000 cm^-1 (aproximadamente una diferencia de 10% entre ellos). Finalmente, se procesan lo datos empleando (6) y (7)

Etapa 8: Para el cálculo de se escoge una línea de Ar I bien aislada y se determina el valor del continuo (background) debajo de la línea, por medio de promediar sus dos valores mínimos. Para este caso se aplica las ecuaciones (16) y (17).

Etapa 9: La densidad electrónica del plasma es obtenida empleando el método de relación entre intensidades de líneas de diferente grado de ionización, empleando la relación (23), obtenida en la sección 3.3.

4. RESULTADOS Y ANÁLISIS

En la figura 5 se muestra la gráfica de Boltzmann obtenida por medio del software diseñado. El valor obtenido fue T_exc= 4778.32 K. Este resultado es similar a aquellos reportados para plasmas generados con la técnica de magnetron sputtering [30, 31]. Esto confirma el buen desempeño del software.

Figura 5. Cálculo de T_exc empleando la gráfica de Boltzmann
Figure 5. Determination of T_exc employing Boltzmann plot

Posteriormente, se determinó la temperatura electrónica T_e empleando la relación de línea continuo [22]. En este caso se elige una línea bien aislada (es decir que se pueda identificar totalmente su máximo y sus dos mínimos). La intensidad de la línea es el área bajo la curva. El continuo se determina midiendo el fondo ubicado directamente debajo del máximo como se muestra en la figura 6.

Figura 6 Línea empleada para el cálculo de T_e
Figure 6 Line employed for determining T_e

En la tabla 1 se presentan las variables empleadas y sus correspondientes valores. Estos valores se reemplazan en la ecuación (11). Debido a la dificultad para resolver analíticamente esta ecuación (4), esta se divide en dos partes, y . En la figura 7 se presenta la gráfica de estas dos funciones y su intercepto, cuyo valor en el eje de las abscisas corresponde a K que es similar al reportado en la literatura [36].

Tabla 1. Parámetros empleados en el cálculo de T_e
Table 1. Parameters used for determining T_e

Figura 7 Gráfica de las ecuaciones y , en donde se puede observar el intercepto T_e
Figure 7 Plot of equations y , where it is possible to observe the intercept T_e

5. CONCLUSIONES

Se implementó un software para el análisis de espectros ópticos de emisión capturados en un proceso de producción de recubrimientos por técnicas asistidas por plasma. Con este fin se emplearon técnicas de procesamiento digital de señales para realizar tareas de filtrado, suavizado, identificación de máximos y mínimos y comparación con bases de datos.

Por medio de este software se redujo el tiempo de análisis de los espectros, el cual normalmente se realiza de una forma manual, que de acuerdo a la pericia de la persona y a la complejidad del espectro, puede tomar desde un día hasta una semana.

Empleando métodos propios del análisis de espectros como son la gráfica de Boltzmann, relación entre líneas de diferente grado de ionización, relación entre intensidades de líneas y su continuo, se obtuvieron propiedades físicas de los plasmas producidos, tales como , , , que arrojaron valores similares a los reportados en la literatura.

6. AGRADECIMIENTOS

Los autores agradecen A la Universidad Nacional de Colombia Sede Manizales por su apoyo para la realización de este artículo, a través de la División para el apoyo a la investigación, (DIMA) bajo el proyecto 10719.

[1] DASKALOVA N., L. ALESKA l.; GENTSHEVA G.; ELICHKOV S. V., Spectrochimica Acta Part B 57 (2002) 755–768.         [ Links ]
[2] ZENG X., MAO S. S.;LIU C.; MAO X.; GREIF R.; RUSO R. E., Spectrochimica Acta Part B 58 (2003) 867–877         [ Links ]
[3] DENAUD J.; HOWES A.; POUSSEL E.; MERMET J.-M., Spectrochimica Acta Part B 56 (2001) 101-112         [ Links ]
[4] TRISTAND P.; DING Z.; TRANG VINH Q.B.; HIDALGO H.; JAUBERTEAU J. L.; DESMAISONA J.; DONG C., Thin Solid Films 390 (2001) 51-58         [ Links ]
[5] VANDEVELDE T; NESLADEK M.; QUAEYHAEGENS C.; STALS L., Thin Solid Films 308–309 (1997) 154–158         [ Links ]
[6] CALI C.; MACALUSO R.; MOSCA M., Spectrochimica Acta Part B 56 (2001)743-751         [ Links ]
[7] KIM Y. M., JUNG M. J.; OH S. G.; HAN J. G., Thin Solid Films 475 (2005) 91         [ Links ]
[8] ZAMBRANO G.; RIASCOS H.; PRIETO P.; RESTREPO E.; DEVIA A.; RINCON C.; Surface and Coating Tecnology 172 (2003) 144         [ Links ]
[9] PETER S.; RICHTER F.; TABERSKY R.; KÖNING, Thin Solid Films, 377-378 (2000) 430-435.         [ Links ]
[10] RAVEH A.; WEISS M.; SCHNECK R., Surface and Coating Technology, 111 (1999) 263-268.         [ Links ]
[11] CHEN R.; HUANG H.; SPANOS C. J., Plasma etch modeling using optical emission spectroscopy, University of California at Berkeley, Berkeley, California 94 720-1770         [ Links ]
[12] GATTO M., Journal of Vacuum Science and Technology A, 14 (1996) 1901-1906.         [ Links ]
[13] SALOUM S.; NADDAF M., Vacuum 82 (2008) 66–71         [ Links ]
[14] WENGER Ch. ; CHAMPION J. P.; BOUDON V; Journal of Quantitative Spectroscopy Radiative Transfer 109 (2008) 2697–2706         [ Links ]
[15] LARIJANI M. M.; LE NORMAND F.; CRÉGUT O., Applied Surface Science 253 (2007) 4051–4059         [ Links ]
[16] AGUILERA J.A.; ARAGÓN C., Spectrochimica Acta Part B 59 (2004) 1861–1876.         [ Links ]
[17] STAVROPOULOS P.; PALAGAS C.; ANGELOPOLOS G. N.; PAPAMANTELLOS D. N.; COURIS S., Spectrochimica Acta Part B 59 (2004) 1885–1892.         [ Links ]
[18] LE DROGOFF B.; MARGOT J.; CHAKER M.; SABSABI M.; BARTHÉLEMY O.; JOHNSTON T.W.; LAVILLE S.; VIDAL F.; VON KAENEL Y., Spectrochimica Acta Part B 56 (2001) 987–1002.         [ Links ]
[19] JIN F.; ZENG J.; YUAN J., Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2707– 2714         [ Links ]
[20] RESTREPO E.; DEVIA A., Journal of Vacuum Science Technology A. 22. (2) (2004) 377-382.         [ Links ]
[21] GARCÍA L. A.; PULZARA A. O.; DEVIA A.; RESTREPO E., Journal of Vacuum Science and Technology A 23 (2005) 551-553.         [ Links ]
[22] SOLA A.; CALZADA M. D.; GAMERO A., Journal of Physics D: Applied Physics 28 (1995) 1099- 1110.         [ Links ]
[23] IVKOVIC M.; JOVICÉVIC S.; KONJEVIC N., Spectrochimica Acta Part B 59 (2004) 591–605         [ Links ]
[24] VASEGHI S., Advanced Digital Signal Processing and Noise Reduction Jhon Wiley and Son, 2000 Bafins Lane (Ingland).         [ Links ]
[25] SPANN M.; NIEMINEM A.; Pattern Recognition Letters, 8 (1988) 251-255         [ Links ]
[26] PASTOR D., Computational Statistics Data Analysis, 52 (2008) 3167-3186.         [ Links ]
[27] LUKAC R.; SMOLKA B.; PLATANIOTIS K. N., Signal Processing 87 (2007) 2085-2099         [ Links ]
[28] TANG H.; QI L.-H.; LI H.-J., Carbon, 45 (2007) 2323-2324.         [ Links ]
[29] Atomic Spectra Databases, Version 2.0, National Institute of Standards Technology (NIST), USA , 1999, in /http://www.nist.gov.         [ Links ]
[30] RESTREPO E.; DEVIA A., Journal of Vacuum Science and Technology A. 22. (2) (2004) 377-382.         [ Links ]
[31] NOZAKI T.; SASAKI K.; OGINO T.; ASAHI D.; OKAZAKI K., Journal of Thermal Science and Technology 2 (2007) 192-199.         [ Links ]
[32] Kurucz RL. Atomic spectral line database from CD-ROM 23. Harvard-Smithsonian Center for Astrophysics, 1995, /http://cfa-www.harvard.edu.         [ Links ]
[33] Environmental chemistry and Hazardous materials information and resources, http://environmentalchemistry.com/yogi/periodic.[citado noviembre 15 de 2008].         [ Links ]
[34] FISHER V.; BERNSHTON V.; GOLTEL H.; MARON Y., Physical Review A, 53 (1996) 2425.         [ Links ]
[35] GRIEM H., Plasma Spectroscopy. New York: MacGraw-Hill; 1964.         [ Links ]
[36] RESTREPO E.; DEVIA A.; COGOLLO R.; MARIÑO A.; SANCHEZ H.; Revista Colombiana de Física 34 (2002) 488-492.         [ Links ] ]]> 1 2002 57 755-768 2 2003 58 867-877 3 2001 56 101-112 4 2001 390 51-58 5 1997 308-309 154-158 6 2001 56 743-751 7 2005 475 91 8 2003 172 144 9 2000 377-378 430-435 10 1999 111 263-268 11 94 720-1770 12 1996 14 1901-1906 13 2008 82 66-71 14 2008 109 2697-2706 15 2007 253 4051-4059 16 2004 59 1861-1876 17 2004 59 1885-1892 18 2001 56 987-1002 19 2008 109 2707- 2714 20 2004 22 377-382 21 2005 23 551-553 22 1995 28 1099- 1110 23 2004 59 591-605 24 2000 25 1988 8 251-255 26 2008 52 3167-3186 27 2007 87 2085-2099 28 2007 45 2323-2324 29 1999 30 2004 22 377-382 31 2007 2 192-199 32 1995 33 34 1996 53 2425 35 1964 36 2002 34 488-492