C  @(#)avwndw.for	17.1.1.1 (ESO-DMD) 01/25/02 17:40:56 
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 AVWNDW
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  program AVWNDW			version 1.00	890123
C  K. Banse                  		ESO - Garching
C  1.20	      890529	1.30	900130
C
C.KEYWORDS
C  bulk data frames, average
C
C.PURPOSE
C  average a no. of image frames, take only values which don't differ too much
C  from the median value (at each pixel), the result will either be a frame with
C  size equal to common area of all frames involved or size equal to union of
C  all frame areas
C
C.ALGORITHM
C  extract input frames either from given parameter string or from catalogue,
C  add up all frames + divide in the end
C  max. 20 frames may be averaged at one time !
C
C.INPUT/OUTPUT
C  the following keys are used:
C
C  ACTION/C/1/3			= WIN, MED, MIN or MAX
C  OUT_A/C/1/60			result frame
C  P3/C/1/80			list of frames to be averaged
C  INPUTR/R/1/2			BGERR,SNOISE - Poisson + Gaussian errors
C				for WINDOW option
C
C  the following descriptors are used with WINDOW option:
C  LHCUTS/R/5/2			low, high cuts for calculations
C  FLAT_BKG/R/1/1		background value
C  O_TIME/D/7/1			exposure time
C	
C.VERSIONS
C 
C 010802		last modif
C	
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER     BEGIN,FRMCNT,IAV,ISEQ
      INTEGER     L,LL,LP,M,IOFF
      INTEGER     N,NAXIS,NAXISC,NCHECK,NIN,NP,NPERC,NTL
      INTEGER*8   PNTRC,PNTRT
      INTEGER     SIZEC,STAT
      INTEGER     IMNOC,IMN(20)
      INTEGER     NPIX(3),NPIXC(3)
      INTEGER     OFX(20),OFY(20),POFF(20),NPNTS(20)
      INTEGER     IORD(20),NORD(20),REJ(20),TNP(4)
      INTEGER     EC,EL,ED
      INTEGER     UNI(1),NULO,MADRID(1)
      INTEGER     LEN,INDEX            !system functions
C
      CHARACTER   FRAME(20)*60,FRAMEC*60,CATFIL*60
      CHARACTER   CUNIT*64,IDENT*72,OUTPUT*80
      CHARACTER   NOCUTS*60,LINE*80,ACTION*3
      CHARACTER   ERROR1*60,ERROR2*50,ERROR3*40,ERROR4*40,ERROR6*40
      CHARACTER   ERROR7*40,ERROR9*60
C
      DOUBLE PRECISION STEP(2,20),STEPC(2)
      DOUBLE PRECISION START(2,20),STARTC(2),OSTART(2),OEND(2)
C
      REAL        EPS(3),CUTS(4)
      REAL        VAL(20),BGV(20),RMS(20),RFC(20)
      REAL        AJCUL(20),AJCUH(20),AJCU(2),SNOISE,BGERR
      REAL        INPUTR(2),EXPTI(20)
C
      DOUBLE PRECISION      EXPTIM
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
C
      COMMON   /VMR/     MADRID
      COMMON   /MYREAL/  VAL,BGV,RMS,RFC,AJCUL,AJCUH,SNOISE,BGERR,EXPTI
      COMMON   /MYINT/   IORD,NORD,REJ,TNP
C
      DATA    IDENT    /'average frame '/
      DATA    ERROR1   /'operands do not match in stepsize... '/
      DATA    ERROR2   /'operands do not overlap...'/
      DATA    ERROR3   /'stepsizes of different sign...'/
      DATA    ERROR4   /'catalogue empty...'/
      DATA    ERROR6   /'no. of input frames < 2 ...'/
      DATA    ERROR7   /'more than 20 entries in catalogue...'/
      DATA    ERROR9   /'no. of axes do not match... '/
      DATA    NOCUTS   /'LHCUTS(5,6) not there - use LHCUTS(3,4)... '/
      DATA    STARTC   /2*0.D0/,    STEPC  /2*1.D0/
C
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
C  set up MIDAS environment + enable automatic error abort
      CALL STSPRO('AVWNDW')
C 
C  init COMMON variables
      DO 50 N=1,20
         REJ(N) = 0
50    CONTINUE
      DO 60 N=1,4
         TNP(N) = 0
60    CONTINUE
C
C  get necessary keys
      CALL STKRDC('ACTION',1,1,3,IAV,ACTION,UNI,NULO,STAT)    !get option
      CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEC,UNI,NULO,STAT)    !get result frame
      CALL STKRDC('P3',1,1,80,IAV,LINE,UNI,NULO,STAT)         !list of frames
      CALL UPCAS(ACTION,ACTION)
      IF (ACTION(1:1).EQ.'W') THEN                        !only for real window option
         CALL STKRDR('INPUTR',1,2,IAV,INPUTR,UNI,NULO,STAT)
         BGERR = INPUTR(1) * INPUTR(1)                   !take square of errors
         SNOISE = INPUTR(2) * INPUTR(2)
      ENDIF
