C @(#)tdcmap.for	17.1.1.1 (ESO-IPG) 01/25/02 17:47:13
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
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  17:58 - 19 NOV 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C
C.IDENTIFICATION
C
C  program TDCMAP
C
C.PURPOSE
C
C  Execute the command
C  CONV/TAB image = table colx[,coly] [z] refima method
C         for method = FREQ
C                    = POLY deg1[,deg2]
C                    = PLOT min,max
C
C
C.KEYWORDS
C
C  table, image, conversion
C
C.ALGORITHM
C
C  use table interface routines
C
C.COMMENTS
C
C  One has to include ths internal routines into the TD library
C
C.MODIF : make it work for FREQ method MPeron 050990
C 
C  010110		last modif
C 
C-----------------------------------------------------------
C
      SUBROUTINE TDCMAP 
      IMPLICIT NONE
C
      INTEGER MADRID
      INTEGER PARVAL, ID
      INTEGER ICOL(3),NPIX(2)
      INTEGER N,ISTAT,STATUS,NPAR,IAV,IMETH,INDX,NDIM,NAXIS
      INTEGER*8   PNTR
      INTEGER TID,NCOL,NROW,NAC,NAR,NSC,I,I1,I2
      INTEGER ICOLV,NIND,K,L,NTROW,N1
      INTEGER INDEX, KUN, KNUL, IREF, NN,TLEN,CLEN1,CLEN2
C
      DOUBLE PRECISION START(2), STEP(2)
      REAL AXMIN(2),BIN(2),RMIN,RMAX,XMIN,XMAX,YMIN,YMAX
      REAL CUTS(4),RPAR(2),ERROR
C
      DOUBLE PRECISION G, XC
C
      CHARACTER*16 MSG
      CHARACTER*80 TABLE
      CHARACTER*80 IMAGE,REFIMA
      CHARACTER*2 METH
      CHARACTER*17 COLREF(2),COLREV
      CHARACTER*60 WORK
      CHARACTER*72 IDENT
      CHARACTER*80 CUNIT
C
       INCLUDE 'MID_INCLUDE:TABLES.INC'
      COMMON /VMR/MADRID(1)
      COMMON /LS/G(50,50),XC(50),N
       INCLUDE 'MID_INCLUDE:TABLED.INC'
C
      DATA PARVAL/8/
      DATA NPIX/1,1/
      DATA START/0.,0./, STEP /1.,1./
      DATA MSG/'ERR:TMAPITBLxxxx'/
C
C ... get input parameters + default
C
      CALL TDPGET(PARVAL,NPAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 30
      IMAGE  = TPARBF(1)
      TABLE  = TPARBF(3)
      CALL GENLEN(TABLE,TLEN)
      WORK   = TPARBF(4)
      METH   = TPARBF(6)
      CALL FORUPC(METH,METH)
      IF (METH(1:2).EQ.'FR' .AND. TPARBF(7)(1:1).EQ.'+') THEN
          IMETH  = 1
          REFIMA = TPARBF(5)
          GO TO 10
      END IF
      METH   = TPARBF(7)
      CALL FORUPC(METH,METH)
      IF (METH(1:2).EQ.'PO') THEN
          IMETH  = 3
          REFIMA = TPARBF(6)
          COLREV = TPARBF(5)
          CALL STKRDR('INPUTR',1,2,IAV,RPAR,KUN,KNUL,ISTAT)
          GO TO 10
      END IF
      IF (METH(1:2).EQ.'PL') THEN
          IMETH  = 2
          REFIMA = TPARBF(6)
          COLREV = TPARBF(5)
          CALL STKRDR('INPUTR',1,2,IAV,RPAR,KUN,KNUL,ISTAT)
          GO TO 10
      END IF
      CALL STTPUT(' Wrong method ... ',ISTAT)
      GO TO 30
   10 CONTINUE
      INDX   = INDEX(WORK,',')
      IF (INDX.EQ.0) THEN
          NDIM   = 1
          COLREF(1) = WORK
          CALL GENLEN(COLREF(1),CLEN1)
      ELSE
          NDIM   = 2
          COLREF(1) = WORK(1:INDX-1)
          CALL GENLEN(COLREF(1),CLEN1)
          COLREF(2) = WORK(INDX+1:)
          CALL GENLEN(COLREF(2),CLEN2)
      END IF
      CALL STFOPN(REFIMA,D_R4_FORMAT,0,F_IMA_TYPE,IREF,ISTAT)
      CALL STDRDI(IREF,'NAXIS',1,1, NN,    NAXIS,KUN,KNUL,ISTAT)
      CALL STDRDI(IREF,'NPIX', 1,NAXIS,NN, NPIX, KUN,KNUL,ISTAT)
      CALL STDRDD(IREF,'START',1,NAXIS,NN, START,KUN,KNUL,ISTAT)
      CALL STDRDD(IREF,'STEP', 1,NAXIS,NN, STEP, KUN,KNUL,ISTAT)
      CALL STDRDC(IREF,'IDENT',1,1,72,NN,  IDENT,KUN,KNUL,ISTAT)
      CALL STDRDC(IREF,'CUNIT',1,1,80,NN,  CUNIT,KUN,KNUL,ISTAT)
      AXMIN(1) = START(1)
      AXMIN(2) = START(2)
      BIN(1)   = STEP(1)
      BIN(2)   = STEP(2)
C [2.1]
C IF (NAXIS.LT.NDIM) THEN
C   CALL STTPUT(' Wrong reference image ... ',ISTAT)
C   GO TO 1000
C ENDIF
C
C ... init input table
C
      CALL TBTOPN(TABLE,F_I_MODE,TID,ISTAT)
      IF (ISTAT.NE.0) GO TO 30
      CALL TBIGET(TID,NCOL,NROW,NSC,NAC,NAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 30
C
C ... find input column
C
      DO 20 I = 1,NDIM
          CALL TBCSER(TID,COLREF(I),ICOL(I),ISTAT)
          IF (ICOL(I).EQ.-1) THEN
              ISTAT  = ERRCOL
              GO TO 30
          END IF
          CALL TBUGET(TID,ICOL(I),WORK,ISTAT)
CC PN 8/99: CUNIT has to be put into CUNIT(17:32) and CUNIT(33:48)
          I1     = I*16 + 1
          I2     = I1 + 15
CC         WRITE(CBUF,*) 'WORK=',WORK(1:16)   
CC         CALL STTPUT(CBUF,ISTAT)
          CUNIT(I1:I2) = WORK(1:16)
   20 CONTINUE
C
C... handle the different methods
C   
      IF (IMETH.EQ.1) THEN		! Method = FREQUENCY
CC PN 8/99: this has to be written into CUNIT(1:16)
CC          CUNIT(I2+1:I2+16) = 'FREQUENCY'
          CUNIT(1:16) = 'FREQUENCY'
          IF (NDIM.EQ.1) THEN
              IDENT  = 'FREQUENCY OF '//TABLE(1:TLEN)
     1                 //' COL. :'//COLREF(1)(1:CLEN1)
          ELSE
              IDENT  = 'FREQUENCY OF '//TABLE(1:TLEN)
     +                 //' COLS. :'//COLREF(1)(1:CLEN1)
     +                 //' '//COLREF(2)(1:CLEN2)
          END IF

      ELSE IF (IMETH.EQ.2) THEN		! Method = PLOT

              CALL TBCSER(TID,COLREV,ICOLV,ISTAT)
              IF (ICOLV.EQ.-1) THEN
                  ISTAT  = ERRCOL
                  GO TO 30
              END IF
              CALL TBUGET(TID,ICOLV,WORK,ISTAT)
CC PN 8/99: This also has to be written into CUNIT(1:16) instead
CC              CUNIT(I2+1:I2+16) = WORK(1:16)
              CUNIT(1:16) = WORK(1:16)
              IF (NDIM.EQ.1) THEN
                  IDENT  = 'MAP OF '//TABLE(1:TLEN)
     +                     //' COL. :'//COLREF(1)(1:CLEN1)
              ELSE
                  IDENT  = 'MAP OF '//TABLE(1:TLEN)
     +                     //' COLS. :'//COLREF(1)(1:CLEN1)
     +                     //' '//COLREF(2)(1:CLEN2)
              END IF

      ELSE IF (IMETH.EQ.3) THEN		! Method = POLYNOMIAL
              CALL TBCSER(TID,COLREV,ICOLV,ISTAT)
              IF (ICOLV.EQ.-1) THEN
                  ISTAT  = ERRCOL
                  GO TO 30
              END IF
              CALL TBUGET(TID,ICOLV,WORK,ISTAT)
CC PN 8/99 same procedure again ... CUNIT(1:16) instead
CC              CUNIT(I2+1:I2+16) = WORK(1:16)
              CUNIT(1:16) = WORK(1:16)
              IF (NDIM.EQ.1) THEN
                  IDENT  = 'POLY.FIT : '//TABLE(1:TLEN)
     +                     //' COL. :'//COLREF(1)(1:CLEN1)
              ELSE
                  IDENT  = 'POLY.FIT : '//TABLE(1:TLEN)
     +                     //' COLS. :'//COLREF(1)(1:CLEN1)
     +                     //' '//COLREF(2)(1:CLEN2)
              END IF

      ENDIF
C
C ... map image
C
      CALL STIPUT(IMAGE,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .     NDIM,NPIX,START,STEP,IDENT,CUNIT,PNTR,ID,STATUS)

      IF (IMETH.EQ.1) THEN		! Method = FREQUENCY
C
C ... frequency analysis
C
          IF (NDIM.EQ.1) THEN
              CALL FREQ1(TID,ICOL(1),NROW,NPIX,
     .        MADRID(PNTR),AXMIN,BIN,RMIN,RMAX)
          ELSE
              CALL FREQ2(TID,ICOL(1),ICOL(2),NROW,
     .        NPIX(1),NPIX(2),MADRID(PNTR),AXMIN,BIN,RMIN,RMAX)
          END IF

      ELSE IF (IMETH.EQ.2) THEN		! Method = PLOT
C
C ... mapping in plot format
C
           IF (NDIM.EQ.1) THEN
               CALL COUNT1(TID,ICOL(1),ICOLV,NROW,
     .         NPIX,MADRID(PNTR),AXMIN,BIN,RPAR,RMIN,RMAX)
           ELSE
              CALL COUNT2(TID,ICOL(1),ICOL(2),ICOLV,NROW,
     .        NPIX(1),NPIX(2),MADRID(PNTR),AXMIN,BIN,RPAR,RMIN,RMAX)
           END IF

      ELSE IF (IMETH.EQ.3) THEN		! Method = POLYNOMIAL
C
C ... polynomial fit
C
              NTROW  = NAR
C [2.1] START
              IF (NDIM.EQ.1) THEN
                  NIND   = 3
                  ICOL(2) = ICOL(1)
                  ICOL(3) = ICOLV
              ELSE
C [2.1] END
                  NIND   = NDIM + 1
                  ICOL(NIND) = ICOLV
              END IF
              K      = NINT(RPAR(1))
              L      = 0
              IF (NDIM.GT.1) L      = NINT(RPAR(2))
              N      = (K+1)* (L+1)
              CALL LSQFIT(TID,NROW,NTROW,NCOL,NIND,ICOL,K,L,N1,
     +                    XMIN,XMAX,YMIN,YMAX,ERROR)
              IF (NDIM.EQ.1) THEN
                  CALL COMFT2(NPIX(1),1,MADRID(PNTR),AXMIN,BIN,K,0,
     +                         RMIN,RMAX)
              ELSE
                  CALL COMFT2(NPIX(1),NPIX(2),MADRID(PNTR),AXMIN,BIN,
     +                         K,L,RMIN,RMAX)
              END IF

      ENDIF	! End of test on IMETH
C
C ... write cuts
C
      CUTS(1) = RMIN
      CUTS(2) = RMAX
      CUTS(3) = RMIN
      CUTS(4) = RMAX
      CALL STDWRR(ID,'LHCUTS',CUTS,1,4,KUN,ISTAT)
C
C ... end
C
      CALL DSCUPT(ID,ID,' ',ISTAT)
      CALL TBTCLO(TID,ISTAT)
   30 IF (ISTAT.NE.0) THEN
          WRITE (MSG(13:16),9000) ISTAT
          CALL TDERRR(ISTAT,MSG,STATUS)
      END IF
      RETURN
 9000 FORMAT (I4)
      END

      SUBROUTINE FREQ1(TID,ICOL,NROW,NPIX,F,START,STEP,RMIN,RMAX)
C
C compute the frequency
C
      IMPLICIT NONE
      INTEGER  TID
      INTEGER  ICOL
      INTEGER  NROW
      INTEGER  NPIX
      REAL     F(NPIX)
      REAL     START
      REAL     STEP
      REAL     RMIN
      REAL     RMAX
C
      REAL    X
      INTEGER I, STATUS, J
      LOGICAL INULL, ISEL
C
      DO 10 I = 1,NPIX
          F(I)   = 0.
   10 CONTINUE
      RMIN   = 0.
      RMAX   = 0.
      DO 20 I = 1,NROW
          CALL TBSGET(TID,I,ISEL,STATUS)
          CALL TBERDR(TID,I,ICOL,X,INULL,STATUS)
          IF (ISEL.AND.(.NOT.INULL)) THEN
              J      = (X-START)/STEP + 1
              IF (J.GE.1 .AND. J.LE.NPIX) THEN
                  F(J)   = F(J) + 1
                  RMAX   = AMAX1(RMAX,F(J))
              END IF
          END IF
   20 CONTINUE
      RETURN
      END

      SUBROUTINE FREQ2(TID,ICOL1,ICOL2,NROW,
     .                 NPIX1,NPIX2,F,START,STEP,RMIN,RMAX)
C
C compute the frequency
C
      IMPLICIT NONE
      INTEGER  TID
      INTEGER  ICOL1
      INTEGER  ICOL2
      INTEGER  NROW
      INTEGER  NPIX1
      INTEGER  NPIX2
      REAL     F(NPIX1,NPIX2)
      REAL     START(2)
      REAL     STEP(2)
      REAL     RMIN
      REAL     RMAX
C
      INTEGER I, J, STATUS, I1, J1
      LOGICAL INULL1, INULL2, ISEL
      REAL STARTX, STARTY, STEPX, STEPY, X, Y
C
C
      DO 20 I = 1,NPIX2
          DO 10 J = 1,NPIX1
              F(J,I) = 0.
   10     CONTINUE
   20 CONTINUE
      RMIN   = 0.
      RMAX   = 0.
      STARTX = START(1)
      STARTY = START(2)
      STEPX  = STEP(1)
      STEPY  = STEP(2)
      DO 30 I = 1,NROW
          CALL TBSGET(TID,I,ISEL,STATUS)
          CALL TBERDR(TID,I,ICOL1,X,INULL1,STATUS)
          CALL TBERDR(TID,I,ICOL2,Y,INULL2,STATUS)
          IF (ISEL.AND.(.NOT.INULL1).AND.(.NOT.INULL2)) THEN
              J1     = (X-STARTX)/STEPX + 1
              I1     = (Y-STARTY)/STEPY + 1
              IF (J1.GE.1 .AND. J1.LE.NPIX1 .AND. I1.GE.1 .AND. I1.LE.
     +            NPIX2) THEN
                  F(J1,I1) = F(J1,I1) + 1
                  RMAX   = AMAX1(RMAX,F(J1,I1))
              END IF
          END IF
   30 CONTINUE
      RETURN
      END

      SUBROUTINE COUNT2
     .(TID,ICOL1,ICOL2,ICOLR,N,NPIX1,NPIX2,A,
     . START,STEP,RFLAG,RMIN,RMAX)
C
C compute distribution
C
      IMPLICIT NONE
      INTEGER  TID
      INTEGER  ICOL1
      INTEGER  ICOL2
      INTEGER  ICOLR
      INTEGER  N
      INTEGER  NPIX1
      INTEGER  NPIX2
      REAL     A(NPIX1, NPIX2)
      REAL     START(2)
      REAL     STEP(2)
      REAL     RFLAG(2)
      REAL     RMIN
      REAL     RMAX
C
      REAL    X1,Y1,Z1
      LOGICAL ISEL, INULL1, INULL2, INULL3
      REAL    RLOW, RUPP
      INTEGER I, J, STATUS, J1, I1
      CHARACTER*80 STRING
      INTEGER      ISTAT,OVER
C
C
C ... clear counts
C
      RLOW   = RFLAG(1)
      RUPP   = RFLAG(2)
      DO 20 I = 1,NPIX2
          DO 10 J = 1,NPIX1
              A(J,I) = RLOW
   10     CONTINUE
   20 CONTINUE
      RMIN   = RLOW
      RMAX   = RLOW
C
C ... iteration on rows
C
      OVER = 0		! Set overwriting flag to 0
      DO 30 I = 1,N
          CALL TBSGET(TID,I,ISEL,STATUS)
          IF (ISEL) THEN
              CALL TBERDR(TID,I,ICOL1,X1,INULL1,STATUS)
              CALL TBERDR(TID,I,ICOL2,Y1,INULL2,STATUS)
              CALL TBERDR(TID,I,ICOLR,Z1,INULL3,STATUS)
              IF (.NOT.INULL1.AND..NOT.INULL2.AND..NOT.INULL3)THEN
                  J1     = NINT((X1-START(1))/STEP(1)) + 1
                  I1     = NINT((Y1-START(2))/STEP(2)) + 1
                  IF (J1.GE.1 .AND. J1.LE.NPIX1 .AND. I1.GE.1 .AND.
     +                I1.LE.NPIX2) THEN

                      IF (A(J1,I1).NE.RLOW) THEN
                          A(J1,I1) = RUPP
                          OVER = 1
                      ELSE
                          A(J1,I1) = Z1
                      END IF
                      RMAX   = AMAX1(RMAX,A(J1,I1))
                  END IF
              END IF
          END IF
   30 CONTINUE

      IF (OVER.NE.0) THEN
         STRING = 'Attempted to overwrite pixels. Change START or STEP'
         CALL STTPUT (STRING,ISTAT)
         WRITE (STRING,500) RUPP
500      FORMAT ('Overwritten value is : ',E15.6)
         CALL STTPUT (STRING,ISTAT)
      ENDIF

      RETURN
      END

      SUBROUTINE COUNT1
     . (TID,ICOL,ICOLR,N,NPIX1,A,START,STEP,RFLAG,RMIN,RMAX)
C
C compute distribution
C
      IMPLICIT NONE
      INTEGER  TID
      INTEGER  ICOL
      INTEGER  ICOLR
      INTEGER  N
      INTEGER  NPIX1
      REAL     A(NPIX1)
      REAL     START
      REAL     STEP
      REAL     RFLAG(2)
      REAL     RMIN
      REAL     RMAX
C
      LOGICAL  ISEL, INULL, INULLR
      INTEGER  I, J, STATUS, J1
      REAL     RLOW, RUPP, X1, Z1
      CHARACTER*80 STRING
      INTEGER      ISTAT,OVER
C
C
C ... clear counts
C
      RLOW   = RFLAG(1)
      RUPP   = RFLAG(2)
      DO 10 J = 1,NPIX1
          A(J)   = RLOW
   10 CONTINUE
      RMIN   = RLOW
      RMAX   = RLOW
C
C ... iteration on rows
C
      OVER = 0		! Set overwriting flag to 0
      DO 20 I = 1,N
          CALL TBSGET(TID,I,ISEL,STATUS)
          IF (ISEL) THEN
              CALL TBERDR(TID,I,ICOL,X1,INULL,STATUS)
              CALL TBERDR(TID,I,ICOLR,Z1,INULLR,STATUS)
              IF (.NOT.INULL.AND..NOT.INULLR)THEN
                  J1     = (X1-START)/STEP + 1
                  IF (J1.GE.1 .AND. J1.LE.NPIX1) THEN
                      IF (A(J1).NE.RLOW) THEN
                          A(J1)  = RUPP
                      ELSE
                          A(J1)  = Z1
                      END IF
                      RMAX   = AMAX1(RMAX,A(J1))
                  END IF
              END IF
          END IF
   20 CONTINUE

      IF (OVER.NE.0) THEN
         STRING = 'Attempted to overwrite pixels. Change START or STEP'
         CALL STTPUT (STRING,ISTAT)
         WRITE (STRING,500) RUPP
500      FORMAT ('Overwritten value is : ',E15.6)
         CALL STTPUT (STRING,ISTAT)
      ENDIF

      RETURN
      END

      SUBROUTINE LSQFIT
     .(TID,NROW,NRALL,NCOL,NIND,IVAR,K,L,N1,XMIN,XMAX,YMIN,YMAX,ERROR)
C
C SOLVE LEAST SQUARE PROBLEM
C
      IMPLICIT NONE
      INTEGER  IVAR(3)
      INTEGER TID, NROW, NRALL, NCOL, NIND, K, L, N1
      REAL    XMIN, XMAX, YMIN, YMAX, ERROR
C
      DOUBLE PRECISION G,XC
C
      LOGICAL ISEL, INULL1, INULL2, INULLY
      REAL    X1, X2, Y
      INTEGER IX1, IX2, IY, LL, I, STATUS, N, K1, KMIN
      COMMON /LS/G(50,50),XC(50),N
C

      XMIN   = 1.E30
      XMAX   = -XMIN
      YMIN   = XMIN
      YMAX   = XMAX
C
C ITERATION
C
      LL     = 0
      N1     = 0
      IX1    = IVAR(1)
      IX2    = IVAR(2)
      IY     = IVAR(3)
      DO 20 I = 1,NROW
          CALL TBSGET(TID,I,ISEL,STATUS)
          IF (ISEL) THEN
             CALL TBERDR(TID,I,IX1,X1,INULL1,STATUS) 
             CALL TBERDR(TID,I,IX2,X2,INULL2,STATUS) 
             CALL TBERDR(TID,I,IY, Y ,INULLY,STATUS) 
             IF (.NOT.INULL1 .AND. .NOT.INULL2 
     .           .AND..NOT.INULLY) THEN 
                XMIN   = AMIN1(XMIN,X1)
                XMAX   = AMAX1(XMAX,X1)
                YMIN   = AMIN1(YMIN,X2)
                YMAX   = AMAX1(YMAX,X2)
                N1     = N1 + 1
                CALL SETUP(LL+1,X1,X2,Y,K,L)
                IF (LL.NE.0) THEN
                    KMIN   = MIN(LL,N+1)
                    DO 10 K1 = 1,KMIN
                        CALL HT(K1,LL+1)
   10               CONTINUE
                END IF
                LL     = MIN(LL+1,N+1)
            END IF
          END IF
   20 CONTINUE
      CALL SOLVE
      ERROR  = ABS(G(N+1,N+1))
      RETURN
      END

      SUBROUTINE SETUP(LL,X1,X2,Y,K,L)
C
C PREPARE DATA IN WORK SPACE
C
      IMPLICIT NONE
      INTEGER LL, K, L
      REAL X1, X2, Y
C
      INTEGER IP, N, L1, K1
C
      DOUBLE PRECISION G,X,DC
C
      COMMON /LS/G(50,50),X(50),N
C
      IP     = 0
      DC     = 1.D0
      DO 20 L1 = 0,L
          IP     = IP + 1
          G(LL,IP) = DC
          DO 10 K1 = 1,K
              IP     = IP + 1
              G(LL,IP) = G(LL,IP-1)*X1
   10     CONTINUE
          DC     = DC*X2
   20 CONTINUE
      G(LL,N+1) = Y
      RETURN
      END

      SUBROUTINE HT(K,L1)

C     FUNCTION : Apply a Householder transformation to the
C                two rows K and L1 of G.
C
C     INPUT :    * matrix G(1:L1,1:n+1)
C                * numbers L1 and K of the two rows
C                * n : rank of the matrix
C
C     OUTPUT :   modified rows K and L1 of the matrix G.
      IMPLICIT NONE
      INTEGER K, L1
C
      INTEGER KH, J, N
      DOUBLE PRECISION G,X,S,H,B
      COMMON /LS/G(50,50),X(50),N

C     Construct the transformation

      S      = DSQRT(G(K,K)**2+G(L1,K)**2)
      IF (G(K,K).GE.0) THEN
          S      = -S
      END IF
      H      = G(K,K) - S
      G(K,K) = S
      B      = G(K,K)*H
C
C     Apply the transformation to the rows K and L1 of G.
C
      KH     = N + 1 - K
      IF ((KH.NE.0) .AND. (B.NE.0)) THEN
          DO 10 J = 1,KH
              S      = (G(K,K+J)*H+G(L1,K+J)*G(L1,K))/B
              G(K,K+J) = G(K,K+J) + S*H
              G(L1,K+J) = G(L1,K+J) + S*G(L1,K)
   10     CONTINUE
      END IF
      END

      SUBROUTINE SOLVE
C     FUNCTION : Solve the set of linear equations which are
C                represented by the upper triangular part
C                of the matrix G(1:n,1:n+1).
C
C     INPUT :    * matrix G(1:n,1:n+1)
C                * n : rank of the matrix
C
C     OUTPUT :   solution vector x(1:n)
      IMPLICIT NONE
      INTEGER  I, II, N, J
      DOUBLE PRECISION G,X,S
      COMMON /LS/G(50,50),X(50),N

      DO 10 I = 1,N
          X(I)   = G(I,N+1)
   10 CONTINUE
      DO 30 II = 1,N
          I      = N - II + 1
          S      = 0.D0
          IF (I.NE.N) THEN
              DO 20 J = 2,N + 1 - I
                  S      = S + G(I,I+J-1)*X(I+J-1)
   20         CONTINUE
          END IF
          X(I)   = (X(I)-S)/G(I,I)
   30 CONTINUE
      END

      SUBROUTINE COMFT2(NPIX1,NPIX2,OUT,START,STEP,K,L,RMIN,RMAX)
C
      IMPLICIT NONE
C
      INTEGER NPIX1,NPIX2,K,L,IY,IX,L1,N,K1,IP
      REAL OUT(NPIX1,NPIX2),RMIN,RMAX
      REAL START(*),STEP(*)
      DOUBLE PRECISION WORK(50),DC,YVAL,G,A
      DOUBLE PRECISION X1,X2,XSTART,YSTART,XSTEP,YSTEP
      COMMON /LS/G(50,50),A(50),N
C
C  initialize variables
      XSTART = START(1)
      YSTART = START(2)
      XSTEP  = STEP(1)
      YSTEP  = STEP(2)
C
C  calculate 1st value already for RMIN, RMAX
      IP = 0
      DC = 1.D0
      YVAL = 0.D0 
      DO 1020, L1 = 0,L
         IP = IP + 1
         WORK(IP) = DC
         YVAL = YVAL + WORK(IP)*A(IP)
         DO 1010, K1 = 1,K
            IP = IP + 1
            WORK(IP) = WORK(IP-1)*XSTART
            YVAL = YVAL + WORK(IP)*A(IP)
1010     CONTINUE
         DC = DC*YSTART
1020  CONTINUE
      RMIN = YVAL
      RMAX = RMIN
C 
C  start loop
      DO 40, IY=1,NPIX2
          X2     = (IY-1)*YSTEP + YSTART
          DO 30, IX = 1,NPIX1
              X1     = (IX-1)*XSTEP + XSTART
              IP     = 0
              DC     = 1.D0
              YVAL   = 0.D0
              DO 20, L1 = 0,L
                  IP     = IP + 1
                  WORK(IP) = DC
                  YVAL   = YVAL + WORK(IP)*A(IP)
                  DO 10, K1 = 1,K
                      IP     = IP + 1
                      WORK(IP) = WORK(IP-1)*X1
                      YVAL   = YVAL + WORK(IP)*A(IP)
   10             CONTINUE
                  DC     = DC*X2
   20         CONTINUE
              OUT(IX,IY) = YVAL
              RMIN   = AMIN1(RMIN,OUT(IX,IY))
              RMAX   = AMAX1(RMAX,OUT(IX,IY))
   30     CONTINUE
   40 CONTINUE
      RETURN
C 
      END