C @(#)pos1b.for	19.1 (ESO-IPG) 02/25/03 13:23:14
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
       PROGRAM POS1b
C---------------------------------------------------------------------------
C      THIS PROGRAM COMPUTES (ALFA,DELTA) FOR OBJECTS 
C      MEASURED ON PHOTOGRAPHIC PLATES.
C      
C---------------------------------------------------------------------------
C    Performs the RA,DEC <--> X,Y conversions. The transformation
C    parameters must have been previously computed by POS1A
C
C    Mode No. 
C      (1) STANDARD REDUCTION OF SCHMIDT PLATE, i.e. measurement
c          file contains the four crosses, as described above.
c          the coordinates of the plate center are given. 
c          after interactive reduction, the computed (alpha,
c          delta) values for a number of objects are printed
c          at the end. it is also possible to get the mean
c          (alpha,delta) of two eachother following positions 
c          (e.g. mp and comet trails).
c      (2) MEASUREMENT FILE DOES NOT CONTAIN THE FOUR CROSSES,
C          but otherwise as mode (1). this mode is the most common.
C      (3) CALCULATION OF (ALPHA,DELTA)1950 FOR SINGLE (X,Y)-VALUES,
C          which are input from terminal one by one.
C      (4) Inverse of mode (3). CALCULATION OF (X,Y)-VALUE FOR GIVEN
C          (ALPHA,DELTA)1950 VALUE, input from terminal. this part
c          uses a iteration procedure, and result is better than
c          1.0 micron. this mode is useful for finding objects
c          on a schmidt plate still on-place in measuring machine.
c      (5) Same as (4), except that the (A,D) INPUT IS TAKEN FROM 
C          A DISC FILE that contains a number of (a,d)-values for 
c          a series of objects. This file is created by program 
c          POS15 (for details, consult that program). This mode is
c          particularly useful for identification of moving objects 
c          within the plate limits. 
c      (6) REPEAT <POS1>
C      (7) END 
C 
C     Refer to pos1a.f for the history...
C #Version: Sat Sep 17 18:11:03 1994
C.VERSION   1995-Mar-26 : Use PARAMETER and increase buffer to 10001, PJG
c----------------------------------------------------------------------------
C
      IMPLICIT REAL*8 (A-H,O-Z), INTEGER*4 (I-N)
C 
      PARAMETER   (MAXOBJ=100001)
C 
      CHARACTER*1 YESNO,SCHMIDT,BLINK
      CHARACTER*2 ITEGN1,FLAGS
      CHARACTER*19 STRING_TERM
      CHARACTER*16 DFILE,ADXYFILE,XYFILE,OUTFILE
      CHARACTER*50 HEADER
      CHARACTER*16 ACUR, AFOR
      CHARACTER*250 MSG

      INTEGER TIDD, NCOLD,NROWD,NSC,ACOL,AROW,ISTAT,NCDIDENT
      INTEGER ILOOP, TIDO

      LOGICAL  NULL

      DIMENSION MNOST(MAXOBJ),NOST1(MAXOBJ)
      DIMENSION KK(9),X(10),AR(MAXOBJ,4),BX(9),BY(9),IDX(9),IDY(9) 
      DIMENSION A0(3),D0(3),XY(MAXOBJ,2),ALFA(3),DELTA(3),
     >SSX(10,10),S(10,10),R(10,10),SLSSC(9),SDC(9),PRXY(MAXOBJ,2), 
     >VBI(9),B(9),SV(10),SIGMA(10),XBAR(10)

      DIMENSION XTEMP(20),YTEMP(20),RESA(3),RESB(3) 
      DIMENSION IA1(2),ID1(2),KICK(100),XEST(2,2),YEST(2,2)


      EQUIVALENCE 
     >     (SSX,S,R),
     >     (S(2,1),B),
     >     (S(3,2),DEV),
     >     (S(4,2),SW,V),
     >     (S(5,2),SSR,RSW),
     >     (S(6,2),CONST),
     >     (S(7,2),SSDR,SUM),
     >     (S(8,2),RSQD),
     >     (S(9,2),CHG,ROFF),
     >     (S(4,3),XMSDR),
     >     (S(5,3),XMSDD),
     >     (S(6,3),VE),
     >     (S(7,3),SDE),
     >     (S(1,4),SLSSC),
     >     (S(1,5),TSQDC),
     >     (S(1,6),SDC),
     >     (S(1,7),VBI),
     >     (SV,SIGMA),
     >     (X,XBAR)
 
c---------------------------------------------------------------------
C---  MIDAS...
c---------------------------------------------------------------------
      INCLUDE          'MID_INCLUDE:ST_DEF.INC'
      COMMON           /VMR/MADRID
      INCLUDE          'MID_INCLUDE:ST_DAT.INC'
c---------------------------------------------------------------------
c--   init
      CALL STSPRO('POS1B')
c---------------------------------------------------------------------
      LU = 30
      DO0 = 0.0D0 
      DO1 = 1.0D0 
      DO12 = 1.2D1
      DO24 = 2.4D1
      DO60 = 6.0D1
      DO180 = 1.8D2 
      DO3600 = 3.6D3
      PI = 4.0D0*DATAN(1.0D0)
      DOPI18 = PI/DO180 
      FOC = 3.054D2/3.6D3
      DFILE = ' '
      ADXYFILE = ' '
      XYFILE = ' '
      OUTFILE = ' '
C 
      DO 5, J=1,100 
    5   KICK(J) = 0 
      NKICK = 0 

      ITRY = 0

c     initialize the ASCII log file
c      LU_LOG = 30
c      OPEN(UNIT=LU_LOG,FILE='pos1.log')

C     Option selected?
      CALL STKRDI('OPT',1,1,N1,IOTYP,IDUM1,IDUM2,ISTAT)

C     Serach the input table
      CALL STKRDC('IN_MES',1,1,20,N1,DFILE,IDUM1,IDUM2,ISTAT)
C     open the table      
      CALL TBTOPN(DFILE,F_I_MODE,TIDD,ISTAT)

c---------------------------------------------------------------------
c---  open the input files
c---------------------------------------------------------------------
C---  Reads the descriptors: 
*     Read the instrument flags:
      CALL STDRDC(TIDD,'FLAGS',1,1,2,N1,FLAGS,IDUM1,IDUM2,ISTAT)
      IF (FLAGS(1:1) .EQ. 'S' .OR. FLAGS(1:1) .EQ. 's' ) THEN
         SCHMIDT = 'Y'
      ELSE
         SCHMIDT = 'N'
      ENDIF
      IF (FLAGS(2:2) .EQ. 'y' .OR. FLAGS(2:2) .EQ. 'Y' 
     >     .OR. FLAGS(2:2) .EQ. 'B' .OR. FLAGS(2:2) .EQ. 'b' ) THEN
         BLINK = 'Y'
      ELSE
         BLINK = 'N'
      ENDIF

*     Read the plate and catalogue epochs:
      CALL STDRDD(TIDD,'EPO_PLA',1,1,N1,EPOCH,IDUM1,IDUM2,ISTAT)
      CALL STDRDD(TIDD,'EPO_CAT',1,1,N1,EP1,IDUM1,IDUM2,ISTAT)

*     Get the permitted permitted RMS
      CALL STDRDD(TIDD,'SIG',1,1,N1,SIG,IDUM1,IDUM2,ISTAT)

*     Get the xy transformation terms used:
      CALL STDRDC(TIDD,'TERMS',1,1,19,N1,STRING_TERM,IDUM1,IDUM2,ISTAT)
*     *decode the xy transformation terms flags
 150  NX = 0
      NY = 0
      DO I = 1,9
         READ(STRING_TERM(I:I),*)IDX(I)
         READ(STRING_TERM(I+10:I+10),*)IDY(I)
         IF (IDX(I).NE.0) NX=NX+1
         IF (IDY(I).NE.0) NY=NY+1
      ENDDO
  165 FORMAT(1X,'  CHECK ',9I2,1X,9I2) 

*     Reads the transformation parameter from the descriptor
      CALL STDRDD(TIDD,'BX',1,9,N1,BX,IDUM1,IDUM2,ISTAT)
      CALL STDRDD(TIDD,'BY',1,9,N1,BY,IDUM1,IDUM2,ISTAT)
      CALL STDRDD(TIDD,'CX',1,1,N1,CX,IDUM1,IDUM2,ISTAT)
      CALL STDRDD(TIDD,'CY',1,1,N1,CY,IDUM1,IDUM2,ISTAT)
      CALL STDRDI(TIDD,'NX',1,1,N1,NX,IDUM1,IDUM2,ISTAT)
      CALL STDRDI(TIDD,'NY',1,1,N1,NY,IDUM1,IDUM2,ISTAT)
      CALL STDRDD(TIDD,'AL_DE0',1,1,N1,ALFA0,IDUM1,IDUM2,ISTAT)
      CALL STDRDD(TIDD,'AL_DE0',2,1,N1,DELTA0,IDUM1,IDUM2,ISTAT)
      CALL STDRDD(TIDD,'AL_DE0',3,1,N1,COSD0,IDUM1,IDUM2,ISTAT)
      CALL STDRDD(TIDD,'AL_DE0',4,1,N1,SIND0,IDUM1,IDUM2,ISTAT)

    
c---------------------------------------------------------------------
C---  Reads the measurements (X,Y)
c---------------------------------------------------------------------
*     needs the measurements only for option 1&2 
*     (conversion of measured XY -> alpha,delta)
      IF (IOTYP.EQ.1 .OR. IOTYP.EQ.2 ) THEN

C        check the numbers of rows and columns of the table
         CALL TBIGET(TIDD,NCOLD,NROWD,NSC,ACOL,AROW,ISTAT)      

C        check if the columns exists
*        PPM identifier
         CALL TBCSER(TIDD,':IDENT',NCDIDENT,ISTAT)
         IF (NCDIDENT.LE.0) THEN
            CALL STTPUT('*error* Input column IDENT not found',ISTAT)
            STOP
         ENDIF

*        x position
         CALL TBCSER(TIDD,':XCEN',NCDX,ISTAT)
         IF (NCDX.LE.0) THEN
            CALL STTPUT('*error* Input column XCEN not found',ISTAT)
            STOP
         ENDIF

*        y position
         CALL TBCSER(TIDD,':YCEN',NCDY,ISTAT)
         IF (NCDY.LE.0) THEN
            CALL STTPUT('*error* Input column YCEN not found',ISTAT)
            STOP
         ENDIF



 55      DO I = 1,3
            RESA(I) = 0.0
            RESB(I) = 0.0
         ENDDO

         IN = 0
         ITOT = 0
         IRES = 0
         IRESNO = 0
         
         DO 62 ILOOP = 1, NROWD
            CALL TBERDC(TIDD,ILOOP,NCDIDENT,ACUR,NULL,ISTAT)

*           *(x/ycen are r*4 read as r*8)
            CALL TBERDD(TIDD,ILOOP,NCDX,XCUR,NULL,ISTAT)
            CALL TBERDD(TIDD,ILOOP,NCDY,YCUR,NULL,ISTAT)

            IF (BLINK.EQ.'Y'.OR.BLINK.EQ.'y') THEN
*           *blink funny format: conversion to (+,-) microns
               XBLINK = XCUR
               YBLINK = YCUR
               CALL BLINKXY(XBLINK,YBLINK,XCUR,YCUR,1)
            ENDIF

            MEND = 0
            IN = IN+ 1


C---        Compute mean (x,y) for each object and R.M.S. in X,Y,R
            IF (IN .LE. 1) THEN
C           *first (X,Y)
               AFOR = ACUR
               XTEMP(1) = XCUR
               YTEMP(1) = YCUR
               IX = 1
               GO TO 62
            ENDIF

            IF (ACUR .EQ. AFOR) THEN
C           *still same object
               IX = IX + 1
               XTEMP(IX) = XCUR
               YTEMP(IX) = YCUR
               GO TO 62
            ENDIF


c-          new object read, do means for the one before
            CALL MEAN(XTEMP,YTEMP,IX,XMEAN,YMEAN,NOUT,RESA,SIG)
            IF (IX.GT.2) THEN
               DO I = 1,3
                  RESB(I) = RESB(I) + RESA(I)
               ENDDO
               IRES = IRES + IX
               IRESNO = IRESNO + 1
            ENDIF
            XTEMP(1) = XCUR
            YTEMP(1) = YCUR


C           *Maximum # of objects is MAXOBJ
            ITOT = ITOT + 1
            IF (ITOT.GT.MAXOBJ) THEN
               WRITE(MSG,9074)
               CALL STTPUT(MSG,ISTAT)
 9074          FORMAT('*error* too many objects ! ')
               CALL STSEPI
               STOP
            ENDIF
 
            PRXY(ITOT,1) = XMEAN/1.0D3
            PRXY(ITOT,2) = YMEAN/1.0D3

*           *conversion Text ident --> Nr ident
            IID = 0
 628        IF (IID .LT. 16) THEN
               IID = IID+1
            ELSE
               WRITE(MSG,*) '*warning* identiers must contain numbers'
               CALL STTPUT(MSG,ISTAT)
               WRITE(MSG,*) '          identifier set to 0'
               CALL STTPUT(MSG,ISTAT)
               NFOR = 0
               GOTO 629
            ENDIF
            READ(AFOR(IID:16),*,err=628) NFOR

 629        MNOST(ITOT) = NFOR
            AFOR = ACUR
            IX = 1

*           *go to read next line
 62      CONTINUE

*        *file is finished
 63      NOBJ = ITOT
         WRITE(MSG,9080) NOBJ
         CALL STTPUT(MSG,ISTAT)
 9080    FORMAT('   No. of objects in (X,Y)-file: ',I5)


C        CALCULATE MEAN R.M.S. FOR MEASURED (X,Y) IN MICRONS

         IF (IRES.GT.0) THEN
            DO I = 1,3
               RESB(I) = RESB(I)/FLOAT(IRESNO)
            ENDDO
            WRITE(MSG,9082) RESB
            CALL STTPUT(MSG,ISTAT)
 9082       FORMAT('   Mean R.M.S. (Microns) in (X,Y,R) :',3f8.2)
            WRITE(MSG,9182) IRES,IRESNO
            CALL STTPUT(MSG,ISTAT)
 9182       FORMAT('   (pix or microns) based on ',i3,
     $           ' (X,Y) measurements of ',
     $           i3,' objects ')
         ELSE
            WRITE(MSG,9083)
 9083       FORMAT('   Not enough measurements of same object to ',
     >           ' calculate mean r.m.s. values ')
         ENDIF
      ENDIF


C----------------------------------------------------------------------
C---  Start the real work
C----------------------------------------------------------------------
C      write(6,*) 'start the real work'
C      ERX = SNGL(SUMD1) 
C      ERY = SNGL(SUMD2) 


c mode:  IOTYP 
c (1) Calculate (A,D) for plate center and several (X,Y)-values'/ 
c (2) Calculate (A,D) for several (X,Y)-values only'/ 
c (3) Calculate (A,D) for one (X,Y)-value'/ 
c (4) Calculate (X,Y) for one (A,D)-position'/
c (5) Calculate (X,Y) for several (A,D)s in disc file'/
c (6) NEW RUN OF POS1'/ 
c (7) EXIT PROGRAMME'// 

      IF (IOTYP.EQ.1) GO TO 1040 
      IF (IOTYP.EQ.2) GO TO 1040
      IF (IOTYP.EQ.3) GO TO 1070
      IF (IOTYP.EQ.4) GO TO 1500 
      IF (IOTYP.EQ.5) GO TO 1600 
      IF (IOTYP.GE.6) GO TO 2000
      GO TO 1010

C-----ENTRY FOR MODE 1
C-----ENTRY FOR MODE 2
C     * Mean positions (e.g. of trail ends)? 
 1040 CALL STKRDC('FLAG',1,1,1,N1,YESNO,IDUM1,IDUM2,ISTAT)

      LMEAN = 0 
      IF (YESNO.EQ.'Y'.OR.YESNO.EQ.'y') LMEAN = 1
      IF (IOTYP-1) 1120,1100,1120 

C-----ENTRY FOR MODE 3
 1070 IF (BLINK.EQ.'Y'.OR.BLINK.EQ.'y')  THEN
         WRITE(MSG,1075)
         CALL STTPUT(MSG,ISTAT)
 1075    FORMAT('  Use *BLINK* Coordinates. ',
     $        'Give identification No., X, Y: ')
         READ(*,*) NOST1(1),XBLINK,YBLINK
         CALL BLINKXY(XBLINK,YBLINK,XCUR,YCUR,1)
      ELSE
         WRITE(MSG,1080) 
         CALL STTPUT(MSG,ISTAT)
 1080    FORMAT('  Give identification Number, X, Y: ') 
         READ(*,*) NOST1(1),XCUR,YCUR
      ENDIF
      XY(1,1) = XCUR/1.0D3 
      XY(1,2) = YCUR/1.0D3 
      PRXY(1,1) = XY(1,1)
      PRXY(1,2) = XY(1,2)
      IST = 0 
      IEX = 1 
      GO TO 1270 

C---- entry
 1100 WRITE(MSG,9110)
      CALL STTPUT(MSG,ISTAT)
      WRITE(MSG,9111)
      CALL STTPUT(MSG,ISTAT)
      WRITE(MSG,9112)
      CALL STTPUT(MSG,ISTAT)
 9110 FORMAT(' *info* Standard Schmidt plate reduction: ')
 9111 FORMAT('        First 4 non-standard stars(X,Y) are supposed ')
 9112 FORMAT('        to be those of the four crosses.')
C---- entry
 1120 IST = 0 
*     Read the output file name
      CALL STKRDC('OUT_MES',1,1,16,N1,OUTFILE,IDUM1,IDUM2,ISTAT)
*     create the output table
      CALL TBTINI(OUTFILE,F_TRANS,F_O_MODE,6,1000,TIDO,ISTAT)
*     create the columns
      CALL TBCINI(TIDO,D_I4_FORMAT,1,'I7'   ,' ','PPM',NSCPPM,ISTAT)
      CALL TBCINI(TIDO,D_R8_FORMAT,1,'F5.2' ,' ','MAG',NSCMAG,ISTAT)
      CALL TBCINI(TIDO,D_R8_FORMAT,1,'R12.7',' ','R_A',NSCR_A,ISTAT)
      CALL TBCINI(TIDO,D_R8_FORMAT,1,'S12.6',' ','DEC',NSCDEC,ISTAT)
      CALL TBCINI(TIDO,D_R8_FORMAT,1,'F8.4' ,' ','PMA',NSCPMA,ISTAT)
      CALL TBCINI(TIDO,D_R8_FORMAT,1,'F8.4' ,' ','PMD',NSCPMD,ISTAT)


 1200 DO 1240 I = 1,NOBJ
            XY(I,1) = PRXY(I,1) 
            XY(I,2) = PRXY(I,2) 
            NOST1(I) = MNOST(I) 
 1240    CONTINUE
      IEX = NOBJ
 1250 WRITE(MSG,9255) IEX
      WRITE(MSG,*) IEX, nobj
      CALL STTPUT(MSG,ISTAT)
 9255 FORMAT(2X,I5,'    (X,Y)-values read ')

      IF (IOTYP.EQ.1) THEN 
*        plate center
         SUMD1=(XY(3,1)-XY(1,1))*(XY(2,2)-XY(4,2))-
     1         (XY(3,2)-XY(1,2))*(XY(2,1)-XY(4,1)) 
         SUMD2=(XY(2,1)-XY(1,1))*(XY(2,2)-XY(4,2))-
     1         (XY(2,2)-XY(1,2))*(XY(2,1)-XY(4,1)) 
         SUMD3=(XY(3,1)-XY(1,1))*(XY(2,2)-XY(1,2))-
     1         (XY(3,2)-XY(1,2))*(XY(2,1)-XY(1,1)) 
         XCEN = XY(1,1)+(SUMD2/SUMD1)*(XY(3,1)-XY(1,1)) 
         YCEN = XY(2,2)+(SUMD3/SUMD1)*(XY(4,2)-XY(2,2)) 
         NCEN = 0 
         XXCEN = XCEN*1.0D3
         YYCEN = YCEN*1.0D3
         WRITE(MSG,9260) XXCEN,YYCEN 
         CALL STTPUT(MSG,ISTAT)
 9260    FORMAT('   Plate center: X = ',F7.3,'  Y = ',F7.3,'  microns') 
         IEX = IEX + 1
         XY(IEX,1) = XCEN
         XY(IEX,2) = YCEN
      ENDIF

*---entry
 1270 DO 1280 I = 1,IEX
         AR(I,1) = PRXY(I,1)
         AR(I,2) = PRXY(I,2)
         AR(I,3) = DO0 
         AR(I,4) = DO0 
 1280    CONTINUE


*     Conversion of the coord
*     X:
      IXY = 1 
      DO 1300 I16= 1,IEX 
         N = NX+1
         CALL RAR(I16,0,N,KK,IXY,IDX,IDY,X,AR)
         V = DO0 
         DO 1290 I17 = 1,NX
 1290       V = V+BX(I17)*X(I17)
         XY(I16,IXY) = V + CX
 1300    CONTINUE

*     Y:
      IXY = 2 
      DO 1320 I16 = 1,IEX
         N = NY+1
         CALL RAR(I16,0,N,KK,IXY,IDX,IDY,X,AR)
         L=NY
         V=DO0 
         DO 1310 I17=1,L
 1310       V=V+BY(I17)*X(I17)
         XY(I16,IXY) = V + CY
 1320    CONTINUE




      DO 1440 I = 1,IEX
         AKSI = XY(I,1)
         ANU = XY(I,2) 
         CALL XYAD(AKSI,ANU,IA1,ID1,DA1S,DD1S,ITEGN1,ALFAS,DELTS,
     *        COSD0,SIND0,ALFA0)

         IF (I.EQ.1) THEN
            CALL STTPUT(' ',ISTAT)
            CALL STTPUT(' ',ISTAT) 
            CALL STTPUT(' ',ISTAT)
            Write(MSG,9340)
            CALL STTPUT(MSG,ISTAT)
            WRITE(MSG,9341)
            CALL STTPUT(MSG,ISTAT)
            WRITE(MSG,9342)
            CALL STTPUT(MSG,ISTAT)
            WRITE(MSG,9343)
            CALL STTPUT(MSG,ISTAT)
            WRITE(MSG,9344)
            CALL STTPUT(MSG,ISTAT)
         ENDIF
 9340 FORMAT('   The *EQUINOX* of the computed coordinates is that of')
 9341 FORMAT('               the standard stars catalogue')
 9342 FORMAT('   Computed positions for measured objects:')
 9343 FORMAT('      No.      R.A.           Dec. ') 
 9344 FORMAT('    <----> <---------->   <---------->') 

C        IF((I.EQ.1.AND.IRUN.EQ.1).AND.IOTYP.EQ.1) THEN
         IF (I.EQ.1 .AND. IOTYP.EQ.1) THEN
            WRITE(MSG,9352) 
            CALL STTPUT(MSG,ISTAT)
         ENDIF
 9352    FORMAT('Crosses N(1),W(2),S(3),E(4)')

C        IF((I.EQ.1.AND.IRUN.EQ.1).AND.IOTYP.EQ.2) THEN
         IF (I.EQ.1 .AND. IOTYP.EQ.2) THEN
            WRITE(MSG,9360) 
            CALL STTPUT(MSG,ISTAT)
         ENDIF

C        IF((I.EQ.5.AND.IRUN.EQ.1).AND.IOTYP.EQ.1) THEN
         IF (I.EQ.5 .AND. IOTYP.EQ.1) THEN
            WRITE(MSG,9360) 
            CALL STTPUT(MSG,ISTAT)
         ENDIF
 9360    FORMAT('Other objects:')
         
         IF(I.EQ.IEX.AND.IOTYP.EQ.1) THEN
            WRITE(MSG,9370)
            CALL STTPUT(MSG,ISTAT)
         ENDIF
 9370    FORMAT('Plate center coordinates:')
         
         WRITE(MSG,9380) NOST1(I),IA1(1),IA1(2),DA1S,ITEGN1,ID1(1),
     >        ID1(2),DD1S
         CALL STTPUT(MSG,ISTAT)
 9380    FORMAT(I10,2I3,F7.3,2X,A2,I2,I3,F6.2)


*---     OPTION 3
         IF (IOTYP.EQ.3) THEN
*           *only one xy
            WRITE(MSG,9381) NOST1(I),IA1(1),IA1(2),DA1S,ITEGN1,
     >           ID1(1),ID1(2),DD1S 
            CALL STTPUT(MSG,ISTAT)
 9381       FORMAT(' (R.A., DEC) for object',I8,' :',2I3,F7.3,2X,
     >           A2,I2,I3,F6.2)
         ELSE
*        *NOT options 3:
*           *a whole set of xy
            IF (LMEAN.NE.1) THEN
*              *not  measurement of trails
               CALL TBEWRI(TIDO, I, NSCPPM, NOST1(I), ISTAT)
               ALFAS = ALFAS/DOPI18
               CALL TBEWRD(TIDO, I, NSCR_A, ALFAS, ISTAT)
               DELTS = DELTS/DOPI18
               CALL TBEWRD(TIDO, I, NSCDEC, DELTS, ISTAT)
            ELSE
*              *measurement of trails
               IF (MOD(I,2).NE.1)THEN 
                  STAL = (ALFAS + STAL)/2D0 
                  STDE = (DELTS + STDE)/2D0 
                  CALL ADCON(STAL,STDE,ALFA,DELTA)
                  
                  DO J = 1,2
                     IA1(J) = IDINT(ALFA(J) + 1.0D-3)
                     ID1(J) = IDINT(DELTA(J) + DSIGN(1.0D-3,DELTA(J))) 
                  ENDDO
                  
                  DA1S = ALFA(3)
                  DD1S = DELTA(3) 
                  ITEGN1 = ' +' 
                  IF (STDE.LT.DO0) ITEGN1 = ' -'
                  ID1(1) = IABS(ID1(1)) 
                  WRITE(MSG,1420) 
     >                 IA1(1),IA1(2),DA1S,ITEGN1,ID1(1),ID1(2),DD1S
                  CALL STTPUT(MSG,ISTAT)
 1420             FORMAT(20X,'Mean',1X,2I3,F7.3,1X,A2,I2,I3,F6.2) 
               ELSE                  
                  STAL = ALFAS
                  STDE = DELTS
               ENDIF
            ENDIF
         ENDIF       

 1440 CONTINUE

      GO TO 1010



C-----ENTRY FOR MODE 4
 1500 WRITE(MSG,1501)
      CALL STTPUT(MSG,ISTAT)
 1501 FORMAT('  Give (alpha(H,M,S),delta(DEG,M,S)): ')
      READ(5,*) A0,D0 

      CALL ADXY(A0,D0,XEST,YEST,NITER,INEX,COSD0,SIND0,ALFA0, 
     *     NX,NY,KK,IXY,IDX,IDY,BX,BY,CX,CY,SCHMIDT)


      IF (INEX.EQ.1) GO TO 1010
      IF (BLINK.EQ.'Y'.OR.BLINK.EQ.'y') THEN
         XCUR = XEST(2,1)
         YCUR = YEST(2,1)
         CALL BLINKXY(XBLINK,YBLINK,XCUR,YCUR,2)
         XEST(2,1) = XBLINK
         YEST(2,1) = YBLINK
         WRITE(MSG,1502)
         CALL STTPUT(MSG,ISTAT)
 1502    FORMAT(' NB! Blink coordinates !')
      ENDIF
      DO 1505 J = 1,2
         IA1(J) = IDINT(A0(J) + 1.0D-3)
 1505    ID1(J) = IDINT(D0(J) + DSIGN(1.0D-3,D0(J))) 
      WRITE(MSG,1510) NITER,XEST(2,1),YEST(2,1)
      CALL STTPUT(MSG,ISTAT)
 1510 FORMAT(I5,' iterations; ( X , Y ) = ',2F12.2) 
C     WRITE(MSG,1520) IA1,A0(3),ID2,D0(3),XEST(2,1),YEST(2,1)
      WRITE(MSG,1520) IA1,A0(3),ID1,D0(3),XEST(2,1),YEST(2,1) 
      CALL STTPUT(MSG,ISTAT)
 1520 FORMAT(1X,2I3,F7.3,2I4,F6.2,' : (X,Y) = ',2F12.2) 
      WRITE(MSG,1530)
      CALL STTPUT(MSG,ISTAT)
 1530 FORMAT(' More positions ? ') 
      READ(*,9016) YESNO
      IF (YESNO.EQ.'Y'.OR.YESNO.EQ.'y') GO TO 1500
      GO TO 1010




C-----ENTRY FOR MODE 5
 1600 WRITE(*,1630)
 1630 FORMAT('***SOOORRRY*** this part works on plain ASCII, not MIDAS'/
     $     '  GIVE NAME OF DISC FILE WITH (A,D)-VALUES: ',$)
      READ(*,1640) ADXYFILE
 1640 FORMAT(A16) 
      OPEN(UNIT=50,FILE=ADXYFILE,STATUS='old',ERR=1600)
      WRITE(*,1650)
 1650 FORMAT('  GIVE MAX.DISTANCE FROM PLATE CENTER IN MICRONS: ',$) 
      READ(*,*) DMAX
      IST = 0 
      I1 = 0
      WRITE(*,1670)
 1670 FORMAT('  STORE THE (X,Y)-VALUES IN A DISC FILE? ',$)
      READ(*,9016) YESNO
      IF (YESNO.EQ.'Y'.OR.YESNO.EQ.'y') IST = 1
      IF (IST.NE.1) GO TO 1710 
 1675 WRITE(*,1680)
 1680 FORMAT('  GIVE FILE NAME: ',$) 
      READ(*,1640) XYFILE
      OPEN(UNIT=60,FILE=XYFILE,STATUS='NEW',ERR=1675)
      ITOT = 0
 1710 READ(50,1720) HEADER
 1720 FORMAT(A50)
C      WRITE(LU_LOG,1730) ADXYFILE,HEADER
      WRITE(*,1730) ADXYFILE,HEADER
 1730 FORMAT(/,1X,'  (A,D)-POSITIONS FROM FILE: ',A16,
     1/,1X,'  IDENT: ',A50       // 
     21X,'    NO.       X            Y') 
 1740 READ(50,*,END=1010) IANO,A0,D0
 1760 CALL ADXY(A0,D0,XEST,YEST,NITER,INEX,COSD0,SIND0,ALFA0, 
     *NX,NY,KK,IXY,IDX,IDY,BX,BY,CX,CY,SCHMIDT)
      IF (BLINK.EQ.'Y'.OR.BLINK.EQ.'y') THEN
         XCUR = XEST(2,1)
         YCUR = YEST(2,1)
         CALL BLINKXY(XBLINK,YBLINK,XCUR,YCUR,2)
         XEST(2,1) = XBLINK
         YEST(2,1) = YBLINK
         WRITE(*,1502)
      ENDIF
      WRITE(*,1775) IANO,XEST(2,1),YEST(2,1)
C      WRITE(LU_LOG,1775) IANO,XEST(2,1),YEST(2,1) 
 1775 FORMAT(1X,I6,2F12.2) 
 1785 IF (IST.NE.1) GO TO 1740 
      WRITE(60,1775) IANO,XEST(2,1),YEST(2,1)
 1790 FORMAT(4F8.0) 
      GO TO 1740 

C---- Error messages
C 8018 WRITE(MSG,*) '*ERROR* cannot open standard file ',PSFILE
C      CALL STTPUT(MSG,ISTAT)
C      STOP
C
C 8019 WRITE(MSG,*) '*ERROR* cannot open positions file ',DFILE
C      CALL STTPUT(MSG,ISTAT)
C      STOP

C-----ENTRY FOR MODE 6
C-----ENTRY FOR MODE 7
 2000 CONTINUE 
 3000 WRITE(MSG,3001) 
      CALL STTPUT(MSG,ISTAT)
 3001 FORMAT('   ')
 1010 CALL STSEPI
 9045 FORMAT(A50)
 9050 FORMAT(2X,A50)
 9999 Format(80a)
 9076 FORMAT(I8,2F12.3,I5,3F8.3,I5)
 9085 FORMAT(I4,I4,I3,F7.3,I5,I4,F7.2,2F8.2)
 9115 FORMAT(I6,F6.2,2I3,F7.3,A2,I2,I3,F6.2,2I6)
 9016 FORMAT(A1)

      END

