C===========================================================================
C Copyright (C) 1995,2003 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
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987,2003 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  18:08 - 11 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.IDENTIFICATION
C    TDREBIN.FOR
C
C
C.AUTHOR: M.ROSA
C REBPOL modified 900117 M Peron
C remove call to NAG 900903 M Peron
C 
C 030925		last modif
C--------------------------------------------------------------------

      SUBROUTINE REBMET(NIN,XI,YI,WI,NOUT,XO,WO,IRFUN,NPAR,NPACT,DPAR,
     +                  IOPT,MINT,YO,RMIN,RMAX)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.KEYWORDS
C        NonLINEAR rebin
C
C.PURPOSE
C        NonLINEAR rebin
C
C.ALGORITHM
C
C        TBF
C
C.USAGE
C        CALL REBMET(NIN,XI,YI,WI,NOUT,XO,WO,IRFUN,NPAR,NPACT,DPAR,
C     IOPT,MINT,YO,RMIN,RMAX)
C
C.PARAMETERS
C        Inputs:
C   NIN,XI,YI,WI    no. of (input) points, indep, dep and bin inputs
C   NOUT,XO,WO      no, of (output) points, indep, bin outputs
C   IRFUN,NPAR,NPACT,PAR  func. no., size par., act no par., parameters
C
C        Note, function relates XO --> XI  i.e.  xi(f(xo))  !!!!!!
C
C        IRFUN = 1: linear    x (old) = A + Bx
C              = 2: poly.     x       = A + Bx + Cxx + (<= 16 coeffs)
C              = 3: inv.      x       = A + B/x
C              = 4: exp.      x       = A + B*EXP(Cx)
C              = 5: dexp.     x       = A + B*10.**(Cx)
C              = 6: log       x       = A + B*LOG(Cx)
C              = 7: dlog      x       = A + B*LOG10(Cx)
C              = 8: inv.pol.  x       = A + B/(x-x0) + C/(x-x0)**2 + ...
C              = 9: user defined function
C
C        IOPT,MINT     interpol/integrate option, no. of
C                                      points to be interpolated
C        IOPT = 1: spline interpolation and gaussian integration (most
C     accurate but most time consuming option)
C             = 2: intensity is simply relevant portion of pixel's
C                  intensity; for heavy oversampling of input, use this.
C             = 3: linear interpolation, followed by staightforward
C                  integration (use if i/p and o/p binwidths similar).
C
C        Outputs:
C        YO,RMIN,RMAX   dep and maxmin of output
C
C.RESTRICTIONS
C
C.MODIFICATIONS
C
C.LONG DESCRIPTION
C       The subroutines below are well commented.  For background
C       information on the use of this routine, see the Users Guide
C       for the modelling system in USG.MEM in [STMODEL.DOC]
C       (available through M. Rosa, ST-ECF).
C
C--------------------------------------------------------------------
C
      INTEGER NIN,NOUT,IRFUN,NPAR,NPACT,IOPT,MINT
      INTEGER NA,NZ,I,J,LOP,LIP,IN,KO,KF,LL,LH,ML,MH
C 
      DOUBLE PRECISION XI(NIN),YI(NIN),WI(NIN)
      DOUBLE PRECISION XO(NOUT),WO(NOUT),DPAR(NPAR),PL,PH
      DOUBLE PRECISION XL,XH,XW,XINCR,XVLS,DLOP,DLIP,DL,DH
C 
      REAL YO(NOUT),RMIN,RMAX,VAL,VALI
C 
C 
C 
C     No. of intervals
      XVLS   = DBLE(MINT-1)
