Programa para el cálculo de los coeficientes de un modelo de apuntado para un telescopio. Leonel Gutiérrez

Octubre, 2005

1. Introducción Generalmente, cuando se instala un telescopio, su alineación no es perfecta. Esta falta de alineación se traduce siempre en errores, tanto en el apuntado como en el seguimiento de los objetos estelares con el telescopio. Por tanto, en un telescopio profesional es necesario alinearlo de la mejor manera posible, moviendo la montura del telescopio hasta que los errores se reducen al mínimo. Esta tarea se ha realizado ya en los telescopios del Observatorio Astronómico Nacional en San Pedro Mártir. Sin embargo, aún en esas circunstancias, quedan errores residuales muy difíciles de corregir mecánicamente, por lo que es necesario contar con un modelo que permita corregir previamente las coordenadas para conseguir un óptimo apuntado del telescopio. Aunque aquí nos referiremos sólo a los telescopios con montura ecuatorial, la idea también se aplica a otros tipos de monturas. Actualmente, los telescopios del Observatorio Astronómico Nacional en San Pedro Mártir cuentan con un modelo que hace un ajuste en base a polinomios. Sin embargo, debido a diversos factores como el hecho de que la distribución de cargas del telescopio no sea constante en el tiempo, pues están sujetos a cambios de instrumentos, es necesario recalcular con cierta frecuencia los coeficientes del modelo para optimizar la precisión de apuntado de telescopio. Pero esto requiere de un método de fácil implementación que permita el cálculo de los coeficientes por parte del personal técnico del Observatorio, a fin de contar siempre con el mejor modelo de apuntado. En este trabajo se describe un modelo de ajuste trigonométrico y, se propone un programa con el que se ha conseguido encontrar los coeficientes para la corrección de los errores de alineación al polo, falta de perpendicularidad de los ejes, flexiones mecánicas y falta de alineación del eje óptico. El programa está escrito en Fortran y hace uso de las nuevas características del Fortran 95, lo que permite hacer una gran cantidad de cálculos de manera muy simple. Existe también una versión para IDL y en breve estará lista una versión del programa en Fortran77. El caso que se presenta como ejemplo usa datos que

corresponden al telescopio de 1 m de diámetro del OAN en Tonantzintla, pero el método puede aplicarse directamente a cualquier tipo de telescopios con montura ecuatorial. Aplicando el modelo de correcciones a una base de datos de errores de apuntado, los errores rms se han reducido en Declinación de 49.8” a 8.9” y en Ángulo Horario de 99.7” a 30.1”.

2. Entendiendo el problema Un objeto en el cielo se identifica por sus coordenadas: ascensión recta y declinación (AR/Dec o α/δ). Las declinaciones se miden hacia el norte y hacia el sur del ecuador celeste, mientras que la ascensión recta se mide hacia el este del equinoccio (la intersección entre el ecuador celeste y la eclíptica, siendo esta última el plano de la órbita de la Tierra alrededor del Sol). Pero, conocer AR y Dec no es el final de la historia. Hay una serie de correcciones y transformaciones astronómicas para hacer antes de estar listos para apuntar el telescopio. Luego, es necesario traducir estas coordenadas a otras coordenadas que sirvan para mover el telescopio. Pero en un telescopio ecuatorial (ver figuras 1, 2 y 3) la situación es, en teoría, simple, ya que si el eje polar del telescopio es paralelo al eje polar de la Tierra, la declinación del telescopio coincide con la declinación celeste; el otro eje del telescopio, denominado ángulo horario (AH) y medido hacia el este y el oeste del meridiano, se relaciona con la ascensión recta mediante la ecuación AR = TS – AH, donde TS es el tiempo sideral local (el tiempo sideral local corre un poco mas rápido que el tiempo regular, de tal manera que se adelanta casi 4 minutos diarios).

Figura 1. Dos tipos de monturas ecuatoriales simétricas.

Figura 2. Dos tipos de monturas ecuatoriales asiméticas.

Pero el problema no ha terminado. Aunque ya tenemos las coordenadas del telescopio y estamos listos para moverlo, hemos olvidado que vivimos en un planeta que gira y se balancea, que vemos a través de una atmósfera y que el telescopio no es una máquina perfecta. Todos estos factores han de tomarse en cuenta antes de apuntar el telescopio. La secuencia de cálculos y ajustes que debemos hacer se muestra en la figura 4.

Figura 3. Telescopio de 2.1m de SPM. Su montura es ecuatorial de yugo (simétrica).

Los tres pasos principales son: pasar de coordenadas medias a aparentes, de aparentes a observadas y de observadas a instrumentales. El primer paso corrige el hecho de que el eje de la Tierra y, por tanto, el ecuador celeste están en constante movimiento debido a la precesión y a

la nutación, así como el hecho de que nuestro movimiento alrededor del Sol hace que las estrellas parezcan desplazarse debido a la aberración anual. La transformación de aparentes a observadas se debe a la rotación de la Tierra, a la ubicación geográfica del telescopio y a la refracción atmosférica. El último paso consiste en corregir las imperfecciones instrumentales del telescopio.

Figura 4. Diagrama que muestra la secuencia de las correcciones.

Las primeras trasformaciones se realizan con ecuaciones que se aplican de igual manera en cualquier telescopio y no dependen de la construcción del mismo. Sin embargo, la última es una transformación única para cada telescopio, pues sólo depende de sus propias características.

La precesión del polo terrestre alrededor del polo de la eclíptica, junto con una inclinación gradual de la eclíptica, es la componente de largo plazo de un movimiento complejo causado por los efectos gravitacionales de los miembros del Sistema Solar sobre la Tierra que gira. El periodo de este movimiento es de 26000 años y produce cambios de hasta 50 segundos de arco por año en las coordenadas de las estrellas. El “cabeceo” residual del eje terrestre, que tiene un período de 18.6 años, es la nutación y afecta el apuntado del telescopio en el orden de los 10 segundos de arco. La aberración anual es culpable de hasta 20 segundos de error en el apuntado de un telescopio. Luego, dado que el telescopio está en la Tierra y dentro de la atmósfera, es necesario hacer las otras correcciones para lo cual es necesario conocer la latitud, la longitud y altura sobre el nivel del mar del sitio donde está el observatorio, así como la presión y la temperatura ambiente y el tiempo exacto. Los errores que pueden cometerse por el efecto de la refracción son hasta del orden de 1 minuto de arco para una elevación de 45º al nivel del mar. Todas estas correcciones proporcionan las coordenadas observadas, que, si todo fuera perfecto, nos permitiría ver la estrella perfectamente centrada en el telescopio. Pero el telescopio es real y, por tanto, sus codificadores pueden tener un ligero desplazamiento en sus lecturas, puede haber un ligero juego en los engranes (backlash), sus ejes pueden no ser perfectamente perpendiculares y el eje polar del telescopio puede no estar apuntando exactamente al polo. Corregir físicamente estas imperfecciones puede ser una tarea larga y costosa, por lo que en ocasiones es mejor aceptarlas y tenerlas en cuenta al calcular las coordenadas del objeto.

