C @(#)spec3b.for	19.1 (ESO-DMD) 02/25/03 13:24:48
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 SPEC3B
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C
C  Program SPEC3B       version 1.00        860906
C                       version 1.01        870609
C                       (Any no. of rows in i/p table).
C                       version 2.00        890424 (pMIDAS)
C.AUTHOR
C  F. Murtagh                  STECF              
C.MODIFICATIONS
C  M. Peron                    IPG           890901  add include files
C  P. Ballester                IPG           910313  Compil. option -u
C                                                    and more
C
C.KEYWORDS
C  Simulated images, test images.
C  Box function taken as PSF - with an option for the delta function.
C
C.PURPOSE
C  Create spectrum of more than 1 line, each a "box" in form.
C
C.OUTPUT
C
C  Keys:  OUT_A/C/1/60       output data array
C         INPUTR/R/1/1       centring relative to pixel centre
C         IN_TAB/C/1/60      table of positions of lines
C         INPUTR/R/2/1       box (PSF) width
C
C-----------------------------------------------------------
C
C
       IMPLICIT NONE
C
       REAL                     RMIN,RMAX
       REAL                     STEPO(3),STARTO(3),CUTS(4)
       REAL                     CENTR,BOX
C 
       DOUBLE PRECISION         DSTEP(3),DSTART(3)
C 
       INTEGER                  NPIXO(3),IMNO
       INTEGER                  NDIMO
       INTEGER                  MADRID,KUN,KNUL,TID
       INTEGER                  IACT,ISTAT,IDIM,IACTV,NC,NR,NS,NAC,NAR
       INTEGER                  KUNIT
C 
       INTEGER*8                JPNTR,IPTR
C 
       CHARACTER*60             OUTIMA,INTAB
       CHARACTER*72             IDENTO
       CHARACTER*80             CUNITO
C
       INCLUDE                 'MID_INCLUDE:ST_DEF.INC/NOLIST'
       COMMON                  /VMR/MADRID(1)
       INCLUDE                 'MID_INCLUDE:ST_DAT.INC/NOLIST'
C
C
C ... get into MIDAS
C
       CALL STSPRO('SPEC3B')
C
C ... get name of output frame
C
        CALL STKRDC('OUT_A',1,1,60,IACT,OUTIMA,KUN,KNUL,ISTAT)
C
C ... get centring.
C
        CALL STKRDR('INPUTR',1,1,IACT,CENTR,KUN,KNUL,ISTAT)
C
C ... map output image
C
       NDIMO = 1
       IDIM  = 660
       NPIXO(1)  = IDIM
       STARTO(1) = 1
       STEPO(1)  = 1.0
       CUNITO    = ' NONE'
       IDENTO    = ' ARTIFICIAL SPECTRUM' 
        DSTART(1) = STARTO(1)
        DSTEP(1)  = STEPO(1)
        CALL STIPUT(OUTIMA,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .                  NDIMO,NPIXO,DSTART,
     .                  DSTEP,IDENTO,CUNITO,JPNTR,IMNO,ISTAT)
C
C ... get line positions from table.
C
        CALL STKRDC('IN_TAB',1,1,60,IACTV,INTAB,KUN,KNUL,ISTAT)
        CALL TBTOPN(INTAB,0,TID,ISTAT)
        CALL TBIGET(TID,NC,NR,NS,NAC,NAR,ISTAT)
       IF (NC.NE.1.OR.NR.LT.1) THEN
          CALL STTPUT(' Invalid nos. of rows or columns.',ISTAT)
          CALL STTPUT(' in table specifying the lines.',ISTAT)
          CALL STTPUT(' Aborting.',ISTAT)
          GOTO 1000
       ENDIF
        CALL TBCMAP(TID,1,IPTR,ISTAT)
C
C ... get box width.
C
        CALL STKRDR('INPUTR',2,1,IACT,BOX,KUN,KNUL,ISTAT)
C
C ... now do the work
C
       IF (BOX.GT.0.001) THEN      ! PSF is a box of given width.
            CALL PATTERN1(MADRID(JPNTR),RMIN,RMAX,NR,MADRID(IPTR),
     +                    CENTR,BOX,ISTAT)
       ELSE                        ! PSF is a delta function.
            CALL PATTERN2(MADRID(JPNTR),RMIN,RMAX,NR,MADRID(IPTR),
     +                    CENTR,ISTAT)
       ENDIF
C
C ... write cuts
C
       CUTS(1) = RMIN
       CUTS(2) = RMAX
       CUTS(3) = RMIN
       CUTS(4) = RMAX
       CALL STDWRR(IMNO,'LHCUTS',CUTS,1,4,KUNIT,ISTAT)
C
C ... end
C
1000   CONTINUE
       CALL STFCLO(IMNO,ISTAT)
       CALL STSEPI
       END

C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C
C  Program PATTERN1    
C  F. MURTAGH                                 ST-ECF  Version 1.0 860411
C
C.KEYWORDS
C
C  Test patterns, simulated images.
C  A box function is taken as PSF.
C
C
C.OUTPUT PARAMETERS
C 
C  ARR  = the frame,
C  RMIN, RMAX = cut values (max. and min. flux values).
C
C----------------------------------------------------------------------
       SUBROUTINE PATTERN1(ARR,RMIN,RMAX,N,VALS,CENTR,BOX,
     +                          ISTAT)
C
       IMPLICIT NONE
C
C
       INTEGER           N,ISTAT,I,II,III,IVAL,IBEG,IEND
       REAL              ARR(660),VALS(N),VAL,XL,XU
       REAL              RMIN,RMAX,CENTR,BOX,X,ZL,ZU,V
C
C
C
       DO I=1,660
          ARR(I) = 0.0
       ENDDO
C
C
       DO II = 1,N
       IVAL = VALS(II)           
C Just a check on integer values,...
C ...since this is impt. for pixel addressing.
       VAL = IVAL
C
C ...   a few calculations:
C
C ...   XL and XU are extremities of the box:
C
       XL = VAL - 0.5*BOX
       XU = VAL + 0.5*BOX
       XL = MAX(0.5,XL)
       XU = MIN(660.5,XU)
C
C
          DO I =1, 660      ! Do, crudely, for all pixels.
             X = I
             ZL = X-0.5     ! Lower limit of pixel,
             ZU = X+0.5     ! upper limit of pixel.
             IF (XL.GE.ZL.AND.XL.LE.ZU) THEN  ! Lr. end of box is in this pix.
                IF (XU.LE.ZU) THEN            ! All of box is in same pixel.
                   ARR(I) = ARR(I) + XU*10000.0-XL*10000.0
                  GOTO 350 
                ENDIF
                ARR(I) = ARR(I) + ZU*10000.0-XL*10000.0
                IBEG = I+1
                GOTO 200
             ENDIF
          ENDDO
  200          CONTINUE
          DO I = IBEG-1,660
             X = I
             ZL = X-0.5
             ZU = X+0.5
             IF (XU.GT.ZL.AND.XU.LE.ZU) THEN
                ARR(I) = ARR(I) + XU*10000.0-ZL*10000.0             
                IEND = I-1
                GOTO 300
             ENDIF
          ENDDO
  300          CONTINUE
          IF (IEND.GE.IBEG) THEN
             DO III = IBEG,IEND
                ARR(III) = ARR(III) + 10000.0
             ENDDO
          ENDIF
  350      CONTINUE
C
C ...   Now offset input value by +0.5
C
C ...   XL and XU are extremities of the box:
C
       V  = VAL + 0.5
       XL = V - 0.5*BOX
       XU = V + 0.5*BOX
       XL = MAX(0.5,XL)
       XU = MIN(660.5,XU)
C
C
          DO I =1, 660      ! Do, crudely, for all pixels.
             X = I
             ZL = X-0.5     ! Lower limit of pixel,
             ZU = X+0.5     ! upper limit of pixel.
             IF (XL.GE.ZL.AND.XL.LE.ZU) THEN  ! Lower end of box is in this pix.
                IF (XU.LE.ZU) THEN            ! All of box is in same pixel.
                  ARR(I+220) = ARR(I+220) + XU*10000.0-XL*10000.0
                  GOTO 550 
                ENDIF
                ARR(I+220) = ARR(I+220) + ZU*10000.0-XL*10000.0
                IBEG = I+1
                GOTO 400
             ENDIF
          ENDDO
  400          CONTINUE
          DO I = IBEG-1,660
             X = I
             ZL = X-0.5
             ZU = X+0.5
             IF (XU.GT.ZL.AND.XU.LE.ZU) THEN
              ARR(I+220) = ARR(I+220) + XU*10000.0-ZL*10000.0             
              IEND = I-1
              GOTO 500
             ENDIF
          ENDDO
  500          CONTINUE
          IF (IEND.GE.IBEG) THEN
             DO III = IBEG,IEND
                ARR(III+220) = ARR(III+220) + 10000.0
             ENDDO
          ENDIF
  550      CONTINUE
C
C ...   Now offset input value by +CENTR
C
C ...   XL and XU are extremities of the box:
C
       V  = VAL + CENTR
       XL = V - 0.5*BOX
       XU = V + 0.5*BOX
       XL = MAX(0.5,XL)
       XU = MIN(660.5,XU)
C
C
          DO I =1, 660      ! Do, crudely, for all pixels.
             X = I
             ZL = X-0.5     ! Lower limit of pixel,
             ZU = X+0.5     ! upper limit of pixel.
             IF (XL.GE.ZL.AND.XL.LE.ZU) THEN  ! Lower end of box is in this pix.
                IF (XU.LE.ZU) THEN            ! All of box is in same pixel.
                   ARR(I+440) = ARR(I+440) + XU*10000.0-XL*10000.0
                  GOTO 750 
                ENDIF
                ARR(I+440) = ARR(I+440) + ZU*10000.0-XL*10000.0
                IBEG = I+1
                GOTO 600
             ENDIF
          ENDDO
  600          CONTINUE
          DO I = IBEG-1,660
             X = I
             ZL = X-0.5
             ZU = X+0.5
             IF (XU.GT.ZL.AND.XU.LE.ZU) THEN
                ARR(I+440) = ARR(I+440) + XU*10000.0-ZL*10000.0             
                IEND = I-1
                GOTO 700
             ENDIF
          ENDDO
  700          CONTINUE
          IF (IEND.GE.IBEG) THEN
             DO III = IBEG,IEND
                ARR(III+440) = ARR(III+440) + 10000.0
             ENDDO
          ENDIF
  750      CONTINUE
C
C ...   that's all:
C
       ENDDO
C
C
C ------DETERMINE CUTS (I.E. MAX AND MIN VALUES)---------------------
C
       RMIN = 1000.0
       RMAX = 0.0
       DO I = 1, 660
          ARR(I) = ARR(I)/100.0
          IF (ARR(I).GT.RMAX) RMAX = ARR(I) 
          IF (ARR(I).LT.RMIN) RMIN = ARR(I) 
       ENDDO
C
C
       RETURN
C
       END

C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C
C  Program PATTERN2    
C  F. MURTAGH                                 ST-ECF  Version 1.0 860411
C
C.KEYWORDS
C
C  Test patterns, simulated images.
C  Delta function taken as PSF of lines.
C
C
C.OUTPUT PARAMETERS
C 
C  ARR  = the frame,
C  RMIN, RMAX = cut values (max. and min. flux values).
C
C----------------------------------------------------------------------
       SUBROUTINE PATTERN2(ARR,RMIN,RMAX,N,VALS,CENTR,
     X                          ISTAT)
C
       IMPLICIT NONE
C
       REAL                     ARR(660),VALS(1)
C
       INTEGER                  N,ISTAT,I
       REAL                     RMIN,RMAX,CENTR
       INTEGER                  LVAL             ! From Real to INT 25.04.91
C
C
C
       DO I = 1, 660
          ARR(I) = 0.0
       ENDDO
C
C ...   first line is located at centre of pixel:
        DO I = 1, N
          LVAL = VALS(I)
          ARR(LVAL) = 100.0
       ENDDO
C ...   next, line is located at pixel boundary - by convention,
C ...   with delta function, all goes to one pixel:
       DO I = 1, N
          LVAL = VALS(I)
          ARR(LVAL+220) = 100.0
       ENDDO
C ...   finally, line is located at offset rel. to pix. centre:
       DO I = 1, N
          LVAL = VALS(I)
          IF (CENTR.GT.0.5.AND.CENTR.LT.1.0) LVAL = LVAL+1
          IF (LVAL+440.LE.660) ARR(LVAL+440) = 100.0
       ENDDO
C
C
C ------DETERMINE CUTS (I.E. MAX AND MIN VALUES)---------------------
C
       RMIN = 0.0
       RMAX = 100.0
C
C
       RETURN
C
       END
