C @(#)comsky.for 17.1.1.1 (ES0-DMD) 01/25/02 17:19:19 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 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding 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 COMSKY C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.IDENT program COMSKY A. Lauberts ESO 890707 C revised from COMB C.KEYWORDS C C C bulk data frames, restitution, sky background C C.PURPOSE C C Given two image frames A and B (= subset of A) and a table file with C preliminary calibration coefficients. C Get: C C SKYBGR - Sky backgr plane through intensity mean levels left of highest C bins x SKYF in histogram of frame A, at 8 rectangles surrounding C frame B and inside two concentric circles, C GALPOP - Start PDS density of upper 5% and 10% populations in histogram C of B, C GALCEN - Centroid of region 10% to 5% highest PDS density population in B. C C Replace (option) intersection of A and B with B. C C Note: if A = B then only part of these operations are done C C.USAGE MIDAS applic program C execute as @@ COMSKY P1 P2 P3 P4 P5 C C.PARAMETERS C Input: C C P1 FRAMEA C*8 first frame A C P2 FRAMEB C*8 second frame B (= centred subset of A) C P3 TABLEC C*8 table file (one row, 9 columns) with preliminary C calibration coefficients for the characteristic curve C (derived in DTOM according to Llebaria from spots only) C P4 METH I*4 method 2 (PDS scan before 14.4.1984, C method 1 (PDS scan after 13.4.1984, C descriptor O_TIME set = 1984.5) C method 0 (same as 1, but no pixel replacement) C method 3 (only one frame, 2nd frame set = 1st frame, C only write new descriptor SKYBGR to A, C default) C P5 SKYF R*4 sky factor, default=1.4 C C Output: Restituted image files + new descriptors C C.COMMENTS Note PDS was not working well at slow speed (<=10) 1987-88. C----------------------------------------------------------------------- IMPLICIT INTEGER (A-Z) C CHARACTER*60 FRAMEA,FRAMEB,TABLEC CHARACTER*24 CUNITA,CUNITB INTEGER HIST(4096),NPIXA(2),NPIXB(2),ISEQ(8),COL(10) INTEGER MADRID REAL T(4096),STEPA(2),STEPB(2),STARTA(2),STARTB(2), + SKY(4),GAL(2),CEN(9),STARTO(2), + A,B,C,D,DI,D1,D2,C1,C2,C3,C4,CONST(9), + X(8),Y(8),Z(8),SKYF, + XX(8),YY(8),ZZ(8),DEV(8),DLIM DOUBLE PRECISION OTIME CHARACTER*21 IDENTA,IDENTB CHARACTER HISTORY*19, SKYVAL*80 LOGICAL NULL(10) C COMMON /VMR/MADRID(1) COMMON /COM/T,HIST,MAXD,I00,J00,IJ1,IJ2 C C set up MIDAS environment C CALL SXSPRO('COMSKY') DO I=1,10 COL(I)=I ENDDO C C get keywords C CALL SXKGET('IN_A','C',1,60,IAV,FRAMEA,STAT) CALL SXKGET('IN_B','C',1,60,IAV,FRAMEB,STAT) CALL SXKGET('INPUTC','C',1,60,IAV,TABLEC,STAT) CALL SXKGET('INPUTI','I',1,1,IAV,METH,STAT) CALL SXKGET('INPUTR','R',1,1,IAV,SKYF,STAT) C C read calibration coefficients from one row table file C CALL TXTOPN(TABLEC,STAT) CALL TXRGET(TABLEC,'R*4',1,9,COL,CONST,NULL,STAT) C+ C C get method of pixel restoration C METH=0: no pix replacement (new amplifier) C METH=1: pixel replacement (new amplifier on PDS) C METH=2: pixel replacement + interpolation C METH=3: (default) FRAMEB = FRAMEA (only frame A needed) C DLIM upper limit of PDS density units C SKYF sky background computation < approx sky*SKYF (default=1.4) C C- IF(METH.EQ.2.OR.METH.EQ.3) THEN CALL SXDGET(FRAMEA,'O_TIME','R*8',1,1,IAV,OTIME,STAT) ELSE OTIME=1984.5D0 ENDIF IF(OTIME.GT.1983.5D0) THEN IF(METH.NE.0.AND.METH.NE.3) METH=1 MAXD=4096 DLIM=4096. DI=.00125 ELSE IF(METH.NE.0.AND.METH.NE.3) METH=2 MAXD=1024. DI=.005 DLIM=800. ENDIF C C look-up table D to I according to the Llebaria formula C with coefficients C1-C5 derived in DTOM with spots only: C C log I = C3 log(D-C1) + C4 log(C2-D) + C5, C C where C C1 = fog density C C2 = saturation density C C3 = A factor C C4 = B factor C C5 = zeropoint (not used here) C C1=CONST(2) C2=CONST(3) C3=CONST(4) C4=CONST(5) D1=C1+.05 D2=C2-.05 DO I=1,MAXD D=(I-1)*DI D=AMIN1(D2,AMAX1(D1,D)) T(I)=(D-C1)**C3*(C2-D)**C4 ENDDO IF(METH.EQ.3.OR.FRAMEA.EQ.FRAMEB) THEN METH=3 FRAMEB=FRAMEA ENDIF WRITE(6,101) METH,MAXD,DLIM,SKYF 101 FORMAT(' METH=',I2,' MAXD=',I5,' DLIM=',F6.0,' SKYF=',F4.1) CALL SXIGET(FRAMEA,'IO',2,NAXISA,NPIXA,STARTA,STEPA,IDENTA, + CUNITA,PNTRA,STAT) CALL SXTPUT(IDENTA,STAT) NDIMA=NPIXA(1)*NPIXA(2) IF(METH.EQ.0.OR.METH.EQ.3) THEN DO I=1,2 NPIXB(I)=NPIXA(I)/3 STEPB(I)=STEPA(I) GAL(I)=MAXD/2 ENDDO IF(METH.EQ.3) GOTO 30 ENDIF HISTORY(1:8)=FRAMEA IF(METH.NE.0) THEN HISTORY(9:11)=' + ' HISTORY(12:19)=FRAMEB ELSE HISTORY(9:19)=' ' ENDIF CALL SXIGET(FRAMEB,'I',2,NAXISB,NPIXB,STARTB,STEPB, + IDENTB,CUNITB,PNTRB,STAT) NDIMB=NPIXB(1)*NPIXB(2) IF(IDENTA.EQ.IDENTB) GOTO 20 CALL SXTPUT(' *** NOT SAME OBJECTS',STAT) GOTO 999 20 IF(NPIXA(1).GT.NPIXB(1)) GOTO 30 CALL SXTPUT(' *** A <= B',STAT) GOTO 999 C+ C C map output file C centre coords in STARTO C set temporarily new origo at centre by redefining STARTA and STARTB C C- 30 X(1)=STARTA(1) Y(2)=STARTA(2) STARTA(1)=-FLOATJ(NPIXA(1)-1)*STEPA(1)/2. STARTA(2)=-FLOATJ(NPIXA(2)-1)*STEPA(2)/2. STARTO(1)=X(1)-STARTA(1) STARTO(2)=Y(2)-STARTA(2) STARTB(1)=-FLOATJ(NPIXB(1)-1)*STEPB(1)/2. STARTB(2)=-FLOATJ(NPIXB(2)-1)*STEPB(2)/2. C update descriptors 40 CALL SXDPUT(FRAMEA,'HISTORY','C',HISTORY,-1,30,STAT) IF (METH.NE.2) THEN CALL SXDPUT(FRAMEA,'O_TIME','R*8',OTIME,1,1,STAT) ENDIF C+ C C compute sky background in 8 rectangles surrounding frame B, C also inside two concentric circles with diams (NPIXB+3*NPIXA)/4 C and NPIXA C C I0J3 I1J3 I2J3 I3J3 C X7Y7 X6Y6 X5Y5 C I0J2 I1J2 I2J2 I3J2 C X8Y8 X0Y0 X4Y4 C I0J1 I1J1 I2J1 I3J1 C X1Y1 X2Y2 X3Y3 C I0J0 I1J0 I2J0 I3J0 C C- 45 I0= 1 I1=(NPIXA(1)-NPIXB(1))/2+1 I2= NPIXA(1)-I1 I3= NPIXA(1) J0= 1 J1=(NPIXA(2)-NPIXB(2))/2+1 J2= NPIXA(2)-J1 J3= NPIXA(2) C centre I00=I3/2+1 J00=J3/2+1 C squared radii of concentric circles IJ1=((I2+I3)**2+(J2+J3)**2)/32 IJ2=(I3**2+J3**2)/8 C compute mean x,y coordinates in rectangles X(1)=(STARTA(1)+STARTB(1))/2. Y(1)=(STARTA(2)+STARTB(2))/2. X(2)=0. Y(2)=(STARTA(2)+Y(1))/2. X(3)=-X(1) Y(3)=Y(1) X(4)=-(STARTA(1)+X(1))/2. Y(4)=0. X(5)=X(3) Y(5)=-Y(3) X(6)=0. Y(6)=-Y(2) X(7)=X(1) Y(7)=Y(5) X(8)=-X(4) Y(8)=0. CALL SKYBAK(NDIMA,NPIXA,MADRID(PNTRA),Z(1),SKYF,I0,J0,I1,J1) CALL SKYBAK(NDIMA,NPIXA,MADRID(PNTRA),Z(2),SKYF,I1,J0,I2,J1) CALL SKYBAK(NDIMA,NPIXA,MADRID(PNTRA),Z(3),SKYF,I2,J0,I3,J1) CALL SKYBAK(NDIMA,NPIXA,MADRID(PNTRA),Z(4),SKYF,I2,J1,I3,J2) CALL SKYBAK(NDIMA,NPIXA,MADRID(PNTRA),Z(5),SKYF,I2,J2,I3,J3) CALL SKYBAK(NDIMA,NPIXA,MADRID(PNTRA),Z(6),SKYF,I1,J2,I2,J3) CALL SKYBAK(NDIMA,NPIXA,MADRID(PNTRA),Z(7),SKYF,I0,J2,I1,J3) CALL SKYBAK(NDIMA,NPIXA,MADRID(PNTRA),Z(8),SKYF,I0,J1,I1,J2) WRITE(SKYVAL,201) Z 201 FORMAT(' SKYVALUES: ',8F8.1) CALL SXTPUT(SKYVAL,STAT) M=8 DO I=1,8 ISEQ(I)=1 ENDDO 50 II=0 DO I=1,8 IF(ISEQ(I).EQ.1) THEN II=II+1 XX(II)=X(I) YY(II)=Y(I) ZZ(II)=Z(I) ENDIF ENDDO CALL PLANE(XX,YY,ZZ,SKY,M,DEV) IF(SKY(4).LT.SKY(1)*0.01) GOTO 53 IF(M.LE.6) GOTO 52 DEV1=0. J=1 DO I=1,8 IF(ISEQ(I).EQ.1.AND.DEV(I).GT.DEV1) THEN J=I DEV1=DEV(I) ENDIF ENDDO ISEQ(J)=0 WRITE(SKYVAL,203) Z(J) 203 FORMAT(' *** BAD SKY, EXCLUDE VALUE',F8.1) CALL SXTPUT(SKYVAL,STAT) M=M-1 GOTO 50 52 CALL SEQUENCE(ZZ,ISEQ,M) SKY(1)=ZZ(ISEQ(3)) SKY(2)=0. SKY(3)=0. WRITE(SKYVAL,202) SKY(1) 202 FORMAT(' *** TAKE ONLY MEDIAN SKY =',F8.1) CALL SXTPUT(SKYVAL,STAT) 53 CALL SXDPUT(FRAMEA,'SKYBGR','R',SKY,1,3,STAT) IF(METH.EQ.3) GOTO 999 C C get location of central galaxy and replace inner region of A with B C CALL EXCHANGE(NDIMA,NDIMB,NPIXA,NPIXB,STARTA,STARTB,STEPA,STEPB, + MADRID(PNTRA),MADRID(PNTRB),SKY,GAL,CEN,DLIM,METH) DO I=1,2 CEN(I)=CEN(I)+STARTO(I) ENDDO CALL SXDPUT(FRAMEA,'GALPOP','R',GAL,1,2,STAT) CALL SXDPUT(FRAMEA,'GALCEN','R',CEN,1,9,STAT) C 999 CALL SXSEPI END C SUBROUTINE SKYBAK(NDIMA,NPIXA,A,SKY,SKYF,I1,J1,I2,J2) C+ C C compute sky background in region A[I1,J1,I2,J2] C also inside two concentric circles with squared radii IJ1 and IJ2 C C- COMMON /COM/T,HIST,MAXD,I00,J00,IJ1,IJ2 INTEGER HIST(4096),NPIXA(2),HJ,HN REAL A(NDIMA),T(4096) C+ C form histogram of subframe C- DO 11 N=1,MAXD 11 HIST(N)=0 DO 3 J=J1,J2 IJ0=(J-J00)**2 LOCA=(J-1)*NPIXA(1)+I1 DO 2 I=I1,I2 IJ=IJ0+(I-I00)**2 IF(IJ.LT.IJ1.OR.IJ.GT.IJ2) GOTO 2 N=MIN0(MAXD,JIFIX(A(LOCA)+1.)) HIST(N)=HIST(N)+1 2 LOCA=LOCA+1 3 CONTINUE C+ C find highest bin C- MAX=0 DO 4 J=1,MAXD HJ=HIST(J) IF(MAX.GT.HJ) GOTO 4 MAX=HJ HN=J 4 CONTINUE C+ C sky = intensity mean of bins left of highest bin x SKYF C- X=T(HN)*SKYF DO L=1,MAXD IF(T(L).GT.X) GOTO 5 ENDDO 5 SKY=0. N=0 DO J=1,L HJ=HIST(J) SKY=SKY+HJ*T(J) N=N+HJ ENDDO SKY=SKY/FLOATJ(N) C C convert back to density C DO J=1,L IF(T(J).GT.SKY) GOTO 6 ENDDO J=L 6 IF(T(J).NE.T(J-1)) THEN SKY=J-(T(J)-SKY)/(T(J)-T(J-1))-1. ELSE SKY=AMAX1(0.,J-1.5) ENDIF RETURN END C SUBROUTINE EXCHANGE(NDIMA,NDIMB,NPIXA,NPIXB,STARTA,STARTB, + STEPA,STEPB,A,B,SKY,GAL,CEN,DLIM,METH) COMMON /COM/T,HIST,MAXD CHARACTER OFFC*80 INTEGER HIST(4096),NPIXA(2),NPIXB(2),H(7),HN(7),COUNT,STAT REAL A(NDIMA),B(NDIMB),T(4096),STARTA(2),STARTB(2),STEPA(2), + STEPB(2),CEN(9),GAL(2),SKY(4) C C form histogram of B C DO 8 N=1,MAXD 8 HIST(N)=0 DO 22 I=1,NDIMB N=MIN0(MAXD,JIFIX(B(I)+1.)) HIST(N)=HIST(N)+1 22 CONTINUE C C find max non-zero bin C DO MAX=MAXD,1,-1 IF(HIST(MAX).NE.0) GOTO 24 ENDDO C C find start position of upper 10% population in histogram of B C 24 J=NDIMB*0.9 K=0 DO 9 N=1,MAXD K=K+HIST(N) IF(K.GT.J) GOTO 10 9 CONTINUE 10 GAL1=FLOATJ(N) N1=N+1 C find start position of upper 5% population in histogram of B J=NDIMB*0.95 DO 91 N=N1,MAXD K=K+HIST(N) IF(K.GT.J) GOTO 101 91 CONTINUE 101 GAL2=FLOATJ(N) GAL(1)=GAL2 GAL(2)=GAL1 C C find unweighted centroid of pixels of values GAL1 to GAL2 in B C MX=0 MY=0 N=0 LOCB=1 DO 13 NY=1,NPIXB(2) DO 12 NX=1,NPIXB(1) BL=B(LOCB) IF(BL.LT.GAL1.OR.BL.GT.GAL2) GOTO 12 MX=MX+NX MY=MY+NY N=N+1 12 LOCB=LOCB+1 13 CONTINUE C C world coordinates in microns (populations 1,2) related to frame A C DO I=1,9 CEN(I)=0. ENDDO IF(N.GT.9) THEN GAL1 = FLOATJ(MX)/FLOATJ(N) GAL2 = FLOATJ(MY)/FLOATJ(N) J1 = (GAL1-1.)*STEPB(1)+STARTB(1) J2 = (GAL2-1.)*STEPB(2)+STARTB(2) IF(J1*J1+J2*J2.GT.10000.OR.GAL(1).LT.SKY(1)+50.) THEN WRITE(OFFC,200) J1,J2,MAX,GAL 200 FORMAT(' OFF CENTER',I5,',',I5,' MAX PDS',I5,' GALPOP',2F8.1, + 17X,'*****') ELSE CEN(1)=J1 CEN(2)=J2 WRITE(OFFC,201) J1,J2,MAX,GAL 201 FORMAT(' OFF CENTER',I5,',',I5,' MAX PDS',I5,' GALPOP',2F8.1) ENDIF CALL SXTPUT(OFFC,STAT) ELSE WRITE(OFFC,202) 202 FORMAT(' CENTROID ILL DEFINED',54X,'*****') CALL SXTPUT(OFFC,STAT) ENDIF IF(METH.EQ.0.OR.METH.EQ.3) RETURN LSKIP=4 LSKIP1=LSKIP+1 GOTO (25,26) METH C+ C C replace inner square region of A with B C no interpolation across pixels C skip four first (usually bad) scan lines in frame B C also skip last (bad) line in B C C- 25 LOCB=1+NPIXB(1)*LSKIP J1=(NPIXA(2)-NPIXB(2))/2+LSKIP1 J2=J1+NPIXB(2)-LSKIP1 I1=(NPIXA(1)-NPIXB(1))/2+1 LOC1=(J1-1)*NPIXA(1)+I1 LOC2=LOC1+NPIXB(1)-1 DO J=J1,J2 DO LOCA=LOC1,LOC2 A(LOCA)=B(LOCB) LOCB=LOCB+1 ENDDO LOC1=LOC1+NPIXA(1) LOC2=LOC2+NPIXA(1) ENDDO RETURN C+ C C Restore part of high densities asymmetrically displaced C at too high scan speed; fill in valleys under straight C lines connecting two maxima, again replacing A with B. C Weak point: interference with nearby bright stars C C- 26 L1=1+NPIXB(1)*LSKIP L2=NPIXB(1)*LSKIP1 LOC1=L1-NPIXB(1) LOC2=L1+NPIXB(1) J1=(NPIXA(2)-NPIXB(2))/2+LSKIP1 J2=J1+NPIXB(2)-LSKIP1-1 I1=(NPIXA(1)-NPIXB(1))/2+1 DO J=J1,J2 LOCA=(J-1)*NPIXA(1)+I1 BM1=0. BM2=0. LM1=0 LM2=0 DO L=L1,L2 BL=B(L) IF(BL.GT.DLIM) THEN C find 1st max on original line IF(BL.GT.BM1) THEN BM1=BL LM1=LOCA ENDIF B1=B(LOC1) B2=B(LOC2) IF(BL.LT.AMIN1(B1,B2)) THEN C interpolate between scan lines BL=(B1+B2)/2. C find 2nd max on interpolated line IF(BL.GT.BM2) THEN BM2=BL LM2=LOCA ENDIF ENDIF ENDIF A(LOCA)=BL LOC1=LOC1+1 LOC2=LOC2+1 LOCA=LOCA+1 ENDDO C if valley exists, fill in depression IF(LM1.NE.0.AND.LM2.NE.0.AND.LM1.NE.LM2 + .AND.JIABS(LM1-LM2).LT.15) THEN R=(BM2-BM1)/(LM2-LM1) I=1 IF(LM1.GT.LM2) I=-1 DO L=LM1+I,LM2-I,I A(L)=BM1+R*(L-LM1) ENDDO ENDIF L1=L2+1 L2=L2+NPIXB(1) ENDDO RETURN END C C SUBROUTINE PLANE(X,Y,Z,A,M,DEV) DOUBLE PRECISION B(4,4) REAL X(8),Y(8),Z(8), A(8),DEV(8) CHARACTER*44 APPROX C C LSQ fitting of plane C Z = A1 + A2*X + A3*Y C through M points C N=3 N1=N+1 DO 2 I=1,4 DO 1 J=1,4 1 B(I,J)=0. 2 CONTINUE DO 3 K=1,M A(1)=1. A(2)=X(K) A(3)=Y(K) A(4)=Z(K) DO 31 I=1,4 DO 32 J=1,4 32 B(I,J)=B(I,J)+A(I)*A(J) 31 CONTINUE 3 CONTINUE B44=B(4,4) B34=B(3,4) B24=B(2,4) B14=B(1,4) CALL SIMUL(N,N1,B,A) A(4)=SQRT(ABS(B44-B34*A(3)-B24*A(2)-B14*A(1))/FLOATJ(M)) WRITE(APPROX,200) (A(I), I=1,4) 200 FORMAT(4E11.3) CALL SXTPUT(' PLANE APPROX'//APPROX,ISTAT) DO I=1,M DEV(I)=ABS(Z(I)-A(1)-A(2)*X(I)-A(3)*Y(I)) ENDDO RETURN END C C SUBROUTINE SIMUL(N,N1,A,X) C THE GAUSS-JORDAN COMPLETE ELIMINATION METHOD IS EMPLOYED C WITH THE MAXIMUM PIVOT STRATEGY. DOUBLE PRECISION A, EPS, AIJCK, PIVOT DIMENSION IROW(6), JCOL(6), Y(6), A(4,4), X(N) MAX=N1 EPS=1.E-20 C BEGIN ELIMINATION PROCEDURE 5 DETER=1. DO 18 K=1,N KM1=K-1 C SEARCH FOR THE PIVOT ELEMENT PIVOT=0. DO 10 I=1,N DO 11 J=1,N C SEARCH IROW AND JCOL ARRAYS FOR INVALID PIVOT SUBSCRIPTS IF(K.EQ.1)GOTO 9 DO 8 ISCAN=1,KM1 DO 7 JSCAN=1,KM1 IF(I.EQ.IROW(ISCAN))GOTO 11 IF(J.EQ.JCOL(JSCAN))GOTO 11 7 CONTINUE 8 CONTINUE 9 CONTINUE IF(DABS(A(I,J)).LE.DABS(PIVOT))GOTO 11 PIVOT=A(I,J) IROW(K)=I JCOL(K)=J 11 CONTINUE 10 CONTINUE C INSURE THAT SELECTED PIVOT IS LARGER THAN EPS IF(DABS(PIVOT).GT.EPS)GOTO 13 CALL SXTPUT(' WARNING! SINGULAR MATRIX',ISTAT) RETURN C UPDATE THE DETERMINANT VALUE 13 IROWK=IROW(K) JCOLK=JCOL(K) DETER=DETER*PIVOT C NORMALIZE PIVOT ROW ELEMENTS DO 14 J=1,MAX 14 A(IROWK,J)=A(IROWK,J)/PIVOT C CARRY OUT ELIMINATION A(IROWK,JCOLK)=1./PIVOT DO 18 I=1,N AIJCK=A(I,JCOLK) IF(I.EQ.IROWK)GOTO 18 A(I,JCOLK)=-AIJCK/PIVOT DO 17 J=1,MAX IF(J.NE.JCOLK) A(I,J)=A(I,J)-AIJCK*A(IROWK,J) 17 CONTINUE 18 CONTINUE C ORDER SOLUTION VALUES DO 20 I=1,N IROWI=IROW(I) JCOLI=JCOL(I) X(JCOLI)= A(IROWI,MAX) 20 CONTINUE RETURN END C SUBROUTINE SEQUENCE(SEQ,ISEQ,N) C C returns integer array ISEQ of indices for C N first values according to descending order C in real*4 array SEQ C REAL SEQ(N), T INTEGER ISEQ(N),JSEQ(8) C IF(N.GT.8) RETURN DO I=1,N JSEQ(I)=0 ENDDO DO J=1,N T=-1.E10 DO I=1,N IF(JSEQ(I).EQ.0.AND.T.LT.SEQ(I)) THEN T=SEQ(I) ISEQ(J)=I ENDIF ENDDO JSEQ(ISEQ(J))=J ENDDO RETURN END