3. El modelo de apuntado En una montura ecuatorial, estas imperfecciones causan un efecto que puede modelarse mediante seis términos puramente geométricos, complementados con tres debidos a flexiones (Wallace y Tritton,1975; Trueblood y Genet, 1998; Wallace, 2001). En la tabla 1 se muestran estos términos. Los términos HI y DI son simplemente errores de inicialización en las lecturas de las coordenadas instrumentales y pueden calcularse con un apuntado inicial del telescopio. El término CO describe qué tan precisa es la alineación entre el eje mecánico del telescopio y el eje óptico y si el instrumento de detección está bien alineado (la componente que afecta a declinación no depende del ángulo, por lo que queda integrada en DI). El término NP depende de qué tan lejos del ángulo recto está el eje de declinación con respecto al eje polar. Término

Descripción

Δh

Δδ

HI

Error constante de posición

HI

DI

Error constante de posición

CO

Error de colimación

CO*secδ

NP

δ y el eje polar no perpendiculares

NP*tanδ

PV

Efecto en elevación de la no coincidencia con el polo

PV*senh*tanδ

PV*cosh

PH

Efecto perpendicular a PV

-PH*cosh*tanδ

PH*senh

FT

Flexión del tubo

FT*cosφ*senh*senδ

FT*(cosφ*cosh*senδsenφ*cosδ)

FY

Flexión del yugo

FD

Flexión del eje de declinación

DI

FY*cosh -FD*(cosφ*cosh+ senφ*tanδ)

Tabla 1. Términos que afectan en una montura ecuatorial. h y δ son el ángulo horario y la declinación y φ es la latitud del observatorio.

Los términos PH y PV describen qué tan apartado está el eje polar del eje verdadero en la dirección horizontal en el caso de PH y en la vertical para PV. FT describe cuánto cae el telescopio al acercarse más a una posición horizontal. FY es la flexión del yugo (o tenedor) y FD es un término de flexiones para el eje de declinación, que sólo se aplica en el caso de los telescopios como los de montura alemana, en la que el eje de declinación no está soportado en ambos extremos. Además de estos términos, generalmente existen términos menores que presentan un comportamiento periódico, causados por defectos de maquinado y de alineación de los diferentes ejes y acoplamientos, que generalmente son de magnitud menor. Estos no se incluyen en el modelo aquí descrito. En la Tabla 1, cada uno de los coeficientes es en general un ángulo pequeño, usualmente del orden de segundos de arco. Las fórmulas para corregir cada uno de los ejes están dadas por la suma de los elementos en sus respectivas columnas. Así, despreciando el término FD, las fórmulas de corrección para el modelo de apuntado son: Δh = HI + CO*secδ + NP*tanδ - PH*cosh*tanδ + PV*senh*tanδ + FT*cosφ*senh*senδ Δδ = DI + PH*senh + (PV+FY)*cosh + FT*(cosφ*cosh*senδ-senφ*cosδ)

Aplicando este modelo puede mejorarse notablemente el apuntado del telescopio.

4. ¿Cómo encontrar los coeficientes? Encontrar los coeficientes del modelo no es una tarea fácil (excepto por HI y DI). Además, los coeficientes varían ligeramente con el tiempo y un modelo no puede durar toda la vida: la fatiga de los elementos del telescopio y el esquema instrumental provocarán ligeras variaciones en los coeficientes de flexión, el desplazamiento del eje óptico durante los procesos de mantenimiento provocarán una variación en el coeficiente CO y los asentamientos del suelo provocarán variaciones en la alineación del eje polar del telescopio. Esto hace que deban calcularse los coeficientes con cierta periodicidad. Por esta razón, se ha creado el programa encuentra_modelo. Este programa ha sido escrito en Fortran, usando las nuevas implementaciones incorporadas en el Fortran 95 y se debe usar en ambiente Linux. También es importante que se encuentre instalada una versión del Gnuplot. Existe también una versión para IDL. Para usarlo, es necesario apuntar el telescopio a diferentes estrellas brillantes en el cielo, cuyas coordenadas observadas hayan sido previamente calculadas. Es importante incluir diferentes posiciones en el cielo, anotando la corrección manual que haya tenido que hacerse, tanto en ángulo horario como en declinación. Los datos deben tabularse en un archivo de texto de la siguiente manera: AH

dAH

Dec

dDec

donde AH y Dec son el ángulo horario (en horas con respecto al meridiano) y la declinación (en grados) de las coordenadas calculadas del objeto y dAH y dDec son las correcciones que tuvieron que hacerse para centrar el objeto en el telescopio, dadas en segundos de arco. Por ejemplo: -2.473

-33

-16.686

-50.2

-1.799

42

22.515

16.0

(AH)

(dAH)

(Dec)

(dDec)

Es importante incluir el mayor número de puntos en la muestra para garantizar un mejor ajuste del modelo. Se recomiendan unos 40 puntos, para tener del orden de 10 en cada cuadrante del cielo.

5. Descripción del programa A este texto se anexa en el Apéndice A un listado del programa en Fortran y en el Apéndice B el correspondiente listado de la versión para IDL. El programa escrito en Fortran funciona de la siguiente manera: a) Busca el archivo con los datos. El programa solicita al usuario que proporcione el nombre del archivo (el que se ha usado para hacer las pruebas se llama ‘epos.txt’). b) Solicita al usuario la latitud del lugar en grados. c) Cuenta cuántas líneas no vacías tiene el archivo, a fin de alojar en memoria un array suficientemente grande para contener los datos. Los elementos de este array son de tipo derivado y contienen cuatro elementos de tipo real: AH, Dec, dAH y dDec. d) Grafica los errores dDec contra dAH para que el usuario tenga una idea de la magnitud de los errores. e) Calcula el término cosφ*cosh*senδ-senφ*cosδ y ajusta la recta tomando como eje vertical a dDec. La pendiente encontrada corresponde al termino FT. Presenta al usuario los datos (pendiente, ordenada al origen y correlación del ajuste), quien puede aceptar o no la corrección. f) Una vez corregidos los errores por este término, se calcula el término cosh y repite el proceso anterior. Esto da como resultado la suma de los términos PV y FY. Nuevamente, pregunta al usuario si acepta la corrección. g) Corrige los errores, si así lo ha decidido el usuario, y calcula el término senh. Al repetir el proceso se tiene el término PH. De esta manera se han encontrado ya los coeficientes que corrigen declinación. h) Grafica, llamando al “gnuplot”, los datos mostrando cómo estaban distribuidos los datos durante los ajustes.

