← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to prim/general/src/modcol.for

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
prim/general/src/modcol.for
 
1
C @(#)modcol.for        13.1.1.2 (ESO-DMD) 02/12/99 19:14:48
 
2
C===========================================================================
 
3
C Copyright (C) 1995 European Southern Observatory (ESO)
 
4
C
 
5
C This program is free software; you can redistribute it and/or 
 
6
C modify it under the terms of the GNU General Public License as 
 
7
C published by the Free Software Foundation; either version 2 of 
 
8
C the License, or (at your option) any later version.
 
9
C
 
10
C This program is distributed in the hope that it will be useful,
 
11
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
C GNU General Public License for more details.
 
14
C
 
15
C You should have received a copy of the GNU General Public 
 
16
C License along with this program; if not, write to the Free 
 
17
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Correspondence concerning ESO-MIDAS should be addressed as follows:
 
21
C       Internet e-mail: midas@eso.org
 
22
C       Postal address: European Southern Observatory
 
23
C                       Data Management Division 
 
24
C                       Karl-Schwarzschild-Strasse 2
 
25
C                       D 85748 Garching bei Muenchen 
 
26
C                       GERMANY
 
27
C===========================================================================
 
28
C
 
29
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
30
C
 
31
C.IDENTIFICATION
 
32
C  program MODCOL               version 1.00    881028
 
33
C  K. Banse                     ESO - Garching
 
34
C  M.Peron                      910417  modif in the call of TCLSER
 
35
C  S. Wolf                      980113  not selected table-rows has not
 
36
C                                       to be processed (swolf@eso.org)!
 
37
C.KEYWORDS
 
38
C  bulk data frame, bad pixels
 
39
C
 
40
C.PURPOSE
 
41
C
 
42
C  to replace a column or row of undefined or bad pixels values by quadratic
 
43
C  interpolation of adjacent columns or rows.The bad columns or rows are
 
44
C  defined via the cursor or coordinates or taken from a table.
 
45
C
 
46
C.ALGORITHM
 
47
C
 
48
C  The adjacent pixels are used to compute the coefficients of a least-squares
 
49
C  approximating second order polynomial
 
50
C
 
51
C.INPUT/OUTPUT
 
52
C  the following keywords are used:
 
53
C  P1/C/1/60                    input frame, table
 
54
C                               or
 
55
C                               CURSOR
 
56
C                               or
 
57
C                               input frame
 
58
C  IN_A/C/1/60                  input frame, if P1 does not indicate table use
 
59
C  OUT_A/C/1/60                 result frame
 
60
C  ACTION/C/1/2                 C or R, for column or row
 
61
C                               V or C for variable or constant offset
 
62
C                               of column/row in question
 
63
C  P4/C/1/80                    column or row world coords. if P1 holds just
 
64
C                               input frame
 
65
C
 
66
C.VERSIONS
 
67
C
 
68
C  1.00         from program MODICOL version 2.20 as of 850212
 
69
C               use FORTRAN 77 + new ST interfaces
 
70
C
 
71
C-------------------------------------------------------------------
 
72
C
 
73
      PROGRAM MODCOL
 
74
C 
 
75
      IMPLICIT NONE
 
76
C
 
77
      INTEGER      N,NCOUNT,IAV,ISTAT
 
78
      INTEGER*8    PNTRA,PNTRB,PNTRX,PNTRY
 
79
      INTEGER      IMNOA,IMNOB,STAT
 
80
      INTEGER      INDXFL,SIZE,NAXIS,SLEN,SLEN1,TBLFLG
 
81
      INTEGER      NCOL,NROW
 
82
      INTEGER      IBUF(2),KINDEX,OFF,KROW,KROW1,KROW2
 
83
      INTEGER      TIDI,TABCLN(1)
 
84
      INTEGER      UNI(1),NULO,MADRID(1)
 
85
      INTEGER      NPIXA(2)
 
86
C
 
87
      CHARACTER    CUNITA*48,IDENTA*72,ACTION*2,COLROW*80
 
88
      CHARACTER    FRAMEA*80,FRAMEB*80,TABLE*80
 
89
      CHARACTER    OUTPUT*80,COLNAM*20
 
90
      CHARACTER    MESS(2)*50,ERRMS1*40,ERRMS2*40
 
91
      CHARACTER    SUBPAR(40)*30
 
92
C
 
93
      REAL         TEMP,CUTS(4),FMIN,FMAX,RMIN,RMAX,RR
 
94
C
 
95
      DOUBLE PRECISION STARTA(2),STEPA(2),DD
 
96
C
 
97
      LOGICAL      TABNUL(1),SELFLG
 
98
C
 
99
      COMMON   /VMR/  MADRID
 
100
C
 
101
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
102
C
 
103
      DATA      FRAMEA   /' '/,    FRAMEB   /' '/
 
104
      DATA      IDENTA   /' '/,    CUNITA    /' '/
 
105
      DATA      OUTPUT   /' '/,    TABLE    /' '/
 
106
      DATA      MESS     /'column/row too close to image boundary... ',
 
107
     +                   'column/row too large...'/
 
108
      DATA      ERRMS1   /'col/row no.s missing...'/
 
109
      DATA      ERRMS2   /'invalid syntax for col/row no.s ...'/
 
110
C
 
111
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
112
C
 
113
C  get into MIDAS environment
 
114
      CALL STSPRO('MODCOL')
 
115
C
 
116
C  get input option (frame or frame,table) + read keys
 
117
      CALL STKRDC('IN_A',1,1,80,IAV,OUTPUT,UNI,NULO,STAT)
 
118
      CALL STKRDC('ACTION',1,1,2,IAV,ACTION,UNI,NULO,STAT)
 
119
      CALL UPCAS(ACTION,ACTION)
 
120
      IF (ACTION(1:1).EQ.'C') THEN
 
121
         INDXFL = 1                  !index = 1 for columns...
 
122
      ELSE
 
123
         INDXFL = 2                  !index = 2 for rows
 
124
      ENDIF
 
125
C
 
126
      N = INDEX(OUTPUT,',')
 
127
      IF (N.GT.0) THEN
 
128
         TBLFLG = 2
 
129
         FRAMEA(1:) = OUTPUT(1:N-1)
 
130
         CALL CLNFRA(FRAMEA,FRAMEA,0)             !because it's not isolated
 
131
         TABLE(1:) = OUTPUT(N+1:)
 
132
         N = INDEX(TABLE,',')
 
133
         IF (N.GT.0) THEN
 
134
            IAV = N
 
135
            IF (TABLE(N+1:N+1).EQ.':') IAV = N + 1
 
136
            COLNAM(1:) = TABLE(IAV+1:)//' '
 
137
            TABLE(N:) = ' '
 
138
         ELSE
 
139
            IF (ACTION(1:1).EQ.'C') THEN         !default column labels
 
140
               COLNAM(1:) = 'X '
 
141
            ELSE
 
142
               COLNAM(1:) = 'Y '
 
143
            ENDIF
 
144
         ENDIF
 
145
         CALL TBTOPN(TABLE,F_I_MODE,TIDI,STAT)
 
146
         CALL TBIGET(TIDI,NCOL,NROW,N,N,N,STAT)
 
147
         CALL TBLSER(TIDI,COLNAM,TABCLN(1),STAT)      !find columns...
 
148
         IF (TABCLN(1).LE.0) THEN
 
149
            N = INDEX(COLNAM,' ')
 
150
            OUTPUT(1:) = 'column '//COLNAM(1:N)//
 
151
     +                   ' not found in input table...'
 
152
            CALL STETER (9,OUTPUT)
 
153
         ENDIF
 
154
      ELSE
 
155
         TBLFLG = 1
 
156
         FRAMEA(1:) = OUTPUT(1:)
 
157
         CALL STKRDC('P4',1,1,80,IAV,COLROW,UNI,NULO,STAT)      !get coords.
 
158
      ENDIF
 
159
      CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEB,UNI,NULO,STAT)
 