C
C  test, if we have list of frames or catalogue
      LL = INDEX(LINE,'.cat')
      IF (LL.LE.0) THEN
         LL = INDEX(LINE,'.CAT')
         IF (LL.LE.0) GOTO 500                     !No. Treat list input
      ENDIF
C
C  here we handle input from a catalog
C  get name of first image file in catalog
C
      ISEQ = 0
      CATFIL(1:) = LINE(1:)
      CALL STCGET(CATFIL,0,FRAME(1),OUTPUT,ISEQ,STAT)
      IF (FRAME(1)(1:1).EQ.' ')                 !there must be at least 1 image
     +   CALL STETER(4,ERROR4)
C
C  now loop through catalogue
      DO 200 N=2,20                                !max. 20 frames
         CALL STCGET(CATFIL,0,FRAME(N),OUTPUT,ISEQ,STAT)
         IF (FRAME(N)(1:1).EQ.' ') THEN            !indicates the end ...
            FRMCNT = N - 1
            GOTO 1000
         ENDIF
200   CONTINUE
C
      CALL STCGET(CATFIL,0,LINE,OUTPUT,ISEQ,STAT)
      IF (LINE(1:1).NE.' ')
     +   CALL STTPUT
     +   ('Max. 20 frames are used, following frames ignored...',STAT)
      FRMCNT = 20
      GOTO 1000
C
C  here we handle input from single inout line
C  pull out operands      (max. 20)
C
500   BEGIN = 1
      DO 800 N=1,21                             !max. of 20 frames...!
         CALL EXTRSS(LINE,',',BEGIN,FRAME(N),LL)
         IF (LL.LE.0) THEN
            FRMCNT = N - 1
            GOTO 1000
         ELSE
            CALL CLNFRA(FRAME(N),FRAME(N),0)
         ENDIF
800   CONTINUE
C
C  check no. of input frames
1000  IF (FRMCNT.LT.2) CALL STETER(6,ERROR6)
C
      DO 1200 N=1,FRMCNT
         DO 1100 M=1,2                        !init values...
            START(M,N) = 0.D0
            STEP(M,N) = 1.D0
            NPIX(M) = 1
1100     CONTINUE
1200  CONTINUE
C
C  get initial area from 1. frame
      DO 1400 M=1,2                  !init values...
         OSTART(M) = 0.D0
         OEND(M) = 0.D0
         NPIX(M) = 1
1400  CONTINUE
C
      CALL STFOPN(FRAME(1),D_R4_FORMAT,0,F_IMA_TYPE,IMN(1),STAT)
      CALL STDRDI(IMN(1),'NAXIS',1,1,IAV,NAXISC,UNI,NULO,STAT)
      IF (NAXISC.GT.2) THEN                !we only work on max 2-dim. frames
         NAXISC = 2
         OUTPUT(1:) = 'currently only 1 or 2-dim frames supported... '
         CALL STTPUT(OUTPUT,STAT)
      ENDIF
