C @(#)subsky.for 17.1.1.1 (ESO-DMD) 01/25/02 17:19:22 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