C @(#)modcol.for	13.1.1.2 (ESO-DMD) 02/12/99 19:14:48
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
C.IDENTIFICATION
C  program MODCOL		version 1.00	881028
C  K. Banse			ESO - Garching
C  M.Peron                      910417  modif in the call of TCLSER
C  S. Wolf                      980113  not selected table-rows has not
C                                       to be processed (swolf@eso.org)!
C.KEYWORDS
C  bulk data frame, bad pixels
C
C.PURPOSE
C
C  to replace a column or row of undefined or bad pixels values by quadratic
C  interpolation of adjacent columns or rows.The bad columns or rows are
C  defined via the cursor or coordinates or taken from a table.
C
C.ALGORITHM
C
C  The adjacent pixels are used to compute the coefficients of a least-squares
C  approximating second order polynomial
C
C.INPUT/OUTPUT
C  the following keywords are used:
C  P1/C/1/60			input frame, table
C				or
C				CURSOR
C				or
C				input frame
C  IN_A/C/1/60			input frame, if P1 does not indicate table use
C  OUT_A/C/1/60			result frame
C  ACTION/C/1/2			C or R, for column or row
C  				V or C for variable or constant offset
C				of column/row in question
C  P4/C/1/80			column or row world coords. if P1 holds just
C				input frame
C
C.VERSIONS
C
C  1.00		from program MODICOL version 2.20 as of 850212
C  		use FORTRAN 77 + new ST interfaces
C
C-------------------------------------------------------------------
C
      PROGRAM MODCOL
C 
      IMPLICIT NONE
C
      INTEGER      N,NCOUNT,IAV,ISTAT
      INTEGER*8    PNTRA,PNTRB,PNTRX,PNTRY
      INTEGER      IMNOA,IMNOB,STAT
      INTEGER      INDXFL,SIZE,NAXIS,SLEN,SLEN1,TBLFLG
      INTEGER      NCOL,NROW
      INTEGER      IBUF(2),KINDEX,OFF,KROW,KROW1,KROW2
      INTEGER      TIDI,TABCLN(1)
      INTEGER      UNI(1),NULO,MADRID(1)
      INTEGER      NPIXA(2)
C
      CHARACTER    CUNITA*48,IDENTA*72,ACTION*2,COLROW*80
      CHARACTER    FRAMEA*80,FRAMEB*80,TABLE*80
      CHARACTER    OUTPUT*80,COLNAM*20
      CHARACTER    MESS(2)*50,ERRMS1*40,ERRMS2*40
      CHARACTER    SUBPAR(40)*30
C
      REAL         TEMP,CUTS(4),FMIN,FMAX,RMIN,RMAX,RR
C
      DOUBLE PRECISION STARTA(2),STEPA(2),DD
C
      LOGICAL      TABNUL(1),SELFLG
C
      COMMON   /VMR/  MADRID
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
C
      DATA      FRAMEA   /' '/,    FRAMEB   /' '/
      DATA      IDENTA   /' '/,    CUNITA    /' '/
      DATA      OUTPUT   /' '/,    TABLE    /' '/
      DATA      MESS     /'column/row too close to image boundary... ',
     +                   'column/row too large...'/
      DATA      ERRMS1   /'col/row no.s missing...'/
      DATA      ERRMS2   /'invalid syntax for col/row no.s ...'/
C
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
C  get into MIDAS environment
      CALL STSPRO('MODCOL')
C
C  get input option (frame or frame,table) + read keys
      CALL STKRDC('IN_A',1,1,80,IAV,OUTPUT,UNI,NULO,STAT)
      CALL STKRDC('ACTION',1,1,2,IAV,ACTION,UNI,NULO,STAT)
      CALL UPCAS(ACTION,ACTION)
      IF (ACTION(1:1).EQ.'C') THEN
         INDXFL = 1                  !index = 1 for columns...
      ELSE
         INDXFL = 2                  !index = 2 for rows
      ENDIF
C
      N = INDEX(OUTPUT,',')
      IF (N.GT.0) THEN
         TBLFLG = 2
         FRAMEA(1:) = OUTPUT(1:N-1)
         CALL CLNFRA(FRAMEA,FRAMEA,0)             !because it's not isolated
         TABLE(1:) = OUTPUT(N+1:)
         N = INDEX(TABLE,',')
         IF (N.GT.0) THEN
            IAV = N
            IF (TABLE(N+1:N+1).EQ.':') IAV = N + 1
            COLNAM(1:) = TABLE(IAV+1:)//' '
            TABLE(N:) = ' '
         ELSE
            IF (ACTION(1:1).EQ.'C') THEN         !default column labels
               COLNAM(1:) = 'X '
            ELSE
               COLNAM(1:) = 'Y '
            ENDIF
         ENDIF
         CALL TBTOPN(TABLE,F_I_MODE,TIDI,STAT)
         CALL TBIGET(TIDI,NCOL,NROW,N,N,N,STAT)
         CALL TBLSER(TIDI,COLNAM,TABCLN(1),STAT)      !find columns...
         IF (TABCLN(1).LE.0) THEN
            N = INDEX(COLNAM,' ')
            OUTPUT(1:) = 'column '//COLNAM(1:N)//
     +                   ' not found in input table...'
            CALL STETER (9,OUTPUT)
         ENDIF
      ELSE
         TBLFLG = 1
         FRAMEA(1:) = OUTPUT(1:)
         CALL STKRDC('P4',1,1,80,IAV,COLROW,UNI,NULO,STAT)      !get coords.
      ENDIF
      CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEB,UNI,NULO,STAT)