C
      CALL STDRDI(IMN(1),'NPIX',1,2,IAV,NPIX,UNI,NULO,STAT)
      NPNTS(1) = NPIX(1)
      CALL STDRDD(IMN(1),'START',1,2,IAV,START(1,1),UNI,NULO,STAT)
      OSTART(1) = START(1,1)
      OSTART(2) = START(2,1)
      CALL STDRDD(IMN(1),'STEP',1,2,IAV,STEP(1,1),UNI,NULO,STAT)
      STEPC(1) = STEP(1,1)
      STEPC(2) = STEP(2,1)
      CALL STDRDC(IMN(1),'CUNIT',1,1,64,IAV,CUNIT,UNI,NULO,STAT)
C
C  for "real" window option get exposure-time, low + high cuts and background
      IF (ACTION(1:1).EQ.'W') THEN
         CALL STECNT('GET',EC,EL,ED)            !disable error abort
         CALL STECNT('PUT',1,0,0)
C
         CALL STDRDD(IMN(1),'O_TIME',7,1,IAV,EXPTIM,UNI,NULO,STAT)
         IF (STAT.NE.0) THEN
            EXPTI(1) = 1.0                      !default to 1 second
         ELSE
            EXPTI(1) = EXPTIM
         ENDIF
C
C  get LHCUTS(5,6) - if not there use LHCUTS(3,4)
         CALL STDRDR(IMN(1),'LHCUTS',5,2,IAV,AJCU,UNI,NULO,STAT)
         IF (STAT.NE.0) THEN
            CALL STTPUT(NOCUTS,STAT)
            CALL STDRDR(IMN(1),'LHCUTS',3,2,IAV,AJCU,UNI,NULO,STAT)
         ENDIF
         IF (AJCU(2).LE.AJCU(1))
     +      CALL STETER(43,'invalid LHCUTS: LHCUTS(6).LE.LHCUTS(5)...')
         AJCUL(1) = AJCU(1)
         AJCUH(1) = AJCU(2)
C
         CALL STDRDR(IMN(1),'FLAT_BKG',1,1,IAV,BGV(1),UNI,NULO,STAT)
         IF (STAT.NE.0) BGV(1) = 0.             !default to background = 0.
         CALL STECNT('PUT',EC,EL,ED)
      ENDIF
C
      DO 2000 M=1,2
         EPS(M) = 0.0001 * ABS(STEPC(M))        !0.01% of stepsize as epsilon
         OEND(M) = OSTART(M) + (NPIX(M)-1)*STEPC(M)
2000  CONTINUE
C
C  now loop through the other input frames
      DO 3000 N=2,FRMCNT
         CALL STFOPN(FRAME(N),D_R4_FORMAT,0,F_IMA_TYPE,IMN(N),STAT)
C
         CALL STDRDI(IMN(N),'NAXIS',1,1,IAV,NAXIS,UNI,NULO,STAT)
         IF (NAXISC.NE.NAXIS) CALL STETER(1,ERROR9)
C
         CALL STDRDI(IMN(N),'NPIX',1,2,IAV,NPIX,UNI,NULO,STAT)
         NPNTS(N) = NPIX(1)
         CALL STDRDD(IMN(N),'START',1,2,IAV,START(1,N),UNI,NULO,STAT)
         CALL STDRDD(IMN(N),'STEP',1,2,IAV,STEP(1,N),UNI,NULO,STAT)
C
         IF (ACTION(1:1).EQ.'W') THEN
            CALL STECNT('GET',EC,EL,ED)            !disable error abort
            CALL STECNT('PUT',1,0,0)
C
            CALL STDRDD(IMN(N),'O_TIME',7,1,IAV,EXPTIM,UNI,NULO,STAT)
            IF (STAT.NE.0) THEN
               EXPTI(N) = 1.0                      !default to 1 second
            ELSE
               EXPTI(N) = EXPTIM
            ENDIF
