C PROGRAM VES C C Para compilar: C f77 -o ves ves.f -L$PGPLOT_DIR -lpgplot -lX11 C---------------------------------------------------------------------- C IMPLICIT NONE C INTEGER I,J,N,NN INTEGER NWIN ! Numero de ventanas activas en el plot. PARAMETER (NWIN=19) INTEGER NOBS ! Cuantos observatorios PARAMETER (NOBS=10) INTEGER NNOBS ! Observatorio elegido INTEGER PGOPEN INTEGER INDATE(3),INTIME(3) INTEGER IDAY,IMONTH,IYEAR INTEGER IHOUR,IMIN,ISEC INTEGER NSTAR,NSTAROLD ! Numero de objetos trazados. INTEGER TESTN INTEGER DEV(2) C REAL WIN(4,NWIN) REAL ANGLE ! Valores para PGPTXT REAL XX,YY ! Posicion del cursor REAL X,Y REAL X1,X2,Y1,Y2 REAL*8 TS0 ! Tiempo sidereo en Greenwich a 0h TU REAL*8 TS1 ! TS0 al dia siguiente. REAL*8 FJ0 ! Fecha Juliana a O TU REAL*8 FJ1 ! FJ0+1. Fecha Juliana de "maniana" REAL*8 FJ ! Fecha Juliana para dia IDAY y hora RHOUR REAL*8 ARSUN,DECSUN ! AR Y DEC del Sol en la fech de observ.(FJ) REAL*8 SEMIA,SEMIB ! semiejes del elipsoide WGS-84 REAL*8 REARTH ! radio terrestre para el observatorio REAL*8 ANGHOA ! angulo depresion horizonte debido a la altura REAL*8 PI PARAMETER(PI=3.141592654D0) REAL*8 DISZEN ! distancia zenital del sol REAL*8 LAT,LONG,ALTOBS ! Latitud, longitud y altura del obser.(DD.dddd) REAL*8 HORTO,HOCASO ! angulos horarios de ortos y ocasos REAL*8 AORTO,AOCASO ! azimuts de ortos y ocasos REAL*8 TUORTO1,TUORTO2 ! Tiempo universal de orto y ocaso REAL*8 TUOCASO1,TUOCASO2 ! utilizados para la iteracion REAL*8 TSL,TSL0 ! Tiempo sidereo local, y a 0h TU REAL*8 TSLORTO,TSLOCASO ! TSL para orto y ocaso REAL*8 TUCREPM1,TUCREPM2 ! Tiempo universal de crepusculos matutino REAL*8 TUCREPV1,TUCREPV2 ! y vespertino utilizados para la iteracion REAL*8 TSLCREPM,TSLCREPV ! TSL inicio crep. matutino y fin vespertino REAL XMIN,XMAX ! Limites de la grafica en segundos. REAL YMIN,YMAX REAL YYY REAL XCREPV,XCREPM ! Crepusculos en segundos. REAL*8 ARII(100),DECII(100) ! AR y DEC del equinocio dado de entrada REAL*8 ARFF(100),DECFF(100) ! AR y DEC de la fecha elegida REAL*8 AZI(100,100),ALT(100,100)! Azimut y altura " " REAL*8 TII,TFF ! Equinocio de entrada y fecha (salida) REAL*8 MONTH REAL TUPH(100),TUPS(100,100)!T.U. para plot estrella, en horas y seg REAL TUPSU(100),ALTU(100) ! T.U. y Alt para la ultima estrella. REAL AZIU(100) ! Azimut para la ultima estrella. REAL SECZ(10) ! Valores de sec Z, masa aire que utiliz. REAL XAZI(100,100),YAZI(100,100) ! Valores de X e Y para el plot REAL RADIO(100) ! azimut-altura REAL XMAZI(48),YMAZI(48) ! Lo mismo para la luna REAL*8 NIGHTL ! Duracion de la noche REAL*8 NIGHTLF ! " " " " sin crepusculos REAL*8 ALTCUL,TUCUL ! Altura y T.U. de la culminacion de un objeto REAL*8 TVIS(100) ! Tiempo que la estrella esta por REAL*8 TVISINI(100),TVISFIN(100) ! encima horizonte en segundos REAL*8 TUORSTAR,TUOCSTAR ! T.U. del orto y el ocaso de una estrella. REAL*8 DAYS,MALFA,MDELTA REAL*8 MOONALT(48),MOONAZI(48) ! Altura y Azimut de Luna, cada hora REAL*8 FASE(48) ! Fase de la Luna C CHARACTER*1 CH CHARACTER*2 CH2 CHARACTER*11 CHDATE CHARACTER*13 LABEL(NWIN),VALUE(NWIN) CHARACTER*11 STARNAME(100),STARRA(100),STARDEC(100) CHARACTER*9 STARRAF(100),STARDECF(100) CHARACTER*6 EQUINOX(100) CHARACTER*10 OBSNAME(10),OBSLONG(10) CHARACTER*9 OBSLAT(10) CHARACTER*5 OBSHEI(10) CHARACTER*11 TEXTBCK CHARACTER*1 ESTADO ! Objeto circumpolar, invisible o con orto/ocaso CHARACTER*7 CHORTOTU,CHOCASOTU ! Hora en formato XXh XXm CHARACTER*7 CHORTOTSL,CHOCASOTSL ! " " " CHARACTER*22 COPYR CHARACTER*3 CHSECZ(10) CHARACTER*11 INFILE ! Listado de estrellas,N, Nombre, AR, DEC, Equin CHARACTER*41 LINE ! con formato I2,A11,A9,A9,A6 CHARACTER*6 CHORA CHARACTER*10 CFJ0 CHARACTER*7 CNIGHTL,CNIGHTLF,CTFF CHARACTER*1 DIRCUL CHARACTER*5 CALTCUL CHARACTER*7 CTUCUL,CTVIS CHARACTER*4 CHF ! Fase de la Luna C LOGICAL LELEV,LSECZ,LUT,LAZI LOGICAL LUNPLOT ! TRUE -> pgsci(0) -> borra el plot anterior LOGICAL LNSTAR ! TRUE -> Hay alguna estrella en el plot actual. LOGICAL LMOON(48) ! TRUE -> La Luna es visible LOGICAL VESIN ! TRUE -> Hay un fichero "vesinput" en ./ que se usara de input C COMMON/BLKOBS/LAT,LONG COMMON/BLKOOC/HORTO,HOCASO,AORTO,AOCASO c COMMON/BLKDATE/FJ,IYEAR,IMONTH,IDAY,IHOUR,IMIN,ISEC C C---------------------------------------------------------------------- C Definimos las ventanas C DATA WIN /0.02,0.25,0.86,0.90, ! 1 Star Name > 0.02,0.25,0.76,0.80, ! 2 Ascension recta > 0.02,0.25,0.66,0.70, ! 3 Declinacion > 0.02,0.13,0.56,0.60, ! 4 Equinocio > 0.02,0.19,0.20,0.24, ! 5 Draw > 0.02,0.19,0.10,0.14, ! 6 Clear > 1.73,1.99,0.86,0.90, ! 7 Nombre del Observatorio > 1.73,1.99,0.76,0.80, ! 8 Longitud del Obs. > 1.73,1.99,0.66,0.70, ! 9 Latitud del Obs. > 1.73,1.90,0.56,0.60, ! 10 Altura del Obs. > 1.73,2.00,0.46,0.50, ! 11 Fecha > 1.73,1.97,0.26,0.30, ! 12 Crea listado estrellas, coord... > 1.73,1.97,0.18,0.22, ! 13 Crea fichero postscript > 1.73,1.90,0.10,0.14, ! 14 Exit > 0.30,0.36,0.50,0.60, ! 15 Altura de la estrella > 0.30,0.36,0.30,0.40, ! 16 Secante de Z, Masa de aire > 0.75,0.90,0.10,0.14, ! 17 Tiempo Universal > 1.18,1.35,0.10,0.14, ! 18 Azimut > 0.60,1.50,0.00,0.06/ ! 19 I/O y ayuda. C DATA LABEL /'Name','R.A.','Dec.','Equinox',' ',' ', > 'Observatory','Longitude','Latitude','Height','Date', > ' ',' ',' ',' ',' ',' ',' ',' '/ C DATA VALUE /' -----',' 20 00 00','+20 00 00','2000', > ' DRAW',' CLEAR', > 'Calar Alto','- 02 32 46','+37 13 25','2168m',' ', > 'file: ves.lis','file: ves.ps',' EXIT', > 'Elev.','sec z',' U.T.',' Azimut',' ' / C DATA OBSNAME/'Calar Alto',' La Palma ',' La Silla ',' Madrid ', > ' Tenerife ',' Lick ',' UKIRT ',' KPNO ', > ' McDonald ','S.Pedro M.'/ C DATA OBSLONG/'- 02 32 46','- 17 52 47','- 70 44 00','- 03 41 15', > '- 16 29 45','-123 38 12','-155 28 24','-111 36 00', > '-104 01 18','-115 27 49'/ C DATA OBSLAT /'+37 13 25','+28 45 40','-29 15 00','+40 24 30', > '+28 17 32','+37 20 36','+19 49 32','+31 57 48', > '+30 40 18','+31 02 39'/ C DATA OBSHEI /'2168m','2334m','2400m',' 656m', > '2400m','1283m','4194m','2120m', > '2075m','2800m'/ C DATA SECZ /1.0,1.1,1.2,1.3,1.4,1.5,2.0,3.0,4.0,5.0/ DATA CHSECZ /'1.0','1.1','1.2','1.3','1.4','1.5', > '2.0','3.0','4.0','5.0'/ C C---------------------------------------------------------------------- C NNOBS=1 DO I=1,48 LMOON(I)=.FALSE. END DO C DEV(1)=PGOPEN('/XWINDOW') CALL PGSLCT(DEV(1)) CALL PGPAGE ! para que no se ponga negra la pantalla CALL PGSVP(0.0,1.0,0.0,1.0) CALL PGSWIN(0.0,2.0,0.0,1.0) CALL PGQWIN(X1,X2,Y1,Y2) CALL PGSCF(2) CALL PGSFS(1) C C---------------------------------------------------------------------- C Le pedimos la fecha a la maquina C CALL IDATE(INDATE) IDAY=INDATE(1) IMONTH=INDATE(2) IYEAR=INDATE(3) WRITE(CHDATE,"(I2,'-',I2,'-',I4)") IDAY,IMONTH,IYEAR WRITE(VALUE(11),'(A11)')CHDATE CALL ITIME(INTIME) IHOUR=INTIME(1) IMIN=INTIME(2) ISEC=INTIME(3) WRITE(CHORA,"(3(I2))")IHOUR,IMIN,ISEC C C---------------------------------------------------------------------- C Dibujamos las ventanitas. C DO I=1,NWIN CALL PGSCI(1) CALL PGRECT(WIN(1,I),WIN(2,I),WIN(3,I),WIN(4,I)) X=WIN(1,I) Y=WIN(4,I)+0.01 CALL PGSCH(0.8) CALL PGPTXT(X,Y,0.,0.0,LABEL(I)) CALL PGSCI(4) CALL PGSCH(1.0) IF(I.EQ.15.OR.I.EQ.16) THEN ANGLE=90. Y=WIN(3,I)+0.02 X=(WIN(1,I)+WIN(2,I)+0.02)/2. ELSE ANGLE=0. Y=WIN(3,I)+0.01 END IF CALL PGPTXT(X,Y,ANGLE,0.0,VALUE(I)) END DO C La ventanita de I/O y ayuda tendra un fondo gris CALL PGSCR(16,0.5,0.5,0.5) CALL IOWIN1(' Click center mouse button for HELP ') C C Al arrancar estaran activadas las opciones de Elevacion en eje Y C y Tiempo Universal en eje X. CALL PGSCI(4) CALL PGSFS(2) CALL PGSLW(12) CALL PGRECT(WIN(1,15),WIN(2,15),WIN(3,15),WIN(4,15)) CALL PGRECT(WIN(1,17),WIN(2,17),WIN(3,17),WIN(4,17)) CALL PGSLW(1) CALL PGQWIN(X1,X2,Y1,Y2) LELEV=.TRUE. LSECZ=.FALSE. LUT=.TRUE. LAZI=.FALSE. LNSTAR=.FALSE. ! Al empezar no hay ninguna estrella en el plot C C El observatorio inicial es Calar Alto, pasamos sus coordenadas C a valor decimal. CALL CH2DEC(VALUE(8),LONG,1) LONG=LONG/15.D0 ! Pasamos de grados a horas. CALL CH2DEC(VALUE(9),LAT,2) READ(VALUE(10)(1:4),*) ALTOBS C WRITE(COPYR,'(A22)')'\(0274) pedraz@caha.es' CALL PGSCI(1) CALL PGSCH(0.55) CALL PGPTXT(1.99,0.01,0.,1.,COPYR) CALL PGSCH(1.0) C C---------------------------------------------------------------------- C Calculamos los momentos de anochecer, amanecer y crepusculos para C delimitar la grafica. C 50 LUNPLOT=.FALSE. C C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C Calculamos la Fecha Juliana (FJ) y TS0 a las 0h TU en Greenwich. C A partir de estas se calcula la posicion del sol. Iteramos para C obtener la posicion del sol en el ocaso. Como valor inicial C tomamos la posicion a las 0h TU CALL FJ_TS0(IYEAR,IMONTH,IDAY,0.D0,FJ,TS0) FJ0=FJ TS1=TS0+24.*2.737909E-3 IF(TS1.GT.24.D0) TS1=TS1-24.D0 CALL POSSOL(FJ,ARSUN,DECSUN) C C Calculamos la distancia zenital en el ocaso teniendo en cuenta la C refraccion (34'), semidiametro del Sol (16') y la depresion del C horizonte(elipsoide WGS-84) dada la altura del observatorio. C SEMIA=6378.137D0 SEMIB=6356.752D0 REARTH=1.D0/ > SQRT(COS(LAT*PI/180.D0)*COS(LAT*PI/180.D0)/SEMIA/SEMIA+ > SIN(LAT*PI/180.D0)*SIN(LAT*PI/180.D0)/SEMIB/SEMIB) ANGHOA=ACOS(REARTH/(REARTH+ALTOBS/1000.D0)) ANGHOA=ANGHOA*180.D0/PI DISZEN=90.D0+34.D0/60.D0+16.D0/60.D0+ANGHOA C C- - - - - - - - - - - OCASO - - - - - - - - - - - - - - - - - - - C Con la declinacion y la distancia zenital calculamos angulo C horario y azimut del sol para el ocaso CALL ORCASO(DECSUN,DISZEN,ESTADO) IF(ESTADO.EQ.'C')THEN CALL IOWIN1(' Circumpolar SUN -- NO orto-ocaso ') GOTO 70 ELSE IF(ESTADO.EQ.'I')THEN CALL IOWIN1(' Invisible SUN -- NO orto-ocaso ') GOTO 70 END IF C Iterando, calculamos el ocaso que se produce despues de las 0h TU C del dia de la observacion. CALL TIEMPOU(ARSUN,HOCASO,LONG,TS0,TSL,TUOCASO1) 52 CALL FJ_TS0(IYEAR,IMONTH,IDAY,TUOCASO1,FJ,TS0) CALL POSSOL(FJ,ARSUN,DECSUN) CALL ORCASO(DECSUN,DISZEN,ESTADO) CALL TIEMPOU(ARSUN,HOCASO,LONG,TS0,TSL,TUOCASO2) TSLOCASO=TSL IF(ABS(TUOCASO1-TUOCASO2).GT.0.0003D0) THEN TUOCASO1=TUOCASO2 GOTO 52 END IF C C - - - - - - - Crepusculo ASTRONOMICO vespertino - - - - - - - - - C Como valor inicial tomamos el ocaso. TUCREPV1=TUOCASO2 C 53 CALL FJ_TS0(IYEAR,IMONTH,IDAY,TUCREPV1,FJ,TS0) CALL POSSOL(FJ,ARSUN,DECSUN) CALL ORCASO(DECSUN,108.0+ANGHOA,ESTADO) IF(ESTADO.EQ.'C') THEN CALL IOWIN1(' Twilight all night long ') GOTO 70 END IF CALL TIEMPOU(ARSUN,HOCASO,LONG,TS0,TSL,TUCREPV2) TSLCREPV=TSL IF(TUCREPV2.LT.TUCREPV1) TUCREPV2=TUCREPV2+24.0 IF(ABS(TUCREPV2-TUCREPV1).GT.0.0003D0) THEN TUCREPV1=TUCREPV2 GOTO 53 END IF C C - - - - - - - - - - ORTO - - - - - - - - - - - - - - - - - - - - C Calculamos el orto del dia siguiente. Partimos de que a las 0h TU C la fecha juliana sera un dia mas FJ1=FJ0+1. FJ1=FJ0+1.D0 IDAY=IDAY+1 54 CALL POSSOL(FJ1,ARSUN,DECSUN) CALL ORCASO(DECSUN,DISZEN,ESTADO) CALL TIEMPOU(ARSUN,HORTO,LONG,TS0,TSL,TUORTO1) C Iteramos 55 CALL FJ_TS0(IYEAR,IMONTH,IDAY,TUORTO1,FJ,TS0) CALL POSSOL(FJ,ARSUN,DECSUN) CALL ORCASO(DECSUN,DISZEN,ESTADO) CALL TIEMPOU(ARSUN,HORTO,LONG,TS0,TSL,TUORTO2) TSLORTO=TSL IF(ABS(TUORTO1-TUORTO2).GT.0.0003D0) THEN TUORTO1=TUORTO2 GOTO 55 END IF C C - - - - - - - Crepusculo ASTRONOMICO matutino - - - - - - - - - C Como valor inicial tomamos el orto. TUCREPM1=TUORTO2 C 56 CALL FJ_TS0(IYEAR,IMONTH,IDAY,TUCREPM1,FJ,TS0) CALL POSSOL(FJ,ARSUN,DECSUN) CALL ORCASO(DECSUN,108.0+ANGHOA,ESTADO) CALL TIEMPOU(ARSUN,HORTO,LONG,TS0,TSL,TUCREPM2) TSLCREPM=TSL IF(TUCREPM2.GT.TUCREPM1) TUCREPM2=TUCREPM2-24.0 IF(ABS(TUCREPM2-TUCREPM1).GT.0.0003D0) THEN TUCREPM1=TUCREPM2 GOTO 56 END IF IF(FJ1.EQ.FJ0) GOTO 57 IDAY=IDAY-1 IF(TUORTO2.GT.TUOCASO2) THEN ! Puede ser que debido a su longitud FJ1=FJ0 ! en un observatorio anochezca y GOTO 54 ! amanezca en el mismo dia juliano. END IF 57 CONTINUE C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C Los limites de la grafica seran el comienzo de la noche C ocaso=TUOCASO2 y el final orto=TUORTO2. El final del crepusculo C por la tarde es TUCREPV2 y el inicio del crepusculo por la C maniana TUCREPM2. C XMIN=TUOCASO2*3600. XCREPV=TUCREPV2*3600. IF(XCREPV.LT.XMIN) XCREPV=XCREPV+(24.*3600.) IF(FJ1.EQ.FJ0) THEN XMAX=TUORTO2*3600. XCREPM=TUCREPM2*3600. ELSE XMAX=(TUORTO2+24.)*3600. XCREPM=(TUCREPM2+24.)*3600. END IF YMIN=0. YMAX=90. 60 CALL PGSVP(0.22,0.85,0.22,0.90) !Definimos el area del plot que C !utilizaremos para la grafica CALL PGSWIN(XMIN,XMAX,YMIN,YMAX) C IF(LUT) THEN IF(LELEV) THEN CALL PGSCI(0) ! Borramos CALL PGSLS(1) CALL PGTBOX('ABCNTSZYXHO',0.0,0,'ABC',0.0,0) CALL PGSLS(4) DO I=1,10 Y=ACOS(1./SECZ(I)) ! Lineas horizontales para marcar las Y=Y*180./PI ! sec Z (masas de aire) Y=90.-Y CALL PGMOVE(XMIN,Y) CALL PGDRAW(XMAX,Y) CALL PGPTXT(XMIN-300.,Y-1.,0.0,1.0,CHSECZ(I)) END DO CALL PGSCI(1) ! Pintamos IF(LUNPLOT) CALL PGSCI(0) CALL PGSLS(1) CALL PGTBOX('ABCNTSZYXHO',0.0,0,'ABCNTS',0.0,0) CALL PGSLS(4) DO I=10,90,10 ! lineas horizontales cada 10 grados de elevac. CALL PGMOVE(XMIN,REAL(I)) CALL PGDRAW(XMAX,REAL(I)) END DO ELSE IF(LSECZ) THEN CALL PGSCI(0) ! Borramos CALL PGSLS(1) CALL PGTBOX('ABCNTSZYXHO',0.0,0,'ABCNTS',0.0,0) CALL PGSLS(4) DO I=10,90,10 ! lineas horizontales cada 10 grados de elevac. CALL PGMOVE(XMIN,REAL(I)) CALL PGDRAW(XMAX,REAL(I)) END DO CALL PGSCI(1) ! Pintamos IF(LUNPLOT) CALL PGSCI(0) CALL PGSLS(1) CALL PGTBOX('ABCNTSZYXHO',0.0,0,'ABC',0.0,0) CALL PGSLS(4) DO I=1,10 Y=ACOS(1./SECZ(I)) ! Lineas horizontales para marcar las Y=Y*180.D0/PI ! sec Z (masas de aire) Y=90.-Y CALL PGMOVE(XMIN,Y) CALL PGDRAW(XMAX,Y) CALL PGPTXT(XMIN-300.,Y-1.,0.0,1.0,CHSECZ(I)) END DO END IF IF(LUNPLOT) CALL PGSCI(0) CALL PGSLS(4) DO I=0,172800,3600 ! Lineas de puntos verticales cada hora IF(I.GT.XMIN.AND.I.LT.XMAX) THEN CALL PGMOVE(REAL(I),0.) CALL PGDRAW(REAL(I),90.) END IF END DO CALL PGSLS(2) ! Lineas verticales a trazos para marcar los crepusc CALL PGMOVE(XCREPV,0.) CALL PGDRAW(XCREPV,90.) CALL PGMOVE(XCREPM,0.) CALL PGDRAW(XCREPM,90.) CALL PGSLS(1) END IF C Labels con tiempos de orto y ocaso en TU y TSL IF(LUNPLOT) THEN CALL PGSCI(0) ELSE CALL PGSCI(1) END IF CALL PGSCH(0.8) CALL PGPTXT(XMIN,91.,0.0,0.0,'Sunset') CALL R2CH(TUOCASO2,CHOCASOTU) CALL PGPTXT(XMIN,94.,0.0,0.0,CHOCASOTU//' UT') CALL R2CH(TSLOCASO,CHOCASOTSL) CALL PGPTXT(XMIN,97.,0.0,0.0,CHOCASOTSL//' LST') CALL PGPTXT(XMAX,91.,0.0,1.0,'Sunrise') CALL R2CH(TUORTO2,CHORTOTU) CALL PGPTXT(XMAX,94.,0.0,1.0,CHORTOTU//' UT ') CALL R2CH(TSLORTO,CHORTOTSL) CALL PGPTXT(XMAX,97.,0.0,1.0,CHORTOTSL//' LST') CALL PGSCH(1.0) C IF(LUNPLOT) GOTO 50 C C---------------------------------------------------------------------- C Calculamos la posicion de la Luna para cada hora C TSL0=TS0+LONG ! Tiempo sidereo local a las 0h TU C DO I=1,48 CALL NDAYS(IYEAR,IMONTH,IDAY,I,0,0,DAYS) !Dias transcurr del anio CALL MOONPOS(FJ0,IYEAR,DAYS,MALFA,MDELTA,FASE(I)) TSL=TSL0+DBLE(I)+DBLE(I)*2.737909D-3 IF(TSL.GT.24.D0) TSL=TSL-24.D0 IF(TSL.GT.24.D0) TSL=TSL-24.D0 CALL ECU2HOR(TSL,LAT,MALFA,MDELTA,MOONAZI(I),MOONALT(I)) IF(MOONALT(I).GT.0.D0) THEN IF(I*3600.GT.XMIN.AND.I*3600.LT.XMAX) THEN LMOON(I)=.TRUE. RADIO(I)=(90.-MOONALT(I))/100.D0 IF(MOONAZI(I).GT.90.D0) THEN MOONAZI(I)=360.D0+90.D0-MOONAZI(I) ELSE MOONAZI(I)=90.D0-MOONAZI(I) END IF XMAZI(I)=RADIO(I)*COS(MOONAZI(I)*PI/180.D0) YMAZI(I)=RADIO(I)*SIN(MOONAZI(I)*PI/180.D0) ELSE LMOON(I)=.FALSE. END IF ELSE LMOON(I)=.FALSE. END IF END DO C CALL PGSCI(5) IF(LUT) THEN DO I=1,48 IF(LMOON(I)) THEN CALL PGSCH(1.5) YYY=MOONALT(I) CALL PGPT1(I*3600.,YYY,2291) WRITE(CHF,'(F4.2)')FASE(I) CALL PGSCH(0.6) CALL PGTEXT(I*3600.,YYY-0.5,CHF) END IF END DO ELSE IF(LAZI) THEN CALL PGSVP(0.24,0.82,0.20,0.94) CALL PGSWIN(-1.,1.,-1.,1.) DO I=1,48 IF(LMOON(I)) THEN CALL PGSCH(1.5) CALL PGPT1(XMAZI(I),YMAZI(I),2291) WRITE(CHF,'(F4.2)')FASE(I) CALL PGSCH(0.6) CALL PGTEXT(XMAZI(I),YMAZI(I)-0.01,CHF) END IF END DO END IF C CALL PGSCI(1) CALL PGSCH(1.) C C---------------------------------------------------------------------- C Miramos en el directorio ./ a ver si hay un fichero "vesinput" C que se utilizara como lista de objetos para pintar en la C grafica al lanzar ves INQUIRE(FILE='vesinput',EXIST=VESIN) C---------------------------------------------------------------------- C Loop para elegir ventanita con el cursor. C IF(LNSTAR) GOTO 500 C XX=1.0 YY=0.5 NSTAR=0 C 70 IF(VESIN) GOTO 100 CALL PGSVP(0.0,1.0,0.0,1.0) CALL PGSWIN(0.0,2.0,0.0,1.0) CALL PGSLW(8) CALL PGSCI(4) DO 90 N=1,1000 CALL PGBAND(0,0,XX,YY,XX,YY,CH) IF(ICHAR(CH).EQ.0) GOTO 2000 CALL PGQWIN(X1,X2,Y1,Y2) DO 80 NN=1,NWIN IF (WIN(1,NN).LE.XX .AND. WIN(2,NN).GE.XX .AND. > WIN(3,NN).LE.YY .AND. WIN(4,NN).GE.YY) THEN IF(CH.EQ.'D') THEN CALL HELP(NN) GOTO 70 END IF IF(NN.EQ.1) GOTO 100 IF(NN.EQ.2) GOTO 200 IF(NN.EQ.3) GOTO 300 IF(NN.EQ.4) GOTO 400 IF(NN.EQ.5) GOTO 500 IF(NN.EQ.6) THEN LNSTAR=.FALSE. GOTO 600 END IF IF(NN.EQ.7) GOTO 700 IF(NN.EQ.8) GOTO 800 IF(NN.EQ.9) GOTO 900 IF(NN.EQ.10) GOTO 1000 IF(NN.EQ.11) GOTO 1100 IF(NN.EQ.12) GOTO 1200 IF(NN.EQ.13) GOTO 1300 IF(NN.EQ.14) GOTO 1400 IF(NN.EQ.15) GOTO 1500 IF(NN.EQ.16) GOTO 1600 IF(NN.EQ.17) GOTO 600 IF(NN.EQ.18) GOTO 600 END IF 80 CONTINUE CALL IOWIN1('No window selected ') 90 CONTINUE C---------------------------------------------------------------------- C Si pinchamos con el boton izquierdo introducimos el nombre de una C nueva estrella. Si es el boton derecho podemos dar el nombre de C un fichero que contenga varias estrellas. C 100 IF(VESIN) WRITE(CH,'(A1)')'X' IF(CH.EQ.'A') THEN ! Boton izquierdo del raton CALL PGRECT(WIN(1,1),WIN(2,1),WIN(3,1),WIN(4,1)) CALL PGRSTR(WIN(1,1),WIN(3,1)+0.01,VALUE(1),TEXTBCK,NN,0) CALL PGSCI(1) CALL PGSLW(8) CALL PGRECT(WIN(1,1),WIN(2,1),WIN(3,1),WIN(4,1)) CALL PGSLW(8) GOTO 200 ELSE IF(CH.EQ.'D') THEN ! Boton central del raton CALL IOWIN1('AYUDA ') ELSE IF(CH.EQ.'X') THEN ! Boton derecho del raton IF(VESIN) THEN WRITE(INFILE,'(A11)') 'vesinput ' VESIN=.FALSE. GOTO 160 END IF CALL IOWIN2('File format : ', > 'Name(4:14),R.A.(16:24),DEC.(26:34),EQUI.(36:41) ') GOTO 150 140 CALL IOWIN1('No file with this name, type again ... ') 150 CALL PGSCI(4) CALL PGSFS(2) CALL PGSLW(8) CALL PGRECT(WIN(1,1),WIN(2,1),WIN(3,1),WIN(4,1)) CALL PGRSTR(WIN(1,1),WIN(3,1)+0.01,VALUE(1),TEXTBCK,NN,0) CALL PGSCI(1) CALL PGSLW(8) CALL PGRECT(WIN(1,1),WIN(2,1),WIN(3,1),WIN(4,1)) CALL PGSLW(1) VALUE(1)=TEXTBCK INFILE=TEXTBCK IF(INFILE.EQ.' ') THEN CALL IOWIN1(' No filename typed ... OK ') GOTO 70 END IF 160 OPEN(21,FILE=INFILE,STATUS='OLD',FORM='FORMATTED',ERR=140) DO I=1,100 NSTAR=NSTAR+1 READ(21,'(A41)',END=190) LINE STARNAME(NSTAR)=LINE(4:14) STARRA(NSTAR)=LINE(16:24) STARDEC(NSTAR)=LINE(26:34) EQUINOX(NSTAR)=LINE(36:41) CALL CH2DEC(STARRA(I),ARII(NSTAR),2) CALL CH2DEC(STARDEC(I),DECII(NSTAR),2) END DO CALL IOWIN1(' More than 100 stars in this file ') ELSE GOTO 70 END IF 190 CLOSE(21) NSTAROLD=NSTAR-1 NSTAR=0 GOTO 500 C---------------------------------------------------------------------- C Rellenamos el valor de la ascension recta de la estrella. C 200 CALL PGSCI(4) CALL PGRECT(WIN(1,2),WIN(2,2),WIN(3,2),WIN(4,2)) CALL PGRSTR(WIN(1,2),WIN(3,2)+0.01,VALUE(2),TEXTBCK,NN,1) VALUE(2)=TEXTBCK CALL PGSCI(1) CALL PGSLW(8) CALL PGRECT(WIN(1,2),WIN(2,2),WIN(3,2),WIN(4,2)) C IF(VALUE(2)(1:1).NE.'-'.AND.VALUE(2)(1:1).NE.'+' > .AND.VALUE(2)(1:1).NE.' ') GOTO 250 IF(VALUE(2)(4:4).NE.' '.OR.VALUE(2)(7:7).NE.' ') GOTO 250 DO J=2,9 IF(TESTN(VALUE(2)(J:J)).EQ.0.AND. > VALUE(2)(J:J).NE.' ') GOTO 250 END DO C IF(NN.EQ.1) THEN GOTO 300 ELSE GOTO 70 END IF C 250 CALL IOWIN1(' Please use format : 00 00 00 ') GOTO 200 C---------------------------------------------------------------------- C Rellenamos el valor de la declinacion de la estrella. C 300 CALL PGSCI(4) CALL PGRECT(WIN(1,3),WIN(2,3),WIN(3,3),WIN(4,3)) CALL PGRSTR(WIN(1,3),WIN(3,3)+0.01,VALUE(3),TEXTBCK,NN,0) VALUE(3)=TEXTBCK CALL PGSCI(1) CALL PGSLW(8) CALL PGRECT(WIN(1,3),WIN(2,3),WIN(3,3),WIN(4,3)) C IF(VALUE(3)(1:1).NE.'-'.AND.VALUE(3)(1:1).NE.'+' > .AND.VALUE(3)(1:1).NE.' ') THEN C ! Borramos el texto CALL PGPTXT (WIN(1,3),WIN(3,3)+0.01,0.0,0.0,VALUE(3)) DO J=9,2,-1 VALUE(3)(J:J)=VALUE(3)(J-1:J-1) END DO VALUE(3)(1:1)='+' CALL PGSCI(4) CALL PGSLW(1) CALL PGPTXT (WIN(1,3),WIN(3,3)+0.01,0.0,0.0,VALUE(3)) CALL PGSCI(1) END IF IF(VALUE(3)(4:4).NE.' '.OR.VALUE(3)(7:7).NE.' ') GOTO 350 DO J=2,9 IF(TESTN(VALUE(3)(J:J)).EQ.0.AND. > VALUE(3)(J:J).NE.' ') GOTO 350 END DO C IF(NN.EQ.1) THEN NN=4 GOTO 400 ELSE GOTO 70 END IF C 350 CALL IOWIN1('Please use format : +00 00 00 ') GOTO 300 C---------------------------------------------------------------------- C Rellenamos el valor del equinocio para la A.R y DEC de la nueva estrella C 400 CALL PGSCI(4) CALL PGRECT(WIN(1,4),WIN(2,4),WIN(3,4),WIN(4,4)) CALL PGRSTR(WIN(1,4),WIN(3,4)+0.01,VALUE(4),TEXTBCK,NN,0) CALL PGSCI(1) CALL PGSLW(8) CALL PGRECT(WIN(1,4),WIN(2,4),WIN(3,4),WIN(4,4)) CALL PGSLW(1) VALUE(4)=TEXTBCK C DO J=1,4 IF(TESTN(VALUE(4)(J:J)).EQ.0.AND. > VALUE(4)(J:J).NE.' ') GOTO 450 END DO C GOTO 70 C 450 CALL IOWIN1('Please use format : 0000 ') GOTO 400 C---------------------------------------------------------------------- C Trazamos la posicion de la nueva estrella que tiene las coordenadas C que se ven en las correspondientes ventanas. La grafica corresponde C a los ejes que estan resaltados. 500 LNSTAR=.TRUE. IF(NN.EQ.5) THEN CALL PGSLW(12) CALL PGSCI(4) CALL PGRECT(WIN(1,5),WIN(2,5),WIN(3,5),WIN(4,5)) NSTAR=NSTAR+1 STARNAME(NSTAR)=VALUE(1) STARRA(NSTAR)=VALUE(2) STARDEC(NSTAR)=VALUE(3) EQUINOX(NSTAR)=VALUE(4) CALL CH2DEC(VALUE(2),ARII(NSTAR),2) CALL CH2DEC(VALUE(3),DECII(NSTAR),2) END IF C READ(VALUE(4)(1:4),*) TII IF(IMONTH.LE.2)THEN MONTH=AINT((IMONTH-1)*63./2.) ELSE MONTH=AINT((IMONTH+1)*30.6)-63. END IF MONTH=MONTH+IDAY TFF=REAL(IYEAR+MONTH/365.25) WRITE(CTFF,'(F7.2)') TFF C IF(NSTAR.EQ.0) THEN ! Hemos cambiado alguna ventana de la 7 a la 11 CALL PGSCI(4) ! o hemos metido un listado de estrellas nuevas CALL PGSLW(2) DO J=1,NSTAROLD c TVISINI(J)=DREAL(XMIN) c TVISFIN(J)=DREAL(XMAX) TVISINI(J)=XMIN TVISFIN(J)=XMAX CALL CHANEQUI(TII,TFF,ARII(J),ARFF(J), > DECII(J),DECFF(J)) DO I=1,100 TUPS(J,I)=XMIN+(XMAX-XMIN)*REAL(I-1)/99. TUPH(I)=TUPS(J,I)/3600. TSL=TS0+TUPH(I)+TUPH(I)*2.737909E-3+LONG CALL ECU2HOR(TSL,LAT,ARFF(J),DECFF(J), > AZI(J,I),ALT(J,I)) TUPSU(I)=TUPS(J,I) ALTU(I)=ALT(J,I) AZIU(I)=AZI(J,I) RADIO(I)=(90.-ALTU(I))/100. IF(AZIU(I).GT.90.) THEN AZIU(I)=360.+90.-AZIU(I) ELSE AZIU(I)=90.-AZIU(I) END IF XAZI(J,I)=RADIO(I)*COS(AZIU(I)*PI/180.) YAZI(J,I)=RADIO(I)*SIN(AZIU(I)*PI/180.) END DO WRITE(CH2,'(I2)') J CALL PGSCH(0.5) IF(LUT) THEN CALL PGSVP(0.22,0.85,0.22,0.90) CALL PGSWIN(XMIN,XMAX,YMIN,YMAX) CALL PGLINE(100,TUPSU,ALTU) DO I=1,9 IF(ALTU(10*I+J).GT.0.AND.ALTU(10*I+J).LT.89) THEN CALL PGTEXT(TUPSU(10*I+J),ALTU(10*I+J),CH2) END IF END DO ELSE IF(LAZI) THEN CALL PGSVP(0.24,0.82,0.20,0.94) CALL PGSWIN(-1.,1.,-1.,1.) DO I=1,100 IF(I.GT.1.AND.ALTU(I-1).GT.0.AND.ALTU(I).GT.0.)THEN CALL PGMOVE(XAZI(J,I-1),YAZI(J,I-1)) CALL PGDRAW(XAZI(J,I),YAZI(J,I)) END IF END DO DO I=1,9 IF(ALTU(10*I+J).GT.0.AND.ALTU(10*I+J).LT.89) THEN CALL PGTEXT(XAZI(J,10*I+J),YAZI(J,10*I+J),CH2) END IF END DO END IF CALL PGSCH(1.0) END DO CALL PGSVP(0.0,1.0,0.0,1.0) CALL PGSWIN(0.0,2.0,0.0,1.0) NSTAR=NSTAROLD GOTO 70 ELSE c TVISINI(NSTAR)=DREAL(XMIN) c TVISFIN(NSTAR)=DREAL(XMAX) TVISINI(NSTAR)=XMIN TVISFIN(NSTAR)=XMAX DO I=1,100 TUPS(NSTAR,I)=XMIN+(XMAX-XMIN)*REAL(I-1)/99. TUPH(I)=TUPS(NSTAR,I)/3600. TSL=TS0+TUPH(I)+TUPH(I)*2.737909D-3+LONG CALL CHANEQUI(TII,TFF,ARII(NSTAR),ARFF(NSTAR), > DECII(NSTAR),DECFF(NSTAR)) CALL ECU2HOR(TSL,LAT,ARFF(NSTAR),DECFF(NSTAR), > AZI(NSTAR,I),ALT(NSTAR,I)) TUPSU(I)=TUPS(NSTAR,I) ALTU(I)=ALT(NSTAR,I) AZIU(I)=AZI(NSTAR,I) IF(ALTU(I).GT.0.) THEN RADIO(I)=(90.-ALTU(I))/100. IF(AZIU(I).GT.90.) THEN AZIU(I)=360.+90.-AZIU(I) ELSE AZIU(I)=90.-AZIU(I) END IF XAZI(NSTAR,I)=RADIO(I)*COS(AZIU(I)*PI/180.) YAZI(NSTAR,I)=RADIO(I)*SIN(AZIU(I)*PI/180.) END IF END DO C CALL PGSCI(4) CALL PGSLW(2) CALL PGSVP(0.22,0.85,0.22,0.90) WRITE(CH2,'(I2)') NSTAR CALL PGSCH(0.5) IF(LUT) THEN CALL PGSVP(0.22,0.85,0.22,0.90) CALL PGSWIN(XMIN,XMAX,YMIN,YMAX) CALL PGLINE(100,TUPSU,ALTU) DO I=1,8 IF(ALTU(10*I+NSTAR).GT.0.) THEN CALL PGTEXT(TUPSU(10*I+NSTAR),ALTU(10*I+NSTAR),CH2) END IF END DO ELSE IF(LAZI) THEN CALL PGSVP(0.24,0.82,0.20,0.94) CALL PGSWIN(-1.,1.,-1.,1.) DO I=1,100 IF(I.GT.1.AND.ALTU(I-1).GT.0.AND.ALTU(I).GT.0.)THEN CALL PGMOVE(XAZI(NSTAR,I-1),YAZI(NSTAR,I-1)) CALL PGDRAW(XAZI(NSTAR,I),YAZI(NSTAR,I)) END IF END DO DO I=1,9 IF(ALTU(10*I+J).GT.0.) THEN CALL PGTEXT(XAZI(NSTAR,10*I+J),YAZI(NSTAR,10*I+J),CH2) END IF END DO END IF CALL PGSCH(1.0) CALL PGSVP(0.0,1.0,0.0,1.0) CALL PGSWIN(0.0,2.0,0.0,1.0) END IF C CALL PGSCI(1) CALL PGSLW(12) CALL PGRECT(WIN(1,5),WIN(2,5),WIN(3,5),WIN(4,5)) CALL PGSLW(1) GOTO 70 C---------------------------------------------------------------------- C Quitamos las trazas de las estrellas calculadas hasta ahora. C 600 CALL PGSLW(12) CALL PGSCI(4) CALL PGSLS(1) CALL PGRECT(WIN(1,6),WIN(2,6),WIN(3,6),WIN(4,6)) C CALL PGSCI(0) CALL PGSLW(2) CALL PGSVP(0.22,0.85,0.22,0.90) CALL PGSWIN(XMIN,XMAX,YMIN,YMAX) DO J=1,NSTAR DO I=1,100 TUPSU(I)=TUPS(J,I) ALTU(I)=ALT(J,I) END DO WRITE(CH2,'(I2)') J CALL PGSCH(0.5) IF(LUT) THEN CALL PGLINE(100,TUPSU,ALTU) DO I=1,9 IF(ALTU(10*I+J).GT.0.AND.ALTU(10*I+J).LT.89) THEN CALL PGTEXT(TUPSU(10*I+J),ALTU(10*I+J),CH2) END IF END DO ELSE CALL PGSVP(0.24,0.82,0.20,0.94) CALL PGSWIN(-1.,1.,-1.,1.) DO I=1,100 CALL PGDRAW(XAZI(J,I),YAZI(J,I)) END DO DO I=1,9 IF(ALTU(10*I+J).GT.0.AND.ALTU(10*I+J).LT.89) THEN CALL PGTEXT(XAZI(J,10*I+J),YAZI(J,10*I+J),CH2) END IF END DO END IF CALL PGSCH(1.0) END DO C Borramos la Luna IF(LUT) THEN DO I=1,48 IF(LMOON(I)) THEN CALL PGSCH(1.5) YYY=MOONALT(I) CALL PGPT1(I*3600.,YYY,2291) WRITE(CHF,'(F4.2)')FASE(I) CALL PGSCH(0.6) CALL PGTEXT(I*3600.,YYY-0.5,CHF) END IF END DO ELSE IF(LAZI) THEN CALL PGSVP(0.24,0.82,0.20,0.94) CALL PGSWIN(-1.,1.,-1.,1.) DO I=1,48 IF(LMOON(I)) THEN CALL PGSCH(1.5) CALL PGPT1(XMAZI(I),YMAZI(I),2291) WRITE(CHF,'(F4.2)')FASE(I) CALL PGSCH(0.6) CALL PGTEXT(XMAZI(I),YMAZI(I)-0.01,CHF) END IF END DO END IF CALL PGSCH(1.) IF(NN.EQ.6) THEN DO I=1,48 LMOON(I)=.FALSE. END DO END IF C CALL PGSLW(12) CALL PGSCI(1) CALL PGSVP(0.0,1.0,0.0,1.0) CALL PGSWIN(0.0,2.0,0.0,1.0) CALL PGRECT(WIN(1,6),WIN(2,6),WIN(3,6),WIN(4,6)) CALL PGSLW(1) C IF(NN.EQ.17) GOTO 1700 IF(NN.EQ.18) GOTO 1800 NSTAROLD=NSTAR NSTAR=0 IF(NN.EQ.7.AND.CH.EQ.'A') GOTO 800 IF(NN.EQ.7.AND.CH.EQ.'X') GOTO 60 IF(NN.GE.8.AND.NN.LE.11) GOTO 60 GOTO 70 C C---------------------------------------------------------------------- C Si pinchamos con el boton izquierdo podemos meter los datos de un C nuevo observatorio. Hemos de rellenar las cuatro ventanas corres- C pondientes. Nombre, longitud, latitud y altura. C Si pinchamos con el boton derecho vamos elegiendo otro observatorio C de los diez que hay en el programa (DATA). C 700 IF(CH.EQ.'A') THEN ! Boton izquierdo del raton CALL PGRECT(WIN(1,7),WIN(2,7),WIN(3,7),WIN(4,7)) CALL PGRSTR(WIN(1,7),WIN(3,7)+0.01,VALUE(7),TEXTBCK,NN,0) CALL PGSCI(1) CALL PGSLW(8) CALL PGRECT(WIN(1,7),WIN(2,7),WIN(3,7),WIN(4,7)) VALUE(7)=TEXTBCK IF(NSTAR.GT.0) GOTO 600 GOTO 800 ELSE IF(CH.EQ.'X') THEN ! Boton derecho del raton CALL PGSLW(1) CALL PGSCI(1) CALL PGPTXT(WIN(1,7),WIN(3,7)+0.01,0.0,0.0,VALUE(7)) CALL PGPTXT(WIN(1,8),WIN(3,8)+0.01,0.0,0.0,VALUE(8)) CALL PGPTXT(WIN(1,9),WIN(3,9)+0.01,0.0,0.0,VALUE(9)) CALL PGPTXT(WIN(1,10),WIN(3,10)+0.01,0.0,0.0,VALUE(10)) NNOBS=NNOBS+1 IF(NNOBS.GT.10) NNOBS=1 VALUE(7)=OBSNAME(NNOBS) VALUE(8)=OBSLONG(NNOBS) CALL CH2DEC(VALUE(8),LONG,1) LONG=LONG/15.D0 VALUE(9)=OBSLAT(NNOBS) CALL CH2DEC(VALUE(9),LAT,2) VALUE(10)=OBSHEI(NNOBS) READ(VALUE(10)(1:4),*)ALTOBS CALL PGSCI(4) CALL PGPTXT(WIN(1,7),WIN(3,7)+0.01,0.0,0.0,VALUE(7)) CALL PGPTXT(WIN(1,8),WIN(3,8)+0.01,0.0,0.0,VALUE(8)) CALL PGPTXT(WIN(1,9),WIN(3,9)+0.01,0.0,0.0,VALUE(9)) CALL PGPTXT(WIN(1,10),WIN(3,10)+0.01,0.0,0.0,VALUE(10)) CALL PGSCI(1) LUNPLOT=.TRUE. GOTO 600 ELSE GOTO 70 END IF C C---------------------------------------------------------------------- C Rellenamos el valor de la longitud del Observatorio C 800 CALL PGSCI(4) CALL PGSLW(8) CALL PGRECT(WIN(1,8),WIN(2,8),WIN(3,8),WIN(4,8)) CALL PGRSTR(WIN(1,8),WIN(3,8)+0.01,VALUE(8),TEXTBCK,NN,0) CALL PGSLW(8) CALL PGSCI(1) CALL PGRECT(WIN(1,8),WIN(2,8),WIN(3,8),WIN(4,8)) VALUE(8)=TEXTBCK C IF(VALUE(8)(1:1).NE.'-'.AND.VALUE(8)(1:1).NE.'+' > .AND.VALUE(8)(1:1).NE.' ') GOTO 850 IF(VALUE(8)(5:5).NE.' '.OR.VALUE(8)(8:8).NE.' ') GOTO 850 DO J=2,10 IF(TESTN(VALUE(8)(J:J)).EQ.0.AND. > VALUE(8)(J:J).NE.' ') GOTO 850 END DO C CALL CH2DEC(VALUE(8),LONG,1) LONG=LONG/15.D0 C IF(NN.EQ.7) THEN GOTO 900 ELSE LUNPLOT=.TRUE. GOTO 600 ENDIF C 850 CALL IOWIN1(' Please use format : +000 00 00 ') GOTO 800 C---------------------------------------------------------------------- C Rellenamos el valor de la latitud del Observatorio C 900 CALL PGSCI(4) CALL PGRECT(WIN(1,9),WIN(2,9),WIN(3,9),WIN(4,9)) CALL PGRSTR(WIN(1,9),WIN(3,9)+0.01,VALUE(9),TEXTBCK,NN,0) CALL PGSCI(1) CALL PGSLW(8) CALL PGRECT(WIN(1,9),WIN(2,9),WIN(3,9),WIN(4,9)) VALUE(9)=TEXTBCK C IF(VALUE(9)(1:1).NE.'-'.AND.VALUE(9)(1:1).NE.'+' > .AND.VALUE(9)(1:1).NE.' ') GOTO 950 IF(VALUE(9)(4:4).NE.' '.OR.VALUE(9)(7:7).NE.' ') GOTO 950 DO J=2,9 IF(TESTN(VALUE(9)(J:J)).EQ.0.AND. > VALUE(9)(J:J).NE.' ') GOTO 950 END DO C CALL CH2DEC(VALUE(9),LAT,2) C IF(NN.EQ.7) THEN GOTO 1000 ELSE LUNPLOT=.TRUE. GOTO 600 END IF C 950 CALL IOWIN1(' Please use format : +00 00 00 ') GOTO 900 C---------------------------------------------------------------------- C Rellenamos el valor de la altura del Observatorio C 1000 CALL PGSCI(4) CALL PGRECT(WIN(1,10),WIN(2,10),WIN(3,10),WIN(4,10)) CALL PGRSTR(WIN(1,10),WIN(3,10)+0.01,VALUE(10),TEXTBCK,NN,0) CALL PGSCI(1) CALL PGSLW(8) CALL PGRECT(WIN(1,10),WIN(2,10),WIN(3,10),WIN(4,10)) CALL PGSLW(1) VALUE(10)=TEXTBCK C DO J=1,4 IF(TESTN(VALUE(10)(J:J)).EQ.0.AND. > VALUE(10)(J:J).NE.' ') GOTO 1050 END DO C IF(VALUE(10)(1:4).EQ.' ') THEN ALTOBS=0. ELSE READ(VALUE(10)(1:4),*) ALTOBS END IF C LUNPLOT=.TRUE. NN=10 GOTO 600 C 1050 CALL IOWIN1(' Please use format : 0000m ') GOTO 1000 C---------------------------------------------------------------------- C Cambiamos la fecha C 1100 CALL PGSCI(4) CALL PGSLW(8) CALL PGRECT(WIN(1,11),WIN(2,11),WIN(3,11),WIN(4,11)) CALL PGRSTR(WIN(1,11),WIN(3,11)+0.01,VALUE(11),TEXTBCK,NN,0) CALL PGSCI(1) CALL PGSLW(8) CALL PGRECT(WIN(1,11),WIN(2,11),WIN(3,11),WIN(4,11)) CALL PGSLW(1) VALUE(11)=TEXTBCK C DO J=1,2 IF(TESTN(VALUE(11)(J:J)).EQ.0.) GOTO 1150 END DO DO J=4,5 IF(TESTN(VALUE(11)(J:J)).EQ.0.) GOTO 1150 END DO DO J=7,10 IF(TESTN(VALUE(11)(J:J)).EQ.0.) GOTO 1150 END DO C READ(VALUE(11)(1:2),'(I2)')IDAY READ(VALUE(11)(4:5),'(I2)')IMONTH READ(VALUE(11)(7:10),'(I4)')IYEAR C LUNPLOT=.TRUE. GOTO 600 C 1150 CALL IOWIN1(' Please use format : DD-MM-YYYY ') GOTO 1100 C---------------------------------------------------------------------- C Escribimos el fichero con el listado de las estrellas. C No Orden, Nombre, A.R., Dec. y Equinocio. C 1200 CALL PGSLW(12) CALL PGSCI(4) CALL PGRECT(WIN(1,12),WIN(2,12),WIN(3,12),WIN(4,12)) OPEN(20,FILE='ves.lis',STATUS='UNKNOWN', > FORM='FORMATTED') DO I=1,NSTAR WRITE(20,'(I2,A1,A11,A1,A9,A1,A9,A1,A6)') > I,',',STARNAME(I),',',STARRA(I),',',STARDEC(I),',',EQUINOX(I) END DO C CLOSE(20) CALL PGSCI(1) CALL PGRECT(WIN(1,12),WIN(2,12),WIN(3,12),WIN(4,12)) GOTO 70 C---------------------------------------------------------------------- C Generamos un fichero postscript con la informacion que tenemos C actualmente en pantalla 1300 CALL PGSLW(12) CALL PGSCI(4) CALL PGRECT(WIN(1,13),WIN(2,13),WIN(3,13),WIN(4,13)) CALL PGSLW(1) CALL PGSCI(1) C DEV(2)=PGOPEN('ves.ps/vps') CALL PGSLCT(DEV(2)) CALL PGSVP(0.,1.,0.,1.) CALL PGSWIN(0.,1.,0.,1.) CALL PGSCF(2) CALL PGSCH(0.55) CALL PGTEXT(0.85,0.016,COPYR) CALL PGSFS(2) ! Dibujamos un doble borde exterior CALL PGRECT(.01,.99,.01,.99) CALL PGRECT(.014,.986,.013,.987) C Escribimos la informacion del observatorio CALL PGSCH(0.8) CALL PGTEXT(.125,.960,'Observatory :') CALL PGTEXT(.125,.942,'Longitud :') CALL PGTEXT(.125,.924,'Latitud :') CALL PGTEXT(.125,.906,'Height :') CALL PGPTXT(.400,.960,0.0,1.0,VALUE(7)) CALL PGPTXT(.400,.942,0.0,1.0,VALUE(8)) CALL PGPTXT(.400,.924,0.0,1.0,VALUE(9)) CALL PGPTXT(.400,.906,0.0,1.0,VALUE(10)) C Escribimos info de la fecha CALL PGTEXT(.600,.960,'Date :') CALL PGPTXT(.875,.960,0.0,1.0,VALUE(11)) CALL PGTEXT(.600,.942,'JD (0h UT) :') WRITE(CFJ0,'(F10.2)')FJ0 CALL PGPTXT(.875,.942,0.0,1.0,CFJ0) CALL PGTEXT(.600,.924,'Night length :') CALL PGTEXT(.600,.906,'Twilight free :') NIGHTL=24.-TUOCASO2+TUORTO2 IF(NIGHTL.GT.24.) NIGHTL=NIGHTL-24. NIGHTLF=24.-TUCREPV2+TUCREPM2 IF(NIGHTLF.GT.24.) NIGHTLF=NIGHTLF-24. CALL R2CH(NIGHTL,CNIGHTL) CALL R2CH(NIGHTLF,CNIGHTLF) CALL PGPTXT(.875,.924,0.0,1.0,CNIGHTL) CALL PGPTXT(.875,.906,0.0,1.0,CNIGHTLF) CALL PGMOVE(.025,.900) CALL PGDRAW(.975,.900) CALL PGMOVE(.025,.897) CALL PGDRAW(.975,.897) C Dibujamos el plot. CALL PGSVP(.075,.925,.340,.840) CALL PGSWIN(XMIN,XMAX,YMIN,YMAX) CALL PGSCH(0.6) CALL PGPTXT(XMIN,91.,0.0,0.0,'Sunset') CALL R2CH(TUOCASO2,CHOCASOTU) CALL PGPTXT(XMIN,93.,0.0,0.0,CHOCASOTU//' UT') CALL R2CH(TSLOCASO,CHOCASOTSL) CALL PGPTXT(XMIN,95.,0.0,0.0,CHOCASOTSL//' LST') CALL PGPTXT(XMAX,91.,0.0,1.0,'Sunrise') CALL R2CH(TUORTO2,CHORTOTU) CALL PGPTXT(XMAX,93.,0.0,1.0,CHORTOTU//' UT ') CALL R2CH(TSLORTO,CHORTOTSL) CALL PGPTXT(XMAX,95.,0.0,1.0,CHORTOTSL//' LST') CALL PGSCH(1.0) C IF(LUT) THEN ! Altura frente a T.U. CALL PGTBOX('ABCNTSZYXHO',0.0,0,'ABCN',0.0,0) CALL PGPTXT(XMIN-2200.,42.,90.,0.0,'Elevation') CALL PGPTXT(XMIN-2000.,-4.5,0.,0.0,'U.T.') CALL PGSLS(4) DO I=10,90,10 ! lineas horizontales cada 10 grados de elevac. CALL PGMOVE(XMIN,REAL(I)) CALL PGDRAW(XMAX,REAL(I)) END DO DO I=0,172800,3600 ! Lineas de puntos verticales cada hora IF(I.GT.XMIN.AND.I.LT.XMAX) THEN CALL PGMOVE(REAL(I),0.) CALL PGDRAW(REAL(I),90.) END IF END DO CALL PGSLS(2) ! Lineas verticales a trazos para marcar los crepusc CALL PGMOVE(XCREPV,0.) CALL PGDRAW(XCREPV,90.) CALL PGMOVE(XCREPM,0.) CALL PGDRAW(XCREPM,90.) CALL PGSLS(1) DO I=1,10 ! Ticks para marcar las sec Z (masas de aire) Y=ACOS(1./SECZ(I)) Y=Y*180./PI Y=90.-Y CALL PGMOVE(XMAX-900.,Y) CALL PGDRAW(XMAX,Y) CALL PGPTXT(XMAX+300.,Y-1.,0.0,0.0,CHSECZ(I)) END DO CALL PGPTXT(XMAX+2000.,72.,90.,0.0,'sec(Z)') DO J=1,NSTAR ! Pintamos las trazas de las estrellas DO I=1,100 TUPSU(I)=TUPS(J,I) ALTU(I)=ALT(J,I) END DO CALL PGSLW(3) CALL PGLINE(100,TUPSU,ALTU) CALL PGSLW(1) WRITE(CH2,'(I2)') J CALL PGSCH(.5) DO I=1,9 IF(ALTU(10*I+J).GT.0.AND.ALTU(10*I+J).LT.89) THEN CALL PGTEXT(TUPSU(10*I+J),ALTU(10*I+J),CH2) END IF END DO CALL PGSCH(1.) END DO C Pintamos la Luna si se ve DO I=1,48 IF(LMOON(I)) THEN CALL PGSCH(1.5) YYY=MOONALT(I) CALL PGPT1(I*3600.,YYY,2291) WRITE(CHF,'(F4.2)')FASE(I) CALL PGSCH(0.6) CALL PGTEXT(I*3600.,YYY-0.5,CHF) END IF END DO ELSE IF(LAZI) THEN CALL PGSVP(.100,.900,.325,.875) CALL PGSWIN(-1.,1.,-1.,1.) CALL PGCIRC(0.,0.,0.9) CALL PGSLS(2) ! Trazamos cruz con lineas de guiones CALL PGMOVE(-.9,0.) CALL PGDRAW(.9,0.) CALL PGMOVE(0.,-0.9) CALL PGDRAW(0.,0.9) CALL PGSCH(2.) CALL PGTEXT(-.05,0.92,'S') CALL PGTEXT(-.05,-1.01,'N') CALL PGTEXT(-1.,0.,'E') CALL PGTEXT(0.92,0.,'W') CALL PGSLS(4) DO I=1,5 ! Circulos de puntos cada 15 grados CALL PGCIRC(0.,0.,REAL(I)*0.15) END DO CALL PGSCH(0.7) CALL PGTEXT(0.0,-.74,'15') CALL PGTEXT(0.0,-.59,'30') CALL PGTEXT(0.0,-.44,'45') CALL PGTEXT(0.0,-.29,'60') CALL PGTEXT(0.0,-.14,'75') CALL PGSCH(1.) CALL PGSLS(1) DO J=1,NSTAR WRITE(CH2,'(I2)') J CALL PGSLW(3) DO I=1,100 RADIO(I)=SQRT(XAZI(J,I)**2+YAZI(J,I)**2) IF(I.GT.1.AND.RADIO(I).LT.0.9.AND.RADIO(I).GT.0.)THEN CALL PGMOVE(XAZI(J,I-1),YAZI(J,I-1)) CALL PGDRAW(XAZI(J,I),YAZI(J,I)) END IF END DO CALL PGSLW(1) CALL PGSCH(0.5) DO I=1,9 IF(RADIO(10*I).LT.0.9.AND.RADIO(10*I).GT.0.) THEN CALL PGTEXT(XAZI(J,10*I),YAZI(J,10*I),CH2) END IF END DO CALL PGSCH(1.0) END DO C Pintamos la Luna si se ve DO I=1,48 IF(LMOON(I)) THEN CALL PGSCH(1.5) CALL PGPT1(XMAZI(I),YMAZI(I),2291) WRITE(CHF,'(F4.2)')FASE(I) CALL PGSCH(0.6) CALL PGTEXT(XMAZI(I),YMAZI(I)-0.015,CHF) END IF END DO END IF C CALL PGSVP(0.,1.,0.,1.) CALL PGSWIN(0.,1.,0.,1.) CALL PGSLS(1) CALL PGMOVE(.025,.300) CALL PGDRAW(.975,.300) CALL PGMOVE(.025,.297) CALL PGDRAW(.975,.297) C Escribimos la lista de objetos IF(NSTAR.GT.0) THEN CALL PGSCH(.6) CALL PGTEXT(.050,.280,' Object') CALL PGTEXT(.235,.280,'\ga Initial \gd ') CALL PGTEXT(.450,.280,'\ga '//CTFF//' \gd ') CALL PGTEXT(.640,.280,'Culmination U.T.') CALL PGTEXT(.840,.280,'Visible') CALL PGMOVE(.050,.275) CALL PGDRAW(.950,.275) DO I=1,NSTAR IF(I.EQ.17) THEN CALL IOWIN1(' More than 16 objects ') GOTO 1350 END IF WRITE(CH2,'(I2)') I CALL PGPTXT(.050,.275-I*.015,0.,0.,CH2) CALL PGPTXT(.070,.275-I*.015,0.,0.,STARNAME(I)) CALL PGPTXT(.280,.275-I*.015,0.,1.,STARRA(I)) CALL PGPTXT(.375,.275-I*.015,0.,1.,STARDEC(I)) CALL RR2CH(ARFF(I),STARRAF(I)) CALL PGPTXT(.500,.275-I*.015,0.,1.,STARRAF(I)) CALL RR2CH(DECFF(I),STARDECF(I)) CALL PGPTXT(.600,.275-I*.015,0.,1.,STARDECF(I)) C Calculamos la altura, direccion y T.U. de la culminacion. CALL CULMINA(LONG,LAT,TS0,ARFF(I),DECFF(I), > ALTCUL,DIRCUL,TUCUL) IF(DIRCUL.NE.'I') THEN WRITE(CALTCUL,'(F5.1)') ALTCUL CALL PGPTXT(.640,.275-I*.015,0.,0.,CALTCUL) CALL PGPTXT(.685,.275-I*.015,0.,0.,DIRCUL) CALL R2CH(TUCUL,CTUCUL) CALL PGPTXT(.785,.275-I*.015,0.,1.,CTUCUL) END IF C Calculamos la visibilidad del objeto CALL ORCASO(DECFF(I),90.D0,ESTADO) IF(ESTADO.EQ.'I') THEN CALL PGTEXT(.835,.275-I*.015,'Invisible') ELSE IF(ESTADO.EQ.'C') THEN CALL PGTEXT(.835,.275-I*.015,'Circumpolar') ELSE CALL TIEMPOU(ARFF(I),HORTO,LONG,TS0,TSL,TUORSTAR) CALL TIEMPOU(ARFF(I),HOCASO,LONG,TS0,TSL,TUOCSTAR) IF(TUOCSTAR.GT.TUORSTAR) THEN IF((TUORSTAR+24.)*3600.LT.TVISFIN(I)) THEN TUORSTAR=TUORSTAR+24. TUOCSTAR=TUOCSTAR+24. END IF TVISINI(I)=DMAX1(TVISINI(I),TUORSTAR*3600.D0) TVISFIN(I)=DMIN1(TVISFIN(I),TUOCSTAR*3600.D0) ELSE TUOCSTAR=TUOCSTAR+24. TVISINI(I)=DMAX1(TVISINI(I),TUORSTAR*3600.D0) TVISFIN(I)=DMIN1(TVISFIN(I),TUOCSTAR*3600.D0) END IF TVIS(I)=(TVISFIN(I)-TVISINI(I))/3600.D0 CALL R2CH(TVIS(I),CTVIS) CALL PGPTXT(.900,.275-I*.015,0.,1.,CTVIS) END IF END DO END IF C 1350 CALL PGSLS(1) CALL PGSCH(1.) C CALL PGCLOS(DEV(2)) CALL PGSLCT(DEV(1)) CALL PGSLW(12) CALL PGRECT(WIN(1,13),WIN(2,13),WIN(3,13),WIN(4,13)) CALL PGSLW(1) GOTO 70 C---------------------------------------------------------------------- C Salimos del programa. C 1400 CALL PGSLW(12) CALL PGSCI(4) CALL PGSLS(1) CALL PGRECT(WIN(1,14),WIN(2,14),WIN(3,14),WIN(4,14)) CALL PGSLW(1) CALL PGSCI(1) GOTO 2000 C C---------------------------------------------------------------------- C Selecionamos la altura sobre el horizonte como coordenada Y. C 1500 CALL PGSCI(1) CALL PGSLW(12) CALL PGSLS(1) CALL PGRECT(WIN(1,16),WIN(2,16),WIN(3,16),WIN(4,16)) CALL PGSCI(4) CALL PGRECT(WIN(1,15),WIN(2,15),WIN(3,15),WIN(4,15)) CALL PGSLW(1) C LSECZ=.FALSE. LELEV=.TRUE. C IF(LUT) THEN GOTO 60 ELSE ! Azimut CALL PGSLS(4) CALL PGSVP(0.24,0.82,0.20,0.94) CALL PGSWIN(-1.,1.,-1.,1.) CALL PGSCI(0) ! Borramos los circulos de puntos que marcan las secZ CALL PGSCH(.7) DO I=1,10 Y=ACOS(1./SECZ(I)) Y=Y*1.8/PI IF(I.EQ.1) CALL PGTEXT(0.,Y,'1.0') IF(I.EQ.2) CALL PGTEXT(0.,Y,'1.1') IF(I.EQ.7) CALL PGTEXT(0.,Y,'2.0') CALL PGCIRC(0.,0.,Y) END DO CALL PGSCI(1) ! Pintamos los circulos de puntos de altura CALL PGSCH(1.) DO I=1,5 CALL PGCIRC(0.,0.,REAL(I)*0.150) END DO CALL PGSCH(0.7) CALL PGTEXT(0.0,-.74,'15') CALL PGTEXT(0.0,-.59,'30') CALL PGTEXT(0.0,-.44,'45') CALL PGTEXT(0.0,-.29,'60') CALL PGTEXT(0.0,-.14,'75') CALL PGSCH(1.) CALL PGSVP(0.0,1.0,0.0,1.0) CALL PGSWIN(0.0,2.0,0.0,1.0) END IF C GOTO 70 C---------------------------------------------------------------------- C Selecionamos sec Z (aprox. masa de aire) como coordenada Y. C 1600 CALL PGSCI(1) CALL PGSLS(1) CALL PGSLW(12) CALL PGRECT(WIN(1,15),WIN(2,15),WIN(3,15),WIN(4,15)) CALL PGSCI(4) CALL PGRECT(WIN(1,16),WIN(2,16),WIN(3,16),WIN(4,16)) CALL PGSLW(1) C LELEV=.FALSE. LSECZ=.TRUE. C IF(LUT) THEN GOTO 60 ELSE ! Azimut CALL PGSLS(1) CALL PGSVP(0.24,0.82,0.20,0.94) CALL PGSWIN(-1.,1.,-1.,1.) CALL PGSCI(0) ! Borramos los circulos de puntos de altura CALL PGSLW(10) CALL PGSCH(1.) DO I=1,5 CALL PGCIRC(0.,0.,REAL(I)*0.150) END DO CALL PGSLW(1) CALL PGSCH(0.7) CALL PGTEXT(0.0,-.74,'15') CALL PGTEXT(0.0,-.59,'30') CALL PGTEXT(0.0,-.44,'45') CALL PGTEXT(0.0,-.29,'60') CALL PGTEXT(0.0,-.14,'75') CALL PGSLS(4) CALL PGSCI(1) ! Pintamos los circulos de puntos que marcan las secZ DO I=1,10 Y=ACOS(1./SECZ(I)) ! Circulos para marcar las Y=Y*1.8/PI ! sec Z (masas de aire) IF(I.EQ.1) CALL PGTEXT(0.,Y,'1.0') IF(I.EQ.2) CALL PGTEXT(0.,Y,'1.1') IF(I.EQ.7) CALL PGTEXT(0.,Y,'2.0') CALL PGCIRC(0.,0.,Y) END DO CALL PGSCH(1.) CALL PGSVP(0.0,1.0,0.0,1.0) CALL PGSWIN(0.0,2.0,0.0,1.0) END IF C GOTO 70 C---------------------------------------------------------------------- C Selecionamos Tiempo Universal como coordenada X C 1700 CALL PGSCI(1) CALL PGSLW(12) CALL PGSLS(1) CALL PGRECT(WIN(1,18),WIN(2,18),WIN(3,18),WIN(4,18)) CALL PGSCI(4) CALL PGRECT(WIN(1,17),WIN(2,17),WIN(3,17),WIN(4,17)) CALL PGSLW(1) C LAZI=.FALSE. LUT=.TRUE. C CALL PGSCI(0) ! Borramos el grafico de azimut CALL PGSVP(0.24,0.82,0.20,0.94) CALL PGSWIN(-1.,1.,-1.,1.) CALL PGSLS(1) CALL PGCIRC(0.,0.,0.9) CALL PGSLS(2) CALL PGMOVE(-.9,0.) CALL PGDRAW(.9,0.) CALL PGMOVE(0.,-0.9) CALL PGDRAW(0.,0.9) CALL PGSCH(2.) CALL PGTEXT(-.05,0.92,'S') CALL PGTEXT(-.05,-1.01,'N') CALL PGTEXT(-1.,0.,'E') CALL PGTEXT(0.92,0.,'W') CALL PGSCH(1.) C CALL PGSLS(4) CALL PGSLW(1) IF(LELEV) THEN DO I=1,5 CALL PGCIRC(0.,0.,REAL(I)*0.15) END DO CALL PGSCH(0.7) CALL PGTEXT(0.0,-.74,'15') CALL PGTEXT(0.0,-.59,'30') CALL PGTEXT(0.0,-.44,'45') CALL PGTEXT(0.0,-.29,'60') CALL PGTEXT(0.0,-.14,'75') CALL PGSCH(1.) ELSE CALL PGSCH(0.7) DO I=1,10 Y=ACOS(1./SECZ(I)) ! Lineas horizontales para marcar las Y=Y*1.8/PI ! sec Z (masas de aire) IF(I.EQ.1) CALL PGTEXT(0.,Y,'1.0') IF(I.EQ.2) CALL PGTEXT(0.,Y,'1.1') IF(I.EQ.7) CALL PGTEXT(0.,Y,'2.0') CALL PGCIRC(0.,0.,Y) END DO CALL PGSCH(1.) END IF C Borramos la Luna DO I=1,48 IF(LMOON(I)) THEN CALL PGSCH(1.5) CALL PGPT1(XMAZI(I),YMAZI(I),2291) WRITE(CHF,'(F4.2)')FASE(I) CALL PGSCH(0.6) CALL PGTEXT(XMAZI(I),YMAZI(I)-0.01,CHF) END IF END DO CALL PGSCH(1.) C NSTAROLD=NSTAR NSTAR=0 C CALL PGSCI(1) CALL PGSLW(1) c GOTO 60 C---------------------------------------------------------------------- C Selecionamos Azimut como coordenada X C 1800 CALL PGSCI(1) CALL PGSLW(12) CALL PGRECT(WIN(1,17),WIN(2,17),WIN(3,17),WIN(4,17)) CALL PGSCI(4) CALL PGRECT(WIN(1,18),WIN(2,18),WIN(3,18),WIN(4,18)) CALL PGSLW(1) C LUT=.FALSE. LAZI=.TRUE. C CALL PGSCI(0) ! Borramos CALL PGSVP(0.22,0.85,0.22,0.90) CALL PGSWIN(XMIN,XMAX,YMIN,YMAX) C IF(LSECZ) THEN CALL PGSLS(1) CALL PGTBOX('ABCNTSZYXHO',0.0,0,'ABC',0.0,0) CALL PGSLS(4) DO I=1,10 Y=ACOS(1./SECZ(I)) ! Lineas horizontales para marcar las Y=Y*180./PI ! sec Z (masas de aire) Y=90.-Y CALL PGMOVE(XMIN,Y) CALL PGDRAW(XMAX,Y) CALL PGPTXT(XMIN-300.,Y-1.,0.0,1.0,CHSECZ(I)) END DO ELSE IF(LELEV) THEN CALL PGSLS(1) CALL PGTBOX('ABCNTSZYXHO',0.0,0,'ABCNTS',0.0,0) CALL PGSLS(4) DO I=10,90,10 ! lineas horizontales cada 10 grados de elevac. CALL PGMOVE(XMIN,REAL(I)) CALL PGDRAW(XMAX,REAL(I)) END DO END IF C CALL PGSLS(4) DO I=0,172800,3600 ! Lineas de puntos verticales cada hora IF(I.GT.XMIN.AND.I.LT.XMAX) THEN CALL PGMOVE(REAL(I),0.) CALL PGDRAW(REAL(I),90.) END IF END DO C CALL PGSLS(2) ! Lineas verticales a trazos para marcar los crepusc CALL PGMOVE(XCREPV,0.) CALL PGDRAW(XCREPV,90.) CALL PGMOVE(XCREPM,0.) CALL PGDRAW(XCREPM,90.) C CALL PGSCI(1) ! Pintamos la trama CALL PGSVP(0.24,0.82,0.20,0.94) CALL PGSWIN(-1.,1.,-1.,1.) CALL PGSLS(1) CALL PGCIRC(0.,0.,0.9) CALL PGSLS(2) ! Trazamos cruz con lineas de guiones CALL PGMOVE(-.9,0.) CALL PGDRAW(.9,0.) CALL PGMOVE(0.,-0.9) CALL PGDRAW(0.,0.9) CALL PGSCH(2.) CALL PGTEXT(-.05,0.92,'S') CALL PGTEXT(-.05,-1.01,'N') CALL PGTEXT(-1.,0.,'E') CALL PGTEXT(0.92,0.,'W') C CALL PGSLS(4) CALL PGSLW(1) CALL PGSCH(1.) IF(LELEV)THEN DO I=1,5 ! Circulos de puntos cada 15 grados CALL PGCIRC(0.,0.,REAL(I)*0.15) END DO CALL PGSCH(0.7) CALL PGTEXT(0.0,-.74,'15') CALL PGTEXT(0.0,-.59,'30') CALL PGTEXT(0.0,-.44,'45') CALL PGTEXT(0.0,-.29,'60') CALL PGTEXT(0.0,-.14,'75') CALL PGSCH(1.) ELSE CALL PGSCH(0.7) DO I=1,10 Y=ACOS(1./SECZ(I)) ! Circulos de puntos para marcar las Y=Y*1.8/PI ! sec Z (masas de aire) IF(I.EQ.1) CALL PGTEXT(0.,Y,'1.0') IF(I.EQ.2) CALL PGTEXT(0.,Y,'1.1') IF(I.EQ.7) CALL PGTEXT(0.,Y,'2.0') CALL PGCIRC(0.,0.,Y) END DO CALL PGSCH(1.0) END IF C CALL PGSCI(4) !Pintamos las trazas de los objetos CALL PGSLS(1) CALL PGSLW(2) DO J=1,NSTAR WRITE(CH2,'(I2)') J CALL PGSVP(0.24,0.82,0.20,0.94) CALL PGSWIN(-1.,1.,-1.,1.) DO I=1,100 RADIO(I)=SQRT(XAZI(J,I)**2+YAZI(J,I)**2) IF(I.GT.1.AND.RADIO(I).LT.0.9.AND.RADIO(I).GT.0.)THEN IF(ABS(XAZI(J,I-1)*YAZI(J,I-1)).GT.0) THEN CALL PGMOVE(XAZI(J,I-1),YAZI(J,I-1)) ELSE CALL PGMOVE(XAZI(J,I),YAZI(J,I)) END IF CALL PGDRAW(XAZI(J,I),YAZI(J,I)) END IF END DO CALL PGSLW(1) CALL PGSCH(0.5) DO I=1,9 IF(RADIO(10*I+J).LT.0.9.AND.RADIO(10*I+J).GT.0.) THEN CALL PGTEXT(XAZI(J,10*I+J),YAZI(J,10*I+J),CH2) END IF END DO CALL PGSCH(1.0) END DO C Ponemos la Luna CALL PGSCI(5) DO I=1,48 IF(LMOON(I)) THEN CALL PGSCH(1.5) CALL PGPT1(XMAZI(I),YMAZI(I),2291) WRITE(CHF,'(F4.2)')FASE(I) CALL PGSCH(0.6) CALL PGTEXT(XMAZI(I),YMAZI(I)-0.01,CHF) END IF END DO C CALL PGSCH(1.) CALL PGSCI(1) CALL PGSLW(1) CALL PGSVP(0.0,1.0,0.0,1.0) CALL PGSWIN(0.0,2.0,0.0,1.0) C GOTO 70 C---------------------------------------------------------------------- C 2000 CALL PGASK(.FALSE.) CALL PGEND C STOP END C C********************************************************************** C Para modificar el texto que haya en una ventana. C SUBROUTINE PGRSTR(X, Y, TEXT,TEXTBCK,NN,SP) C IMPLICIT NONE INTEGER NN ! Numero de ventana INTEGER LSTR ! Longitud del texto. INTEGER SP ! SP=1 pone un espacio REAL X, Y ! Posicion donde empieza el texto. REAL XBOX(4),YBOX(4) ! Coordenadas de la caja que enmarca el texto. REAL XX,YY CHARACTER*1 CH CHARACTER*(*) TEXT CHARACTER*11 TEXTOLD CHARACTER*11 TEXTBCK C TEXTOLD=TEXT LSTR=0 CALL PGSLW(1) CALL PGSCI(1) ! Borramos el texto viejo CALL PGSCH(1.) CALL PGPTXT (X,Y,0.0,0.0,TEXT) CALL PGSCI(4) C IF(SP.EQ.1) THEN LSTR=1 TEXT=' ' END IF C 10 IF(LSTR.GT.0) THEN CALL PGPTXT(X,Y,0.0,0.0,TEXT(1:LSTR)) CALL PGQTXT(X,Y,0.0,0.0,TEXT,XBOX,YBOX) XX=XBOX(4) YY=YBOX(4) ELSE XX=X YY=Y END IF C CALL PGBAND(0,0,XX,YY,XX,YY,CH) C CALL PGSCI(1) ! Borramos el texto viejo IF(LSTR.GT.0) CALL PGPTXT (X,Y,0.0,0.0,TEXT(1:LSTR)) CALL PGSCI(4) C IF(CH.EQ.CHAR(92)) CH = '*' ! Para evitar problemas con "\" y PGPLOT C IF(ICHAR(CH).EQ.8.) THEN ! Backspace o ctrl H borran ultimo caracter IF(LSTR.GT.0) THEN TEXT(LSTR:LSTR)= ' ' LSTR=LSTR-1 END IF ELSE IF (ICHAR(CH).EQ.127) THEN ! Del borra todo el texto TEXT(1:LSTR)=' ' LSTR=0 IF(NN.EQ.2) LSTR=1 ELSE IF (ICHAR(CH).LT.32) THEN IF(LSTR.GT.0) THEN IF(NN.EQ.2.AND.TEXT(3:3).EQ.' ') THEN CALL PGPTXT(X,Y,0.0,0.0,TEXTOLD) TEXTBCK=TEXTOLD CALL IOWIN1(' Please use format : 00 00 00 ') ELSE CALL PGPTXT(X,Y,0.0,0.0,TEXT(1:LSTR)) TEXTBCK=TEXT(1:LSTR) END IF ELSE CALL PGPTXT(X,Y,0.0,0.0,TEXTOLD) TEXTBCK=TEXTOLD END IF GOTO 20 ELSE IF (LSTR.LT.LEN(TEXT)) THEN ! Aniadimos un caracter al texto LSTR=LSTR+1 TEXT(LSTR:LSTR)=CH END IF C GOTO 10 C 20 RETURN END C C********************************************************************** INTEGER FUNCTION TESTN(CH) C Funcion que chequea si un caracter es un numero del 0 al 9. C CHARACTER*1 CH C IF(ICHAR(CH).GT.47.AND.ICHAR(CH).LT.58) THEN TESTN = 1 ELSE TESTN=0 END IF C RETURN END C C ********************************************************************** C SUBROUTINE FJ_TS0(IYEAR,IMONTH,IDAY,RHOUR,FJ,TS0) C C Calcula la fecha juliana y el Tiempo Sidereo en Greenwich a 0h UT a C partir de los datos de una fecha concreta. En la fecha juliana se C tiene en cuenta la fraccion del dia correspondiente a la hora. C Los calculos se realizan en doble precision para no perder informacion C por redondeo. C IMPLICIT NONE C INTEGER IYEAR,IMONTH,IDAY ! Entrada REAL*8 RHOUR ! Entrada REAL*8 TS0 ! Salida REAL*8 FJ ! Salida C INTEGER*4 FECHA REAL*8 M,A,B,Y,DT,DTS0 REAL*8 DYEAR,DMONTH,DDAY C DYEAR=DBLE(IYEAR) DMONTH=DBLE(IMONTH) DDAY=DBLE(IDAY) FECHA=IYEAR*10000+IMONTH*100+IDAY IF(FECHA.GE.15821015)THEN A=DINT(DYEAR/1.D2) B=2.D0-A+DINT(A/4.D0) ELSE B=0.D0 END IF IF(IYEAR.GE.0)THEN A=0.D0 ELSE A=-.75D0 END IF IF(IMONTH.GE.3)THEN Y=DYEAR M=DMONTH ELSE Y=DYEAR-1.D0 M=DMONTH+1.2D1 END IF FJ=DINT(365.25D0*Y+A)+DINT(30.6001D0*(M+1.D0))+DDAY+B+ > 1720994.5D0 DT=(FJ-2451545.D0)/36525.D0 DTS0=6.D0*3.6D3+41.D0*6.D1+50.54851D0+8640184.812866D0*DT > +.093104D0*DT*DT-6.2D-6*DT*DT*DT DTS0=DMOD(DTS0,8.64D4) DTS0=DTS0/8.64D4 DTS0=DTS0*2.4D1 TS0=REAL(DTS0) IF(TS0.LT.0.)TS0=TS0+24. IF(TS0.GT.24.)TS0=TS0-24. FJ=FJ+DBLE(RHOUR/24.) END C C ********************************************************************** C SUBROUTINE POSSOL(FJ,ARSUN,DECSUN) C Calculamos AR y DEC del Sol para la FJ dada C Astronomical Formulae for Calculators C J. Meeus Cap. 18, pag. 79 y siguientes C IMPLICIT NONE C REAL*8 FJ ! Entrada REAL*8 ARSUN,DECSUN ! Salida REAL*8 PI,RAG REAL*8 ALAMBDA,PAR REAL*8 AR,DEC C PI=3.141592654D0 RAG=180.0D0/PI CALL PSOL(FJ,ALAMBDA,PAR) CALL ECUAT(FJ,ALAMBDA,0.D0,AR,DEC) ARSUN=AR*RAG/15.D0 DECSUN=DEC*RAG END C ******************************************************************** C SUBROUTINE PSOL(FJ,ALON,PAR) C C AL : Long. media C AM : Anomalia media C EXC: Excentricidad C EC : Ecuacion del centro C V : Anomalia verdadera C R : Radiovector C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PI=3.141592654D0 RAG=180.0D0/PI DPI=2.0D0*PI T=(FJ-2415020.D0)/36525.D0 C AL=279.69668D0+36000.76892D0*T+.0003025D0*T*T AM=358.47583D0+35999.04975D0*T-0.000150D0*T*T > +.0000033D0*T*T*T AL=DMOD(AL,360.D0)/RAG AM=DMOD(AM,360.D0)/RAG EXC=.016751014D0-.0000418D0*T-.000000126D0*T*T C1=1.919460D0-.004789D0*T-.000014D0*T*T C2=.020094D0-.000100D0*T EC=(C1*DSIN(AM)+C2*DSIN(2.*AM)+.000293D0*DSIN(3.*AM))/RAG C ALON=AL+EC V=AM+EC R=1.0000002D0*(1.D0-EXC*EXC)/(1.D0+EXC*DCOS(V)) C C Perturbaciones A=(153.23D0+22518.7541D0*T)/RAG B=(216.57D0+45037.5082D0*T)/RAG C=(312.69D0+32964.3577D0*T)/RAG D=(350.74D0+445267.1142D0*T-.00144D0*T*T)/RAG E=(231.19D0+20.20D0*T)/RAG H=(353.40D0+65928.7155D0*T)/RAG C A=DMOD(A,DPI) B=DMOD(B,DPI) C=DMOD(C,DPI) D=DMOD(D,DPI) H=DMOD(H,DPI) E=DMOD(E,DPI) C ALON=ALON+(.00134D0*DCOS(A)+.00154D0*DCOS(B) > +.00200D0*DCOS(C)+.00179D0*DSIN(D) > +.00178D0*DSIN(E))/RAG R=R+.00000543D0*DSIN(A)+.00001575D0*DSIN(B) > +.00001627D0*DSIN(C)+.00003076D0*DCOS(D) > +.00000927D0*DSIN(H) C c IF(IAP .EQ. 1) THEN ! Por si se necesitan c OM=(259.18D0-1934.142D0*T)/RAG ! coordenadas aparentes. c OM=DMOD(OM,DPI) c IF(OM .LT. 0.) OM=OM+DPI c ALON=ALON-(0.00569D0+0.00479D0*DSIN(OM))/RAG c ENDIF RR=R*1.49597870D11/6.37814D6 PAR=ASIN(1.D0/RR)*RAG C RETURN END C C C C ******************************************************************** C Calculamos Ascension Recta y Declinacion C SUBROUTINE ECUAT(FJ,AL,B,A,D) C IMPLICIT DOUBLE PRECISION(A-H,O-Z) C DPI=6.283185307D0 CALL EPSIL(FJ,EPSI) C X0=DCOS(B)*DCOS(AL) Y0=DCOS(B)*DSIN(AL) Z0=DSIN(B) C X=X0 Y=Y0*DCOS(EPSI)-Z0*DSIN(EPSI) Z=Y0*DSIN(EPSI)+Z0*DCOS(EPSI) C A=ATAN2(Y,X) IF(A .LT. 0.) A=A+DPI D=ASIN(Z) C RETURN END C********************************************************************** C Calculo de la oblicuidad de la ecliptica C SUBROUTINE EPSIL(FJ,EPSI) C IMPLICIT DOUBLE PRECISION(A-H,O-Z) C U=(FJ-2451545.D0)/3652500.D0 C A1=2.18D0-3375.70D0*U+0.36D0*U**2 A2=3.51D0+125666.39D0*U+0.10D0*U**2 EPSI=0.4090928D0-0.0226938D0*U-75.D-7*U**2+96926.D-7*U**3 > -2491.D-7*U**4-12104.D-7*U**5+(446*DCOS(A1)+28*DCOS(A2))*1.D-7 RETURN END C C ********************************************************************** C Calculamos para un objeto (sol) el orto y ocaso, condicion en la C cual obtiene una distancia cenital igual a DISZEN.Si el objeto C es circumpolar estado=C, si es invisible estado=I y sino estado=O C SUBROUTINE ORCASO(DECSUN,DISZEN,ESTADO) C IMPLICIT NONE C REAL*8 DECSUN,DISZEN ! Entrada REAL*8 LAT,LONG ! Entrada REAL*8 HORTO,HOCASO ! Salida REAL*8 AORTO,AOCASO ! Salida C REAL*8 PI PARAMETER(PI=3.141592654D0) REAL*8 LATR,DECR REAL*8 COSENOH,COSENOA REAL*8 ALTR CHARACTER*1 ESTADO C COMMON/BLKOBS/LAT,LONG COMMON/BLKOOC/HORTO,HOCASO,AORTO,AOCASO C IF(LAT.GE.0)THEN IF(DECSUN.GE.90.-LAT)THEN ESTADO='C' ELSE IF(DECSUN.LE.LAT-90.)THEN ESTADO='I' ELSE ESTADO='O' END IF ELSE IF(DECSUN.LE.-(90.+LAT))THEN ESTADO='C' ELSE IF(DECSUN.GE.90.+LAT)THEN ESTADO='I' ELSE ESTADO='O' END IF END IF C LATR=LAT*PI/180.D0 DECR=DECSUN*PI/180.D0 ALTR=(90.D0-DISZEN)*PI/180.D0 COSENOH=(SIN(ALTR)-SIN(LATR)*SIN(DECR))/(COS(LATR)*COS(DECR)) IF(COSENOH.LT.-1.) COSENOH=-1. IF(COSENOH.GT.+1.) COSENOH=+1. COSENOA=(SIN(ALTR)*SIN(LATR)-SIN(DECR))/(COS(ALTR)*COS(LATR)) IF(COSENOA.LT.-1.) COSENOA=-1. IF(COSENOA.GT.+1.) COSENOA=+1. HOCASO=ACOS(COSENOH) HOCASO=HOCASO*180.D0/PI AOCASO=ACOS(COSENOA) AOCASO=AOCASO*180.D0/PI HORTO=360.D0-HOCASO AORTO=360.D0-AOCASO HORTO=HORTO/15.D0 HOCASO=HOCASO/15.D0 C END C ********************************************************************** C SUBROUTINE TIEMPOU(ARSUN,H,LONG,TS0,TSL,TU) C C Calculo del Tiempo Sidereo Local y del Tiempo Universal. C IMPLICIT NONE C REAL*8 ARSUN ! Entrada REAL*8 H,LONG,TS0 ! Entrada REAL*8 TSL,TU ! Salida C REAL*8 DTS,TSG C TSL=ARSUN+H IF(TSL.GT.24.)TSL=TSL-24. TSG=TSL-LONG IF(TSG.GT.24.)TSG=TSG-24. IF(TSG.LT.0.)TSG=TSG+24. DTS=TSG-TS0 IF(DTS.LT.0.)DTS=DTS+24. TU=.997270*DTS IF(TU.LT.0.)TU=TU+24. C END C C ********************************************************************** C Pasamos un valor escrito como character formato +000 00 00 (II=1) C o +00 00 00 (II=2) a decimal DD.dddd C SUBROUTINE CH2DEC(CH11,DD,II) C IMPLICIT NONE INTEGER II INTEGER A1,A2,A3 REAL*8 DD CHARACTER*11 CH11 C IF(II.EQ.1) THEN IF(CH11(2:4).EQ.' ') THEN A1=0 ELSE READ(CH11(2:4),*) A1 ! Grados END IF IF(CH11(6:7).EQ.' ') THEN A2=0 ELSE READ(CH11(6:7),*) A2 ! Minutos END IF IF(CH11(9:10).EQ.' ') THEN A3=0 ELSE READ(CH11(9:10),*) A3 ! Segundos END IF ELSE IF (II.EQ.2) THEN IF(CH11(2:3).EQ.' ') THEN A1=0 ELSE READ(CH11(2:3),*) A1 ! Grados END IF IF(CH11(5:6).EQ.' ') THEN A2=0 ELSE READ(CH11(5:6),*) A2 ! Minutos END IF IF(CH11(8:9).EQ.' ') THEN A3=0 ELSE READ(CH11(8:9),*) A3 ! Segundos END IF END IF DD=DBLE(A1)+DBLE(A2)/60.D0+DBLE(A3)/3600.D0 IF(CH11(1:1).EQ.'-') DD=-DD END C C ********************************************************************** C Pasamos la hora de formato XX.xxxx (real) a XXh XXm (character) C SUBROUTINE R2CH(RR,CH) C INTEGER RR1,RR2 REAL*8 RR CHARACTER*7 CH C RR1=INT(RR) RR2=(RR-RR1)*60. WRITE(CH(1:2),'(I2)')RR1 WRITE(CH(3:4),'(A1)')'h ' WRITE(CH(5:6),'(I2)')RR2 WRITE(CH(7:7),'(A1)')'m' C END C C ********************************************************************** C Pasamos coordenadas de formato XX.xxxx (real) a +gg mm ss (character) C SUBROUTINE RR2CH(RR,CH) C INTEGER RR1,RR2,RR3 REAL*8 RR CHARACTER*9 CH C RR1=INT(RR) RR2=INT(ABS((RR-RR1)*60.)) RR3=INT((ABS((RR-RR1)*60.)-RR2)*60.) WRITE(CH(1:3),'(I3)')RR1 WRITE(CH(4:4),'(A1)') ' ' WRITE(CH(5:6),'(I2)')RR2 WRITE(CH(7:7),'(A1)')' ' WRITE(CH(8:9),'(I2)')RR3 C END C C********************************************************************** C Pasamos de coordenadas ecuatoriales a horizontales C SUBROUTINE ECU2HOR(TSL,LAT,AR,DEC,AZI,ALT) C IMPLICIT NONE REAL*8 TSL,LAT !Entrada REAL*8 AR,DEC ! Entrada REAL*8 HOR REAL*8 AZI,ALT ! Salida REAL*8 PI,G2R,R2G PARAMETER(PI=3.141592654D0) REAL*8 SENOAZ,COSEAZ REAL*8 HORR,DECR,LATR ! Angulos en radianes C G2R=PI/180.D0 R2G=180.D0/PI HOR=TSL-AR ! Angulo horario (en horas) HORR=HOR*G2R*15.D0 DECR=DEC*G2R LATR=LAT*G2R SENOAZ=COS(DECR)*SIN(HORR) COSEAZ=-(SIN(DECR)*COS(LATR))+COS(DECR)*SIN(LATR)*COS(HORR) AZI=ATAN2(SENOAZ,COSEAZ) AZI=AZI*R2G IF(AZI.LT.0.)AZI=AZI+360.D0 ALT=ASIN(SIN(LATR)*SIN(DECR)+COS(LATR)*COS(DECR)*COS(HORR)) ALT=ALT*R2G C END C C********************************************************************** C Cambio de coordenadas ecuatoriales debido a la precesion de los C equinocios. SUBROUTINE CHANEQUI(TII,TFF,ARII,ARFF,DECII,DECFF) C REAL*8 TII,TFF,TI,TF REAL*8 ARII,ARFF,DECII,DECFF REAL*8 X(3),X0(3) REAL*8 M(3,3) REAL*8 GIO,ZETA,TETA REAL*8 PI PARAMETER(PI=3.141592654D0) C TI=(TII-2000.)/100. TF=(TFF-2000.-100.*TI)/100. GIO=((2306.2181D0+1.39656D0*TI-.000139D0*TI*TI)*TF+ > (.30188D0-.000344D0*TI)*TF*TF+.017998D0*TF*TF*TF)/3600.D0 ZETA=GIO+((.7928+.00041*TI)*TF*TF+.000205*TF*TF*TF)/3600.D0 TETA=((2004.3109D0-.8533D0*TI-.000217D0*TI*TI)*TF- > (.42665D0+.000217D0*TI)*TF*TF-.041833*TF*TF*TF)/3600.D0 C GIO=GIO*PI/180.D0 ZETA=ZETA*PI/180.D0 TETA=TETA*PI/180.D0 C M(1,1)=-SIN(GIO)*SIN(ZETA)+COS(GIO)*COS(TETA)*COS(ZETA) M(1,2)=-COS(GIO)*SIN(ZETA)-SIN(GIO)*COS(ZETA)*COS(TETA) M(1,3)=-SIN(TETA)*COS(ZETA) M(2,1)=SIN(GIO)*COS(ZETA)+COS(GIO)*COS(TETA)*SIN(ZETA) M(2,2)=COS(GIO)*COS(ZETA)-SIN(GIO)*COS(TETA)*SIN(ZETA) M(2,3)=-SIN(TETA)*SIN(ZETA) M(3,1)=COS(GIO)*SIN(TETA) M(3,2)=-SIN(GIO)*SIN(TETA) M(3,3)=COS(TETA) C X0(1)=COS(DECII*PI/180.D0)*COS(ARII*15.*PI/180.D0) X0(2)=COS(DECII*PI/180.D0)*SIN(ARII*15.*PI/180.D0) X0(3)=SIN(DECII*PI/180.D0) C DO 505,K=1,3 X(K)=0. DO 500,KK=1,3 X(K)=X(K)+X0(KK)*M(K,KK) 500 CONTINUE 505 CONTINUE ARFF=ATAN2(X(2),X(1)) DECFF=ASIN(X(3)) ARFF=ARFF*180.D0/PI/15.D0 IF(ARFF.LT.0.)ARFF=ARFF+24.D0 DECFF=DECFF*180.D0/PI END C C********************************************************************** C SUBROUTINE CULMINA(LONG,LAT,TS0,AR,DEC,ALTCUL,DIRCUL,TUCUL) C C Se calcula la altura, direccion y tiempo universal de la culminacion C de una estrella a partir de sus AR. y DEC y la latitud del observador C IMPLICIT NONE C REAL*8 LONG,LAT REAL*8 AR,DEC REAL*8 ALTCUL,TUCUL REAL*8 TS0,TSL C CHARACTER*1 DIRCUL C C Comprobamos que el objeto sea visible IF(LAT.GE.0.) THEN IF(DEC.LE.LAT-90.) THEN DIRCUL='I' RETURN END IF ELSE IF(DEC.GE.90.+LAT) THEN DIRCUL='I' RETURN END IF END IF C Calculamos la altura y direccion de la culminacion IF(LAT.GE.0)THEN IF(DEC.LE.LAT)THEN ALTCUL=90.-LAT+DEC DIRCUL='S' ELSE ALTCUL=90.-DEC+LAT DIRCUL='N' END IF ELSE IF(DEC.LE.LAT)THEN ALTCUL=90.-ABS(DEC-LAT) DIRCUL='S' ELSE ALTCUL=90.-ABS(LAT-DEC) DIRCUL='N' END IF END IF C Calculamos T.U. del momento de la ocultacion CALL TIEMPOU(AR,0.D0,LONG,TS0,TSL,TUCUL) C END C C********************************************************************** C SUBROUTINE IOWIN1(IOTXT) C C Escribimos el texto IOTXT en la ventanita de I/O y ayuda C CHARACTER*40 IOTXT C CALL PGSVP(0.0,1.0,0.0,1.0) CALL PGSWIN(0.0,2.0,0.0,1.0) CALL PGSLW(1) CALL PGSCI(16) CALL PGSFS(1) CALL PGRECT(0.60,1.50,0.00,0.06) CALL PGSCI(1) CALL PGSCF(3) CALL PGTEXT(0.61,0.02,IOTXT) CALL PGSCF(2) CALL PGSFS(2) C END C********************************************************************** C SUBROUTINE IOWIN2(IOTXT1,IOTXT2) C C Escribimos el texto IOTXT en la ventanita de I/O y ayuda C CHARACTER*50 IOTXT1,IOTXT2 C CALL PGSVP(0.0,1.0,0.0,1.0) CALL PGSWIN(0.0,2.0,0.0,1.0) CALL PGSLW(1) CALL PGSCI(16) CALL PGSFS(1) CALL PGRECT(0.60,1.50,0.00,0.06) CALL PGSCI(1) CALL PGSCF(3) CALL PGSCH(.8) CALL PGTEXT(0.62,0.035,IOTXT1) CALL PGTEXT(0.62,0.010,IOTXT2) CALL PGSCF(2) CALL PGSFS(2) CALL PGSCH(1.) C END C********************************************************************** C SUBROUTINE MOONPOS(FJ,IYEAR,DAYS,MALFA,MDELTA,FASE) C C Calculamos la posicion de la Luna siguiendo los pasos que se dan C 'Practical Astronoy with your calculator' Peter Duffet-Smith, C cap.40 pg 91 y sgtes. C Tenemos como entrada un momento dado por la fecha juliana y como C salida las coord. ecuatoriales de la Luna en ese momento. C Elegimos como epoca de partida 1975.0 en la que : C Longitud de la Luna MLONG = 124.8756 grados C Longitud del perigeo MPER = 145.9601 " C Longitud del nodo MLONGN = 248.6441 " C Inclinacion de orbita de Luna MOINC = 5.1453 " C Excentricidad orbita Luna MOEXC = 0.0549 C Tamanio angular medio de Luna MSIZE = 0.5181 " C Semieje mayor de orbita Luna MSEJE = 384401 Km C Paralaje a su distancia media MPAR = 0.9507 grados C Longitud ecliptica del Sol SELONG = 279.041470 " C Longitud ecliptica de su perigeo SPELONG = 282.510396 " C Excentricidad de su orbita SOEXC = 0.016720 C IMPLICIT NONE C INTEGER I INTEGER IYEAR C REAL*8 MLONG,MPER,MLONGN,MOINC,MOEXC,MSIZE,MSEJE,MPAR REAL*8 DAYS,PI,G2R,R2G REAL*8 SELONG,SPELONG,SOEXC REAL*8 NN,EC REAL*8 SUNLONG ! Longitud ecliptica del Sol en un momento dado REAL*8 SUNANOM ! Anomalia media del Sol " " " " REAL*8 MOONLONG ! Longitud ecliptica de la Luna en un momento dado REAL*8 MOONLC,MOONLV ! " " " corregida Y verdadera REAL*8 MOONANOM ! Anomalia media " " " " " " REAL*8 MOONAC ! Anomalia media Luna, corregida REAL*8 MOONLN ! Longitud media del nodo ascendente REAL*8 MOONLNC ! " " " " corregida REAL*8 EVEC ! Evection REAL*8 ANECU ! Ecuacion anual REAL*8 VARIA ! Variacion REAL*8 ECC ! Ecuacion del centro REAL*8 CCC REAL*8 MLANDA,MBETA ! Latitud y longitud eclipticas de la Luna REAL*8 FJ,MALFA,MDELTA REAL*8 KK1,KK2 REAL*8 FASE ! Fraccion iluminada de cara Luna que vemos REAL*8 DANG ! Angulo Tierra-Luna-Sol C C Inicializamos las variables PI=3.141592654D0 G2R=PI/180.D0 R2G=180.D0/PI MLONG = 124.8756D0 MPER = 145.9601D0 MLONGN = 248.6441D0 MOINC = 5.1453D0 MOEXC = 0.0549D0 MSIZE = 0.5181D0 MSEJE = 384401D0 MPAR = 0.9507D0 SELONG = 279.041470D0 SPELONG = 282.510396D0 SOEXC = 0.016720D0 C Sumamos los dias de los anios desde 1975 hasta el actual. DO I=1975,IYEAR-1 IF((I/4.-INT(I/4.)).EQ.0) THEN DAYS=DAYS+366.D0 ELSE DAYS=DAYS+365.D0 END IF END DO C NN=360.D0*DAYS/365.25D0 30 IF(NN.GT.360.D0) THEN NN=NN-360.D0 GOTO 30 END IF C Calculamos la anomalia media del Sol SUNANOM=NN+SELONG-SPELONG IF(SUNANOM.LT.0.) SUNANOM=SUNANOM+360.D0 C Calculamos la longitud ecliptica del Sol EC=(360.D0/PI)*SOEXC*DSIN(SUNANOM*PI/180.D0) SUNLONG= NN+EC+SELONG IF(SUNLONG.GT.360.D0) SUNLONG=SUNLONG-360.D0 C Calculamos la longitud media de la Luna MOONLONG=360.D0*DAYS/27.3217D0+MLONG C Quitamos multiplos de 360, para que 0 > MOONLONG > 360 50 IF(MOONLONG.GT.360.D0) THEN MOONLONG=MOONLONG-360.D0 GOTO 50 ENDIF C Calculamos la anomalia media de la luna MOONANOM=MOONLONG-(360.D0/365.25D0)*(DAYS/8.85D0)-MPER C Aniadimos multiplos de 360, para que 0 > MOONANOM > 360 60 IF(MOONANOM.LT.0.) THEN MOONANOM=MOONANOM+360.D0 GOTO 60 ENDIF C Calculamos la longitud media del nodo ascendente MOONLN=MLONGN-(360.D0/365.25D0)*(DAYS/18.61D0) C Aniadimos multiplos de 360, para que 0 > MOONLN > 360 DO I=1,100 IF(MOONLN.LT.0.) THEN MOONLN=MOONLN+360.D0 ELSE GOTO 70 ENDIF END DO 70 CONTINUE C Calculamos la variacion de la excentricidad aparente de la orbita C de la Luna (evection). Esta se le aplicara a la longitud y anomalia C medias de la Luna EVEC=1.274D0*DSIN((2.*(MOONLONG-SUNLONG)-MOONANOM)*G2R) C Calculamos la variacion de la distancia Sol-Tierra a lo largo del C anio (ecuacion anual). Esta se le restara a la longitud y anomalia C medias de la Luna ANECU=0.18581D0*DSIN(SUNANOM*G2R) C Una tercera correccion que restar a la anomalia media de la Luna CCC=0.37D0*DSIN(SUNANOM*G2R) C Con estas tres correciones calculamos la anomalia media corregida C de la Luna MOONAC=MOONANOM+EVEC-ANECU-CCC C La ecuacion del centro sera: ECC=6.289D0*DSIN(MOONAC*G2R) C Calculamos la longitud corregida MOONLC=MOONLONG+EVEC-ANECU+ECC C Calculamos lo que afecta el que la Luna se encuentre a un lado u otro C de la Tierra, es decir mas cerca o mas lejos del Sol, y por tanto C la variacion de la atraccion gravitatorio del Sol sobre la Luna C (variacion) Esta se le restara a la anomalia media de la Luna. VARIA=0.65831D0*DSIN(2.D0*(MOONLC-SUNLONG)*G2R) C Se la sumamos a la longitud corregida y nos da la longitud verdadera MOONLV=MOONLC+VARIA C Con esta podemos calcular la longitud y latitud eclipticas. MOONLNC=MOONLN-0.16D0*DSIN(SUNANOM*G2R) C Para evitar ambiguidad en arcotangente, kk1 y kk2 deben ser muy parecidos KK2=MOONLV-MOONLNC KK1=DATAN(DTAN((KK2)*G2R)*DCOS(MOINC*G2R)) KK1=KK1*R2G 75 IF(ABS(KK2-KK1).GT.5.) THEN IF(KK2.GT.KK1) THEN KK1=KK1+180. ELSE KK1=KK1-180. END IF GOTO 75 END IF C MLANDA=KK1+MOONLNC MBETA=DASIN(DSIN(KK2*G2R)*DSIN(MOINC*G2R)) MBETA=MBETA*R2G C Pasamos a coord. ecuatoriales CALL ECLI2ECU(FJ,MLANDA,MBETA,MALFA,MDELTA) C Calculamos la fase de la Luna DANG=180.D0-SUNLONG+MOONLV FASE=0.5D0*(1.D0+COS(DANG*G2R)) C END c C********************************************************************** C SUBROUTINE NDAYS(IIYEAR,IIMONTH,IIDAY,IIHOUR,IIMIN,IISEC,DAYS) C C Calculamos el numero de dias y fracc. (DAYS) transcurridos del anio. C Se considera Enero 0.0 a la medianoche entre 30 y 31 Dec del C anio anterior. Asi mediodia del dia 1 de Enero sera 1.5 C IMPLICIT NONE C INTEGER IIYEAR,IIMONTH,IIDAY,IIHOUR,IIMIN,IISEC INTEGER IMONTH,IDAYS C REAL LEAP REAL*8 DAYS C LOGICAL LLEAP !TRUE -> Anio bisiesto C LLEAP=.FALSE. LEAP=(IIYEAR/4.-INT(IIYEAR/4.)) IF(LEAP.EQ.0.) LLEAP=.TRUE. IF(IIMONTH.GT.2) GOTO 80 IMONTH=IIMONTH-1 IF(LLEAP) THEN DAYS=IMONTH*62.D0 ELSE DAYS=IMONTH*63.D0 END IF DAYS=DAYS/2.D0 IDAYS=INT(DAYS) GOTO 120 80 IMONTH=IIMONTH+1 DAYS=IMONTH*30.6D0 IDAYS=INT(DAYS) IF(LLEAP) THEN IDAYS=IDAYS-62 ELSE IDAYS=IDAYS-63 END IF 120 DAYS=IDAYS+IIDAY+IIHOUR/24.D0+IIMIN/1440.D0+IISEC/86400.D0 C END C C********************************************************************** C Pasamos de coordenadas eclipticas a ecuatoriales C SUBROUTINE ECLI2ECU(FJ,ECLONG,ECLAT,ALFA,DELTA) C IMPLICIT NONE C REAL*8 ECLONG,ECLAT ! Latitud y longitud eclipticas (entrada) REAL*8 ALFA,DELTA ! Coordenadas ecuatoriales A.R. Y DEC (salida) REAL*8 FJ,EPSI REAL*8 G2R,R2G,PI C PI=3.141592654D0 G2R=PI/180.D0 R2G=180.D0/PI C Calculamos la oblicuidad de la ecliptica (epsil nos la da en radianes) CALL EPSIL(FJ,EPSI) C Calculamos la declinacion DELTA=DASIN(DSIN(ECLAT*G2R)*DCOS(EPSI) > +DCOS(ECLAT*G2R)*DSIN(EPSI)*DSIN(ECLONG*G2R)) DELTA=DELTA*R2G C Calculamos la ascension recta ALFA=DATAN(DTAN(ECLONG*G2R)*DCOS(EPSI) > -(DTAN(ECLAT*G2R)*DSIN(EPSI)/DCOS(ECLONG*G2R))) ALFA=ALFA*R2G C Para evitar ambiguedad al tomar el inverso de la tangente, nos C aseguramos de que ALFA y ECLONG esten en el mismo cuadrante+ IF(ECLONG.GT.90.) ALFA=ALFA+180. IF(ECLONG.GT.270.) ALFA=ALFA+180. C ALFA=ALFA/15.D0 C END C C********************************************************************** C Mostramos la ayuda correspondiente a la ventana que se haya elegido (NN) C SUBROUTINE HELP(NN) C IMPLICIT NONE C INTEGER NN C IF(NN.EQ.1) THEN CALL IOWIN2('Click left button to give the name of a new star ', > 'or right button for a file with a list of objects ') ELSEIF(NN.EQ.2) THEN CALL IOWIN2('Enter right ascension for a new star ', > 'with format : 00 00 00 ') ELSEIF(NN.EQ.3) THEN CALL IOWIN2('Enter declination for a new star ', > 'with format : +00 00 00 ') ELSEIF(NN.EQ.4) THEN CALL IOWIN2('Enter equinox for a new star ', > 'with format : 0000.00 ') ELSEIF(NN.EQ.5) THEN CALL IOWIN2('Draw the position of the star with the coordinates', > 'shown for the chosen observatory ') ELSEIF(NN.EQ.6) THEN CALL IOWIN2('Delete the traces of the stars and the Moon ', > 'currently drawed ') ELSEIF(NN.EQ.7) THEN CALL IOWIN2('Click left mouse button to change the name or the ', > 'right one to select another observatory ') ELSEIF(NN.EQ.8) THEN CALL IOWIN2('Enter the longitude of a new observatory ', > 'with format : +000 00 00 ') ELSEIF(NN.EQ.9) THEN CALL IOWIN2('Enter the latitude of a new observatory ', > 'with format : +00 00 00 ') ELSEIF(NN.EQ.10) THEN CALL IOWIN2('Enter the height of a new observatory ', > 'with format : 0000. ') ELSEIF(NN.EQ.11) THEN CALL IOWIN2('Enter a new date ', > 'with format : DD-MM-YYYY ') ELSEIF(NN.EQ.12) THEN CALL IOWIN2('Create a file "ves.lis" that list the drawed stars', > 'This can be used later as input-file ') ELSEIF(NN.EQ.13) THEN CALL IOWIN2('Create a postscript file with the current plot and', > 'data like coord., culmination, visibility, etc ...') ELSEIF(NN.EQ.14) THEN CALL IOWIN2('Yes, exactly. Exit = Quit = Bye ', > 'Adios, auf Wiedersehen .... ') ELSEIF(NN.EQ.15) THEN CALL IOWIN1(' Set Elevation for Y-axis ') ELSEIF(NN.EQ.16) THEN CALL IOWIN1(' Set sec-Z ~ airmass for Y-axis ') ELSEIF(NN.EQ.17) THEN CALL IOWIN1(' Set U.T. for X-axis ') ELSEIF(NN.EQ.18) THEN CALL IOWIN1(' Plot elevation or sec-Z vs azimut ') END IF C END