Cubierta de configuración dinámica para iluminación solar en hangar

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo de Fin de Grado Ingeniería Aeroespacial Cubierta de configuración dinámica para iluminación solar en hangar Autor: Jesús María Alberto Romero Tutor: Johan Wideberg Departamento de Ingeniería y Ciencia de los Materiales y del Transporte Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 Trabajo de Fin de Grado Ingeniería Aeroespacial Cubierta de configuración dinámica para iluminación solar en hangar Autor: Jesús María Alberto Romero Tutor: Johan Wideberg Profesor Titular Departamento de Ingeniería y Ciencia de los Materiales y del Transporte Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 Trabajo de Fin de Grado: Cubierta de configuración dinámica para iluminación solar en hangar Autor: Jesús María Alberto Romero Tutor: Johan Wideberg El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos El día que decidí estudiar esta preciosa carrera no era totalmente consciente de las implicaciones que eso conllevaría a nivel personal y profesional. A nivel personal, me obligaría a tomar decisiones complicadas en pro de la responsabilidad, hipotecando por ello los mejores años de la juventud de una persona. Me llevaría a límites de desesperación insospechados, porque el esfuerzo nunca estaba ligado al éxito. De hecho, el éxito siempre llegaba después de muchos fracasos. Poco a poco el cuerpo y la mente se empiezan a adaptar a un nivel de exigencia psicológica jamás experimentado, y empieza a entender que la clave del éxito no es solo el esfuerzo, sino que es la perseverancia y la constancia. La mente se va moldeando y acostumbrando a trabajar en situaciones muy complicadas, permitiendo cada vez ver las situaciones con más perspectiva, calma y análisis. Nada de lo que parece ser tan terrible, lo acaba siendo. Quizás eso sea lo que me haya enseñado la carrera, una forma de ver la vida que sólo otro compañero que haya pasado por la misma escuela podrá entender. Toda mi dedicatoria de agradecimiento solo puede ir dirigida a mi familia. Mis padres, junto con mis dos hermanas son los que han hecho posible que hoy haya podido finalizar este camino. Sin su apoyo incondicional permanente, sin su confianza plena hacia mi, y sin su amor infinito, no habría tenido el suficiente valor para haber podido superar todos los escollos que me he ido encontrando. A mis amigos, los que he ido creando en Sevilla y los que dejé en mi tierra natal en Manilva, los cuales no me han soltado en todo el camino, a pesar de mis dificultades. A Carlos B. y María S., por darme la oportunidad en su día de empezar una bonita carrera profesional aeronáutica. Antes jefes, ahora amigos. A Johan Wideberg, por depositar en mí toda su confianza desde el primer día que le propuse el proyecto, y que me ha ayudado todo lo que se puede uno imaginar. I Resumen En este proyecto se tratará la forma de permitir la entrada de luz natural en un hangar que se encuentra totalmente tapado en su cara interior. La luz sólo podrá entrar cumpliendo ciertas restricciones. Por eso, primero se determinará la posición del Sol, para luego conocer la sombra que arroja el edificio y por tanto, cuál es la superficie del techo que ha de ser tapada para garantizar sombra en todo su interior. III XNotación Cmax Longitud total del hangar en el punto x0 0,y0 0 según la dirección ,Az(t) 1 Introducción La sociedad en la que vivimos está generación tras generación, cada vez más concienciada por el impacto ambiental, la limitación de los recursos naturales y el calentamiento global. Es por eso que cada vez más personas, y también más empresas se comprometen en hacer un uso responsable de los recursos, y optimizarlos para evitar pérdidas y derroches. Es también desde la administración pública, no solo nacional sino a nivel global, quienes realizan campañas de concienciación para dar valor a los recursos que disponemos, o incluso permiten mediante ayudas directas la adquisición de elementos respetuosos con el medio ambiente, como pueden ser las ayudas para la compra de vehículos eléctricos o las ayudas para la adquisición de placas solares y fotovoltaicas. En el caso que nos ocupa, trataremos con nuestro proyecto estudiar la forma de mejorar la iluminación interior de un hangar haciendo un uso óptimo de la luz natural, basado en ciertas restricciones que se deben de cumplir. Actualmente en el edifico en el que se realizará el estudio no entra absolutamente nada de luz natural. Eso quiere decir que toda la iluminación que se usa en su interior es artificial, y que teniendo en cuenta las dimensiones que posee, se puede presuponer un altísimo coste energético. Por tanto, cualquier adición de luz natural que entre en el hangar de forma controlada, será un hecho beneficioso en términos ambientales. Por otro lado, la implementación de un sistema para la entrada de luz natural no sólo tendrá beneficios ecológicos, sino también de salud laboral. Aunque la iluminación artificial instalada satisface los requisitos mínimos de iluminación en espacios cerrados, según indica la norma de Salud Laboral del centro, es cierto que la actividad laboral prolongada en un espacio cerrado sin consciencia temporal a partir de la iluminación ambiental, promueve estados de ánimo de apatía, deficiencia visual y fatiga muscular [ 7 ]. Existen numerosos estudios publicados acerca de cómo mejorar las condiciones de iluminación en interiores, a partir de luz natural. Si bien es cierto que aunque nuestro proyecto tendrá algunos puntos en común con algunas referencias que se comentan a continuación, ninguno de los estudios a los que se ha podido tener acceso tratan el tema de la misma manera. En primera lugar, en la referencia [ 2 ] se estudia el diseño de una fachada dinámica de dos grados de libertad. Este proyecto consiste en construir una estructura por delante de la propia fachada del edificio, la cual se orienta de forma autónoma para garantizar las condiciones de iluminación óptimas en el interior del edificio. Para ello, usa diversos sensores, entre ellos lumínicos. 1 2Capítulo 1. Introducción En segundo lugar, en el proyecto [ 5 ] se desarrolla un sistema automático de control de persianas, las cuales se abren y se cierran a conveniencia controladas por un microcontrolador de tipo Arduino. La lógica que sigue el Arduino para abrir o cerrar las persianas es garantizar en la estancia unas condiciones de iluminación predefinidas previamente por el usuario. Por último, en la referencia [ 3 ] se trata el diseño del control automático de un sistema híbrido de iluminación. En este proyecto se usa tanto luz artificial como luz natural para proporcionar en conjunto la iluminación necesaria en el espacio interior. Sin embargo, la forma de usar la luz natural es llamativa, y a priori podría ser un método alternativo al que se propone en este proyecto. En esta referencia, la luz es reconducida desde el exterior hacia el habitáculo interior a través de tubos de luz. Es un dispositivo que tiene una parte que comunica al exterior y otra al interior, y mediante lentes reconduce la luz hacia el interior. El problema que indica durante el proyecto es que este sistema no tiene un buen desempeño en latitudes altas como las de Europa (el proyecto se realiza en Ecuador). El proyecto que se desarrolla en este documento consiste en determinar con la precisión de minutos (o segundos, si fuera necesario) la posición del Sol en coordenadas topocéntricas. Luego, se va a hacer uso de la altura de las paredes del hangar para determinar cómo de útiles son las propias paredes del edificio para generar sombra en la zona útil del interior. Por último, para aquellas zonas del interior que no sea posible cubrir con la sombra de las paredes, se calculará la porción de techo que necesariamente hace falta cubrir para que garantice sombra en la parte restante, dejando el excedente no cubierto del techo libre para que pueda entrar la luz del sol, en algunos casos luz directa y en otro casos luz difusa. Esa porción de techo que necesariamente ha de ser tapada se hará con persianas domóticas gobernadas por algún tipo de controlador, que podría ser perfectamente Arduino. La principal diferencia que radica en este proyecto con los comentados anteriormente es que no necesita de sensores para poder determinar cuando se deben abrir o cerrar las persianas, sino que ya está a priori definido, tomando como única variable de entrada el día y la hora exactos. Otro de los puntos por el que hace diferente este proyecto frente a los otros es por el uso de las paredes como método de sombra. En general, esto no puede usarse a nivel doméstico (que suele ser la aplicación a la que están principalmente enfocados las referencias anteriores) puesto que las dimensiones de las paredes son considerablemente más pequeñas. El uso de sensores que se usa en los proyectos anteriormente comentados no es un aspecto negativo en absoluto. Sin embargo, pueden dar lugar a comportamientos indeseables en aquellos días de tiempo inestable, como puede ser los días parcialmente nublados. Sí vemos de buen interés el uso de sensores de iluminación para controlar la iluminación interior, por lo que el proyecto [3] podría ser un buen documento de referencia. En el presente documento, el Capítulo 2 presenta el lugar donde se ubica el hangar en estudio. Se muestran imágenes vía satélite del lugar y se obtienen a través de Google Earth datos de relevancia para capítulos posteriores. Acto seguido, en el Capítulo 3 se plantea el problema con mayor detalle, y se propone una estrategia de resolución. Se establecen las pautas que se van a seguir para poder resolver el problema, consiguiendo orden y estructura de secciones sucesivas. En el Capítulo 4 se establecen las ecuaciones de la mecánica orbital, se determina la 3 órbita de la Tierra alrededor del Sol, y se particulariza todo para las coordenadas del lugar en estudio. Se obtienen gráficas del comportamiento del Sol en cuanto a elevación y azimuth, datos claves para capítulos posteriores. Una vez definido el movimiento del Sol para cada día del año, se procede en el Capítulo 5 a calcular la sombra producida por las paredes del hangar, y luego cuál es la porción de techo que es necesario tapar para cada instante de cada día del año. Por último, en el Capítulo 6 se comenta a modo de conclusiones los resultados obtenidos y se proponen opciones para la realización del proyecto. 2 Lugar del proyecto El lugar del cual se desarrolla el proyecto que se trata en este documento está emplazado en la factoría de Airbus San Pablo Norte, ubicado en las inmediaciones del Aeropuer- to de Sevilla. Sevilla tiene históricamente una relevancia importante en la fabricación y puesta en servicio de aeronaves. Ya en los años 40 se establecía como un lugar de referencia aeronáutico, principalmente motivado por la base militar de Tablada, junto a la cual se estableció una factoría de la compañía C.A.S.A. (Construcciones Aeronáuticas Sociedad Anónima), dedicada a la fabricación de aviones emblemáticos como el Heikel He-111. Años más tarde, esta misma compañía sitúa junto al actual Aeropuerto de Sevilla factorías para la fabricación y mantenimiento de diversas aeronaves militares, como son el CN-235 y el C-295. Poco antes de entrar el siglo XXI, Airbus absorbe a la compañía C.A.S.A., y crea la nueva factoría de San Pablo Sur, donde actualmente se fabrica el avión de transporte medio y ligero C-295 y el avión táctico y de transporte pesado A400-M. Figura 2.1 FAL del avión C295 en San Pablo Sur [15]. 5 6Capítulo 2. Lugar del proyecto 2.1 Trabajo que se desarrolla en San Pablo Norte En la actualidad, la factoría de San Pablo Norte continua dedicada al mantenimiento en base de los aviones militares CN-235 y C-295, así como el avión de combate Harrier. San Pablo Norte es, por tanto, un centro MRO cuya normativa se rige según la EASA Parte 145 - Organizaciones de Mantenimiento. Este centro consta de 3 hangares, cada uno de ellos destinados a un propósito específico. A la llegada de un avión para su correspondiente labor de mantenimiento, el primer hangar que visita es el de Línea de Vuelo. Allí se realiza la recepción del avión y se hacen sendas pruebas de rodaje de motores para determinar el estado funcional de los elementos de la aeronave. Una vez creada la relación de defectos de entrada, pasa al hangar de mantenimiento en base. Este será el hangar objeto de estudio de este proyecto. Se trata de un hangar en el que entran 6 aviones C295, 3 de ellos alineados en un lado, y otros 3 alineados en el otro. Figura 2.2 Disposición de los aviones en el hangar en estudio. Imagen propia, manteniendo la relación de aspecto en la escala. Dibujo de avión C295 obtenido de [9]. Al llegar el avión a este hangar, se empiezan a realizar todas las tareas de mantenimiento que indica la CAMO del cliente y que previamente ha contratado. La CAMO es la figura encargada de controlar el estado de mantenimiento de la aeronave, en base al manual de mantenimiento que le suministra el fabricante. La CAMO puede ser llevada por el propio cliente, o puede ser subcontratado el servicio a terceros para que le mantengan al día todas 2.2 Obtención de datos: Google Earth Pro 7 las revisiones obligatorias de su aeronave. Es el fabricante del aeronave el encargado de generar el libro de mantenimiento que acompaña al avión cuando se vende, en el cual indica cuándo debe realizarse una comprobación del estado de servicio de según qué elementos, de forma obligatoria. Dependiendo de la criticidad de cada uno de los sistemas del avión y de si puede ser reparable o no, el fabricante determina cada cuantos ciclos/horas de vuelo debe revisarse el elemento, o por el contrario, cada cuantos ciclos/horas de vuelo hay que cambiar el elemento, aún estando serviciable (vida límite de los materiales). Por tanto, los trabajos de mantenimiento que se acometen en este centro son totalmente irrepetibles de un avión a otro. Cada avión tiene su paquete de trabajo personalizado definido por el departamento de Ingeniería de Fabricación en base a lo dictaminado por la CAMO, y es el departamento de Producción el encargado de llevar a cabo todas esas tareas: desmontajes, inspecciones, reparaciones, pruebas y montajes. 2.2 Obtención de datos: Google Earth Pro En esta sección nos centraremos en encontrar la localización exacta del hangar en el globo terráqueo, así como sus dimensiones, puesto que serán necesarias para capítulos posteriores. Para ello, usaremos la herramienta Google Earth Pro [ 1 ]. Esta aplicación de ordenador es gratuita y permite localizar cualquier punto de la Tierra, proporcionando datos de geolocalización muy precisas, y permite medición de distancia y altura. Además, en la versión más reciente Google utiliza la tecnología de superposición de imágenes para poder generar espacios tridimensionales de muy buena calidad, siendo una herramienta perfecta para describir con suficiente detalle el proyecto. Figura 2.3 Captura de inicio del programa [1]. 8Capítulo 2. Lugar del proyecto En la siguiente figura puede verse una imagen general de la factoría, a vista de satélite. Se ha marcado en rojo el hangar propio del estudio. Fíjese como es la disposición de la cubierta del tejado respecto de las naves aledañas, sin dejar de lado la ubicación del Norte, que se puede ver en la esquina superior derecha. Puede verse tras comparar este hangar con el que hay inmediatamente superior a él (marcado en azul) que los dientes de sierra del tejado están girados 90º. Sabiendo que el amanecer se produce en el Este y el atardecer en el Oeste, y que el Sol pasa por el sur en todo su recorrido, puede decirse que la factoría de San Pablo Norte está preparada para generar sombra durante el amanecer (por la orientación de sus dientes de sierra), momento en el cual el Sol tiene su menor intensidad. Por otro lado, la nave marcada en azul tiene los dientes de sierra provistos para generar sombra en el momento en el que el Sol llega a su elevación máxima. Figura 2.4 Captura general de San Pablo Norte y hangar del proyecto. Debido a los trabajos que se desarrollan en el interior del hangar, no puede entrar nada de luz directa, puesto que podría deslumbrar a cualquier persona y eso desencadenar en un accidente en el peor de los casos. Por esta razón, se tapó todo el techo del hangar del estudio con unas planchas opacas de espuma, de unas dimensiones de 1.2m x 1.5m, no permitiendo la entrada absoluta de luz natural. Ver Figura 2.5. 2.2.1 Azimuth del hangar, latitud y longitud La posición geográfica del hangar es un dato relevante para poder saber cuál es el movimiento que describe el Sol en la esfera celeste. La esfera celeste es una esfera imaginaria, donde se proyectan radialmente todos los cuerpos que hay en el espacio [ 10 ]. En nuestro caso, obtener estos datos es algo inmediato usando Google Earth. 2.2 Obtención de datos: Google Earth Pro 9 Figura 2.5 Fotografía del interior del hangar en estudio. Cortesía de Airbus San Pablo Norte. •Latitud y longitud: Situamos el cursor en el centro del hangar, y en la esquina inferior derecha aparece la posición geográfica de ese punto. Por tanto, obtenemos: λ=−5°52059.5400 =−5.883221°φ=37°25023.8900 =37.423292°(2.1) Figura 2.6 Captura en Google Earth de la situación del hangar. Latitud y longitud. 3 Planteamiento del problema y estrategia de resolución En el punto 2.2 ya se ha introducido parcialmente el origen del problema. Como se ha podido ver, la orientación de los dientes de sierra del tejado de la nave no están dispuestos en su forma óptima (ver Figura 2.4) para cancelar los rayos de Sol en las horas centrales del día, que es cuando alcanza su elevación máxima. Es por ello que todo el interior del hangar tuvo que ser tapado con planchas opacas que no permitieran nada de luz (ver Figura 2.5). Sin embargo, esto genera unos problemas secundarios y que son los que se intentan minimizar con este proyecto. •Fatiga visual y cansancio: El hecho de que no pueda entrar nada de luz natural difusa obliga a disponer de potentes focos de iluminación artificial por todo el interior de la nave (ver Figura 2.5), para alcanzar los niveles mínimos de iluminación recomendados por el departamento de Seguridad y Salud Laboral. Aún así, después de varias horas durante la jornada laboral, es notoria la fatiga visual que se acumula, produciendo mayor cansancio. Permitir que entre luz natural, por poca que sea, proporciona un plus de iluminación que actualmente no existe. •Pobre uso de los recursos naturales : Por otro lado, no hay que dejar de lado el enorme consumo que supone mantener la iluminación artificial a pleno rendimiento durante los turnos de trabajo. Esto supone coste elevado de consumo de electricidad y coste en mantenimiento de los focos. Las lámparas de cualquier elemento de iluminación tienen provistas por el fabricante un tiempo de vida límite que está directamente relacionado con las horas de funcionamiento. Por tanto, a mayor uso, mayor será la frecuencia con la que hay que cambiar las lámparas, y con ello, mayor coste. Es cierto que con este proyecto no se van a solucionar los problemas comentados anteriormente, puesto que la única manera de encontrar la solución óptima al aprovechamiento de los recursos naturales sería rehacer el diseño de la cubierta del tejado hacia otro con una geometría óptima. Sin embargo, sí va a minimizar los problemas que existen actualmente. 17 18 Capítulo 3. Planteamiento del problema y estrategia de resolución 3.1 Estrategia de resolución Los pasos que se van a seguir para encontrar solución al problema son los siguientes1: 1. Determinar la posición del Sol en cada instante del año, con la precisión que se estime oportuno. Esto es, encontrar los elementos orbitales de posicionamiento del Sol en el sistema de referencia adecuado, para cada instante del año. De este modo, se consigue saber cuando se produce el amanecer y el atardecer de cada día, así como cual es la trayectoria que describe el Sol, siendo totalmente necesario para conocer cuál es la sombra que proyecta. 2. Una vez determinada la posición del Sol, se va a calcular cuál es la sombra que arrojan las paredes y puertas del hangar sobre el interior de él. La proyección de la sombra de las paredes será función de la elevación y Azimuth del Sol, para cada instante de tiempo. Ver apartado 3.1.1. 3. Sabida cuál es la sombra que las propias paredes de la estructura proyectan sobre el interior del hangar, se calcula cuál es la porción de techo que debe taparse para que arroje sombra sobre el resto de superficie interior por la que entraría luz directa (ver apartado 3.1.2). Esta porción de techo que sería necesario tapar se calcula en planta, a la altura de Hre f (ver 2.6). Se hace así para proporcionar un margen de seguridad a los resultados, y así garantizar que siempre se va a proporcionar sombra en el interior del hangar. La forma de proceder más exacta sería proyectando dicha sombra sobre la geometría inclinada de los dientes de sierra. Sin embargo, el procedimiento de cálculo se complica notoriamente, sin proporcionar un mayor valor añadido. En este caso, la imprecisión que se produce en el cálculo de la porción de techo a tapar según la manera que se sugiere, es un hecho positivo, puesto que nos proporciona el factor de seguridad que necesitamos a priori. 4. Sabiendo ya qué fracción de techo debe ser tapado en cada instante t del año, se indica cuales son las persianas domotizadas que deben abrirse o cerrarse. 3.1.1 Cálculo de la sombra arrojada por las paredes En la Figura 3.1 se ve a la izquierda una representación esquemática del hangar, con sus ejes absolutos y locales. Se ha representado en amarillo el Sol, junto con los rayos que provienen de él. Se toma como Az el ángulo formado entre el Norte y el rayo de Sol que pasa por el centro del origen de coordenadas absoluto. Se ha hecho la hipótesis de que todos los rayos de Sol son paralelos a éste, puesto que: • La escala sobre la que estamos trabajando es mucho menor que las dimensiones del radio de la Tierra. •El Sol se considera una estrella "infinitamente lejos" (Ver [12]). En dicha representación se ha realizado el corte A-A a modo de ejemplo, el cual muestra cuál es la sombra que proyecta la pared de la nave en el punto (x0,y0) = (x0 0,0) , para los valores Az(t)yh(t). 1 Se hace la hipótesis de que se quitan todas las planchas opacas comentadas anteriormente, permitiendo la entrada de luz directa por toda la cubierta del tejado. 3.1 Estrategia de resolución 19 N E X Y X’ Y’ AzN A A Az(t)CORTE A-A h(t) x0 0 Csombra (h(t)) Cmax (x0 0,y0 0),Az(t) Hre f Figura 3.1 Sombra proyectada por las paredes. El procedimiento para calcular la sombra de la nave es como sigue: para cada uno de los puntos (x0 0,y0 0) que conforman el contorno del hangar (en ejes locales), se trazan lineas paralelas según la dirección del azimuth del Sol Az , de modo que realicen cortes de las paredes del hangar como el A-A. Luego, se calculará la sombra usando la elevación del Sol h. Por tanto, Csombra (h(t)) = Hre f tan (h(t)) (3.1) Como puede verse en la ecuación 3.1, la longitud de la sombra solo depende de la elevación del Sol, y no del punto del contorno de la nave. Lo único que depende del punto de la nave es la proyección de la sombra sobre el suelo. 3.1.2 Cálculo del techo necesario Como ya se ha comentado, el objeto del proyecto es determinar cuál debe ser la porción de techo óptima que debe taparse para que genere sombra en el interior del hangar. En la Figura 3.2 se ha representado en verde la porción de techo que debe estar necesariamente tapada para generar sombra en el interior a cierta altura Hs . Es importante darse cuenta que la sombra debe estar a partir de cierta altura Hs respecto del suelo de la nave. Esto es así para salvar la altura típica de las personas, con cierto margen de seguridad. Por eso, se determina como valor inicial el que se muestra a continuación, aunque es un parámetro que puede variar según las necesidades. Hs=3m (3.2) 20 Capítulo 3. Planteamiento del problema y estrategia de resolución Por tanto, sabiendo la elevación del Sol y la altura mínima Hs a la cual debe existir sombra, el cálculo de la longitud de techo es inmediata, según la siguiente expresión: Ltecho (x0 0,y0 0),Az(t)=Cmax (x0 0,y0 0),Az(t)−Llibre (h(t)) (3.3) Primero, es necesario saber la longitud de techo que quedaría libre de ser tapada. De la figura, se deduce que: Llibre (h(t)) = Hre f −Hs tan (h(t)) (3.4) CORTE A-A h(t) Cmax (x0 0,y0 0),Az(t) Hre f Hs Ltecho (x0 0,y0 0),Az(t) Llibre Figura 3.2 Techo necesario para generar sombra en el interior del hangar a cierta altura Hs. Una vez sabido esto, es necesario conocer la longitud máxima del corte producido en el hangar según la dirección de Az. Esto es, Cmax (x0 0,y0 0),Az(t). Cálculo de Cmax (x0 0,y0 0),Az(t) Así como Llibre se calcula de forma inmediata, con Cmax no ocurre lo mismo. Al depender de la posición del contorno del hangar y del azimuth del Sol, es necesario tratarlo de forma particular. En primer lugar, trataremos la dependencia que tiene con Az(t) . Si nos fijamos en la Figura 3.1 y pensamos el movimiento que describiría el Sol, de Este a Oeste y pasando por el Sur, la línea roja que corta cada una de las paredes de la nave estaría continuamente girando a medida que aumenta Az(t) . Durante ese recorrido, existen ángulos para los que el corte A-A es máximo, y otros ángulos para los que el corte es mínimo o incluso inexistente. Por ello, es necesario definir las ecuaciones que caracterizan a Cmax (x0 0,y0 0),Az(t) para cuando Az(t) se encuentra en cada uno de los intervalos acotados. Existen 4 intervalos (o cuadrantes), en los cuales el Sol iluminará siempre 2 caras del hangar, hasta llegar al punto de cambio. El punto de cambio es aquel instante en el que Az(t) pasa de un cuadrante al 3.1 Estrategia de resolución 21 siguiente. Cuando Az(t) es exactamente igual al punto de cambio, sólo en ese momento el Sol iluminará una de las caras, pasando los rayos de sol paralelos a las colindantes. •Cuadrante 1: Az(t)∈[AzN−90°,AzN) –Az(t)<α1 En la Figura 3.3 se ha representado el caso para el cual el Sol se encuentra con Az(t)∈[Az−90°,Az) . Se ha denotado cada lado del hangar con los nombres Lado a (Este), Lado b (Sur), Lado c (Oeste) y Lado d (Norte). A su vez, cada esquina del hangar se ha llamado Vértice y un número que aumenta en el sentido de las agujas del reloj, siendo el Vértice 1 el que se encuentra en el Noreste, y el Vértice 4 en el Noroeste. Como ya se ha comentado, el Sol ilumina siempre 2 caras, por lo que en el caso que se representa son los Lados a yLado d. Se han pintado estos lados de amarillo. Por otro lado, puede verse como estas caras tienen ciertos puntos equidistribuidos desde los que parten ciertas líneas rojas, azules y verdes. El hecho de que existan esos puntos están relacionados con la programación en Matlab (ver sección 5.3). Se han discretizado los lados del hangar en un número finito y manejable de divisiones, de modo que las representaciones en Matlab puedan ser entendibles. En general, los lados paralelos al eje local x0 (Lados b yd) se han dividido un número total de nx divisiones, mientras que los lados paralelos al eje local y0 (Lados a yc) se han dividido un número total de ny divisiones. El sentido positivo de avance en el recorrido de estos puntos es siempre a favor de las agujas del reloj, puesto que también así es el movimiento del Sol. De este modo, la división d1 que se ve en la figura tendrá una longitud Lx/nx , mientras que la división que corresponde al penúltimo punto (nx−1) tendrá una longitud desde el Vértice 4 de (Lx/nx)(nx−1) . Dicho esto, la longitud de cualquier punto intermedio ktendrá una longitud de: Ld,k=Lx nx k k =1,2,...,nx(3.5) De forma análoga ocurre con el Lado a. Existen sendos puntos equidistribuidos que dividen el lado en ny subdivisiones, siendo q un punto cualquiera de las subdivisiones. En este caso, la longitud de los lados se mide respecto del Vértice 2, aunque el contador qavance en sentido positivo desde el Vértice 1 al Vértice 2, según el movimiento de las agujas del reloj. Esto es así porque interesa para poder usar el triángulo derecho que se ve en la parte inferior de la figura, y que contiene el Vértice 2. Por tanto, en este caso, la definición de longitud de lado para la subdivisión qserá: La,q=Ly−Ly ny qq=1,2,...,ny(3.6) Las líneas azules y rojas del dibujo son así para delimitar dos zonas bien diferenciadas en las que las longitudes son diferentes. Las líneas verdes por otro lado indican uno de los puntos genéricos de las subdivisiones. Aún así, todas ellas son paralelas a la dirección marcada por el azimuth del Sol. 22 Capítulo 3. Planteamiento del problema y estrategia de resolución N E X Y x’ y’ AzN Az(t) (AzN−Az(t)) Lado d Lado a Vértice 1 Vértice 2 Vértice 3 Vértice 4 Lado b Lado c d1 d2 dk d4 d5 d... d,(nx−1) a1 a2 a3 a4 aq a,... a,(ny−1) Ly Lx Cmax,dk Vértice 4 Cmax,dk (AzN−Az(t)) C1max dk Lx nx k Cmax,ak Vértice 2 aq (AzN−Az(t)) Ly−Ly ny q Cmax,ak Xcambio α1 Figura 3.3 Cuadrante 1. Az<α1. 3.1 Estrategia de resolución 23 Si analizamos las partes del dibujo que engloban líneas azules (del Lado d, por ejemplo), éstas van aumentando su longitud desde el Vértice 4 en su recorrido hasta el Vértice 2, hasta llegar a un punto en el que la programación en Matlab se ha llamado Xcambio2 . Si nos damos cuenta, esta línea roja cae directamente en el Vértice 3, y a partir de ahí todas las subdivisiones k (según la expresión 3.5) mayores que Xcambio tendrán una longitud fija e igual a esta distancia, hasta llegar al Vértice 1, que delimita la segunda línea roja a partir de la cual cambia de nuevo el comportamiento. En la figura, esta distancia es C1max . Por tanto, el Lado d tendrá dos zonas bien diferenciadas para las longitudes: una variable y otra fija y máxima. En el Lado a, sin embargo, no existe punto de cambio a medida que se avanza con el contador qdesde el Vértice 1 hacia el 2. Por tanto, todas las líneas serán azules e irán decreciendo en longitud hasta llegar al Vértice 2. Este razonamiento es común para todos los lados y para todos los ángulos de azimuth del Sol. Lo único que cambiará serán las definiciones de la longitud de esas líneas azules y rojas dependiendo de en qué cuadrante se encuentre el Sol. El dibujo que se ha representado en la Figura 3.3 será válido para cuando Az(t)∈[Az−90°,Az) y además Az<α1 . Se define el ángulo α1 como el ángulo formado entre el Norte y el Vértice 1. Así, cada vértice tendrá un αi , definido de la misma manera. Ver sección 5.1. Una vez definidos todos los términos de relevancia, pasamos a definir las expresiones matemáticas de ellos. C1max =Ly sin (AzN−Az(t)),Lx≥Ld,k≥Xcambio (3.7) Xcambio =C1max cos (AzN−Az(t)) (3.8) Cmax,d,k= Lx nx k cos (AzN−Az(t)) Ld,k<Xcambio (3.9) Cmax,a,q= Ly−Ly ny q sin (AzN−Az(t)) (3.10) –Az(t)≥α1 Como ya se ha comentado con suficiente detalle en el anterior punto cuál es la forma de proceder, se va a limitar únicamente a actualizar el 2 En este caso se ha llamado Xcambio porque el Lado d es paralelo al eje x0 . Si el punto de cambio hubiera caído en un lado paralelo al eje y0, se llamaría Ycambio. 24 Capítulo 3. Planteamiento del problema y estrategia de resolución valor que deben tomar las variables. C1max =Lx cos (AzN−Az(t)),Ly≥La,q≥Ycambio (3.11) Ycambio =C1max sin (AzN−Az(t)) (3.12) Cmax,d,k= Lx nx k cos (AzN−Az(t)) (3.13) Cmax,a,q= Ly−Ly ny q sin (AzN−Az(t)) La,q<Ycambio (3.14) d1 d2 dk d4 d5 d... d,(nx−1) a1 a2 a3 a4 aq a,... a,(ny−1) Lado d Lado a Vértice 1 Vértice 2 Vértice 3 Vértice 4 Lado b Lado c Ly Lx Ycambio (AzN−Az(t)) N E x’ y’ X Y AzN Az(t) Cmax,dk C1max Cmax,ak α1 Figura 3.4 Cuadrante 1. AzN>Az≥α1. 3.1 Estrategia de resolución 25 •Cuadrante 2: Az(t)∈[AzN,AzN+90°) –Az(t)<α2 C2max =Lx cos (Az(t)−AzN),La,q≤Ycambio (3.15) Ycambio =Ly−C2max sin (Az(t)−AzN)(3.16) Cmax,a,q= Ly ny q sin (Az(t)−AzN)La,q>Ycambio (3.17) Cmax,b,k= Lx−Lx nx k cos (Az(t)−AzN)(3.18) –Az(t)≥α2 Lb,k=Lx nx k k =1,2,...,nx(3.19) C2max =Ly sin (Az(t)−AzN)(3.20) Xcambio =Lx−C2max cos (Az(t)−AzN)(3.21) Cmax,a,q= Ly ny q sin (Az(t)−AzN)(3.22) Cmax,b,k= Lx−Lx nx k cos (Az(t)−AzN),Lb,k<Xcambio (3.23) •Cuadrante 3: Az(t)∈[AzN+90°,AzN+180°) –Az(t)<α3 32 Capítulo 4. Definición de la posición del Sol observador (en nuestro caso, Sevilla). Este sistema de referencia se conoce como Sistema Topocéntrico. Figura 4.4 Sistema de referencia Topocénctrico. Dibujo sacado de [10]. Tiene su origen de coordenadas en el lugar del observador, por lo que también está ligado al movimiento de la Tierra, y es útil para mediciones astronómicas. El eje X coincide con el Este y el eje Y con el Norte. El eje Z está ubicado en la vertical del lugar, también conocido como cénit. Sus variables de posicionamiento son ρ (distancia al cuerpo), la cual no se va a utilizar; Azimuth Az , que es el ángulo que hay entre el Norte y la proyección del cuerpo sobre el plano terrestre; y por último la elevación h , que es el ángulo que forma el objeto con el plano horizontal. El dominio de estos ángulos son: Az ∈[0,360°] y h∈[−90°,90°] . 4.1.2 Tiempos de uso en astronáutica Tan importante es utilizar el sistema de referencia apropiado según el problema en el que nos encontremos, como ubicarlo bien en el tiempo. Al estar tratando con cuerpos móviles y masivos como son los cuerpos celestes, el movimiento de ellos viene impuesto, sin posibilidad alguna de poder alterarlos. Por ello, es necesario tomar tiempos de referencia para ubicarnos correctamente en el punto adecuado de la cinemática de los cuerpos. Existen muchas definiciones de tiempos usados en este campo de la ciencia. Unos se basan en la inmovilidad de las estrellas lejanas. Estos se conocen como sidéreos. Otros sin embargo, son artificios matemáticos creados para hacer el uso de los tiempos en la vida cotidiana un poco más fácil. Estos serán los que usaremos en el proyecto. Todos los tiempos que aparecen a continuación han sido obtenidos de [10]. •Sol Medio: es un "sol ficticio", que se mueve a velocidad constante y en el plano del ecuador. Está definido de forma que coincide con el Sol real (si el Ecuador y la eclíptica coincidieran) cuando la Tierra se encuentra en su perihelio y afelio. El perihelio y el afelio son aquellos puntos de la órbita elíptica de la Tierra alrededor del Sol en el que la Tierra está más próxima y más alejada, respectivamente, del Sol. •HSM (Hora Solar Media): Esta es la hora que marcaría un reloj de Sol, si el Sol fuera el Sol Medio. La hora solar media en Greenwich es la base de UT. •UT (Universal Time): Previamente llamado GMT, se basa en el día solar medio. 4.1 Introducción general a la Mecánica Orbital 33 •Día Solar Medio: En base al Sol Medio, se considera que el día tiene 24h. Se define el día Solar Medio como 1/365.2421988 de un año tropical. De esta manera, se consigue que el día conste de 24h exactas. •Año Tropical: Es el tiempo transcurrido entre 2 Equinoccios sucesivos. Este es el período del Sol Medio. Por tanto, a modo de corolario, los tiempos que se van a usar en el proyecto son: •Duración de 1 día: 24h •Duración de 1 año:365.25días •Hora de referencia: La que tiene el meridiano de Greenwich en UT. 4.1.3 Elementos orbitales Si usamos el Sistema de Referencia Geocéntrico Ecuatorial (ver Figura 4.2), la órbita del Sol se comportaría a ojos de la Tierra como "una órbita más", como si fuera la de cualquier cuerpo orbitando alrededor de la Tierra. Es por ello necesario definirla. Para definir cualquier órbita en el espacio, es necesaria la definición de ciertos conceptos previos. Todos los términos que se definen a continuación han sido obtenidos de [11]. Figura 4.5 Definición de elementos orbitales. Dibujo sacado de [11]. •Plano de referencia: En nuestro caso, el plano del Ecuador. • Plano orbital: Plano en el cual se mueve el cuerpo del que se quiere determinar la posición, en nuestro caso, el Sol. Este plano sería el plano de la eclíptica. 34 Capítulo 4. Definición de la posición del Sol • Línea de nodos: Es la línea resultante de la intersección entre el plano orbital (eclíptica) y el de referencia (plano ecuatorial). En nuestro caso, la línea de nodos coincide con  , tal y como se definió al comienzo del apartado 4.1.1. • Nodo ascendente: La órbita corta la línea de nodos en dos puntos. Por uno de ellos asciende, y por el otro desciende. Pues bien, el nodo por el cual el cuerpo asciende se conoce como nodo ascendente, y se denota por el símbolo k. • Ascensión Recta del Nodo Ascendente: Es el ángulo que existe entre  y k . Se denota por Ω . En nuestro caso, Ωk=0 , puesto que el nodo ascendente se sitúa en la línea de nodos, y esta a su vez coincide con la dirección de  , por definición. Además, el Sol empieza a ascender en la órbita a partir del Equinoccio de Primavera. • Argumento de perihelio ω y Anomalía verdadera θ :El argumento de perihelio ω es el ángulo formado por  y la línea de ápsides (por donde pasa el perihelio y el afelio). La anomalía verdadera θ(t) es el ángulo entre la posición del cuerpo a posicionar y la línea de ápsides. En nuestro caso, al hacer la hipótesis de que la órbita de la Tierra alrededor del Sol es circular, ninguno de los dos elementos orbitales está definido. En ese caso, se define un término adicional tal y como sigue: u(t) = ω+θ(t) . Éste término tiene su u=0 en k , es decir, en el Equinoccio de Primavera. Luego, recorre toda la órbita a medida que va avanzando el cuerpo, en este caso el Sol. • Inclinación de la órbita: Es el ángulo que forma el plano orbital con el plano de referencia. Se denota por i. En nuestro caso, i=ε=23.5°. 4.1.4 Trigonometría esférica Una vez realizada una breve introducción a los conceptos más relevantes para nuestro proyecto, podemos introducir los procedimientos y conceptos para el cálculo de las variables buscadas. En Mecánica Orbital es habitual trabajar con distancias muy grandes y en cuerpos esféricos. En ese caso, la trigonometría plana deja de tener aplicación, puesto que el error que se comete es inadmisible. Por eso, para solventar el problema se recurre al uso de Trigonometría esférica. Este tipo de trigonometría hace lo que la plana, pero en entornos esféricos.Véase [8] en el que se realiza un estudio para seguimiento del Sol. La condición que se deben cumplir para el uso de la trigonometría esférica es que todos los lados (rectas esféricas o geodésicas) del triángulo esférico deben pertenecer a un círculo máximo. Como se puede ver en la Figura 4.6, todos los círculos que cortan la esfera y que generan el triángulo contienen el centro de la esfera. Un círculo máximo es aquel que contiene el centro de la esfera. Para este tipo de trigonometría, la ley de los cosenos queda 4.2 Cálculo de hyAz35 Figura 4.6 Triángulo esférico. Dibujo sacado de [12]. como sigue, apoyándose en el triángulo esférico de la Figura 4.6: cos a=cos bcos c+sin bsin ccos α(4.1) cos b=cos acos c+sin asin ccos β(4.2) cos c=cos bcos a+sin bsin acos γ(4.3) cos α=−cos βcos γ+sin βsin γcos a(4.4) cos β=−cos αcos γ+sin αsin γcos b(4.5) cos γ=−cos βcos α+sin βsin αcos c(4.6) 4.2 Cálculo de hyAz Vista la manera de calcular las rectas esféricas y los ángulos de un triángulo esférico, podemos finalmente plantear cómo se calcula la elevación y el azimuth del Sol, desde el punto en estudio. Para ello, se recurre al llamado triángulo astronómico. Este triángulo es de aplicación sobre cuerpos infinitamente distantes de la Tierra, con declinación δ y cuyas coordenadas topocéntricas respecto a un observador son su elevación h y su azimuth Az . Por otro lado, el observador debe encontrarse en un punto de la Tierra de latitud φ , y con ángulo horario del cuerpo respecto al observador llamado HS . El ángulo horario es por definición la diferencia entre el meridiano del observador y el meridiano en el que se encuentra el cuerpo al que se 36 Capítulo 4. Definición de la posición del Sol observa. Este ángulo se define matemáticamente como sigue: HS=LST −ARS=GST0+λ+ω⊕t−ARS(4.7) donde •LST :Local Sidereal Time. Es el tiempo sidéreo local, o dicho de otro modo, la ascensión recta del meridiano del observador. •ARS: Esta es la ascensión recta del Sol. •GST0 : Esta es la ascensión recta del meridiano de Greenwich a las 0:00 UT del día de referencia. Este dato se puede calcular o mirar en tablas. En nuestro caso, lo calcularemos como se indica más adelante. Como puede verse, para poder saber LST del observador, es necesario conocer GST0 en el instante deseado. Nosotros lo haremos usando Días Julianos. Este concepto consigue ubicar cualquier día de cualquier año, sea bisiesto o no, a partir únicamente del día, el mes y el año en cuestión. El resultado es un número bastante largo (puesto que la fórmula es un conteo de días desde el 1 de Enero del año 4713 AC a las 12:00 UT). Cuando el número resultante a través de la fórmula es exacto, se ubica el día a las 12:00 UT. Cuando el número resultante tiene 0.5 decimales, entonces nos encontramos a las 0:00 UT. La fórmula que permite calcular el día Juliano a partir del año A , el mes M y el día D es así (datos sacados de [10]): JD =367A−$7A+7$M+9 12 % 4%+$275M 9%+D+1721013.5(4.8) Dicho esto, la forma de calcular GST0es como sigue (datos sacados de [10]): 1. Calcular la "Centuria Juliana" T0 : Aquí J0 es el día Juliano calculado con la ecuación anterior para el día de referencia. T0=J0−2451545 36525 (4.9) 2. Calcular GST0: GST0=100.4606184 +36000.77004T0+0.000387933T2 0−2.583 ×10−8T3 0 (4.10) Aplicando este procedimiento al 1 de Enero de 2021 a las 0:00 UT, se llega a: GST0=100.8688°(4.11) Para el cálculo de la velocidad angular de la Tierra partimos de la duración del día sidéreo (T⊕=23 h56 m4s=86164 s). Por tanto, ω⊕=2π T⊕ =0.000072921158rad s(4.12) 4.2 Cálculo de hyAz37 Lo único que nos falta por calcular de la ecuación 4.7 es ARS . Para ello, se van a usar las leyes de senos, cosenos, la identidad trigonométrica habitual y la Figura 4.7. Hay que recordar que en nuestro caso se está considerando que la órbita que describe la Tierra alrededor del Sol es circular, por lo que u=ω+θen la Figura 4.7. Por un lado, de la ley de senos habitual (que también es de aplicación en trigonometría esférica) se deduce que: sin δ=sin usin i(4.13) Por otro lado, de la expresión 4.3 particularizada para γ=90 °, a=δ , c=u , b=λu y α=i, se llega a: cos u=cos λucos δ⇒cos δ=cos u cos λu (4.14) Sustituimos en la identidad trigonométrica este último resultado y 4.13, sin2δ+cos2δ=1⇒sin2usin2i+cos2u cos2λu =1(4.15) También de la identidad trigonométrica se deduce que, sin2λu+cos2λu=1⇒tan2λu+1=1 cos2λu (4.16) Sustituyendo este resultado en 4.15, se llega a, sin2usin2i+cos2utan2λu+1=1(4.17) de donde, si despejamos tan2λunos queda, tan2λu=1−cos2u−sin2usin2i cos2u=1−sin2isin2u cos2u=cos2itan2u(4.18) llegando finalmente al resultado tan λu=cos itan u(4.19) Figura 4.7 Triángulo esférico. Dibujo modificado de [11]. 38 Capítulo 4. Definición de la posición del Sol El término λu es el ángulo que va desde k hasta el meridiano en el que se encuentra el cuerpo, en este caso el Sol. Como ya se ha comentado varias veces en este documento, el nodo ascendente k coincide con  por definición para el Sol, por lo que la ascensión recta del nodo ascendente es cero, dando lugar a la identidad AR=λu . Por tanto, conocido λu conocemos ARS. Sin embargo, para poder aplicar la expresión 4.19 necesitamos primero conocer u . El ángulo u es el que recorre la órbita junto con el objeto en estudio, desde k . Como se considera que la órbita es circular, la velocidad de movimiento del Sol en la eclíptica será constante. Es por ello que el movimiento del Sol será descrito por la siguiente expresión: u(t) = u0+˙u(t−t0)(4.20) donde u0 será el valor que tome u en el instante t0 . Para comodida, se exige que t0 coincida con el momento exacto del Equinoccio de Primavera, puesto que es en este momento cuando u=0 °. Para saber cuando ocurre el Equinoccio de Primavera en 2021 se consulta la referencia [ 6 ]. Se ve por tanto que t0 ocurre el 20 de Marzo, a las 09:37 UT. Sin embargo, Figura 4.8 Equinoccios y Solsticios tabulados en la serie 2020-2025 [6]. en nuestro proyecto nos interesa que el tiempo empiece el 1 de Enero a las 0:00 UT (en lo que sigue, esta fecha será t1 ), para entender mejor el movimiento del Sol en un año natural. En ese caso, hay que propagar u0 hasta el instante t1 . Hay que tener en cuenta que la variable t en todas las ecuaciones de la mecánica orbital están en la escala de segundos. Es por ello que hay que saber el número de días que hay desde el 1 de Enero a las 0:00 UT hasta el 20 de Marzo a las 0:00 UT. Por simple conteo, se infiere que hay 78 días. Entonces, 4.2 Cálculo de hyAz39 el tiempo que hay desde el Equinoccio de Primavera hasta el 1 de Enero a las 0:00 UT es: t1=78dias·24 h 1dia ·60min 1h·60s 1min +9·60min 1h·60s 1min +37 ·60s 1min =6739200 s+34620 s= =6773820 s (4.21) Antes de poder propagar u0, necesitamos saber ˙u. ˙u=360° 365.25dias ·1dia 24h·1h 60min ·1min 60s=1 87660 ° s(4.22) Por tanto, usando la expresión 4.20, sabiendo que ∆t=−t1, obtenemos u0=−t1˙u=−77.2738°(4.23) Así, la expresión 4.20 queda perfectamente determinada, siendo t0 el 1 de Enero de 2021 a las 0:00 UT. Habiendo definido u(t) , la obtención de λu es inmediata a partir de la expresión 4.19, sabiendo que i=ε=23.5°. En efecto, ARS=λu=tan−1(cos εtan u)(4.24) Por tanto, sustituyendo λ de 2.1, GST0 de 4.11, ω⊕ de 4.12 y ARS de 4.24 en 4.7, queda perfectamente determinado el ángulo horario entre ambos meridianos. Conocido entonces el ángulo horario, conocer la elevación y el azimuth es inmediato. En la Figura 4.9 se ha representado el triángulo astronómico para un observador con latitud φ y longitud λ . Notar que tiene representado los ejes cardinales en el plano que está inclinado φ respecto de la estrella polar (punto P). El punto S representa el Sol, el cual tiene declinación δS respecto del ecuador celeste. El mismo punto S tiene elevación h , esta vez respecto del sistema de referencia topocéntrico. Por tanto, para generar el triángulo astronómico se hace pasar un círculo máximo por el Sol y por P. Dicho círculo máximo coincide con el meridiano que pasa por el Sol. Por otro lado, se hace pasar un segundo círculo máximo que contenga al Sol y el cénit del observador (punto Z). Por último, se traza un último círculo máximo que contenga el cénit y el punto P. De este modo, se obtiene el triángulo de la figura. Usando la ecuación 4.3 para obtener δS y la ecuación 4.1 para obtener h , donde γ=Az , c=90°−δS,b=90°−φ,α=−HSya=90°−h, llegamos a las expresiones: sin δS=sin φsin h+cos φcos hcos Az (4.25) sin h=sin φsin δS+cos φcos δScos Hs(4.26) 40 Capítulo 4. Definición de la posición del Sol Aplicando la ley de senos en el triángulo de la Figura 4.7, obtenemos sin δ sin i=sin u sin 90°⇒sin δ=sin usin i⇒δ=sin−1(sin usin i) (4.27) Figura 4.9 Triángulo Astronómico. Dibujo obtenido de [14]. Sustituyendo 4.27 en 4.25, obtenemos: sin usin i=sin φsin h+cos φcos hcos Az⇒ ⇒cos Az=sin usin i−sin φsin h cos φcos h (4.28) Del mismo modo, sustituyendo 4.27 en 4.26, obtenemos: h=sin−1(sin φ(sin usin i))+cos φcos sin−1(sin usin i)cos HS(4.29) 4.3 Representación gráfica del movimiento del Sol A continuación se representan en diferentes gráficas las curvas a las que dan lugar las ecuaciones encontradas en el apartado anterior. 4.3 Representación gráfica del movimiento del Sol 41 Figura 4.10 Elevación del Sol para todos los días del año. Figura 4.11 Azimuth del Sol para todos los días del año. 48 Capítulo 5. Obtención computacional de sombra usando Matlab Vértice 1 Vértice 2 Vértice 3 Vértice 4 Ly Lx N E X Y η x’ y’ O O’ EJES LOCALES EJES GLOBALES o ABSOLUTOS AzN α1 Figura 5.1 Ángulos de los vértices del hangar. donde (x0,y0) son las coordenadas del punto expresado en ejes locales, y (x,y) las coordenadas del mismo punto en ejes globales. Como suele ser habitual con las matrices de giro, la matriz G es ortogonal, y por tanto su inversa es igual que su traspuesta. Nótese en este caso que se cumple que G=GT Sin embargo, el cambio de base no está completamente resuelto, puesto que el sistema de coordenadas local no solo está rotado sino que además está desplazado respecto del origen de coordenadas global (ver Figura 5.2). Por tanto, después de efectuar el giro, es necesario realizar una traslación según el vector ~rO0|global, definido como: ~rO0|global =1 2qL2 x+L2 ysin (180°−α2) −cos (180°−α2)(5.7) Por tanto, el cambio de base quedaría como sigue:   x0 y0 =  −cos (αN−90°)sin (αN−90°) sin (αN−90°)cos (αN−90°)  | {z } G   x y +1 2qL2 x+L2 ysin (180°−α2) −cos (180°−α2) | {z } ~rO0|global (5.8) 5.2 Sombra arrojada por las paredes del hangar 49 De este modo, cualquier punto k podrá expresarse en un sistema de coordenadas o en otro, según convenga, a partir de la siguiente relación (ver Figura 5.3): ~rkO|global =~rkO0|global +~rO|global (5.9) donde ~rO|global =−~rO0|global Vértice 1 Vértice 2 Vértice 3 Vértice 4 Ly Lx N E X Y x’ y’ O O’ EJES LOCALES EJES GLOBALES o ABSOLUTOS AzN−90° AzN−90° ~rO0|global α2 AzN Figura 5.2 Desplazamiento y giro del eje de coordenadas local. Sabiendo esto, estamos en disposición de poder calcular la posición de los vértices del hangar. Para ello, los definiremos en ejes locales. Luego, sería necesario aplicar 5.6 para convertirlos en ejes globales y por último aplicar 5.9. ~r1|local =0 Ly~r2|local =0 0~r3|local =Lx 0~r4|local =Lx Ly(5.10) 5.2 Sombra arrojada por las paredes del hangar El cálculo de la sombra que arrojan las paredes se hará según la expresión 3.1. Si nos fijamos en esta expresión, la elevación se encuentra en el denominador, por lo que cuando el Sol se encuentra cercano al amanecer o al atardecer (elevación con valores cercanos a cero o incluso cero), el cociente tiende a infinito. En ese caso, la sombra proyectada también tendería a infinito, lo que es sinónimo de que toda la superficie del hangar estaría cubierta por sombra. Por tanto, vamos a definir un nuevo término llamado hmin que será la elevación mínima que debe tener el Sol para que deje de arrojar sombra en toda la superficie 50 Capítulo 5. Obtención computacional de sombra usando Matlab Vértice 1 Vértice 2 Vértice 3 Vértice 4 Ly Lx N E X Y x’ y’ O O’ EJES LOCALES EJES GLOBALES o ABSOLUTOS ~rO0|global ~rkO|global ~rkO0|global k Figura 5.3 Punto kexpresado en ambos ejes de coordenadas. -100 -50 0 50 100 -100 -80 -60 -40 -20 0 20 40 60 80 100 X 32.8442 Y -42.6293 X 45.8415 Y 28.1878 X -45.8415 Y -28.1878 X -32.8442 Y 42.6293 Vértice 1 Vértice 2 Vértice 3 Vértice 4 Figura 5.4 Representación gráfica en Matlab de los lados y los vértices del hangar. 5.2 Sombra arrojada por las paredes del hangar 51 interior del hangar. La distancia más restrictiva de sombra proyectada por las paredes sería la diagonal de la nave. Por tanto, hmin =tan−1  Hre f −Hs qL2 x+L2 y =6.887°(5.11) CORTE A-A qL2 x+L2 y Hre f Hs hmin Figura 5.5 Cálculo de hmin. Habiendo definido la condición h≥hmin , hacemos un bucle que avance en el tiempo, lo cual dará de forma inequívoca pares de Az y h para cada instante de tiempo. El valor de h nos permite resolver la expresión 3.1, mientras que el valor de Az nos permite conocer la dirección hacia donde se proyecta la sombra. Esta dirección vendrá dada por el vector unitario en ejes globales~esombra: ~esombra =~e=−sin Az cos Az(5.12) En la Figura 5.6 se ha representado a modo de ejemplo el vector sombra del vértice 4, en ejes globales, aunque se define de la misma manera para el resto de vértices. Sería como sigue: ~rsombra,i=~ri|global +Csombra ~ei=1,2,3,4(5.13) Sabiendo cómo resulta el procedimiento del cálculo de la sombra arrojada por las paredes, se representa a continuación los resultados obtenidos en dos días del año: el 2 de Enero y el 1 de Julio, ambos de 2021. Se han elegido estos días por estar muy distantes entre sí. De hecho, el primero está muy cercano al Solsticio de Invierno, mientras que el segundo está muy cercano al Solsticio de Verano. Se han hecho las simulaciones a la misma hora del día (hora local del lugar, sin tener en cuenta el horario de verano) para poder comparar entre ambas fechas. Las Figuras que van desde la5.7 hasta la 5.10 son las dedicadas al 2 de Enero, mientras que las que van desde la Figura 5.11 hasta la 5.14 son las dedicadas al 1 de Julio. Puede verse la gran diferencia que existe en la sombra arrojada en la comparativa de ambas fechas a las mismas horas. Mientras que el 2 de Enero a las 10h de la mañana la sombra arrojada por las paredes cubre prácticamente toda la superficie del hangar, el 1 de Julio sin embargo sólo cubre una mínima parte. Téngase en cuenta en ambos casos la elevación del Sol. En 52 Capítulo 5. Obtención computacional de sombra usando Matlab Vértice 1 Vértice 2 Vértice 3 Vértice 4 Ly Lx N E X Y x’ y’ O O’ EJES LOCALES EJES GLOBALES o ABSOLUTOS Csombra ~e ~r4|global ~rsombra,4 AzN ~e ~esombra Figura 5.6 Sombra arrojada por las paredes. el primero, la elevación está en torno a los 11.26 °, mientras que en el segundo alcanza los 41.23 °. Véase también en ambas figuras donde se encuentra comparativamente el Sol, representado con un punto amarillo. En el primer caso, está en el Sureste, mientras que en el segundo caso está coincidente con el Este. A modo general comparando todas las gráficas de un mismo día con las del otro, se puede ver cómo el 2 de Enero está gran parte del día cubierto por la sombra que arrojan las paredes, mientras que en el caso del 1 de Julio la porción de superficie cubierta por la sombra de las paredes es mínima. 5.2 Sombra arrojada por las paredes del hangar 53 Figura 5.7 Sombra proyectada por las paredes en el 2 de Enero de 2021, a las 10h, hora local. Figura 5.8 Sombra proyectada por las paredes en el 2 de Enero de 2021, a las 14h, hora local. 54 Capítulo 5. Obtención computacional de sombra usando Matlab Figura 5.9 Sombra proyectada por las paredes en el 2 de Enero de 2021, a las 16h, hora local. Figura 5.10 Sombra proyectada por las paredes en el 2 de Enero de 2021, a las 17.5h, hora local. 5.2 Sombra arrojada por las paredes del hangar 55 Figura 5.11 Sombra proyectada por las paredes en el 1 de Julio de 2021, a las 10h, hora local. Figura 5.12 Sombra proyectada por las paredes en el 1 de Julio de 2021, a las 14h, hora local. 56 Capítulo 5. Obtención computacional de sombra usando Matlab Figura 5.13 Sombra proyectada por las paredes en el 1 de Julio de 2021, a las 16h, hora local. Figura 5.14 Sombra proyectada por las paredes en el 1 de Julio de 2021, a las 17.5h, hora local. 5.3 Cálculo de la porción de techo tapado 57 5.3 Cálculo de la porción de techo tapado Llegados a este punto, sólo falta calcular cuál es la porción de techo que necesariamente debe ser tapada en cada punto y para cada instante del año, para que cumpla con la condición que se mostraba en la Figura 3.2. Para este cálculo, el procedimiento que se ha seguido ha sido el siguiente, apoyándonos en la Figura 5.15. Hay que notar que en el procedimiento que se explica a continuación sólo se están tratando los lados d ya, pero es exactamente el mismo para el resto de lados del hangar. Lo único que haría falta sería cambiar la notación propia del lado por el lado que se esté evaluando. d1 d2 dk d4 d5 d... d,(nx−1) a1 a2 a3 a4 aq a,... a,(ny−1) Lado d Lado a Vértice 1 Vértice 2 Vértice 3 Vértice 4 Lado b Lado c Ly Lx ~rtecho,d,k|global N E x’ y’ X Y ~rtecho,d,k|global ~rd,k|global ~rd,k|local ~rO0|global Ltecho ~esombra AzN ~e ~esombra Figura 5.15 Cálculo de la porción de techo necesaria. 1. Calcular el vector~rd,k|global(vector rojo de la figura) Este vector se define como ~rd,k|global =G~rd,k|local +~rO0|global (5.14) donde G es la matriz de cambio de base definida en la expresión 5.6 y ~rO0|global el vector para la traslación del cambio de base, definido en la expresión 5.7. En definitiva, lo que se ha hecho ha sido expresar en ejes globales la posición del punto k del lado 64 Capítulo 5. Obtención computacional de sombra usando Matlab Figura 5.24 Porción de techo tapado necesario para el 1 de Julio de 2021, a las 19.5h, hora local. 5.4 Consideraciones adicionales Ya se ha definido para cada día del año cuál es la necesidad de techo tapado para que garantice la sombra deseada en el interior del hangar. Sin embargo, es necesario tomar algunas consideraciones adicionales que mejoran el diseño. Todas las representaciones que se han hecho hasta ahora han sido necesidades de techo en planta, pero no hay que perder de vista que la cubierta del tejado real tiene forma de diente de sierra, y cuyo fenómeno puede jugar a nuestro favor si se usa de la forma apropiada. En efecto, gracias a este diseño, se determinan 3 zonas de funcionamiento del sistema dinámico de persianas. 1. Az(t)<(AzN+90°) Siempre y cuando el azimuth del Sol cumpla dicha condición, los laterales totalmente verticales de los dientes de sierra que se ven en la Figura 5.26 pueden permanecer abiertos, independientemente del cálculo de techo necesario. Es decir, todos los laterales verticales de los dientes de sierra permanecerán destapados hasta que Az<(AzN+90°) . En la Figura 5.26 se ha pintado en verde aquellos laterales que pueden quedar abiertos mientras se cumpla la condición, y en rojo el techo que se taparía en el caso de que el cálculo de porción de techo hecho en el apartado anterior lo requiriera. En la Figura 5.25 se ha marcado en amarillo toda la zona de acción en la que puede moverse el Sol para que se cumpla la condición. 5.4 Consideraciones adicionales 65 ZONA EN LA QUE SE CUMPLE QUE Az(t)<(AzN+90°) Figura 5.25 Zona en la que se cumple la condición para Az(t). Figura 5.26 En verde, la zona que quedaría abierta mientras se cumpla la condición. En rojo, la zona que se taparía si fuera necesario. 66 Capítulo 5. Obtención computacional de sombra usando Matlab También se ha señalado en la Figura 5.27 todo el rango de azimuth en cada día del año para el cual se garantizaría la condición. Puede verse que prácticamente la mitad del día permanecerían abiertos estos laterales. 19-Mar-2021 20-Jun-2021 21-Sep-2021 20-Dec-2021 Día del año 50 100 150 200 250 300 350 Azsol (t)[º] Azimuth del Sol para cada día del año, desde el 1 de Enero de 2021 X 14673180 Y 300.1228 X 14620680 Y 59.9527 X 6590220 Y 91.177 X 6633060 Y 269.1165 X 22703700 Y 269.1283 X 22660740 Y 90.6495 X 30475320 Y 239.796 X 30441480 Y 120.1764 X 47700 Y 190.392 Figura 5.27 En rojo, aquella zona donde se garantiza la condición del azimuth, para cada día del año. 2. Az(t)≥(AzN+90°)yh(t)≥hsierra Se define un nuevo ángulo, en este caso asociado a la geometría del diente de sierra de la cubierta del tejado. Teniendo en cuenta la Figura 2.16, podemos calcular el ángulo de los dientes de sierra de la cubierta respecto de la horizontal como: hsierra =tan119 −16 8=20.55°(5.20) Cuando el Sol tenga un azimuth mayor que el anterior y su elevación sea mayor que hsierra,todas las caras de los dientes de la cubierta deben estar tapadas, en caso de que el cálculo de techo lo exija. El Sol se encuentra en una posición en el cual todos los rayos de Sol inciden en todas las caras de la cubierta. 3. Az(t)≥(AzN+90°)yh(t)<hsierra En este caso, la situación se invierte. Ahora el Sol tiene una elevación menor que el ángulo que forman los dientes de sierra de la cubierta (ver Figura 5.28, por lo que la cara inclinada de dichos dientes puede permanecer abierta, con la garantía de que no entrarán los rayos de Sol directos (marcado en verde en la Figura 5.29). Lo que sí debe permanecer cerrado en todo momento son los lados verticales (marcados en rojo en la Figura 5.29), hasta donde indique el cálculo de porción de techo. 5.4 Consideraciones adicionales 67 hsierra h(t) Figura 5.28 Sol con una elevación inferior a hsierra. Figura 5.29 En verde, la zona que quedaría abierta mientras se cumpla la condición. En rojo, la zona que se taparía si fuera necesario. Por último, es interesante marcar dichas zonas de funcionamiento en la gráfica h(t)− Az(t) , Figura 5.30. Como ya se ha comentado, la zona de funcionamiento 1 es la más amplia de todas, ocupando más de la mitad del día. Por otro lado, la zona de funcionamiento 3 también permite en el atardecer abrir de forma segura las partes de la cubierta que ya se han comentado. Sin embargo, la más deseable de las dos es la zona de funcionamiento 1, no sólo porque tiene más tiempo de aplicación, sino porque parte desde el amanecer hasta horas centrales del día donde la intensidad de luz es mayor, y por tanto entraría más claridad. 68 Capítulo 5. Obtención computacional de sombra usando Matlab Por otro lado, la zona de funcionamiento 2 no permite abrir superficies con seguridad, sino que hay que atanerse a los resultados arrojados por el cálculo de porción de techo. Figura 5.30 Sombreado en amarillo la zona de funcionamiento 1 ( Az(t)<(AzN+90°) ); en verde la zona de funcionamiento 3 ( Az(t)≥(AzN+90°) y h(t)<hsierra );sin sombreado, la zona de funcionamiento 2 (Az(t)≥(AzN+90°)yh(t)≥hsierra) . 6 Conclusiones En este proyecto se ha evaluado con todo detalle y de forma personalizada la casuística que ocurre en cuanto al Sol, para el hangar de San Pablo Norte. Se han representado las ecuaciones que determinan la posición del Sol para la latitud y longitud en estudio, para cada instante del año. Gracias a eso, se ha podido analizar de forma cualitativa por donde se produce el amanecer y el atardecer, y cuál es la elevación máxima y mínima respecto del horizonte. Luego, se ha establecido un procedimiento de resolución para tratar de encontrar cuál es la sombra que las propias paredes del hangar arrojan sobre el interior del edificio. Esto era necesario para saber cuánto de techo efectivo hacía falta para generar sombra en aquellas partes donde las paredes no llegaban. Se ha podido ver en las Figuras 5.16-5.24 la comparativa que existe entre dos días totalmente opuestos en el ciclo orbital de la Tierra alrededor del Sol. Ambos días estaban muy cercanos a los Solsticios de Invierno y Verano, pudiendo dar una visión bastante fiel de cómo sería el comportamiento en los casos extremos. En esas figuras pudo verse que la diferencia es bastante notoria. En el caso más cercano al Solsticio de Invierno (2 de Enero), era poca la superficie del techo que necesitaba ser cubierta, puesto que el Sol tenía una elevación pequeña. En el caso más cercano al Solsticio de Verano (1 de Julio), ocurría justo al contrario. Desde muy temprano en el día hasta muy entrada la tarde era necesario tapar prácticamente la totalidad del techo. Sin embargo, volviendo a tener en cuenta el movimiento que realizad el Sol y la orientación que tiene la nave respecto del Norte, se llegó a determinar 3 zonas de funcionamiento del sistema dinámico de persianas, que complementan profundamente al cálculo de techo necesario. En esas zonas de funcionamiento se tenía en cuenta la geometría de dientes de sierra del techo, en el cual uno de los lados del diente es absolutamente vertical. Por esa razón, el Sol no incidiría sobre él hasta que no tuviera un azimuth suficientemente elevado como para que los rayos de Sol pudieran incidir sobre dicho lado. Afortunadamente, hasta llegar a dicho azimuth en el que se produce el cambio de funcionamiento, transcurre la mitad del día, desde el amanecer hasta las horas centrales, cuando hay más luz, lo cual es un hecho muy positivo. Otra de las zonas de funcionamiento aprovecha el ángulo que tiene los dientes de sierra de la cubierta del techo, permitiendo que el lado inclinado se abra libremente siempre que el Sol tenga una elevación inferior al ángulo que posee dicha cuña. En relación a todo lo expuesto, se sugieren diferentes vías de actuación: 1. Situación actual En la situación actual del hangar, todo el techo interior está cubierto por las planchas 69 70 Capítulo 6. Conclusiones opacas que se ha comentado, aunque no todas ellas están dispuestas para bloquear la luz directa. Viendo el hangar desde fuera con las imágenes desde el satélite, puede verse que sólo las caras verticales de los dientes de sierra del techo tienen cristales. Las caras inclinadas de dichos dientes están cubiertas de chapa, por lo que, aunque se quitaran las planchas opacas interiores de esos lados, no podría entrar luz. Figura 6.1 Puede verse a través de las cristaleras la cara interior de las planchas opacas que apoyan sobre los cristales. Por tanto, se sugiere para esta situación: • Colocar persianas domotizadas sólo en los laterales por donde puede entrar la luz. Permanecerán abiertas mientras el Sol se encuentre en la zona de funcionamiento 1, y luego se cerrarán a conveniencia según el cálculo de techo. Se obtendrá una solución parcial al estudio realizado, aunque si bien será efectivo gracias al rango amplio de la zona de funcionamiento 1. Esta sería la solución más barata. 2. Situación Óptima En caso de que se dispongan de los recursos económicos y se quiera acometer el proyecto en su totalidad, sería necesario: a) Determinar si las cubiertas de chapa constituyen un elemento estructural del hangar, y en caso de que no lo sean, estudiar si pueden ser sustituidas por algún elemento translúcido. b) En caso de que se pudieran sustituir, instalar persianas domotizadas en ambas caras del diente de sierra, tanto en el vertical como el inclinado. El funcionamiento del hangar será pues el que se ha tratado en el proyecto. 3. Persianas domotizadas Existen numerosas empresas dedicadas a la instalación de persianas, y aún más 71 en auge aquellas que permiten un control electrónico absoluto sobre ellas. Sin embargo, el control de cada una de las persianas podría realizarse de forma sencilla y no muy costosa a través de electrónica de software libre, como puede ser Arduino o Raspberry Pi. En cualquiera de los dos casos, ambos deben: • Permitir que el funcionamiento de apertura y cierre de las persianas venga determinado por los tiempos deducidos en las ecuaciones. Pueden almacenar a modo de tabla los instantes de todo el año de apertura y cierre, o si tiene suficiente capacidad de computación, que lo calcule a través de las ecuaciones. Por tanto, la variable de entrada que necesitarán el microcontrolador o microprocesador será únicamente la variable tiempo. •Deberán permitir la apertura/cierre manual en caso de fallo. •Ventajas e inconvenientes de usar microcontroladores o proveedor: –Ventajas: ∗Microcontroladores: Más barato. ∗ Microcontroladores: Máxima flexibilidad. Diseñado para cumplir exactamente con lo requerido. ∗Proveedor: son expertos en la materia. –Inconvenientes: ∗ Microcontroladores: Es necesario tiempo de desarrollo para programar todo el sistema. ∗Proveedor: Más costoso. ∗ Proveedor: Probablemente use un desarrollo comercial que no permita flexibilidad a la hora de cumplir con lo requerido. • Las persianas serán preferentemente enrollables, aislantes, que discurran suavemente por los raíles sin hacer ruido y sin trabarse. Las dimensiones serán las más apropiadas según indique el fabricante. •El mantenimiento de las persianas deberá ser accesible y sencillo. • Las persianas deberán cerrarse tiempo antes de la necesidad establecida por el cálculo de techo, para asegurar que la persiana esté cerrada cuando el Sol esté en la posición adecuada para entrar. 4. Luces interiores de intensidad dinámica Puesto que el techo estará en un estado de continuo cambio, es necesario que la iluminación interior esté monitorizada por sensores de iluminación, que se traduzca en generar la intensidad adecuada en las luces interiores y con la suficiente rapidez como para que en todo momento se garanticen en el interior las condiciones mínimas de iluminación. Apéndice A Código de Matlab La forma de utilizar el código es como sigue: ubicar todos los códigos que se adjuntan aquí en archivos .m diferentes y en la misma carpeta. Luego, debe saber que el código A.2 es un código listo para ejecutar, y depende de A.1 y A.3. Por otro lado, A.4 también es un código listo para ejecutar, y depende de A.1, A.5 y A.6 A.1 Posición del Sol: Cálculo de Az,h 1% clc; clear all; clf; close all; 2function [Az_fix,h_sol_fix, T]=PosicionSol(n) % n representa ... el paso en 3% minutos. 4 5% POSICIÓN GEOGRÁFICA DE LA NAVE 6phi_Sev=37.423292*pi/180; % rad 7lambda_Sev=-5.883221*pi/180; % rad 8 9% DATOS ORBITALES Y PLANETARIOS 10 epsilon=23.5*pi/180; 11 mu_sol=132712439935.5; % km^3/s^2 12 L_Tierra=149.598e6; % km 13 e=0.0167; 14 w_argPerigeo=114.2*pi/180; 15 16 17 % CÁLCULO DE GST0 USANDO DÍAS JULIANOS 18 % Día del Equinoccio de Primavera: 20 Marzo 2021, 09:37 UT 19 A=2021; M=1; D=1; 20 % Día Juliano a las 0:00 UT 21 J_0= 367*A-floor((7*A+7*floor((M+9)/12))/4)+floor(275*M/9)+D+ ... 1721013.5; 22 % Cálculo de la centuria Juliana 23 T_0=(J_0-2451545)/36525; 24 GST0=(100.4606184+36000.77004*T_0+0.000387933*T_0^2-2.583e- ... 8*T_0^3)*pi/180; 25 GST0 = mod(GST0, 2*pi); 26 w_Tierra=0.004178*pi/180; % rad/s 73 80 Capítulo A. Código de Matlab 150 L_abierto(i,j)=(H_Pared-Hs)/tan(h_sol(i,j)); 151 152 % Dirección unitaria del sol, en ejes globales 153 e_sol(i,j,1)=-sin(Az_sol(i,j)); % En la primera capa ... de la matriz tridimensional se almacena la ... coordenada X del vector e_sol 154 e_sol(i,j,2)=-cos(Az_sol(i,j)); % En la segunda capa ... de la matriz tridimensional se almacena la ... coordenada Y del vector e_sol 155 156 % Techo lado a 157 for k=1:length(r_k_local_a(:,1)) 158 techo=(C_a(i,j,k)-L_abierto(i,j)); 159 if techo≥0 160 % En ejes absolutos 161 r_techo_k_a_X(i,j,k)=r_k_a(1,k)+ ... techo*e_sol(i,j,1); 162 r_techo_k_a_Y(i,j,k)=r_k_a(2,k)+ ... techo*e_sol(i,j,2); 163 164 % En ejes locales 165 r_techo_k_a_local=G_la'*([r_techo_k_a_X(i,j,k); ... r_techo_k_a_Y(i,j,k)]-despl_local); 166 r_techo_k_a_X_local(i,j,k)=r_techo_k_a_local(1); 167 r_techo_k_a_Y_local(i,j,k)=r_techo_k_a_local(2); 168 169 end 170 end 171 172 % Techo lado b 173 for k=1:length(r_k_local_b(:,1)) 174 techo=(C_b(i,j,k)-L_abierto(i,j)); 175 if techo≥0 176 % En ejes absolutos 177 r_techo_k_b_X(i,j,k)=r_k_b(1,k)+ ... techo*e_sol(i,j,1); 178 r_techo_k_b_Y(i,j,k)=r_k_b(2,k)+ ... techo*e_sol(i,j,2); 179 180 % En ejes locales 181 r_techo_k_b_local=G_la'*([r_techo_k_b_X(i,j,k); ... r_techo_k_b_Y(i,j,k)]-despl_local); 182 r_techo_k_b_X_local(i,j,k)=r_techo_k_b_local(1); 183 r_techo_k_b_Y_local(i,j,k)=r_techo_k_b_local(2); 184 185 end 186 end 187 188 % Techo lado d 189 for k=1:length(r_k_local_d(:,1)) 190 techo=(C_d(i,j,k)-L_abierto(i,j)); 191 if techo≥0 192 % En ejes absolutos 193 r_techo_k_d_X(i,j,k)=r_k_d(1,k)+ ... techo*e_sol(i,j,1); A.2 Representación de sombra y techo 81 194 r_techo_k_d_Y(i,j,k)=r_k_d(2,k)+ ... techo*e_sol(i,j,2); 195 196 % En ejes locales 197 r_techo_k_d_local=G_la'*([r_techo_k_d_X(i,j,k); ... r_techo_k_d_Y(i,j,k)]-despl_local); 198 r_techo_k_d_X_local(i,j,k)=r_techo_k_d_local(1); 199 r_techo_k_d_Y_local(i,j,k)=r_techo_k_d_local(2); 200 201 end 202 end 203 204 % Techo lado c 205 for k=1:length(r_k_local_c(:,1)) 206 techo=(C_c(i,j,k)-L_abierto(i,j)); 207 if techo≥0 208 % En ejes absolutos 209 r_techo_k_c_X(i,j,k)=r_k_c(1,k)+ ... techo*e_sol(i,j,1); 210 r_techo_k_c_Y(i,j,k)=r_k_c(2,k)+ ... techo*e_sol(i,j,2); 211 212 % En ejes locales 213 r_techo_k_c_local=G_la'*([r_techo_k_c_X(i,j,k); ... r_techo_k_c_Y(i,j,k)]-despl_local); 214 r_techo_k_c_X_local(i,j,k)=r_techo_k_c_local(1); 215 r_techo_k_c_Y_local(i,j,k)=r_techo_k_c_local(2); 216 end 217 end 218 219 220 221 hold on; 222 % Movimiento del Sol 223 figure(1); ... plot(r_sol_X(i,j),r_sol_Y(i,j),'*y','LineWidth',3); 224 % Indicamos la fecha y hora del día 225 dim = [0.63 0.72 0.1 0.1]; 226 str = {['HORA LOCAL = ... ',num2str(1+(T(i,j)-T(1,j))/3600), ' h'],['DIA= ... ',datestr(datetime(2021,1,j+1))],['h_{sol}= ... ',num2str(h_sol(i,j)*180/pi), 'º']}; 227 Caja=annotation('textbox',dim, 'String', str, ... 'FitBoxToText','on'); 228 sz=Caja.BackgroundColor; 229 Caja.BackgroundColor='white'; 230 %% SOMBRA POR LAS PAREDES 231 % Sombra del vértice 1 y 2 232 if Az_sol(i,j)<Az_Nave+pi/2 233 Sombra_a_X=[r_1(1) r_2(1) r_sombra_2_X(i,j) ... r_sombra_1_X(i,j)]; 234 Sombra_a_Y=[r_1(2) r_2(2) r_sombra_2_Y(i,j) ... r_sombra_1_Y(i,j)]; 235 fill(Sombra_a_X,Sombra_a_Y,[0.44 0.8 0.8]); ... alpha(0.5); 236 end 82 Capítulo A. Código de Matlab 237 % Sombra del vértice 2 y 3 238 if Az_sol(i,j)≥Az_Nave && Az_sol(i,j)≤Az_Nave+pi 239 Sombra_b_X=[r_2(1) r_3(1) r_sombra_3_X(i,j) ... r_sombra_2_X(i,j)]; 240 Sombra_b_Y=[r_2(2) r_3(2) r_sombra_3_Y(i,j) ... r_sombra_2_Y(i,j)]; 241 fill(Sombra_b_X,Sombra_b_Y,[0.44 0.8 0.8]); ... alpha(0.5); 242 end 243 % Sombra del vértice 3 y 4 244 if Az_sol(i,j)≥Az_Nave+pi/2 && ... Az_sol(i,j)≤Az_Nave+3*pi/2 245 Sombra_c_X=[r_3(1) r_4(1) r_sombra_4_X(i,j) ... r_sombra_3_X(i,j)]; 246 Sombra_c_Y=[r_3(2) r_4(2) r_sombra_4_Y(i,j) ... r_sombra_3_Y(i,j)]; 247 fill(Sombra_c_X,Sombra_c_Y,[0.44 0.8 0.8]); ... alpha(0.5); 248 end 249 250 % Sombra del vértice 4 y 1 251 if Az_sol(i,j)≥Az_Nave+pi || Az_sol(i,j)≤Az_Nave 252 Sombra_d_X=[r_4(1) r_1(1) r_sombra_1_X(i,j) ... r_sombra_4_X(i,j)]; 253 Sombra_d_Y=[r_4(2) r_1(2) r_sombra_1_Y(i,j) ... r_sombra_4_Y(i,j)]; 254 fill(Sombra_d_X,Sombra_d_Y,[0.44 0.8 0.8]); ... alpha(0.5); 255 end 256 257 258 %% REPRESENTACIÓN DEL TECHO 259 % Representación de los puntos de las paredes 260 plot(r_k_d(1,:),r_k_d(2,:),'o'); 261 plot(r_k_a(1,:),r_k_a(2,:),'o'); 262 plot(r_k_c(1,:),r_k_c(2,:),'o'); 263 plot(r_k_b(1,:),r_k_b(2,:),'o'); 264 % Techo d 265 if abs(sum(r_techo_k_d_X(i,j,:)))>0 || ... abs(sum(r_techo_k_d_Y(i,j,:)))>0 266 for k=1:length(r_k_d(1,:)) % Recorremos todos ... los puntos k de la pared d 267 if r_techo_k_d_X(i,j,k)6=0||... r_techo_k_d_Y(i,j,k)6=0 268 plot([r_k_d(1,k) r_techo_k_d_X(i,j,k)], ... [r_k_d(2,k) r_techo_k_d_Y(i,j,k)],'b'); 269 end 270 end 271 end 272 % Techo a 273 if abs(sum(r_techo_k_a_X(i,j,:)))>0 || ... abs(sum(r_techo_k_a_Y(i,j,:)))>0 274 for k=1:length(r_k_a(1,:)) 275 if r_techo_k_a_X(i,j,k)6=0||... r_techo_k_a_Y(i,j,k)6=0 A.2 Representación de sombra y techo 83 276 plot([r_k_a(1,k) r_techo_k_a_X(i,j,k)], ... [r_k_a(2,k) r_techo_k_a_Y(i,j,k)],'b'); 277 end 278 end 279 end 280 % Techo b 281 if abs(sum(r_techo_k_b_X(i,j,:)))>0 || ... abs(sum(r_techo_k_b_Y(i,j,:)))>0 282 for k=1:length(r_k_b(1,:)) 283 if r_techo_k_b_X(i,j,k)6=0||... r_techo_k_b_Y(i,j,k)6=0 284 plot([r_k_b(1,k) r_techo_k_b_X(i,j,k)], ... [r_k_b(2,k) r_techo_k_b_Y(i,j,k)],'b'); 285 end 286 end 287 end 288 % Techo c 289 if abs(sum(r_techo_k_c_X(i,j,:)))>0 || ... abs(sum(r_techo_k_c_Y(i,j,:)))>0 290 for k=1:length(r_k_c(1,:)) 291 if r_techo_k_c_X(i,j,k)6=0||... r_techo_k_c_Y(i,j,k)6=0 292 plot([r_k_c(1,k) r_techo_k_c_X(i,j,k)], ... [r_k_c(2,k) r_techo_k_c_Y(i,j,k)],'b'); 293 end 294 end 295 end 296 alpha(0.5); hold off; pause(0.4); 297 end 298 end 299 end 84 Capítulo A. Código de Matlab A.3 Cmax 1function [C_a,r_local_a,C_b,r_local_b,C_c,r_local_c, C_d, ... r_local_d]=C_Nave_Toda_longitud(Az_Nave, Lx, Ly, nx, ny, ... Az_sol, h_sol, hmin, alpha1, alpha2, alpha3, alpha4) 2%% Definición de los lados respecto del sistema de referencia local 3% Lado d 4r_local_d(1:nx+1,1)=Lx-(Lx/nx)*(0:nx); 5r_local_d(1:nx+1,2)=Ly; 6% Lado a 7r_local_a(1:ny+1,1)=0; 8r_local_a(1:ny+1,2)=Ly-(Ly/ny)*(0:ny); 9% Lado b 10 r_local_b(1:nx+1,1)=(Lx/nx)*(0:nx); 11 r_local_b(1:nx+1,2)=0; 12 % Lado c 13 r_local_c(1:ny+1,1)=Lx; 14 r_local_c(1:ny+1,2)=(Ly/ny)*(0:ny); 15 16 %% 17 alpha43=Az_Nave-pi/2; 18 alpha12=alpha43; 19 [f,c]=size(Az_sol); 20 C_a=zeros(f,c,ny); 21 C_b=zeros(f,c,nx); 22 C_c=zeros(f,c,ny); 23 C_d=zeros(f,c,nx); 24 25 for j=1:length(Az_sol(1,:)) 26 for i=1:length(Az_sol(:,1)) 27 if h_sol(i,j)>0 && h_sol(i,j)≥hmin 28 %% CASO Az_sol EN EL CUADRANTE QUE CONTIENE AL ... VÉRTICE 1 29 % Az_sol menor que alpha1 30 if Az_sol(i,j)<alpha1 && Az_sol(i,j)≥Az_Nave-pi/2 31 C3max=Ly/cos(Az_sol(i,j)-alpha43); 32 X_cambio_d=C3max*sin(Az_sol(i,j)-alpha43); 33 % LADOS DE INTERÉS: A Y D 34 % Lado d 35 for k=1:length(r_local_d(:,1)) 36 if r_local_d(k,1)<X_cambio_d 37 C_d(i,j,k)=Lx*k/(nx*sin(Az_sol(i,j)- ... alpha43)); 38 else 39 C_d(i,j,k)=C3max; 40 end 41 end 42 43 % Lado a 44 for k=1:length(r_local_a(:,1)) 45 C_a(i,j,k)=r_local_a(k,2)/cos(Az_sol(i,j)- ... alpha12); 46 end 47 end A.3 Cmax 85 48 49 % Az_sol entre alpha1 y Az_Nave 50 if Az_sol(i,j)<Az_Nave && Az_sol(i,j)≥alpha1 51 C1max=Lx/cos(Az_Nave-Az_sol(i,j)); 52 Y_cambio_a=C1max*sin(Az_Nave-Az_sol(i,j)); 53 % LADOS DE INTERÉS: A Y D 54 % Lado d 55 for k=1:length(r_local_d(:,1)) 56 C_d(i,j,k)=(Lx-r_local_d(k,1))/cos(Az_Nave- ... Az_sol(i,j)); 57 end 58 % Lado a 59 for k=1:length(r_local_a(:,1)) 60 if r_local_a(k,2)≥Y_cambio_a 61 C_a(i,j,k)=C1max; 62 else 63 C_a(i,j,k)=r_local_a(k,2)/sin(Az_Nave- ... Az_sol(i,j)); 64 end 65 end 66 end 67 68 %% CASO Az_sol EN EL CUADRANTE QUE CONTIENE AL ... VÉRTICE 2 69 % Az_sol menor que alpha2 70 if Az_sol(i,j)<alpha2 && Az_sol(i,j)≥Az_Nave 71 C2max=Lx/cos(Az_sol(i,j)-Az_Nave); 72 Y_cambio_a=Ly-C2max*sin(Az_sol(i,j)-Az_Nave); 73 % LADOS DE INTERÉS: A Y B 74 % Lado a 75 for k=1:length(r_local_a(:,2)) 76 if r_local_a(k,2)≥Y_cambio_a 77 C_a(i,j,k)=(Ly-r_local_a(k,2))/sin( ... Az_sol(i,j)-Az_Nave); 78 else 79 C_a(i,j,k)=C2max; 80 end 81 end 82 83 % Lado b 84 for k=1:length(r_local_b(:,1)) 85 C_b(i,j,k)=(Lx-r_local_b(k,1))/cos( ... Az_sol(i,j)-Az_Nave); 86 end 87 end 88 89 % Az_sol mayor que alpha2 90 if Az_sol(i,j)≥alpha2 && Az_sol(i,j)<Az_Nave+pi/2 91 C2max=Ly/sin(Az_sol(i,j)-Az_Nave); 92 X_cambio_b=Lx-C2max*cos(Az_sol(i,j)-Az_Nave); 93 % LADOS DE INTERÉS: A Y B 94 % Lado a 95 for k=1:length(r_local_a(:,2)) 96 C_a(i,j,k)=(Ly-r_local_a(k,2))/sin( ... Az_sol(i,j)-Az_Nave); 97 end 86 Capítulo A. Código de Matlab 98 99 % Lado b 100 for k=1:length(r_local_b(:,1)) 101 if r_local_b(k,1)>X_cambio_b 102 C_b(i,j,k)=(Lx-r_local_b(k,1))/cos( ... Az_sol(i,j)-Az_Nave); 103 else 104 C_b(i,j,k)=C2max; 105 end 106 end 107 end 108 %% CASO Az_sol EN EL CUADRANTE QUE CONTIENE AL ... VÉRTICE 3 109 % Az_sol menor que alpha3 110 if Az_sol(i,j)<alpha3 && Az_sol(i,j)≥Az_Nave+pi/2 111 beta=pi-(Az_sol(i,j)-Az_Nave); 112 C3max=Ly/sin(beta); 113 X_cambio_b=C3max*cos(beta); 114 % LADOS DE INTERÉS: B Y C 115 % Lado b 116 for k=1:length(r_local_b(:,1)) 117 if r_local_b(k,1)<X_cambio_b 118 C_b(i,j,k)=r_local_b(k,1)/cos(beta); 119 else 120 C_b(i,j,k)=C3max; 121 end 122 end 123 124 % Lado c 125 for k=1:length(r_local_c(:,2)) 126 C_c(i,j,k)=(Ly-r_local_c(k,2))/sin(beta); 127 end 128 end 129 % Az_sol mayor que alpha3 130 if Az_sol(i,j)≥alpha3 && Az_sol(i,j)<Az_Nave+pi 131 beta=Az_sol(i,j)-Az_Nave-pi/2; 132 C3max=Lx/sin(beta); 133 Y_cambio_c=Ly-C3max*cos(beta); 134 % LADOS DE INTERÉS: B Y C 135 % Lado b 136 for k=1:length(r_local_b(:,1)) 137 C_b(i,j,k)=r_local_b(k,1)/sin(beta); 138 end 139 140 % Lado c 141 for k=1:length(r_local_c(:,2)) 142 if r_local_c(k,2)>Y_cambio_c 143 C_c(i,j,k)=(Ly-r_local_c(k,2))/cos(beta); 144 else 145 C_c(i,j,k)=C3max; 146 end 147 end 148 end 149 150 %% CASO Az_sol EN EL CUADRANTE QUE CONTIENE AL ... VÉRTICE 4 A.3 Cmax 87 151 % Az_sol menor que alpha4 152 if Az_sol(i,j)<alpha4 && Az_sol(i,j)≥Az_Nave+pi 153 beta=-Az_sol(i,j)+Az_Nave+3*pi/2; 154 C4max=Lx/sin(beta); 155 Y_cambio_c=C4max*cos(beta); 156 % LADOS DE INTERÉS: C Y D 157 % Lado c 158 for k=1:length(r_local_c(:,2)) 159 if r_local_c(k,2)<Y_cambio_c 160 C_c(i,j,k)=r_local_c(k,2)/cos(beta); 161 else 162 C_c(i,j,k)=C4max; 163 end 164 end 165 % Lado d 166 for k=1:length(r_local_d(:,1)) 167 C_d(i,j,k)=(Lx-r_local_d(k,1))/sin(beta); 168 end 169 end 170 end 171 end 172 end 88 Capítulo A. Código de Matlab A.4 Representación posición del Sol 1clc; clear all; 2clf; close all; 3%%%% Representación de Az_Sol y h_Sol %%%% 4n=1; 5[Az_fix,h_sol_fix, T]=PosicionSol(n); 6 7inicio=0; 8final=0; 9% REPRESENTACIÓN DE TODOS LOS DÍAS DEL AÑO 10 for i=1:length(T(1,:)) 11 % Todos los días del año 12 % Elevación del Sol para cada día del año 13 figure(1); plot(T(:,i),h_sol_fix(:,i)*180/pi,'b'); hold on; 14 % Azimuth del Sol para cada día del año 15 for t=2:length(T(:,1)) 16 if Az_fix(t,i)==0 17 if Az_fix(t-1,i)6=0 18 if i==79 && t==1115 19 final=t-2; 20 else 21 final=t-1; 22 end 23 if inicio6=0 24 break 25 end 26 else 27 continue 28 end 29 else 30 if inicio==0 31 inicio=t; 32 end 33 end 34 end 35 figure(3);plot(T(inicio:final,i), ... Az_fix(inicio:final,i)*180/pi, 'b'); hold on; 36 inicio=0; final=0; 37 end 38 %% ELECCIÓN DE COLORES 39 paso=30; %Paso de días 40 C=[1 0 0;0 0 1]; 41 Saturacion=0.5; 42 Brillo=1; 43 NumeroCurvas=2+4+floor(length(T(1,:))/paso); %+2 porque se ... representa el día 1 y el día 365; +4 porque se representan ... los solsticios y los equinoccios 44 Colores=CambiaColor(C,NumeroCurvas,Saturacion,Brillo); 45 46 %% REPRESENTACIÓN SÓLO DE AQUELLOS DÍAS DEL AÑO QUE INDICA LA ... ETIQUETA 47 k=1; 48 for i=[1 paso:paso:paso*floor(length(T(1,:))/paso) length(T(1,:))] A.4 Representación posición del Sol 89 49 50 % Elevación del Sol para cada día i del año 51 figure(5); 52 plot(T(:,i),h_sol_fix(:,i)*180/pi,'color',Colores(k,:), ... 'LineWidth', 1.75); hold on;%pause; 53 54 % Azimuth del Sol para cada día del año 55 [inicio,final]=inicio_final_dia(T,Az_fix,i); 56 figure(6);plot(T(inicio:final,i), ... Az_fix(inicio:final,i)*180/pi, 'color',... Colores(k,:),'LineWidth',1.75); hold on; 57 58 % En un solo día de 24h 59 figure(7); h_Sol_Unitario=plot(T(:,1), ... h_sol_fix(:,i)*180/pi, 'color',... Colores(k,:),'LineWidth',1.75); hold on; 60 h_Sol_Unitario_legend{k}=['Día: ',datestr(datetime(2021,1,i))]; 61 [inicio,final]=inicio_final_dia(T,Az_fix,i); 62 figure(8); ... plot(T(inicio:final,1),Az_fix(inicio:final,i)*180/pi, ... 'color',Colores(k,:),'LineWidth',1.75); hold on; 63 Az_Sol_Unitario_legend{k}=['Día: ... ',datestr(datetime(2021,1,i))]; 64 k=k+1; 65 end 66 67 %% REPRESENTACIÓN DE LOS EQUINOCCIOS Y SOLSTICIOS 68 % Equinoccio de Primavera 2021: 20 de Marzo 69 Equi_Sols(1)=caldays(between(datetime(2021,1,1), ... datetime(2021,3,20), 'days'))-1; 70 % Solsticio de Verano: 21 de Junio 71 Equi_Sols(2)=caldays(between(datetime(2021,1,1), ... datetime(2021,6,21),'days'))-1; 72 % Equinoccio de Otoño: 22 de Septiembre 73 Equi_Sols(3)=caldays(between(datetime(2021,1,1), ... datetime(2021,9,22), 'days'))-1; 74 % Solsticio de Invierno: 21 de Diciembre 75 Equi_Sols(4)=caldays(between(datetime(2021,1,1), ... datetime(2021,12,21), 'days'))-1; 76 77 for i=1:4 78 % Elevación del Sol para cada día i del año 79 figure(1); 80 plot(T(:,Equi_Sols(i)),h_sol_fix(:,Equi_Sols(i))*180/pi,'.-',... 'color', Colores(k,:), 'LineWidth',1); hold on;%pause; 81 figure(5); 82 plot(T(:,Equi_Sols(i)),h_sol_fix(:,Equi_Sols(i))*180/pi,'.-',... 'color', Colores(k,:), 'LineWidth',1); hold on;%pause; 83 84 85 % Azimuth del Sol para cada día del año 86 figure(2);plot(T(:,Equi_Sols(i)),Az_fix(:,Equi_Sols(i))*180/pi, ... '.-','color', Colores(k,:), 'LineWidth',1); hold on; 87 [inicio,final]=inicio_final_dia(T,Az_fix,Equi_Sols(i)); 96 Capítulo A. Código de Matlab A.6 Función auxiliar: Colores graduales para las gráficas 1function Colores=CambiaColor(C,n,Saturacion,Brillo) 2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3%%%%%%%%FUNCIÓN PARA GENERAR COLORES CORRELATIVOS A ... CONVENIENCIA %%%%%%%%% 4%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 5%% ARGUMENTOS DE ENTRADA ... %% 6%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 7% C: Vector en el que se incluye el color inicial y final 8% en su forma RGB entre los que se va a interpolar para 9% obtener el resto de colores. 10 % 11 % n: Número de colores que queremos obtener. Si los colores ... inicial y 12 % final son muy parecidos y además n es muy grande, la ... diferencia entre 13 % colores será imperceptible. 14 % 15 % Saturacion: Número entre 0 y 1. Si no se incluye, se toma el ... valor 1 por 16 % defecto. 17 % 18 % Brillo: Número entre 0 y 1. Si no se incluye, se toma el ... valor 1 por 19 % defecto. 20 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 21 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 22 %% ARGUMENTOS DE SALIDA ... %% 23 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 24 % Colores: Vector nx3 con los colores interpolados entre ... C(1,:) y C(2,:), 25 % en formato RGB, esta vez RGB [0, 1] para que "plot" lo reconozca. 26 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 27 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 28 %% EJEMPLO ... %% 29 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 30 % [X,Y,Z] = peaks(30); 31 % surfc(X,Y,Z) 32 % axis([-3 3 -3 3 -10 5]) 33 % 34 % C=[1 0 0;0 0 1]; %Color Inicial: Rojo; Color Final: Azul. 35 % Saturacion=0.5; 36 % Brillo=1; 37 % n=20; 38 % Colores=CambiaColor(C,n,Saturacion,Brillo); 39 % colormap(Colores) 40 % colorbar 41 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 42 43 A.6 Función auxiliar: Colores graduales para las gráficas 97 44 45 if nargin==2 46 Saturacion=1; Brillo=1; 47 elseif nargin==3 48 Brillo=1; 49 end 50 51 52 % Comprobamos que C esté comprendido entre [0, 1] 53 m=max(max(C)); 54 if m>1 55 C=C/255; 56 end 57 58 % Convertimos de RGB al formato HSV. Este formato está formado ... por las 59 % siglas H = Hue (Matiz o Color), S = Saturation, V = Valor ... (Brillo). 60 % Este formato devolverá una fila y tres columnas (por cada ... color RGB), 61 % donde cada elemento estará en [0, 1]. La primera columna ... corresponde al 62 % Matiz, la segunda a la Saturación y la tercera al Brillo. 63 C_HSV = rgb2hsv(C); 64 65 % Interpolamos el matiz (C_HSV(:,1)) entre el color inicial ... (que correspon- 66 % de al punto 0) y el color final (que corresponde al punto n), ... para n 67 % puntos. 68 C_HSV_interp = interp1([1 n], C_HSV(:, 1), 1:n); 69 70 % Una vez interpolado el matiz de cada color para los puntos n ... elegidos, le 71 % incorporamos a cada matiz la Saturación y el Brillo elegidos. ... En la 72 % siguiente matriz, la primera columna es el matiz interpolado. ... La segunda 73 % y tercera columna serán la saturación y el brillo, ... respectivamente. La 74 % orden repmat(A,m,n) lo que hace es tratar la matriz A como si ... fuera el 75 % elemento de una matriz, de modo que forma una "matriz" mxn ... donde cada 76 % elemento es la propia matriz A. 77 if length(Saturacion)==1 && length(Brillo)==1 78 C_HSV = [C_HSV_interp(:), repmat([Saturacion Brillo], n, 1)]; 79 elseif length(Saturacion)==1 && length(Brillo)>1 80 C_HSV = [C_HSV_interp(:), ones(n,1)*Saturacion, interp1([1 ... n], Brillo,1:n)']; 81 elseif length(Saturacion)>1 && length(Brillo)==1 82 C_HSV = [C_HSV_interp(:), interp1([1 n],Saturacion,1:n)', ... ones(n,1)*Brillo]; 83 elseif length(Saturacion)>1 && length(Brillo)>1 84 C_HSV = [C_HSV_interp(:), interp1([1 n],Saturacion,1:n)', ... interp1([1 n], Brillo,1:n)']; 98 Capítulo A. Código de Matlab 85 end 86 87 % Convertimos los colores de HSV a RGB 88 Colores = hsv2rgb(C_HSV); Bibliografía [1] Google Earth Pro ,https://earth.google.com/static/9.130.0.0/app_min__es.html, Version: 7.3.3.7786. [2] Nicolás Fernando Marrugo Cardenas, Tfg: Control de una fachada dinámica de 2 grados de libertad con implementación de paneles solares , http://repository.unipiloto.edu.co/bitstream/handle/20.500.12277/1070/ 00001752.pdf?sequence=1&isAllowed=y, 2014, Universidad Piloto de Colombia. Facultad Ingeniería Mecatrónica. Bogotá D.C. [3] Josué Alberto Alonzo Chavarría, Tesis: Diseño del control automático de un sistema híbrido de iluminación , 2012, Universidad Autónoma de Yucatán.Facultad de Ingeniería. [4] Howard Curtis, Orbital mechanics for engineering students , Elsevier Butterworth- Heinemann, 2005. [5] Alberto Barres Martínez, Tfg: Diseño e implementación de un sistema de control de iluminación natural basado en arduino , 2018, Escola Tècnica Superior dÉnginyeria Informàtica. Universitat Politècnica de València. [6] National Oceanic and Atmospheric Administration (NOAA), Equinoccios y Solsticios, serie 2020-2025 ,https://www.weather.gov/dvn/Climate_ Astronomical_Seasons. [7] Quirón Prevención, Trabaja con iluminación perfecta, 5 de mayo de 2017 ,https:// www.quironprevencion.com/blogs/es/prevenidos/trabaja-iluminacion-perfecta. [8] A. A. Rizvi, K. Addoweesh, A. El-Leathy, and H. Al-Ansary, Sun position algorithm for sun tracking applications , IECON 2014 - 40th Annual Conference of the IEEE Industrial Electronics Society, 2014, pp. 5595–5598. [9] Sorruno, Dibujo de avión C295 ,https://upload.wikimedia.org/wikipedia/commons/ d/d3/EADS_CASA_C-295.svg, Fecha: 13-Abril-2009. [10] Rafael Vázquez Valenzuela, Apuntes Mecánica Orbital/Astronáutica y Vehículos Espaciales. Tema 1: Introducción ,http://www.aero.us.es/move/files/Tema1.pdf , Septiembre 2015, Departamento de Ingeniería Aeroespacial. Escuela Superior de Ingenieros. Universidad de Sevilla. 99 100 Bibliografía [11] , Apuntes Mecánica Orbital/Astronáutica y Vehículos Espaciales. Tema 2: Mecánica Orbital Básica ,http://www.aero.us.es/move/files/t2.pdf , Septiembre 2015. [12] , Apuntes Mecánica Orbital/Astronáutica y Vehículos Espaciales. Tema 3: Análisis y Diseño de Misiones Geocéntricas ,http://aero.us.es/move/files/t3p1.pdf , Marzo 2015, Departamento de Ingeniería Aeroespacial. Escuela Superior de Ingenieros. Universidad de Sevilla. [13] , Formulario ,http://aero.us.es/move/files/formulario1516.pdf , 2015, Departamento de Ingeniería Aeroespacial. Escuela Superior de Ingenieros. Universidad de Sevilla. [14] , Apuntes Mecánica Orbital y Vehículos Espaciales. Tema 5: Trigonometría Esférica. Círculos Esféricos. Aplicaciones: traza, cobertura, visibilidad. ,http://aero. us.es/move/files1920/Tema5_1920.pdf , Septiembre 2019, Departamento de Ingeniería Aeroespacial. Escuela Superior de Ingenieros. Universidad de Sevilla. [15] Esteban Villarejo, FAL C295 en Airbus San Pablo Sur ,https://static-abcblogs.abc. es/wp-content/uploads/sites/46/2017/11/FALC295-2-516x315.jpg, Fecha: 15-Nov- 2017.