C
            CALL STDRDR(IMN(N),'LHCUTS',5,2,IAV,AJCU,UNI,NULO,STAT)
            IF (STAT.NE.0)
     +         CALL STDRDR(IMN(N),'LHCUTS',3,2,IAV,AJCU,UNI,NULO,STAT)
            AJCUL(N) = AJCU(1)
            AJCUH(N) = AJCU(2)
C
            CALL STDRDR(IMN(N),'FLAT_BKG',1,1,IAV,BGV(N),UNI,NULO,STAT)
            IF (STAT.NE.0) BGV(N) = 0.            !default to background = 0.
            CALL STECNT('PUT',EC,EL,ED)
         ENDIF
C
C  create new result frame with dimension as intersection of input frames
C  stepsizes should have same sign and not differ too much...
         DO 2400 M=1,NAXISC
            IF (STEPC(M)*STEP(M,N).LE.0.) CALL STETER(1,ERROR3)
            IF (ABS(STEP(M,N)-STEPC(M)).GT.EPS(M))
     +         CALL STTPUT(ERROR1,STAT)                 !just warning message
2400     CONTINUE
C
C  get intersection of image areas
         DO 2600 M=1,NAXISC
            IF (STEPC(M).GT.0.) THEN
               OSTART(M) = MAX(OSTART(M),START(M,N))
               OEND(M) = MIN(OEND(M),START(M,N) + (NPIX(M)-1)*STEP(M,N))
            ELSE
               OSTART(M) = MIN(OSTART(M),START(M,N))
               OEND(M) = MAX(OEND(M),START(M,N) + (NPIX(M)-1)*STEP(M,N))
            ENDIF
2600     CONTINUE
C
C  compute individual offsets
         DO 2800 NIN=1,N
            OFX(NIN) = NINT((OSTART(1)-START(1,NIN))/STEP(1,NIN)) + 1
            OFY(NIN) = NINT((OSTART(2)-START(2,NIN))/STEP(2,NIN)) + 1
2800     CONTINUE
C
3000  CONTINUE
C
C  test, if something is left...
C
      DO 3300 M=1,NAXISC
         IF (STEPC(M)*(OEND(M)-OSTART(M)).LT.0.) CALL STETER(2,ERROR2)
3300  CONTINUE
C
C  set up standard stuff for result frame
      SIZEC = 1
      NPIXC(2) = 1
      DO 3600 M=1,NAXISC
         STARTC(M) = OSTART(M)
         NPIXC(M) = NINT((OEND(M)-OSTART(M))/STEPC(M)) + 1
         SIZEC = SIZEC * NPIXC(M)
3600  CONTINUE
      LP = NPIXC(2)
      NP = NPIXC(1)
C
C  map result frame
      CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +            NAXISC,NPIXC,STARTC,STEPC,
     +            IDENT,CUNIT,PNTRC,IMNOC,STAT)
C
C  initialize some working variables
      DO 3900 N=1,FRMCNT
         IF (ACTION(1:1).EQ.'W') RFC(N) = 1./EXPTI(N)
         POFF(N) = (OFY(N)-1)*NPNTS(N) + OFX(N)      !determine 1.pixel no.
3900  CONTINUE
C
C  and allocate working space
      NTL = NP * FRMCNT
      IF (NTL.LT.512) NTL = 512
      CALL STFXMP(NTL,D_R4_FORMAT,PNTRT,STAT)
C
C  set up counters for calculation of progress
      NTL = INT(FLOAT(LP)*.2)                  !20% of total no. of lines
      IF (NTL.LE.2) THEN
         NCHECK = LP + 1                       !not worth it...
      ELSE
         NPERC = 20
         NCHECK = NTL
      ENDIF
C
C  now, loop through lines
      DO 5000 L=1,LP
C
C  and loop through frames
         IOFF = 1
         DO 4200 N=1,FRMCNT
            CALL COPULA(IOFF,IMN(N),POFF(N),MADRID(PNTRT),NP)
            POFF(N) = POFF(N) + NPNTS(N)                  !move to next line
            IOFF = IOFF + NP
