C @(#)subsky.for	19.1 (ESO-DMD) 02/25/03 13:31:24
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 SUBSKY
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  Program    SUBSKY                        version 1.00           881010
C  Principle proposed by: L. Lucy        ESO - Garching     1988 October 
C  MIDAS implementation:  D. Baade       ESO - Garching
C                                                   1.10           890518
C                                                   1.20           900220
C
C.KEYWORDS
C  Noise reduction, background subtraction 
C
C.PURPOSE
C  Subtract sky from images (model sky by extrapolating reference sky 
C  histogram to full frame). Reduce number of pixels with negative fluxes. 
C
C.ALGORITHM
C  Get the names of input + output frames from IN_A + OUT_A. Use the 
C  supplied histogram (in descriptor HISTOGRAM; descriptor STATISTIC 
C  must also be set as by MIDAS command STATISTIC/IMAGE) as a model of 
C  the background. Scale up this sky reference histogram to area of 
C  full frame. Bin-wise ratio to actual image histogram provides probability 
C  that a pixel in a given bin is due to background. If this probability 
C  is unity, set all pixels in the bin concerned to zero. Otherwise rank 
C  pixels in the bin considered by the mean flux in the ambient m x n image 
C  pixels. If NZERO pixels are probably due to background and therefore 
C  to be zeroed, select NZERO lowest ranking pixels and set them to zero. 
C  Subtract the expectation value computed for all lower bins excluding the 
C  current one from the remaining pixels in the bin. Flux conservation is 
C  not fully maintained. 
C
C.INPUT/OUTPUT
C  the following keys are used:
C
C  IN_A/C/1/60                  input frame (with STAT/IMA previously applied) 
C  OUT_A/C/1/60                 output frame 
C  INPUTI/I/1/2                 dimensions of neighbourhood box  
C
C.REFERENCES TO EXTERNAL SUBROUTINES OTHER THAN THE MIDAS STANDARD INTERFACES
C
C  HSTVLS                       to compute histogram (link with ugenlib)
C
C.VERSIONS
C  1.10  ???    KB   adapted to Portable MIDAS 
C  1.11  900807 RHW  bug fixed in routine INDEXX (sorting of 1 element failed)
C  1.20  910220 DB   a) made box size free parameter, 
C                    b) changed some working buffers from REAL to INTEGER 
C                    c) introduced finer resolution of flux values before 
C                       sorting (subroutine AMBINT)
C                    d) re-introduced original version of subroutine INDEXX
C                       but now paying special attention to case N=1 
C        
C 
C 010201		last modif
C 
C --------------------------------------------------------------------
C
      IMPLICIT NONE
C      
      INTEGER     STAT,NAXIS,NPIX(3),IAV,SHISTO(1000),WTO(2)
      INTEGER     WFR(2),NBIN,ISIZ
      INTEGER     INNO,OUTNO,BOX(2)
      INTEGER     AMBNO,SBINNO,INDXNO,SCR1NO,SCR2NO,POINNO
      INTEGER     UNI(1),NULO,MADRID(1)
      INTEGER*8   INPT,OUTPT,AMBPT,SBINPT,INDXPT,SCR1PT,SCR2PT,POINPT
C      
      DOUBLE PRECISION START(3),STEP(3)
C 
      REAL        STATS(11),CUTS(2)
C
      CHARACTER   INFRM*60,OUTFRM*60,CUNIT*64,IDENT*72
      CHARACTER   INCFLG*1

C 
      INCLUDE     'MID_INCLUDE:ST_DEF.INC'
C 
      COMMON      /VMR/  MADRID
C      
      INCLUDE     'MID_INCLUDE:ST_DAT.INC'
C 
C  set up MIDAS environment + enable automatic error abort
      CALL STSPRO('SUBSKY')