i) Incorpora los términos calculados (FT y PH) y corrige los errores en ángulo horario. Revisa si existe algún posible backlash en el mecanismo y lo corrige. j) Calcula el término senh*tanδ y ajusta de nuevo una recta, tomando ahora como eje vertical a dAH. La pendiente del ajuste es ahora PV. Esto permite discriminar ya también el término FY, pues sólo se conocía la suma PV+FY. k) Si el usuario acepta la corrección, se corrigen los errores y se calcula ahora el término tanδ. Esto permite encontrar el término NP. l) Finalmente, calcula el término secδ y encuentra el término CO. m) Grafica ahora los datos de ángulo horario mostrando cómo estaban distribuidos durante los ajustes. También muestra los errores corregidos graficando nuevamente dDec contra dAH y presenta los valores de los coeficientes encontrados. El programa escrito para IDL funciona de manera parecida, con algunas pequeñas variantes: a) El programa recibe como argumentos lel nombre del archive con los datos, la latitud del observatorio y un arreglo de “unos” y “ceros” que indican qué correcciones se harán. Por ejemplo: ENCUENTRA_MODELO, archivo="epos.txt", latitud=19.5, corrige=[1,1,1,0,1,0,1]

indica que se harán las correcciones correspondientes a PH, PV, FT, NP y BL y que no se aplicarán las correcciones por el término combinado PV_FY ni la de colimación óptica. El orden es: PH, PV, FT, PV_FY, NP, CO, BL La ausencia de este arreglo en el argumento hará que se calculen todas las correcciones. De hecho, se recomienda habilitar todas las

correcciones en una primera ejecución y determinar de allí los términos que pueden o deben omitirse. La latitud deberá darse en grados. b) Busca el archivo con los datos. c) Grafica los errores dDec contra dAH para que el usuario tenga una idea de la magnitud de los errores. d) Calcula el término cosφ*cosh*senδ-senφ*cosδ y ajusta la recta tomando como eje vertical a dDec. La pendiente encontrada corresponde al termino FT. e) Una vez corregidos los errores por este término, se calcula el término cosh y repite el proceso anterior. Esto da como resultado la suma de los términos PV y FY. f) Corrige los errores, si así lo ha decidido el usuario con los valores del argumeno, y calcula el término senh. Al repetir el proceso se tiene el término PH. De esta manera se han encontrado ya los coeficientes que corrigen declinación. g) Grafica los datos usando las propias utilerías de IDL, mostrando cómo estaban distribuidos durante los ajustes. h) Incorpora los términos calculados (FT y PH) y corrige los errores en ángulo horario. Revisa si existe algún posible backlash en el mecanismo y lo corrige. i) Calcula el término senh*tanδ y ajusta de nuevo una recta, tomando ahora como eje vertical a dAH. La pendiente del ajuste es ahora PV. Esto permite discriminar ya también el término FY, pues sólo se conocía la suma PV+FY. j) Se corrigen los errores y se calcula ahora el término tanδ. Esto permite encontrar el término NP. k) Finalmente, calcula el término secδ y encuentra el término CO.

l) Grafica ahora los datos de ángulo horario mostrando cómo estaban distribuidos durante los ajustes. También muestra los errores corregidos graficando nuevamente dDec contra dAH y presenta los valores de los coeficientes encontrados. Los términos HI y DI encontrados durante este ajuste no deben ser definitivos, ya que dependen del proceso de inicialización al empezar una sesión de observación. Estos valores deben recalcularse cada vez que se inicia una sesión de observación, apuntando una estrella después de haber aplicado el modelo calculado. El error en este primer apuntado corresponde directamente a los términos HI y DI. Si se desea, pueden apuntarse varias estrellas y promediar los errores encontrados. Pero este hecho puede obviarse en los telescopios del OAN en SPM, pues el sistema de control de los telescopios calcula estos términos cuando se corrigen las coordenadas después del primer apuntado. Los coeficientes encontrados pueden sustituirse ahora en las fórmulas de la sección 3 para hacer las correcciones a las coordenadas. Para obtener mejores resultados de este modelo trigonométrico, en el caso de los telescopios del OAN en San Pedro Mártir, se sugiere inhibir el modelo polinomial que actualmente corrige los errores mecánicos de apuntado, calcular los coeficientes del modelo trigonométrico y, después de aplicarlo, recalcular el modelo polinomial, pues con este último pueden modelarse los demás errores residuales (errores de los acoplamientos de los codificadores, errores de los engranes, etc.) que prácticamente no cambian con el tiempo.

6. Resultados Para probar el algoritmo, hemos aplicado el programa a una base de datos con errores de apuntado del telescopio del Observatorio Astronómico Nacional en Tonantzintla y hemos encontrado los siguientes coeficientes para el modelo (en este caso se usó descartaron las correcciones para CO y para el término combinado PV_FY:

Término

Valor en segundos de arco

Offset en AH

HI

9.328

Offset en Dec

DI

9.550

Error de colimación

CO

0.000

Falta de perpendicularidad

NP

67.520

Error horizontal en el polo

PH

-17.693

Error horizontal en el polo

PV

-117.582

Flexión del tubo

FT

137.124

Flexión del yugo

FY

117.582

Backlash en AH

BL

171.394

Estos valores permiten disminuir el valor rms de los errores desde 99.7” y 49.8” para AH y Dec, respectivamente, hasta 30.1” y 8.9”. En las figuras 5, 6, 7 y 8 se reproducen las gráficas entregadas por el programa. El código de ambas versiones del programa presentadas aquí se encuentra en la dirección: http://www.astrosen.unam.mx/~leonel/progs/modelo/

Figura 4. Gráfica con los errores de Dec y de AH, antes de ser corregidos. Los aparentes 2 grupos de datos se deben al error por juego mecánico (backlash).

Figura 5. Gráfica de los errores de la figura 5 después de ser corregidos. Nótese la escala de los ejes.

Figura 6. Gráficas de los errores de Declinación durante el proceso de corrección.

Figura 7. Gráfica de los errores en Ángulo Horario durante el proceso de corrección. Nótese el escalón en la primera gráfica al pasar por cero en AH debido al backlash.

7. Referencias 1. Wallace, P. T, y Tritton, K. P., 1979, MNRAS, 189, 115. 2. Trueblood, M. y Genet, R. M., Telescope Control, Willman-Bell, Inc., Richmond, Virginia, USA, 1998. 3. Wallace, P., “Telescope Pointing”, http://www.tpsoft.demon.co.uk/pointing.htm (2001).

Apéndice A. Listado del programa en Fortran95 ! Listado del programa !******************************************************************************************! ! Programa 'encuentra_modelo' elaborado por Leonel Gutierrez ! ! Este programa lee un archivo con una serie de datos que muestran los errores de ! apuntado de un telescopio con montura ecuatorial y calcula los coeficientes de ! un modelo trigonometrico para corregir dicho apuntado. ! !******************************************************************************************! ! ! Lo que sigue es un modulo con los elementos necesarios para el programa ! Contiene las siguientes funciones y subrutinas ! ! 1) ! Funcion para contar los renglones del archivo ! function cuantas_lineas(archivo) result(cuantas) ! ! 2) ! Subrutina para traer todos los elementos a sus lugares en el arreglo ! subroutine trae_lineas(archivo, arreglo) ! ! 3) ! Subrutina para guardar en un archivo provisional los datos a graficar ! subroutine guarda_lineas(archivo, matriz) ! ! 4) ! Subrutina para ajustar por minimos cuadrados una recta a una distribucion. Calcula ! ! la pendiente, la ordenada al origen y la correlacion ! subroutine ajusta_recta(x, y, m, b, corr) ! ! 5) ! Funcion que calcula el valor rms de un vector. ! function calcula_rms(x) ! ! 6) ! Subrutina que prepara el archivo que se le va a entregar al 'gnuplot' ! subroutine prepara_archivo_de_grafica() ! ! 7) ! Subrutina que agrega las diferentes graficas al archivo que usara el 'gnuplot'. ! ! Podemos poner hasta 4 graficas en una ventana. (El 'gnuplot' permite muchas ! ! mas, pero en este programa nos hemos limitado a 4.) ! subroutine agrega_grafica(x, y, sx, sy, posx, posy, etiqueta, archivo_prov) !