C
C  map input + result frame
      CALL GENEQF(FRAMEA,FRAMEB,STAT)
      IF (STAT.EQ.1) THEN                !same name
         CALL STIGET(FRAMEA,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,
     +               2,NAXIS,NPIXA,STARTA,STEPA,IDENTA,
     +               CUNITA,PNTRA,IMNOA,STAT)
         IF (NAXIS.LT.2) THEN
            CALL STETER(1,'input frame must be a 2-dim frame...')
         ELSE IF (NAXIS.GT.2) THEN
            CALL STTPUT
     +      ('only first plane of 3-dim input frame processed',STAT)
            NAXIS = 2
         ENDIF
         PNTRB = PNTRA
         IMNOB = IMNOA
         SIZE = NPIXA(1)*NPIXA(2)
      ELSE
         CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +               2,NAXIS,NPIXA,STARTA,STEPA,IDENTA,
     +               CUNITA,PNTRA,IMNOA,STAT)
         IF (NAXIS.LT.2) THEN
            CALL STETER(1,'input frame must be a 2-dim frame...')
         ELSE IF (NAXIS.GT.2) THEN
            CALL STTPUT
     +      ('only first plane of 3-dim input frame processed',STAT)
            NAXIS = 2
         ENDIF
         CALL STIPUT(FRAMEB,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +               NAXIS,NPIXA,STARTA,STEPA,IDENTA,
     +               CUNITA,PNTRB,IMNOB,STAT)
         SIZE = NPIXA(1)*NPIXA(2)
         CALL COPYF(MADRID(PNTRA),MADRID(PNTRB),SIZE)            !copy data...
      ENDIF
      CALL FRAMOU(FRAMEA)
C
      IF (ACTION(1:1).EQ.'C') THEN                !allocate workspace
         SIZE = NPIXA(2)
      ELSE
         SIZE = NPIXA(1)
      ENDIF
      CALL STFXMP(SIZE,D_R4_FORMAT,PNTRX,STAT)
      CALL STFXMP(SIZE,D_R4_FORMAT,PNTRY,STAT)
C 
C  save min,max + cuts of image
      CALL STDRDR(IMNOA,'LHCUTS',1,4,IAV,CUTS,UNI,NULO,STAT)
      RMIN = CUTS(3)
      RMAX = CUTS(4)
C
C  fork according to input mode
      IF (TBLFLG.EQ.2) THEN
         KROW2 = 0
         ISTAT = 2                           !thus we start with KROW = 0
         GOTO 3000                           !and ISTAT != 1
      ENDIF