C
C  get name of input image  + map frame 
      CALL STKRDC('IN_A',1,1,60,IAV,INFRM,UNI,NULO,STAT)
      CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            3,NAXIS,NPIX,START,STEP,
     +            IDENT,CUNIT,INPT,INNO,STAT)
      IF ((NAXIS.GT.2) .AND. (NPIX(3).GT.1)) THEN
         CALL STETER
     +   (1,'Maximum image dimension supported is 2...')
      ELSE IF (NAXIS.EQ.1) THEN 
         NPIX(2) = 1
         NAXIS = 2
      ENDIF
C
      ISIZ = NPIX(1)*NPIX(2)
C
C  Read descriptors STATISTIC, HISTOGRAM, WINDOW_FROM, WINDOW_TO, and HIST_BINS
      CALL STDRDR(INNO,'STATISTIC',1,11,IAV,STATS,UNI,NULO,STAT)
      NBIN = STATS(10)
      IF (NBIN.GT.1000) 
     +   CALL STETER
     +   (2,'Cannot handle more than 1000 histogram bins.')
      IF (NBIN.LT.5) 
     +   CALL STETER
     +   (3,'Histograms with fewer than 5 bins do not make sense.')
C 
      CALL STDRDI(INNO,'HISTOGRAM',1,NBIN,IAV,SHISTO,UNI,NULO,STAT) 
      CALL STDRDR(INNO,'HIST_BINS',3,2,IAV,CUTS,UNI,NULO,STAT)
      CALL STDRDI(INNO,'WINDOW_FROM',1,2,IAV,WFR,UNI,NULO,STAT)
      CALL STDRDI(INNO,'WINDOW_TO',1,2,IAV,WTO,UNI,NULO,STAT)
C
C  get name of output frame, map it
      CALL STKRDC('OUT_A',1,1,60,IAV,OUTFRM,UNI,NULO,STAT)
      CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,3,
     2             NPIX,START,STEP,IDENT,CUNIT,OUTPT,OUTNO,STAT)
C 
C  read and test dimensions of neighbourhood box
      CALL STKRDI ('INPUTI',1,2,IAV,BOX,UNI,NULO,STAT)
      IF ((BOX(1)*BOX(2)) .LT. 0) THEN
         CALL STETER(4,'Dimensions of neighbourhood box must be '//
     +                 'non-negative') 
      ENDIF
      IF ((BOX(1)*BOX(2)) .GT. 20) THEN                ! box suspiciously large?
         CALL STTPUT('WARNING: Size of neighbourhood box appears '//
     +                         'very large',STAT)
      ENDIF
C
C  read flag for inclusion or not of central pixel in calculation
C  of flux in box
      CALL STKRDC('INPUTC',1,1,1,IAV,INCFLG,UNI,NULO,STAT)      
C
C *** We need six temporary arrays identical to inframe:
C  1. for the flux in the neighbourhood of each pixel 
      CALL STIPUT('DUMMY1',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +             NAXIS,NPIX,START,STEP,IDENT,CUNIT,AMBPT,
     +             AMBNO,STAT)
C 
C  2. for the number of the sky histogram bin corresponding to each pixel 
      CALL STIPUT('DUMMY2',D_I4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +             NAXIS,NPIX,START,STEP,IDENT,CUNIT,SBINPT,
     +             SBINNO,STAT)
C 
C  3. for the index determined during the sorting procedure 
      CALL STIPUT('DUMMY3',D_I4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +             NAXIS,NPIX,START,STEP,IDENT,CUNIT,INDXPT,
     +             INDXNO,STAT)
C 
C  4. for intermediate storage 
      CALL STIPUT('DUMMY4',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +             NAXIS,NPIX,START,STEP,IDENT,CUNIT,SCR1PT,
     +             SCR1NO,STAT)
C 
C  5. also for intermediate storage 
      CALL STIPUT('DUMMY5',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +             NAXIS,NPIX,START,STEP,IDENT,CUNIT,SCR2PT,
     +             SCR2NO,STAT)