C
C  We start here doing some brain work in order to speed up
C  the program --- all this should be ok but who knows (MR/19/4/85)
C
C  There are 4 cases for the indep. arrays
C XI(i) < XI(i+1)
C XI(i) > XI(i+1)
C XO(i) < XO(i+1)
C XO(i) > XO(i+1)
C
C  and two for the mapping of XO onto XI
C       XO small onto XI small
C       XO small onto XI large
C
C  lets make a few tests, so we can optimize the DO LOOPS
C       XO increasing
      IF (XO(1).LE.XO(2)) THEN
          KO     = 1
          XL     = XO(1)
          XH     = XO(NOUT)
      ELSE
          KO     = -1
          XH     = XO(1)
          XL     = XO(NOUT)
      END IF

      CALL REBFNC(IRFUN,XL,NPAR,NPACT,DPAR)
C  the function is inverting direction
      CALL REBFNC(IRFUN,XH,NPAR,NPACT,DPAR)  
C  1: normal, -1: inverted direction
      KF     = 1  
      IF (XL.GT.XH) KF     = -1
      LOP    = KF*KO
C
C  So we set up the loops in such a way that the input indep. will be
C  increasing for all cases. That means we can keep track of search
C  limits (large gain in search for long arrays).
C
C  REBIN with flux conservation means that we have to integrate the flux
C  in the input data between the projected flanks of our new pixel.
C   Remember, that the binwidth is negative if the arrays are
C  monotonicaly decreasing !!! Since we will have to look for the
C  projection of the flanks of each pixel XO in the XI array (XL,XR),
C  where Left and Right referes to the direction of indexing in XI,
C  we have to be careful with the directions (4 cases):
C
C  XO incr, LOP pos, WO pos, xl = xo - 0.5*wo , xr = xl + wo
C  XO incr, LOP neg, WO pos, xl = xo + 0.5*wo , xr = xl - wo
C  XO decr, LOP pos, WO neg, xl = xo + 0.5*wo , xr = xl - wo
C  XO decr, LOP neg, WO neg, xl = xo - 0.5*wo , xr = xl + wo
C
C  Hence, XL = XO - 0.5*LOP*WO, XR = XL + LOP*WO  (in REAL*8)
C
C  The new indep. may not map completely onto the old and vice versa.
C  We will NOT extrapolate and we will fill undefined output with
C  NULLVALUE. Now we need
C  to obtain the applicable limits in XO, the direction in XI and the
C  search boundaries for XI and fill YO with 0 outside, in order to
C  avoid all sorts of limit testing in the actual run through XO.
C     XI increasing  !        decreasing
      LIP    = 1
      IF (XI(1).GT.XI(2)) LIP    = -1
      DLIP   = 0.5D0*DBLE(LIP)
      DLOP   = 0.5D0*DBLE(LOP)
C
C Check for the active area in XO
C
      IF (LOP.GT.0) THEN
          NA     = 1
          NZ     = NOUT
      ELSE
          NA     = NOUT
          NZ     = 1
      ENDIF
C
      DL = XI(1)
      DH = XI(NIN)
C
      XL = XO(NA) - DLOP*WO(NA)
      XH = XO(NZ) + DLOP*WO(NZ)
      CALL REBFNC(IRFUN,XL,NPAR,NPACT,DPAR)
      CALL REBFNC(IRFUN,XH,NPAR,NPACT,DPAR)
C  complete mapping of XO
      IF (XL.GE.DL .AND. XH.LE.DH) GO TO 50
C  no mapping of XO
      IF (XL.GE.DH .OR. XH.LE.DL) THEN
          WRITE (*,*) ' No overlap between input',
     +      ' and output independent !','/'
      END IF
C  check lower DX bound
      IN = NA
      DO 10, I=NA,NZ,LOP
          XL = XO(I) - DLOP*WO(I)
          IN = I
          CALL REBFNC(IRFUN,XL,NPAR,NPACT,DPAR)
          IF (XL.GE.DL) GO TO 20
          YO(I) = 0.00
10    CONTINUE
20    NA = IN  
C  check upper DX bound
      IN     = NZ
      DO 30, I=NZ,NA,-LOP
          XH = XO(I) + DLOP*WO(I)
          IN = I
          CALL REBFNC(IRFUN,XH,NPAR,NPACT,DPAR)
          IF (XH.LT.DH) GO TO 40
          YO(I)  = 0.00
