C @(#)comsky.for	19.1 (ES0-DMD) 02/25/03 13:31:21
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