160
C
 
161
C  map input + result frame
 
162
      CALL GENEQF(FRAMEA,FRAMEB,STAT)
 
163
      IF (STAT.EQ.1) THEN                !same name
 
164
         CALL STIGET(FRAMEA,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,
 
165
     +               2,NAXIS,NPIXA,STARTA,STEPA,IDENTA,
 
166
     +               CUNITA,PNTRA,IMNOA,STAT)
 
167
         IF (NAXIS.LT.2) THEN
 
168
            CALL STETER(1,'input frame must be a 2-dim frame...')
 
169
         ELSE IF (NAXIS.GT.2) THEN
 
170
            CALL STTPUT
 
171
     +      ('only first plane of 3-dim input frame processed',STAT)
 
172
            NAXIS = 2
 
173
         ENDIF
 
174
         PNTRB = PNTRA
 
175
         IMNOB = IMNOA
 
176
         SIZE = NPIXA(1)*NPIXA(2)
 
177
      ELSE
 
178
         CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
 
179
     +               2,NAXIS,NPIXA,STARTA,STEPA,IDENTA,
 
180
     +               CUNITA,PNTRA,IMNOA,STAT)
 
181
         IF (NAXIS.LT.2) THEN
 
182
            CALL STETER(1,'input frame must be a 2-dim frame...')
 
