GALAC v1.02

        Como respuesta a las diferentes practicas planteadas en la 
asignatura optativa del octavo quatrimestre de la carrera de Fisica de 
la Universidad de Barcelona,Astrofisica Galactica,GALAC es la evolucion 
logica de GALAC1.

NUESTRAS INTENCIONES

        En la primera practica de esta asignatura se nos pide que hagamos
un estudio de la Galaxia proxima al sistema solar a partir de unos datos 
contenidos en un fichero (velocidades,posiciones,movimientos propios,
estadisticas,etc.);por lo tanto,GALAC esta pensado para leer ficheros con 
el formato del fichero de la practica (lo que en el programa llamamos es-
tandar TOTAL2.dat).
        Este formato es:

200     FORMAT(I4,A1,I2,I5,2I3,I2,I3,I4,A1,I2,I4,3I2,I3,A1,2I2,I4,I4,
     +  A1,I2,I4,A1,I2,I4)
300     FORMAT(2I5,A1,I3,I4,A1,I3,I4,A1,I3,I4,2I4,I2,I3,I2)

donde los datos de una estrella se introducen en dos filas.
        En la primera fila se nos dan los numeros de catalogo BS (4 digitos),
DM (un signo,dos digitos para la zona y cinco para el numwero),HD (dos numeros 
de tres digitos),GC (dos numeros de dos y tres digitos),FK4 (4 digitos),AGK2
(un signo,dos digitos y cuatro digitos),la ascension recta (A,alfa;h,m,s,mi-
lesimas de segundo),la declinacion (D,delta;signo,,'," y centesimas de segun-
do),epoca (4 digitos),movimiento propio en alfa (AM,por el cosinus de alfa;
signo,"/ao,diezmilesimas de "/ao) y movimiento propio en delta (DM;idem);
en la segunda fila se da paralaje (IP;diezmilesimas de "),velocidad radial
(VR o Vr;en decimas de Km/s),U,V,W (las tres con signo,tres digitos,cuatro 
digitos;en Km/s,diezmilesimas de Km/s),magnitud visual (m;centesimas de mag-
nitud),B-V (idem),tipo espectral (TS;dos digitos para el tipo y tres para el
subtipo),clase de luminosidad (CL o LC;dos digitos).
        Como este formato es complicado,muchos de los datos no son de inte-
res y nos interesa un formato y unas unidades mas manejables para operar
GALAC utiliza el formato estandar en GALAC1

500     FORMAT(1X,I3,1X,I3,1X,F11.7,1X,F10.6,1X,F8.4,1X,F8.4,1X,F9.4,
     +  1X,I5,1X,F9.4,1X,F9.4,1X,F9.4,1X,I4,1X,I4,1X,I2,1X,I3,1X,I2,1X,
     +  F10.6,1X,F10.6,1X,F9.4,1X,F8.4,1X,F10.4,1X,F10.4,1X,F10.4)

formato utilizado en todo el programa.
        Este formato es el que corresponde para dar numero de catalogo HD  
(puede haber estrellas en las que no aparezca;las dos primeras columnas),
A (en h),D (en ),AM (en "/ao),DM (idem),R (o r,distancia radial;en pc),
Vr (en decimas de Km/s),U,V,W (en Km/s),m (en centesimas de magnitud),B-V 
(idem),TS (tipo y subtipo,ver mas abajo;en dos columnas),CL (ver mas abajo),
l (longitud galactica;en ),b (latitud galactica;en ),movimientos propios en
l y b (lM,bM;en "/ao),X,Y y Z (en pc).
        Por otra parte,en la segunda practica de la asignatura se nos pide
que hagamos un estudio cinematico del entorno solar;para ello,GALAC incorpora
a las prestaciones de GALAC1 una opcion para realizar este estudio,ademas de
otras mejoras y opciones.
        Asi pues,nuestra intencion ha sido hacer un programa que nos solu-
cione lo que necesitamos para la susodicha practica,pero intentando a la 
vez crear un programa mucho mas general.


MANUAL DE USO
                   
        1)CALCULOS
                
                Escogiendo la opcion 1) del programa entramos en la
        parte del programa que nos calcula l (longitud galactica),b
        (latitud galactica),movimientos propios en l y b,posiciones
        cartesianas (X,Y y Z) y velocidades referidas a este sitema 
        cartesiano (U,V y W).
                Para este apartado es imprescindible introducir un 
        fichero con el formato TOTAL2.dat;si no es asi,el programa 
        no funcionara.
                A la salida nos proporcionara un fichero en el forma-
        to estandar de GALAC1 donde aparecera:numero de catalogo HD 
        (las dos primeras columnas),alfa(en h),delta (en ),movimientos 
        propios en alfa y delta (en "/ao),distancia radial (en pc),
        velocidad radial (en 1e-1 Km/s),U,V,W (en Km/s),magnitud visual 
        (en 1e-2 mag.),indice de color (en 1e-2 mag.),tipo espectral,
        subtipo espectral (segun codigos indicados mas abajo),clase de 
        luminosidad (segun codigos indicados mas abajo),l,b (en ),mo-
        vimientos propios en l y b (en "/ao),X,Y y Z (en pc).
                Los codigos de tipos espectrales (y subtipos) son:

          TS     W  O  B  A  F  G  K  M  N  R  S  C
                 0  1  2  3  4  5  6  7  8  9  10 11

          sTS    0  1   2   3   4   5   6   7   8   9   no se conoce
                 0  10  20  30  40  50  60  70  80  90      100

                Y los codigos para clase de luminosodad son:

                 I   II  III IV  V   VI  wd  no se conoce
                 10  20  30  40  50  60  61       99
                   15  25  35  45  55

                Un  999 en el indice de color quiere decir que no esta 
        determinado.
                Para realizar los calculos se han empleado constantes
        (posicion del polo norte galactico,...) referidas al 1950.0;por
        lo tanto,segun la epoca en que se tomaran los datos puede ser ne-
        cesario coregirlos (mirar seccion EL PROGRAMA o ponerse en con-
        tacto con nosotros).

        2)ORDENACIONES

                La opcion 2) nos clasifica las estrellas de ficheros 
        creados con la opcion 1),con la opcion 2) o con el formato 
        estandar en GALAC1,en ficheros segun tipo espectral (O,B,A,F,
        G,K,M,OTROS con extension .DAT), clase deluminosidad (I,I-II,II,
        II-III,III,III-IV,IV,IV-V,V,V-IV,VI,OTRAS con extension .DAT) o
        indice de color (10 ficheros donde se nos presentan las estre-
        llas que cumplen B-Vmin+paso*numero del fichero-1< B-V* >B-Vmin+
        paso*numero del fichero ;los ficheros son B-V1,B-V2,B-V3,B-V4,
        B-V5,B-V6,B-V7,B-V8,B-V9,B-V10 con extension .DAT).
                Devemos advertir que los nombres de estos ficheros son 
        fijos,el usuario no puede definirlos;por lo tanto,si se quieren 
        conservar,el usuario devera cambiarles el nombre "a mano" (por
        ejemplo:supongamos que hemos creado los ficheros O.DAT,etc;ahora
        creamos los ficheros I.DAT,etc. a partir del G.DAT -por ejemplo-;
	  si despues queremos calcular los ficheros de clase de lumi-
        nosidad para otro fichero,estos se llamaran tambien I.DAT...con 
        lo cual los ficheros existentes se nos borraran).

        3)LIMPIEZA

                Con esta opcion GALAC nos permite eliminar estrellas
        que no tienen definido un parametro de un fichero para despues
        trabajar con el.
                Estos parametros no definidos pueden ser el subtipo es-
        pectral,la clase de luminosidad o el indice de color.
                Solo se limpia un parametro cada vez que se usa la op-
        cion;si se desean eliminar mas parametros del  fichero se deven
        hacer ejecuciones iteradas de la opcion.
                Ahora bien,puede no resultar interesante eliminar 2 o
        3 paramtros porque nos podriamos quedar con muy pocas estrellas
        (estrellas que no tienen indice de color conocido pueden tener
        subtipo espectral y/o clase de luminosidad,por ejemplo).

        4)ESTADISTICA

                La opcion 1) nos permite contar el numero de estrellas
        que responden a una condicion concreta (n de estrellas de tipo
        espectral O,a distancia R,con velocidad V...).
                Obviamente,segun que dato queramos estudiar es conve-
        niente antes "limpiar" el fichero para evitar errores (GALAC
        calcula N( ) dividiendo el intervalo de valores que nosotros le 
        damos en N1 subintervalos -ver NOTA-; aun asi,en algunos casos 
        puede ser necesario).Por lo tanto,GALAC nos da N( ) discretizado;
        de todas formas,podemos hacer la discretizacion tan fina como 
        queramos (o nos permita el dato).
                NOTA:Cada iteracion de este apartado tarda aproximadamente
        1 segundo en un PENTIUM MMX a 200 MHz;por lo tanto,recomendamos 
        tener cuidado al escoger los limites y el paso.
                Tambien incorpora una opcion 2) que nos permite calcular
        la media y la desviacion estandar de r,Vr,U,V o W.El programa tam-
	  bien nos permite restringir nuestro intervalo de estudio;pero hacer 
	  esto nos puede reducir el numero de estrellas consideradas;por eso
	  GALAC nos dice el numero de estrellas que ha considerado.
                Debemos advertir que en esta opcion el programa realiza 
        dos lecturas del fichero de entrada (por eso la segunda vez que 
        aparece el indicador TRABAJANDO... el numero total de iteraciones
        sera el doble del numero de entradas indicado -ver para mas acla-
        racion la parte del programa correspondiente mas abajo-).

        5)CALCULO DEL POLO

                Con esta opcion el programa calcula la posicion del polo
        norte galactico mediante el metodo de Newcomb.
                Como es de esperar tratandose de un ajuste,los resultados
        asi encontrados seran mejores como mejores sean los datos utili-
        zados.

        6)ESTUDIO CINEMATICO DEL ENTORNO SOLAR
                
                En esta opcion GALAC nos calcula los valores de las varian-
        cias (o las raices cuadradas de estas para los terminos diagonales 
        de la matriz) con sus errores.
                Tambien calcula la desviacion del vertice del elipsoide de
        velocidades con su error.
		    Como en la opcion de calculo de medias del apartado 4),GALAC
	  tambien nos da el numero de estrellas consideradas al hacer los cal-
	  culos.
                En esta version de GALAC introducimos la posibilidad de ob-
        tener los valores de los semiejes principales de este elipsoide,si
        bien todabia no nos proporciona su error;aun asi,hemos considerado 
        interesante darlos,ya que puede ser de utilidad a la hora de repre-
        sentaciones graficas o similares.