! 8) ! ! ! 9) ! ! !10) ! ! !11) !

! Subrutina que inicia una grafica llamando al 'gnuplot' subroutine inicia_grafica() ! Subrutina que grafica los errores de Declinacion contra los de AH subroutine grafica_errores(etiqueta) ! Subrutina que calcula los coeficientes del modelo usando los errores de Declinacion subroutine encuentra_coeficientes_dec() ! Subrutina que calcula los demas coeficientes del modelo usando los errores de AH subroutine encuentra_coeficientes_ah()

module utilerias ! variables y constantes para usarse en el modulo integer, parameter character(len = 256) integer

:: dp = selected_real_kind(8,30) ! La clase de doble precision :: linea ! define una cadena para ir leyendo el archivo :: ioss ! variable para almacenar los errores del READ

! Voy a poner los reales de doble precision para no tener problemas de redondeo ! en las funciones trigonometricas real(kind=dp), parameter :: pi = acos(-1.) real(kind=dp), parameter :: gra_rad = pi/180. ! Factor multiplicativo para convertir ! grados en radianes real(kind=dp), parameter :: horas_rad = gra_rad*15.0 ! Factor multiplicativo para convertir ! horas de cielo en radianes integer :: d1, d2, d3,d4 ! dummies real(kind=dp), dimension(4) :: correlacion ! para guardar las correlaciones real(kind=dp), dimension(4) :: m ! las pendientes real(kind=dp), dimension(4) :: b ! y las ordenadas character(len=10) :: resp ! para los dialogos real(kind=dp) real(kind=dp) real(kind=dp) real(kind=dp) real(kind=dp) real(kind=dp) real(kind=dp) real(kind=dp)

:: :: :: :: :: :: :: ::

PH PV FT FY NP CO PV_FY DAF

= = = = = = = =

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

! ! ! ! ! ! ! ! !

Error horizontal en el polo Error vertical en el polo Flexiones del tubo Flexiones del yugo Falta de perpendicularidad Colimacion Flexiones del yugo + error vertical en polo Flexiones del eje de declinacion; No se aplica en este caso

real(kind=dp) real(kind=dp) real(kind=dp)

:: BL :: HI :: DI

= 0.0 = 0.0 = 0.0

real(kind=dp) real(kind=dp) real(kind=dp)

:: phi :: senphi :: cosphi

! Error de backlash ! Offset inicial en AH ! Offset inicial en Dec

= 19.5*gra_rad

! La latitud del observatorio ! Esto es el seno de la latitud ! Esto es el coseno de la latitud

! variables derivadas type coordenadas real(kind=dp) real(kind=dp) real(kind=dp) real(kind=dp) end type coordenadas

:: AH :: Dec :: dAH :: dDec

! Angulo horario del telescopio ! Declinacion ! Error en AH ! Error en Dec

! variables globales ! En este arreglo pongo todos los datos del archivo type(coordenadas), dimension(:), allocatable :: arreglo ! ! ! ! ! !

Esta es una matriz de n x 9 donde pondremos los ejes a graficar; las columnas quedaran asignadas como sigue: 1) para el eje x 2) para los errores antes de corregirse; estos se graficaran en el eje y En las demas columnas se pondran los mismos errores conforme se les vaya aplicando las correcciones

real(kind=dp), dimension(:,:), allocatable

:: matriz

integer

:: ii

! un indice

contains !***************************************************! ! Funcion para contar los renglones no vacios del archivo. Esto permite reservar el numero de lugares ! necesario para los arreglos usados en el programa. Basicamente son 'arreglo' y 'matriz', asi como ! algunos vectores auxiliares dentro de las subrutinas ! function cuantas_lineas(archivo) result(cuantas) character(*), intent(in) :: archivo integer :: cuantas

logical :: existe ! variable para verificar que el archivo existe inquire(FILE = archivo, EXIST = existe) cuantas = 0 if (existe) then open (unit=1,file=archivo,status='old') i = 0 do read(1,'(a)',IOSTAT=ioss) linea if(IOSS == 0) then if(len_trim(linea) > 0) i = i + 1 else exit endif enddo cuantas = i close(1) else print*,"No existe el archivo" end if end function cuantas_lineas !***************************************************! ! Subrutina que lee todos los renglones del archivo de entrada y pone todos los elementos en ! sus lugares correspondientes en el arreglo. ! Subroutine trae_lineas(archivo, arreglo) character(*), intent(in) :: archivo type(coordenadas), dimension(*) :: arreglo logical

:: existe

! variable para verificar que el archivo existe

inquire(FILE = archivo, EXIST = existe) if (existe) then open (unit=1,file=archivo,status='old') i = 1 do read(1,'(a)',IOSTAT=ioss) linea if(IOSS == 0) then if(len_trim(linea) > 0) then

read(linea,*), arreglo(i)%AH, arreglo(i)%dAH, & arreglo(i)%Dec, arreglo(i)%dDec i = i + 1 endif else exit endif end do close(1) else print*,"No existe el archivo" end if end Subroutine trae_lineas !***************************************************! ! Esta subrutina guarda la matriz en un archivo provisional para que sea de alli donde tome el 'gnuplot' los ! datos a graficar. El archivo provisional no se borrar al terminar el programa, a fin de que el usuario ! pueda tener los datos como referencia. El archivo para graficar los errores de Declinacion en cada etapa ! de la correccion, quedan guardados en el archivo 'prov.txt'; los de AH quedan en 'prov_ah.txt'; los errores ! ya corregidos por el modelo quedan en 'prov_i.txt'. Al reiniciarse el programa, los archivos se sobreescriben. ! Subroutine guarda_lineas(archivo, matriz) character(*), intent(in) :: archivo real(kind=dp), dimension(:,:), intent(in) :: matriz integer n = size(matriz,1)

:: n