183
         ELSE IF (NAXIS.GT.2) THEN
 
184
            CALL STTPUT
 
185
     +      ('only first plane of 3-dim input frame processed',STAT)
 
186
            NAXIS = 2
 
187
         ENDIF
 
188
         CALL STIPUT(FRAMEB,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
189
     +               NAXIS,NPIXA,STARTA,STEPA,IDENTA,
 
190
     +               CUNITA,PNTRB,IMNOB,STAT)
 
191
         SIZE = NPIXA(1)*NPIXA(2)
 
192
         CALL COPYF(MADRID(PNTRA),MADRID(PNTRB),SIZE)            !copy data...
 
193
      ENDIF
 
194
      CALL FRAMOU(FRAMEA)
 
195
C
 
196
      IF (ACTION(1:1).EQ.'C') THEN                !allocate workspace
 
197
         SIZE = NPIXA(2)
 
198
      ELSE
 
199
         SIZE = NPIXA(1)
 
200
      ENDIF
 
201
      CALL STFXMP(SIZE,D_R4_FORMAT,PNTRX,STAT)
 
202
      CALL STFXMP(SIZE,D_R4_FORMAT,PNTRY,STAT)
 
203
C 
 
204
C  save min,max + cuts of image
 
205
      CALL STDRDR(IMNOA,'LHCUTS',1,4,IAV,CUTS,UNI,NULO,STAT)
 
206
      RMIN = CUTS(3)
 
207
      RMAX = CUTS(4)
 
208
C
 
209
C  fork according to input mode
 
210
      IF (TBLFLG.EQ.2) THEN
 
211
         KROW2 = 0
 
212
         ISTAT = 2                           !thus we start with KROW = 0
 
213
         GOTO 3000                           !and ISTAT != 1
 
214
      ENDIF
 
215
C
 
216
C  extract columns/rows from COLROW
 
217
      ISTAT = 1  
 
218
      KINDEX = 0
 
219
      NCOUNT = 0
 
220
      OFF = 1
 
221
      DO 600, N=1,40
 
222
         CALL EXTRSS(COLROW,',',OFF,SUBPAR(N),SLEN)
 
223
         IF (SLEN.LE.0) THEN
 
224
            IF (NCOUNT.LE.0) CALL STETER(2,ERRMS1)
 
225
            GOTO 1000                                   !we reached the end...
 
226
         ELSE
 
227
            NCOUNT = N
 
228
         ENDIF
 
229
600   CONTINUE
 
230
C
 
231
C  finally loop through columns or rows
 
232
C
 
233
1000  KINDEX = KINDEX + ISTAT
 
234
      IF (KINDEX.GT.NCOUNT) GOTO 4000
 
235
C 
 
236
      N = KINDEX
 
237
      IF (SUBPAR(N)(1:1).EQ.'@') THEN               !input = @pixel
 
238
         CALL GENCNV(SUBPAR(N)(2:),1,1,IBUF,RR,DD,SLEN)
 
239
         IF (SLEN.LE.0) CALL STETER(3,ERRMS2)
 
240
      ELSE                                          !input = world coord
 
241
         CALL GENCNV(SUBPAR(N),2,1,IAV,TEMP,DD,SLEN)
 
242
         IF (SLEN.LE.0) CALL STETER(3,ERRMS2)
 
243
         IBUF(1) = NINT((TEMP - STARTA(INDXFL))/STEPA(INDXFL)) + 1
 
244
      ENDIF
 
245
      IF (KINDEX.EQ.NCOUNT) THEN
 
246
         IBUF(2) = 0
 
247
      ELSE
 
248
         N = N + 1
 
249
         IF (SUBPAR(N)(1:1).EQ.'@') THEN               !input = @pixel
 
250
            CALL GENCNV(SUBPAR(N)(2:),1,1,IBUF(2),TEMP,DD,SLEN)
 
251
            IF (SLEN.LE.0) CALL STETER(3,ERRMS2)
 
252
         ELSE                                          !input = world coord
 
253
            CALL GENCNV(SUBPAR(N),2,1,IAV,TEMP,DD,SLEN)
 
254
            IF (SLEN.LE.0) CALL STETER(3,ERRMS2)
 
255
            IBUF(2) = NINT((TEMP - STARTA(INDXFL))/STEPA(INDXFL)) + 1
 
256
         ENDIF
 
257
      ENDIF
 
258
C 
 
259
1010  FMIN = RMIN
 
260
      FMAX = RMAX
 
261
      CALL DOIT(ACTION,MADRID(PNTRB),MADRID(PNTRX),MADRID(PNTRY),
 
262
     +          STARTA,STEPA,NPIXA,
 
263
     +          IBUF,FMIN,FMAX,ISTAT)
 