EL PROGRAMA

        Por si algun usuario deseara completar o mejorar el programa 
(los valores de algunas constantes deven actualizarse con el tiempo,las
opciones pueden ampliarse...) a continuacion ofrecemos el programa para
su ejecucion en FORTRAN.

c       Programa para las practica de la asignatura de Astrofisica 
c       Galactica
        PROGRAM GALAC
        IMPLICIT REAL*8 (A-H,O-Z)
        PARAMETER (PI=3.14159)
        CHARACTER*1 DMCS,AGK2S,DS,AMS,DMS,US,VS,WS,RESP
        CHARACTER*11 NOM1,NOM2

10      CONTINUE
c       Preguntamos que parte del programa queremos usar
        WRITE (*,*) 'QUE QUIERE HACER?'
        WRITE (*,*) '1)Calcular l,b,movimientos propios...;2)Ordenar en  
     +ficheros segun clase de     luminosidad (CL),tipo espectral (TS) o
     + indice de color B-V (IC);3)Limpiar un   fichero;4)Estadistica;5)C
     +alculo del polo norte galactico;6)Estudio cinematico  del entorno 
     +solar'
        READ (*,*) IRESP
        IF (IRESP.EQ.2) GO TO 20
        IF (IRESP.EQ.3) GO TO 30
        IF (IRESP.EQ.4) GO TO 40
        IF (IRESP.EQ.5) GO TO 50
        IF (IRESP.EQ.6) GO TO 60
c       Esta es la parte de calculo
c       Pedimos el monbre con extension del fichero de entrada
        WRITE (*,*) 'ESCRIBA NOMBRE DEL FICHERO DE ENTRADA (con ext.):'
        READ (*,100) NOM1
100     FORMAT (A)
c       Pedimos el numero de datos que vamos a estudiar (n de estre-       
c       llas)
        WRITE(*,*) 'NUMERO DE ENTRADAS:'
        READ (*,*) N
c       Abrimos el fichero de datos y el de salida        
        OPEN (1,FILE=NOM1)
        WRITE (*,*) 'ESCRIBA NOMBRE DEL FICHERO DE SALIDA (con ext.):'
        READ (*,100) NOM2
        OPEN (2,FILE=NOM2)
c       Definimos una variable para controlar el paso de iteracion        
        ICOUNT=0
c       Iniciamos el bucle de lectura        
        DO I=1,N
c       Leemos la primera linea de la estrella i                
                READ(1,200) IBS,DMCS,IDMC1,IDMC2,IHD1,IHD2,IGC1,IGC2,
     +          IFK4,AGK2S,IAGK21,IAGK22,IA1,IA2,IA3,IA4,DS,ID1,ID2,
     +          ID3,IEPO,AMS,IAM1,IAM2,DMS,IDM1,IDM2
c       Leemos la segunda linea de la estrella i                
                READ(1,300) IP,IVR,US,IU1,IU2,VS,IV1,IV2,WS,IW1,IW2,
     +          IVM,IC,ITS1,ITS2,LC
c       Llamamos a la subrutina de calculo de l y b                
                CALL COORD(IA1,IA2,IA3,IA4,DS,ID1,ID2,ID3,GL,B,AR,DR,
     +          SF,CF)
c       Llamamos a la subrutina de calculo de movimientos propios en
c       l y b
                CALL MOV(GL,B,AR,DR,AMS,IAM1,IAM2,DMS,IDM1,IDM2,AM,
     +          DM,GLM,BM)
c       Llamamos a la subrutina de calculo de X,Y y Z                
                CALL XYZ(IP,GL,B,X,Y,Z,R)
c       Llamamos a la subrutina de calculo de U,V y W                
                CALL UVW(IVR,GL,B,R,AM,DM,SF,CF,U,V,W)
c       Pasamos las coordenadas angulares a un sistema de unidades mas                
c       convencional
                A=(AR*360.0/(2*PI))/15.0
                D=DR*360.0/(2*PI)
                GL=GL*360.0/(2*PI)
                B=B*360.0/(2*PI)