C
C  extract columns/rows from COLROW
      ISTAT = 1  
      KINDEX = 0
      NCOUNT = 0
      OFF = 1
      DO 600, N=1,40
         CALL EXTRSS(COLROW,',',OFF,SUBPAR(N),SLEN)
         IF (SLEN.LE.0) THEN
            IF (NCOUNT.LE.0) CALL STETER(2,ERRMS1)
            GOTO 1000                                   !we reached the end...
         ELSE
            NCOUNT = N
         ENDIF
600   CONTINUE
C
C  finally loop through columns or rows
C
1000  KINDEX = KINDEX + ISTAT
      IF (KINDEX.GT.NCOUNT) GOTO 4000
C 
      N = KINDEX
      IF (SUBPAR(N)(1:1).EQ.'@') THEN               !input = @pixel
         CALL GENCNV(SUBPAR(N)(2:),1,1,IBUF,RR,DD,SLEN)
         IF (SLEN.LE.0) CALL STETER(3,ERRMS2)
      ELSE                                          !input = world coord
         CALL GENCNV(SUBPAR(N),2,1,IAV,TEMP,DD,SLEN)
         IF (SLEN.LE.0) CALL STETER(3,ERRMS2)
         IBUF(1) = NINT((TEMP - STARTA(INDXFL))/STEPA(INDXFL)) + 1
      ENDIF
      IF (KINDEX.EQ.NCOUNT) THEN
         IBUF(2) = 0
      ELSE
         N = N + 1
         IF (SUBPAR(N)(1:1).EQ.'@') THEN               !input = @pixel
            CALL GENCNV(SUBPAR(N)(2:),1,1,IBUF(2),TEMP,DD,SLEN)
            IF (SLEN.LE.0) CALL STETER(3,ERRMS2)
         ELSE                                          !input = world coord
            CALL GENCNV(SUBPAR(N),2,1,IAV,TEMP,DD,SLEN)
            IF (SLEN.LE.0) CALL STETER(3,ERRMS2)
            IBUF(2) = NINT((TEMP - STARTA(INDXFL))/STEPA(INDXFL)) + 1
         ENDIF
      ENDIF
C 
1010  FMIN = RMIN
      FMAX = RMAX
      CALL DOIT(ACTION,MADRID(PNTRB),MADRID(PNTRX),MADRID(PNTRY),
     +          STARTA,STEPA,NPIXA,
     +          IBUF,FMIN,FMAX,ISTAT)

      IF (ISTAT.LT.0) THEN
         CALL STTPUT(MESS(-ISTAT),STAT)
         ISTAT = 1            !on return, ISTAT = 1 or 2
      ELSE                    !depending on no. of columns/rows used
         SLEN = INDEX(SUBPAR(KINDEX),' ')            !cut off trailing blanks
         IF (SLEN.LE.0) SLEN = LEN(SUBPAR(KINDEX))
C
         IF (ISTAT.EQ.1) THEN
            OUTPUT(1:) = 'column/row at x/y = '
     +                   //SUBPAR(KINDEX)(1:SLEN)//
     +                   'cleaned... '
         ELSE
            SLEN1 = INDEX(SUBPAR(KINDEX+1),' ')      !cut off trailing blanks
            IF (SLEN1.LE.0) SLEN1 = LEN(SUBPAR(KINDEX+1))
            OUTPUT(1:) = 'columns/rows at x/y = '//
     +                   SUBPAR(KINDEX)(1:SLEN)//', '//
     +                   SUBPAR(KINDEX+1)(1:SLEN1)//'cleaned...'
         ENDIF
         CALL STTPUT(OUTPUT,STAT)
         RMIN = MIN(RMIN,FMIN)
         RMAX = MAX(RMAX,FMAX)
      ENDIF
      GOTO (1000,3000),TBLFLG                  !look for more...
C
C  here we get our input data from a table
3000  KROW = KROW2 + ISTAT - 1                 !KROW -> KROW2 or -> KROW2+1
      IF (KROW.LE.NROW) THEN
         KROW1 = KROW
         IF (ISTAT.EQ.1) THEN                  !only one row used last time...
            IBUF(1) = IBUF(2)
            SUBPAR(1)(1:) = SUBPAR(2)(1:)
            GOTO 3100                          !so reuse last value
         ENDIF