264
 
 
265
      IF (ISTAT.LT.0) THEN
 
266
         CALL STTPUT(MESS(-ISTAT),STAT)
 
267
         ISTAT = 1            !on return, ISTAT = 1 or 2
 
268
      ELSE                    !depending on no. of columns/rows used
 
269
         SLEN = INDEX(SUBPAR(KINDEX),' ')            !cut off trailing blanks
 
270
         IF (SLEN.LE.0) SLEN = LEN(SUBPAR(KINDEX))
 
271
C
 
272
         IF (ISTAT.EQ.1) THEN
 
273
            OUTPUT(1:) = 'column/row at x/y = '
 
274
     +                   //SUBPAR(KINDEX)(1:SLEN)//
 
275
     +                   'cleaned... '
 
276
         ELSE
 
277
            SLEN1 = INDEX(SUBPAR(KINDEX+1),' ')      !cut off trailing blanks
 
278
            IF (SLEN1.LE.0) SLEN1 = LEN(SUBPAR(KINDEX+1))
 
279
            OUTPUT(1:) = 'columns/rows at x/y = '//
 
280
     +                   SUBPAR(KINDEX)(1:SLEN)//', '//
 
281
     +                   SUBPAR(KINDEX+1)(1:SLEN1)//'cleaned...'
 
282
         ENDIF
 
283
         CALL STTPUT(OUTPUT,STAT)
 
284
         RMIN = MIN(RMIN,FMIN)
 
285
         RMAX = MAX(RMAX,FMAX)
 
286
      ENDIF
 
287
      GOTO (1000,3000),TBLFLG                  !look for more...
 
288
C
 
289
C  here we get our input data from a table
 
290
3000  KROW = KROW2 + ISTAT - 1                 !KROW -> KROW2 or -> KROW2+1
 
291
      IF (KROW.LE.NROW) THEN
 
292
         KROW1 = KROW
 
293
         IF (ISTAT.EQ.1) THEN                  !only one row used last time...
 
294
            IBUF(1) = IBUF(2)
 
295
            SUBPAR(1)(1:) = SUBPAR(2)(1:)
 
296
            GOTO 3100                          !so reuse last value
 
297
         ENDIF
 
298
C 
 
299
3030     CALL TBSGET(TIDI,KROW1,SELFLG,STAT)       !only use selected rows
 
300
         IF (.NOT.SELFLG) THEN
 
301
            KROW1 = KROW1 + 1                      !get next row
 
302
            IF (KROW1.LE.NROW) THEN
 
303
               GOTO 3030
 
304
            ELSE
 
305
               GOTO 3300                           !end of table reached
 
306
            ENDIF
 
307
         ENDIF
 
308
C 
 
309
         CALL TBRRDR(TIDI,KROW1,1,TABCLN,TEMP,TABNUL,STAT)
 
310
         IBUF(1) = NINT((TEMP - STARTA(INDXFL))/STEPA(INDXFL)) + 1
 
311
         WRITE(SUBPAR(1),10000) TEMP
 
312
         DO 3060, N=1,30
 
313
            IF (SUBPAR(1)(N:N).NE.' ') THEN
 
314
               IF (N.GT.1) SUBPAR(1)(1:) = SUBPAR(1)(N:)//' '
 
315
               GOTO 3100
 
316
            ENDIF
 
317
3060     CONTINUE
 
318
C 
 
319
3100     KROW2 = KROW1 + 1                      !point to next possible row
 
320
         IF (KROW1.EQ.NROW) THEN   
 
321
            IBUF(2) = 0
 
322
         ELSE
 
323
3130        CALL TBSGET(TIDI,KROW2,SELFLG,STAT)       !only use selected rows
 
324
            IF (.NOT.SELFLG) THEN
 
325
               KROW2 = KROW2 + 1                      !get next row
 
326
               IF (KROW2.LE.NROW) THEN
 
327
                  GOTO 3130
 
328
               ELSE
 
329
                  IBUF(2) = 0
 
330
                  GOTO 3200 
 
331
               ENDIF
 
332
            ENDIF
 
333
C 
 
334
            CALL TBRRDR(TIDI,KROW2,1,TABCLN,TEMP,TABNUL,STAT)
 
335
            IBUF(2) = NINT((TEMP - STARTA(INDXFL))/STEPA(INDXFL)) + 1
 
336
            WRITE(SUBPAR(2),10000) TEMP
 
337
            DO 3160, N=1,30
 
338
               IF (SUBPAR(2)(N:N).NE.' ') THEN
 
339
                  IF (N.GT.1) SUBPAR(2)(1:) = SUBPAR(2)(N:)//' '
 