C 
C  6. for pointers to selected image elements
      CALL STIPUT('DUMMY6',D_I4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +             NAXIS,NPIX,START,STEP,IDENT,CUNIT,POINPT,
     +             POINNO,STAT)
C
C  Now comes the real thing:
      CALL REMSKY(MADRID(INPT),MADRID(OUTPT),
     +            MADRID(AMBPT),MADRID(SBINPT),
     +            MADRID(INDXPT),MADRID(POINPT),
     +            MADRID(SCR1PT),MADRID(SCR2PT),
     +            STATS,SHISTO,NPIX,CUTS,WFR,WTO,BOX,INCFLG)
C
C  Copy all descriptors (leave also LHCUTS as it was since no large 
C  changes are expected), update HISTORY:
      CALL STDCOP(INNO,OUTNO,3,' ',STAT)
C
C  Free data
      CALL STSEPI
      END

      SUBROUTINE REMSKY(IN,OUT,AMB,SBIN,INDX,POINT,SCR1,SCR2,
     2                  STATS,SHISTO,NPIX,CUTS,WFR,WTO,BOX,INCFLG)
C
C  
      IMPLICIT NONE
C
      INTEGER   NPICS,SHISTO(1),HISTO(1000),I,J,K,NBIN,STAT
      INTEGER   NSKY,WFR(2),WTO(2),NZERO,NPIX(3),SUBLO(2),BOX(2)
      INTEGER   SBIN(1),INDX(1),POINT(1)
C
      REAL      IN(1),OUT(1),AMB(1),SCR1(1),SCR2(1)
      REAL      VCENT(1000),PBIN(1000)
      REAL      STATS(11),CUTS(1)
      REAL      EXPVAL,V1,V2,SCL,STEP,MIN,MEAN
C
      CHARACTER WARNNG*80,INCFLG*1
C
      MIN=CUTS(1)
      MEAN=STATS(3)
      NBIN=STATS(10)
      STEP=STATS(11)
C
C  Issue warning if area from which sky histogram has been determined appears 
C  rather small:
      NPICS=NPIX(1)*NPIX(2)                        ! full frame
      NSKY=(WTO(1)-WFR(1)+1)*(WTO(2)-WFR(2)+1)     ! sky reference subframe
      SCL=NSKY*100./NPICS
      IF (SCL.LT.5) THEN 
         WRITE (WARNNG,300) SCL
         CALL STTPUT(WARNNG,STAT)
      ENDIF
C
C  Construct identically structured histogram of entire frame
      SUBLO(1)=1
      SUBLO(2)=1
      CALL HSTVLS(IN,2,NPIX,SUBLO,NPIX,CUTS,STEP,NBIN,HISTO)
C
C  Determine central values of the non-excess bins (see command 
C  STATISTIC/IMAGE for definition of excess bins)
      DO 50 J=2,NBIN-1
         VCENT(J)=CUTS(1)+(J-1.5)*STEP
50    CONTINUE
C
C  Compute on the basis of a comparison of the reference sky histogram and 
C  the histogram of the complete frame the probability that the flux 
C  of pixels in a given histogram bin is due to sky (background in general) 
C  only. To this end, scale up sky histogram to area of full frame. 
C  Excess bins are special.
      DO 80 I=2,NBIN-1
         PBIN(I)=0                              ! minimum
         IF (HISTO(I).NE.0) PBIN(I)=SHISTO(I)*100./HISTO(I)/SCL  ! normally
         IF (PBIN(I).GT.1) PBIN(I)=1            ! maximum
80    CONTINUE
      PBIN(1)=0
      PBIN(NBIN)=0
C
C  Determine for each pixel the histogram bin to which its flux corresponds:
C  In order to save time, subtract from pixels in the high excess bin 
C  the expectation value of the sky reference histogram and from the lower 
C  excess bin the lower cutoff flux of the histogram :
      CALL BINS(IN,OUT,SBIN,CUTS,STEP,NPICS,NBIN,MEAN)