C 
3030     CALL TBSGET(TIDI,KROW1,SELFLG,STAT)       !only use selected rows
         IF (.NOT.SELFLG) THEN
            KROW1 = KROW1 + 1                      !get next row
            IF (KROW1.LE.NROW) THEN
               GOTO 3030
            ELSE
               GOTO 3300                           !end of table reached
            ENDIF
         ENDIF
C 
         CALL TBRRDR(TIDI,KROW1,1,TABCLN,TEMP,TABNUL,STAT)
         IBUF(1) = NINT((TEMP - STARTA(INDXFL))/STEPA(INDXFL)) + 1
         WRITE(SUBPAR(1),10000) TEMP
         DO 3060, N=1,30
            IF (SUBPAR(1)(N:N).NE.' ') THEN
               IF (N.GT.1) SUBPAR(1)(1:) = SUBPAR(1)(N:)//' '
               GOTO 3100
            ENDIF
3060     CONTINUE
C 
3100     KROW2 = KROW1 + 1                      !point to next possible row
         IF (KROW1.EQ.NROW) THEN   
            IBUF(2) = 0
         ELSE
3130        CALL TBSGET(TIDI,KROW2,SELFLG,STAT)       !only use selected rows
            IF (.NOT.SELFLG) THEN
               KROW2 = KROW2 + 1                      !get next row
               IF (KROW2.LE.NROW) THEN
                  GOTO 3130
               ELSE
                  IBUF(2) = 0
                  GOTO 3200 
               ENDIF
            ENDIF
C 
            CALL TBRRDR(TIDI,KROW2,1,TABCLN,TEMP,TABNUL,STAT)
            IBUF(2) = NINT((TEMP - STARTA(INDXFL))/STEPA(INDXFL)) + 1
            WRITE(SUBPAR(2),10000) TEMP
            DO 3160, N=1,30
               IF (SUBPAR(2)(N:N).NE.' ') THEN
                  IF (N.GT.1) SUBPAR(2)(1:) = SUBPAR(2)(N:)//' '
                  GOTO 3200
               ENDIF
3160        CONTINUE
         ENDIF
C 
3200     KINDEX = 1                              !force index to 1
         GOTO 1010                               !now continue as above
      ENDIF
C 
C  close table
3300  CALL TBTCLO(TIDI,STAT)
C
C  get new dynamic range + write descriptor LHCUTS
4000  CUTS(3) = MIN(RMIN,CUTS(3))
      CUTS(4) = MAX(RMAX,CUTS(4))
C
C  copy descriptors + append HISTORY
      CALL DSCUPT(IMNOA,IMNOB,' ',STAT)
      CALL STDWRR(IMNOB,'LHCUTS',CUTS,1,4,UNI,STAT)      !use old cuts...
C
C  That's it folks...
      CALL STSEPI
C
C  format statements
10000 FORMAT(F20.10)
      END

      SUBROUTINE DOIT(ACTION,A,TEMPA,TEMPB,
     +                START,STEP,NPIX,COORD,FMIN,FMAX,ISTAT)
C
      IMPLICIT NONE
C
      INTEGER      COORD(2),ISTAT
      INTEGER      JOF,JOF1,N,M,NFLAG,KFLAG
      INTEGER      NEXT,NP,JNDX,NOLINS
      INTEGER      NPIX(2),INDX(4),IOF(4)
C
      CHARACTER    ACTION*2
C
      REAL         A(*),TEMPA(*),TEMPB(*),FMIN,FMAX
      REAL         V,VV,W,WW,XX,XXX
      REAL         RMIN,RMAX
C
      DOUBLE PRECISION START(2),STEP(2)
      DOUBLE PRECISION   S(4),Z(4),X(4),F(4)
      DOUBLE PRECISION   RR1,RR2,RR3,P1,P2,P3,Q
      DOUBLE PRECISION   SUM1,SUM2
      DOUBLE PRECISION   A0,A1,A2
C
      IF (ACTION(1:1).EQ.'C') THEN
         NFLAG = 1                              !columns...
      ELSE
         NFLAG = 2                              !rows...
      ENDIF
      KFLAG = 3 - NFLAG 
      INDX(1) = COORD(1) - 2                    !take 2 pixels less
      NEXT = COORD(2) - 2