340
                  GOTO 3200
 
341
               ENDIF
 
342
3160        CONTINUE
 
343
         ENDIF
 
344
C 
 
345
3200     KINDEX = 1                              !force index to 1
 
346
         GOTO 1010                               !now continue as above
 
347
      ENDIF
 
348
C 
 
349
C  close table
 
350
3300  CALL TBTCLO(TIDI,STAT)
 
351
C
 
352
C  get new dynamic range + write descriptor LHCUTS
 
353
4000  CUTS(3) = MIN(RMIN,CUTS(3))
 
354
      CUTS(4) = MAX(RMAX,CUTS(4))
 
355
C
 
356
C  copy descriptors + append HISTORY
 
357
      CALL DSCUPT(IMNOA,IMNOB,' ',STAT)
 
358
      CALL STDWRR(IMNOB,'LHCUTS',CUTS,1,4,UNI,STAT)      !use old cuts...
 
359
C
 
360
C  That's it folks...
 
361
      CALL STSEPI
 
362
C
 
363
C  format statements
 
364
10000 FORMAT(F20.10)
 
365
      END
 
366
 
 
367
      SUBROUTINE DOIT(ACTION,A,TEMPA,TEMPB,
 
368
     +                START,STEP,NPIX,COORD,FMIN,FMAX,ISTAT)
 
369
C
 
370
      IMPLICIT NONE
 
371
C
 
372
      INTEGER      COORD(2),ISTAT
 
373
      INTEGER      JOF,JOF1,N,M,NFLAG,KFLAG
 
374
      INTEGER      NEXT,NP,JNDX,NOLINS
 
375
      INTEGER      NPIX(2),INDX(4),IOF(4)
 
376
C
 
377
      CHARACTER    ACTION*2
 
378
C
 
379
      REAL         A(*),TEMPA(*),TEMPB(*),FMIN,FMAX
 
380
      REAL         V,VV,W,WW,XX,XXX
 
381
      REAL         RMIN,RMAX
 
382
C
 
383
      DOUBLE PRECISION START(2),STEP(2)
 
384
      DOUBLE PRECISION   S(4),Z(4),X(4),F(4)
 
385
      DOUBLE PRECISION   RR1,RR2,RR3,P1,P2,P3,Q
 
386
      DOUBLE PRECISION   SUM1,SUM2
 
387
      DOUBLE PRECISION   A0,A1,A2
 
388
C
 
389
      IF (ACTION(1:1).EQ.'C') THEN
 
390
         NFLAG = 1                              !columns...
 
391
      ELSE
 
392
         NFLAG = 2                              !rows...
 
393
      ENDIF
 
394
      KFLAG = 3 - NFLAG 
 
395
      INDX(1) = COORD(1) - 2                    !take 2 pixels less
 
396
      NEXT = COORD(2) - 2
 
397
C
 
398
      INDX(2) = INDX(1) + 1
 
399
      JNDX = INDX(2) + 1                        !keep index of bad line...
 
400
      IF (NEXT .EQ. INDX(1)+1) THEN
 
401
         NOLINS = 2                             !double lines...
 
402
      ELSE
 
403
         NOLINS = 1
 
404
      ENDIF
 
405
      INDX(3) = INDX(2) + NOLINS + 1                  !single line...
 
406
      INDX(4) = INDX(3) + 1
 
407
C 
 
408
C  check input pars
 
409
      IF ((INDX(1).LT.1).OR.(INDX(4).GT.NPIX(NFLAG))) THEN
 
410
         ISTAT = -1
 
411
         RETURN
 
412
      ELSE
 
413
         ISTAT = NOLINS
 
414
      ENDIF
 
415
C
 
416
C  get "x" values + related sums for later calculations
 
417
      DO 200, N=1,4                              !get x values...
 
418
         X(N) = START(NFLAG) + (INDX(N)-1)*STEP(NFLAG)
 
419
         Z(N) = 1.D0
 
420
200   CONTINUE
 
421
      XX = START(NFLAG) + (JNDX-1)*STEP(NFLAG)      !x of bad row/column...
 
422
      XXX = XX + STEP(NFLAG)
 
423
C
 
424
      DO 1000, M=1,4
 
425
         DO 400, N=1,4                               !get x**2,x**3,x**4 ...
 
426
            Z(N) = Z(N)*X(N)
 
427
400      CONTINUE
 
428
         S(M) = 0.D0                                !compute sums
 
429
         DO 600, N=1,4
 
430
            S(M) = S(M) + Z(N)
 
431
600      CONTINUE
 
432
1000  CONTINUE
 
433
C
 