C
C  Now start replacing the probable sky pixels. The bins with a probability 
C  of unity of their contributing pixel values being entirely due to 
C  background are easy to handle - signals are completely zeroed. 
      DO 100 I=1,NPICS
         K=SBIN(I)
         IF (PBIN(K).EQ.1) OUT(I)=0            ! all signal due to sky
100   CONTINUE
C
C  For the remainder we procede histogram bin-wise rather than pixel-wise. 
C  Since a criterion for the order in which pixels in a given bin will be 
C  treated is the ambient flux level, we first compute the mean flux in 
C  a (2*BOX(1)+1)*(2*BOX(2)+1) box centered on each pixel:
C 
      CALL AMBINT(AMB,IN,NPIX,STEP,BOX,INCFLG)
C
      DO 400 J=2,NBIN-1                        ! exclude excess bins
         IF (PBIN(J).EQ.1) GOTO 400            ! already handled above 
         NZERO=NINT(HISTO(J)*PBIN(J))          ! number of pixels to be zeroed
         IF (NZERO.NE.0) THEN                  ! some pixels due to sky 
C
            CALL SELECT(SBIN,AMB,INDX,SCR1,SCR2,POINT,J,NPICS)
C
C  Zero the pixels selected, subtract from the remaining pixels in bin J 
C  the expectation value
            DO 360 I=1,NZERO
               K = POINT(I)
               OUT(K)=0                        ! set selected pixels to zero
360         CONTINUE
C
C  Compute expectation value as the histogram-weighted mean over all 
C  bins (excluding the low excess bin) up to the one with the observed flux:
            V1=0                        ! initialize
            V2=0                        ! dito
            DO 370 K=2,J-1
               V1=V1+VCENT(K)*SHISTO(K)
               V2=V2+SHISTO(K)
370         CONTINUE
C 
            EXPVAL=MIN                  ! lower limit to true exp. value
            IF (V2.NE.0) EXPVAL=V1/V2   ! expectation value 
C 
            DO 380 I=NZERO+1,HISTO(J)
               K = POINT(I)
               OUT(K)=IN(K) - EXPVAL    ! subtract expect. value elsewhere 
380         CONTINUE
         ELSE                           ! no pixel to be zeroed 
C
C  If the probability of a certain signal being entirely due to background 
C  is zero, just subtract the mean background flux:
            DO 390 I=1,NPICS
               K=SBIN(I)
               IF (ABS(K-J).LT.0.1) OUT(I)=IN(I)-MEAN
390         CONTINUE
         ENDIF
400   CONTINUE                          ! replacement probability was 1
C
      RETURN
