C @(#)minmax.for	17.1.1.1 (ESO-DMD) 01/25/02 17:40:35
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
      SUBROUTINE MNMX(ARRAY,SIZE,CUTS,PIX)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine MNMX            version 1.50      850813
C             MYMX            version 1.50      850813
C             MNVAL            version 1.1      830407
C             MOMENT                    1.30	840116
C  K. Banse                        ESO - Garching
C  source code on minmax.for
C
C.KEYWORDS
C  minimum, maximum, mean value, standard deviation
C
C.PURPOSE
C  calculate min, max and mean + standard deviation of an array
C
C.ALGORITHM
C  std is calculated as SQRT of VAR(X), with VAR(X) = E(X)**2 - E(X**2)
C
C.INPUT/OUTPUT
C  call as MNMX(ARRAY,SIZE,CUTS,PIX)
C          DMYMX(ARRAY,SIZE,CUTS)	double precision
C          MYMX(ARRAY,SIZE,CUTS)
C          MNVAL(ARRAY,NAXIS,NPIX,PIXELS,FMEAN,STD,FMIN,FMAX)
C          MOMENT(ARRAY,NAXIS,NPIX,PIXELS,FMOM,STD,FMIN,FMAX,MINPIX,MAXPIX)
C
C  input par:
C  ARRAY:	R*4 array      input data array
C  SIZE:	I*4            size of ARRAY
C  NAXIS:	I*4            no. of dimensions of ARRAY
C  NPIX:	I*4 array      no. of pixels/rows in complete frame
C  PIXELS:	I*4 array      start + end pixels in each dimension
C
C  output par:
C  CUTS:	R*4 array      minimum,maximum of ARRAY
C  FMEAN:	R*4            mean value of ARRAY
C  STD:		R*4            standard deviation from FMEAN
C  FMOM:	R*4 array      1., 2., 3. + 4. moments of ARRAY
C  PIX:		I*4 array      index of min. + max. pixel
C
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER      SIZE,N
      INTEGER      MINPIX,MAXPIX,PIX(2)
C
      REAL         ARRAY(SIZE)
      REAL         CUTS(2),FMIN,FMAX
C
C  initialize
      FMIN = ARRAY(1)
      MINPIX = 1
      FMAX = FMIN
      MAXPIX = MINPIX
C
C  now loop through array
      DO 100, N=1,SIZE
         IF (ARRAY(N).LT.FMIN) THEN
            FMIN = ARRAY(N)
            MINPIX = N
         ELSE IF (ARRAY(N).GT.FMAX) THEN
            FMAX = ARRAY(N)
            MAXPIX = N
         ENDIF
100   CONTINUE
C
C  set up output parms
      CUTS(1) = FMIN
      CUTS(2) = FMAX
      PIX(1) = MINPIX
      PIX(2) = MAXPIX
C
      RETURN
      END

      SUBROUTINE DMYMX(ARRAY,SIZE,CUTS)
C
      IMPLICIT NONE
C
      INTEGER      SIZE,N
C
      DOUBLE PRECISION   ARRAY(SIZE)
      DOUBLE PRECISION   CUTS(2),FMIN,FMAX
C
C  initialize
      FMIN = ARRAY(1)
      FMAX = FMIN
C
C  now loop through array
      DO 100, N=1,SIZE
         IF (ARRAY(N).LT.FMIN) THEN
            FMIN = ARRAY(N)
         ELSE IF (ARRAY(N).GT.FMAX) THEN
            FMAX = ARRAY(N)
         ENDIF
100   CONTINUE
C
C  set up output parms
      CUTS(1) = FMIN
      CUTS(2) = FMAX
C
      RETURN
      END

      SUBROUTINE MYMX(ARRAY,SIZE,CUTS)
C
      IMPLICIT NONE
C
      INTEGER      SIZE,N
C
      REAL         ARRAY(SIZE)
      REAL         CUTS(2),FMIN,FMAX
C
C  initialize
      FMIN = ARRAY(1)
      FMAX = FMIN
C
C  now loop through array
      DO 100, N=1,SIZE
         IF (ARRAY(N).LT.FMIN) THEN
            FMIN = ARRAY(N)
         ELSE IF (ARRAY(N).GT.FMAX) THEN
            FMAX = ARRAY(N)
         ENDIF
100   CONTINUE
C
C  set up output parms
      CUTS(1) = FMIN
      CUTS(2) = FMAX
C
      RETURN
      END

      SUBROUTINE MNVAL(ARRAY,NAXIS,NPIX,PIXELS,FMEAN,STD,FMIN,FMAX)
C
      IMPLICIT NONE
C 
      INTEGER      NAXIS,NPIX(*),PIXELS(2,2)
      INTEGER      NOPIX,LOWX,HIX,SIZE
      INTEGER      LOWY,HIY,IOFF,NX,NY,N
C 
      REAL     ARRAY(*),FMIN,FMAX,FMEAN,STD,PIX
C 
      DOUBLE PRECISION    SUM,SUMS,S1,S2,EXPECT,SIGMA
C
C  initialize
      NOPIX = NPIX(1)                     !points per line in father frame...
      SUM = 0.D0
      SUMS = 0.D0
      LOWX = PIXELS(1,1)
      HIX = PIXELS(1,2)
      SIZE = HIX - LOWX + 1
C
      IF (NAXIS.GE.2) THEN
         LOWY = PIXELS(2,1)
         HIY = PIXELS(2,2)
         SIZE = SIZE * (HIY - LOWY + 1)
      ELSE
         LOWY = 1
         HIY = 1
      ENDIF
C
      N = (LOWY-1)*NOPIX + LOWX           !get index of 1. element in sybframe
      FMIN = ARRAY(N)
      FMAX = FMIN
C
C  now loop through array
      DO 800, NY=LOWY,HIY
         IOFF = (NY-1)*NOPIX
         DO 500, NX=LOWX,HIX
            N = IOFF + NX
            PIX = ARRAY(N)
            IF(PIX.LT.FMIN) THEN
               FMIN = PIX
            ELSE
               IF(PIX.GT.FMAX) FMAX = PIX
            ENDIF
            SUM = SUM + PIX                     !sum of values
            SUMS = SUMS + (PIX*PIX)            	!sum of squares
500      CONTINUE
800   CONTINUE
C
C  now calculate mean + standard deviation
      EXPECT = SUM/SIZE
      S1 = SUMS/SIZE
      S2 = EXPECT * EXPECT
      SIGMA = SQRT(MAX(0.D0,(S1-S2)))
      FMEAN = EXPECT                         !get single precision values
      STD = SIGMA
C 
      RETURN
      END

      SUBROUTINE MOMENT(ARRAY,NAXIS,NPIX,PIXELS,FMOM,STD,FMIN,FMAX,
     +                     MINPIX,MAXPIX)
C
      IMPLICIT NONE
C 
      INTEGER      NAXIS,NPIX(*),PIXELS(2,2)
      INTEGER      NOPIX,LOWX,HIX,SIZE
      INTEGER      LOWY,HIY,IOFF,NX,NY
      INTEGER      MINPIX,MAXPIX,N
C 
      REAL            ARRAY(*),FMIN,FMAX,FMOM(4),STD
      REAL            PIX,PIX2,PIX3,PIX4
C 
      DOUBLE PRECISION SUM,SUM2,SUM3,SUM4,SIGMA
C
      DATA      SUM /0.D0/, SUM2 /0.D0/, SUM3 /0.D0/, SUM4 /0.D0/
C
C  initialize
      NOPIX = NPIX(1)                   !points per line in father frame...
      LOWX = PIXELS(1,1)
      HIX = PIXELS(1,2)
      SIZE = HIX - LOWX + 1
C
      IF (NAXIS.GE.2) THEN
         LOWY = PIXELS(2,1)
         HIY = PIXELS(2,2)
         SIZE = SIZE * (HIY - LOWY + 1)
      ELSE
         LOWY = 1
         HIY = 1
      ENDIF
C
      N = (LOWY-1)*NOPIX + LOWX          !get index of 1. element in subframe
      FMIN = ARRAY(N)
      MINPIX = N
      FMAX = FMIN
      MAXPIX = MINPIX
C
C  now loop through array
      IOFF = (LOWY-1) * NOPIX            		!init offset...
      DO 800, NY=LOWY,HIY
C
         DO 600, NX=LOWX,HIX
            N = IOFF + NX
            PIX = ARRAY(N)
            IF (PIX.LT.FMIN) THEN
               FMIN = PIX
             MINPIX = N
            ELSE
               IF (PIX.GT.FMAX) THEN
                FMAX = PIX
                MAXPIX = N
             ENDIF
            ENDIF
            PIX2 = PIX*PIX
            PIX3 = PIX2*PIX
            PIX4 = PIX3*PIX
            SUM = SUM + PIX                     !sum of values
            SUM2 = SUM2 + PIX2                  !sum of squares
            SUM3 = SUM3 + PIX3                  !sum of PIX**3
            SUM4 = SUM4 + PIX4                  !sum of PIX**4
600      CONTINUE
C
         IOFF = IOFF + NOPIX                    !increment offset...
800   CONTINUE
C
C  now calculate moments + standard deviation
      SUM = SUM/SIZE                            !1. moment
      SUM2 = SUM2/SIZE                          !2. moment
      SUM3 = SUM3/SIZE                          !3. moment
      SUM4 = SUM4/SIZE                          !4. moment
      SIGMA = SQRT(MAX(0.D0,(SUM2-(SUM*SUM))))  !standard deviation
      FMOM(1) = SUM                             !get single precision values
      FMOM(2) = SUM2
      FMOM(3) = SUM3
      FMOM(4) = SUM4
      STD = SIGMA
C
C  return
      RETURN
      END