C
      INDX(2) = INDX(1) + 1
      JNDX = INDX(2) + 1                        !keep index of bad line...
      IF (NEXT .EQ. INDX(1)+1) THEN
         NOLINS = 2                             !double lines...
      ELSE
         NOLINS = 1
      ENDIF
      INDX(3) = INDX(2) + NOLINS + 1                  !single line...
      INDX(4) = INDX(3) + 1
C 
C  check input pars
      IF ((INDX(1).LT.1).OR.(INDX(4).GT.NPIX(NFLAG))) THEN
         ISTAT = -1
         RETURN
      ELSE
         ISTAT = NOLINS
      ENDIF
C
C  get "x" values + related sums for later calculations
      DO 200, N=1,4                              !get x values...
         X(N) = START(NFLAG) + (INDX(N)-1)*STEP(NFLAG)
         Z(N) = 1.D0
200   CONTINUE
      XX = START(NFLAG) + (JNDX-1)*STEP(NFLAG)      !x of bad row/column...
      XXX = XX + STEP(NFLAG)
C
      DO 1000, M=1,4
         DO 400, N=1,4                               !get x**2,x**3,x**4 ...
            Z(N) = Z(N)*X(N)
400      CONTINUE
         S(M) = 0.D0                                !compute sums
         DO 600, N=1,4
            S(M) = S(M) + Z(N)
600      CONTINUE
1000  CONTINUE
C
C  branch according to columns + rows
      IF (ACTION(1:1).EQ.'C') GOTO 5000
C
C  handle rows...
      DO 1200, N=1,4                                 !get y values...
         IOF(N) = (INDX(N)-1)*NPIX(1)
1200  CONTINUE
      JOF = (JNDX-1)*NPIX(1)
      JOF1 = JOF + NPIX(1)
C
C  go through the row
      DO 3300, NP=1,NPIX(1)
         DO 1400, M=1,4                              !get function values
            F(M) = A(IOF(M)+NP)
            Z(M) = 0.D0
1400     CONTINUE
C
         DO 1600, M=1,4                              !compute sums on right side
            Z(1) = Z(1) + F(M)
            Z(2) = Z(2) + F(M)*X(M)
            Z(3) = Z(3) + F(M)*X(M)*X(M)
1600     CONTINUE
C
C  reduce + solve the linear system  4.*A0 + S1*A1 + S2*A2 = Z1
C                                    S1*A0 + S2*A1 + S3*A2 = Z2
C                                    S2*A0 + S3*A1 + S4*A2 = Z3
         Q = S(1) * 0.25D0
         RR1 = S(2) - S(1)*Q
         RR2 = S(3) - S(2)*Q
         RR3 = Z(2) - Z(1)*Q
         Q = S(2) * 0.25D0
         P1 = S(3) - S(1)*Q
         P2 = S(4) - S(2)*Q
         P3 = Z(3) - Z(1)*Q
C                   we now have:          RR1*A1 + RR2*A2 = RR3
C                                         P1*A1 + P2*A2 = P3
         Q = P1/RR1
         IF (ABS(P2-RR2*Q).LE.10.D-30) THEN
            IF (ABS(P1 - RR1*Q).LE.10.D-30) THEN      !linearly dependant system
               A0 = 0.D0
               A1 = 1.D0
               A2 = 0.D0
               GOTO 3000
            ENDIF
            Q = P2/RR2
            A1 = (P3 - RR3*Q) / (P1 - RR1*Q)
            A2 = (RR3 - A1*RR1)/RR2
         ELSE
            A2 = (P3 - RR3*Q) / (P2 - RR2*Q)
            A1 = (RR3 - A2*RR2)/RR1
         ENDIF
C
         A0 = (Z(1) - S(2)*A2 - S(1)*A1) * .25D0
C
C  finally calculate interpolated value...!
3000     TEMPA(NP) = A0 + A1*XX + A2*XX*XX
         IF (NOLINS.EQ.2)
     +      TEMPB(NP) = A0 + A1*XXX + A2*XXX*XXX
3300  CONTINUE
C
C  check, if we work on single line or not
C
      IF (NOLINS.EQ.2) GOTO 4000