434
C  branch according to columns + rows
 
435
      IF (ACTION(1:1).EQ.'C') GOTO 5000
 
436
C
 
437
C  handle rows...
 
438
      DO 1200, N=1,4                                 !get y values...
 
439
         IOF(N) = (INDX(N)-1)*NPIX(1)
 
440
1200  CONTINUE
 
441
      JOF = (JNDX-1)*NPIX(1)
 
442
      JOF1 = JOF + NPIX(1)
 
443
C
 
444
C  go through the row
 
445
      DO 3300, NP=1,NPIX(1)
 
446
         DO 1400, M=1,4                              !get function values
 
447
            F(M) = A(IOF(M)+NP)
 
448
            Z(M) = 0.D0
 
449
1400     CONTINUE
 
450
C
 
451
         DO 1600, M=1,4                              !compute sums on right side
 
452
            Z(1) = Z(1) + F(M)
 
453
            Z(2) = Z(2) + F(M)*X(M)
 
454
            Z(3) = Z(3) + F(M)*X(M)*X(M)
 
455
1600     CONTINUE
 
456
C
 
457
C  reduce + solve the linear system  4.*A0 + S1*A1 + S2*A2 = Z1
 
458
C                                    S1*A0 + S2*A1 + S3*A2 = Z2
 
459
C                                    S2*A0 + S3*A1 + S4*A2 = Z3
 
460
         Q = S(1) * 0.25D0
 
461
         RR1 = S(2) - S(1)*Q
 
462
         RR2 = S(3) - S(2)*Q
 
463
         RR3 = Z(2) - Z(1)*Q
 
464
         Q = S(2) * 0.25D0
 
465
         P1 = S(3) - S(1)*Q
 
466
         P2 = S(4) - S(2)*Q
 
467
         P3 = Z(3) - Z(1)*Q
 
468
C                   we now have:          RR1*A1 + RR2*A2 = RR3
 
469
C                                         P1*A1 + P2*A2 = P3
 
470
         Q = P1/RR1
 
471
         IF (ABS(P2-RR2*Q).LE.10.D-30) THEN
 
472
            IF (ABS(P1 - RR1*Q).LE.10.D-30) THEN      !linearly dependant system
 
473
               A0 = 0.D0
 
474
               A1 = 1.D0
 
475
               A2 = 0.D0
 
476
               GOTO 3000
 
477
            ENDIF
 
478
            Q = P2/RR2
 
479
            A1 = (P3 - RR3*Q) / (P1 - RR1*Q)
 
480
            A2 = (RR3 - A1*RR1)/RR2
 
481
         ELSE
 
482
            A2 = (P3 - RR3*Q) / (P2 - RR2*Q)
 
483
            A1 = (RR3 - A2*RR2)/RR1
 
484
         ENDIF
 
485
C
 
486
         A0 = (Z(1) - S(2)*A2 - S(1)*A1) * .25D0
 
487
C
 
488
C  finally calculate interpolated value...!
 
489
3000     TEMPA(NP) = A0 + A1*XX + A2*XX*XX
 
490
         IF (NOLINS.EQ.2)
 
491
     +      TEMPB(NP) = A0 + A1*XXX + A2*XXX*XXX
 
492
3300  CONTINUE
 
493
C
 
494
C  check, if we work on single line or not
 
495
C
 
496
      IF (NOLINS.EQ.2) GOTO 4000
 
497
C
 
498
C  for constant offset, get average difference to orignal bad row
 
499
      IF (ACTION(2:2).EQ.'C') THEN
 
500
         SUM1 = 0.D0
 
501
         DO 3500, NP=1,NPIX(1)
 
502
            SUM1 = SUM1 + DBLE(A(JOF+NP) - TEMPA(NP))  !sum differences
 
503
3500     CONTINUE
 
504
         V = SNGL(SUM1/NPIX(1))            !get average offset
 
505
C
 
506
C  and correct pixels in bad row with that offset
 
507
         DO 3600, NP=1,NPIX(1)
 
508
            W = A(JOF+NP)
 
509
            A(JOF+NP) = W - V
 
510
            IF (W.LT.FMIN) THEN
 
511
               FMIN = W
 
512
            ELSE
 
513
               IF (W.GT.FMAX) FMAX = W
 
514
            ENDIF
 
515
3600     CONTINUE
 
516
C
 
517
      ELSE
 
518
C
 
519
C  store approximated pixels in bad row
 
520
         DO 3800, NP=1,NPIX(1)
 
521
            W = TEMPA(NP)
 
522
            A(JOF+NP) = W
 
523
            IF (W.LT.FMIN) THEN
 
524
               FMIN = W
 
525
            ELSE
 