open (unit=1,file=archivo,status='replace') i = 1 do write(1,'(9(f10.5,3X))',IOSTAT=ioss) matriz(i,1), matriz(i,2), matriz(i,3), matriz(i,4), & matriz(i,5), matriz(i,6), matriz(i,7), matriz(i,8), matriz(i,9) if(IOSS == 0) then i = i + 1 if (i == n) exit else exit endif end do close(1) end Subroutine guarda_lineas

!***************************************************! ! Esta subrutina ajusta una recta usando el metodo de minimos cuadrados ! subroutine ajusta_recta(x, y, m, b, corr) real(kind = dp), dimension(:), intent(in) :: x real(kind = dp), dimension(:), intent(in) :: y real(kind = dp), intent(out) :: m real(kind = dp), intent(out) :: b real(kind = dp), intent(out) :: corr real, dimension(:), allocatable integer :: n

:: va, vb, vc ! vectores auxiliares ! tamanio de x o y

! variables para calcular la estadistica real(kind = dp) :: sumax, sumay, sumaxy, sumax2, & cociente, varx, vary, covxy, Xm, Ym ! Procedemos con los minimos cuadrados canonicos n

= size(x)

allocate(va(n)) allocate(vb(n)) allocate(vc(n)) sumay sumax va sumax2 va sumaxy cociente b m

! Aqui calculamos el tamanio del arreglo ! Reserva n lugares para va, vb y vc

= = = = = = =

sum(y) ! Suma de todos los elementos en y sum(x) ! Suma de todos los elementos en x x*x sum(va) ! Suma de todas las x^2 x*y sum(va) ! Suma de todas los elementos cruzados n*sumax2 - sumax*sumax ! La ordenada al origen = (sumay*sumax2-sumax*sumaxy)/cociente ! La pendiente = (n*sumaxy-sumax*sumay)/cociente

! Vamos a sacar ahora la correlacion Xm = sumax/n ! La media de x Ym = sumay/n ! La media de y va vb

= x - Xm = y - Ym

! Calcula los errores en x ! Calcula los errores en y

vc

= vb*vb

sumay2 vc

= sum(vc) = va*va

sumax2 vc sumaxy

= sum(vc) = va*vb = sum(vc)

varx vary covxy corr

= = = =

! Calcula el vector de los cuadrados de las ! diferencias en y ! Aqui la suma de todo ese vector ! Calcula el vector de los cuadrados de las ! diferencias en x ! Aqui la suma de todo ese vector ! Calcula el vector de los productos cruzados ! Aqui la suma de todo ese vector

sumax2/n ! La varianza en x sumay2/n ! La varianza en y sumaxy/n ! La covarianza covxy/(sqrt(varx*vary)) ! Finalmente ya tenemos el coeficiente de correlacion ! Libera la memoria

deallocate(va) deallocate(vb) deallocate(vc) end subroutine ajusta_recta

!***************************************************! ! Esta funcion calcula el valor rms de un vector. En este caso, se trata de un ! vector de errores ! function calcula_rms(x) result(rms) real(kind = dp), dimension(:), intent(in) :: x real(kind = dp) :: rms real, dimension(:), allocatable integer :: n

:: va, vb, vc ! vectores auxiliares ! tamanio de x o y

! variables para calcular la estadistica real(kind = dp) :: sumax, sumax2, varx, Xm, Ym n

= size(x)

! Aqui calculamos el tamanio del arreglo

allocate(va(n)) ! Vamos a calcular el RMS sumax = sum(x) Xm = sumax/n va = x - Xm va = va*va

! Reserva n lugares para va

! ! ! ! !

Suma de todos los elementos en x La media de x Calcula los errores en x Calcula el vector de los cuadrados de las diferencias en

sumax2

= sum(va)

varx rms

= sumax2/n = sqrt(varx)

deallocate(va) end function calcula_rms

! Aqui la suma de todo ese vector ! La varianza en x ! ! Finalmente ya tenemos el coeficiente de correlacion ! Libera la memoria

!***************************************************! ! Esta subrutina prepara el archivo que luego se le entregara al 'gnuplot' para hacer las ! graficas ! subroutine prepara_archivo_de_grafica() open (unit=3,file="pp.prg",status='replace') write(3,'(A)') "set origin 0,0" write(3,'(A)') "set size 1,1" write(3,'(A)') "set multiplot" close(3) end subroutine prepara_archivo_de_grafica !***************************************************! ! Esta subrutina ordena al sistema que arranque el 'gnuplot' usando el archivo 'pp.prg' ! como archivo de lotes para hacer las graficas ! subroutine inicia_grafica() character(len = 25) :: salida salida = "gnuplot -persist pp.prg" call system(salida) end subroutine inicia_grafica !***************************************************! ! Esta subrutina agrega al archivo 'pp.prg' los mandos necesarios para insertar una grafica en ! la ventana del 'gnuplot' ! subroutine agrega_grafica(x, y, sx, sy, posx, posy, etiqueta, archivo_prov) ! 'x' es el numero de columna que se usara en el eje x; lo mismo con 'y' ! 'sx' es la etiqueta que ira en el eje x; lo mismo con y integer, intent(in) :: x, y character(len = *), intent(in) :: sx, sy ! etiquetas para la grafica integer, intent(in) :: posx, posy ! pueden valer 0 y 1, para ! indicar en que cuadrante se van a poner

! 'archivo_prov' es donde se han puesto los datos a graficar character(len = *), intent(in) :: etiqueta character(len = *), intent(in) :: archivo_prov integer :: mx, my, mmx, mmy ! Calcula los valores maximos y minimos mx = floor(minval(matriz(:,x))) my = floor(minval(matriz(:,y))) mmx = ceiling(maxval(matriz(:,x))) mmy = ceiling(maxval(matriz(:,y)))

! para poner los max y min

open (unit=3,file="pp.prg",status='old', position='append') write(3,'(A12,I4,A1,I4,A1)') "set xrange [", mx, ":", mmx, "]" write(3,'(A12,I4,A1,I4,A1)') "set yrange [", my, ":", mmy, "]" write(3,'(A)') "set size 0.5,0.5" write(3,'(A,F3.2,A1,F3.2)') "set origin ", 0.5*posx,",", 0.5*posy write(3,'(A12,A)') "set xlabel ", sx write(3,'(A12,A)') "set ylabel ", sy write(3,'(A6,A,A7,I2,A1,I2,A3,A)') "plot '",archivo_prov, "' using ", x, ":", y, " t ",etiqueta close(3) end Subroutine agrega_grafica !***************************************************! ! Esta subrutina grafica los errores de Declinacion contra los errores de Angulo horario para que el ! usuario tenga una idea de la magnitud de los errores. Se llama al inicio del programa y al final ! para poder comparar los errores iniciales contra los finales ! subroutine grafica_errores(etiqueta, rmsdec, rmsah) character(len = *), intent(in) :: etiqueta real, intent(out) :: rmsdec, rmsah integer

:: mx, my, mmx, mmy

