C @(#)optop.for	19.1 (ES0-DMD) 02/25/03 14:27:53
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 OPTOP
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENTIFICATION: program OPTOP 
C.PURPOSE: to provide X and Y coordinates for plates to be used at the
C          Cassegrain focus of the 3.6mt telescope.
C.AUTHOR:  A. Gemmo         Padova Department of Astronomy
C.VERSION: 861101 G. Lund       Creation 
C          900625 S. Cristiani  Last update on HP machine
C          910301 A. Gemmo      New portable version running in MIDAS env.
C          910524 A. Gemmo      Modified for OPTOPUS context.
C          910717 A. Gemmo      Modified to handle "-00" in coordinates
C                               of the center. DEC read as a string.  
C------------------------------------------------------------------------------
       IMPLICIT NONE
C
C ***
C
       INTEGER		MADRID(1)
       INTEGER		NRINP
       INTEGER		NCINP
       INTEGER		NROUT
       INTEGER		NCOUT
       PARAMETER	(NCINP=7)
       PARAMETER	(NCOUT=9)
       PARAMETER        (NROUT=300)
       INTEGER		ISTAT,IAC,IAV
       INTEGER		LFIEL1,LFIEL2,LFIEL3,LFIEL4
       INTEGER          LFIEL5,LFIEL6,LFIEL7
       INTEGER          TYP1,TYP2,TYP3,TYP4
       INTEGER          TYP5,TYP6,TYP7
       INTEGER		KUN,KNUL
       INTEGER		NACINP,NARINP,NSINP
       INTEGER		OUTTYP,OUTCOL(10),COLNUM(10),CN(10)
       INTEGER		TIDINP
       INTEGER		TIDOUT
       INTEGER		NBIG
       INTEGER		I1,I7,I8,I9,I10,NCOL,II
       INTEGER       	I12,I13,I14,I15,I16,I17,I18,J4,J5,J7
       INTEGER          K1,K2,NCINP1,IND1,IND2,IND3
C
C ***  
C
       INTEGER		ICOUN0,ICOUN1,ICOUN2,ICOUN3
       INTEGER		FLAG(NROUT),FLAG1(NROUT),FLAG2(NROUT)
C
C ***
C
       REAL             CDEC(3)
C
C ***
C
       DOUBLE PRECISION		FRAD,SGSEP,BGSEP,OOSEP,OGSEP
       DOUBLE PRECISION		SCALFA,EXDGUI,EXDOBJ
       DOUBLE PRECISION 	XFAREA,YFAREA
       DOUBLE PRECISION         CRA(3),DCDEC(3),CRA1,CDEC1
       DOUBLE PRECISION		DRA(NROUT),DDEC(NROUT)
       DOUBLE PRECISION 	X(NROUT),Y(NROUT),Z(NROUT)
       DOUBLE PRECISION         PI  
       DOUBLE PRECISION		RARAD,DECRAD,CRARAD,CDECRA
       DOUBLE PRECISION 	DISTAN
       DOUBLE PRECISION         FAR1,FAR2  
       DOUBLE PRECISION		DBGUID,DSGUID
       DOUBLE PRECISION         XX(360),YY(360)
       DOUBLE PRECISION         RAXRAD(360),DECYRAD(360)
       DOUBLE PRECISION         RAX(360),DECY(360)
       DOUBLE PRECISION         XF(4),YF(4),XFRAD(4),YFRAD(4)
       DOUBLE PRECISION         FX(4),FY(4)
       DOUBLE PRECISION         EQUI0,EQUI1,PCRARA,PCDECR
       DOUBLE PRECISION         PCRA,PCDEC
C
C ***
       CHARACTER*80	STRING(25)
       CHARACTER*60	INPFIL
       CHARACTER*60	OUTFIL
       CHARACTER*16	ILAB(NCINP),RLAB(NCINP),OUTLAB
       CHARACTER*16	OUTUNI
       CHARACTER*16	OUTFOR
       CHARACTER*8      FORM
       CHARACTER*16     FORMC16,FORMC1,FORMR8
       CHARACTER*16     IDENT(NROUT)
       CHARACTER*1      TYPE(NROUT)
       CHARACTER*1      CHECK(NROUT)
       CHARACTER*80     CDECCH
       CHARACTER*80     DECCHA(3)
       CHARACTER*1      SIGN
       CHARACTER*1      ACFLAG,PFLAG
C
C ***
       LOGICAL          NULL
C
C ***
       INCLUDE          'MID_INCLUDE:TABLES.INC'
       COMMON		/VMR/ MADRID
       INCLUDE          'MID_INCLUDE:TABLED.INC'
C
C ***
       DATA 	ILAB/'IDENT   ','TYPE    ','CHECK   ',
     2               'OLDRA   ','OLDDEC  ','RA      ','DEC     '/
C
C ***
C
       DATA    FORMC16/'A16'/
       DATA    FORMC1 /'A1'/
       DATA    FORMR8 /'F7.0'/
       DATA    PI/3.14159265358979D0/
       DATA    ICOUN0,ICOUN1,ICOUN2,ICOUN3/0,0,0,0/
       DATA    DBGUID,DSGUID/0.0D0,0.0D0/
       DATA    FX/-20010.,-20010.,20010.,20010./
       DATA    FY/137000.,118100.,118100.,137000./
       DATA    SIGN/'+'/
C
C *** start definition of constraints of starplates
C
       DATA	FRAD   /1.37D5/  !radius of the field                (microns)
       DATA	SGSEP  /3.45D3/  !min. sep. between small guidestars (microns)
       DATA	BGSEP  /1.25D4/  !min. sep. between big guidestars   (microns)
       DATA	NBIG   /2/       !number of big guidestars
       DATA	OOSEP  /3.45D3/  !min. object-object separation      (microns)
       DATA	OGSEP  /9.00D3/  !min. object-guidestar separation   (microns)
       DATA 	SCALFA /1.56D-7/ !scal. fact. for depth of the holes (microns)
       DATA	EXDGUI /-6.00D2/ !extra-depth for guidestars         (microns)
       DATA     EXDOBJ /-1.49D3/ !extra-depth for objects            (microns)
       DATA	XFAREA /2.001D4/ !half X-dimension of forbidden area (microns)
       DATA	YFAREA /1.890D4/ !Y-dimension of forbidden area      (microns)
C
C *** start the code
C
         CALL STSPRO('OPTOP')
C
       DO II=1,NROUT
          FLAG (II) = 0
          FLAG1(II) = 0
          FLAG2(II) = 0
       ENDDO
C
C *** get the input table and the output table
C
       CALL STKRDC('INPUTFI',1,1,60,IAC,INPFIL,KUN,KNUL,ISTAT)
       CALL STKRDC('OUTPUTF',1,1,60,IAC,OUTFIL,KUN,KNUL,ISTAT)
C
C *** open input table
C
       CALL TBTOPN(INPFIL,F_I_MODE,TIDINP,ISTAT)
C
C *** get information about the input table
C
       CALL TBIGET(TIDINP,NCOL,NRINP,NSINP,NACINP,NARINP,ISTAT)
       IF(NRINP.EQ.0)THEN      !no rows in table
         STRING(1) = '*** FATAL: There are no data in the input TABLE'
         CALL STETER(9,STRING(1))
       ENDIF
C
       CALL TBLSER(TIDINP,'IDENT',COLNUM(1),ISTAT)
       CALL TBLSER(TIDINP,'TYPE',COLNUM(2),ISTAT)
       CALL TBLSER(TIDINP,'CHECK',COLNUM(3),ISTAT)
       CALL TBLSER(TIDINP,'OLDRA',COLNUM(4),ISTAT)
       CALL TBLSER(TIDINP,'OLDDEC',COLNUM(5),ISTAT)
       CALL TBLSER(TIDINP,'RA',COLNUM(6),ISTAT)
       CALL TBLSER(TIDINP,'DEC',COLNUM(7),ISTAT)
C
       DO I1=1,NCINP
         IF(COLNUM(I1).GE.1)THEN
           CALL TBLGET(TIDINP,COLNUM(I1),RLAB(I1),ISTAT)
         ELSE
           GOTO 1902
         ENDIF
         IF(RLAB(I1).NE.ILAB(I1))THEN   
            STRING(2) = '*** FATAL: Wrong column label'
           CALL STETER(9,STRING(2))
         ENDIF
1902   ENDDO
C
       IF (COLNUM(1).GE.1)THEN
          CALL TBFGET(TIDINP,COLNUM(1),FORM,LFIEL1,TYP1,ISTAT)
          IF (TYP1.NE.D_C_FORMAT)THEN
             STRING(3) = '*** FATAL: Wrong `IDENT` column format'
             CALL STETER(9,STRING(3))
          ENDIF   
       ENDIF

       IF (COLNUM(2).GE.1)THEN
          CALL TBFGET(TIDINP,COLNUM(2),FORM,LFIEL2,TYP2,ISTAT)
          IF (TYP2.NE.D_C_FORMAT)THEN
             STRING(4) = '*** FATAL: Wrong `TYPE` column format'
             CALL STETER(9,STRING(4))
          ENDIF
       ENDIF

       IF (COLNUM(3).GE.1)THEN
          CALL TBFGET(TIDINP,COLNUM(3),FORM,LFIEL3,TYP3,ISTAT)
          IF (TYP3.NE.D_C_FORMAT)THEN
             STRING(5) = '*** FATAL: Wrong `CHECK` column format'
             CALL STETER(9,STRING(5))
          ENDIF
       ENDIF

       IF (COLNUM(4).GE.1)THEN
          CALL TBFGET(TIDINP,COLNUM(4),FORM,LFIEL4,TYP4,ISTAT)
          IF (TYP4.NE.D_R8_FORMAT)THEN
             STRING(6) = '*** FATAL: Wrong `OLDRA` column format'
             CALL STETER(9,STRING(6))
          ENDIF
       ENDIF

       IF (COLNUM(5).GE.1)THEN
          CALL TBFGET(TIDINP,COLNUM(5),FORM,LFIEL5,TYP5,ISTAT)
          IF (TYP5.NE.D_R8_FORMAT)THEN
            STRING(7) = '*** FATAL: Wrong `OLDDEC` column format'
            CALL STETER(9,STRING(7))	  
          ENDIF
       ENDIF

       IF (COLNUM(6).GE.1)THEN
          CALL TBFGET(TIDINP,COLNUM(6),FORM,LFIEL6,TYP6,ISTAT)
          IF (TYP6.NE.D_R8_FORMAT)THEN
             STRING(8) = '*** FATAL: `RA` column format'
             CALL STETER(9,STRING(8))	  
          ENDIF
       ENDIF

       IF (COLNUM(7).GE.1)THEN
          CALL TBFGET(TIDINP,COLNUM(7),FORM,LFIEL7,TYP7,ISTAT)
          IF (TYP7.NE.D_R8_FORMAT)THEN
             STRING(9) = '*** FATAL: Wrong `DEC` column format'
             CALL STETER(9,STRING(9))	  
          ENDIF
       ENDIF
C
C *** create the output table
C
       CALL TBTINI(OUTFIL,0,F_O_MODE,NCOUT,NROUT,TIDOUT,ISTAT)
C
       OUTTYP = TYP1     !initialize ident column
       OUTFOR = FORMC16
       OUTUNI = '        '
       OUTLAB = 'IDENT   '
       CALL TBCINI(TIDOUT,OUTTYP,16,OUTFOR,OUTUNI,
     2                     OUTLAB,OUTCOL(1),ISTAT)      
C
       OUTTYP = D_C_FORMAT    !initialize type column
       OUTFOR = FORMC1
       OUTUNI = '        '
       OUTLAB = 'TYPE    '
       CALL TBCINI(TIDOUT,OUTTYP,1,OUTFOR,OUTUNI,
     2                     OUTLAB,OUTCOL(2),ISTAT)
C
       OUTTYP = D_C_FORMAT    !initialize check column
       OUTFOR = FORMC1
       OUTUNI = '        '
       OUTLAB = 'CHECK   '
       CALL TBCINI(TIDOUT,OUTTYP,1,OUTFOR,OUTUNI,
     2                     OUTLAB,OUTCOL(3),ISTAT)
C
       OUTTYP = D_R8_FORMAT    !initialize x column
       OUTFOR = FORMR8
       OUTUNI = 'MICRONS'
       OUTLAB = 'X      '
       CALL TBCINI(TIDOUT,OUTTYP,1,OUTFOR,OUTUNI,
     2                     OUTLAB,OUTCOL(4),ISTAT)
C 
       OUTTYP = D_R8_FORMAT    !initialize y column
       OUTFOR = FORMR8
       OUTUNI = 'MICRONS'
       OUTLAB = 'Y      '
       CALL TBCINI(TIDOUT,OUTTYP,1,OUTFOR,OUTUNI,
     2                     OUTLAB,OUTCOL(5),ISTAT)
C
       OUTTYP = D_R8_FORMAT    !initialize z column
       OUTFOR = FORMR8
       OUTUNI = 'MICRONS'
       OUTLAB = 'Z      '
       CALL TBCINI(TIDOUT,OUTTYP,1,OUTFOR,OUTUNI,
     2                     OUTLAB,OUTCOL(6),ISTAT)
C
C *** read center coordinates (if present)
C
       CALL STKRDD('PLATECEN',1,1,IAV,CRA(1),KUN,KNUL,ISTAT)
       CALL STKRDD('PLATECEN',2,1,IAV,CRA(2),KUN,KNUL,ISTAT)
       CALL STKRDD('PLATECEN',3,1,IAV,CRA(3),KUN,KNUL,ISTAT)
       CALL STKRDC('ACFLAG',1,1,1,IAC,ACFLAG,KUN,KNUL,ISTAT)
       IF(ACFLAG.EQ.'Y'.OR.ACFLAG.EQ.'y')THEN
         CALL STKRDD('PLATECEN',4,1,IAV,DCDEC(1),KUN,KNUL,ISTAT)
         CALL STKRDD('PLATECEN',5,1,IAV,DCDEC(2),KUN,KNUL,ISTAT)
         CALL STKRDD('PLATECEN',6,1,IAV,DCDEC(3),KUN,KNUL,ISTAT)
         GOTO 8888
       ELSE IF(ACFLAG.EQ.'N'.OR.ACFLAG.EQ.'n')THEN
         CALL STKRDC('PLATECHA',1,1,80,IAV,CDECCH,KUN,KNUL,ISTAT)
       ENDIF
C
       IND1 = INDEX(CDECCH,',')
       IND2 = INDEX(CDECCH(IND1+1:80),',') + IND1
       DECCHA(1) = CDECCH(1:IND1-1)
       IND3 = INDEX(CDECCH,'-')
       DECCHA(2) = CDECCH(IND1+1:IND2-1)
       DECCHA(3) = CDECCH(IND2+1:80)
C
C *** decode character DECCHA to real array
C
      CALL USRINP(CDEC(1),80,'R',DECCHA(1))
      CALL USRINP(CDEC(2),80,'R',DECCHA(2))
      CALL USRINP(CDEC(3),80,'R',DECCHA(3))
C
C *** write decoded dec in keyword
C
      DCDEC(1) = DBLE(CDEC(1))
      DCDEC(2) = DBLE(CDEC(2))
      DCDEC(3) = DBLE(CDEC(3))
C
      IF(DCDEC(1).LT.-90..OR.DCDEC(1).GT.90.)THEN
        CALL STETER(ISTAT,'*** FATAL: Wrong DEC-degrees !!')
      ELSE IF(DCDEC(2).GT.60.)THEN
        CALL STETER(ISTAT,'*** FATAL: Wrong DEC-arcmin !!')
      ELSE IF(DCDEC(3).GT.60.)THEN
        CALL STETER(ISTAT,'*** FATAL: Wrong DEC-arcsec !!')
      ENDIF
C
      CALL STKWRD('PLATECEN',DCDEC(1),4,1,KUN,ISTAT)
      CALL STKWRD('PLATECEN',DCDEC(2),5,1,KUN,ISTAT)
      CALL STKWRD('PLATECEN',DCDEC(3),6,1,KUN,ISTAT)
C
C *** convert RA and DEC to decimal format
C
8888  CRA1 = CRA(1)+CRA(2)/60.+CRA(3)/3600.
      IF(DCDEC(1).LT.0..OR.
     1    DCDEC(2).LT.0..OR.
     2    DCDEC(3).LT.0..OR.
     3    IND3.NE.0) THEN
       SIGN(1:1) = '-'
       DCDEC(1) = DABS(DCDEC(1))
       DCDEC(2) = DABS(DCDEC(2))
       DCDEC(3) = DABS(DCDEC(3))
      ENDIF 
       CDEC1 = DCDEC(1)+DCDEC(2)/60.+DCDEC(3)/3600.
       IF(SIGN(1:1).EQ.'-')CDEC1 = -CDEC1
C  
C *** do the work
C
C
C *** start computation of X and Y positions
C
       DO I7=1,NRINP      !read in table elements
          CALL TBERDC(TIDINP,I7,COLNUM(1),IDENT(I7),KNUL,ISTAT)
          CALL TBERDC(TIDINP,I7,COLNUM(2),TYPE(I7),KNUL,ISTAT)
C
C *** convert TYPE to uppercase
C
          IF(TYPE(I7).EQ.'o')TYPE(I7) = 'O'
          IF(TYPE(I7).EQ.'s')TYPE(I7) = 'S'
          IF(TYPE(I7).EQ.'b')TYPE(I7) = 'B'
C
          CALL TBERDD(TIDINP,I7,COLNUM(6),DRA(I7),KNUL,ISTAT)
          CALL TBERDD(TIDINP,I7,COLNUM(7),DDEC(I7),KNUL,ISTAT)
       ENDDO
C
C *** convert center RA and DEC to radians
C
       CRARAD     = CRA1*1.5D1*(PI/1.8D2)
       CDECRA     = CDEC1*(PI/1.8D2)
C
C *** precess center coordinates
C
       CALL STKRDC('PFLAG',1,1,1,IAC,PFLAG,KUN,KNUL,ISTAT)
       CALL STKRDD('EQUINOX',1,1,IAC,EQUI0,KUN,KNUL,ISTAT)
       CALL STKRDD('EQUINOX',2,1,IAC,EQUI1,KUN,KNUL,ISTAT)
C
       IF(PFLAG.EQ.'N'.OR.PFLAG.EQ.'n')THEN
         CALL STTPUT('     ',ISTAT)
      CALL STTPUT('No precession applied to center coordinates',
     1ISTAT)
         CALL STTPUT('     ',ISTAT)
         PCRARA = CRARAD
         PCDECR = CDECRA
       ELSE IF(PFLAG.EQ.'Y'.OR.PFLAG.EQ.'y')THEN
         CALL PRE(CRARAD,CDECRA,PCRARA,PCDECR,EQUI0,EQUI1)
       ENDIF
C
C *** convert precessed RA and DEC to hours and degrees
C
       PCRA = PCRARA*(1.8D2/PI)/1.5D1
       PCDEC = PCDECR*(1.8D2/PI)
       CALL STKWRD('PLATEC1',PCRA,1,1,KUN,ISTAT)
       CALL STKWRD('PLATEC1',PCDEC,2,1,KUN,ISTAT)
C
C *** subroutine ADXY taken from Yanling Fang's astrometric package
C     (modified by A. Gemmo)
       DO J7=1,NRINP
          RARAD  = DRA(J7)*1.5D1*(PI/1.8D2)
          DECRAD  = DDEC(J7)*(PI/1.8D2)
          CALL ADXY(RARAD,DECRAD,PCRARA,PCDECR,X(J7),Y(J7))
       ENDDO
C
C *** fill in table elements
C
       DO J4=1,NRINP
           CALL TBEWRC(TIDOUT,J4,OUTCOL(1),IDENT(J4),ISTAT)
           CALL TBEWRC(TIDOUT,J4,OUTCOL(2),TYPE(J4),ISTAT)
           CALL TBEWRD(TIDOUT,J4,OUTCOL(4),X(J4),ISTAT)
           CALL TBEWRD(TIDOUT,J4,OUTCOL(5),Y(J4),ISTAT)
       ENDDO
C
C *** check maximum distance from the center
C       
       DO I8=1,NRINP
           DISTAN = DSQRT(X(I8)**2+Y(I8)**2)
           IF(DISTAN.GE.FRAD)THEN
            ICOUN1 = ICOUN1+1
            WRITE(STRING(7),222)IDENT(I8)
            WRITE(STRING(8),333)TYPE(I8),X(I8),Y(I8)
222    FORMAT(1X,'Object: ',A16)
333    FORMAT(1X,'Type: ',A1,',',' X = ',F10.0,' Y = ',F10.0,
     2                ' is too far from the center')
            CALL STTPUT(STRING(7),ISTAT)
            CALL STTPUT(STRING(8),ISTAT)
            CHECK(I8) = 'N'
          ENDIF
       ENDDO
       IF(ICOUN1.NE.0)THEN
         WRITE(STRING(9),444)ICOUN1,NRINP
444    FORMAT(1X,'You are going to loose',I3,' object(s) out of ',I3)
         CALL STTPUT(STRING(9),ISTAT)
         STRING(10)=' You can delete it/them or change the center positi
     2on'
          CALL STTPUT(STRING(10),ISTAT)
       ENDIF
C
C *** check if there is any object in the forbidden area
C
       FAR1 = XFAREA
       FAR2 = FRAD-YFAREA
       DO I9=1,NRINP
          IF(DABS(X(I9)).LE.FAR1.AND.Y(I9).GE.FAR2)THEN
            ICOUN2 = ICOUN2+1
            WRITE(STRING(11),555)IDENT(I9)
            WRITE(STRING(12),666)TYPE(I9),X(I9),Y(I9)
555    FORMAT(1X,'Object: ',A16)
666    FORMAT(1X,'Type: ',A1,',',' X = ',F10.0,' Y = ',F10.0,
     2                ' is in the forbidden area')
            CALL STTPUT(STRING(11),ISTAT)
            CALL STTPUT(STRING(12),ISTAT)
            CHECK(I9) = 'N'
          ENDIF
       ENDDO
       IF(ICOUN2.NE.0)THEN
          WRITE(STRING(13),777)ICOUN2
777    FORMAT(1X,I3,' object(s) is/are in the forbidden area')
          CALL STTPUT(STRING(13),ISTAT)
          WRITE(STRING(13),778)
778       FORMAT(1X,' You MUST delete it/them or change the center posit
     2ion')
          CALL STTPUT(STRING(13),ISTAT)
       ENDIF
C
C *** check the number of big guidestars 
C
       DO I10=1,NRINP
          IF(TYPE(I10).EQ.'B')THEN
            ICOUN3 = ICOUN3+1
          ENDIF
       ENDDO
       IF(ICOUN3.LT.NBIG)THEN
         STRING(14)= ' *** FATAL: Your `BIG` guidestars are not enough'
         CALL STETER(9,STRING(14))
       ENDIF
C
C *** check the reciprocal superpositions of big guidestars
C
       DO I12=1,NRINP
           IF(TYPE(I12).EQ.'B'.AND.FLAG(I12).EQ.0)THEN
               DO I13=1,NRINP
                IF(TYPE(I13).EQ.'B'.AND.I12.NE.I13.AND.
     2              FLAG(I13).EQ.0)THEN
                    DISTAN = DSQRT((X(I12)-X(I13))**2
     2                            +(Y(I12)-Y(I13))**2)
                    DBGUID = DMAX1(DBGUID,DISTAN)
                   IF(DISTAN.LT.BGSEP)THEN
       	              WRITE(STRING(16),888)IDENT(I12),IDENT(I13)
888    FORMAT(1X,' Big guidestars: ',A16,' and ',A16,' are too close')
                      CALL STTPUT(STRING(16),ISTAT)
                      CHECK(I12) = 'N'
                      CHECK(I13) = 'N'
                      FLAG(I12) = FLAG(I12)+1
                      FLAG(I13) = FLAG(I13)+1
                    ENDIF
       	         ELSE 
                    GOTO 1111
                 ENDIF
1111           ENDDO
            ELSE 
               GOTO 1112
            ENDIF
1112      ENDDO
C
C *** check the reciprocal superpositions of small guidestars and objects
C
       DO I14=1,NRINP
           IF(TYPE(I14).NE.'B'.AND.FLAG1(I14).EQ.0)THEN
               DO I15=1,NRINP
                  IF(TYPE(I15).NE.'B'.AND.I14.NE.I15.AND.
     2              FLAG1(I15).EQ.0)THEN
                    DISTAN = DSQRT((X(I14)-X(I15))**2
     2                            +(Y(I14)-Y(I15))**2)
                    DSGUID = DMAX1(DSGUID,DISTAN)
                    IF (DISTAN.LT.SGSEP)THEN
                      IF(TYPE(I14).EQ.'S'.AND.TYPE(I15).EQ.'S')THEN
                        WRITE(STRING(18),2222)IDENT(I14),IDENT(I15)
2222   FORMAT(1X,' Small guidestars: ',A16,' and ',A16,
     2           ' are too close')
                        CALL STTPUT(STRING(18),ISTAT)
                        CHECK(I14) = 'N'
                        CHECK(I15) = 'N'
                        FLAG1(I14) = FLAG1(I14)+1
                        FLAG1(I15) = FLAG1(I15)+1

       	              ELSE IF(TYPE(I14).EQ.'S'.AND.TYPE(I15).EQ.'O')THEN
                        WRITE(STRING(18),2223)IDENT(I14),IDENT(I15)
2223   FORMAT(1X,' Small guidestar: ',A16,' and object:', A16)
                        CALL STTPUT(STRING(18),ISTAT)
                        WRITE(STRING(18),3333)
3333   FORMAT(1X,' are too close')
                        CALL STTPUT(STRING(18),ISTAT)
                        CHECK(I14) = 'N'
                        CHECK(I15) = 'N'
                        FLAG1(I14) = FLAG1(I14)+1
                        FLAG1(I15) = FLAG1(I15)+1

                      ELSE IF(TYPE(I14).EQ.'O'.AND.TYPE(I15).EQ.'S')THEN
                        WRITE(STRING(18),2224)IDENT(I14),IDENT(I15)
2224   FORMAT(1X,' Object: ',A16,' and small guidestar:', A16)
                        CALL STTPUT(STRING(18),ISTAT)
                        WRITE(STRING(18),3334)
3334                    FORMAT(1X,' are too close')
                        CALL STTPUT(STRING(18),ISTAT)
                        CHECK(I14) = 'N'
                        CHECK(I15) = 'N'
       		        FLAG1(I14) = FLAG1(I14)+1
                        FLAG1(I15) = FLAG1(I15)+1

                      ELSE IF(TYPE(I14).EQ.'O'.AND.TYPE(I15).EQ.'O')THEN
                        WRITE(STRING(18),2225)IDENT(I14),IDENT(I15)
2225   FORMAT(1X,' Objects: ',A16,' and ',A16,' are too close')
                        CALL STTPUT(STRING(18),ISTAT)
                        CHECK(I14) = 'N'
                        CHECK(I15) = 'N'
                        FLAG1(I14) = FLAG1(I14)+1
                        FLAG1(I15) = FLAG1(I15)+1
                      ENDIF
                    ENDIF
       	         ELSE 
                    GOTO 1113
                 ENDIF
1113           ENDDO
            ELSE 
               GOTO 1114
            ENDIF
1114      ENDDO
C
C *** check the reciprocal superpositions of small and big guidestars
C
       DO I16=1,NRINP
           IF(TYPE(I16).EQ.'B'.AND.FLAG2(I16).EQ.0)THEN
               DO I17=1,NRINP
                  IF(TYPE(I17).NE.'B'.AND.I16.NE.I17.AND.
     2              FLAG2(I17).EQ.0)THEN
                    DISTAN = DSQRT((X(I16)-X(I17))**2
     2                            +(Y(I16)-Y(I17))**2)
                    DSGUID = DMAX1(DSGUID,DISTAN)
                   IF(DISTAN.LT.OGSEP)THEN
       	             IF (TYPE(I16).EQ.'B'.AND.TYPE(I17).EQ.'S')THEN
                        WRITE(STRING(20),4444)IDENT(I16),IDENT(I17)
4444   FORMAT(1X,' Big guidestar: ',A16, ' and small guidestar: ',A16)
                        CALL STTPUT(STRING(20),ISTAT)
                        WRITE(STRING(21),5555)
5555   FORMAT(1X,' are too close')
                        CALL STTPUT(STRING(21),ISTAT)
                        CHECK(I16) = 'N'
                        CHECK(I17) = 'N'

       	             ELSE IF(TYPE(I16).EQ.'B'.AND.TYPE(I17).EQ.'O')THEN
                        WRITE(STRING(20),4441)IDENT(I16),IDENT(I17)
4441   FORMAT(1X,' Big guidestar: ',A16,' and object: ',A16)
                        WRITE(STRING(21),5551)
5551   FORMAT(1X,' are too close')
                        CALL STTPUT(STRING(21),ISTAT)
                        CHECK(I16) = 'N'
                        CHECK(I17) = 'N'
                        FLAG2(I16) = FLAG2(I16)+1
                        FLAG2(I17) = FLAG2(I17)+1
                      ENDIF
                    ENDIF
       	    ELSE 
                      GOTO 4442
                   ENDIF
4442            ENDDO
            ELSE 
               GOTO 4443
             ENDIF
4443      ENDDO
C
C *** compute depth of the hole for big guidestars, small guidestars 
C     and objects
C
       DO I18=1,NRINP
          IF(TYPE(I18).EQ.'B')THEN
            Z(I18) = -(X(I18)*X(I18)+Y(I18)*Y(I18))*SCALFA+EXDGUI
           ELSE IF(TYPE(I18).NE.'B')THEN
             Z(I18) = -(X(I18)*X(I18)+Y(I18)*Y(I18))*SCALFA+EXDOBJ
           ENDIF
          X(I18)=-X(I18)
       ENDDO
C
C *** fill in the output table
C
       DO J5=1,NRINP
          CALL TBEWRC(TIDOUT,J5,OUTCOL(3),CHECK(J5),ISTAT)
          CALL TBEWRD(TIDOUT,J5,OUTCOL(6),Z(J5),ISTAT)
       ENDDO
C
C *** calculate circle to show plate limits
C
       CALL TBTOPN('circle',F_IO_MODE,TIDINP,ISTAT)   !open table
C
       CALL TBIGET(TIDINP,NCINP1,NRINP,NSINP,NACINP,NARINP,ISTAT) !table info
C
       DO K2=1,NRINP    !fill in RAX and DECY columns in circle.tbl
C
C *** subroutine XYAD taken from Yanling Fang's astrometric package
C              (God help us!)
         CALL TBLSER(TIDINP,'XX',CN(1),ISTAT)
         CALL TBLSER(TIDINP,'YY',CN(2),ISTAT)
         CALL TBERDD(TIDINP,K2,CN(1),XX(K2),NULL,ISTAT)
         CALL TBERDD(TIDINP,K2,CN(2),YY(K2),NULL,ISTAT)
         CALL XYAD(XX(K2),YY(K2),PCRARA,PCDECR,
     2             RAXRAD(K2),DECYRAD(K2))
         RAX(K2) = RAXRAD(K2)*360./(2.*15.*PI)
         DECY(K2) = DECYRAD(K2)*360./(2.*PI)
         CALL TBLSER(TIDINP,'RAX',CN(3),ISTAT)
         CALL TBLSER(TIDINP,'DECY',CN(4),ISTAT) 
         CALL TBEWRD(TIDINP,K2,CN(3),RAX(K2),ISTAT)
         CALL TBEWRD(TIDINP,K2,CN(4),DECY(K2),ISTAT)
       ENDDO
C
       CALL STKWRD('CIRDIM',RAX(1),1,1,KUN,ISTAT)
       CALL STKWRD('CIRDIM',DECY(90),2,1,KUN,ISTAT)
C
C *** conversion of forbidden area coordinates
C
       DO K1=1,4
         CALL XYAD(FX(K1),FY(K1),PCRARA,PCDECR,XFRAD(K1),YFRAD(K1))
         XF(K1) = XFRAD(K1)*360./(2.*15.*PI)
         CALL STKWRD('XF',XF(K1),K1,1,KUN,ISTAT)
         YF(K1) = YFRAD(K1)*360./(2.*PI)
         CALL STKWRD('YF',YF(K1),K1,1,KUN,ISTAT)
       ENDDO
C
C *** initialize row selection flag
C
       CALL TBSINI(TIDOUT,ISTAT)
C
C *** close the table
C
       CALL TBTCLO(TIDOUT,ISTAT)
C
C *** over and out
C
       CALL STSEPI
       END
C 
C *** Subroutines used in the code:
C
      SUBROUTINE PRE(A0,D0,A1,D1,EQX0,EQX1)
C
C ... calculate general precession for two given epochs
C
C ... References: STUMPFF P.,ASTRON. ASTROPHYS. SUPPL. 41, P. 1 (1980)
C ...             EXPLAN. SUPPL. ASTRON. EPHEREMERIS, P. 28 (1961)
C
      IMPLICIT NONE             !REAL*8 (A-H,O-Z)
C
      DOUBLE PRECISION A0,D0,D1,EQX0,EQX1
      DOUBLE PRECISION DT0,DT,DTS,DTC,DC1900,DC1M2
      DOUBLE PRECISION DC1,DC2,DC3,DC4,DC5,DC6,DC7,DC8,DC9
      DOUBLE PRECISION DZED,DZETA,DCSAR,DTHETA
      DOUBLE PRECISION PI,SINR,COSR,A1,R
C
      DATA DCSAR/4.848136812D-6/,DC1900/1900.0D0/,DC1M2/0.01D0/,
     * DC1/2304.25D0/,DC2/1.396D0/,DC3/0.302D0/,DC4/0.018D0/,
     * DC5/0.791D0/,DC6/2004.683D0/,DC7/-0.853D0/,DC8/-0.426D0/,
     * DC9/-0.042D0/
C
      PI = 3.1415926535897928
C
C ... first determine the precession angles zeta, z and theta
C
      DT0=(EQX0-DC1900)*DC1M2
      DT=(EQX1-EQX0)*DC1M2
      DTS=DT*DT
      DTC=DTS*DT
      DZETA=((DC1+DC2*DT0)*DT+DC3*DTS+DC4*DTC)*DCSAR
      DZED=DZETA + DC5*DTS*DCSAR
      DTHETA=((DC6+DC7*DT0)*DT+DC8*DTS+DC9*DTC)*DCSAR
C
C ... with these and A0,D0 compute the precessed coordinates A1,D1
C
      D1 = ASIN(DCOS(DTHETA)*DSIN(D0)+DSIN(DTHETA)*DCOS(D0)*
     1                                              DCOS(A0+DZETA))
      SINR = (DCOS(D0)*DSIN(A0+DZETA))/DCOS(D1)
      COSR = (DCOS(DTHETA)*DCOS(D0)*DCOS(A0+DZETA)-
     1                         DSIN(DTHETA)*DSIN(D0))/DCOS(D1)
      R = DASIN(SINR)
      IF (COSR.LT.0) R=PI-R                ! determine quadrant of R
      A1 = R+DZED                      ! this is the precessed R.A.
      IF (A1.GT.2.0D0*PI) A1=A1-2.0D0*PI
      IF (A1.LT.0.0D0) A1=A1+2.0D0*PI
C
      RETURN
      END
C
C *************** XYAD *******
C
C  The subroutine XYAD transforms X and Y coordinates in RAs and DECs
C
      SUBROUTINE XYAD(XX,YY,RA0RAD,D0RAD,RARAD,DECRAD)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     XX        (INPUT)  Value of X coordinate to be transformed (microns)
C     YY        (INPUT)  Value of Y coordinate to be transformed (microns)
C     RA0RAD    (INPUT)  Center RA                               (radians)
C     D0RAD     (INPUT)  Center DEC                              (radians)
C     RARAD     (OUTPUT) RA                                      (radians)
C     DECRAD    (OUTPUT) DEC                                     (radians)
C------------------------------------------------------------------------------
      IMPLICIT NONE
C
C ***
C
      DOUBLE PRECISION       PI,FOC
      DOUBLE PRECISION       XX,YY,AKSI,ANU
      DOUBLE PRECISION       COSD0,SIND0
      DOUBLE PRECISION       CX,CY,BX(9),BY(9)
      DOUBLE PRECISION       TERM1,TERM2,XX1,YY1
      DOUBLE PRECISION       RH,SINPH,COSPH,SINRH,COSRH
      DOUBLE PRECISION       RARAD,DECRAD
      DOUBLE PRECISION       RA0RAD,D0RAD
C
C ***
C
      DATA         PI/3.14159265358979D0/
      DATA         BX/0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,
     2                0.0D0,0.0D0,0.0D0,0.0D0/
      DATA         BY/0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,
     2                0.0D0,0.0D0,0.0D0,0.0D0/
      DATA         CX,CY/0.0D0,0.0D0/
C
C *** do the work
C
      FOC   = 1.0D0/3.6D3
      BX(1) = 9.9973D-1/0.14006
      BY(2) = BX(1)
C
      COSD0 = DCOS(D0RAD)
      SIND0 = DSIN(D0RAD)
C
      TERM1 = BX(1)*BY(2)-BX(2)*BY(1)
      TERM2 = BX(2)*BY(1)-BX(1)*BY(2)
C
      XX1 = XX*FOC/1.0D3
      YY1 = YY*FOC/1.0D3
C
      ANU = XX1*BY(1)-(YY1*BY(2)*TERM2/TERM1)
     2      -(CY*BX(2)*BY(1)/TERM1)+(CY*BX(1)*BY(2)/TERM1)
      AKSI = (XX1*TERM1+ANU*BX(2)-CY*BX(2))/BY(2)+CX
C
      RH = DSQRT(AKSI**2 + ANU**2)
C
      SINPH = AKSI/RH 
      COSPH = ANU/RH
C
      RH = RH*PI/1.8D2
C
      SINRH= DSIN(RH) 
      COSRH = DCOS(RH)
C
      DECRAD = COSRH*SIND0+SINRH*COSD0*COSPH
      DECRAD = DATAN(DECRAD/DSQRT(1.0D0-DECRAD**2))
C
      RARAD = SINRH*SINPH/DCOS(DECRAD) 
      RARAD = RA0RAD + DATAN(RARAD/DSQRT(1.0D0-RARAD**2))
C      
      RETURN
      END
      
C
C
C *************** ADXY *******
C
C   The subroutine ADXY transforms RAs and DECs in X and Y coordinates 
C 
      SUBROUTINE ADXY(ALFAS,DELTS,ALFA0,DELT0,XX,YY)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     ALFAS (INPUT) Value of RA to be transformed  (radians)
C     DELTS (INPUT) Value of DEC to be transformed (radians)
C     ALFA0 (INPUT) RA of plate center             (radians)
C     DELT0 (INPUT) DEC of plate center            (radians)
C     XX    (OUTPUT) X coordinate                  (microns)
C     YY    (OUTPUT) Y coordinate                  (microns)
C-------------------------------------------------------------------------- 
      IMPLICIT NONE	
C
C ***
C
      INTEGER	IDX(9),IDY(9)
      INTEGER   NX,NY,IXY,NITER
      INTEGER   I16,I17,N,L
C
C ***
C 
      DOUBLE PRECISION  	XX,YY
      DOUBLE PRECISION  	DELTS,ALFAS
      DOUBLE PRECISION          RH,COSRH,SINRH,COSPH,SINPH
      DOUBLE PRECISION          V,AKSI1,ANU1,AKSI,ANU
      DOUBLE PRECISION          DKSI,DNU,DDIST,DCRIT
      DOUBLE PRECISION          DKDX,DKDY,DNDX,DNDY 
      DOUBLE PRECISION          XBIS(10),AR(50,4),CX,CY,BX(9),BY(9) 
      DOUBLE PRECISION          ALFA0,DELT0,COSD0,SIND0
      DOUBLE PRECISION          XEST(2,2),YEST(2,2)
      DOUBLE PRECISION          PI,FOC
C 
C *** peculiar Cassegrain focus parameters
C 
      DATA     NX/3/
      DATA     NY/3/
      DATA     CX/0.0D0/
      DATA     CY/0.0D0/
      DATA     IDX/1,1,1,0,0,0,0,0,0/
      DATA     IDY/1,1,1,0,0,0,0,0,0/
      DATA     BX/0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,
     2           0.0D0,0.0D0,0.0D0,0.0D0/
      DATA     BY/0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,
     2           0.0D0,0.0D0,0.0D0,0.0D0/
C
C ***
C
      DATA     PI/3.14159265358979D0/
C
C *** do the work
C 
      NITER = 0 
      FOC   = 1.0D0/3.6D3
      BX(1) = 9.9973D-1/0.14006
      BY(2) = BX(1)
C
      COSD0 = DCOS(DELT0) 
      SIND0 = DSIN(DELT0) 
C
      IF (ALFA0.LT.PI/2D0.AND.ALFAS.GT.1.5D0*PI) ALFAS = ALFAS-2D0*PI 
      IF (ALFA0.GT.1.5D0*PI.AND.ALFAS.LT.PI/2D0) ALFAS = ALFAS+2D0*PI 
C
      COSRH = COSD0*DCOS(DELTS)*DCOS(ALFAS-ALFA0) + SIND0*DSIN(DELTS) 
      SINRH = DSQRT(1.0D0-COSRH**2) 
      SINPH = DCOS(DELTS)*DSIN(ALFAS-ALFA0)/SINRH 
      COSPH = (DSIN(DELTS)-SIND0*COSRH)/(COSD0*SINRH) 
      RH = DATAN(SINRH/COSRH) 
C 
C ***  estimate X and Y position 
C 
      AKSI = RH*SINPH/(PI/1.8D2)
      ANU = RH*COSPH/(PI/1.8D2) 
C
      XEST(1,1) = ((AKSI-CX)*BY(2)-(ANU-CY)*BX(2))/ 
     2            (BX(1)*BY(2)-BX(2)*BY(1)) 
      YEST(1,1) = ((AKSI-CX)*BY(1)-(ANU-CY)*BX(1))/ 
     2            (BX(2)*BY(1)-BX(1)*BY(2)) 
30    XEST(2,1) = XEST(1,1)*1.0D3/FOC 
      YEST(2,1) = YEST(1,1)*1.0D3/FOC 
C
C60    NITER = NITER + 1 
      NITER = NITER + 1 
C
C70   AR(1,1) = XEST(1,1) 
      AR(1,1) = XEST(1,1) 
      AR(1,2) = YEST(1,1) 
      AR(1,3) = 0.0D0 
      AR(1,4) = 0.0D0 
C
      I16 = 1 
      N = NX+1
      IXY = 1 
      CALL RAR(XBIS,N,IXY,IDX,IDY,AR,I16,0)
      L = NX
      V = 0.0D0 
      DO I17 = 1,L 
         V = V + BX(I17)*XBIS(I17)
      ENDDO
      AKSI1 = V + CX
C
      N = NY + 1
      IXY = 2 
      CALL RAR(XBIS,N,IXY,IDX,IDY,AR,I16,0)
      L = NY
      V = 0.0D0 
      DO I17 = 1,L 
         V = V + BY(I17)*XBIS(I17)
      ENDDO
      ANU1 = V + CY 
      DKSI = AKSI1 - AKSI 
      DNU = ANU1 - ANU
      DDIST = DSQRT(DKSI**2 + DNU**2) 
      DCRIT = 1.0D-6
      IF (DDIST.LE.DCRIT ) GOTO 1000
C 
C100  N = NX + 1
      N = NX + 1
      IXY = 1 
      CALL RAR(XBIS,N,IXY,IDX,IDY,AR,1,1)
      L = NX
      V = 0.0D0 
      DO I17 = 1,L
         V = V + BX(I17)*XBIS(I17)
      ENDDO
      DKDX = V
C
      N = NX + 1
      IXY = 1 
      CALL RAR(XBIS,N,IXY,IDX,IDY,AR,1,2)
      L = NX
      V = 0.0D0 
      DO I17 = 1,L
         V = V + BX(I17)*XBIS(I17)
      ENDDO
      DKDY = V
C
      N = NY + 1
      IXY = 2 
      CALL RAR(XBIS,N,IXY,IDX,IDY,AR,1,1)
      L = NY
      V = 0.0D0 
      DO I17 = 1,L
         V = V + BY(I17)*XBIS(I17)
      ENDDO
      DNDX = V
C
      N = NY + 1
      IXY = 2 
      CALL RAR(XBIS,N,IXY,IDX,IDY,AR,1,2)
      L = NY
      V = 0.0D0 
      DO I17 = 1,L
         V = V + BY(I17)*XBIS(I17)
      ENDDO
      DNDY = V
C
      YEST(1,2) = (DKSI*DNDX-DNU*DKDX)/(DNDX*DKDY-DKDX*DNDY)
      XEST(1,2) = (DKSI*DNDY-DNU*DKDY)/(DKDX*DNDY-DNDX*DKDY)
      XEST(1,1) = XEST(1,1) - XEST(1,2) 
      YEST(1,1) = YEST(1,1) - YEST(1,2) 
      GOTO 30
1000  XX = -XEST(2,1)
      YY = YEST(2,1)
      RETURN
      END 
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C 
C  The subroutine RAR comes directly from the POS package (cfr Astrom. prgs.) 
C 
C----------------------------------------------------------------------------
      SUBROUTINE RAR(XBIS,N,IXY,IDX,IDY,AR,I,K)
C
C ***
C
      IMPLICIT NONE
C
C ***
C
      INTEGER IDX(1),IDY(1)
      INTEGER K,I,I1,IXY,J,N
C
C ***
C
      DOUBLE PRECISION 	  X1(10) 
      DOUBLE PRECISION    XBIS(10),AR(50,4)
C
C ***
C
      CHARACTER STRING*80
C
C *** do the work
C
      IF (K.EQ.0) GOTO 1 
      IF (K.EQ.1) GOTO 2 
      IF (K.EQ.2) GOTO 3 
C
1     X1(1) = AR(I,1) 
      X1(2) = AR(I,2) 
      X1(3) = AR(I,1)*AR(I,2) 
      X1(4) = AR(I,1)**2
      X1(5) = AR(I,2)**2
      X1(6) = AR(I,1)**3
      X1(7) = AR(I,2)**3
      X1(8) = AR(I,1)*AR(I,2)**2
      X1(9) = AR(I,2)*AR(I,1)**2
      GOTO 5
C 
2     X1(1) = 1.0D0 
      X1(2) = 0.0D0 
      X1(3) = AR(I,2) 
      X1(4) = 2.0D0*AR(I,1) 
      X1(5) = 0.0D0 
      X1(6) = 3.0D0*AR(I,1)**2
      X1(7) = 0.0D0 
      X1(8) = AR(I,2)**2
      X1(9) = 2.0D0*AR(I,2)*AR(I,1) 
      GOTO 5
C 
3     X1(1) = 0.0D0 
      X1(2) = 1.0D0 
      X1(3) = AR(I,1) 
      X1(4) = 0.0D0 
      X1(5) = 2.0D0*AR(I,2) 
      X1(6) = 0.0D0 
      X1(7) = 3.0D0*AR(I,2)**2
      X1(8) = 2.0D0*AR(I,1)*AR(I,2) 
      X1(9) = AR(I,1)**2
C
5     IF(IXY.EQ.1) X1(10) = AR(I,3) 
      IF(IXY.EQ.2) X1(10) = AR(I,4) 
      IF(IXY.EQ.2) GOTO 20
      I1 = 0
      DO J=1,9 
         IF(IDX(J).EQ.0) GOTO 10 
         I1 = I1+1 
         XBIS(I1) = X1(J) 
10    ENDDO
      IF (N-I1.NE.1) GOTO 50 
      XBIS(I1+1) = X1(10)
      GOTO 40
20    I1= 0 
      DO J=1,9 
         IF (IDY(J).EQ.0) GOTO 30 
         I1 = I1+1 
         XBIS(I1) = X1(J) 
30    ENDDO
      IF(N-I1.NE.1) GOTO 50
      XBIS(I1+1) = X1(10)
40    RETURN
50    STRING = '*** FATAL: Coeff. trouble'
      CALL STETER(9,STRING)
      RETURN
      END 
C---------------------------------------------------------------------------