30    CONTINUE
40    NZ     = IN
C  the initial lower limit DX
50    IF (LIP.EQ.1) THEN
          LL = 1
          LH = NIN
      ELSE
          LL = NIN
          LH = 1
      ENDIF
C
C ... and finally we go through the XO array
C
      DO 130, I = NA,NZ,LOP
          VAL = 0.0
          XW  = DLOP*WO(I)
          XL  = XO(I) - XW
          XH  = XL + XW + XW
          CALL REBFNC(IRFUN,XL,NPAR,NPACT,DPAR)
          CALL REBFNC(IRFUN,XH,NPAR,NPACT,DPAR)
C
C   This is overdoing it for arrays with equal binwidth and pixel
C   filling factors of 1, (i.e. next run through loop XR could be
C   copied into the new XL), but we want to be fancy (general).
C
C   Now search for the XI(j's) falling between XL and XR
C   Recall that we set up everything in such a way that we are
C   increasing in the value of XI and that the first XL will fall
C   after XI(LL) and the last XR before XI(LH). We will update
C   LL each time (we can do so because XO and XI were assumed to be
C   monotonic (or at least constant, but that doesn't make sense
C   at all).
C
C   Fast loop here to find full pixels falling into XL,XR
C   A complete pix is defined GE left, LT right border
C
C   DL is the distance between the lower border of pix XI(j) and
C   XL - if positive or 0 then we found the first applic. pixel.
C
C   found the pix
          DO 60, J = LL,LH,LIP  
C   this is the low pix
              DL = XL - XI(J) - DLIP*WI(J)  
              IF (DL.LT.0.0D0) THEN
                  ML = J
C   all in one pix
                  DH = XI(J) + DLIP*WI(J) - XH  
C   goto one-pix handling
                  IF (DH.GE.0.0D0) THEN  
                      MH = ML
                      GO TO 110
                  ENDIF
                  GOTO 70
              ENDIF
60        CONTINUE
C
          CALL STETER(41,' ERROR - no lower bound of pixel')
C
C   DR is the distance between the left border of pix XI(j) and
C   XR - if positive or 0 then we found the last applic. pixel.
C   Note that ML+1 is always .LE.LH because of boundary search
C
C  dist. to right border
CMvK   70     DO 80, J = ML + 1,LH,LIP  
   70     DO 80, J = ML+LIP,LH,LIP  
              MH = J
              DH = XI(J) + DLIP*WI(J) - XH  
C  JMP to MH integration
              IF (DH.GT.0.0D0) THEN
                  GO TO 90
              ELSE  
C  JMP to ML integration
                  VAL = VAL + SNGL(YI(J))
C  instead of  :  IF (DH.EQ.0.0D0) GO TO 100
                  IF (DABS(DH) .LT. 1.D-30) GOTO 100
              ENDIF
80        CONTINUE
C
          CALL STETER(41,' ERROR - no upper bound of pixel')
C
C  handle the integrations (MH.NE.ML), in MH first
C
90        PL = XI(MH) - DLIP*WI(MH)
          XINCR  = (XH-PL)/XVLS
          CALL REBPIX(IOPT,NIN,XI,YI,WI,PL,XH,MINT,XINCR,MH,LIP,VALI)
          VAL = VAL + VALI/WI(MH)
C
C  and look into ML in any case
C
100       DL     = XL - XI(ML) + DLIP*WI(ML)
          IF (DL.EQ.0.0D0) THEN
              VAL    = VAL + SNGL(YI(ML))

          ELSE
              PH     = XI(ML) + DLIP*WI(ML)
              XINCR  = (PH-XL)/XVLS
              CALL REBPIX(IOPT,NIN,XI,YI,WI,XL,PH,MINT,XINCR,ML,LIP,
     +                    VALI)
              VAL    = VAL + VALI/WI(ML)
          END IF
C
          GO TO 120
C
C  handle the single pix case (MH.EQ.ML)
C
C  only one full pixel
110       DL     = XL - XI(ML) + DLIP*WI(ML)  
          IF (DL.EQ.0.0D0 .AND. DH.EQ.0.0D0) THEN
              VAL    = SNGL(YI(MH))

          ELSE
              XINCR  = (XH-XL)/XVLS
              CALL REBPIX(IOPT,NIN,XI,YI,WI,XL,XH,MINT,XINCR,ML,LIP,
     +                    VALI)
              VAL    = VALI/WI(ML)
          END IF
C
C  update cuts and result array, CONSERVE FLUX IN PIXEL DOMAIN
C
120       YO(I)  = VAL
          IF (I.EQ.1) THEN
             RMIN = YO(1)
             RMAX = RMIN
          ELSE
             RMIN   = AMIN1(RMIN,YO(I))
             RMAX   = AMAX1(RMAX,YO(I))
          ENDIF
C
C  update the loop start LL
C
CMvK       LL     = MH
          LL = ML
130   CONTINUE
C
C  Here is the end of REBMET
C
      RETURN
      END
C 
      SUBROUTINE REBPIX(JOP,NP,DX,DY,WI,XA,XB,NPOI,DINCR,ISTA,IDIR,F)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  Determine flux in fractions of pixel using various options JOP
C
C.INPUT
C JOP I*4 option
C NP      I*4 number of array elements DX,DY
C DX      R*8 independent array
C       DY      R*8     dependent array
C     WI R*8 bin (or pixel) widths
C XA      R*8     lower boundary of area of interest
C       XB      R*8     upper boundary
C       NPOI    I*4     number of points for integration/interpolation
C                       to be used
C       DINCR   R*8     interval between succesive points NPOI
C       ISTA    I*4     start index or search loop in DX
C       IDIR    I*4     Direction of DX array (X(1) < X(2) => -1)
C
C.OUTPUT
C       F       R*4     intensity between XA and XB under DY = f(DX)
C
C -------------------------------------------------------------------
C
      INTEGER JOP,NP,NPOI,ISTA,IDIR
      REAL F
      DOUBLE PRECISION DX(NP),DY(NP),WI(NP),XA,XB
      DOUBLE PRECISION DINCR
C
      GO TO (10,20,30),JOP
CD WRITE (6,*) ' The option requested for interp./integ. is not'
CD.                'possible'
CD WRITE (6,*) ' Simple proportion of pixel will be used'
      GO TO 20
C
C  Spline interpolation and gaussian integration
C
   10 CALL REBISP(NP,DX,DY,XA,XB,NPOI,DINCR,ISTA,IDIR,F)
      RETURN
C
C   Intensity is simply the corresponding proportion of the pixel's
C   intensity
   20 F      = DY(ISTA)* (XB-XA)
      RETURN
C
C  linear interpolation and integration
   30 CALL REBLNR(NP,DX,DY,XA,XB,ISTA,IDIR,F)
      RETURN
C
      END

      SUBROUTINE REBISP(NP,DX,DY,XA,XB,NPOI,DINCR,IST,IDI,F)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  F. Murtagh ST-ECF  Version 1.0 15-AUG-1985
C  M. Rosa      ST-ECF                  2.0     15-OCT-1985
C
C  Set up for REBIPL and integrate --> F
C
      INTEGER I,NP,NPOI,IST,IDI,IFAIL,IS,ID,FLAG
      REAL F
      DOUBLE PRECISION DX(NP),DY(NP),XA,XB
      DOUBLE PRECISION VECI(10),VECD(10)
      DOUBLE PRECISION DINCR,DVAL,DER,XP
C
      IFAIL  = 0
      IS     = IST
      ID     = IDI
C
      XP     = XA
      FLAG=1
      CALL REBIPL(FLAG,XP,DVAL,DX,DY,NP,IS,ID)
      VECI(1) = XP
      VECD(1) = DVAL
C
      DO 10 I = 2,NPOI
          XP     = XP + DINCR
          FLAG=0
          CALL REBIPL(FLAG,XP,DVAL,DX,DY,NP,IS,ID)
          VECI(I) = XP
          VECD(I) = DVAL
   10 CONTINUE
C
      DVAL = 0
      DER = 0
      CALL FINDIFF(VECI,VECD,NPOI,DVAL,DER,IFAIL)
      F      = SNGL(DVAL)
      RETURN

      END

      SUBROUTINE REBIPL(FLAG,XP,P,X,F,N,IS,ID)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Interpolate a function P at a given value XP using
C the table of values (X,F).
C Spline interpolation scheme based on Hermite polynomials.
C
C Ref.: U.S. Airforce Surveys in Geophysics, no.272
C       Publ of the Dominion Astrophys. Obs. 204, 1982
C
C.AUTHOR  J.D. PONZ (ESO)
C
C.MODIFIED  M. ROSA (STECF) 02-OCT-1985
C
C                 DOUBLE PRECISION calculation
C                 Stepping through X is always increasing X and XP
C                 All checking for validity range done outside
C                 Loop parameters IS (start) and ID (direction)
C           supplied from outside
C
C Arguments :
C XP      Input argument
C P Output interpolated value
C X Input vector of independent values
C F Input vector of dependent values
C N Input number of points in (X,F) vector
C       IS      Input   index of respective pixel
C       ID      Input   direction of X array
C
      INTEGER N,IS,ID,KA,KZ,IRUN,I,FLAG
      DOUBLE PRECISION F(N),X(N)
      DOUBLE PRECISION XP,P
      DOUBLE PRECISION LP1,LP2,L1,L2,FP1,FP2,XPI1,XPI,L1S,L2S
      SAVE KA,KZ,IRUN,FP1,FP2,LP1,LP2
C
      IF (FLAG.EQ.1) THEN
      IF (ID.GT.0) THEN
          KA     = MAX0(IS-ID,1)
          KZ     = N

      ELSE
          KA     = MIN0(IS-ID,N)
          KZ     = 1
      END IF

      IRUN   = 0
      ENDIF
C
C
      DO 10 I = KA,KZ,ID
          IF (XP.LT.X(I)) GO TO 20
   10 CONTINUE
      P      = 10.**36
      RETURN
C
   20 I      = I - ID
      IF (I.EQ.KA-1 .AND. IRUN.NE.0) GO TO 30
      KA     = I + ID
      IRUN   = 1
      LP1    = 1./ (X(I)-X(I+1))
      LP2    = 1./ (X(I+1)-X(I))
C
      IF (I.EQ.1) THEN
          FP1    = (F(2)-F(1))/ (X(2)-X(1))

      ELSE
          FP1    = (F(I+1)-F(I-1))/ (X(I+1)-X(I-1))
      END IF
C
      IF (I.EQ.N-1) THEN
          FP2    = (F(N)-F(N-1))/ (X(N)-X(N-1))

      ELSE
          FP2    = (F(I+2)-F(I))/ (X(I+2)-X(I))
      END IF
C
   30 XPI1   = XP - X(I+1)
      XPI    = XP - X(I)
      L1     = XPI1*LP1
      L2     = XPI*LP2
      L1S    = L1*L1
      L2S    = L2*L2
      P      = F(I)* (1.D0-2.D0*LP1*XPI)*L1S +
     +         F(I+1)* (1.D0-2.D0*LP2*XPI1)*L2S + FP2*XPI1*L2S +
     +         FP1*XPI*L1S
C
      RETURN

      END

      SUBROUTINE REBFNC(NF,D,NPA,NPAC,PA)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.AUTHOR M.ROSA (STECF)  VERSION 1.0     02-OCT-1985
C        F.MURTAGH (STECF)              1.1     11-NOV-1985
C              (more options, check on existence)
C
C  Dummy routine to call the requested function f[NF]
C  Returns D[new] as f(D[old])
C
C  Functions:  1 LIN    D --> a + b*D
C              2 POL    D --> a + b*D + c*DD + ...
C              3 INV    D --> a + b/D
C              4 EXP    D --> a + b*exp(D*c)
C              5 DEX    D --> a + b*10.**(D*c)
C              6 LOG    D --> a + b*LOG(D*c)
C              7 DLG    D --> a + b*LOG10(D*c)
C              8 IPO    D --> a + b/(D-q) + c/((D-q)**2) + ...
C               where q = PA(NPAC) and a,b,c... = PA(1,2,...(NPAC-1))
C              9 U01    User defined function REBU01
C
      INTEGER NF,NPA,NPAC
      DOUBLE PRECISION D,PA(NPA)
C
      GO TO (10,20,30,40,50,60,70,80,90) NF
C
      WRITE (*,*) ' Function not existent, sorry - look up manual '
      D      = 0.
      RETURN
C
   10 CALL REBLIN(D,PA,NPA,NPAC)
      RETURN

   20 CALL REBPOL(D,PA,NPA,NPAC)
      RETURN

   30 CALL REBINV(D,PA,NPA,NPAC)
      RETURN

   40 CALL REBEXP(D,PA,NPA,NPAC)
      RETURN

   50 CALL REBDEX(D,PA,NPA,NPAC)
      RETURN

   60 CALL REBLOG(D,PA,NPA,NPAC)
      RETURN

   70 CALL REBDLG(D,PA,NPA,NPAC)
      RETURN

   80 CALL REBIPO(D,PA,NPA,NPAC)
      RETURN

   90 CALL REBU01(D,PA,NPA,NPAC)
      RETURN
C
      END

      SUBROUTINE REBLIN(DV,P,NP,NPC)
      INTEGER NP,NPC
      DOUBLE PRECISION DV,P(NP)

      DV     = P(1) + P(2)*DV
      RETURN

      END
C
      SUBROUTINE REBPOL(DV,P,NP,NPC)
      INTEGER NP,NPC,I
      DOUBLE PRECISION DV,P(NP),VD  !MP 900117
      VD     = 0.0D0
      DO 10 I = NPC,1,-1
          VD     = VD*DV + P(I)
   10 CONTINUE
      DV     = VD
      RETURN

      END
C
      SUBROUTINE REBINV(DV,P,NP,NPC)
      INTEGER NP,NPC
      DOUBLE PRECISION DV,P(NP)

      DV     = P(1) + P(2)/DV
      RETURN

      END
C
      SUBROUTINE REBEXP(DV,P,NP,NPC)
      INTEGER NP,NPC
      DOUBLE PRECISION DV,P(NP)

      DV     = P(1) + P(2)*EXP(DV*P(3))
      RETURN

      END
C
      SUBROUTINE REBDEX(DV,P,NP,NPC)
      INTEGER NP,NPC
      DOUBLE PRECISION DV,P(NP)

      DV     = P(1) + P(2)*10.0D0** (DV*P(3))
      RETURN

      END
C
      SUBROUTINE REBLOG(DV,P,NP,NPC)
      INTEGER NP,NPC
      DOUBLE PRECISION DV,P(NP)

      DV     = P(1) + P(2)*LOG(DV*P(3))
      RETURN

      END
C
      SUBROUTINE REBDLG(DV,P,NP,NPC)
      INTEGER NP,NPC
      DOUBLE PRECISION DV,P(NP)

      DV     = P(1) + P(2)*LOG10(DV*P(3))
      RETURN

      END
C
      SUBROUTINE REBIPO(DV,P,NP,NPC)
C
C INVERSE POLYNOMIAL WITH OFFSET
C Implement the PRISM formula for dispersion relations
C
C y = a + bz + cz**2 + dz**3 + .... P(npc-1)z**(npc-1)
C
C where z = x - p(npc)
C
      INTEGER NP,NPC,I
      DOUBLE PRECISION P(NP),DV,VD

      VD     = 0.0D0
      DV     = 1./ (DV-P(NPC))
      IF (NPC.EQ.1) RETURN
      DO 10 I = NPC - 1,1,-1
          VD     = VD*DV + P(I)
   10 CONTINUE
      DV     = VD
      RETURN

      END
C

      SUBROUTINE REBLNR(NP,DX,DY,XA,XB,ISTA,IDIR,F)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C   Subroutine linear
C
C   F. Murtagh   ST-ECF  Version 1.0 17-Oct-1985
C
C   linear interpolation and subsequent integration.
C
C.INPUT
C  NP I*4 Number of array elements DX, DY
C  DX R*8 Independent array
C DY R*8 Dependent array
C XA R*8 Lower boundary of area of interest
C XB R*8 Upper boundary
C ISTA I*4 Start index for search loop in DX;
C   index number of pixel in which both
C   XA and XB are located
C IDIR I*4 Direction of DX array (X(1) < X(2) => -1)
C
C.OUTPUT
C F R*4 Intensity between XA and XB under DY = f(DX)
C
C------------------------------------------------------------------
      INTEGER NP,ISTA,IDIR
      REAL F,YA,YB,XLIN
      DOUBLE PRECISION DX(NP),DY(NP),XA,XB
      DOUBLE PRECISION X,Y,X2,Y2
C
C For convenience:
      X      = DX(ISTA)
      Y      = DY(ISTA)
C
C Four cases arise: XA or XB = X; otherwise, XA and XB are on the same
C side of the pixel-centre, or they are on different sides.
      IF (XA.EQ.X) THEN
          IF (XB-X.GE.0) THEN
              Y2     = DY(ISTA+IDIR)
              X2     = DX(ISTA+IDIR)

          ELSE
              Y2     = DY(ISTA-IDIR)
              X2     = DX(ISTA-IDIR)
          END IF

          YB     = XLIN(Y2,Y,X,XB,X2)
          F      = ABS(X-XB)* (Y+YB)/2.

      ELSE IF (XB.EQ.X) THEN
          IF (XA-X.GE.0) THEN
              Y2     = DY(ISTA+IDIR)
              X2     = DX(ISTA+IDIR)

          ELSE
              Y2     = DY(ISTA-IDIR)
              X2     = DX(ISTA-IDIR)
          END IF

          YA     = XLIN(Y2,Y,X,XA,X2)
          F      = ABS(X-XA)* (Y+YA)/2.

      ELSE IF (((XA-X)/ (XB-X)).GT.0) THEN
C       (If XA and XB are both on the same side of the pixel-centre,
C       we get num. and denom. both of same sign, hence above is true:
          IF (XA-X.GT.0) THEN
              Y2     = DY(ISTA+IDIR)
              X2     = DX(ISTA+IDIR)

          ELSE
              Y2     = DY(ISTA-IDIR)
              X2     = DX(ISTA-IDIR)
          END IF

          YA     = XLIN(Y2,Y,X,XA,X2)
          YB     = XLIN(Y2,Y,X,XB,X2)
          F      = (XB-XA)* (YA+YB)/2.

      ELSE IF (((XA-X)/ (XB-X)).LE.0) THEN
          IF (XA-X.GT.0) THEN
              YA     = XLIN(DY(ISTA+IDIR),Y,X,XA,DX(ISTA+IDIR))
              YB     = XLIN(DY(ISTA-IDIR),Y,X,XB,DX(ISTA-IDIR))

          ELSE
              YA     = XLIN(DY(ISTA-IDIR),Y,X,XA,DX(ISTA-IDIR))
              YB     = XLIN(DY(ISTA+IDIR),Y,X,XB,DX(ISTA+IDIR))
          END IF

          F      = ABS(X-XA)* (Y+YA)/2. + ABS(X-XB)* (Y+YB)/2.
      END IF
C
      RETURN

      END

      REAL FUNCTION XLIN(Y2,Y1,X1,X,X2)
C
C Function which returns linearly interpolated value corresponding to
C X, where we have X1 <= X <= X2 or vice versa, and Y1 = f(X1) and
C Y2 = f(X2).
C
      DOUBLE PRECISION Y1,Y2,X,X1,X2
C
      XLIN   = (Y2-Y1)* (X-X1)/ (X2-X1) + Y1
C
      END

      SUBROUTINE FINDIFF(X,Y,N,INT,ERR,IFAIL)
C       This subroutine integrates a function which is specified numerically
C       at four or more points, using third-order finite-difference formulae 
C       according to a method due to Gill and Miller. (Computer Journal 15,
C       page 80-83,1972).
C       Author: Michele Peron, ESO-IPG. August 1990



        INTEGER N,I,J,K,IFAIL
        DOUBLE PRECISION X(N),Y(N),INT,ERR,E,H1,H2,H3,H4,R1,R2,R3,R4
        DOUBLE PRECISION D1,D2,D3,C,S,H3S
        IFAIL = 0
        J  = 3
        K  =  N-1
        ERR = 0
        E = 0
        INT=0
        H1 = 0
        H2 = 0
        H3 = 0
        H4 = 0
        R1 = 0
        R2 =0
        R3 = 0
        R4 = 0
        D1  = 0
        D2 = 0
        D3 = 0
        C = 0
        S = 0
        DO 10 I=J,K
         IF (I.EQ.J) THEN        
             H2 = X(J-1)-X(J-2)
             D3 = (Y(J-1)-Y(J-2))/H2
             H3 = X(J)-X(J-1)
             D1  = (Y(J)-Y(J-1))/H3
             H1 = H2+H3
             D2 = (D1-D3)/H1
             H4 = X(J+1)-X(J)
             R1 = (Y(J+1)-Y(J))/H4
             R2 = (R1-D1)/(H4+H3)
             H1 = H1+H4
             R3 = (R2-D2)/H1
             INT = H2*(Y(1)+H2*(D3/2.0-H2*
     1               (D2/6.-(H2+2.*H3)*R3/12.)))
             S  = -H2**3*(H2*(3.*H2+5.*H4)+10.*H3*H1)/60.
         ELSE
             H4 = X(I+1)-X(I)
             R1 = (Y(I+1)-Y(I))/H4
             R4 = H4+H3
             R2 = (R1-D1)/R4
             R4 = R4+H2
             R3 = (R2-D2)/R4
             R4 = R4+H1
             R4 = (R3-D3)/R4
         ENDIF
         H3S = H3**2
        INT = INT+H3*((Y(I)+Y(I-1))/2.-H3S*
     1         (D2+R2+(H2-H4)*R3)/12.)
        C = H3*H3S*(2.*H3S+5.*(H3*(H4+H2)+
     2             2.*H4*H2))/120.
        E = E+(C+S)*R4
        IF (I.EQ.J) THEN 
           S = 2*C+S
        ELSE
           S = C
        ENDIF
        IF (I.EQ.K) THEN
            INT = INT+H4*(Y(N)-H4*(R1/2.+H4*(R2/6.+
     1                (2.*H3+H4)*R3/12.)))
            E = E-H4**3*R4*(H4*(3.*H4+5.*H2)+10.
     2            *H3*(H2+H3+H4))/60. 
            E = E+S*R4
        ELSE
            H1 = H2
            H2 = H3
            H3 = H4
            D1 = R1
            D2 = R2
            D3 = R3
        ENDIF
10      CONTINUE
        ERR = E
        INT = INT+E
        END
           