matriz(:,1) = arreglo%dDec matriz(:,2) = arreglo%dAH ! Calcula los valores maximos y minimos mx = floor(minval(matriz(:,1))) my = floor(minval(matriz(:,2))) mmx = ceiling(maxval(matriz(:,1))) mmy = ceiling(maxval(matriz(:,2)))

! para poner los max y min

rmsdec = calcula_rms(matriz(:,1)) rmsah = calcula_rms(matriz(:,2)) call guarda_lineas("prov_i.txt", matriz) open (unit=3,file="pp.prg",status='replace') write(3,'(A)') "set origin 0,0" write(3,'(A)') "set size 1,1" write(3,'(A12,I4,A1,I4,A1)') "set xrange [", mx, ":", mmx, "]" write(3,'(A12,I4,A1,I4,A1)') "set yrange [", my, ":", mmy, "]" write(3,'(A,F5.1,A)') "set xlabel 'dDec: ", rmsdec, " rms'" write(3,'(A,F5.1,A)') "set ylabel 'dAH: ", rmsah, " rms'" write(3,'(A,A)') "plot 'prov_i.txt' using 1:2 t ", etiqueta close(3) call inicia_grafica() end subroutine grafica_errores !***************************************************! ! Esta subrutina se encarga de hacer los calculos necesarios para encontrar los coeficientes del ! modelo de correccion correspondientes a Declinacion. Solicita al usuario la confirmacion de ! las correcciones. El usuario debe decidir si se hacen las correcciones o no, en funcion de la ! informacion presentada, en particular la correlacion del ajuste ! subroutine encuentra_coeficientes_dec() ! Pone en la primera columna los errores en declinacion antes de ninguna correccion integer :: tam, kk matriz(:,1) = arreglo%dDec ! En la segunda columna pone el termino de flexiones del tubo; en los calculos, va ! multiplicando a Dec, que esta dada en grados, por el termino 'gra_rad' para ponerla ! en radianes y a AH, que esta en horas, por 'horas_rad' matriz(:,2) = -senphi*cos(arreglo%Dec*gra_rad) + & cosphi*sin(arreglo%Dec*gra_rad)*cos(arreglo%AH*horas_rad); !!!!! Primer ajuste (para el termino FT) !!!!! ! Calcula la mejor recta que se ajusta a esta distribucion call ajusta_recta(matriz(:,2), matriz(:,1), m(1), b(1), correlacion(1)) ! Imprime los valores de m, b y la correlacion para que el usuario decida si se

! aplica la correccion write(*,'(A,3F14.6)') " Flexiones del tubo: m, b y corr=", m(1), b(1), correlacion(1) ! Pregunta si aplica la correccion write(*,'(/a)', advance = "no") "Aplico esta correccion (s,n)? " read(*,*) resp if(resp == "s") then ! Corrige los errores usando la recta encontrada arreglo%dDec = arreglo%dDec-( m(1)*matriz(:,2) + b(1) ); FT = m(1) DI = DI + b(1) endif ! En la tercera columna de la matriz pone los errores con esta primera correccion; si no hubo ! correccion se dejan como estaban pero se ponen en esa columna para poder graficarlos matriz(:,3) = arreglo%dDec ! En la cuarta columna pone el termino de flexiones del yugo + alineacion vertical al polo matriz(:,4) = cos(arreglo%AH*horas_rad); !!!!! Segundo ajuste (para el termino PV_FY) !!!!! ! Calcula la mejor recta que se ajusta a esta distribucion call ajusta_recta(matriz(:,4), matriz(:,3), m(1), b(1), correlacion(1)) write(*,'(A,3F14.6)') " Flexiones del yugo y no alin. vertical: m, b y corr=", m(1), b(1), correlacion(1) ! De nuevo, dejamos que el usuario decida write(*,'(/a)', advance = "no") "Aplico esta correccion (s,n)? " read(*,*) resp if(resp == "s") then ! Corrige los errores usando la recta encontrada arreglo%dDec = arreglo%dDec-( m(1)*matriz(:,4) + b(1) ); PV_FY = m(1) DI = DI + b(1) endif ! En la columna 5 de la matriz pone los errores para ser graficados despues matriz(:,5) = arreglo%dDec ! Y en la columna 6 pone el termino de alineacion horizontal al polo matriz(:,6) = sin(arreglo%AH*horas_rad)

!!!!! Tercer ajuste (termino PH) !!!!! ! Calcula la mejor recta que se ajusta a esta distribucion call ajusta_recta(matriz(:,6), matriz(:,5), m(1), b(1), correlacion(1)) write(*,'(A,3F14.6)') " No alin. horiz.: m, b y corr=", m(1), b(1), correlacion(1) ! Igual. Dejamos que el usuario decida write(*,'(/a)', advance = "no") "Aplico esta correccion (s,n)? " read(*,*) resp if(resp == "s") then ! Corrige los errores usando la recta encontrada arreglo%dDec = arreglo%dDec-( m(1)*matriz(:,6) + b(1) ); PH = m(1) DI = DI + b(1) endif ! En la columna 7 de la matriz pone los errores para ser graficados despues matriz(:,7) = arreglo%dDec ! Aqui los datos estan acomodados como sigue: ! Col 1 | 2 | 3 | 4 | 5 | 6 | 7 ! ! --------------------------------------------------------------! ! dDec | FT | dDec_FT | PV_FY | dDec_PV_FY | PH | dDec_PH ! call guarda_lineas("prov.txt", matriz) call call call call call call

prepara_archivo_de_grafica() agrega_grafica(2,1,'"FT"', '"dDec"', 0, agrega_grafica(4,3,'"PV+FY"', '"dDec"', agrega_grafica(6,5,'"PH"', '"dDec"', 0, agrega_grafica(2,7,'"FT"', '"dDec"', 1, inicia_grafica()

1,'"Flexiones del tubo"', "prov.txt") 1, 1,'"Flex. yugo + No alin. vert."', "prov.txt") 0,'"No alin. horiz."', "prov.txt") 0,'"Errores corregidos vs. FT"', "prov.txt")

tam = size(matriz(:,7)) do kk = 1, tam if( matriz(kk,7) < -30) then write(*,*) "El valor es: ",matriz(kk,7), "con kk = ", k endif enddo end subroutine encuentra_coeficientes_dec

!***************************************************! ! Como en el caso de Declinacion, esta subrutina se encarga de hacer los calculos necesarios para ! encontrar los coeficientes del modelo de correccion correspondientes al Angulo Horario. Tambien ! solicita al usuario la confirmacion de las correcciones, para que sea el usuario el que decida ! si se hacen las correcciones o no. ! subroutine encuentra_coeficientes_ah() real :: media_mas, media_menos integer :: cuantos_mas integer :: cuantos_menos media = 0.0 cuantos_mas = 0 cuantos_menos = 0 media_mas media_menos cuantos_mas cuantos_menos media_mas media_menos

= = = = = =