C
C  for constant offset, get average difference to orignal bad row
      IF (ACTION(2:2).EQ.'C') THEN
         SUM1 = 0.D0
         DO 3500, NP=1,NPIX(1)
            SUM1 = SUM1 + DBLE(A(JOF+NP) - TEMPA(NP))  !sum differences
3500     CONTINUE
         V = SNGL(SUM1/NPIX(1))            !get average offset
C
C  and correct pixels in bad row with that offset
         DO 3600, NP=1,NPIX(1)
            W = A(JOF+NP)
            A(JOF+NP) = W - V
            IF (W.LT.FMIN) THEN
               FMIN = W
            ELSE
               IF (W.GT.FMAX) FMAX = W
            ENDIF
3600     CONTINUE
C
      ELSE
C
C  store approximated pixels in bad row
         DO 3800, NP=1,NPIX(1)
            W = TEMPA(NP)
            A(JOF+NP) = W
            IF (W.LT.FMIN) THEN
               FMIN = W
            ELSE
               IF (W.GT.FMAX) FMAX = W
            ENDIF
3800     CONTINUE
      ENDIF
      RETURN
C
C  here to handle two adjacent "bad" rows
C
4000  IF (ACTION(2:2).EQ.'C') THEN
         SUM1 = 0.D0                  !get average difference to adjacent lines
         SUM2 = 0.D0
         DO 4200, NP=1,NPIX(1)
            SUM1 = SUM1 + DBLE(A(JOF+NP) - TEMPA(NP))   !sum differences
            SUM2 = SUM2 + DBLE(A(JOF1+NP) - TEMPB(NP))  !sum differences
4200     CONTINUE
         V = SNGL(SUM1/NPIX(1))                         !get average offset
         VV = SNGL(SUM2/NPIX(1))
C
C  and correct pixels in bad rows with these offsets
         DO 4400, NP=1,NPIX(1)
            W = A(JOF+NP)
            A(JOF+NP) = W - V
            WW = A(JOF1+NP)
            A(JOF1+NP) = WW - VV
            RMIN = MIN(W,WW)
            RMAX = MAX(W,WW)
            IF (RMIN.LT.FMIN) FMIN = RMIN
            IF (RMAX.GT.FMAX) FMAX = RMAX
4400     CONTINUE
C
      ELSE
C
C  store approximated pixels in bad rows
         DO 4600, NP=1,NPIX(1)
            W = TEMPA(NP)
            WW = TEMPB(NP)
            A(JOF+NP) = W
            A(JOF1+NP) = WW
            RMIN = MIN(W,WW)
            RMAX = MAX(W,WW)
            IF (RMIN.LT.FMIN) FMIN = RMIN
            IF (RMAX.GT.FMAX) FMAX = RMAX
4600     CONTINUE
      ENDIF
      RETURN
C
C  handle columns...
C
C  go through the column
5000  JOF = 0
      DO 6600, NP=1,NPIX(2)
         DO 5200, M=1,4                        !get function values
            F(M) = A(JOF+INDX(M))
            Z(M) = 0.D0
5200     CONTINUE
C
         DO 5400, M=1,4                        !compute sums on right side
            Z(1) = Z(1) + F(M)
            Z(2) = Z(2) + F(M)*X(M)
            Z(3) = Z(3) + F(M)*X(M)*X(M)
5400     CONTINUE
C
C  reduce + solve the linear system     4.*A0 + S1*A1 + S2*A2 = Z1
C                                       S1*A0 + S2*A1 + S3*A2 = Z2
C                                       S2*A0 + S3*A1 + S4*A2 = Z3
         Q = S(1)*.25D0
         RR1 = S(2) - S(1)*Q
         RR2 = S(3) - S(2)*Q
         RR3 = Z(2) - Z(1)*Q
         Q = S(2)*.25D0
         P1 = S(3) - S(1)*Q
         P2 = S(4) - S(2)*Q
         P3 = Z(3) - Z(1)*Q