526
               IF (W.GT.FMAX) FMAX = W
 
527
            ENDIF
 
528
3800     CONTINUE
 
529
      ENDIF
 
530
      RETURN
 
531
C
 
532
C  here to handle two adjacent "bad" rows
 
533
C
 
534
4000  IF (ACTION(2:2).EQ.'C') THEN
 
535
         SUM1 = 0.D0                  !get average difference to adjacent lines
 
536
         SUM2 = 0.D0
 
537
         DO 4200, NP=1,NPIX(1)
 
538
            SUM1 = SUM1 + DBLE(A(JOF+NP) - TEMPA(NP))   !sum differences
 
539
            SUM2 = SUM2 + DBLE(A(JOF1+NP) - TEMPB(NP))  !sum differences
 
540
4200     CONTINUE
 
541
         V = SNGL(SUM1/NPIX(1))                         !get average offset
 
542
         VV = SNGL(SUM2/NPIX(1))
 
543
C
 
544
C  and correct pixels in bad rows with these offsets
 
545
         DO 4400, NP=1,NPIX(1)
 
546
            W = A(JOF+NP)
 
547
            A(JOF+NP) = W - V
 
548
            WW = A(JOF1+NP)
 
549
            A(JOF1+NP) = WW - VV
 
550
            RMIN = MIN(W,WW)
 
551
            RMAX = MAX(W,WW)
 
552
            IF (RMIN.LT.FMIN) FMIN = RMIN
 
553
            IF (RMAX.GT.FMAX) FMAX = RMAX
 
554
4400     CONTINUE
 
555
C
 
556
      ELSE
 
557
C
 
558
C  store approximated pixels in bad rows
 
559
         DO 4600, NP=1,NPIX(1)
 
560
            W = TEMPA(NP)
 
561
            WW = TEMPB(NP)
 
562
            A(JOF+NP) = W
 
563
            A(JOF1+NP) = WW
 
564
            RMIN = MIN(W,WW)
 
565
            RMAX = MAX(W,WW)
 
566
            IF (RMIN.LT.FMIN) FMIN = RMIN
 
567
            IF (RMAX.GT.FMAX) FMAX = RMAX
 
568
4600     CONTINUE
 
569
      ENDIF
 
570
      RETURN
 
571
C
 
572
C  handle columns...
 
573
C
 
574
C  go through the column
 
575
5000  JOF = 0
 
576
      DO 6600, NP=1,NPIX(2)
 
577
         DO 5200, M=1,4                        !get function values
 
578
            F(M) = A(JOF+INDX(M))
 
579
            Z(M) = 0.D0
 
580
5200     CONTINUE
 
581
C
 
582
         DO 5400, M=1,4                        !compute sums on right side
 
583
            Z(1) = Z(1) + F(M)
 
584
            Z(2) = Z(2) + F(M)*X(M)
 
585
            Z(3) = Z(3) + F(M)*X(M)*X(M)
 
586
5400     CONTINUE
 
587
C
 
588
C  reduce + solve the linear system     4.*A0 + S1*A1 + S2*A2 = Z1
 
589
C                                       S1*A0 + S2*A1 + S3*A2 = Z2
 
590
C                                       S2*A0 + S3*A1 + S4*A2 = Z3
 
591
         Q = S(1)*.25D0
 
592
         RR1 = S(2) - S(1)*Q
 
593
         RR2 = S(3) - S(2)*Q
 
594
         RR3 = Z(2) - Z(1)*Q
 
595
         Q = S(2)*.25D0
 
596
         P1 = S(3) - S(1)*Q
 
597
         P2 = S(4) - S(2)*Q
 
598
         P3 = Z(3) - Z(1)*Q
 
599
C  we now have:                         RR1*A1 + RR2*A2 = RR3
 
600
C                                       P1*A1 + P2*A2 = P3
 
601
         Q = P1/RR1
 
602
         IF (ABS(P2-RR2*Q).LE.10.D-30) THEN
 
603
            IF (ABS(P1 - RR1*Q).LE.10.D-30) THEN      !linearly dependant system
 
604
               A0 = 0.D0
 
605
               A1 = 1.D0
 
606
               A2 = 0.D0
 
607
               GOTO 6000
 
608
            ENDIF
 
609
            Q = P2/RR2
 
610
            A1 = (P3 - RR3*Q) / (P1 - RR1*Q)
 
611
            A2 = (RR3 - A1*RR1)/RR2
 
612
         ELSE
 
613
            A2 = (P3 - RR3*Q) / (P2 - RR2*Q)
 
614
            A1 = (RR3 - A2*RR2)/RR1
 
615
         ENDIF
 
616
C
 