4200     CONTINUE
C
C  and do it ...
         IF (ACTION(1:1).EQ.'W') THEN
            CALL AVWIN1(MADRID(PNTRT),FRMCNT,NP)
         ELSE
            CALL AVWIN2(ACTION,MADRID(PNTRT),FRMCNT,NP)
         ENDIF
         CALL COPUL(MADRID(PNTRT),MADRID(PNTRC),L,NP)    !save resulting line
C
C  display progress
         IF (L.GE.NCHECK) CALL PROGRS(20,NTL,NPERC,NCHECK)
5000  CONTINUE
C
C  display more info
      IF (ACTION(1:1).EQ.'W') THEN
         DO 5500 N=1,FRMCNT
            WRITE(OUTPUT,10000) AJCUL(N),AJCUH(N)
            L = INDEX(FRAME(N),' ') - 1
            IF (L.LE.0) L = LEN(FRAME(N))
            LINE(1:) = FRAME(N)(1:L)//OUTPUT(1:)
            CALL STTPUT(LINE,STAT)
            WRITE(OUTPUT,10001) EXPTI(N),BGV(N),REJ(N)
            CALL STTPUT(OUTPUT,STAT)
5500     CONTINUE
         WRITE(OUTPUT,10002) TNP(1),TNP(2),TNP(3),TNP(4)
         CALL STTPUT(OUTPUT,STAT)
      ENDIF
C
C  copy descriptors + calculate new dynamic range of result frame
      CALL DSCUPT(IMN(1),IMNOC,' ',STAT)
      CALL MYMX(MADRID(PNTRC),SIZEC,CUTS(3))
      CUTS(1) = CUTS(3)
      CUTS(2) = CUTS(4)
      CALL STDWRR(IMNOC,'LHCUTS',CUTS,1,4,UNI,STAT)
C
      CALL STSEPI
C
C  Formats
10000 FORMAT(' low and high cuts',G12.5,' ,',G12.5)
10001 FORMAT('         EXPTIM',G12.5,', BGV',G12.5,', rej_pix',I8)
10002 FORMAT('pixels with N=0:',I8,', N=1:',I8,', N=2:',I8,', N>2:',I8)
      END

      SUBROUTINE COPULA(IOFF,IMNO,OFFSET,A,NP)
C
      IMPLICIT NONE
C
      INTEGER     IOFF,IMNO,OFFSET,NP
      INTEGER     IAV,STAT
C
      REAL        A(1)
C
      CALL STFGET(IMNO,OFFSET,NP,IAV,A(IOFF),STAT)
      RETURN
      END

      SUBROUTINE COPUL(A,B,NLINE,NP)
C
      IMPLICIT NONE
C
      INTEGER     IOFF,N,NLINE,NP
C
      REAL        A(1),B(1)
C
      IOFF = (NLINE-1) * NP
      DO 800 N=1,NP
         B(N+IOFF) = A(N)
800   CONTINUE
C
      RETURN
      END

      SUBROUTINE AVWIN1(T,TOTAL,NP)
C
      IMPLICIT NONE
C
      REAL        T(1)
      REAL        V,VM,ER,ERR,WG,AM,ALP,F
      REAL        VAL(20),BGV(20),RMS(20),RFC(20)
      REAL        AJCUL(20),AJCUH(20),SNOISE,BGERR
      REAL        EXPTI(20)
      REAL        SQRT                        !system function
C
      LOGICAL     NOK(20)
C
      INTEGER     TOTAL,NP
      INTEGER     IORD(20),NORD(20),REJ(20),TNP(4)
      INTEGER     I,IDX,II,IM1,IO
      INTEGER     N1,N2,NH,NN,NX
      INTEGER     TNP1,TNP2,TNP3,TNP4
C
      COMMON  /MYREAL/  VAL,BGV,RMS,RFC,AJCUL,AJCUH,SNOISE,BGERR,EXPTI
      COMMON  /MYINT/   IORD,NORD,REJ,TNP
