C @(#)creastar.for	19.1 (ES0-DMD) 02/25/03 13:19:04
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
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C.IDENTIFICATION: CRSTAR
C.LANGUAGE:  F77+ESOext
C.AUTHOR:    Rein H. Warmels
C.PURPOSE:   Create an image file pixels containing the profile
C            of the standard reference star, by adding and averaging
C            some selected  stars, after recentering
C.VERSION:   840524 ChO ???
C.VERSION:   931122 RHW rewritten 
C ---------------------------------------------------------------------
        PROGRAM CRSTAR
C
        IMPLICIT NONE
        REAL             PARAM(7)
        REAL             XYZ(3,1000)
        REAL             CUTS(4)
        REAL             DMIN, DMAX
C
        INTEGER          NAXIS
        INTEGER          NPIX1(2), NPIX(2), NPX
        INTEGER          ICOL(3)
        INTEGER          IAC, I, IZONE
        INTEGER          TID,IMF
        INTEGER*8        IP,IPOUT
        INTEGER          IMF1
        INTEGER          KUN,KNUL,ISTAT
        INTEGER          NCOL,NROW,NSC,NAC,NAR
        INTEGER          MADRID(1)
C
        DOUBLE PRECISION START(2), STEP(2)
        DOUBLE PRECISION START1(2), STEP1(2)
C
        CHARACTER*16     LABEL(3)
        CHARACTER*60     FRAMIN,TBIN,FRAMOUT
        CHARACTER        IDENT*72,CUNIT*64
C     
        LOGICAL          NULL(3)
C
        INCLUDE          'MID_INCLUDE:ST_DEF.INC/NOLIST'
        COMMON           /VMR/MADRID
        INCLUDE          'MID_INCLUDE:ST_DAT.INC/NOLIST'
C
        DATA LABEL       /'X_COORD','Y_COORD','VALUE'/
C
C ***   start the code
        CALL STSPRO('CSTAR')
        CALL STKRDC('IN_A',1,1,60,IAC,FRAMIN,KUN,KNUL,ISTAT)
        CALL STKRDC('IN_B',1,1,60,IAC,TBIN,KUN,KNUL,ISTAT)
        CALL STKRDC('OUT_A',1,1,60,IAC,FRAMOUT,KUN,KNUL,ISTAT)
        CALL STKRDI('INPUTI',1,1,IAC,NPX,KUN,KNUL,ISTAT)
        CALL STKRDR('INPUTR',1,7,IAC,PARAM,KUN,KNUL,ISTAT)
C
C ***   open the image
        CALL STIGET(FRAMIN,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,2,NAXIS,
     2              NPIX,START,STEP,IDENT,CUNIT,IP,IMF,ISTAT)
C
C *** prepare the tables
        CALL TBTOPN(TBIN,F_I_MODE,TID,ISTAT)
        CALL TBIGET(TID,NCOL,NROW,NSC,NAC,NAR,ISTAT)               ! read table
C
C *** search of colonne indices by name
        DO I=1,3
           CALL TBLSER(TID,LABEL(I),ICOL(I),ISTAT)
        ENDDO
        DO I = 1,NROW
           CALL TBRRDR(TID,I,3,ICOL,XYZ(1,I),NULL,ISTAT)
        ENDDO
        CALL TBTCLO(TID,ISTAT)
C
C *** open the output frame
        START1(1) = PARAM(1)
        START1(2) = PARAM(2)
        STEP1(1)  = PARAM(3)
        STEP1(2)  = PARAM(4)
        NPIX1(1)  = NPX
        NPIX1(2)  = NPX
        CALL STIPUT(FRAMOUT,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,NAXIS,
     2              NPIX1,START1,STEP1,IDENT,CUNIT,IPOUT,IMF1,ISTAT)
C
C *** get minimum and maximum
        DMIN  = PARAM(5)
        DMAX  = PARAM(6)
        IZONE = INT(PARAM(7))
        IF (DMIN .EQ. DMAX) THEN
            CALL STDRDR(IMF,'LHCUTS',1,4,IAC,CUTS,KUN,KNUL,ISTAT) ! min and max
            DMIN = CUTS(3)
            DMAX = CUTS(4)
        ENDIF
C
C *** do the stellar creation
        CALL SPPHOT(MADRID(IP),NPIX,START,STEP,XYZ,NROW,
     2              NPX,IZONE,DMIN,DMAX,MADRID(IPOUT))
C
C *** finish
        CALL STSEPI
        END




        SUBROUTINE SPPHOT(XIN,NPIX,START,STEP,XYZ,NET,
     *                    NPNT,IZONE,DMIN,DMAX,XOUT)
C+++
C       routine for stellar photometry; generates a stellar image by averaging
C---

        IMPLICIT NONE
 
        INTEGER          NPN, NPN1
        PARAMETER        (NPN  = 1024)
        PARAMETER        (NPN1 = NPN + 1)

        REAL             XIN(*)
        INTEGER          NPIX(*)
        DOUBLE PRECISION START(*),STEP(*)
        REAL             XYZ(3,1000)
        INTEGER          NET
        INTEGER          NPNT
        INTEGER          IZONE
        REAL             DMIN, DMAX
        REAL             XOUT(*)
C            
        INTEGER          NPNT1
        INTEGER          J, LL
        INTEGER          NETOT
        INTEGER          I, IE
        INTEGER          IVAL, ISTAT
        INTEGER          KL, KLT, KP, KTOT
        INTEGER          IP, IPMIN, IPMAX
        INTEGER          IL, ILMIN, ILMAX
        INTEGER          IN1, IN2
        INTEGER          NPL, NL, NP
        REAL             AIL, AIP
        REAL             D, DPHI
        REAL             ETOT
        REAL             X0, Y0, XD, YD
        REAL             DX, DY
        REAL             DP, DL
        REAL             P1, P2, P3, P4
        REAL             RMES(NPN1,2),RMESS(NPN)
        REAL             SMOY(NPN,NPN), SMOYB
        REAL             TOT
        CHARACTER        TEXT*80
C
 1100   FORMAT('Number of stars:',I5)
C
        NPNT1 = NPNT + 1
        DX    = STEP(1)
        DY    = STEP(2)
        NPL   = NPIX(1)
        NL    = NPIX(2)
        X0    = START(1)
        Y0    = START(2)
        IVAL  = IZONE/2
C
C *** 
        NETOT = 0
        DO 40 KL = 1,NPNT
           DO 41 KP = 1,NPNT
              SMOY(KP,KL) = 0.
   41      CONTINUE
   40   CONTINUE            
C
        DO 100 IE = 1,NET
           XD   = XYZ(1,IE)
           YD   = XYZ(2,IE)
           D    = XYZ(3,IE)
           DPHI = D
           IF (DPHI.GT.DMAX.OR.DPHI.LT.DMIN) THEN
              GO TO 100
           ENDIF
C
C ***   calculate zone limits
           AIP   = (XD-X0)/DX+1.
           AIL   = (YD-Y0)/DY+1.
           IP    = AIP
           IL    = AIL
           IPMIN = IP-IVAL
           IPMAX = IP+IVAL
           ILMIN = IL-IVAL
           ILMAX = IL+IVAL
C
C *** exclusion of bordwer points
           IF (IPMIN.LT.1) GO TO 100
           IF (IPMAX.GT.NPL) GO TO 100
           IF (ILMIN.LT.1) GO TO 100
           IF (ILMAX.GT.NL) GO TO 100
           DP = AIP-IP
           DL = AIL-IL
C
C *** interpolation 
           P1 = (1.-DP)*(1.-DL)
           P2 = DP*(1.-DL)
           P3 = (1.-DP)*DL
           P4 = DP*DL
           J  = 1
           CALL LIRE(ILMIN,NPL,IPMIN,IPMAX,1,XIN,RMES(1,J))
           KTOT = IVAL*2
           KL   = 0
C
C ** added after interpolation
           DO 50 LL=ILMIN+1,ILMAX
              KL= KL+1
              I = J
              J = MOD(J,2)+1
              CALL LIRE(LL,NPL,IPMIN,IPMAX,1,XIN,RMES(1,J))
              DO 51 KP = 1,KTOT
                 SMOY(KP,KL) = SMOY(KP,KL) +
     2                         P1*RMES(KP,I) + P2*RMES(KP+1,I) +
     3                         P3*RMES(KP,J) + P4*RMES(KP+1,J)
   51         CONTINUE
   50      CONTINUE
           NETOT = NETOT+1 
  100    CONTINUE
C
C *** 
        WRITE (TEXT,1100) NETOT
        CALL STTPUT(TEXT,ISTAT)
        IF (NETOT.GT.0) GO TO 110
           CALL STTPUT('*** WARNING: Stars are beyound limits',ISTAT)
           RETURN
  110   CONTINUE
C
C *** normalize
        ETOT = FLOAT(NETOT)
        DO 120 KLT = 1,KL
           DO 121 KP = 1,KTOT
              SMOY(KP,KLT) = SMOY(KP,KLT)/ETOT
  121      CONTINUE
  120   CONTINUE
C
C *** prepare output table; calculate the pixel values
        SMOYB = 0.                
        DO 401 I = 1,KTOT
           SMOYB = SMOYB+SMOY(I,1)
           SMOYB = SMOYB+SMOY(I,KL)
  401   CONTINUE
        DO 402 I = 2,KL-1
           SMOYB = SMOYB+SMOY(1,I)
           SMOYB = SMOYB+SMOY(KTOT,I)
  402   CONTINUE
        TOT   = KTOT*2+(KL-2)*2
        SMOYB = SMOYB/TOT
        IN1   = INT(NPNT/2) - 1 - IVAL
        IN2   = IN1+IZONE-1
        IF (IN1.EQ.0) THEN
           IN2 = NPNT
           GO TO 2
        END IF
C
C *** write the otput tabler after re-centering
        DO 291 KP = 1,NPNT
           RMESS(KP) = SMOYB
  291   CONTINUE
        DO 301 NL = 1,IN1
           CALL ECRIT(NL,NPNT,RMESS,XOUT)
  301   CONTINUE
C
        J=1
    2   CONTINUE
        DO 302 NL = IN1+1,IN2
           I=1
           DO 292 NP = 1,IN1
              RMESS(NP) = SMOYB
  292      CONTINUE
           DO 293 NP = IN1+1,IN2
              RMESS(NP) = SMOY(I,J)
              I=I+1
  293      CONTINUE
           DO 294 NP = IN2+1,NPNT
              RMESS(NP) = SMOYB
  294      CONTINUE
           CALL ECRIT(NL,NPNT,RMESS,XOUT)
           J = J+1
  302   CONTINUE
C         
        IF (IN2.EQ.NPNT) THEN
           RETURN
        ENDIF
        DO 295 KP=1,NPNT
           RMESS(KP) = SMOYB
  295   CONTINUE        
        DO 303 NL = IN2+1,NPNT
           CALL ECRIT(NL,NPNT,RMESS,XOUT)
  303   CONTINUE
C
        RETURN
        END


        SUBROUTINE ECRIT(NLI,NPOINT,RMES,FMES)
C
        IMPLICIT NONE
        INTEGER  NLI, NPOINT
        REAL     RMES(1),FMES(1)
C
        INTEGER  KI, I
C
        DO I = 1,NPOINT
              KI = I + (NLI-1)*NPOINT
              FMES(KI) = RMES(I)
        ENDDO
C
        RETURN
        END



        SUBROUTINE LIRE(NL,NPL,NPL1,NPL2,NPL3,FMES,RMES)
C+++
C.PURPOSE:  Write part of a frame into an array
C.AUTHOR:   ???
C.VERSION:  ?????? ??? created
C.VERSION:  890117 RHW documented
C.COMMENTS: none
C---
        IMPLICIT   NONE
        INTEGER    NL
        INTEGER    NPL
        INTEGER    NPL1
        INTEGER    NPL2
        INTEGER    NPL3
        REAL       FMES(*)
        REAL       RMES(*)                         

        INTEGER    NPD, NPF, K, I 

C ***
        NPD = NPL* (NL-1) + NPL1    
        NPF = NPD + NPL2 - NPL1
        K   = 0 
C 
        DO 10 I = NPD,NPF,NPL3
           K       = K + 1  
           RMES(K) = FMES(I)
   10   CONTINUE 
C       RETURN
        END