c       Escrivimos los resultados y datos de interes en el fichero de                 
c       salida
                WRITE (2,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
c       Indicamos el paso en el que se encuentra el programa        
                ICOUNT=ICOUNT+1
                WRITE (*,400) 'TRABAJANDO...',ICOUNT
400     FORMAT ('+',A13,1X,I5)
        END DO
c       Formatos del fichero de entrada (Segun estandar TOTAL2.DAT)
200     FORMAT(I4,A1,I2,I5,2I3,I2,I3,I4,A1,I2,I4,3I2,I3,A1,2I2,I4,I4,
     +  A1,I2,I4,A1,I2,I4)
300     FORMAT(2I5,A1,I3,I4,A1,I3,I4,A1,I3,I4,2I4,I2,I3,I2)
c       Formato del fichero de salida (estandar usado por GALAC1)
500     FORMAT(1X,I3,1X,I3,1X,F11.7,1X,F10.6,1X,F8.4,1X,F8.4,1X,F9.4,
     +  1X,I5,1X,F9.4,1X,F9.4,1X,F9.4,1X,I4,1X,I4,1X,I2,1X,I3,1X,I2,1X,
     +  F10.6,1X,F10.6,1X,F9.4,1X,F8.4,1X,F10.4,1X,F10.4,1X,F10.4)
c       Cerramos los ficheros         
        CLOSE (1)
        CLOSE (2)
        GO TO 70
20      CONTINUE
c       Parte del programa que ordena las estrellas en ficheros segun        
c       tipo espectral o clase de luminosodad
        CTL=1
        CALL ORD(CTL,CTL1)
        GO TO 70
30      CONTINUE
c       Parte del programa que elimina de un fichero las estrellas que        
c       no tienen un determinado dato definido
        CTL=1
        CALL LIMPIA(CTL,CTL1)
        GO TO 70
40      CONTINUE
c       Parte del programa que cuenta el numero de estrellas con una         
c       cierta condicion en un fichero dado
        CTL=1
        CALL CONTA(CTL,CTL1)
        GO TO 70
50      CONTINUE
c       Parte del programa que calcula la posicion del polo norte galac-
c       tico (AG y DG) haciendo un ajuste
        CTL=1
        CALL PNG(CTL,CTL1)
        GO TO 70
60      CONTINUE
c       Parte del programa que nos permite hacer un estudio cinematico 
c       del entorno solar
        CTL=1
        CALL ECES(CTL,CTL1)
70      CONTINUE
c       Preguntamos si queremos abandonar el programa        
        WRITE (*,*) 'QUIERE FINALIZAR EL PROGRAMA?(S/N)'
        READ (*,100) RESP
        IF (RESP.EQ.'n'.OR.RESP.EQ.'N') GO TO 10
        STOP
        END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

c       Subrutina de calculo de coordenadas (tambien nos da cosinus        
c       y sinus del angulo paralactico
        SUBROUTINE COORD(IA1,IA2,IA3,IA4,DS,ID1,ID2,ID3,GL,B,AR,DR,
     +  SF,CF)
        IMPLICIT REAL*8 (A-H,O-Z)
c       Estos parametros se deven actualizar si son demasiado viejos        
        PARAMETER (PI=3.141592654,DG=0.478220215,AG=3.355395489,
     +  GLE=2.14675498)
        CHARACTER DS*1

c       Pasamos alfa (A) y delta (D) a radianes        
        A=15.0*(IA1+IA2/60.0+IA3/3600.0+IA4*1D-3/3600.0)
        D=ID1+ID2/60.0+(ID3*1D-2)/3600.0
        IF (DS.EQ.'-') D=-D
        AR=A*2.0*PI/360.0
        DR=D*2.0*PI/360.0
c       Calculamos b y l (GL;en el cuadrante correcto)        
        B=ASIN(COS(DR)*COS(DG)*COS(AR-AG)+SIN(DR)*SIN(DG))
        V1=ACOS((-COS(DR)*SIN(DG)*COS(AR-AG)+SIN(DR)*COS(DG))/
     +  COS(B))
        V2=ASIN(-COS(DR)*SIN(AR-AG)/COS(B))
        IF (V1.GT.0.0.AND.V1.LT.PI/2.AND.V2.GT.0.0.AND.V2.LT.PI/2)
     +  GL=V1+GLE
        IF (V1.GT.PI/2.AND.V1.LT.PI.AND.V2.GT.0.0.AND.V2.LT.PI/2)
     +  GL=V1+GLE
        IF (V1.GT.0.0.AND.V1.LT.PI/2.AND.V2.GT.-PI/2.AND.V2.LT.0.0)
     +  GL=2.0*PI+V2+GLE
        IF (V1.GT.PI/2.AND.V1.LT.PI.AND.V2.GT.-PI/2.AND.V2.LT.0.0)
     +  GL=PI-V2+GLE
        IF (V1.EQ.0.0.AND.V2.EQ.0.0) GL=V1+GLE
        IF (V1.EQ.PI/2.AND.V2.EQ.PI/2) GL=V1+GLE
        IF (V1.EQ.PI.AND.V2.EQ.0.0) GL=V1+GLE
        IF (V1.EQ.PI/2.AND.V2.EQ.-PI/2) GL=2.0*PI+V2+GLE
        IF (GL.GT.2.0*PI) GL=GL-2.0*PI
c       Calculamos el sinus (SF) y el cosinus (CF)del angulo paralactico        
        SF=-COS(DG)*SIN(GL-GLE)/COS(DR)
        CF=(-COS(DG)*SIN(B)*COS(GL-GLE)+SIN(DG)*COS(B))/COS(DR)
        RETURN
        END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

c       Subrutina de calculo de movimientos propios en l y b        
        SUBROUTINE MOV(GL,B,AR,DR,AMS,IAM1,IAM2,DMS,IDM1,IDM2,AM,DM,
     +  GLM,BM)
        IMPLICIT REAL*8 (A-H,O-Z)
        PARAMETER (DG=0.4734772828,AG=3.366032541,GLE=2.133871541)
        CHARACTER*1 AMS,DMS

c       Pasamos movimiento propio en alfa (AM) y en delta (DM) a         
c       unidades mas convenientes
        AM=IAM1+IAM2*1D-4
        IF (AMS.EQ.'-') AM=-AM
        DM=IDM1+IDM2*1D-4
        IF (DMS.EQ.'-') DM=-DM
c       Calculamos el movimiento propio en b (BM) y en l (GLM)        
        BM=(DM*(COS(DR)*SIN(DG)-SIN(DR)*COS(DG)*COS(AR-AG))-
     +  AM*COS(DG)*SIN(AR-AG))/COS(B)
        GLM=(DM*SIN(DR)*SIN(AR-AG)-AM*COS(AR-AG)-BM*COS(AR-AG)-
     +  BM*SIN(B)*SIN(GL-GLE))/COS(GL-GLE)
        RETURN
        END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        
c       Subrutina de calculo de X,Y y Z        
        SUBROUTINE XYZ(IP,GL,B,X,Y,Z,R)
        IMPLICIT REAL*8 (A-H,O-Z)

c       Calculamos la distancia radial (R) a partir de la paralaje (IP)        
        R=1/(IP*1D-4)
c       Calculamos X,Y y Z        
        X=R*COS(B)*COS(GL)
        Y=R*COS(B)*SIN(GL)
        Z=R*SIN(B)
        RETURN
        END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        
c       Subrutina de calculo de U,V y W        
        SUBROUTINE UVW(IVR,GL,B,R,AM,DM,SF,CF,U,V,W)
        IMPLICIT REAL*8 (A-H,O-Z)
c       Parametro K en las unidades adecuadas         
        PARAMETER (XK=4.741)

c       Pasamos la velocidad radial (VR) a Km/s        
        VR=IVR*1D-1
c       Calculamos U,V y W        
        U=VR*COS(GL)*COS(B)+XK*R*AM*(COS(GL)*SIN(B)*SF-SIN(GL)*CF)-
     +  XK*R*DM*(COS(GL)*SIN(B)*CF+SIN(GL)*SF)
        V=VR*SIN(GL)*COS(B)+XK*R*AM*(SIN(GL)*SIN(B)*SF+COS(GL)*CF)-
     +  XK*R*DM*(SIN(GL)*SIN(B)*CF-COS(GL)*SF)
        W=VR*SIN(B)-XK*R*AM*COS(B)*SF+XK*R*DM*COS(B)*CF
        RETURN
        END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        
c       Subrutina de ordenacion de estrellas        
c       Utilizamos una subrutina para marcar claramente que es una 
c       parte independiente del programa (se podria hacer un programa
c       independiente de GALAC1)
        SUBROUTINE ORD(CTL,CTL1)
        IMPLICIT REAL*8 (A-H,O-Z)
        CHARACTER NOM1*11,RESP*1

70      CONTINUE
c       Pedimos el monbre con extension del fichero de entrada
        WRITE (*,*) 'ESCRIBA NOMBRE DEL FICHERO DE ENTRADA (con ext.):'
        READ (*,100) NOM1
100     FORMAT (A)
c       Pedimos el numero de datos que vamos a estudiar (n de estre-       
c       llas)
        WRITE(*,*) 'NUMERO DE ENTRADAS:'
        READ (*,*) N
        OPEN (1,FILE=NOM1)       
c       Pedimos el criterio de clasificacion        
        WRITE(*,*) 'COMO QUIERE LOS FICHEROS DE SALIDA?' 
        WRITE(*,*) '(1 por TS , 2 por CL y 3 por IC)'
        READ (*,*) IRESP
        IF (IRESP.EQ.2) GO TO 20  
        IF (IRESP.EQ.3) GO TO 30
c       Abrimos el fichero de datos y los de salida para TS
        OPEN (2,FILE='O.DAT')
        OPEN (3,FILE='B.DAT')
        OPEN (4,FILE='A.DAT')
        OPEN (7,FILE='F.DAT')
        OPEN (8,FILE='G.DAT')
        OPEN (9,FILE='K.DAT')
        OPEN (10,FILE='M.DAT')
        OPEN (11,FILE='OTROS.DAT')
c       Definimos una variable para controlar el paso de iteracion        
        ICOUNT=0
c       Iniciamos el bucle de lectura        
        DO I=1,N
                READ (1,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
c       Clasificamos la estrella con sus resultados en un TS                
                IF (ITS1.EQ.1) WRITE (2,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                IF (ITS1.EQ.2) WRITE (3,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                IF (ITS1.EQ.3) WRITE (4,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                IF (ITS1.EQ.4) WRITE (7,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                IF (ITS1.EQ.5) WRITE (8,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                IF (ITS1.EQ.6) WRITE (9,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                IF (ITS1.EQ.7) WRITE (10,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                IF (ITS1.NE.1.AND.ITS1.NE.2.AND.ITS1.NE.3.AND.ITS1.NE.
     +          4.AND.ITS1.NE.5.AND.ITS1.NE.6.AND.ITS1.NE.7)
     +          WRITE (11,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z  
c       Indicamos el paso en el que se encuentra el programa        
                ICOUNT=ICOUNT+1
                WRITE (*,400) 'TRABAJANDO...',ICOUNT
        END DO
400     FORMAT ('+',A13,1X,I5)
500     FORMAT(1X,I3,1X,I3,1X,F11.7,1X,F10.6,1X,F8.4,1X,F8.4,1X,F9.4,
     +  1X,I5,1X,F9.4,1X,F9.4,1X,F9.4,1X,I4,1X,I4,1X,I2,1X,I3,1X,I2,1X,
     +  F10.6,1X,F10.6,1X,F9.4,1X,F8.4,1X,F10.4,1X,F10.4,1X,F10.4)
c       Cerramos los ficheros         
        CLOSE (2)
        CLOSE (3)
        CLOSE (4)
        CLOSE (7)
        CLOSE (8)
        CLOSE (9)
        CLOSE (10)
        CLOSE (11)
        IF (IRESP.EQ.1) GO TO 60
20      CONTINUE
c       Clasificamos las estrellas en una CL
c       Abrimos los ficheros para CL
        OPEN (2,FILE='I.DAT')
        OPEN (3,FILE='II.DAT')
        OPEN (4,FILE='III.DAT')
        OPEN (7,FILE='IV.DAT')
        OPEN (8,FILE='V.DAT')
        OPEN (9,FILE='VI.DAT')
        OPEN (10,FILE='VII.DAT')
        OPEN (11,FILE='I-II.DAT')
        OPEN (12,FILE='II-III.DAT')
        OPEN (13,FILE='III-IV.DAT')
        OPEN (14,FILE='IV-V.DAT')
        OPEN (15,FILE='V-VI.DAT')
        OPEN (16,FILE='OTRAS.DAT')
c       Iniciamos el bucle de lectura,etc.        
        ICOUNT=0
        DO I=1,N
                READ (1,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                IF (LC.EQ.10) WRITE (2,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                IF (LC.EQ.20) WRITE (3,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                IF (LC.EQ.30) WRITE (4,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                IF (LC.EQ.40) WRITE (7,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                IF (LC.EQ.50) WRITE (8,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                IF (LC.EQ.60) WRITE (9,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                IF (LC.EQ.61) WRITE (10,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                IF (LC.EQ.15) WRITE (11,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                IF (LC.EQ.25) WRITE (12,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                IF (LC.EQ.35) WRITE (13,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                IF (LC.EQ.45) WRITE (14,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                IF (LC.EQ.55) WRITE (15,500) IHD1,IHD2,A,D,AM,
     +          DM,R,IVR,U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                IF (LC.NE.10.AND.LC.NE.20.AND.LC.NE.30.AND.LC.NE.40.
     +          AND.LC.NE.50.AND.LC.NE.60.AND.LC.NE.61.AND.LC.NE.15.
     +          AND.LC.NE.25.AND.LC.NE.35.AND.LC.NE.45.AND.LC.NE.55)
     +          WRITE (16,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z  
c       Indicamos el paso en el que se encuentra el programa        
                ICOUNT=ICOUNT+1
                WRITE (*,400) 'TRABAJANDO...',ICOUNT
        END DO
        CLOSE (2)
        CLOSE (3)
        CLOSE (4)
        CLOSE (7)
        CLOSE (8)
        CLOSE (9)
        CLOSE (10)
        CLOSE (11)
        CLOSE (12)
        CLOSE (13)
        CLOSE (14)
        CLOSE (15)
        CLOSE (16)
        GO TO 60
c       Clasificamos las estrellas en ficheros segun su indice de color
30      CONTINUE
        ICOUNT=0  
        XMAX=0.0
        XMIN=0.0
        DO I=1,N
                READ (1,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                IF (IC.GT.XMAX) XMAX=IC
                IF (IC.LT.XMIN) XMIN=IC
                ICOUNT=ICOUNT+1
                WRITE (*,400) 'TRABAJANDO...',ICOUNT
        END DO
        REWIND (1)
        WRITE (*,200) 'EL VALOR MAXIMO ES:',XMAX
        WRITE (*,200) 'EL VALOR MINIMO ES:',XMIN
200     FORMAT (1X,A19,1X,F10.4)
c       Pedimos los valores maximo y minimo a estudir 
        WRITE (*,*) 'VALOR MAXIMO A CONSIDERAR:'
        READ (*,*) XMAX
        WRITE (*,*) 'VALOR MINIMO A CONSIDERAR:'
        READ (*,*) XMIN
        PASO=(XMAX-XMIN)/10.0
        WRITE (*,300) 'EL INCREMENTO ENTRE UN FICHERO Y OTRO ES DE:',
     +  PASO
300     FORMAT (1X,A44,1X,F9.3)
        WRITE (*,*) ' '
c       Abrimos los ficheros        
        OPEN (2,FILE='B-V1.DAT')
        OPEN (3,FILE='B-V2.DAT')
        OPEN (4,FILE='B-V3.DAT')
        OPEN (7,FILE='B-V4.DAT')
        OPEN (8,FILE='B-V5.DAT')
        OPEN (9,FILE='B-V6.DAT')
        OPEN (10,FILE='B-V7.DAT')
        OPEN (11,FILE='B-V8.DAT')
        OPEN (12,FILE='B-V9.DAT')
        OPEN (13,FILE='B-V10.DAT')
c       Definimos una variable para controlar el paso de iteracion        
        ICOUNT=0
c       Iniciamos el bucle de lectura        
        DO I=1,N
                READ (1,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
c       Clasificamos la estrella con sus resultados en un intervalo en                
c       el indice de color B-V
                IF (IC.GE.XMIN.AND.IC.LT.XMIN+PASO) THEN
                        WRITE (2,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +                  IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                END IF
                IF (IC.GE.XMIN+PASO.AND.IC.LT.XMIN+2.0*PASO) THEN
                        WRITE (3,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +                  IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                END IF
                IF (IC.GE.XMIN+2.0*PASO.AND.IC.LT.XMIN+3.0*PASO) THEN
                        WRITE (4,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +                  IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                END IF
                IF (IC.GE.XMIN+3.0*PASO.AND.IC.LT.XMIN+4.0*PASO) THEN
                        WRITE (7,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +                  IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                END IF
                IF (IC.GE.XMIN+4.0*PASO.AND.IC.LT.XMIN+5.0*PASO) THEN
                        WRITE (8,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +                  IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                END IF
                IF (IC.GE.XMIN+5.0*PASO.AND.IC.LT.XMIN+6.0*PASO) THEN
                        WRITE (9,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +                  IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                END IF
                IF (IC.GE.XMIN+6.0*PASO.AND.IC.LT.XMIN+7.0*PASO) THEN
                        WRITE (10,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +                  IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                END IF
                IF (IC.GE.XMIN+7.0*PASO.AND.IC.LT.XMIN+8.0*PASO) THEN
                        WRITE (11,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +                  IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                END IF
                IF (IC.GE.XMIN+8.0*PASO.AND.IC.LT.XMIN+9.0*PASO) THEN
                        WRITE (12,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +                  IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                END IF
                IF (IC.GE.XMIN+9.0*PASO.AND.IC.LE.XMIN+10.0*PASO) THEN
                        WRITE (13,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +                  IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z           
                END IF
c       Indicamos el paso en el que se encuentra el programa        
                ICOUNT=ICOUNT+1
                WRITE (*,400) 'TRABAJANDO...',ICOUNT
        END DO
c       Cerramos los ficheros         
        CLOSE (2)
        CLOSE (3)
        CLOSE (4)
        CLOSE (7)
        CLOSE (8)
        CLOSE (9)
        CLOSE (10)
        CLOSE (11)
        CLOSE (12)
        CLOSE (13)
60      CONTINUE
        CLOSE (1)
c       Pedimos si se quiere regresar a GALAC
        WRITE (*,*) 'QUIERE CLASIFICAR OTRO FICHERO?(S/N)'
        READ (*,100) RESP
        IF (RESP.EQ.'S'.OR.RESP.EQ.'s') GO TO 70
c       Damos el valor a la variable de control        
        CTL1=CTL
        RETURN
        END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        
c       Subrutina de limpieza de ficheros        
c       Utilizamos una subrutina para marcar claramente que es una 
c       parte independiente del programa (se podria hacer un programa
c       independiente de GALAC1)
        SUBROUTINE LIMPIA(CTL,CTL1)
        IMPLICIT REAL*8 (A-H,O-Z)
        CHARACTER NOM1*11,NOM2*11,RESP*1

20      CONTINUE
c       Pedimos el fichero a limpiar y el numero de filas que contiene       
c       (numero de estrellas)
        WRITE (*,*) 'NOMBRE DEL FICHERO DE ENTRADA (con ext.):'
        READ (*,100) NOM1
100     FORMAT (A)
        WRITE (*,*) 'NUMERO DE ENTRADAS:'
        READ (*,*) N
c       Pedimos el fichero de salida        
        WRITE (*,*) 'NOMBRE DEL FICHERO DE SALIDA (con ext.):'
        READ (*,100) NOM2
c       Abrimos los ficheros        
        OPEN (1,FILE=NOM1)
        OPEN (2,FILE=NOM2)
c       Pedimos que dato se quiere limpiar de valores no determinados        
        WRITE (*,*) 'QUE VALOR QUIERE LIMPIAR?'
        WRITE (*,*) '(1 para Indice de Color,2 para Subtipo Espectral,3
     +para Clase de Luminosidad)'
        READ (*,*) IRESP
        ICOUNT=0
c       Iniciamos el bucle de lectura        
        DO I=1,N
                READ (1,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
c       Determinamos los criterios de limpieza                
                IF (IRESP.EQ.1) THEN 
                        IX=IC
                        IREF=999
                END IF
                IF (IRESP.EQ.2) THEN
                        IX=ITS2
                        IREF=100
                END IF
                IF (IRESP.EQ.3) THEN
                        IX=LC
                        IREF=99
                END IF
c       Limpiamos la fila;si la estrella i no tiene definido el valor                
c       a limpiar la fila de esta estrella no se considera 
                IF (IX.EQ.IREF) GO TO 10
                WRITE (2,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
10      CONTINUE
        ICOUNT=ICOUNT+1
        WRITE (*,400) 'TRABAJANDO...',ICOUNT
400     FORMAT('+',A13,1X,I5)
        END DO
500     FORMAT(1X,I3,1X,I3,1X,F11.7,1X,F10.6,1X,F8.4,1X,F8.4,1X,F9.4,
     +  1X,I5,1X,F9.4,1X,F9.4,1X,F9.4,1X,I4,1X,I4,1X,I2,1X,I3,1X,I2,1X,
     +  F10.6,1X,F10.6,1X,F9.4,1X,F8.4,1X,F10.4,1X,F10.4,1X,F10.4)
c       Cerramos los ficheros        
        CLOSE (1)
        CLOSE (2)
c       Preguntamos si se desea vober a GALAC1        
        WRITE (*,*) 'QUIERE LIMPIAR OTRO FICHERO?(S/N)'
        READ (*,100) RESP
        IF (RESP.EQ.'S'.OR.RESP.EQ.'s') GO TO 20
        CTL1=CTL
        RETURN
        END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        
c       Subrutina que cuenta el numero de estrellas con una determinada        
c       propiedad y permite hacer medias y desviaciones estandar
c       Utilizamos una subrutina para marcar claramente que es una 
c       parte independiente del programa (se podria hacer un programa
c       independiente de GALAC1)
        SUBROUTINE CONTA(CTL,CTL1)
        IMPLICIT REAL*8 (A-H,O-Z)
        CHARACTER NOM1*11,NOM2*11,RESP*1

10      CONTINUE
c       Las explicaciones estan de mas
        WRITE (*,*) 'ESCOJA:'
        WRITE (*,*) '1)Calcular N();2)Calcular un valor medio con su des
     +biacion'
        READ (*,*) IRESP
        WRITE (*,*) 'NOMBRE DEL FICHERO DE ENTRADA (con ext.):'
        READ (*,100) NOM1
100     FORMAT (A)
        WRITE (*,*) 'NUMERO DE ENTRADAS:'
        READ (*,*) N
        IF (IRESP.EQ.2) GO TO 20
        WRITE (*,*) 'NOMBRE DEL FICHERO DE SALIDA (con ext.):'
        READ (*,100) NOM2
        OPEN (1,FILE=NOM1)
        OPEN (2,FILE=NOM2)
        WRITE (*,*) 'QUE QUIERE?'
        WRITE (*,*) '1 para N(A),2 para N(D),3 para N(AM),4 para N(DM),5
     + para N(R),6 para N(Vr),7   para N(U),8 para N(V),9 para N(W),10 p
     +ara N(m),11 para N(B-V),12 para N(TS),13 para N(LC),14 para N(l),1
     +5 para N(b),16 para N(lM),17 para N(bM),18 para N(X), 19 para N(Y) 
     +,20 para N(Z)'
        READ (*,*) IRESP
        XMAX=0.0
        XMIN=0.0
        ICOUNT=0
        DO I=1,N
                READ (1,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                IF (IRESP.EQ.1) X=A
                IF (IRESP.EQ.2) X=D
                IF (IRESP.EQ.3) X=AM
                IF (IRESP.EQ.4) X=DM
                IF (IRESP.EQ.5) X=R
                IF (IRESP.EQ.6) X=IVR
                IF (IRESP.EQ.7) X=U
                IF (IRESP.EQ.8) X=V
                IF (IRESP.EQ.9) X=W
                IF (IRESP.EQ.10) X=IVM
                IF (IRESP.EQ.11) X=IC
                IF (IRESP.EQ.12) X=ITS1
                IF (IRESP.EQ.13) X=LC
                IF (IRESP.EQ.14) X=GL
                IF (IRESP.EQ.15) X=B
                IF (IRESP.EQ.16) X=GLM
                IF (IRESP.EQ.17) X=BM
                IF (IRESP.EQ.18) X=X
                IF (IRESP.EQ.19) X=Y
                IF (IRESP.EQ.20) X=Z
c       Determinamos el valor maximo y minimo del valor estudiado                
                IF (X.GT.XMAX) XMAX=X
                IF (X.LT.XMIN) XMIN=X
                ICOUNT=ICOUNT+1
                WRITE (*,400) 'TRABAJANDO...',ICOUNT
        END DO
        REWIND (1)
        WRITE (*,200) 'EL VALOR MAXIMO ES:',XMAX
        WRITE (*,200) 'EL VALOR MINIMO ES:',XMIN
200     FORMAT (1X,A19,1X,F10.4)
c       Pedimos los valores maximo y minimo a estudir 
        WRITE (*,*) 'VALOR MAXIMO A CONSIDERAR:'
        READ (*,*) XMAX
        WRITE (*,*) 'VALOR MINIMO A CONSIDERAR:'
        READ (*,*) XMIN
        WRITE (*,*) 'CUAL ES EL PASO DE INTEGRACION?'
        READ (*,*) PASO
c       Determinamos el numero de intervalos en los que dividir el valor        
        N1=((XMAX-XMIN)/PASO)+1
        REF1=XMIN
        ICOUNT=N1
c       Contamos y escrivimos los resultados en el fichero de salida       
        DO I=1,N1
                REF2=REF1+PASO 
                NUM=0.0
                DO J=1,N
                        READ (1,500) IHD1,IHD2,A,D,AM,DM,R,IVR,
     +                  U,V,W,IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                        IF (IRESP.EQ.1) X=A
                        IF (IRESP.EQ.2) X=D
                        IF (IRESP.EQ.3) X=AM
                        IF (IRESP.EQ.4) X=DM
                        IF (IRESP.EQ.5) X=R
                        IF (IRESP.EQ.6) X=IVR
                        IF (IRESP.EQ.7) X=U
                        IF (IRESP.EQ.8) X=V
                        IF (IRESP.EQ.9) X=W
                        IF (IRESP.EQ.10) X=IVM
                        IF (IRESP.EQ.11) X=IC
                        IF (IRESP.EQ.12) X=ITS1
                        IF (IRESP.EQ.13) X=LC
                        IF (IRESP.EQ.14) X=GL
                        IF (IRESP.EQ.15) X=B
                        IF (IRESP.EQ.16) X=GLM
                        IF (IRESP.EQ.17) X=BM
                        IF (IRESP.EQ.18) X=X
                        IF (IRESP.EQ.19) X=Y
                        IF (IRESP.EQ.20) X=Z
                        IF (X.GE.REF1.AND.X.LT.REF2) NUM=NUM+1
                END DO
                REF=(REF1+REF2)/2.0
                WRITE (2,300) REF,NUM
                REF1=REF2
                REWIND (1)
c       Ahora indicamos el numero de iteraciones que faltan para acabar               
                ICOUNT=ICOUNT-1
                WRITE (*,400) 'TRABAJANDO...',ICOUNT
        END DO
300     FORMAT (1X,F9.4,1X,I5)
400     FORMAT ('+',A13,1X,I5)                
500     FORMAT(1X,I3,1X,I3,1X,F11.7,1X,F10.6,1X,F8.4,1X,F8.4,1X,F9.4,
     +  1X,I5,1X,F9.4,1X,F9.4,1X,F9.4,1X,I4,1X,I4,1X,I2,1X,I3,1X,I2,1X,
     +  F10.6,1X,F10.6,1X,F9.4,1X,F8.4,1X,F10.4,1X,F10.4,1X,F10.4)
        CLOSE (1)
        CLOSE (2)
        GO TO 30
20      CONTINUE
c       Pedimos el valor a estudiar
        WRITE (*,*) 'QUE QUIERE ESTUDIAR?'
        WRITE (*,*) '1)r;2)Vr;3)U;4)V;5)W'
        READ (*,*) IRESP
        OPEN (1,FILE=NOM1)
        ICOUNT=0
c       Determinamos el valor maximo ,etc.
        XMAX=0.0
        XMIN=0.0
        DO I=1,N
                READ (1,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                IF (IRESP.EQ.1) X=R
                IF (IRESP.EQ.2) X=IVR
                IF (IRESP.EQ.3) X=U
                IF (IRESP.EQ.4) X=V
                IF (IRESP.EQ.5) X=W
                IF (X.GT.XMAX) XMAX=X
                IF (X.LT.XMIN) XMIN=X
                ICOUNT=ICOUNT+1
                WRITE (*,400) 'TRABAJANDO...',ICOUNT
        END DO
        REWIND (1)
        WRITE (*,200) 'EL VALOR MAXIMO ES:',XMAX
        WRITE (*,200) 'EL VALOR MINIMO ES:',XMIN
        WRITE (*,*) 'VALOR MAXIMO A CONSIDERAR:'
        READ (*,*) XMAX
        WRITE (*,*) 'VALOR MINIMO A CONSIDERAR:'
        READ (*,*) XMIN
        XMED=0.0
        ICOUNT=0
        N1=0
        DO I=1,N        
                READ (1,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                IF (IRESP.EQ.1) X=R
                IF (IRESP.EQ.2) X=IVR
                IF (IRESP.EQ.3) X=U
                IF (IRESP.EQ.4) X=V
                IF (IRESP.EQ.5) X=W
c       Descartamos los valores que no queremos
                IF (X.GT.XMAX.OR.X.LT.XMIN) GO TO 40
                XMED=XMED+X
                N1=N1+1
40      CONTINUE
                ICOUNT=ICOUNT+1
                WRITE (*,400) 'TRABAJANDO...',ICOUNT
        END DO
c       Calculamos la media
        XMED=XMED/N1
        REWIND (1)
c       Idem para la desviacion estandar        
        SIGMA=0.0
        DO I=1,N        
                READ (1,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                IF (IRESP.EQ.1) X=R
                IF (IRESP.EQ.2) X=IVR
                IF (IRESP.EQ.3) X=U
                IF (IRESP.EQ.4) X=V
                IF (IRESP.EQ.5) X=W
                IF (X.GT.XMAX.OR.X.LT.XMIN) GO TO 50
                SIGMA=SIGMA+(X-XMED)**2
50      CONTINUE
                ICOUNT=ICOUNT+1
                WRITE (*,400) 'TRABAJANDO...',ICOUNT
        END DO
        SIGMA=SQRT(SIGMA/(N1-1))
        WRITE (*,900) 'NUMERO DE ESTRELLAS CONSIDERADAS:',N1
	  WRITE (*,600) 'VALOR MEDIO:',XMED
        WRITE (*,700) 'DESVIACION ESTANDAR:',SIGMA
600     FORMAT (1X,A12,1X,F9.3)
700     FORMAT (1X,A20,1X,F9.3)
900	  FORMAT (1X,A33,1X,I5)
        CLOSE (1)
30      CONTINUE
        WRITE (*,*) 'QUIERE HACER OTRA ESTADISTICA?(S/N)'
        READ (*,100) RESP
        IF (RESP.EQ.'S'.OR.RESP.EQ.'s') GO TO 10
        CTL1=CTL
        RETURN
        END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

c       Subrutina que calcula la posicion del polo norte galactico
c       mediante el metodo de Newcomb
c       Utilizamos una subrutina para marcar claramente que es una 
c       parte independiente del programa (se podria hacer un programa
c       independiente de GALAC1)
        SUBROUTINE PNG(CTL,CTL1)
        IMPLICIT REAL*8 (A-H,O-Z)
        PARAMETER (PI=3.141592654)
        CHARACTER NOM1*11,RESP*1

40      CONTINUE
        WRITE (*,*) 'NOMBRE DEL FICHERO DE ENTRADA (con ext.):'
        READ (*,100) NOM1
        WRITE (*,*) 'NUMERO DE ENTRADAS:'
        READ (*,*) N
100     FORMAT (A)
        OPEN (1,FILE=NOM1)
c       Calculamos los coeficientes para el determinante de multi-
c       plicadores de Lagrange
        A11=0.0
        A12=0.0
        A13=0.0
        A22=0.0
        A23=0.0
        A33=0.0
        ICOUNT=0
        DO I=1,N
                READ (1,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                A=15.0*A*2.0*PI/360.0
                D=D*2.0*PI/360.0
                X=COS(A)*COS(D)
                Y=SIN(A)*COS(D)
                Z=SIN(D)
                A11=A11+X**2
                A12=A12+X*Y
                A13=A13+X*Z
                A22=A22+Y**2
                A23=A23+Y*Z
                A33=A33+Z**2
                ICOUNT=ICOUNT+1
                WRITE (*,400) 'TRABAJANDO...',ICOUNT
        END DO 
500     FORMAT(1X,I3,1X,I3,1X,F11.7,1X,F10.6,1X,F8.4,1X,F8.4,1X,F9.4,
     +  1X,I5,1X,F9.4,1X,F9.4,1X,F9.4,1X,I4,1X,I4,1X,I2,1X,I3,1X,I2,1X,
     +  F10.6,1X,F10.6,1X,F9.4,1X,F8.4,1X,F10.4,1X,F10.4,1X,F10.4)
400     FORMAT ('+',A13,1X,I5)       
        CLOSE (1)
        A1=-(A11+A22+A33)
        A2=A11*A22+A11*A33+A22*A33-A13**2-A23**2-A12**2
        A3=-(A11*A22*A33+2.0*A12*A23*A13-A13**2*A22-A23**2*A11-A12**2*
     +  A33)
c       Resolbemos el determinante (ecuacion cubica) y encontramos el
c       resultado correspondiente
        Q=(3.0*A2-A1**2)/9.0
        R=(-9.0*A1*A2+27.0*A3+2.0*A1**3)/27.0
        D=4.0*Q**3+R**2
        IF (D.GT.0.0) THEN
                S=(-R+SQRT(D))/2.0
                IF (S.LT.0.0) THEN
                        S=-S
                        S=-S**(1.0/3.0)
                ELSE
                        S=S**(1.0/3.0)
                END IF
                T=(-R-SQRT(D))/2.0
                IF (T.LT.0.0) THEN
                        T=-T
                        T=-T**(1.0/3.0)
                ELSE
                        T=T**(1.0/3.0)
                END IF
                X1=S+T-A1/3.0
                CALL EQ(A11,A12,A13,A22,A23,A33,X1,AG1,DG1)
                AG=(AG1*360.0)/(30.0*PI)
                DG=(DG1*360.0)/(2.0*PI)
                IAGH=INT(AG)
                IAGM=INT((AG-IAGH)*60.0)
                AGS=((AG-IAGH)*60.0-IAGM)*60.0
                IDGG=INT(DG)
                IDGM=INT((DG-IDGG)*60.0)
                DGS=((DG-IDGG)*60.0-IDGM)*60.0
                WRITE (*,200) 'Alfa=',IAGH,'h',IAGM,'m',AGS,'s'
                WRITE (*,201) 'Delta=',IDGG,'',IDGM,'''',DGS,'"'
        END IF
200     FORMAT (1X,A5,1X,I3,A1,1X,I3,A1,1X,F7.3,A1)
201     FORMAT (1X,A6,I3,A1,1X,I3,A1,1X,F7.3,A1)
        IF (D.LE.0.0) THEN
                THETA=ACOS(-R/(2.0*SQRT(-Q**3)))
                X1=2.0*SQRT(-Q)*COS(THETA/3.0)-A1/3.0
                X2=2.0*SQRT(-Q)*COS(THETA/3.0+2.094395102)-A1/3.0
                X3=2.0*SQRT(-Q)*COS(THETA/3.0+4.188790205)-A1/3.0
                CALL EQ (A11,A12,A13,A22,A23,A33,X1,AG1,DG1)      
                CALL EQ (A11,A12,A13,A22,A23,A33,X2,AG2,DG2) 
                CALL EQ (A11,A12,A13,A22,A23,A33,X3,AG3,DG3)
                CALL FUNCIO (A11,A12,A13,A22,A23,A33,X1,AG1,DG1,F1)
                CALL FUNCIO (A11,A12,A13,A22,A23,A33,X2,AG2,DG2,F2)
                CALL FUNCIO (A11,A12,A13,A22,A23,A33,X3,AG3,DG3,F3)
                IF (F1.LE.F2.AND.F1.LE.F3) THEN
                        AG=(AG1*360.0)/(30.0*PI)
                        DG=(DG1*360.0)/(2.0*PI)
                        IAGH=INT(AG)
                        IAGM=INT((AG-IAGH)*60.0)
                        AGS=((AG-IAGH)*60.0-IAGM)*60.0
                        IDGG=INT(DG)
                        IDGM=INT((DG-IDGG)*60.0)
                        DGS=((DG-IDGG)*60.0-IDGM)*60.0
                        WRITE (*,200)'Alfa=',IAGH,'h',IAGM,'m',AGS,'s'
                        WRITE (*,201)'Delta=',IDGG,'',IDGM,'''',DGS,'"'
                END IF
                IF (F2.LE.F1.AND.F2.LE.F3) THEN
                        AG=(AG2*360.0)/(30.0*PI)
                        DG=(DG2*360.0)/(2.0*PI)
                        IAGH=INT(AG)
                        IAGM=INT((AG-IAGH)*60.0)
                        AGS=((AG-IAGH)*60.0-IAGM)*60.0
                        IDGG=INT(DG)
                        IDGM=INT((DG-IDGG)*60.0)
                        DGS=((DG-IDGG)*60.0-IDGM)*60.0
                        WRITE (*,200)'Alfa=',IAGH,'h',IAGM,'m',AGS,'s'
                        WRITE (*,201)'Delta=',IDGG,'',IDGM,'''',DGS,'"'
                END IF
                IF (F3.LE.F1.AND.F3.LE.F2) THEN
                        AG=(AG3*360.0)/(30.0*PI)
                        DG=(DG3*360.0)/(2.0*PI)
                        IAGH=INT(AG)
                        IAGM=INT((AG-IAGH)*60.0)
                        AGS=((AG-IAGH)*60.0-IAGM)*60.0
                        IDGG=INT(DG)
                        IDGM=INT((DG-IDGG)*60.0)
                        DGS=((DG-IDGG)*60.0-IDGM)*60.0
                        WRITE (*,200)'Alfa=',IAGH,'h',IAGM,'m',AGS,'s'
                        WRITE (*,201)'Delta=',IDGG,'',IDGM,'''',DGS,'"'
                END IF
        END IF
        WRITE (*,*) 'QUIERE VOLVER A HACER LOS CALCULOS (S/N)?'
        READ (*,100) RESP
        IF (RESP.EQ.'S'.OR.RESP.EQ.'s') GO TO 40
        CTL1=CTL
        RETURN
        END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

c       Subrutina de calculo de AG y DG para un valor del multiplica-
c       dor de Lagrange X
        SUBROUTINE EQ(A11,A12,A13,A22,A23,A33,X,AG,DG)
        IMPLICIT REAL*8 (A-H,O-Z)
        PARAMETER (PI=3.141592654)
        
        AG=ATAN(((A23*(A11-X))/A13-A12)/(A22-X-(A12*A23)/A13))
        DG=ATAN(-(A12*SIN(AG)+(A11-X)*COS(AG))/A13)
        F=A13*COS(DG)*COS(AG)+A23*COS(DG)*SIN(AG)+(A33-X)*SIN(DG)
        IF (F.GT.1D-20.AND.F.LT.-1D-20) THEN
                AG=AG+PI
                IF (AG.GT.2.0*PI) AG=AG-2.0*PI
                DG=ATAN(-(A12*SIN(AG)+(A11-X)*COS(AG))/A13)  
        END IF
        IF (AG.LT.0.0) AG=AG+2.0*PI
        IF (DG.LT.0.0) THEN
                DG=-DG
                AG=AG+PI
        END IF
        IF (AG.GT.2.0*PI) AG=AG-2.0*PI
        RETURN
        END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

c       Subrutina para el calculo de la funcion a minimizar en el me-
c       todo de Newcomb
        SUBROUTINE FUNCIO(A11,A12,A13,A22,A23,A33,X,AG,DG,F)
        IMPLICIT REAL*8 (A-H,O-Z)

        RL=COS(DG)*COS(AG)
        RM=COS(DG)*SIN(AG)
        RN=SIN(DG)
        F=(A11-X)*RL**2+(A22-X)*RM**2+(A33-X)*RN**2+2.0*(A12*RL*RM+
     +  A13*RL*RN+A23*RM*RN)+X
        RETURN
        END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

c       Subrutina para el Estudio Cinematico del Entorno Solar
c       Como en las anteriores,usamos una subrutina para marcar que es
c       un bloque diferenciado del programa
        SUBROUTINE ECES(CTL,CTL1)
        IMPLICIT REAL*8 (A-H,O-Z)
        CHARACTER NOM1*11,RESP*1

10      CONTINUE
c       Pedimos el nombre del fichero a estudiar y los datos pertinentes        
        WRITE (*,*) 'NOMBRE DEL FICHERO DE ENTRADA (con ext.):'
        READ (*,100) NOM1
100     FORMAT (A)
        WRITE (*,*) 'NUMERO DE ENTRADAS:'
        READ (*,*) N
        OPEN (1,FILE=NOM1)
c       Calculamos los valores maximos y minimos presentes en el fichero        
c       y ofrecemos la posibilidad de restringir el intervalo de estudio
        UMAX=0.0
        UMIN=0.0
        VMAX=0.0
        VMIN=0.0
        WMAX=0.0
        WMIN=0.0
        ICOUNT=0
        DO I=1,N
                READ (1,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
c       Determinamos el valor maximo y minimo del valor estudiado                
                IF (U.GT.UMAX) UMAX=U
                IF (U.LT.UMIN) UMIN=U
                IF (V.GT.VMAX) VMAX=V
                IF (V.LT.VMIN) VMIN=V
                IF (W.GT.WMAX) WMAX=W
                IF (W.LT.WMIN) WMIN=W
                ICOUNT=ICOUNT+1
                WRITE (*,400) 'TRABAJANDO...',ICOUNT
        END DO
        WRITE (*,300) 'EL VALOR MAXIMO DE U ES:',UMAX
        WRITE (*,300) 'EL VALOR MINIMO DE U ES:',UMIN
        WRITE (*,*) 'VALOR MAXIMO Y MINIMO A CONSIDERAR (max,min):'
        READ (*,*) UMAX,UMIN
        WRITE (*,300) 'EL VALOR MAXIMO DE V ES:',VMAX
        WRITE (*,300) 'EL VALOR MINIMO DE V ES:',VMIN
        WRITE (*,*) 'VALOR MAXIMO Y MINIMO A CONSIDERAR (max,min):'
        READ (*,*) VMAX,VMIN
        WRITE (*,300) 'EL VALOR MAXIMO DE W ES:',WMAX
        WRITE (*,300) 'EL VALOR MINIMO DE W ES:',WMIN
        WRITE (*,*) 'VALOR MAXIMO Y MINIMO A CONSIDERAR (max,min):'
        READ (*,*) WMAX,WMIN                                     
300     FORMAT (1X,A24,1X,F10.4)
        REWIND (1)                
c       Calculamos las medias        
        U0=0.0
        V0=0.0
        W0=0.0
        ICOUNT=2*N
        N1=0
        DO I=1,N
                READ (1,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                IF (U.GT.UMAX.OR.U.LT.UMIN.OR.V.GT.VMAX.OR.V.LT.VMIN.
     +          OR.W.GT.WMAX.OR.W.LT.WMIN) GO TO 40
                U0=U0+U
                V0=V0+V
                W0=W0+W
                N1=N1+1
40      CONTINUE                   
                ICOUNT=ICOUNT-1
                WRITE (*,400) 'TRABAJANDO...',ICOUNT
        END DO
400     FORMAT ('+',A13,1X,I5)
        U0=U0/N1
        V0=V0/N1
        W0=W0/N1
500     FORMAT(1X,I3,1X,I3,1X,F11.7,1X,F10.6,1X,F8.4,1X,F8.4,1X,F9.4,
     +  1X,I5,1X,F9.4,1X,F9.4,1X,F9.4,1X,I4,1X,I4,1X,I2,1X,I3,1X,I2,1X,
     +  F10.6,1X,F10.6,1X,F9.4,1X,F8.4,1X,F10.4,1X,F10.4,1X,F10.4)
        REWIND (1)
c       Y los valores deseados        
        VMU200=0.0
        VMU110=0.0
        VMU101=0.0
        VMU020=0.0
        VMU011=0.0
        VMU002=0.0
        VMU400=0.0
        VMU220=0.0
        VMU202=0.0
        VMU022=0.0
        VMU040=0.0
        VMU004=0.0
        DO I=1,N
                READ (1,500) IHD1,IHD2,A,D,AM,DM,R,IVR,U,V,W,
     +          IVM,IC,ITS1,ITS2,LC,GL,B,GLM,BM,X,Y,Z
                IF (U.GT.UMAX.OR.U.LT.UMIN.OR.V.GT.VMAX.OR.V.LT.VMIN.
     +          OR.W.GT.WMAX.OR.W.LT.WMIN) GO TO 50
                VMU200=VMU200+(U-U0)**2
                VMU110=VMU110+(U-U0)*(V-V0)
                VMU101=VMU101+(U-U0)*(W-W0)
                VMU020=VMU020+(V-V0)**2
                VMU011=VMU011+(V-V0)*(W-W0)
                VMU002=VMU002+(W-W0)**2
                VMU400=VMU400+(U-U0)**4
                VMU220=VMU220+(U-U0)**2*(V-V0)**2
                VMU202=VMU202+(U-U0)**2*(W-W0)**2
                VMU022=VMU022+(V-V0)**2*(W-W0)**2
                VMU040=VMU040+(V-V0)**4
                VMU004=VMU004+(W-W0)**4
50      CONTINUE                
                ICOUNT=ICOUNT-1
                WRITE (*,400) 'TRABAJANDO...',ICOUNT
        END DO
        VMU200=VMU200/N1
        VMU110=VMU110/N1
        VMU101=VMU101/N1
        VMU020=VMU020/N1
        VMU011=VMU011/N1
        VMU002=VMU002/N1
        VMU400=VMU400/N1
        VMU220=VMU220/N1
        VMU202=VMU202/N1
        VMU022=VMU022/N1
        VMU040=VMU040/N1
        VMU004=VMU004/N1
        WRITE (*,900) 'NUMERO DE ESTRELLAS CONSIDERADAS:',N1
	  VMU11=SQRT(VMU400/N1+VMU200**2.0*(2.0/(N1-1)-3.0/N1))
        VMU200=SQRT(VMU200)
        VMU11=VMU11/(2.0*VMU200)
        WRITE (*,200) 'sqrt<u**2>=',VMU200,'+/-',VMU11
        VMU12=SQRT(VMU220/N1+VMU110**2*(1.0/(N1-1)-2.0/N1)+VMU200**2*
     +  VMU020*(1.0/(N1-1)-1.0/N1))
        WRITE (*,201) '<uv>=',VMU110,'+/-',VMU12
        VMU13=SQRT(VMU202/N1+VMU101**2*(1.0/(N1-1)-2.0/N1)+VMU200**2*
     +  VMU002*(1.0/(N1-1)-1.0/N1))
        WRITE (*,201) '<uw>=',VMU101,'+/-',VMU13
        VMU22=SQRT(VMU040/N1+VMU020**2.0*(2.0/(N1-1)-3.0/N1))
        VMU020=SQRT(VMU020)
        VMU22=VMU22/(2.0*VMU020)
        WRITE (*,200) 'sqrt<v**2>=',VMU020,'+/-',VMU22
        VMU23=SQRT(VMU022/N1+VMU011**2*(1.0/(N1-1)-2.0/N1)+VMU020**2*
     +  VMU002*(1.0/(N1-1)-1.0/N1))
        WRITE (*,201) '<vw>=',VMU011,'+/-',VMU23
        VMU33=SQRT(VMU004/N1+VMU002**2.0*(2.0/(N1-1)-3.0/N1))
        VMU002=SQRT(VMU002)
        VMU33=VMU33/(2.0*VMU002)
        WRITE (*,200) 'sqrt<w**2>=',VMU002,'+/-',VMU33
200     FORMAT (1X,A11,D12.4,1X,A3,D10.2)
201     FORMAT (1X,A5,6X,D12.4,1X,A3,D10.2)
900	  FORMAT (1X,A33,1X,I5)
c       Calculamos la desviacion del vertice              
        X=2.0*(VMU110/(VMU200**2-VMU020**2))
        GLV=0.5*ATAN(X)
        EGLV=ABS(X)*SQRT((VMU12/VMU110)**2+(SQRT((VMU11*2.0*VMU200)
     +  **2+(VMU22*2.0*VMU020)**2)/(VMU200**2-VMU020**2))**2)
        EGLV=ABS(0.5/(1.0+X**2))*EGLV
        GLV=GLV*360.0/(2.0*3.141592654)
        EGLV=EGLV*360.0/(2.0*3.141592654)
        WRITE (*,202) 'lv=',GLV,'+/-',EGLV
202     FORMAT (1X,A3,F7.2,1X,A3,F7.2)
c       Ofrecemos la posibilidad de calcular los semiejes del elipsoide         
c       de velocidades (de momento sin error)
        WRITE (*,*) 'QUIERE LOS VALORES DE LOS SEMIEJES DEL ELIPSOIDE D 
     +E VELOCIDADES?(S/N)'
        READ (*,100) RESP
        IF (RESP.EQ.'S'.OR.RESP.EQ.'s') THEN
                A1=-(VMU200**2+VMU020**2+VMU002**2)
                A2=VMU200**2*VMU020**2+VMU200**2*VMU002**2+VMU020**2*
     +          VMU002**2-VMU101**2-VMU011**2-VMU110**2
                A3=-(VMU200**2*VMU020**2*VMU002**2+2.0*VMU110*VMU011*
     +          VMU101-VMU101**2*VMU020**2-VMU011**2*VMU200**2-VMU110**
     +          2*VMU002**2)
c       Resolbemos el determinante (ecuacion cubica) y encontramos el
c       resultado correspondiente
                Q=(3.0*A2-A1**2)/9.0
                R=(-9.0*A1*A2+27.0*A3+2.0*A1**3)/27.0
                D=4.0*Q**3+R**2
                IF (D.GT.0.0) THEN
                        S=(-R+SQRT(D))/2.0
                        IF (S.LT.0.0) THEN
                                S=-S
                                S=-S**(1.0/3.0)
                        ELSE
                                S=S**(1.0/3.0)
                        END IF
                        T=(-R-SQRT(D))/2.0
                        IF (T.LT.0.0) THEN
                                T=-T
                                T=-T**(1.0/3.0)
                        ELSE
                                T=T**(1.0/3.0)
                        END IF
                        X1=S+T-A1/3.0
                        WRITE (*,*) X1
                END IF
                IF (D.LE.0.0) THEN
                        THETA=ACOS(-R/(2.0*SQRT(-Q**3)))
                        X1=2.0*SQRT(-Q)*COS(THETA/3.0)-A1/3.0
                        X2=2.0*SQRT(-Q)*COS(THETA/3.0+2.094395102)-
     +                  A1/3.0
                        X3=2.0*SQRT(-Q)*COS(THETA/3.0+4.188790205)-
     +                  A1/3.0
                        WRITE (*,*) X1,X2,X3
                END IF
        END IF
        CLOSE (1)
c       Como siempre...        
        WRITE (*,*) 'QUIERE HACER OTRO ESTUDIO?(S/N)'
        READ (*,100) RESP
        IF (RESP.EQ.'S'.OR.RESP.EQ.'s') GO TO 10
        CTL1=CTL
        RETURN
        END
      
      
DATOS TECNICOS

        NOMBRE:GALAC
        VERSION:1.02
        ESPACIO NECESARIO EN DISCO:minimo 82K
        CPU:funciona en un 286,por lo tanto no deveria presentar
            problemas en ordenadores mas potentes
        PROGRAMADOR:
                -Jerkos
        E-MAIL DE CONTACTO:
                (seele01_es@yahoo.com)