617
         A0 = (Z(1) - S(2)*A2 - S(1)*A1)*.25D0
 
618
C
 
619
C  finally calculate interpolated value...!
 
620
6000     TEMPA(NP) = A0 + A1*XX + A2*XX*XX
 
621
         IF (NOLINS.EQ.2)
 
622
     +      TEMPB(NP) = A0 + A1*XXX + A2*XXX*XXX
 
623
         JOF = JOF + NPIX(1)                  !update offset to next row...
 
624
6600  CONTINUE
 
625
C
 
626
C  check, if we work on single line or not
 
627
C
 
628
      IF (NOLINS.EQ.2) GOTO 8000
 
629
C
 
630
C  for constant offset, get average difference to original bad column
 
631
      IF (ACTION(2:2).EQ.'C') THEN
 
632
         SUM1 = 0.D0
 
633
         JOF = JNDX
 
634
         DO 6800, NP=1,NPIX(2)
 
635
            SUM1 = SUM1 + DBLE(A(JOF) - TEMPA(NP))      !sum differences
 
636
            JOF = JOF + NPIX(1)
 
637
6800     CONTINUE
 
638
         W = SNGL(SUM1/NPIX(2))                         !get average offset
 
639
C
 
640
C  and correct pixels in bad column with that offset
 
641
         JOF = JNDX
 
642
         DO 7000, NP=1,NPIX(2)
 
643
            V = A(JOF) - W
 
644
            A(JOF) = V 
 
645
            IF (V.LT.FMIN) THEN
 
646
               FMIN = V
 
647
            ELSE
 
648
               IF (V.GT.FMAX) FMAX = V
 
649
            ENDIF
 
650
            JOF = JOF + NPIX(1)
 
651
7000     CONTINUE
 
652
C
 
653
      ELSE
 
654
C
 
655
C  fill approximated pixels into frame
 
656
         JOF = JNDX
 
657
         DO 7400, NP=1,NPIX(2)
 
658
            V = TEMPA(NP)
 
659
            A(JOF) = V
 
660
            IF (V.LT.FMIN) THEN
 
661
               FMIN = V
 
662
            ELSE
 
663
               IF (V.GT.FMAX) FMAX = V
 
664
            ENDIF
 
665
            JOF = JOF + NPIX(1)
 
666
7400     CONTINUE
 
667
      ENDIF
 
668
      RETURN
 
669
C
 
670
C  here to handle two adjacent "bad" columns
 
671
8000  IF (ACTION(2:2).EQ.'C') THEN
 
672
         SUM1 = 0.D0                  !get average difference to adjacent lines
 
673
         SUM2 = 0.D0
 
674
         JOF = JNDX
 
675
C
 
676
         DO 8200, NP=1,NPIX(2)
 
677
            SUM1 = SUM1 + DBLE(A(JOF) - TEMPA(NP))        !sum differences
 
678
            SUM2 = SUM2 + DBLE(A(JOF+1) - TEMPB(NP))
 
679
            JOF = JOF + NPIX(1)
 
680
8200     CONTINUE
 
681
         V = SNGL(SUM1/NPIX(2))                           !get average offset
 
682
         VV = SNGL(SUM2/NPIX(2))
 
683
C
 
684
C  and correct pixels in bad columns with these offsets
 
685
         JOF = JNDX
 
686
         DO 8400, NP=1,NPIX(2)
 
687
            W = A(JOF)
 
688
            WW = A(JOF+1)
 
689
            A(JOF) = W - V
 
690
            A(JOF+1) = WW - VV
 
691
            RMIN = MIN(W,WW)
 
692
            RMAX = MAX(W,WW)
 
693
            IF (RMIN.LT.FMIN) FMIN = RMIN
 
694
            IF (RMAX.GT.FMAX) FMAX = RMAX
 
695
            JOF = JOF + NPIX(1)
 
696
8400     CONTINUE
 
697
C
 
698
      ELSE
 
699
C
 
700
C  fill approximated pixels into frame
 
701
         JOF = JNDX
 
702
         DO 8600, NP=1,NPIX(2)
 
703
            W = TEMPA(NP)
 
704
            WW = TEMPB(NP)
 
705
            A(JOF) = W
 
706
            A(JOF+1) = WW
 
707
            W = A(JOF)
 
708
            RMIN = MIN(W,WW)
 
709
            RMAX = MAX(W,WW)
 
710
            IF (RMIN.LT.FMIN) FMIN = RMIN
 
711
            IF (RMAX.GT.FMAX) FMAX = RMAX
 
712
            JOF = JOF + NPIX(1)
 
713
8600     CONTINUE
 
714
      ENDIF
 
715
      RETURN
 
716
C
 
717
      END