C @(#)meanstar.for	19.1 (ES0-DMD) 02/25/03 13:29:08
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 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding 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 MEANSTAR
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT (c)
C.IDENT         meanstar.for
C.MODULE
C.AUTHOR        Andrew T. Young
C.KEYWORD
C.LANGUAGE      FORTRAN 77
C.PURPOSE       collapse star table with multiple positions to means.
C.COMMENTS
C.VERSION       4.3
C.RETURNS
C.ENVIRONMENT
C.
C-----------------------------------------------------------------------------
C*****************************************************************************
C
C
C       Reads a sorted star table with duplicate positions, and tries to
C       estimate a "best" position for each star.  Complains if the spread
C       in either coordinate exceeds a minute of arc.
C
C       This is done in the command CONVERT/PHOT.
C
C       This program is modified from the "esodcon" program, and may contain
C       fossils from it.
C
C
C*****************************************************************************
C
C
        IMPLICIT NONE
C
C  BEGIN Declarations:
C
C
        INTEGER ITBL, NEWTBL
        INTEGER NCOLS,NROWS,NSORTC,NWPRAL,NROWSAL, ISTAT
C
C
        CHARACTER CARD*78, TBLFIL*80
        CHARACTER*32 OBJECT,OLDOBJ
C
        INTEGER ITEM,ITEMS
        INTEGER KOBJ,KRA,KDEC,KEQ,NROW
        INTEGER JOBJ,JRA,JDEC,JEQ,NROUT
C
        INTEGER MXITEM
        PARAMETER (MXITEM=900)
        REAL RA,DEC,EQUINOX
        REAL RAS(MXITEM),DECS(MXITEM),EQUINOXS(MXITEM)
C
        LOGICAL NULL
C
C    Types for external fcns.:
C
C
C  Set up MIDAS declarations:
C
        INTEGER MADRID(1)
C
        COMMON /VMR/ MADRID
C
        INCLUDE 'MID_INCLUDE:ST_DEF.INC'
C
        INTEGER IUNIT
        INTEGER NACTEL,NULLS
C
C  END Declarations.
C
C
C  BEGIN DATA statements:
C
        INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
C
C  END DATA statements.
C
C
C  ********************  PROLOGUE  ********************
C
      CALL STSPRO ('MEANSTAR')
C
      OLDOBJ=' '
      NROW=0
      NROUT=0
C
C       Set up INPUT table file:
C
      CALL TV('Opening  sdata.tbl')
      CALL TBTOPN('sdata.tbl', 1,  ITBL, ISTAT)
      IF (ISTAT.NE.0) CALL TERROR (ITBL,1,'Could not open "sdata.tbl".')
C
      CALL TBIGET(ITBL,  NCOLS,NROWS,NSORTC,NWPRAL,NROWSAL, ISTAT)
      IF (ISTAT.NE.0) CALL TERROR
     1      (ITBL,1,'Could not get basic table data.')
C
C       Get column pointers...
C
      CALL TBLSER (ITBL, 'OBJECT',  KOBJ,ISTAT)
      IF (ISTAT.NE.0 .OR. KOBJ.EQ.-1)
     1    CALL TERROR(ITBL,2,'Could not find column OBJECT')
      CALL TBLSER (ITBL, 'RA',  KRA,ISTAT)
      IF (ISTAT.NE.0 .OR. KRA.EQ.-1)
     1    CALL TERROR(ITBL,3,'Could not find column RA')
      CALL TBLSER (ITBL, 'DEC',  KDEC,ISTAT)
      IF (ISTAT.NE.0 .OR. KDEC.EQ.-1)
     1    CALL TERROR(ITBL,4,'Could not find column DEC')
      CALL TBLSER (ITBL, 'EQUINOX',  KEQ,ISTAT)
      IF (ISTAT.NE.0 .OR. KEQ.EQ.-1)
     1    CALL TERROR(ITBL,5,'Could not find column EQUINOX')
C
C       Make sure table is sorted on OBJECT column:
C
      CALL TBCSRT (ITBL, 1, KOBJ, 1,  ISTAT)
      IF (ISTAT.NE.0)
     1        CALL TERROR(ITBL,8,'Could not sort OBJECT column')
C
C
C       Set up OUTPUT table file:
C
C     First, get name from tblfil local keyword:
      CALL STKRDC ('TBLFIL', 1, 1, 80,  NACTEL,TBLFIL,IUNIT,NULLS,ISTAT)
C
      CARD='Creating '//TBLFIL(:71)
      CALL TV(CARD)
      CALL TBTINI (TBLFIL, 0, 0, 1, 1,  NEWTBL, ISTAT)
      CARD='Could not create '//TBLFIL                      ! RHW 4/10/1993
      IF (ISTAT.NE.0) CALL TERROR (ITBL,11,CARD)
C
C       create column pointers...
C
      CALL TBCINI(NEWTBL, D_C_FORMAT,32,'A32',' ','OBJECT',JOBJ,ISTAT)
      IF (ISTAT.NE.0)
     1        CALL TERROR(NEWTBL,12,'Could not create OBJECT column')
      CALL TBCINI(NEWTBL, D_R4_FORMAT, 1,'R10.5',' ','RA',JRA,ISTAT)
      IF (ISTAT.NE.0)
     1        CALL TERROR(NEWTBL,13,'Could not create RA column')
      CALL TBCINI(NEWTBL, D_R4_FORMAT, 1,'s9.4',' ','DEC',JDEC,ISTAT)
      IF (ISTAT.NE.0)
     1        CALL TERROR(NEWTBL,14,'Could not create DEC column')
      CALL TBCINI
     1        (NEWTBL, D_R4_FORMAT, 1,'F10.3',' ','EQUINOX',JEQ,ISTAT)
      IF (ISTAT.NE.0)
     1        CALL TERROR(NEWTBL,15,'Could not create EQUINOX column')
C
C     Initialize:
C
   20 ITEM=0
C
   21 NROW=NROW+1
C
C       Read row:
C
      CALL TBERDC (ITBL, NROW, KOBJ,  OBJECT, NULL, ISTAT)
      IF (ISTAT.NE.0)
     1        CALL TERROR(NEWTBL,21,'Could not read OBJECT column')
C
      IF (OBJECT.EQ.OLDOBJ .OR.OLDOBJ.EQ.' ') THEN
C         continue collecting data for this star.
          OLDOBJ=OBJECT
          ITEM=ITEM+1
          IF (ITEM.GT.MXITEM) THEN
             CALL TV('Too many data for this star.')
             CALL STETER(21,'Increase parameter MXITEM and recompile.')
          END IF
      ELSE
C         Summarize this star.
          GO TO 30
      END IF
C
      CALL TBERDR(ITBL, NROW, KRA,  RAS(ITEM),NULL, ISTAT)
      IF (ISTAT.NE.0)
     1        CALL TERROR(NEWTBL,21,'Could not read RA column')
      CALL TBERDR(ITBL, NROW, KDEC,  DECS(ITEM),NULL, ISTAT)
      IF (ISTAT.NE.0)
     1        CALL TERROR(NEWTBL,21,'Could not read DEC column')
      CALL TBERDR(ITBL, NROW, KEQ,  EQUINOXS(ITEM),NULL, ISTAT)
      IF (ISTAT.NE.0)
     1        CALL TERROR(NEWTBL,21,'Could not read EQUINOX column')
C
      IF (NROW.EQ.NROWS) GO TO 30
      GO TO 21
C
C
C       Summarize this star:
C
   30 NROUT=NROUT+1
      ITEMS=ITEM
C
      EQUINOX=EQUINOXS(1)
      IF (ITEM.EQ.1) THEN
         RA=RAS(1)
         DEC=DECS(1)
         CARD='Only one observation of '//OLDOBJ
         CALL TVN(CARD)
         GO TO 80
      ELSE
C        extract robust estimates.
         CARD='Processing '//OLDOBJ
         CALL TVN(CARD)
      END IF
C
C       (here is where we should display scatter.)
C
      CALL SORT1(RAS,ITEMS)
      CALL SORT1(DECS,ITEMS)
C
      RA=0.5*(RAS((ITEMS+1)/2) + RAS(ITEMS/2+1))
      DEC=0.5*(DECS((ITEMS+1)/2) + DECS(ITEMS/2+1))
C
C
C        Write summary row:
C
   80 CALL TBEWRC (NEWTBL, NROUT, JOBJ, OLDOBJ,  ISTAT)
      IF (ISTAT.NE.0)
     1        CALL TERROR(NEWTBL,81,'Could not write OBJECT column')
      CALL TBEWRR (NEWTBL, NROUT, JRA, RA,  ISTAT)
      IF (ISTAT.NE.0)
     1        CALL TERROR(NEWTBL,82,'Could not write RA column')
      CALL TBEWRR (NEWTBL, NROUT, JDEC, DEC,  ISTAT)
      IF (ISTAT.NE.0)
     1        CALL TERROR(NEWTBL,83,'Could not write DEC column')
      CALL TBEWRR (NEWTBL, NROUT, JEQ, EQUINOX,  ISTAT)
      IF (ISTAT.NE.0)
     1        CALL TERROR(NEWTBL,84,'Could not write EQUINOX column')
C
      IF (NROW.EQ.NROWS) GO TO 90
C
C     Prepare for next star:
      OLDOBJ=OBJECT
      NROW=NROW-1
      GO TO 20
C
C
  90  CALL TBTCLO(ITBL,  ISTAT)
      IF (ISTAT.NE.0) CALL TERROR(ITBL,1,'Could not close "sdata.tbl".')
      CALL TV('   sdata.tbl closed.')
C
      CALL TBTCLO(NEWTBL,  ISTAT)
      CARD='Could not close '//TBLFIL(1:60)                    ! RHW 4/10/93
      IF (ISTAT.NE.0) CALL TERROR (NEWTBL,1,CARD)
      CARD='File '//TBLFIL(1:60)//' closed'                    ! RHW 4/10/93
      CALL TV(CARD)
C
C
      CALL STSEPI
C
      END
