C @(#)rctint.for 17.1.1.1 (ES0-DMD) 01/25/02 17:19:01 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 C +++ C.PURPOSE: Interpolation C.ALGORYTHM: Each input pixel is divided into 9 subpixels C Their indices are defined as follows: C 1,3 2,3 3,3 FIRST INDEX: X C 1,2 2,2 3,2 C 1,1 2,1 3,1 SECOND INDEX: Y C C For the interpolation, the following weighting scheme is used: C .1667 .1667 .333 .1667 .1667 C -------------------------------- C .1667 | .5 .667 .5 | .1667 C .333 | .667 (1.) .667 | .333 C .1667 | .5 .667 .5 | .1667 C -------------------------------- C .1667 .1667 .333 .1667 .1667 C C (the box marks the edges of one input pixel) C C.NOTE: This "deconvolution" makes an implicit assumption about C the point spread function, namely that its fwhm is smaller C than or roughly equal to one input pixel. future improvements C could include reading in a more realistic or even exact C point spread function and doing a proper deconvolution. C --- SUBROUTINE RCTINT(IN,NPX,NPY,IPOL) C IMPLICIT NONE REAL IN(1) INTEGER NPX INTEGER NPY REAL IPOL(3*NPX,3*NPY) C REAL SUB(3,3) INTEGER NPERC INTEGER INDXX INTEGER INDXY CHARACTER PROGR*50 INTEGER ISTAT INTEGER IPOS, TWO INTEGER I, J, K, L C 9000 FORMAT ('*** INFO: ',I3,' percent completed ...') C CALL STTPUT('*** INFO: Deconvolution started',ISTAT) NPERC = 0 TWO = 2 C C *** Begin with all pixels on the edge of the input frame C *** lower left corner C DO 20 K = 1,3 DO 10 L = 1,3 SUB(K,L) = IN(1) 10 CONTINUE 20 CONTINUE C SUB(3,1) = .667*IN(1) + .333*IN(TWO) ! first index: x, second: y SUB(3,2) = SUB(3,1) SUB(3,3) = .5*IN(1) + .1667* (IN(TWO)+IN(NPX+1)+IN(NPX+2)) SUB(2,3) = .667*IN(1) + .333*IN(NPX+1) SUB(1,3) = SUB(2,3) CALL RCTSCL(SUB,IN(1)) C DO 40 K = 1,3 DO 30 L = 1,3 INDXX = K INDXY = L IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L) 30 CONTINUE 40 CONTINUE C C *** lower right corner C DO 60 K = 1,3 DO 50 L = 1,3 SUB(K,L) = IN(NPX) 50 CONTINUE 60 CONTINUE C SUB(1,1) = .667*IN(NPX) + .333*IN(NPX-1) SUB(1,2) = SUB(1,1) SUB(1,3) = .5*IN(NPX) + .1667* (IN(NPX-1)+IN(2*NPX-1)+IN(2*NPX)) SUB(2,3) = .667*IN(NPX) + .333*IN(2*NPX) SUB(3,3) = SUB(2,3) C CALL RCTSCL(SUB,IN(NPX)) C DO 80 K = 1,3 DO 70 L = 1,3 INDXX = NPX*3 - 3 + K INDXY = L IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L) 70 CONTINUE 80 CONTINUE C C *** upper left corner C IPOS = NPX*(NPY-1) + 1 DO 100 K = 1,3 DO 90 L = 1,3 SUB(K,L) = IN(IPOS) 90 CONTINUE 100 CONTINUE C SUB(1,1) = .667*IN(IPOS) + .333*IN(IPOS-NPX) SUB(2,1) = SUB(1,1) SUB(3,1) = .5*IN(IPOS) + .1667* (IN(IPOS-NPX)+IN(IPOS-NPX-1)+ + IN(IPOS+1)) SUB(3,2) = .667*IN(IPOS) + .333*IN(IPOS+1) SUB(3,3) = SUB(3,2) CALL RCTSCL(SUB,IN(IPOS)) C DO 120 K = 1,3 DO 110 L = 1,3 INDXX = K INDXY = NPY*3 - 3 + L IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L) 110 CONTINUE 120 CONTINUE C C *** upper right corner C IPOS = NPX*NPY DO 140 K = 1,3 DO 130 L = 1,3 SUB(K,L) = IN(IPOS) 130 CONTINUE 140 CONTINUE C SUB(1,3) = .667*IN(IPOS) + .333*IN(IPOS-1) SUB(1,2) = SUB(1,3) SUB(3,1) = .5*IN(IPOS) + .1667* (IN(IPOS-1)+IN(IPOS-NPX-1)+ + IN(IPOS-NPX)) SUB(2,1) = .667*IN(IPOS) + .333*IN(IPOS-NPX) SUB(3,1) = SUB(2,1) CALL RCTSCL(SUB,IN(IPOS)) C DO 160 K = 1,3 DO 150 L = 1,3 INDXX = NPX*3 - 3 + K INDXY = NPY*3 - 3 + L IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L) 150 CONTINUE 160 CONTINUE C C *** bottom edge C DO 210 I = 2,NPX - 1 IPOS = I DO 180 K = 1,3 DO 170 L = 1,3 SUB(K,L) = IN(IPOS) 170 CONTINUE 180 CONTINUE C SUB(1,1) = .667*IN(IPOS) + .333*IN(IPOS-1) SUB(1,2) = SUB(1,1) SUB(1,3) = .5*IN(IPOS) + .1667* + (IN(IPOS-1)+IN(IPOS+NPX-1)+IN(IPOS+NPX)) SUB(2,3) = .667*IN(IPOS) + .333*IN(IPOS+NPX) SUB(3,3) = .5*IN(IPOS) + .1667* + (IN(IPOS+1)+IN(IPOS+NPX)+IN(IPOS+NPX+1)) SUB(3,2) = .667*IN(IPOS) + .333*IN(IPOS+1) SUB(3,1) = SUB(3,2) CALL RCTSCL(SUB,IN(IPOS)) C DO 200 K = 1,3 DO 190 L = 1,3 INDXX = I*3 - 3 + K INDXY = L IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L) 190 CONTINUE 200 CONTINUE 210 CONTINUE C C *** top edge C DO 260 I = 2,NPX - 1 IPOS = (NPY-1)*NPX + I DO 230 K = 1,3 DO 220 L = 1,3 SUB(K,L) = IN(IPOS) 220 CONTINUE 230 CONTINUE C SUB(1,3) = .667*IN(IPOS) + .333*IN(IPOS-1) SUB(1,2) = SUB(1,3) SUB(1,1) = .5*IN(IPOS) + .1667* + (IN(IPOS-1)+IN(IPOS-NPX-1)+IN(IPOS-NPX)) SUB(2,1) = .667*IN(IPOS) + .333*IN(IPOS-NPX) SUB(3,1) = .5*IN(IPOS) + .1667* + (IN(IPOS-NPX)+IN(IPOS-NPX)+IN(IPOS-NPX+1)) SUB(3,2) = .667*IN(IPOS) + .333*IN(IPOS+1) SUB(3,3) = SUB(3,2) CALL RCTSCL(SUB,IN(IPOS)) C DO 250 K = 1,3 DO 240 L = 1,3 INDXX = I*3 - 3 + K INDXY = NPY*3 - 3 + L IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L) 240 CONTINUE 250 CONTINUE 260 CONTINUE C C *** left edge C DO 310 J = 2,NPY - 1 IPOS = (J-1)*NPX + 1 DO 280 K = 1,3 DO 270 L = 1,3 SUB(K,L) = IN(IPOS) 270 CONTINUE 280 CONTINUE C SUB(1,1) = .667*IN(IPOS) + .333*IN(IPOS-NPX) SUB(2,1) = SUB(1,1) SUB(3,1) = .5*IN(IPOS) + .1667* + (IN(IPOS-NPX)+IN(IPOS-NPX+1)+IN(IPOS+1)) SUB(3,2) = .667*IN(IPOS) + .333*IN(IPOS+1) SUB(3,3) = .5*IN(IPOS) + .1667* + (IN(IPOS+1)+IN(IPOS+NPX+1)+IN(IPOS+NPX)) SUB(2,3) = .667*IN(IPOS) + .333*IN(IPOS+NPX) SUB(1,3) = SUB(2,3) CALL RCTSCL(SUB,IN(IPOS)) C DO 300 K = 1,3 DO 290 L = 1,3 INDXX = K INDXY = J*3 - 3 + L IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L) 290 CONTINUE 300 CONTINUE 310 CONTINUE C C *** right edge C DO 360 J = 2,NPY - 1 IPOS = J*NPX DO 330 K = 1,3 DO 320 L = 1,3 SUB(K,L) = IN(IPOS) 320 CONTINUE 330 CONTINUE C SUB(3,1) = .667*IN(IPOS) + .333*IN(IPOS-NPX) SUB(2,1) = SUB(3,1) SUB(1,1) = .5*IN(IPOS) + .1667* + (IN(IPOS-NPX)+IN(IPOS-NPX-1)+IN(IPOS-1)) SUB(1,2) = .667*IN(IPOS) + .333*IN(IPOS-1) SUB(1,3) = .5*IN(IPOS) + .1667* + (IN(IPOS-1)+IN(IPOS+NPX-1)+IN(IPOS+NPX)) SUB(2,3) = .667*IN(IPOS) + .333*IN(IPOS+NPX) SUB(3,3) = SUB(2,3) CALL RCTSCL(SUB,IN(IPOS)) C DO 350 K = 1,3 DO 340 L = 1,3 INDXX = NPX*3 - 3 + K INDXY = J*3 - 3 + L IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L) 340 CONTINUE 350 CONTINUE 360 CONTINUE C C *** and now the rest is straight forward C DO 400 J = 2,NPY - 1 DO 390 I = 2,NPX - 1 IPOS = (J-1)*NPX + I ! CENTRAL SUBPIXEL SUB(2,2) = IN(IPOS) SUB(1,1) = .5*IN(IPOS) + .1667* + (IN(IPOS-1)+IN(IPOS-NPX-1)+IN(IPOS-NPX)) SUB(2,1) = .667*IN(IPOS) + .333*IN(IPOS-NPX) SUB(3,1) = .5*IN(IPOS) + .1667* + (IN(IPOS-NPX)+IN(IPOS-NPX+1)+IN(IPOS+1)) SUB(3,2) = .667*IN(IPOS) + .333*IN(IPOS+1) SUB(3,3) = .5*IN(IPOS) + .1667* + (IN(IPOS+1)+IN(IPOS+NPX+1)+IN(IPOS+NPX)) SUB(2,3) = .667*IN(IPOS) + .333*IN(IPOS+NPX) SUB(1,3) = .5*IN(IPOS) + .1667* + (IN(IPOS+NPX)+IN(IPOS+NPX-1)+IN(IPOS-1)) SUB(1,2) = .667*IN(IPOS) + .333*IN(IPOS-1) CALL RCTSCL(SUB,IN(IPOS)) C DO 380 K = 1,3 DO 370 L = 1,3 INDXX = I*3 - 3 + K INDXY = J*3 - 3 + L IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L) 370 CONTINUE 380 CONTINUE 390 CONTINUE C C *** feed them progress messages ... C IF (((10*J)/NPY).GE. (NPERC+1)) THEN NPERC = NPERC + 1 WRITE (PROGR,9000) 10*NPERC CALL STTPUT(PROGR,ISTAT) END IF 400 CONTINUE RETURN END