C  we now have:                         RR1*A1 + RR2*A2 = RR3
C                                       P1*A1 + P2*A2 = P3
         Q = P1/RR1
         IF (ABS(P2-RR2*Q).LE.10.D-30) THEN
            IF (ABS(P1 - RR1*Q).LE.10.D-30) THEN      !linearly dependant system
               A0 = 0.D0
               A1 = 1.D0
               A2 = 0.D0
               GOTO 6000
            ENDIF
            Q = P2/RR2
            A1 = (P3 - RR3*Q) / (P1 - RR1*Q)
            A2 = (RR3 - A1*RR1)/RR2
         ELSE
            A2 = (P3 - RR3*Q) / (P2 - RR2*Q)
            A1 = (RR3 - A2*RR2)/RR1
         ENDIF
C
         A0 = (Z(1) - S(2)*A2 - S(1)*A1)*.25D0
C
C  finally calculate interpolated value...!
6000     TEMPA(NP) = A0 + A1*XX + A2*XX*XX
         IF (NOLINS.EQ.2)
     +      TEMPB(NP) = A0 + A1*XXX + A2*XXX*XXX
         JOF = JOF + NPIX(1)                  !update offset to next row...
6600  CONTINUE
C
C  check, if we work on single line or not
C
      IF (NOLINS.EQ.2) GOTO 8000
C
C  for constant offset, get average difference to original bad column
      IF (ACTION(2:2).EQ.'C') THEN
         SUM1 = 0.D0
         JOF = JNDX
         DO 6800, NP=1,NPIX(2)
            SUM1 = SUM1 + DBLE(A(JOF) - TEMPA(NP))      !sum differences
            JOF = JOF + NPIX(1)
6800     CONTINUE
         W = SNGL(SUM1/NPIX(2))                         !get average offset
C
C  and correct pixels in bad column with that offset
         JOF = JNDX
         DO 7000, NP=1,NPIX(2)
            V = A(JOF) - W
            A(JOF) = V 
            IF (V.LT.FMIN) THEN
               FMIN = V
            ELSE
               IF (V.GT.FMAX) FMAX = V
            ENDIF
            JOF = JOF + NPIX(1)
7000     CONTINUE
C
      ELSE
C
C  fill approximated pixels into frame
         JOF = JNDX
         DO 7400, NP=1,NPIX(2)
            V = TEMPA(NP)
            A(JOF) = V
            IF (V.LT.FMIN) THEN
               FMIN = V
            ELSE
               IF (V.GT.FMAX) FMAX = V
            ENDIF
            JOF = JOF + NPIX(1)
7400     CONTINUE
      ENDIF
      RETURN
C
C  here to handle two adjacent "bad" columns
8000  IF (ACTION(2:2).EQ.'C') THEN
         SUM1 = 0.D0                  !get average difference to adjacent lines
         SUM2 = 0.D0
         JOF = JNDX
C
         DO 8200, NP=1,NPIX(2)
            SUM1 = SUM1 + DBLE(A(JOF) - TEMPA(NP))        !sum differences
            SUM2 = SUM2 + DBLE(A(JOF+1) - TEMPB(NP))
            JOF = JOF + NPIX(1)
8200     CONTINUE
         V = SNGL(SUM1/NPIX(2))                           !get average offset
         VV = SNGL(SUM2/NPIX(2))
C
C  and correct pixels in bad columns with these offsets
         JOF = JNDX
         DO 8400, NP=1,NPIX(2)
            W = A(JOF)
            WW = A(JOF+1)
            A(JOF) = W - V
            A(JOF+1) = WW - VV
            RMIN = MIN(W,WW)
            RMAX = MAX(W,WW)
            IF (RMIN.LT.FMIN) FMIN = RMIN
            IF (RMAX.GT.FMAX) FMAX = RMAX
            JOF = JOF + NPIX(1)
8400     CONTINUE
C
      ELSE
C
C  fill approximated pixels into frame
         JOF = JNDX
         DO 8600, NP=1,NPIX(2)
            W = TEMPA(NP)
            WW = TEMPB(NP)
            A(JOF) = W
            A(JOF+1) = WW
            W = A(JOF)
            RMIN = MIN(W,WW)
            RMAX = MAX(W,WW)
            IF (RMIN.LT.FMIN) FMIN = RMIN
            IF (RMAX.GT.FMAX) FMAX = RMAX
            JOF = JOF + NPIX(1)
8600     CONTINUE
      ENDIF
      RETURN
C
      END