sum(arreglo%dAH, mask = arreglo%AH > 0) sum(arreglo%dAH, mask = arreglo%AH 0) count(arreglo%AH 0) arreglo%dAH = arreglo%dAH - (media_mas - media_menos) BL = (media_mas - media_menos) endif ! Ahora corrige con los coeficientes que ya conoce arreglo%dAH = arreglo%dAH - FT*cosphi*sin(arreglo%AH*horas_rad)*sin(arreglo%Dec*gra_rad) arreglo%dAH = arreglo%dAH + PH*cos(arreglo%AH*horas_rad)*tan(arreglo%Dec*gra_rad) ! Pone en la segunda columna los errores en AH despues de las correcciones conocidas

matriz(:,2) = arreglo%dAH ! En la tercera columna pone el termino de falta de alineacion vertical matriz(:,3) = tan(arreglo%Dec*gra_rad)*sin(arreglo%AH*horas_rad); ! Limita los valores pues la tangente puede hacerse muy grande where(matriz(:,3) > 1000) matriz(:,3) = 1000 where(matriz(:,3) < -1000) matriz(:,3) = -1000 ! Calcula la mejor recta que se ajusta a esta distribucion call ajusta_recta(matriz(:,3), matriz(:,2), m(1), b(1), correlacion(1)) write(*,'(A,3F14.6)') " No alineacion vertical: m, b y corr=", m(1), b(1), correlacion(1) ! Pregunta si aplica la correccion write(*,'(/a)', advance = "no") "Aplico esta correccion (s,n)? " read(*,*) resp if(resp == "s") then ! Corrige los errores usando la recta encontrada arreglo%dAH = arreglo%dAH-( m(1)*matriz(:,3) + b(1) ); PV = m(1) HI = HI + b(1) endif ! En la cuarta columna pone los errores con la correccion anterior matriz(:,4) = arreglo%dAH ! En la quinta columna pone el termino de falta de perpendicularidad matriz(:,5) = tan(arreglo%Dec*gra_rad); ! Tambien limitamos los valores de la tangente where(matriz(:,5) > 1000) matriz(:,5) = 1000 where(matriz(:,5) < -1000) matriz(:,5) = -1000 ! Calcula la mejor recta que se ajusta a esta distribucion call ajusta_recta(matriz(:,5), matriz(:,4), m(1), b(1), correlacion(1)) write(*,'(A,3F14.6)') " No perpendicularidad: m, b y corr=", m(1), b(1), correlacion(1) ! Pregunta si aplica la correccion write(*,'(/a)', advance = "no") "Aplico esta correccion (s,n)? " read(*,*) resp if(resp == "s") then ! Corrige los errores usando la recta encontrada arreglo%dAH = arreglo%dAH-( m(1)*matriz(:,5) + b(1) ); NP = m(1) HI = HI + b(1)

endif ! En la columna 6 pone los errores con la correccion anterior matriz(:,6) = arreglo%dAH ! En la columna 7 pone el termino de falta de colimacion matriz(:,7) = 1.0/cos(arreglo%Dec*gra_rad); where(matriz(:,7) > 1000) matriz(:,7) = 1000 where(matriz(:,7) < -1000) matriz(:,7) = -1000 ! Calcula la mejor recta que se ajusta a esta distribucion call ajusta_recta(matriz(:,7), matriz(:,6), m(1), b(1), correlacion(1)) write(*,'(A,3F14.6)') " Colimacion: m, b y corr=", m(1), b(1), correlacion(1) ! Pregunta si aplica la correccion write(*,'(/a)', advance = "no") "Aplico esta correccion (s,n)? " read(*,*) resp if(resp == "s") then ! Corrige los errores usando la recta encontrada arreglo%dAH = arreglo%dAH-( m(1)*matriz(:,7) + b(1) ) CO = m(1) HI = HI + b(1) endif ! En la columna 8 pone los errores con todas las correcciones matriz(:,8) = arreglo%dAH ! Aqui los datos estan acomodados como sigue: ! Col 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 ! ! ---------------------------------------------------------------------------! ! dAH | dAH_ | PV | dAH_PV | NP | dAH_PV_NP | CO | dAH_PV_NP_CO | AH ! call guarda_lineas("prov_ah.txt", matriz) call call call call call call

prepara_archivo_de_grafica() agrega_grafica(9,1,'"AH"', '"dAH"', agrega_grafica(3,2,'"PV"', '"dAH"', agrega_grafica(5,4,'"NP"', '"dAH"', agrega_grafica(7,6,'"CO"', '"dAH"', inicia_grafica()

end subroutine encuentra_coeficientes_ah end module utilerias

0, 1, 0, 1,

1,'"Errores iniciales"', "prov_ah.txt") 1,'"Alineacion vertical"', "prov_ah.txt") 0,'"Falta de perpendicularidad"', "prov_ah.txt") 0,'"Falta de colimacion"', "prov_ah.txt")

!*************************************************************! ! Empieza el programa ! !*************************************************************! Program calcula_modelo ! Trate de poner lo mas posible en subrutinas en el modulo a fin de hacer el programa ! mas compacto. use utilerias

! En este modulo estan las subrutinas

implicit none character(len=40)

:: archivo

integer

:: n_lineas

real

:: rmsdec, rmsah, rmsdec_ant, rmsah_ant, dummy

!

! Para el nombre del archivo que ! contiene los datos ! Numero de lineas en el archivo

Primero pide al usuario el nombre del archivo con los datos

write(unit=*, fmt='(a)', advance="no") "Teclea el nombre del archivo con los datos: " read(*,*) archivo ! do

Pide al usuario la latitude del observatorio y verifica que sea un numero razonable write(unit=*, fmt='(/a,f4.1,a)', advance="no") "Dame la latitud del observatorio en grados (", phi/gra_rad, "): " read(*,*) dummy if(dummy > -80.0 .and. dummy < 80.0) exit

enddo phi senphi cosphi ! !

= dummy*gra_rad = sin(phi) = cos(phi)

! la convierte en radianes ! Esto es el seno de la latitud ! Esto es el coseno de la latitud

Luego calcula el numero de renglones del archivo para definir los arreglos

n_lineas = cuantas_lineas(archivo)

! Verifica que exista y revisa cuantas lineas ! con datos tiene

print*, "Hay", n_lineas, "lineas validas" ! Si no hay lineas validas, termina el proceso

if(n_lineas == 0) return allocate(arreglo(n_lineas)) allocate(matriz(n_lineas,9))

! Acomoda el arreglo de n_lineas ! Acomoda la matriz donde pondra ! los datos a graficar

call trae_lineas(archivo, arreglo)

! Trae todos los datos y los pone en 'arreglo'

