C @(#)intfr.for	19.1 (ES0-DMD) 02/25/03 13:19:39
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
      SUBROUTINE INTFRI(VAL,NPL,IMAGE,SX,XS,NCUR,NDEG)             
C+++
C---
      IMPLICIT      NONE
      INTEGER       NPL
      REAL          VAL(NPL)
      REAL          IMAGE(4)          
      REAL          SX
      REAL          XS
      INTEGER       NCUR
      INTEGER       NDEG
C
      INTEGER       MAXC
      INTEGER       MAXT
      PARAMETER     (MAXC=100)
      PARAMETER     (MAXT=100000)
      INTEGER       IBAR
      INTEGER       BNMOD
      INTEGER       IAC
      INTEGER       ISTAT
      INTEGER       I
      INTEGER       INNER
      INTEGER       KA
      INTEGER       KUN
      INTEGER       LINET2, LINET3
      INTEGER       MADRID(1)
      INTEGER       NPL1,NPL2,NPL3 
      INTEGER       NP1, NP2, NP3, NP4
      INTEGER       NPIX
      INTEGER       NPT
      REAL          AEQW
      REAL          ALICO3, ACONT3
      REAL          ALICO4, ACONT4
      REAL          ALINE, ALICO, ACONT 
      REAL          ALINT, AITKEN
      REAL          FLXFRC
      REAL          OUTP(10)  
      REAL          RC3, XC3, DX3
      REAL          RC4, XC4, DX4
      REAL          SX3, SX4
      REAL          XARG
      REAL          XCUR, YCUR
      REAL          XRANG
      REAL          XPC(MAXC),YPC(MAXC) 
      REAL          X(MAXT),Y(MAXT)
      REAL          WORK(MAXT)
      CHARACTER*80  TEXT,HEADER 
      CHARACTER     BIN*4
      LOGICAL       FIRST 
      INCLUDE       'MID_INCLUDE:PLTDEC.INC/NOLIST' 
      COMMON        /VMR/MADRID 
C
      DATA IBAR/1/
      DATA LINET2/2/
      DATA LINET3/3/
C
 9010 FORMAT(1X,5G13.6)
C
C *** start of executable code **********************************************
C     Now do the looping for cursers and interpolations
C     NCUR positions are asked for
      NPL1  = INT(MIN(IMAGE(1),IMAGE(2))) 
      NPL2  = INT(MAX(IMAGE(1),IMAGE(2))) 
      NPL3  = 1
      NPT   = ABS(NPL1-NPL2) + 1
C
C *** Now do the looping for cursers and integration
C     Actions: NCUR cursors are asked for
C     A space bar terminates this loop
      FIRST = .TRUE.

   30 CONTINUE
      DO 50 I = 1,NCUR
         CALL PTGCUR(XCUR,YCUR,KA,ISTAT)
         IF (KA.EQ.32) THEN                                 ! exit on space bar
            GO TO 120
         ELSE                                 ! correct input; store data point
            XPC(I) = XCUR
            YPC(I) = YCUR
            CALL PTDATA(5,0,0,XPC(I),YPC(I),0.0,1)
         END IF 
   50 CONTINUE
      CALL SORT2(XPC,YPC,NCUR)
C
C *** Determine the range effective for change (= the EXTREME cursors) 
      NP1 = NINT((XPC(1)-XS)/SX) + 1             ! pixel number corr. to cursor
      NP2 = NINT((XPC(NCUR)-XS)/SX) + 1
 
      IF (NP1.LT.NPL1) THEN
         NP1 = NPL1
      ENDIF
C
C *** Do the interpolation here and plot the line 
      IF (NDEG.LT.2) THEN
         DO 70 I = NP1,NP2
            XARG    = XS + FLOAT(I-1)*SX
            WORK(I) = ALINT(XPC,YPC,NCUR,XARG,1,NCUR) 
   70    CONTINUE
      ELSE

         DO 80 I = NP1,NP2 
            XARG    = XS + FLOAT(I-1)*SX
            WORK(I) = AITKEN(XPC,YPC,NCUR,XARG,1,NCUR,NDEG)               
   80    CONTINUE
      END IF
C
C *** get the bin option, plot and write the results
      CALL PTKRDC('BINMOD',4,IAC,BIN)
      IF (BIN(1:2).EQ.'ON') THEN
         BNMOD = 1
      ELSE
         BNMOD = 0
      ENDIF
      NPIX = NP2 - NP1 + 1
      DO 90 I = NP1,NP2
         X(I) = XS + FLOAT(I-1)*SX
         Y(I) = WORK(I)
   90 CONTINUE 
      CALL PTDATA(0,LINET3,BNMOD,X(NP1),Y(NP1),0.0,NPIX)
C
C *** Do the integration
C     FIRST compute the flux density in the first bin.
C     The bin limits are the input cursor position and the middle between
C     the first and second image pixel. We compute the interpolated values
C     at these points and make the integral in the first bin.
      INNER  = NCUR/2                               ! first inner pixel
      XC3    = XPC(INNER)                           ! real cursor position
      RC3    = (XC3-XS)/SX+1.0                      ! real cursor pixel
      NP3    = NINT((XC3-XS)/SX)+1                  ! nearest image pixel
      DX3    = FLOAT(NP3)-RC3
      SX3    = (DX3+0.5)*SX
      ACONT3 = SX3*WORK(NP3)
      ALICO3 = SX3*VAL(NP3)
C
C *** SECOND compute the flux density in the last bin.
C     The bin limits are the last cursor position and the middle between
C     the last and one but last image pixel. We compute the interpolated 
C     values at these points and make the integral in the last bin.
C     and then do the integration.
      XC4    = XPC(INNER+1)
      RC4    = (XC4-XS)/SX+1.0
      NP4    = NINT((XC4-XS)/SX)+1
      DX4    = RC4-FLOAT(NP4)
      SX4    = (DX4+0.5)*SX 
      ACONT4 = SX4*WORK(NP4)
      ALICO4 = SX4*VAL(NP4)
C
C *** do the pixel summation for continuum and line
      ACONT  = 0.
      ALICO  = 0.
      DO 110 I = NP3+1,NP4-1 
          ACONT = WORK(I)*SX + ACONT                                          
          ALICO = VAL(I)*SX + ALICO
  110 CONTINUE                  
      ACONT    = ACONT + ACONT3 + ACONT4
      ALICO    = ALICO + ALICO3 + ALICO4
      XRANG    = XPC(INNER+1)-XPC(INNER)
      ALINE    = ALICO - ACONT                                      
      AEQW     = -1.0*(ALINE/ACONT)*XRANG                          
      FLXFRC   = ALINE/ACONT
C ***
C
      IF (FIRST) THEN
          HEADER = '   X_start (pix/world)          '//
     2             ' X_end (pix/world)      Pixel sep. '
          CALL STTPUT(HEADER,ISTAT)
          HEADER = '  Line+Cont.   Continuum        '//
     2             'Line       Line/Cont    Equiv. w.  '
          CALL STTPUT(HEADER,ISTAT)
          HEADER = ' ---------------------------------------'//
     2             '-------------------------'
          CALL STTPUT(HEADER,ISTAT) 
          FIRST  = .FALSE.
      END IF 
      WRITE(TEXT,9010) XC3,RC3,XC4,RC4,SX
      CALL STTPUT(TEXT,ISTAT)

      WRITE(TEXT,9010) ALICO,ACONT,ALINE,FLXFRC,AEQW
      CALL STTPUT(TEXT,ISTAT)
C ***
C
      OUTP(1)  = RC3                                      ! first NUMBER
      OUTP(2)  = RC4                                      ! last pixel
      OUTP(3)  = XC3                                      ! first cursor 
      OUTP(4)  = XC4                                      ! last cursor 
      OUTP(5)  = SX                                       ! step in x
      OUTP(6)  = ALICO                                    ! line + cont
      OUTP(7)  = ACONT                                    ! continuum
      OUTP(8)  = ALINE                                    ! line
      OUTP(9)  = FLXFRC                                   ! flux fraction
      OUTP(10) = AEQW                                     ! equivalent width
      CALL STKWRR('OUTPUTR',OUTP,1,10,KUN,ISTAT)
      GO TO 30 

  120 CONTINUE
      RETURN
      END


      SUBROUTINE INTFRB(VAL,NPL,IMAGE,SX,XS,NCUR,NDEG,             
     2                  DELTAX,XPC)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.AUTHOR:   F. Murtagh, ST-ECF                                                 
C.VERSION:  860321 FM Creation
C.VERSION:  880607 RHW Restructured and recoded for portable MIDAS
C.VERSION:  900430 RHW Make code similar to INTFRI 
C.VERSION:  900430 RHW Store more output in the keyword OUTPUTR
C.VERSION:  900629 RHW Integration over subpixel included
C.VERSION:  910115 RHW IMPLICIT NONE added
C----------------------------------------------------------------------------
      IMPLICIT  NONE
      INTEGER   MAXC
      PARAMETER (MAXC=100)
C
      INTEGER   NPL
      REAL      VAL(NPL)
      REAL      IMAGE(4)          
      REAL      SX
      REAL      XS
      INTEGER   NCUR
      INTEGER   NDEG
      REAL      DELTAX
      REAL      XPC(MAXC)
C
      INTEGER   MAXT
      PARAMETER (MAXT=100000)
      INTEGER   ISTAT
      INTEGER   KUN
      INTEGER   I, I2, INNER
      INTEGER   NP1, NP2, NP3, NP4
      INTEGER   NPL1, NPL2, NPL3, NPT
      REAL      AITKEN, ALINE, ALINT
      REAL      ACONT, ALICO
      REAL      ACONT3, ACONT4, ALICO3, ALICO4
      REAL      AEQW
      REAL      DX3, DX4
      REAL      FLXFRC
      REAL      OUTP(10)                                            
      REAL      RC3, RC4
      REAL      SX3, SX4
      REAL      WORK(MAXT)
      REAL      XARG
      REAL      Y
      REAL      XC3, XC4
      REAL      XRANG
      REAL      YPC(MAXC) 
C
C *** start of executable code **********************************************
C     Now do the looping for cursers and interpolations.
C     NCUR positions are asked for.
      NPL1  = INT(MIN(IMAGE(1),IMAGE(2)))
      NPL2  = INT(MAX(IMAGE(1),IMAGE(2))) 
      NPL3  = 1 
      NPT   = ABS(NPL1-NPL2) + 1 
C
C *** loop through cursor position and determine nearest pixels
      DO 10 I = 1, NCUR 
         NP1 = NINT((XPC(I)-DELTAX-XS)/SX) + 1                             
         NP2 = NINT((XPC(I)+DELTAX-XS)/SX) + 1                             
         IF (NP1.LT.NPL1 .OR. NP1.GT.NPL2 .OR.
     2       NP2.LT.NPL1 .OR. NP2.GT.NPL2) THEN
             CALL STTPUT('*** FATAL: Cursor position(s) outside'//
     2                   ' frame; try again', ISTAT)
             CALL STSEPI
         ENDIF
         Y   = 0.0
         DO 20 I2 = NP1, NP2 
            Y = Y + VAL(I2)
   20    CONTINUE
         Y     = Y/FLOAT(NP2-NP1+1)
         YPC(I) = Y 
   10 CONTINUE 
      CALL SORT2(XPC,YPC,NCUR)
C
C *** Determine effective range of cursor positions
      NP1 = NINT((XPC(1)-XS)/SX) + 1                    ! first cursor position
      NP2 = NINT((XPC(NCUR)-XS)/SX) + 1                 ! last cursor position
C
C *** Do the interpolation here
      IF (NDEG.LT.2) THEN                              !compute the continuum
         DO 30 I = NP1,NP2
            XARG    = XS + FLOAT(I-1)*SX
            WORK(I) = ALINT(XPC,YPC,NCUR,XARG,1,NCUR)
   30    CONTINUE
     
       ELSE
         DO 40 I = NP1,NP2
            XARG    = XS + FLOAT(I-1)*SX             
            WORK(I) = AITKEN(XPC,YPC,NCUR,XARG,1,NCUR,NDEG)
   40    CONTINUE
      ENDIF
C
C *** Do the integration 
C     FIRST compute the flux density in the first bin.
C     The bin limits are the input cursor position and the middle between
C     the first and second image pixel. We compute the interpolated values
C     at these points and make the integral in the first bin.
      INNER  = NCUR/2                               ! first inner pixel
      XC3    = XPC(INNER)                           ! real cursor position
      RC3    = (XC3-XS)/SX+1.0                      ! real cursor pixel
      NP3    = NINT((XC3-XS)/SX)+1                  ! nearest image pixel
      DX3    = FLOAT(NP3)-RC3
      SX3    = (DX3+0.5)*SX
      ACONT3 = SX3*WORK(NP3)
      ALICO3 = SX3*VAL(NP3)
C
C *** SECOND compute the flux density in the last bin.
C     The bin limits are the last cursor position and the middle between
C     the last and one but last image pixel. We compute the interpolated 
C     values at these points and make the integral in the last bin.
C     and then do the integration.
      XC4    = XPC(INNER+1)
      RC4    = (XC4-XS)/SX+1.0
      NP4    = NINT((XC4-XS)/SX)+1
      DX4    = RC4-FLOAT(NP4)
      SX4    = (DX4+0.5)*SX 
      ACONT4 = SX4*WORK(NP4)
      ALICO4 = SX4*VAL(NP4)
C
C *** do the pixel summation for continuum and line
      ACONT  = 0.
      ALICO  = 0.
      DO 50 I = NP3+1,NP4-1                                                    
          ACONT = WORK(I)*SX + ACONT                                          
          ALICO = VAL(I)*SX + ALICO
   50 CONTINUE                  
      ACONT    = ACONT + ACONT3 + ACONT4
      ALICO    = ALICO + ALICO3 + ALICO4
      XRANG    = XPC(INNER+1)-XPC(INNER)
      ALINE    = ALICO - ACONT                                      
      AEQW     = -1.0*(ALINE/ACONT)*XRANG                          
      FLXFRC   = ALINE/ACONT
C ***
C
      OUTP(1)  = RC3                                      ! first NUMBER
      OUTP(2)  = RC4                                      ! last pixel
      OUTP(3)  = XC3                                      ! first cursor
      OUTP(4)  = XC4                                      ! last cursor 
      OUTP(5)  = SX                                       ! step in x
      OUTP(6)  = ALICO                                    ! line + cont
      OUTP(7)  = ACONT                                    ! continuum
      OUTP(8)  = ALINE                                    ! line
      OUTP(9)  = FLXFRC                                   ! flux fraction
      OUTP(10) = AEQW                                     ! equivalent width
      CALL STKWRR('OUTPUTR',OUTP,1,10,KUN,ISTAT)
C
      RETURN
      END 