C 
300   FORMAT ('WARNING: Histogram pertains to only ',G10.3,'% of
     + total image area!')
      END

      SUBROUTINE BINS(IN,OUT,SBIN,CUTS,STEP,NPICS,NBIN,MEAN)
C
C  Determine for each pixel its corresponding bin number in the histogram 
C  (conform to definition of bin boundaries used in MIDAS routine HISTVALS)
C  Simultaneously subtract from lower and upper excess bin pixels respectively 
C  lower cutoff and expectation value of sky reference histogram.
C
      IMPLICIT NONE
C
      INTEGER      NPICS,NBIN,I,SBIN(1)
C
      REAL         STEP,IN(1),MAX,MIN,MEAN,OUT(1),CUTS(1)
C
      MIN=CUTS(1)
      MAX=MIN+(NBIN-2)*STEP
C
      DO 100 I=1,NPICS
         IF ((IN(I)-MIN).LT.0) THEN 
            SBIN(I)=1
            OUT(I)=IN(I)-MIN
            GOTO 100
         ENDIF
         IF((IN(I)-MAX).GE.0) THEN 
            SBIN(I)=NBIN
            OUT(I)=IN(I)-MEAN
            GOTO 100
         ENDIF
         SBIN(I)=NINT((IN(I)-MIN)/STEP)+2
100   CONTINUE
C
      RETURN 
      END

      SUBROUTINE AMBINT(AMB,IN,NPIX,STP,BOX,INCFLG)
C
C  Compute for each pixel the mean flux in its neighbourhood of 
C  (2*BOX(1)+1)*(2*BOX(2)+1) pixels (not including the central pixel itself). 
C  Purposely degrade the resolution in flux so that later the sorting by 
C  neighbourhood flux values is somewhat faster. However, the degradation is 
C  limited to 0.01 times the histogram bin width in order to avoid a bias 
C  in the result because otherwise the effective background removal would 
C  be different at the beginning and at the end of the frame. 
C
      IMPLICIT NONE
C
      INTEGER     NPIX(3),NPICS,I,J,K,M,N,NX,NY,BOX(2),BX,BY,P,SIZE,FLAG
C
      REAL        IN(1),AMB(1),STP,STEP
C
      CHARACTER   INCFLG*1
C
      STEP=STP/100.
      NX=NPIX(1)
      NY=NPIX(2)
      NPICS=NX*NY
      BX=BOX(1)
      BY=BOX(2)
      SIZE=(2*BX+1)*(2*BY+1)-1
      IF (SIZE.LT.0.5) SIZE=1                        ! avoid zero division
C
      IF ((INCFLG.EQ.'Y').OR.(INCFLG.EQ.'y')) THEN 
         FLAG=0 
      ELSE
         FLAG=1
      ENDIF
C
C  Initialize array AMB to zero
      DO 10 I=1,NPICS
         AMB(I)=0
10    CONTINUE

C
C  Start computation for inner pixels
      DO 200 I=1+BX,NX-BX
         DO 100 J=1+BY,NY-BY
            P=(J-1)*NX+I
            DO 50 M=-BY,BY
               K=(J-1+M)*NX+I
               DO 25 N=-BX,BX
                  AMB(P)=AMB(P)+IN(K+N)
25             CONTINUE
50          CONTINUE
            AMB(P)=(NINT(((AMB(P)-FLAG*IN(P))/SIZE)/STEP)+0.5)*STEP 
100      CONTINUE
200   CONTINUE
C
C  Finish off at edges ...
      DO 300 I=1+BX,NX-BX
         DO 310 J=1,BY                             ! bottom
            P=(J-1)*NX+I
            DO 320 M=1,J+BY
               K=(M-1)*NX+I
               DO 330 N=-BX,BX
                  AMB(P)=AMB(P)+IN(K+N)
330            CONTINUE
320         CONTINUE
            SIZE=(2*BX+1)*(J+BY)-1
            IF (SIZE.LT.0.5) SIZE=1                ! avoid zero division
            AMB(P)=(NINT(((AMB(P)-FLAG*IN(P))/SIZE)/STEP)+0.5)*STEP
310      CONTINUE      
C
         DO 360 J=NY-BY+1,NY                       ! top 
            P=(J-1)*NX+I
            DO 370 M=J-BY,NY
               K=(M-1)*NX+I
               DO 380 N=-BX,BX
                  AMB(P)=AMB(P)+IN(K+N)
380            CONTINUE
370         CONTINUE
            SIZE=(2*BX+1)*(NY-J+1+BY)-1
            IF (SIZE.LT.0.5) SIZE=1                ! avoid zero division 
            AMB(P)=(NINT(((AMB(P)-FLAG*IN(P))/SIZE)/STEP)+0.5)*STEP
360      CONTINUE      
300   CONTINUE
C 
      DO 400 J=1+BY,NY-BY
         DO 410 I=1,BX                             ! left margin
            P=(J-1)*NX+I
            DO 420 M=-BY,BY
               K=(J-1+M)*NX+I
               DO 430 N=1,I+BX
                  AMB(P)=AMB(P)+IN(K+N)
430            CONTINUE
420         CONTINUE
            SIZE=(I+BX)*(2*BY+1)-1        
            IF (SIZE.LT.0.5) SIZE=1                ! avoid zero division
            AMB(P)=(NINT(((AMB(P)-FLAG*IN(P))/SIZE)/STEP)+0.5)*STEP
410      CONTINUE
C
         DO 460 I=NX-BX+1,NX                       ! right margin
            P=(J-1)*NX+I
            DO 470 M=-BY,BY
               K=(J-1+M)*NX+I
               DO 480 N=1,I+BX
                  AMB(P)=AMB(P)+IN(K+N)
480            CONTINUE
470         CONTINUE
            SIZE=(NX-I+1+BX)*(2*BY+1)-1        
            IF (SIZE.LT.0.5) SIZE=1                ! avoid zero division
            AMB(P)=(NINT(((AMB(P)-FLAG*IN(P))/SIZE)/STEP)+0.5)*STEP
460      CONTINUE
400   CONTINUE
C
C  ... and corners
      DO 511 I=1,BX                                ! left
         DO 512 J=1,BY                             ! lower
            P=(J-1)*NX+I
            DO 513 M=1,J+BY
               K=(M-1)*NX+I
               DO 514 N=1,I+BX
                  AMB(P)=AMB(P)+IN(K+N)
514            CONTINUE
513         CONTINUE
            SIZE=(I+BX)*(J+BY)-1
            IF (SIZE.LT.0.5) SIZE=1                ! avoid zero division 
            AMB(P)=(NINT(((AMB(P)-FLAG*IN(P))/SIZE)/STEP)+0.5)*STEP
512      CONTINUE
511   CONTINUE
C
      DO 611 I=1,BX                                ! left
         DO 612 J=NY-BY+1,NY                       ! upper
            P=(J-1)*NX+I
            DO 613 M=J-BY,NY
               K=(M-1)*NX+I
               DO 614 N=1,I+BX
                  AMB(P)=AMB(P)+IN(K+N)
614            CONTINUE
613         CONTINUE
            SIZE=(I+BX)*(NY-J+1+BY)-1
            IF (SIZE.LT.0.5) SIZE=1                ! avoid zero division 
            AMB(P)=(NINT(((AMB(P)-FLAG*IN(P))/SIZE)/STEP)+0.5)*STEP
612      CONTINUE
611   CONTINUE
C
      DO 711 I=NX-BX+1,NX                          ! right
         DO 712 J=1,BY                             ! lower
            P=(J-1)*NX+I
            DO 713 M=1,J+BY
               K=(M-1)*NX+I
               DO 714 N=I-BX,NX
                  AMB(P)=AMB(P)+IN(K+N)
714            CONTINUE
713         CONTINUE
            SIZE=(NX-I+1+BX)*(J+BY)-1
            IF (SIZE.LT.0.5) SIZE=1                ! avoid zero division 
            AMB(P)=(NINT(((AMB(P)-FLAG*IN(P))/SIZE)/STEP)+0.5)*STEP
712      CONTINUE
711   CONTINUE
C
      DO 811 I=NX-BX+1,NX                          ! right
         DO 812 J=NY-BY+1,NY                       ! upper
            P=(J-1)*NX+I
            DO 813 M=J-BY,NY
               K=(M-1)*NX+I
               DO 814 N=I-BX,NX
                  AMB(P)=AMB(P)+IN(K+N)
814            CONTINUE
813         CONTINUE
            SIZE=(NX-I+1+BX)*(NY-J+1+BY)-1
            IF (SIZE.LT.0.5) SIZE=1                ! avoid zero division 
            AMB(P)=(NINT(((AMB(P)-FLAG*IN(P))/SIZE)/STEP)+0.5)*STEP
812      CONTINUE
811   CONTINUE
C
      RETURN
      END

      SUBROUTINE SELECT(SBIN,AMB,INDX,SCR1,SCR2,POINT,IBIN,NPICS)
C
C  Select pixels in bin IBIN whose signal is to be zeroed completely. 
C  Selection criterion is the ambient flux level. Therefore, rank 
C  values of pixels according to the ambient flux (in AMB). 
C
      IMPLICIT NONE
C
      INTEGER    IBIN,NSORT,NPICS,I,J,K
      INTEGER    SBIN(NPICS),INDX(NPICS),POINT(NPICS)
C
      REAL       AMB(NPICS),SCR1(NPICS),SCR2(NPICS)
C
C  Prepare an index array with pointers to all pixels in IBIN
      K=0
      DO 100 I=1,NPICS
         IF (ABS(SBIN(I)-IBIN).LT.0.1) THEN 
            K=K+1
            POINT(K)=I
         ENDIF
100   CONTINUE
      NSORT=K
C
C  Store the relevant NSORT elements of array AMB in the first NSORT 
C  elements of temporary buffer SCR1 so that the sorting is more efficient. 
      DO 200 J=1,NSORT
         K = POINT(J)
         SCR1(J)=AMB(K)
         SCR2(J)=POINT(J)+0.1
200   CONTINUE
C
C  Now sort array SCR2 according to ascending values of the neighbouring 
C  flux level (in AMB). 
      CALL SORT1(NSORT,SCR1,SCR2,INDX)
C
C  Finally, revert from REAL format (array SCR2) to INTEGER (array POINT)
      DO 300 J=1,NSORT
         POINT(J)=INT(SCR2(J))
300   CONTINUE
C
      RETURN
      END

      SUBROUTINE SORT1(N,RA,RB,IWKSP)
C
C  Sort array RB using index IWKSP prepared for array RA in subroutine INDEXX
C  Code adapted from "Numerical Recipes", Chapter 8.3 "Indexing and Ranking", 
C  p. 233-234 (subroutine SORT3).
C
      IMPLICIT NONE
C
      INTEGER     N,J,K
      INTEGER     IWKSP(N)
C
      REAL        RA(N),RB(N)
C
      CALL INDEXX(RA,IWKSP,N)
C
      DO 100 J=1,N
         RA(J)=RB(J)           ! initial contents of buffer RA no longer needed
100   CONTINUE
C
      DO 200 J=1,N
         K = IWKSP(J)
         RB(J)=RA(K)           ! do the desired re-arrangement
200   CONTINUE
C
      RETURN
      END

      SUBROUTINE INDEXX(ARRIN,INDX,N)
C
C  Index array ARRIN according to ascending order of values in SCR1, result in 
C  INDX.
C  Code taken from "Numerical Recipes", Chapter 8.3 "Indexing and Ranking", 
C  p. 233 (Subroutine INDEXX)
C
      IMPLICIT NONE
C
      INTEGER     N,I,J,L,IR,INDXT,INDX(N)
C
      REAL        ARRIN(N),Q
C
      IF (N.LT.1.5) THEN 
         INDX(1)=1
         RETURN
      ENDIF
C
      DO 100 J=1,N
         INDX(J)=J                        ! initialize
100   CONTINUE
C
      L=N/2+1
      IR=N
C
10    CONTINUE
         IF (L.GT.1) THEN
            L = L-1
            INDXT=INDX(L)
            Q=ARRIN(INDXT)
         ELSE
            INDXT=INDX(IR)
            Q = ARRIN(INDXT)
            INDX(IR)=INDX(1)
            IR=IR-1
            IF (IR.LE.1) THEN
               INDX(1)=INDXT
               RETURN                     ! job done 
            ENDIF
         ENDIF
C 
         I=L
         J=L+L
C
20       CONTINUE
         IF (J.LE.IR) THEN
            IF (J.LT.IR) THEN
               IF (ARRIN(INDX(J)).LT.ARRIN(INDX(J+1))) J=J+1
            ENDIF 
            IF (Q.LT.ARRIN(INDX(J))) THEN
               INDX(I) = INDX(J)
               I=J
               J=J+J
            ELSE
               J=IR+1
            ENDIF
C    
            GOTO 20
         ENDIF
C 
         INDX(I)=INDXT
      GOTO 10
C
      END