! Luego verifica que AH este entre -12 y 12 where(arreglo%AH > 12.0) arreglo%AH = arreglo%AH-24.0; call grafica_errores('"Errores iniciales"', rmsdec_ant, rmsah_ant) ! Aqui llama a una subrutina que calcula los coeficientes de correccion para la ! Declinacion call encuentra_coeficientes_dec() ! Luego llama a otra subrutina para calcular los coeficientes de correccion para ! el Angulo Horario call encuentra_coeficientes_ah() call grafica_errores('"Errores finales"', rmsdec, rmsah) ! Muestra los coeficientes calculados write(*,'(/A)') "Los coeficientes calculados son:" write(*,'(/A,F8.3)') " (Alineacion horizontal) - PH: ", PH write(*,'(A,F8.3)') " (Alineacion vertical) - PV: ", PV write(*,'(A,F8.3)') " (Flexion del tubo) - FT: ", FT write(*,'(A,F8.3)') " (Flexion del yugo) - FY: ", PV_FY - PV write(*,'(A,F8.3)') "(Falta de perpendicularidad) - NP: ", NP write(*,'(A,F8.3)') " (Error de colimacion) - CO: ", CO write(*,'(A,F8.3)') " (Backlash en AH) - BL: ", BL write(*,'(A,F8.3)') " (Offset en AH) - HI: ", HI write(*,'(A,F8.3)') " (Offset en Dec) - DI: ", DI write(*,'(/A)') "Con esta correccion se han bajado del errores RMS como sigue:" write(*,'(A,F5.1,A,F5.1,A)') " AH: de ", rmsah_ant, " a ", rmsah, " (arcsec)" write(*,'(A,F5.1,A,F5.1,A)') "Dec: de ", rmsdec_ant, " a ", rmsdec, " (arcsec)" ! Pregunta si cerramos los graficos write(*,'(/a)', advance = "no") "Cierro los graficos (s,n)? " read(*,*) resp

if(resp == "s") then call system("killall gnuplot_x11") ! Cierra todas las ventanas con gnuplot endif deallocate(arreglo) deallocate(matriz) End Program calcula_modelo !************************************************************************

Apéndice B. Listado del programa en IDL FUNCTION rms, vector media = mean(vector) va = vector - media return, sqrt(mean(va*va)) END

PRO encuentra_modelo, ARCHIVO=archivo, LATITUD=latitud, CORRIGE=corrige ;********************************************************************************** ; Programa 'encuentra_modelo.pro' elaborado por Leonel Gutierrez . ; ; Este programa lee un archivo con una serie de datos que muestran ; los errores de apuntado de un telescopio con montura ecuatorial ; y calcula los coeficientes de un modelo trigonometrico para corregir ; dicho apuntado. ;********************************************************************************** ; sintaxis: ENCUENTRA_MODELO, archivo="epos.txt", latitud=19.5, corrige=[1,1,1,0,1,0,1] if (NOT keyword_set(ARCHIVO)) then archivo="epos.txt" if (NOT keyword_set(LATITUD)) then latitud=19.5 if (NOT keyword_set(CORRIGE)) then begin corrigePH=1 corrigePV=1 corrigeFT=1 corrigePV_FY=1 corrigeNP=1 corrigeCO=1 corrigeBL=1 endif else begin if n_elements(corrige) ne 7 then begin print, "Corrige debe tener 7 elementos en el orden [PH,PV,FT,PV_FY,NP,CO,BL]." print, "Poner 1 en la correccion que se desea y 0 en la que no." return endif corrigePH=corrige(0) corrigePV=corrige(1)

corrigeFT=corrige(2) corrigePV_FY=corrige(3) corrigeNP=corrige(4) corrigeCO=corrige(5) corrigeBL=corrige(6) endelse ;Trae los datos del archivo readcol,archivo,f='(f,f,f,f)', ah, dah, dec, ddec ; Inicia una ventana para desplegar las graficas window, 1, xsi=300,ysi=300,retain=2, title="Errores iniciales" white = !D.N_COLORS-1 !P.color = 0 !p.background = white ; Desplegara solo una grafica !p.multi = [0,1,1,0,0] sx = "dDec (rms=" +string(rms(dDec),f='(f6.2)') + ")" sy = "dAH (rms=" +string(rms(dAH),f='(f6.2)') + ")" plot,ddec,dah,psym=4,symsize=0.5,xtitle=sx, ytitle=sy print, "rms de dAH y dDec (antes): ", rms(dAH), rms(dDec) ; Mi definicion de PI pi=acos(-1) ; Pone la latitud y calculas el seno y el coseno phi=latitud senphi = sin(phi*pi/180.0) cosphi = cos(phi*pi/180.0) ; Inicia otra ventana para desplegar el proceso de correccion en declinacion window, 2, xsi=600,ysi=600,retain=2, title="Correcciones en declinacion" ; Desplegara 4 graficas !p.multi = [0,2,2,0,0] ; Para calcular el termino correspondiente a las flexiones del tubo term1 = -senphi*cos(dec*pi/180.0) + cosphi*sin(dec*pi/180.0)*cos(ah*15*pi/180.0) ; Grafica los errores plot,term1,ddec,psym=4,symsize=0.5,xtitle='FT', ytitle='dDec' xyouts, 120,260+300,"Flexiones del tubo", /device ; Ajusta una recta y calcula FT params = fltarr(2)

params = linfit(term1,ddec) DI = params(0) FT = 0.0 if corrigeFT then begin FT = params(1) ; Corrige los errores ddec = ddec - (FT*term1 + DI) endif ; Para corregir las flexiones del yugo y la falta de alineacion vertical term2 = cos(ah*15*pi/180.0) plot, term2,ddec,psym=4,symsize=0.5,xtitle='PV+FY', ytitle='dDec' xyouts, 120+300,260+300,"Flex. yugo + alin. vert.", /device ; Ajusta una recta y calcula PV_FY params = linfit(term2,ddec) DI = params(0) + DI PV_FY = 0.0 if (corrigePV_FY) then begin PV_FY = params(1) ; Corrige los errores ddec = ddec-(PV_FY*term2 + params(0)) endif ; Para corregir la falta de alineacion horizontal term3 = sin(ah*15*pi/180.0) plot, term3,ddec,psym=4,symsize=0.5,xtitle='PH', ytitle='dDec' xyouts, 120,260,"Alin. horizontal", /device ; Ajusta una recta y calcula FT params = linfit(term3,ddec) DI = params(0) + DI PH = 0.0 if (corrigePH) then begin PH = params(1) ; Corrige los errores ddec = ddec-(PH*term3 + params(0)) endif ; Grafica los errores en declinacion una vez corregidos

plot,term1,ddec,psym=4,symsize=0.5,xtitle='FT', ytitle='dDec' xyouts, 120+300,260,"Errores corregidos en Dec", /device

; Inicia otra ventana para desplegar el proceso de correccion en AH window, 3, xsi=600,ysi=600,retain=2, title="Correcciones en AH" ; Grafica los errores plot,AH,dAH,psym=4,symsize=0.5,xtitle='AH', ytitle='dAH' xyouts, 120,260+300,"Errores iniciales", /device ; Para calcular el Backlash, calcula la media de los errores donde AH > 0 crit = where(ah>0) mediamas = mean(dAh(crit)) ; Luego calcula la media del lado negativo crit1 = where(ah