C
      EQUIVALENCE    (NORD(1),N1),(NORD(2),N2)
      EQUIVALENCE    (TNP(1),TNP1),(TNP(2),TNP2),
     +               (TNP(3),TNP3),(TNP(4),TNP4)
C
      ALP = 0.
C
C  loop through pixels in each line
      DO 5000 NX=1,NP
         NN = 0                          !init no. of valid pixels
         IDX = NX                        !init offset in working area
C
         DO 400 I=1,TOTAL
            V = T(IDX)
            NOK(I) = V.LT.AJCUL(I) .OR. V.GT.AJCUH(I)      !set NOK, if outside
C
            IF (.NOT.NOK(I)) THEN
               V = V - BGV(I)
               VAL(I) = RFC(I) * V       !fill value stack
             NN = NN + 1                 !increment no. of valid pixels
             IF (V.LT.0.) THEN
                F = 0.
             ELSE
                F = SNOISE * V
             ENDIF
             RMS(I) = RFC(I) * SQRT(BGERR + F)      !estimate deviation
            ENDIF
C
            IDX = IDX + NP                          !increase offset
400      CONTINUE
C
C  order pixel values + find median
         IF (NN.GE.3) THEN
            GOTO 2000                       !more than 2 pixels o.k.
         ELSE
            GOTO (1400,1410,1450),NN+1      !branch acc. to no. of good pixels
         ENDIF
C
C  no good pixels, take 1. frame
1400     NOK(1) = .FALSE.
         TNP1 = TNP1 + 1
         GOTO 3900
C
C  only one good pixel, take it
1410     TNP2 = TNP2 + 1
         GOTO 3900
C
C  two good pixels, test them
1450     II = 1
         DO 1600 I=1,TOTAL
            IF (.NOT.NOK(I)) THEN
               NORD(II) = I
               II = II + 1
            ENDIF
1600     CONTINUE
C
         IF (ABS(VAL(N1)-VAL(N2)).LE.RMS(N1)+RMS(N2)) THEN
            TNP3 = TNP3 + 1
         ELSE
            TNP2 = TNP2 + 1
            IF (ABS(VAL(N1)-ALP).LE.ABS(VAL(N2)-ALP)) THEN
               NOK(N2) = .TRUE.
            ELSE
               NOK(N1) = .TRUE.
            ENDIF
         ENDIF
         GOTO 3900
C
C  more than 2 pixels, compare them
2000     DO 2500 I=1,TOTAL
            IF (.NOT.NOK(I)) THEN
               V = VAL(I)
               IO = 0
               IF (I.EQ.1) GOTO 2400
               IM1 = I - 1
               DO 2200 II=1,IM1
                  IF (.NOT.NOK(II)) THEN
                     IF (V.GT.VAL(II)) THEN
                        IO = IO + 1
                     ELSE
                        IORD(II) = IORD(II) + 1
                     ENDIF
                  ENDIF
2200           CONTINUE
2400           IORD(I) = IO
            ENDIF
2500     CONTINUE
C
C  update NORD
         DO 2600 I=1,TOTAL
            IF (.NOT.NOK(I)) THEN
               II = IORD(I) + 1
               NORD(II) = I
            ENDIF
2600     CONTINUE
C
C  and get median
         NH = (NN+1) / 2
         IM1 = NH - 1
         I = NORD(NH)
         VM = VAL(I)
         ER = RMS(I)
         NH = NN
C
C  now remove all pixels with high residuals (not close to median...)
         DO 3000 I=1,IM1
            II = NORD(I)                                !start at bottom...
            ERR = RMS(II) + ER
            IF (ABS(VAL(II)-VM).GT.ERR) THEN    !discard, if value differs much
               NN = NN - 1
               NOK(II) = .TRUE.
            ENDIF
            II = NORD(NH)                               !and move to top...
            NH = NH - 1
            ERR = RMS(II) + ER
            IF (ABS(VAL(II)-VM).GT.ERR) THEN
               NN = NN - 1
               NOK(II) = .TRUE.
            ENDIF
3000     CONTINUE
C
         NN = MIN(4,NN+1)
         TNP(NN) = TNP(NN) + 1
C
C  calculate final flux
3900     WG = 0.
         AM = 0.
         DO 4000 I=1,TOTAL
            IF (.NOT.NOK(I)) THEN
               WG = WG + EXPTI(I)
               AM = AM + EXPTI(I)*VAL(I)
            ELSE
               REJ(I) = REJ(I) + 1
            ENDIF
4000     CONTINUE
         ALP = AM / WG
         T(NX) = ALP
C
5000  CONTINUE                               !end of loop through pixels
      RETURN
      END

      SUBROUTINE AVWIN2(ACTION,T,TOTAL,NP)
C
      IMPLICIT NONE
C
      REAL        T(1)
      REAL        VAL(20),V,O(11),RVAL
      REAL        VVAL(122)                         !just to pad the COMMON ...
C
      INTEGER     TOTAL,NP
      INTEGER     I,IDX,NN
      INTEGER     NH,NX,NH1
C
      CHARACTER   ACTION*3
C
      COMMON   /MYREAL/  VAL,VVAL
C
C ***
C
C  handle maximum search
C
C ***
C
      IF (ACTION(1:2).EQ.'MA') THEN
C
C  loop through pixels in each line
         DO 800 NX=1,NP
C
            IDX = NX                           !init offset in working area
            DO 200 I=1,TOTAL
               VAL(I) = T(IDX)
               IDX = IDX + NP                  !increase offset
200         CONTINUE
C
            V = VAL(1)
            DO 400 I=2,TOTAL
               IF (V.LT.VAL(I)) V = VAL(I)
400         CONTINUE
            T(NX) = V
800      CONTINUE
C
C ***
C
C  handle minimum search
C
C ***
C
      ELSE IF (ACTION(1:2).EQ.'MI') THEN
C
C  loop through pixels in each line
         DO 1800 NX=1,NP
C
            IDX = NX                           !init offset in working area
            DO 1200 I=1,TOTAL
               VAL(I) = T(IDX)
               IDX = IDX + NP                  !increase offset
1200        CONTINUE
C
            V = VAL(1)
            DO 1400 I=2,TOTAL
               IF (V.GT.VAL(I)) V = VAL(I)
1400        CONTINUE
            T(NX) = V
1800     CONTINUE
C
C ***
C
C  handle median search
C
C ***
C
      ELSE
C
C  loop through pixels in each line
         NH = (TOTAL+1) / 2
         NH1 = NH + 1
         DO 5000 NX=1,NP
C
            IDX = NX                           !init offset in working area
            DO 2200 I=1,TOTAL
               VAL(I) = T(IDX)
               IDX = IDX + NP                  !increase offset
2200        CONTINUE
C
C  order pixel values
            O(1) = VAL(1)
C
C  build up lower half of ordered array (including median index)
            DO 3600 I=2,NH
               NN = I
               RVAL = VAL(NN)
C
3100           IF ((NN.LE.1) .OR. (O(NN-1).LE.RVAL) ) GOTO 3200
               O(NN) = O(NN-1)
               NN = NN - 1
               GOTO 3100
C
3200           O(NN) = RVAL
3600        CONTINUE
C
C  now only merge upper half of array VAL into ordered array O if necessary
            DO 4000 I=NH1,TOTAL
               RVAL = VAL(I)
               NN = NH
               IF (RVAL.GE.O(NN)) GOTO 4000
C
3700           IF ( (NN.LE.1) .OR. (O(NN-1).LE.RVAL) ) GOTO 3800
               O(NN) = O(NN-1)
               NN = NN - 1
               GOTO 3700
C
3800           O(NN) = RVAL
4000        CONTINUE
C
            T(NX) = O(NH)
5000  CONTINUE
      ENDIF
C
      RETURN
      END