C C DRUTIL.F C C Utility routines for Drizzle and Blot C C These routines should call only themselves and the STSDAS f77/vos C library. C C This version for DRIZZLE 2.2 and BLOT 0.6, January 2001 C C Richard Hook, ST-ECF/ESO C C History: C C The module GOLDH was transferred from drizzle as it is also C now needed by blot. C C WCS routines added (XY2RD & RD2XY), September 1998 C C Modified GTCOEF routine to also handle STIS and NICMOS, C January 1999 C C Added additional code for reversing the cubic transformation added C June 1999 C C Major extensions and some routines transferred from the EIS utilities C collection (for context handling), October 2000 C C Extensive modifications for a more flexible geometric coefficients C scheme, changes in calling sequences, December 2000 C SUBROUTINE SETIM(A,NX,NY,V) C C Set a 2D array to a specified value C IMPLICIT NONE INTEGER NX,NY REAL A(NX,NY),V INTEGER I,J DO J=1,NY DO I=1,NX A(I,J)=V ENDDO ENDDO RETURN END SUBROUTINE MULC(A,NX,NY,V) C C Multiply a 2D array to a specified value C IMPLICIT NONE INTEGER NX,NY REAL A(NX,NY),V INTEGER I,J DO J=1,NY DO I=1,NX A(I,J)=A(I,J)*V ENDDO ENDDO RETURN END SUBROUTINE COPYIM(IN,OUT,NX,NY) C C Simply copy an array into another C IMPLICIT NONE INTEGER NX,NY REAL IN(NX,NY),OUT(NX,NY) INTEGER I,J DO J=1,NY DO I=1,NX OUT(I,J)=IN(I,J) ENDDO ENDDO RETURN END REAL FUNCTION EVALCU(X,Y,CO) C C Evaluate a cubic geometric distortion in 2d C IMPLICIT NONE REAL X,Y,CO(*) EVALCU=CO(1)+ : CO(2)*X+ : CO(3)*Y+ : CO(4)*X*X+ : CO(5)*X*Y+ : CO(6)*Y*Y+ : CO(7)*X*X*X+ : CO(8)*X*X*Y+ : CO(9)*X*Y*Y+ : CO(10)*Y*Y*Y RETURN END REAL FUNCTION EVAL4(X,Y,CO) C C Evaluate a 4th order (quartic) geometric distortion in 2d C IMPLICIT NONE REAL X,Y,CO(*) EVAL4=CO(1)+ : CO(2)*X+ : CO(3)*Y+ : CO(4)*X*X+ : CO(5)*X*Y+ : CO(6)*Y*Y+ : CO(7)*X*X*X+ : CO(8)*X*X*Y+ : CO(9)*X*Y*Y+ : CO(10)*Y*Y*Y+ : CO(11)*X*X*X*X+ : CO(12)*X*X*X*Y+ : CO(13)*X*X*Y*Y+ : CO(14)*X*Y*Y*Y+ : CO(15)*Y*Y*Y*Y RETURN END REAL FUNCTION EVAL5(X,Y,CO) C C Evaluate a 5th order geometric distortion in 2d C IMPLICIT NONE REAL X,Y,CO(*) EVAL5=CO(1)+ : CO(2)*X+ : CO(3)*Y+ : CO(4)*X*X+ : CO(5)*X*Y+ : CO(6)*Y*Y+ : CO(7)*X*X*X+ : CO(8)*X*X*Y+ : CO(9)*X*Y*Y+ : CO(10)*Y*Y*Y+ : CO(11)*X*X*X*X+ : CO(12)*X*X*X*Y+ : CO(13)*X*X*Y*Y+ : CO(14)*X*Y*Y*Y+ : CO(15)*Y*Y*Y*Y+ : CO(16)*X*X*X*X*X+ : CO(17)*X*X*X*X*Y+ : CO(18)*X*X*X*Y*Y+ : CO(19)*X*X*Y*Y*Y+ : CO(20)*X*Y*Y*Y*Y+ : CO(21)*Y*Y*Y*Y*Y RETURN END SUBROUTINE RAD3(X,Y,CO,XO,YO) C C Evaluate a 3rd order radial geometric distortion in 2d C X version. Note that there is no zero order coefficient C as this is physically meaningless. C IMPLICIT NONE REAL X,Y,XO,YO,R,F,CO(*) R=SQRT(X*X+Y*Y) F=1.0+CO(1)+CO(2)*R+CO(3)*R*R XO=F*X YO=F*Y RETURN END SUBROUTINE GETCO(LUN,LAM,COTY,COMAX,CONUM,XCO,YCO,ISTAT) C C Read geometric distortion coefficients in free format C from a text file which has already been opened and return them. C C If there is a problem reading the numbers then the status flag will be C returned as 1. If all is well it will be returned C as 0. C C This routine caused problems with SUN SPARC SC4 compilers and was C modified to avoid them. March 1997 C C Modified to handle higher-order terms in addition. C December 2000 C IMPLICIT NONE INTEGER LUN,ISTAT,I,J CHARACTER*80 LINE CHARACTER*1024 BUFFER REAL LAM,MGF2,N REAL A,B,C INTEGER COTY,COMAX,CONUM REAL XCO(COMAX),YCO(COMAX) C Skip blank lines and those beginning with # until we find a label DO WHILE(.TRUE.) CALL UFGLIN(LUN,LINE,ISTAT) IF(ISTAT.NE.0) GO TO 99 IF(LINE.NE.' ' .AND. LINE(1:1).NE.'#') THEN C First "Trauger" style coefficients - 60 of them IF(LINE(1:7).EQ.'trauger') THEN C Now we need to calculate the coefficients which are a C function of wavelength N=MGF2(LAM) C Now we loop through extracting the coefficients, 3 per line, C and calculating the wavelength dependence J=1 DO WHILE(J.LE.20) CALL UFGLIN(LUN,LINE,ISTAT) IF(ISTAT.NE.0) GO TO 99 IF(LINE.NE.' ' .AND. LINE(1:1).NE.'#') THEN READ(LINE,*,IOSTAT=ISTAT) A,B,C IF(ISTAT.NE.0) GO TO 99 IF(J.LE.10) THEN XCO(J)=A+B*(N-1.5)+C*(N-1.5)**2 ELSE YCO(J-10)=A+B*(N-1.5)+C*(N-1.5)**2 ENDIF J=J+1 ENDIF ENDDO COTY=3 CONUM=10 ISTAT=0 GO TO 99 C Next "cubic" style - 20 coefficients ELSE IF(LINE(1:5).EQ.'cubic' .OR. : LINE(1:5).EQ.'poly3') THEN C Copy the rest of the file into a buffer ready for a free-format C read operation CALL BFILL(LUN,BUFFER,ISTAT) CONUM=10 IF(ISTAT.EQ.0) THEN READ(BUFFER,*,IOSTAT=ISTAT) : (XCO(I),I=1,CONUM),(YCO(I),I=1,CONUM) ENDIF IF(ISTAT.NE.0) GO TO 99 COTY=3 ISTAT=0 GO TO 99 C Now 4th order - 30 coefficients ELSE IF(LINE(1:5).EQ.'poly4' .OR. : LINE(1:7).EQ.'quartic') THEN CALL BFILL(LUN,BUFFER,ISTAT) CONUM=15 IF(ISTAT.EQ.0) THEN READ(BUFFER,*,IOSTAT=ISTAT) : (XCO(I),I=1,CONUM),(YCO(I),I=1,CONUM) ENDIF IF(ISTAT.NE.0) GO TO 99 COTY=4 ISTAT=0 GO TO 99 C Now 5th order - 42 coefficients ELSE IF(LINE(1:5).EQ.'poly5' .OR. : LINE(1:7).EQ.'quintic') THEN CALL BFILL(LUN,BUFFER,ISTAT) CONUM=21 IF(ISTAT.EQ.0) THEN READ(BUFFER,*,IOSTAT=ISTAT) : (XCO(I),I=1,CONUM),(YCO(I),I=1,CONUM) ENDIF IF(ISTAT.NE.0) GO TO 99 COTY=5 ISTAT=0 GO TO 99 C Next third-order radial polynomial ELSE IF(LINE(1:6).EQ.'radial' .OR. : LINE(1:4).EQ.'rad3') THEN C Copy the rest of the file into a buffer ready for a free-format C read operation CALL BFILL(LUN,BUFFER,ISTAT) CONUM=3 IF(ISTAT.EQ.0) THEN READ(BUFFER,*,IOSTAT=ISTAT) : (XCO(I),I=1,CONUM) ENDIF IF(ISTAT.NE.0) GO TO 99 COTY=-3 ISTAT=0 GO TO 99 ELSE CALL UMSPUT('! Unknown coefficient set type', : 1,0,ISTAT) ISTAT=1 GO TO 99 ENDIF ENDIF ENDDO 99 CONTINUE RETURN END REAL FUNCTION MGF2(LAM) C C Calculate the refractive index of MgF2 for a given C wavelength (in nm) using the formula given by Trauger (1995) C IMPLICIT NONE REAL LAM,SIG SIG=1.0E7/LAM MGF2=SQRT(1.0 + 2.590355E10/(5.312993E10-SIG*SIG) : + 4.4543708E9/(11.17083E9-SIG*SIG) : + 4.0838897E5/(1.766361E5-SIG*SIG)) RETURN END SUBROUTINE GETWCS(ID,WCS,CTYPE1,CTYPE2,ISTAT) C C Get the WCS header items (8 values) and CTYPE strings C C The image is assumed to be available with descriptor ID C and the status flag is returned as non-zero if a parameter C is not found C IMPLICIT NONE INTEGER ID,ISTAT,IS CHARACTER*8 CTYPE1,CTYPE2 DOUBLE PRECISION WCS(8) ISTAT=0 CALL UHDGST(ID,'CTYPE1',CTYPE1,ISTAT) IF(ISTAT.NE.0) CALL UMSPUT( : '! Warning, unable to read CTYPE1 from header',1,0,IS) IF(ISTAT.NE.0) RETURN CALL UHDGST(ID,'CTYPE2',CTYPE2,ISTAT) IF(ISTAT.NE.0) CALL UMSPUT( : '! Warning, unable to read CTYPE2 from header',1,0,IS) IF(ISTAT.NE.0) RETURN CALL UHDGSD(ID,'CRPIX1',WCS(1),ISTAT) IF(ISTAT.NE.0) CALL UMSPUT( : '! Warning, unable to read CRPIX1 from header',1,0,IS) IF(ISTAT.NE.0) RETURN CALL UHDGSD(ID,'CRVAL1',WCS(2),ISTAT) IF(ISTAT.NE.0) CALL UMSPUT( : '! Warning, unable to read CRVAL1 from header',1,0,IS) IF(ISTAT.NE.0) RETURN CALL UHDGSD(ID,'CRPIX2',WCS(3),ISTAT) IF(ISTAT.NE.0) CALL UMSPUT( : '! Warning, unable to read CRPIX2 from header',1,0,IS) IF(ISTAT.NE.0) RETURN CALL UHDGSD(ID,'CRVAL2',WCS(4),ISTAT) IF(ISTAT.NE.0) CALL UMSPUT( : '! Warning, unable to read CRVAL2 from header',1,0,IS) IF(ISTAT.NE.0) RETURN CALL UHDGSD(ID,'CD1_1',WCS(5),ISTAT) IF(ISTAT.NE.0) CALL UMSPUT( : '! Warning, unable to read CD1_1 from header',1,0,IS) IF(ISTAT.NE.0) RETURN CALL UHDGSD(ID,'CD2_1',WCS(6),ISTAT) IF(ISTAT.NE.0) CALL UMSPUT( : '! Warning, unable to read CD2_1 from header',1,0,IS) IF(ISTAT.NE.0) RETURN CALL UHDGSD(ID,'CD1_2',WCS(7),ISTAT) IF(ISTAT.NE.0) CALL UMSPUT( : '! Warning, unable to read CD1_2 from header',1,0,IS) IF(ISTAT.NE.0) RETURN CALL UHDGSD(ID,'CD2_2',WCS(8),ISTAT) IF(ISTAT.NE.0) CALL UMSPUT( : '! Warning, unable to read CD2_2 from header',1,0,IS) RETURN END SUBROUTINE UWCS(ID,WCS,CTYPE1,CTYPE2) C C Update the WCS header items using the 8 values C supplied in WCS and the values of CTYPE1 and CTYPE2. C C The image is assumed to be available with descriptor ID C and to already have all the parameters C C Modified to get rid of CDELT and CROTA keywords if present C and also the PC matrix to avoid confusion. October 2000 C IMPLICIT NONE INTEGER ID,ISTAT CHARACTER*8 CTYPE1,CTYPE2 DOUBLE PRECISION WCS(8) CALL UHDAST(ID,'CTYPE1',CTYPE1,' ',0,ISTAT) CALL UHDAST(ID,'CTYPE2',CTYPE2,' ',0,ISTAT) CALL UHDASD(ID,'CRPIX1',WCS(1),' ',0,ISTAT) CALL UHDASD(ID,'CRVAL1',WCS(2),' ',0,ISTAT) CALL UHDASD(ID,'CRPIX2',WCS(3),' ',0,ISTAT) CALL UHDASD(ID,'CRVAL2',WCS(4),' ',0,ISTAT) CALL UHDASD(ID,'CD1_1',WCS(5),' ',0,ISTAT) CALL UHDASD(ID,'CD2_1',WCS(6),' ',0,ISTAT) CALL UHDASD(ID,'CD1_2',WCS(7),' ',0,ISTAT) CALL UHDASD(ID,'CD2_2',WCS(8),' ',0,ISTAT) C Delete CDELT and CROTA keywords C and the PC matrix CALL UHDDSP(ID,'CDELT1',ISTAT) CALL UHDDSP(ID,'CDELT2',ISTAT) CALL UHDDSP(ID,'CROTA1',ISTAT) CALL UHDDSP(ID,'CROTA2',ISTAT) CALL UHDDSP(ID,'PC001001',ISTAT) CALL UHDDSP(ID,'PC002001',ISTAT) CALL UHDDSP(ID,'PC001002',ISTAT) CALL UHDDSP(ID,'PC002002',ISTAT) RETURN END SUBROUTINE PUTFIL(DAT,COU,ONX,ONY,FILVAL) C C Put in the missing filling value where the weight is C zero C IMPLICIT NONE INTEGER ONX,ONY REAL DAT(ONX,ONY),COU(ONX,ONY),FILVAL INTEGER I,J DO J=1,ONY DO I=1,ONX IF(COU(I,J).EQ.0.0) DAT(I,J)=FILVAL ENDDO ENDDO RETURN END SUBROUTINE GTCOEF(ID,COEFFS,LAMBDA,ISTAT) C C This routine is invoked when it has been requested that C the coefficients file name for geometric distortion, and C the associated wavelength is to be extracted from the C header. C C It is currently handles WFPC2, STIS and NICMOS and assumes C coefficients in a directory drizzle$coeffs. C C Extended for STIS and NICMOS, January 1999 C IMPLICIT NONE INTEGER ID,ISTAT,IDET REAL LAMBDA,PLAM CHARACTER*(*) COEFFS CHARACTER*8 INST,DETECT CHARACTER*80 CHARS C First try to get the instrument name CALL UHDGST(ID,'INSTRUME',INST,ISTAT) IF(ISTAT.NE.0) RETURN C Check instrument can be handled IF(INST.EQ.'WFPC2') THEN C Get the other keywords we need CALL UHDGSI(ID,'DETECTOR',IDET,ISTAT) IF(ISTAT.NE.0) RETURN IF(IDET.LT.1 .OR. IDET.GT.4) THEN CALL UMSPUT('! Invalid detector number in header', : 1,0,ISTAT) ISTAT=1 RETURN ENDIF CALL UHDGSR(ID,'PHOTPLAM',PLAM,ISTAT) IF(ISTAT.NE.0) RETURN C We have to convert the wavelength from Angstroms to nm LAMBDA=PLAM/10.0 C Now we build the full "Trauger style" name for the coefficients C file IF(IDET.EQ.1) THEN COEFFS='drizzle$coeffs/pc1-trauger' ELSE IF(IDET.EQ.2) THEN COEFFS='drizzle$coeffs/wf2-trauger' ELSE IF(IDET.EQ.3) THEN COEFFS='drizzle$coeffs/wf3-trauger' ELSE IF(IDET.EQ.4) THEN COEFFS='drizzle$coeffs/wf4-trauger' ENDIF WRITE(CHARS, : '(''-Instrument: '',A8,'' Chip: '''// : ',I1,'' Wavelength: '',F7.2,'' (nm)'')') : INST,IDET,LAMBDA CALL UMSPUT(CHARS,1,0,ISTAT) ELSE IF(INST.EQ.'NICMOS') THEN C Get the other keywords we need CALL UHDGSI(ID,'CAMERA',IDET,ISTAT) IF(ISTAT.NE.0) RETURN IF(IDET.LT.1 .OR. IDET.GT.3) THEN CALL UMSPUT('! Invalid camera number in header', : 1,0,ISTAT) ISTAT=1 RETURN ELSE WRITE(COEFFS,'(''drizzle$coeffs/nic-'',I1)') IDET ENDIF WRITE(CHARS, : '(''-Instrument: NICMOS Camera: '',i1)') IDET CALL UMSPUT(CHARS,1,0,ISTAT) ELSE IF(INST.EQ.'STIS') THEN C Get the other keywords we need CALL UHDGST(ID,'DETECTOR',DETECT,ISTAT) IF(DETECT.EQ.'CCD') THEN COEFFS='drizzle$coeffs/stis-ccd' ELSE IF(DETECT.EQ.'NUV-MAMA') THEN COEFFS='drizzle$coeffs/stis-nuv-mama' ELSE IF(DETECT.EQ.'FUV-MAMA') THEN COEFFS='drizzle$coeffs/stis-fuv-mama' ELSE CALL UMSPUT('! Invalid STIS detector name in header', : 1,0,ISTAT) ISTAT=1 RETURN ENDIF WRITE(CHARS, : '(''-Instrument: STIS Detector: '',A8)') DETECT CALL UMSPUT(CHARS,1,0,ISTAT) ELSE CALL UMSPUT( : '! We can currently only handle WFPC2, NICMOS or STIS headers', : 1,0,ISTAT) ISTAT=1 RETURN ENDIF RETURN END SUBROUTINE GOLDH(ID,USEWCS,IC,LASTPF,LASTSC, : LASTUN,LASTOX,LASTOY,ISTAT) C C This routine finds out details of the last image C to have been drizzled onto the image with identifier C ID and returns then for comparison with the new C ones. C C This version was moved from the source code of V1.2 of C drizzle. C C Small modifications to fit in with changed header scheme. C March 2001 C IMPLICIT NONE INTEGER IC,ISTAT,ID LOGICAL USEWCS CHARACTER*8 LASTUN CHARACTER*4 STEM CHARACTER*80 BUFFER REAL LASTPF,LASTSC,LASTOX,LASTOY C Before doing anything check for old-style headers and C warn the user CALL UHDGST(ID,'D01VER',BUFFER,ISTAT) IF(ISTAT.EQ.0) THEN CALL UMSPUT( : '! Warning - there is an old-style Drizzle header present,', : 1,0,ISTAT) CALL UMSPUT( : '! please consult the header for full processing history', : 1,0,ISTAT) ENDIF C First find out if there are old values C This is very ugly STEM='D000' DO IC=1,999 WRITE(STEM(2:4),'(I3)') IC IF(STEM(2:2).EQ.' ') STEM(2:2)='0' IF(STEM(3:3).EQ.' ') STEM(3:3)='0' CALL UHDGST(ID,STEM//'VER',BUFFER,ISTAT) C Check for a status code indicating that this header C item doesn't exist IF(ISTAT.EQ.40) GO TO 111 ENDDO 111 CONTINUE IC=IC-1 WRITE(STEM(2:4),'(I3)') IC IF(STEM(2:2).EQ.' ') STEM(2:2)='0' IF(STEM(3:3).EQ.' ') STEM(3:3)='0' CALL UHDGSR(ID,STEM//'PIXF',LASTPF,ISTAT) IF(ISTAT.NE.0) RETURN CALL UHDGST(ID,STEM//'OUUN',LASTUN,ISTAT) IF(ISTAT.NE.0) RETURN IF(.NOT.USEWCS) THEN CALL UHDGSR(ID,STEM//'SCAL',LASTSC,ISTAT) IF(ISTAT.NE.0) RETURN CALL UHDGSR(ID,STEM//'OUXC',LASTOX,ISTAT) IF(ISTAT.NE.0) RETURN CALL UHDGSR(ID,STEM//'OUYC',LASTOY,ISTAT) IF(ISTAT.NE.0) RETURN ENDIF ISTAT=0 RETURN END SUBROUTINE INMAT(P,Q,R,S) C C Invert a 2 by 2 double precision matrix C in place. It assumes that the matrix is not C singular. C IMPLICIT NONE DOUBLE PRECISION P,Q,R,S DOUBLE PRECISION A,B,C,D,DET DET=P*S-Q*R A=S/DET B=-Q/DET C=-R/DET D=P/DET P=A Q=B R=C S=D RETURN END SUBROUTINE XY2RD(X,Y,R,D,WCS) C C XY2RD - convert a pixel position to an equatorial (RA,DEC) C position assuming a TAN projection and WCS with the C standard 8 elements. C C It is based on the xy2rd task in STSDAS C C Output angles are in degrees C C Richard Hook, STScI, 11th August 1998 C IMPLICIT NONE DOUBLE PRECISION X,Y,R,D,WCS(8) DOUBLE PRECISION XI,ETA,RA0,DEC0 DOUBLE PRECISION PIBY PARAMETER (PIBY=3.141592653589793/180.0) C First convert pixel coordinates to tangent plane in radians XI=(WCS(5)*(X-WCS(1))+WCS(7)*(Y-WCS(3)))*PIBY ETA=(WCS(6)*(X-WCS(1))+WCS(8)*(Y-WCS(3)))*PIBY C And convert the reference point on the sky to radians RA0=WCS(2)*PIBY DEC0=WCS(4)*PIBY C Now go to equatorial from tangent plane R=ATAN2(XI, COS(DEC0)-ETA*SIN(DEC0)) + RA0 D=ATAN2(ETA*COS(DEC0)+SIN(DEC0), : SQRT((COS(DEC0)-ETA*SIN(DEC0))**2 + XI**2)) C Convert back to degrees and check the range R=R/PIBY D=D/PIBY IF(R.LT.0.0) R=R+360.0 RETURN END SUBROUTINE RD2XY(R,D,X,Y,WCS,ISTAT) C C RD2XY - convert an equatorial (RA,DEC) position to a pixel position C assuming a TAN projection and WCS with the C standard 8 elements. C C It is based on the rd2xy task in STSDAS C C Input angles are in degrees C C Richard Hook, STScI, 13th August 1998 C IMPLICIT NONE INTEGER ISTAT DOUBLE PRECISION X,Y,R,D,WCS(8),RA,DEC DOUBLE PRECISION XI,ETA,RA0,DEC0 DOUBLE PRECISION DET,CDINV(2,2),BOTTOM DOUBLE PRECISION PIBY PARAMETER (PIBY=3.141592653589793/180.0) C First invert the CD matrix DET=WCS(5)*WCS(8)-WCS(7)*WCS(6) IF(DET.EQ.0.0) THEN ISTAT=1 RETURN ENDIF CDINV(1,1)=WCS(8)/DET CDINV(1,2)=-WCS(7)/DET CDINV(2,1)=-WCS(6)/DET CDINV(2,2)=WCS(5)/DET C Translate from RA,Dec to X,Y RA0=WCS(2)*PIBY DEC0=WCS(4)*PIBY RA=R*PIBY DEC=D*PIBY BOTTOM=SIN(DEC)*SIN(DEC0)+COS(DEC)*COS(DEC0)*COS(RA-RA0) IF(BOTTOM.EQ.0) THEN ISTAT=1 RETURN ENDIF C Calculate tangent plane position and convert to degrees XI=COS(DEC)*SIN(RA-RA0)/BOTTOM/PIBY ETA=(SIN(DEC)*COS(DEC0)-COS(DEC)*SIN(DEC0)*COS(RA-RA0))/ : BOTTOM/PIBY C Convert back to pixels using the inverse of the CD matrix X=CDINV(1,1)*XI+CDINV(1,2)*ETA+WCS(1) Y=CDINV(2,1)*XI+CDINV(2,2)*ETA+WCS(3) ISTAT=0 RETURN END SUBROUTINE WCSLIN(WCSIN,WCSOUT, : XCEN,YCEN,COTY,CONUM, : XCO,YCO,XC,YC,XS,YS,XT,YT) C C WCSLIN - derive best linear transformation coefficients to map C pixel positions from one world coordinate system to another. C C Supplied: C C WCSIN - double precision 8 element array. The WCS of the input image. C This includes CRPIX1/2, CRVAL1/2 and a CD matrix. C C WCSOUT - double precision 8 element array. The WCS of the output image. C This includes CRPIX1/2, CRVAL1/2 and a CD matrix. C C XCEN,YCEN - reals, the centre of the chip, for reference C C COTY,CONUM - coefficients type code and number of coefficients C C XCO,YCO - real arrays, the non-linear distortion coefficients C C Returned: C C XC,YC,XS,YS,XT,YT - reals, linear transformation coefficients C C Richard Hook, 18th September 1998 C Modified for more general coefficients and changed calling C sequence, December 2000 C C Further modification to try to handle ACS images with big distortions C and offsets. January 2001. C IMPLICIT NONE C Input parameters DOUBLE PRECISION WCSIN(8),WCSOUT(8) INTEGER COTY,CONUM REAL XCO(CONUM),YCO(CONUM),XCEN,YCEN C Output parameters REAL XC,YC,XS,YS,XT,YT C Local variables DOUBLE PRECISION R,D DOUBLE PRECISION XIN(4),YIN(4),XOUT(4),YOUT(4) REAL EVALCU,EVAL4,EVAL5,XX,YY,PIBY REAL X(4),Y(4),XO(4),YO(4) INTEGER I,ISTAT PARAMETER (PIBY=180.0/3.14159265359) C Set up a square at the reference pixel of the input C image (WCSIN(1)=CRPIX1 and WCSIN(3)=CRPIX2) XIN(1)=WCSIN(1) XIN(2)=WCSIN(1) XIN(3)=WCSIN(1)+1.0D0 XIN(4)=WCSIN(1)+1.0D0 YIN(1)=WCSIN(3) YIN(2)=WCSIN(3)+1.0D0 YIN(3)=WCSIN(3)+1.0D0 YIN(4)=WCSIN(3) C Transform these points onto the sky and then back out again C using the target WCS DO I=1,4 CALL XY2RD(XIN(I),YIN(I),R,D,WCSIN) CALL RD2XY(R,D,XOUT(I),YOUT(I),WCSOUT,ISTAT) ENDDO C Now we apply the geometric distortion to the input points so that C the linear transformation which we derive is appropriate after C the distortion is corrected DO I=1,4 XX=SNGL(XIN(I)) YY=SNGL(YIN(I)) IF(COTY.EQ.3) THEN X(I)=EVALCU(XX-XCEN,YY-YCEN,XCO) Y(I)=EVALCU(XX-XCEN,YY-YCEN,YCO) ELSE IF(COTY.EQ.4) THEN X(I)=EVAL4(XX-XCEN,YY-YCEN,XCO) Y(I)=EVAL4(XX-XCEN,YY-YCEN,YCO) ELSE IF(COTY.EQ.5) THEN X(I)=EVAL5(XX-XCEN,YY-YCEN,XCO) Y(I)=EVAL5(XX-XCEN,YY-YCEN,YCO) ELSE IF(COTY.EQ.-3) THEN CALL RAD3(XX-XCEN,YY-YCEN,XCO,X(I),Y(I)) ELSE X(I)=XX-XCEN Y(I)=YY-YCEN ENDIF ENDDO C Now we have the inputs and outputs and can derive the linear C transform between them C This is now done in a general way using least squares DO I=1,4 XO(I)=SNGL(XOUT(I)) YO(I)=SNGL(YOUT(I)) ENDDO CALL FITLIN(XO,YO,X,Y,4,XT,YT,XC,YS,XS,YC,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to determine mapping from WCS', : 1,0,ISTAT) ISTAT=1 RETURN ENDIF C We change a sign here to fit in with convention later YS=-YS C And now the linear offset XT=SNGL(XOUT(1))-XC*X(1)+YS*Y(1) YT=SNGL(YOUT(1))-XS*X(1)-YC*Y(1) C Before returning these values check that the two axes C have the same scales and are perpendicular after the geometric C distortion has been included CCC XSS=SQRT(((XOUT(4)-XOUT(1))**2+(YOUT(4)-YOUT(1))**2) / CCC : ((XIN(4)-XIN(1))**2+(YIN(4)-YIN(1))**2)) CCC YSS=SQRT(((XOUT(2)-XOUT(1))**2+(YOUT(2)-YOUT(1))**2) / CCC : ((XIN(2)-XIN(1))**2+(YIN(2)-YIN(1))**2)) CCC XA=PIBY*ATAN2(YOUT(4)-YOUT(1),XOUT(4)-XOUT(1)) CCC YA=PIBY*ATAN2(YOUT(2)-YOUT(1),XOUT(2)-XOUT(1)) CCC IF(ABS(1.0-YSS/XSS).GT.1.0E-3) THEN CCC CALL UMSPUT( CCC : '! Warning X and Y scales from WCS differ by more than 0.1%', CCC : 1,0,ISTAT) CCC ENDIF CCC IF(ABS(90.0-YA+XA).GT.0.01) THEN CCC CALL UMSPUT( CCC : '! Warning X and Y axes from WCS skewed by more than 0.01deg', CCC : 1,0,ISTAT) CCC ENDIF RETURN END SUBROUTINE SETWCS(OUTSCL,ORIENT,CRPIX1,CRVAL1, : CRPIX2,CRVAL2,WCS) C C SETWCS - convert a scale and orientation, along with C reference pixel and position on sky, to a WCS C including a CD matrix. C C Note: OUTSCL is in arcsecs/pix, ORIENT and CRVAL1/2 C are in degrees. C C Richard Hook, 22nd September 1998 C IMPLICIT NONE C Supplied: REAL OUTSCL,CRPIX1,CRPIX2,ORIENT DOUBLE PRECISION CRVAL1,CRVAL2 C Returned: DOUBLE PRECISION WCS(8) C Local variables DOUBLE PRECISION PIBY,ROR,DSC PARAMETER (PIBY=3.14159265358979/180.0) C Convert to radians for orientation ROR=DBLE(ORIENT)*PIBY C and to degrees for scale DSC=DBLE(OUTSCL)/3600.0 C Set reference points WCS(1)=DBLE(CRPIX1) WCS(3)=DBLE(CRPIX2) WCS(2)=CRVAL1 WCS(4)=CRVAL2 C and CD matrix (no skew, equal X,Y scales) WCS(5)=-DSC*DCOS(ROR) WCS(6)=DSC*DSIN(ROR) WCS(7)=WCS(6) WCS(8)=-WCS(5) RETURN END SUBROUTINE COPY1D(IN,OUT,N) C C Copy a 1d real array from one place to another. C IMPLICIT NONE INTEGER N REAL IN(*),OUT(*) INTEGER I DO I=1,N OUT(I)=IN(I) ENDDO RETURN END SUBROUTINE LCORN(DNX,DNY,ONX,ONY, : XSH,YSH,ROT,SCALE,ALIGN,ROTFIR, : SECPAR,XSH2,YSH2,ROT2,XSCALE,YSCALE,SHFR2,ROTF2, : USEWCS,WCSIN,WCSOUT, : COTY,CONUM,XCO,YCO,XMIN,XMAX,YMIN,YMAX,ISTAT) C C LCORN - transform the range of input pixel coordinates covered C by an image onto the output and find the extreme values. C C Modified to also pass the secondary geometric parameters, C Richard Hook, STScI, July 2000 C C Added SECPAR parameter, Sept 2000 C C Added check of the border and just 5 pixel margin, October 2000 C C Added more general coefficient support, December 2000 C IMPLICIT NONE INTEGER DNX,DNY,ONX,ONY REAL XSH,YSH,ROT,SCALE LOGICAL ROTFIR,USEWCS,SECPAR CHARACTER*8 ALIGN DOUBLE PRECISION WCSIN(8),WCSOUT(8) INTEGER COTY,CONUM REAL XCO(CONUM),YCO(CONUM) INTEGER XMIN,XMAX,YMIN,YMAX,ISTAT,I REAL XIN(16),YIN(16) REAL XOUT(16),YOUT(16),XMA,YMA,XMI,YMI C Secondary geometrical parameters, added in V1.5 REAL XSH2,YSH2,ROT2,XSCALE,YSCALE CHARACTER*8 SHFR2 LOGICAL ROTF2 XIN(1)=1.0 XIN(2)=1.0 XIN(3)=REAL(DNX) XIN(4)=REAL(DNX) YIN(1)=1.0 YIN(2)=REAL(DNY) YIN(3)=REAL(DNY) YIN(4)=1.0 C Transform onto output coordinate system CALL DRIVAL(XIN,YIN,4,DNX,DNY,ONX,ONY, : XSH,YSH,ROT,SCALE,ALIGN,ROTFIR, : SECPAR,XSH2,YSH2,ROT2,XSCALE,YSCALE,SHFR2,ROTF2, : USEWCS,WCSIN,WCSOUT, : COTY,CONUM,XCO,YCO,XOUT,YOUT) C Calculate the extreme values XMA=MAX(XOUT(1),XOUT(2),XOUT(3),XOUT(4)) YMA=MAX(YOUT(1),YOUT(2),YOUT(3),YOUT(4)) XMI=MIN(XOUT(1),XOUT(2),XOUT(3),XOUT(4)) YMI=MIN(YOUT(1),YOUT(2),YOUT(3),YOUT(4)) C Now check a few more points on the edges DO I=1,4 XIN(I)=REAL(I)*REAL(DNX)/5.0 YIN(I)=1.0 XIN(I+4)=XIN(I) YIN(I+4)=REAL(DNY) XIN(I+8)=1.0 YIN(I+8)=REAL(I)*REAL(DNY)/5.0 XIN(I+12)=REAL(DNX) YIN(I+12)=YIN(I+8) ENDDO CALL DRIVAL(XIN,YIN,16,DNX,DNY,ONX,ONY, : XSH,YSH,ROT,SCALE,ALIGN,ROTFIR, : SECPAR,XSH2,YSH2,ROT2,XSCALE,YSCALE,SHFR2,ROTF2, : USEWCS,WCSIN,WCSOUT, : COTY,CONUM,XCO,YCO,XOUT,YOUT) DO I=1,16 IF(XOUT(I).GT.XMA) XMA=XOUT(I) IF(YOUT(I).GT.YMA) YMA=YOUT(I) IF(XOUT(I).LT.XMI) XMI=XOUT(I) IF(YOUT(I).LT.YMI) YMI=YOUT(I) ENDDO C Calculate limits allowing for a 5 pixel margin of error XMAX=NINT(XMA+5.0) XMIN=NINT(XMI-5.0) YMAX=NINT(YMA+5.0) YMIN=NINT(YMI-5.0) ISTAT=0 RETURN END SUBROUTINE DRIVAL(XIN,YIN,N,DNX,DNY,ONX,ONY, : XSH,YSH,ROT,SCALE,ALIGN,ROTFIR, : SECPAR,XSH2,YSH2,ROT2,XSCALE,YSCALE,SHFR2,ROTF2, : USEWCS,WCSIN,WCSOUT, : COTY,CONUM,XCO,YCO,XOUT,YOUT) C C DRIVAL - apply the standard Drizzle transformation C from input to output pixel coordinates. C This may optional include a polynomial C distortion solution or be specified using C input and output WCS. It replaces the C logic originally in the DOBOX routine within Drizzle. C C Richard Hook, ST-ECF, 5th October 1998 C C Modified to include the secondary geometric transformation, C Richard Hook, STScI, July 2000 C C Added SECPAR logical flag, Sept 2000 C C Added more general coefficients, December 2000 C Used COTY to define the transformation and remove NONLIN C Added EVAL4/5 for higher order options C IMPLICIT NONE INTEGER N,DNX,DNY,ONX,ONY,I REAL XIN(N),YIN(N),XOUT(N),YOUT(N) REAL XCEN,YCEN,XF,YF,XCORN,YCORN REAL XSH,YSH,ROT,SCALE INTEGER COTY,CONUM REAL XCO(CONUM),YCO(CONUM) REAL XS,XC,YS,YC REAL EVALCU,EVAL4,EVAL5,SINTH,COSTH,XOFF,YOFF,XT,YT,XP,YP LOGICAL ROTFIR,USEWCS,SECPAR CHARACTER*8 ALIGN DOUBLE PRECISION WCSIN(8),WCSOUT(8) C Secondary geometrical parameters, added in V1.5 REAL XSH2,YSH2,ROT2,XSCALE,YSCALE CHARACTER*8 SHFR2 LOGICAL ROTF2 REAL XUT,YUT,XC2,XS2,YC2,YS2,XT2,YT2,XP2,YP2 REAL SINTH2,COSTH2 C-- IF(ALIGN.EQ.'corner') THEN XCEN=FLOAT(DNX/2)+0.5 YCEN=FLOAT(DNY/2)+0.5 ELSE XCEN=FLOAT(DNX/2)+1.0 YCEN=FLOAT(DNY/2)+1.0 ENDIF C Calculate some numbers to simplify things later SINTH=SIN(ROT) COSTH=COS(ROT) C The shifts are in units of INPUT pixels. XOFF=XSH/SCALE YOFF=YSH/SCALE XF=1.0/SCALE YF=1.0/SCALE C Some more economies XS=XF*SINTH XC=XF*COSTH YS=YF*SINTH YC=YF*COSTH C Secondary ones IF(SECPAR) THEN SINTH2=SIN(ROT2) COSTH2=COS(ROT2) XS2=SINTH2 XC2=COSTH2 YS2=SINTH2 YC2=COSTH2 ENDIF IF(ALIGN.EQ.'corner') THEN XP=REAL(ONX/2)+0.5 YP=REAL(ONY/2)+0.5 ELSE XP=REAL(ONX/2)+1.0 YP=REAL(ONY/2)+1.0 ENDIF XT=XOFF+XP YT=YOFF+YP C Set the secondary ones IF(SECPAR) THEN XP2=XP YP2=YP XT2=XP2+XSH2 YT2=YP2+YSH2 ENDIF C Here is some new code to support the WCS option C first we work out a linear transformation from input to C out IF(USEWCS) THEN CALL WCSLIN(WCSIN,WCSOUT, : XCEN,YCEN,COTY,CONUM, : XCO,YCO,XC,YC,XS,YS,XT,YT) ROTFIR=.TRUE. ENDIF C We split the cases with and without secondary parameters to C avoid excessive calculations C First the case with secondary parameters IF(SECPAR) THEN DO I=1,N IF(COTY.EQ.3) THEN XCORN=EVALCU(XSCALE*(XIN(I)-XCEN), : YSCALE*(YIN(I)-YCEN),XCO) YCORN=EVALCU(XSCALE*(XIN(I)-XCEN), : YSCALE*(YIN(I)-YCEN),YCO) ELSE IF(COTY.EQ.4) THEN XCORN=EVAL4(XSCALE*(XIN(I)-XCEN), : YSCALE*(YIN(I)-YCEN),XCO) YCORN=EVAL4(XSCALE*(XIN(I)-XCEN), : YSCALE*(YIN(I)-YCEN),YCO) ELSE IF(COTY.EQ.5) THEN XCORN=EVAL5(XSCALE*(XIN(I)-XCEN), : YSCALE*(YIN(I)-YCEN),XCO) YCORN=EVAL5(XSCALE*(XIN(I)-XCEN), : YSCALE*(YIN(I)-YCEN),YCO) ELSE IF(COTY.EQ.-3) THEN CALL RAD3(XSCALE*(XIN(I)-XCEN), : YSCALE*(YIN(I)-YCEN),XCO, : XCORN,YCORN) ELSE XCORN=XSCALE*(XIN(I)-XCEN) YCORN=YSCALE*(YIN(I)-YCEN) ENDIF C Apply the linear transform C There are two ways this can be done - shift then C rotate or rotate then shift IF(ROTFIR) THEN XUT=XC*XCORN-YS*YCORN+XT-XP YUT=XS*XCORN+YC*YCORN+YT-YP ELSE XUT=XC*(XCORN+XSH)-YS*(YCORN+YSH) YUT=XS*(XCORN+XSH)+YC*(YCORN+YSH) ENDIF C Apply the secondary transform IF(ROTF2) THEN XOUT(I)=XC2*XUT-YS2*YUT+XT2 YOUT(I)=XS2*XUT+YC2*YUT+YT2 ELSE XOUT(I)=XC2*(XUT+XSH2)-YS2*(YUT+YSH2)+XP2 YOUT(I)=XS2*(XUT+XSH2)+YC2*(YUT+YSH2)+YP2 ENDIF ENDDO C and now without secondary parameters ELSE DO I=1,N IF(COTY.EQ.3) THEN XCORN=EVALCU(XIN(I)-XCEN,YIN(I)-YCEN,XCO) YCORN=EVALCU(XIN(I)-XCEN,YIN(I)-YCEN,YCO) ELSE IF(COTY.EQ.4) THEN XCORN=EVAL4(XIN(I)-XCEN,YIN(I)-YCEN,XCO) YCORN=EVAL4(XIN(I)-XCEN,YIN(I)-YCEN,YCO) ELSE IF(COTY.EQ.5) THEN XCORN=EVAL5(XIN(I)-XCEN,YIN(I)-YCEN,XCO) YCORN=EVAL5(XIN(I)-XCEN,YIN(I)-YCEN,YCO) ELSE IF(COTY.EQ.-3) THEN CALL RAD3(XIN(I)-XCEN,YIN(I)-YCEN,XCO,XCORN,YCORN) ELSE XCORN=XIN(I)-XCEN YCORN=YIN(I)-YCEN ENDIF C Apply the linear transform C There are two ways this can be done - shift then C rotate or rotate then shift IF(ROTFIR) THEN XOUT(I)=XC*XCORN-YS*YCORN+XT YOUT(I)=XS*XCORN+YC*YCORN+YT ELSE XOUT(I)=XC*(XCORN+XSH)-YS*(YCORN+YSH)+XP YOUT(I)=XS*(XCORN+XSH)+YC*(YCORN+YSH)+YP ENDIF ENDDO ENDIF RETURN END SUBROUTINE DERICU(X,Y,CO,DX,DY) C C Evaluate the derivatives of a cubic polynomial distortion with C respect to X and Y. This is mainly for inverting the distortion C using Newton-Raphson. C C Richard Hook, STScI, May 1999 C IMPLICIT NONE REAL X,Y,CO(*),DX,DY DX=CO(2)+ : 2.0*CO(4)*X+ : CO(5)*Y+ : 3.0*CO(7)*X*X+ : 2.0*CO(8)*X*Y+ : CO(9)*Y*Y DY=CO(3)+ : CO(5)*X+ : 2.0*CO(6)*Y+ : CO(8)*X*X+ : 2.0*CO(9)*X*Y+ : 3.0*CO(10)*Y*Y RETURN END SUBROUTINE INVECU(XOUT,YOUT,XCO,YCO,ERR,XIN,YIN) C C Invert a cubic distortion in 2d using Newton-Raphson C C Supplied: C C Xout,Yout - real - the output position from the geometric C distortion. C C Xco,Yco - real arrays - the cubic distortion coefficients C C Err - real - the accuracy required of the inversion (in X and Y) C C Returned: C C Xin, Yin - real - the position which, when distorted, ends up C at Xout,Yout to accuracy Err. C C This method is interative and relatively slow. C C Richard Hook, STScI, May 1999 C IMPLICIT NONE REAL XOUT,YOUT,XCO(*),YCO(*),XIN,YIN,ERR REAL DXX,DXY,DYX,DYY,XO,YO,X,Y,D REAL EVALCU C First guess X=2.0*XOUT-EVALCU(XOUT,YOUT,XCO) Y=2.0*YOUT-EVALCU(XOUT,YOUT,YCO) DO WHILE(.TRUE.) C Distort current guesses XO=EVALCU(X,Y,XCO) YO=EVALCU(X,Y,YCO) C Check for error criterion IF(ABS(XOUT-XO).LT.ERR .AND. : ABS(YOUT-YO).LT.ERR) GO TO 999 C Calculate derivatives - there are four of these CALL DERICU(X,Y,XCO,DXX,DXY) CALL DERICU(X,Y,YCO,DYX,DYY) C Work out the determinant D=DXX*DYY-DYX*DXY C Improve guess with Newton-Raphson X=X+((XOUT-XO)*DYY-(YOUT-YO)*DXY)/D Y=Y+((YOUT-YO)*DXX-(XOUT-XO)*DYX)/D ENDDO 999 CONTINUE XIN=X YIN=Y RETURN END REAL FUNCTION OVER(I,J,XMIN,XMAX,YMIN,YMAX) C C OVER - calculate overlap between an arbitrary rectangle, aligned C with the axes, and a pixel. C C This is a simplified version of the BOXER code. C C Richard Hook, 6th May 1998 C IMPLICIT NONE INTEGER I,J REAL XMIN,XMAX,YMIN,YMAX REAL DX,DY DX=MIN(XMAX,REAL(I)+0.5)-MAX(XMIN,REAL(I)-0.5) DY=MIN(YMAX,REAL(J)+0.5)-MAX(YMIN,REAL(J)-0.5) IF(DX.GT.0 .AND. DY.GT.0) THEN OVER=DX*DY ELSE OVER=0.0 ENDIF RETURN END SUBROUTINE COPY1I(IN,OUT,N) C C Copy a 1d integer array from one place to another. C IMPLICIT NONE INTEGER N INTEGER IN(*),OUT(*) INTEGER I DO I=1,N OUT(I)=IN(I) ENDDO RETURN END SUBROUTINE UPCON(NCON,II,JJ,OLDCON,NEWCON,DONE,LX,LY, : INTAB,NEN,MAXIM,MAXEN,UNIQID,ISTAT) C C Update the context image C C This routine is called in the heart of Drizzle when the C context image needs to be updated. C C October 2000 (from the EIS version) C Added "bitmask" style context (16 bit), January 2001 C Converted to 32bit integers C IMPLICIT NONE INTEGER ISTAT,K INTEGER MAXEN,MAXIM,NEN,ICON,NN,LX,LY INTEGER UNIQID,INTAB(MAXIM,MAXEN) INTEGER OLDCON,NEWCON INTEGER NEWMA(100),II,JJ INTEGER NCON(LX,LY) INTEGER DONE(LX,LY) LOGICAL MATCH CHARACTER*80 CHARS C Look up the current context value ICON=NCON(II,JJ) C If it is the same as the last one we don't need to C go further IF(ICON.EQ.OLDCON) THEN NCON(II,JJ)=NEWCON GO TO 88 ENDIF C Combine with the new one IF(ICON.EQ.0) THEN NN=0 NEWMA(1)=1 NEWMA(2)=UNIQID ELSE NN=INTAB(1,ICON) C Check for whether this image is already on this pixel DO K=3,NN+2 IF(UNIQID.EQ.INTAB(K,ICON)) GO TO 88 ENDDO C If not create the new context by adding the current image NEWMA(1)=NN+1 DO K=2,NN+1 NEWMA(K)=INTAB(K+1,ICON) ENDDO C Check for too many images at a given context IF(NN.GT.MAXIM-3) THEN CALL UMSPUT( : '! Too many images - context table overloaded', : 1,0,ISTAT) ISTAT=1 RETURN ENDIF NEWMA(NN+2)=UNIQID ENDIF C Before matching sort the context array IF(NN.GT.0) CALL CSORT(NN+1,NEWMA(2)) C See whether we have had this context before DO K=NEN,1,-1 IF(MATCH(NEWMA,INTAB(1,K),MAXIM)) THEN NCON(II,JJ)=K GO TO 88 ENDIF ENDDO C No context match found - make a new one NEN=NEN+1 C Check for full table IF(NEN.EQ.MAXEN) THEN CALL UMSPUT('! Context table full',1,0,ISTAT) ISTAT=1 RETURN ENDIF NCON(II,JJ)=NEN INTAB(1,NEN)=NEWMA(1) INTAB(2,NEN)=0 DO K=3,NN+3 INTAB(K,NEN)=NEWMA(K-1) ENDDO C Tell the user about this new context IF(NN.LE.4) THEN WRITE(CHARS,'(''--New context #'',I5, : '' | '',I4,'' images: '', : 5I7)') NEN,(NEWMA(K),K=1,NN+2) ELSE WRITE(CHARS,'(''--New context #'',I5, : '' | '',I4,'' images: '', : 5I7,''...'')') NEN,(NEWMA(K),K=1,6) ENDIF IF(NEWMA(1).EQ.1) CHARS(34:34)=' ' CALL UMSPUT(CHARS,1,0,ISTAT) 88 CONTINUE C Save the old values for quick comparison OLDCON=ICON NEWCON=NCON(II,JJ) C Lastly we update the counter IF(OLDCON.NE.NEWCON) THEN IF(OLDCON.GT.0) INTAB(2,OLDCON)=INTAB(2,OLDCON)-1 INTAB(2,NEWCON)=INTAB(2,NEWCON)+1 ENDIF C Note that we have been here DONE(II,JJ)=1 ISTAT=0 RETURN END SUBROUTINE GETIMR(ID,NX,NY,SUB,XMIN,XMAX,YMIN,YMAX, : BUFF,P,ISTAT) C C Read a subset of an already opened image into an already C allocated memory space. This version for REAL images. C C A one dimensional buffer is used and must be big enough. C C October 2000 C IMPLICIT NONE INTEGER ID,NX,NY,XMIN,XMAX,YMIN,YMAX,ISTAT REAL P(XMAX-XMIN+1,YMAX-YMIN+1) REAL BUFF(1) LOGICAL SUB INTEGER I C If we have a subset then read just the bit of interest IF(SUB) THEN DO I=YMIN,YMAX CALL UIGL2R(ID,I,BUFF,ISTAT) IF(ISTAT.NE.0) RETURN C Copy the part of the line into the subset image CALL COPY1D(BUFF(XMIN), : P(1,I-YMIN+1), : XMAX-XMIN+1) ENDDO ELSE C Read the whole image, line by line DO I=1,NY CALL UIGL2R(ID,I,P(1,I),ISTAT) IF(ISTAT.NE.0) RETURN ENDDO ENDIF ISTAT=0 RETURN END SUBROUTINE GETIMI(ID,NX,NY,SUB,XMIN,XMAX,YMIN,YMAX, : BUFF,P,ISTAT) C C Read a subset of an already opened image into an already C allocated memory space. This version for INTEGER images. C C A one dimensional buffer is used and must be big enough. C C October 2000 C IMPLICIT NONE INTEGER ID,NX,NY,XMIN,XMAX,YMIN,YMAX,ISTAT INTEGER P(XMAX-XMIN+1,YMAX-YMIN+1) INTEGER BUFF(1) LOGICAL SUB INTEGER I C If we have a subset then read just the bit of interest IF(SUB) THEN DO I=YMIN,YMAX CALL UIGL2I(ID,I,BUFF,ISTAT) IF(ISTAT.NE.0) RETURN C Copy the part of the line into the subset image CALL COPY1I(BUFF(XMIN), : P(1,I-YMIN+1), : XMAX-XMIN+1) ENDDO ELSE C Read the whole image, line by line DO I=1,NY CALL UIGL2I(ID,I,P(1,I),ISTAT) IF(ISTAT.NE.0) RETURN ENDDO ENDIF ISTAT=0 RETURN END SUBROUTINE PUTIMR(ID,NX,NY,SUB,XMIN,XMAX,YMIN,YMAX, : UPDATE,OUTCPS,FS,BUFF,P,ISTAT) C C Write out an output image, possibly a subset. C C REAL version - with image scaling if required. C C October 2000 C IMPLICIT NONE INTEGER ID,NX,NY,XMIN,XMAX,YMIN,YMAX,ISTAT REAL P(XMAX-XMIN+1,YMAX-YMIN+1) REAL BUFF(1),FS LOGICAL SUB,UPDATE,OUTCPS INTEGER I C First the case of a subset image IF(SUB) THEN C First we read the ranges of Y for which the drizzle C has had no effect. This is not needed in the CPS C case except the first time around. IF(.NOT.UPDATE .OR. .NOT.OUTCPS) THEN IF(YMIN.GT.1) THEN DO I=1,YMIN-1 IF(UPDATE) THEN CALL UIGL2R(ID,I,BUFF,ISTAT) IF(FS.NE.1.0) THEN CALL MULC(BUFF,NX,1,FS) ENDIF ELSE CALL SETIM(BUFF,NX,1,0.0) ENDIF CALL UIPL2R(ID,I,BUFF,ISTAT) ENDDO ENDIF IF(YMAX.LT.NY) THEN DO I=YMAX+1,NY IF(UPDATE) THEN CALL UIGL2R(ID,I,BUFF,ISTAT) IF(FS.NE.1.0) THEN CALL MULC(BUFF,NX,1,FS) ENDIF ELSE CALL SETIM(BUFF,NX,1,0.0) ENDIF CALL UIPL2R(ID,I,BUFF,ISTAT) ENDDO ENDIF ENDIF C Now we read in the range of lines which we have already C processed, scale them and copy in the appropriate section DO I=YMIN,YMAX IF(UPDATE) THEN CALL UIGL2R(ID,I,BUFF,ISTAT) IF(FS.NE.1.0) THEN CALL MULC(BUFF,NX,1,FS) ENDIF ELSE CALL SETIM(BUFF,NX,1,0.0) ENDIF CALL COPY1D(P(1,I-YMIN+1),BUFF(XMIN),XMAX-XMIN+1) CALL UIPL2R(ID,I,BUFF,ISTAT) ENDDO IF(ISTAT.NE.0) RETURN C and now the simple case of the whole image ELSE DO I=1,NY CALL UIPL2R(ID,I,P(1,I),ISTAT) IF(ISTAT.NE.0) RETURN ENDDO ENDIF ISTAT=0 RETURN END SUBROUTINE PUTIMI(ID,NX,NY,SUB,XMIN,XMAX,YMIN,YMAX, : UPDATE,OUTCPS,FS,BUFF,P,ISTAT) C C Write out an output image, possibly a subset. C C INTEGER version. C C Values of FS other than 1.0 result in an error exit. C C October 2000 C IMPLICIT NONE INTEGER ID,NX,NY,XMIN,XMAX,YMIN,YMAX,ISTAT INTEGER P(XMAX-XMIN+1,YMAX-YMIN+1) INTEGER BUFF(1) REAL FS LOGICAL SUB,UPDATE,OUTCPS INTEGER I C Check the FS value - if not 1.0 then take an error exit C (this is only here for compatibility with the real version) IF(FS.NE.1.0) THEN ISTAT=2 RETURN ENDIF C Handle subsets first IF(SUB) THEN C First we read the ranges of Y for which the drizzle C has had no effect. This is not needed in the CPS C case except the first time around. IF(.NOT.UPDATE .OR. .NOT.OUTCPS) THEN IF(YMIN.GT.1) THEN DO I=1,YMIN-1 IF(UPDATE) THEN CALL UIGL2I(ID,I,BUFF,ISTAT) ELSE CALL SETIMI(BUFF,NX,1,0.0) ENDIF CALL UIPL2I(ID,I,BUFF,ISTAT) ENDDO ENDIF IF(YMAX.LT.NY) THEN DO I=YMAX+1,NY IF(UPDATE) THEN CALL UIGL2I(ID,I,BUFF,ISTAT) ELSE CALL SETIMI(BUFF,NX,1,0.0) ENDIF CALL UIPL2I(ID,I,BUFF,ISTAT) ENDDO ENDIF ENDIF C Now we read in the range of lines which we have already C processed, scale them and copy in the appropriate section DO I=YMIN,YMAX IF(UPDATE) THEN CALL UIGL2I(ID,I,BUFF,ISTAT) ELSE CALL SETIMI(BUFF,NX,1,0.0) ENDIF CALL COPY1I(P(1,I-YMIN+1),BUFF(XMIN),XMAX-XMIN+1) CALL UIPL2I(ID,I,BUFF,ISTAT) ENDDO IF(ISTAT.NE.0) RETURN C and now the simple case of the whole image ELSE DO I=1,NY CALL UIPL2I(ID,I,P(1,I),ISTAT) IF(ISTAT.NE.0) RETURN ENDDO ENDIF ISTAT=0 RETURN END SUBROUTINE LENSTR(STRING,I1,I2) C C Find the start and end of a string C IMPLICIT NONE CHARACTER*(*) STRING INTEGER I,I1,I2 LOGICAL IN IN=.FALSE. DO I=1,LEN(STRING) IF(STRING(I:I).NE.' ' .AND. : .NOT.IN) THEN I1=I IN=.TRUE. ENDIF IF(STRING(I:I).EQ.' ' .AND. : IN) THEN I2=I-1 GO TO 99 ENDIF ENDDO 99 CONTINUE RETURN END SUBROUTINE PTCOIN(ID,CONTAB,MAXHCN,INTAB,MAXIM,MAXEN,NEN,ISTAT) C C Write a context index table back to an open image header C C Richard Hook, November 1997 C Modified, 24th Nov 97 to included a "used" logical C vector so that only the relevant contexts are saved. C C Revised for new context format and other things, Richard Hook C March 1999 C C Added option to write to a file if there are more than a certain C number of entries. This added the CONTAB,MAXHCN parameters. C Richard Hook, July 1999 C IMPLICIT NONE INTEGER MAXIM,MAXEN,MAXHCN INTEGER ID,NEN INTEGER INTAB(MAXIM,MAXEN) INTEGER I,J,K,ISTAT,KWL CHARACTER*8 KEY CHARACTER*80 CHARS CHARACTER*80 CONTAB C Table of letters CHARACTER*1 LET(26) LET(1)='A' LET(2)='B' LET(3)='C' LET(4)='D' LET(5)='E' LET(6)='F' LET(7)='G' LET(8)='H' LET(9)='I' LET(10)='J' LET(11)='K' LET(12)='L' LET(13)='M' LET(14)='N' LET(15)='O' LET(16)='P' LET(17)='Q' LET(18)='R' LET(19)='S' LET(20)='T' LET(21)='U' LET(22)='V' LET(23)='W' LET(24)='X' LET(25)='Y' LET(26)='Z' C First we clear out all traces of an old header C Initialise header keyword matching for CON... CALL UHDOKL(ID,'CON*',.TRUE.,KWL,ISTAT) C Loop over entries until we get a bad status return DO WHILE(.TRUE.) C Get a keyword CALL UHDGNK(KWL,KEY,ISTAT) IF(ISTAT.NE.0) GO TO 77 C Check it is meaningful - to avoid too many mistakes READ(KEY(4:7),'(I4)',IOSTAT=ISTAT) IF(ISTAT.EQ.0) CALL UHDDSP(ID,KEY,ISTAT) ENDDO 77 CONTINUE C Close the keyword loop CALL UHDCKL(KWL,ISTAT) C Check whether we have lots of entries, if so write them C to a file rather than to the header and record the fact in the header C itself. Added 1/7/99 IF(NEN.GT.MAXHCN) THEN C Write the name of the file to the header CALL UHDAST(ID,'CONFILE',CONTAB, : 'File containing context text table',0,ISTAT) C Write the context information to that file CALL PTGLCO(CONTAB,INTAB,MAXIM,MAXEN,NEN,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to write context table file', : 1,0,ISTAT) ISTAT=1 RETURN ELSE CALL UMSPUT('-Writing context table file: '//CONTAB, : 1,0,ISTAT) ENDIF ELSE C Loop over the different entries DO I=1,NEN C Check for empty contexts and ignore them IF(INTAB(2,I).GT.0) THEN C Put together a name for the keyword WRITE(KEY,'(''CON'',I4)') I DO J=4,7 IF(KEY(J:J).EQ.' ') KEY(J:J)='0' ENDDO C Write the first line WRITE(CHARS,'(I3,I10,7I7)') : (INTAB(J,I),J=1,MIN(9,INTAB(1,I)+2)) C If this is the first item then label it IF(I.EQ.1) THEN CALL UHDAST(ID,KEY,CHARS,'Image context index', : 0,ISTAT) ELSE CALL UHDAST(ID,KEY,CHARS,' ', : 0,ISTAT) ENDIF C Now loop over the additional lines which may be C required IF(INTAB(1,I).GT.7) THEN DO J=2,INT((INTAB(1,I)+1)/9)+1 WRITE(CHARS,'(9I7)') : (INTAB(K,I),K=(J-1)*9+1,MIN(J*9,INTAB(1,I)+2)) C Delete any old keyword of the same name CALL UHDDSP(ID,KEY//LET(J-1),ISTAT) CALL UHDAST(ID,KEY//LET(J-1),CHARS,' ', : 0,ISTAT) IF(ISTAT.NE.0) RETURN ENDDO ENDIF ENDIF ENDDO CALL UMSPUT('-Writing context table to context image header', : 1,0,ISTAT) ENDIF ISTAT=0 RETURN END SUBROUTINE GTCOIN(ID,INTAB,MAXIM,MAXEN,NEN,ISTAT) C C Read a context index table from an open image header C C Richard Hook, November 1997 C Modified, Richard Hook, March 1999 C C Modified to use the file option if the CONFILE keyword is C present - Richard Hook, July 1999 C IMPLICIT NONE INTEGER MAXIM,MAXEN INTEGER ID,NEN INTEGER INTAB(MAXIM,MAXEN) INTEGER I,J,K,L,N,NR,NVAL,ISTAT,KWL CHARACTER*7 KEY CHARACTER*80 CHARS,CONTAB C Table of letters CHARACTER*1 LET(26) LET(1)='A' LET(2)='B' LET(3)='C' LET(4)='D' LET(5)='E' LET(6)='F' LET(7)='G' LET(8)='H' LET(9)='I' LET(10)='J' LET(11)='K' LET(12)='L' LET(13)='M' LET(14)='N' LET(15)='O' LET(16)='P' LET(17)='Q' LET(18)='R' LET(19)='S' LET(20)='T' LET(21)='U' LET(22)='V' LET(23)='W' LET(24)='X' LET(25)='Y' LET(26)='Z' C Initialise the array DO J=1,MAXEN INTAB(1,J)=1 DO I=2,MAXIM INTAB(I,J)=0 ENDDO ENDDO C First check for a CONFILE keyword. If this is found then the C context table is read from there and not the header CALL UHDGST(ID,'CONFILE',CONTAB,ISTAT) IF(ISTAT.EQ.0) THEN CALL GTGLCO(CONTAB,INTAB,MAXIM,MAXEN,NEN,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to read context from table file', : 1,0,ISTAT) ISTAT=1 RETURN ELSE CALL UMSPUT('-Read context table file: '//CONTAB, : 1,0,ISTAT) ENDIF ELSE C Initialise the total in case of error I=0 C Initialise header keyword matching for CON... CALL UHDOKL(ID,'CON*',.TRUE.,KWL,ISTAT) IF(ISTAT.NE.0) GO TO 99 C Loop over entries until we get a bad status return DO WHILE(.TRUE.) C Get a keyword CALL UHDGNK(KWL,KEY,ISTAT) IF(ISTAT.NE.0) GO TO 99 C Check keyword is valid READ(KEY(4:8),*,IOSTAT=ISTAT) I IF(ISTAT.EQ.0) THEN C Read the keyword from the header CALL UHDGST(ID,KEY,CHARS,ISTAT) IF(ISTAT.NE.0) GO TO 99 C First just extract the first number which is the number of C items READ(CHARS,*) NVAL IF(NVAL.LT.8) THEN C Case of one row only (less than 11 values) NR=NVAL+2 READ(CHARS,*) (INTAB(K,I),K=1,NR) ELSE C Case for more than one row C First read the first row READ(CHARS,*) (INTAB(K,I),K=1,9) C Then work out how many remain NR=NVAL-7 C and loop over them in groups of 9 (if possible) L=0 DO WHILE(NR.GT.0) N=MIN(9,NR) L=L+1 CALL UHDGST(ID,KEY//LET(L),CHARS,ISTAT) IF(ISTAT.NE.0) GO TO 99 READ(CHARS,*,IOSTAT=ISTAT) : (INTAB(K,I),K=L*9+1,L*9+N) NR=NR-N ENDDO ENDIF ENDIF ENDDO C If we get here something must be wrong ISTAT=1 RETURN C This is the normal exit point 99 CONTINUE C Close keyword list CALL UHDCKL(KWL,ISTAT) NEN=I ENDIF C Check that we have successfully read some values IF(NEN.GT.0) THEN ISTAT=0 ELSE ISTAT=1 ENDIF RETURN END SUBROUTINE SETIMI(A,NX,NY,V) C C Set a 2d integer image to a constant value C INTEGER NX,NY INTEGER A(NX,NY),V INTEGER I,J DO J=1,NY DO I=1,NX A(I,J)=V ENDDO ENDDO RETURN END LOGICAL FUNCTION MATCH(A,B,N) C C Match up a context against a context table C Note that after drizzle V0.40 (EIS) the context C table itself has an extra second column for the C counter. C IMPLICIT NONE INTEGER I,N INTEGER A(N),B(N),AN,BN AN=A(1) BN=B(1) IF(AN.NE.BN) THEN MATCH=.FALSE. RETURN ELSE MATCH=.TRUE. DO I=2,AN+1 IF(A(I).NE.B(I+1)) THEN MATCH=.FALSE. RETURN ENDIF ENDDO ENDIF RETURN END SUBROUTINE CSORT(N,ARR) C C This routine is modified from the Numerical Recipes one C to work on INTEGER arrays. C C It sorts an array of integers in place. C INTEGER N,M,NSTACK INTEGER ARR(N) PARAMETER (M=7,NSTACK=50) INTEGER I,IR,J,JSTACK,K,L,ISTACK(NSTACK) INTEGER A,TEMP jstack=0 l=1 ir=n 1 if(ir-l.lt.M)then do 12 j=l+1,ir a=arr(j) do 11 i=j-1,1,-1 if(arr(i).le.a)goto 2 arr(i+1)=arr(i) 11 continue i=0 2 arr(i+1)=a 12 continue if(jstack.eq.0)return ir=istack(jstack) l=istack(jstack-1) jstack=jstack-2 else k=(l+ir)/2 temp=arr(k) arr(k)=arr(l+1) arr(l+1)=temp if(arr(l+1).gt.arr(ir))then temp=arr(l+1) arr(l+1)=arr(ir) arr(ir)=temp endif if(arr(l).gt.arr(ir))then temp=arr(l) arr(l)=arr(ir) arr(ir)=temp endif if(arr(l+1).gt.arr(l))then temp=arr(l+1) arr(l+1)=arr(l) arr(l)=temp endif i=l+1 j=ir a=arr(l) 3 continue i=i+1 if(arr(i).lt.a)goto 3 4 continue j=j-1 if(arr(j).gt.a)goto 4 if(j.lt.i)goto 5 temp=arr(i) arr(i)=arr(j) arr(j)=temp goto 3 5 arr(l)=arr(j) arr(j)=a jstack=jstack+2 if(jstack.gt.NSTACK)pause 'NSTACK too small in sort' if(ir-i+1.ge.j-l)then istack(jstack)=ir istack(jstack-1)=i ir=j-1 else istack(jstack)=j-1 istack(jstack-1)=l l=i endif endif goto 1 END SUBROUTINE GTGLCO(CONTAB,INTAB,MAXIM,MAXEN,NEN,ISTAT) C C Read a context index table from global context file C C Richard Hook, November 1997 C Modified to handle missing context values, March 1999 C IMPLICIT NONE INTEGER MAXIM,MAXEN INTEGER NEN,CN INTEGER INTAB(MAXIM,MAXEN) INTEGER I,J,K,NVAL,ISTAT CHARACTER*80 CONTAB CHARACTER*1024 BUFFER C Try to open the context file OPEN(11,FILE=CONTAB,STATUS='OLD',IOSTAT=ISTAT) IF(ISTAT.NE.0) RETURN C Initialise the context table DO J=1,MAXEN INTAB(1,J)=1 DO I=2,MAXIM INTAB(I,J)=0 ENDDO ENDDO C Work through the lines skipping comments and blank ones DO WHILE(.TRUE.) READ(11,'(A)',END=99) BUFFER IF(BUFFER(1:1).NE.'#' .AND. BUFFER.NE.' ') THEN READ(BUFFER,*,IOSTAT=ISTAT) CN,NVAL IF(ISTAT.NE.0 .OR. NVAL.GT.99 .OR. CN.GT.MAXEN) GO TO 100 C Check for empty contexts IF(NVAL.GT.0) THEN READ(BUFFER,*,IOSTAT=ISTAT) CN,(INTAB(K,CN),K=1,NVAL+2) IF(ISTAT.NE.0) GO TO 100 ENDIF ENDIF ENDDO 99 CONTINUE ISTAT=0 100 CONTINUE CLOSE(11) NEN=CN RETURN END SUBROUTINE PTGLCO(CONTAB,INTAB,MAXIM,MAXEN,NEN,ISTAT) C C Write a context index table to global context file C C Richard Hook, November 1997 C Modified not to write empty contexts, March 1999 C IMPLICIT NONE INTEGER MAXIM,MAXEN INTEGER NEN INTEGER INTAB(MAXIM,MAXEN) INTEGER I,K,ISTAT CHARACTER*80 CONTAB C Open the global context file with write access - this will work C whether or not it exists OPEN(12,FILE=CONTAB,IOSTAT=ISTAT) IF(ISTAT.NE.0) RETURN C Write a header of sorts WRITE(12,'(''# GLOBAL context table'')') WRITE(12,'(''#'')') WRITE(12,'(''# Context | Nima | Npix | Unique image ids...'')') DO K=1,NEN IF(INTAB(2,K).NE.0) THEN WRITE(12,'(I7,I6,I10,3X,200I7)') : K,(INTAB(I,K),I=1,INTAB(1,K)+2) ENDIF ENDDO CLOSE(12) RETURN END SUBROUTINE GTUNID(NAME,ID,UNIQID,ISTAT) C C Get a unique image id (an integer) from the header of C an open image. This version is for EMMI and looks for C firstly a keyword EISID and, if that doesn't exist, C it uses the name of the image and takes the value C between the last _ and the first . in the filename. C IMPLICIT NONE INTEGER I,ISTAT,UNIQID,ID,IDOT,IUND CHARACTER*80 NAME C First try to get the keyword CALL UHDGSI(ID,'EISID',UNIQID,ISTAT) IF(ISTAT.EQ.0) RETURN C Then try to find it from the image name IUND=0 IDOT=0 DO I=1,LEN(NAME) IF(NAME(I:I).EQ.'_') IUND=I IF(NAME(I:I).EQ.'.') THEN IDOT=I GO TO 88 ENDIF ENDDO 88 CONTINUE IF(IUND.LT.IDOT-1 .AND. IUND.GT.0) : READ(NAME(IUND+1:IDOT-1),'(I)',IOSTAT=ISTAT) UNIQID RETURN END SUBROUTINE UCLGSD(PNAME,VAL,ISTAT) C C This is a routine to emulate the F77/VOS routine of the same C name which seems to be missing. It gets the value as a character C string from the CL (although the CL itself checks that it is a C valid number) and converts it to a double precision value via C a free-format read. C IMPLICIT NONE CHARACTER*80 PNAME CHARACTER*80 CHARS DOUBLE PRECISION VAL INTEGER ISTAT CALL UCLGST(PNAME,CHARS,ISTAT) READ(CHARS,*,IOSTAT=ISTAT) VAL RETURN END SUBROUTINE BFILL(LUN,BUFFER,ISTAT) C C Read free-format numbers from an open text file and C concatenate them into a character string buffer. C C This is mostly to facilitate free-format numerical C reads from text files. C C Richard Hook, ST-ECF, December 2000 C IMPLICIT NONE INTEGER LUN CHARACTER*(*) BUFFER CHARACTER*256 LINE INTEGER ISTAT INTEGER IS,I1,I2 IS=1 BUFFER=' ' DO WHILE(.TRUE.) CALL UFGLIN(LUN,LINE,ISTAT) IF(ISTAT.NE.0) GO TO 88 IF(LINE.NE.' ' .AND. LINE(1:1).NE.'#') THEN CALL ENDS(LINE,I1,I2) IF(IS+I2-I1.GT.LEN(BUFFER)) THEN ISTAT=1 RETURN ELSE BUFFER(IS:IS+I2-I1)=LINE(I1:I2) IS=IS+I2-I1+2 ENDIF ENDIF ENDDO 88 CONTINUE ISTAT=0 RETURN END SUBROUTINE ENDS(STRING,I1,I2) C C Find the start and end of a string C This differs from LENSTR in that internal whitespace C is retained rather than just the first non-white section C being extracted. C IMPLICIT NONE CHARACTER*(*) STRING INTEGER I,I1,I2 I1=1 I2=1 DO I=1,LEN(STRING) IF(STRING(I:I).NE.' ') THEN I1=I GO TO 88 ENDIF ENDDO 88 CONTINUE DO I=LEN(STRING),1,-1 IF(STRING(I:I).NE.' ') THEN I2=I GO TO 99 ENDIF ENDDO 99 CONTINUE RETURN END SUBROUTINE OPCON(NAME,ONX,ONY,ID,ODCO,UNIQID,ISTAT) C C Open a context image in bitmap, 3d form, check C whether there is room for the bit to be set and, C if not, recreate a bigger one. Return the ID C of the opened channel. C C It is expected that the file exists and isn't being C created for the first time. C C Richard Hook, ST-ECF, January 2001 C IMPLICIT NONE CHARACTER*80 NAME INTEGER ID,UNIQID,ONX,ONY,IDOLD,ISTAT,J,N INTEGER IB LOGICAL ODCO C Local variables INTEGER NDIMS,DIMS(7),DATTYP,I1,I2,NP CHARACTER*80 PNAME C Buffer array INTEGER MEMI(1) COMMON /MEM/MEMI C Assume we fail ODCO=.FALSE. C First try to open what we have CALL UIMOPN(NAME,2,ID,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to open current context image', : 1,0,ISTAT) ISTAT=1 RETURN ENDIF C Look at its size and shape CALL UIMGID(ID,DATTYP,NDIMS,DIMS,ISTAT) C Check it is the right size IF(NDIMS.NE.2 .AND. NDIMS.NE.3) THEN CALL UMSPUT( : '! Input image is not 2 or 3 dimensional', : 1,0,ISTAT) ISTAT=1 RETURN ENDIF IF(NDIMS.EQ.2) DIMS(3)=1 C Close the image CALL UIMCLO(ID,ISTAT) C Calculate the number of planes required NP=(UNIQID-1)/32+1 C Check whether we need a new one IF(NP.GT.DIMS(3)) THEN CALL UMSPUT('! Creating new, more spacious, context file', : 1,0,ISTAT) C Take a copy of the current one CALL UIMDEL('ConTemp',ISTAT) CALL UIMRNM(NAME,'ConTemp',ISTAT) IF(ISTAT.NE.0) RETURN C Open it again CALL UIMOPN('ConTemp',1,IDOLD,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to re-open copy of context image', : 1,0,ISTAT) ISTAT=1 RETURN ENDIF C Create the new one DIMS(3)=NP NDIMS=3 CALL UIMCRE(NAME,4,3,DIMS,ID,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to create new output context image', : 1,0,ISTAT) ODCO=.FALSE. ISTAT=1 RETURN ELSE CALL UMSPUT( : '-Created new output context image: '//NAME, : 1,0,ISTAT) C Write some header information which will otherwise be C neglected CALL UHDAST(ID,'CTYPE3','CHANNEL',' ',0,ISTAT) CALL UHDASD(ID,'CRVAL3',0.0D0,' ',0,ISTAT) CALL UHDASD(ID,'CRPIX3',0.0D0,' ',0,ISTAT) CALL UHDASD(ID,'CD1_3',0.0D0,' ',0,ISTAT) CALL UHDASD(ID,'CD2_3',0.0D0,' ',0,ISTAT) CALL UHDASD(ID,'CD3_3',1.0D0,' ',0,ISTAT) CALL UHDASD(ID,'CD3_2',0.0D0,' ',0,ISTAT) CALL UHDASD(ID,'CD3_1',0.0D0,' ',0,ISTAT) ODCO=.TRUE. ENDIF C Allocate a buffer CALL UDMGET(DIMS(1),4,IB,ISTAT) C Copy the old stuff over DO N=1,NP-1 DO J=1,DIMS(2) IF(NDIMS.EQ.2) THEN CALL UIGL2I(IDOLD,J,MEMI(IB),ISTAT) IF(ISTAT.NE.0) RETURN ELSE CALL UIGL3I(IDOLD,J,N,MEMI(IB),ISTAT) IF(ISTAT.NE.0) RETURN ENDIF CALL UIPL3I(ID,J,N,MEMI(IB),ISTAT) IF(ISTAT.NE.0) RETURN ENDDO ENDDO C Copy the header as well CALL UHDCPY(IDOLD,ID,ISTAT) IF(ISTAT.NE.0) RETURN C Write some header information which will otherwise be C neglected IF(NP.EQ.2) THEN CALL UHDAST(ID,'CTYPE3','CHANNEL',' ',0,ISTAT) CALL UHDASD(ID,'CRVAL3',0.0D0,' ',0,ISTAT) CALL UHDASD(ID,'CRPIX3',0.0D0,' ',0,ISTAT) CALL UHDASD(ID,'CD1_3',0.0D0,' ',0,ISTAT) CALL UHDASD(ID,'CD2_3',0.0D0,' ',0,ISTAT) CALL UHDASD(ID,'CD3_3',1.0D0,' ',0,ISTAT) CALL UHDASD(ID,'CD3_2',0.0D0,' ',0,ISTAT) CALL UHDASD(ID,'CD3_1',0.0D0,' ',0,ISTAT) ENDIF C Close the temporary copy CALL UIMCLO(IDOLD,ISTAT) IF(ISTAT.NE.0) RETURN C Initialise the new plane CALL SET1I(MEMI(IB),DIMS(1),0) DO J=1,DIMS(2) CALL UIPL3I(ID,J,NP,MEMI(IB),ISTAT) IF(ISTAT.NE.0) RETURN ENDDO C Free buffer CALL UDMFRE(IB,4,ISTAT) C Close the new copy CALL UIMCLO(ID,ISTAT) IF(ISTAT.NE.0) RETURN C If we have got here all must be well so we can delete the C old copy CALL UIMDEL('ConTemp',ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Warning, failed to delete copy of old context image', : 1,0,ISTAT) ENDIF ENDIF IF(NDIMS.EQ.2) THEN PNAME=NAME ELSE CALL LENSTR(NAME,I1,I2) WRITE(PNAME,'(A,''[*,*,'',I2,'']'')') NAME(I1:I2),NP IF(PNAME(I2+5:I2+5).EQ.' ') PNAME(I2+5:I2+5)='0' ENDIF CALL UIMOPN(PNAME,2,ID,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to open context image'//PNAME, : 1,0,ISTAT) ODCO=.FALSE. ISTAT=1 RETURN ELSE CALL UMSPUT('-Opening extant context image: '//PNAME, : 1,0,ISTAT) ODCO=.TRUE. ENDIF RETURN END SUBROUTINE SET1I(A,N,V) C C Set a 1d array to a value C IMPLICIT NONE INTEGER N,V,I INTEGER A(N) DO I=1,N A(I)=V ENDDO RETURN END SUBROUTINE FITLIN(XO,YO,X,Y,N,X0,Y0,A,B,C,D,ISTAT) C C Fit a linear transformation (with six parameters, equivalent C to the linear WCS) to two sets of points. C C This uses the standard least-squares linear equations method. C C Richard Hook, ST-ECF, January 31st 2001 C IMPLICIT NONE INTEGER ISTAT,N REAL XO(N),YO(N),X(N),Y(N) REAL X0,Y0,A,B,C,D DOUBLE PRECISION MAT(10,10),XORG,YORG,XOORG,YOORG REAL DET DOUBLE PRECISION SIGXOX,SIGXOY,SIGXO,SIGYOX,SIGYOY,SIGYO INTEGER I,J C Initialize the matrix (3x3) DO J=1,3 DO I=1,3 MAT(I,J)=0.0 ENDDO ENDDO C Also initialise the vectors SIGXOX=0.0D0 SIGXOY=0.0D0 SIGXO=0.0D0 SIGYOX=0.0D0 SIGYOY=0.0D0 SIGYO=0.0D0 C Take off an offset XORG=X(1) YORG=Y(1) XOORG=XO(1) YOORG=YO(1) C Setup the normal equations DO I=1,N MAT(1,1)=MAT(1,1)+(X(I)-XORG)**2 MAT(1,2)=MAT(1,2)+(X(I)-XORG)*(Y(I)-YORG) MAT(1,3)=MAT(1,3)+(X(I)-XORG) MAT(2,2)=MAT(2,2)+(Y(I)-YORG)**2 MAT(2,3)=MAT(2,3)+(Y(I)-YORG) SIGXOX=SIGXOX+(XO(I)-XOORG)*(X(I)-XORG) SIGXOY=SIGXOY+(XO(I)-XOORG)*(Y(I)-YORG) SIGXO=SIGXO+(XO(I)-XOORG) SIGYOX=SIGYOX+(YO(I)-YOORG)*(X(I)-XORG) SIGYOY=SIGYOY+(YO(I)-YOORG)*(Y(I)-YORG) SIGYO=SIGYO+(YO(I)-YOORG) ENDDO C Use symmetry (the matrix is diagonal) MAT(3,3)=DBLE(N) MAT(2,1)=MAT(1,2) MAT(3,1)=MAT(1,3) MAT(3,2)=MAT(2,3) C Invert the matrix (we check it isn't singular) CALL MATINV(MAT,3,DET) IF(DET.EQ.0.0) THEN CALL UMSPUT('! Linear transformation matrix is singular', : 1,0,ISTAT) ISTAT=1 RETURN ENDIF C Multiply the inverse by the vector A=SNGL(SIGXOX*MAT(1,1)+SIGXOY*MAT(1,2)+SIGXO*MAT(1,3)) B=SNGL(SIGXOX*MAT(2,1)+SIGXOY*MAT(2,2)+SIGXO*MAT(2,3)) X0=SNGL(SIGXOX*MAT(3,1)+SIGXOY*MAT(3,2)+SIGXO*MAT(3,3)) C=SNGL(SIGYOX*MAT(1,1)+SIGYOY*MAT(1,2)+SIGYO*MAT(1,3)) D=SNGL(SIGYOX*MAT(2,1)+SIGYOY*MAT(2,2)+SIGYO*MAT(2,3)) Y0=SNGL(SIGYOX*MAT(3,1)+SIGYOY*MAT(3,2)+SIGYO*MAT(3,3)) C Note that X0 and Y0 haven't been corrected for the offsets C Normally they are not used ISTAT=0 RETURN END C SUBROUTINE MATINV.F C C SOURCE C BEVINGTON, PAGES 302-303. C C PURPOSE C INVERT A SYMMETRIC MATRIX AND CALCULATE ITS DETERMINANT C C USAGE C CALL MATINV (ARRAY, NORDER, DET) C C DESCRIPTION OF PARAMETERS C ARRAY - INPUT MATRIX WHICH IS REPLACED BY ITS INVERSE C NORDER - DEGREE OF MATRIX (ORDER OF DETERMINANT) C DET - DETERMINANT OF INPUT MATRIX C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C COMMENT C DIMENSION STATEMENT VALID FOR NORDER UP TO 10 C SUBROUTINE MATINV (ARRAY,NORDER,DET) DOUBLE PRECISION ARRAY,AMAX,SAVE DIMENSION ARRAY(10,10),IK(10),JK(10) C 10 DET=1. 11 DO 100 K=1,NORDER C C FIND LARGEST ELEMENT ARRAY(I,J) IN REST OF MATRIX C AMAX=0. 21 DO 30 I=K,NORDER DO 30 J=K,NORDER 23 IF (DABS(AMAX)-DABS(ARRAY(I,J))) 24,24,30 24 AMAX=ARRAY(I,J) IK(K)=I JK(K)=J 30 CONTINUE C C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K,K) C 31 IF (AMAX) 41,32,41 32 DET=0. GOTO 140 41 I=IK(K) IF (I-K) 21,51,43 43 DO 50 J=1,NORDER SAVE=ARRAY(K,J) ARRAY(K,J)=ARRAY(I,J) 50 ARRAY(I,J)=-SAVE 51 J=JK(K) IF (J-K) 21,61,53 53 DO 60 I=1,NORDER SAVE=ARRAY(I,K) ARRAY(I,K)=ARRAY(I,J) 60 ARRAY(I,J)=-SAVE C C ACCUMULATE ELEMENTS OF INVERSE MATRIX C 61 DO 70 I=1,NORDER IF (I-K) 63,70,63 63 ARRAY(I,K)=-ARRAY(I,K)/AMAX 70 CONTINUE 71 DO 80 I=1,NORDER DO 80 J=1,NORDER IF (I-K) 74,80,74 74 IF (J-K) 75,80,75 75 ARRAY(I,J)=ARRAY(I,J)+ARRAY(I,K)*ARRAY(K,J) 80 CONTINUE 81 DO 90 J=1,NORDER IF (J-K) 83,90,83 83 ARRAY(K,J)=ARRAY(K,J)/AMAX 90 CONTINUE ARRAY(K,K)=1./AMAX 100 DET=DET*AMAX C C RESTORE ORDERING OF MATRIX C 101 DO 130 L=1,NORDER K=NORDER-L+1 J=IK(K) IF (J-K) 111,111,105 105 DO 110 I=1,NORDER SAVE=ARRAY(I,K) ARRAY(I,K)=-ARRAY(I,J) 110 ARRAY(I,J)=SAVE 111 I=JK(K) IF (I-K) 130,130,113 113 DO 120 J=1,NORDER SAVE=ARRAY(K,J) ARRAY(K,J)=-ARRAY(I,J) 120 ARRAY(I,J)=SAVE 130 CONTINUE 140 RETURN END SUBROUTINE UHEAD(ID,VERS,DATA,WEIGHT,EXPTIM,WTSCL,OUTDAT, : OUTCOU,CONTIM,COEFFS,SHFTUN,PFRACT, : LAM,SCALE,ROT,SHFTFR,ALIGN,KERNEL,EXPKEY, : FILSTR,INUN,OUTUN, : SECPAR,XSH2,YSH2,ROT2,XSCALE,YSCALE,SHFR2, : USEWCS,XSH,YSH,DNX,DNY,ONX,ONY) C C Update the header of the output image with all the C parameters of the current DRIZZLE run. C C Added Kernel value and context file name, October 2000 C C Added exposure time, January 2001 C C Corrected bug in CONTIM output, March 2001 C C Modified to be a little more intelligent and not to write C out irrelevant information, March 2001 C IMPLICIT NONE INTEGER ID,DNX,DNY,ONX,ONY,ISTAT,I CHARACTER*80 DATA,WEIGHT CHARACTER*80 OUTDAT,OUTCOU,CONTIM,COEFFS,SHFTUN CHARACTER*80 BUFFER,FILSTR CHARACTER*40 VERS CHARACTER*8 EXPKEY,INUN,OUTUN,ALIGN CHARACTER*4 STEM CHARACTER*8 SHFTFR,KERNEL REAL PFRACT,LAM,SCALE,ROT,EXPTIM,WTSCL REAL DROT,XSH,YSH,XSHI,YSHI,OFF LOGICAL USEWCS C Secondary geometrical parameters, added in V1.5 REAL XSH2,YSH2,ROT2,XSCALE,YSCALE CHARACTER*8 SHFR2 LOGICAL SECPAR STEM='D000' C First find out if there are old values C This is very ugly DO I=1,999 WRITE(STEM(2:4),'(I3)') I IF(STEM(2:2).EQ.' ') STEM(2:2)='0' IF(STEM(3:3).EQ.' ') STEM(3:3)='0' CALL UHDGST(ID,STEM//'VER',BUFFER,ISTAT) C Check for a status code indicating that this header C item doesn't exist IF(ISTAT.EQ.40) GO TO 111 ENDDO 111 CONTINUE C First update the number of images on this output IF(I.EQ.1) THEN CALL UHDASI(ID,'NDRIZIM',I, : 'Drizzle, number of images drizzled onto this output', : 0,ISTAT) ELSE CALL UHDPSI(ID,'NDRIZIM',I,ISTAT) ENDIF C Then the version information CALL UHDAST(ID,STEM//'VER',VERS, : 'Drizzle, task version',0,ISTAT) C Then the source of the geometric information IF(USEWCS) THEN CALL UHDAST(ID,STEM//'GEOM','Header WCS', : 'Drizzle, source of geometric information',0,ISTAT) ELSE CALL UHDAST(ID,STEM//'GEOM','User parameters', : 'Drizzle, source of geometric information',0,ISTAT) ENDIF C Now we continue to add the other items using the same C "stem" CALL UHDAST(ID,STEM//'DATA',DATA, : 'Drizzle, input data image',0,ISTAT) CALL UHDASR(ID,STEM//'DEXP',EXPTIM, : 'Drizzle, input image exposure time (s)',0,ISTAT) CALL UHDAST(ID,STEM//'OUDA',OUTDAT, : 'Drizzle, output data image',0,ISTAT) CALL UHDAST(ID,STEM//'OUWE',OUTCOU, : 'Drizzle, output weighting image',0,ISTAT) CALL UHDAST(ID,STEM//'OUCO',CONTIM, : 'Drizzle, output context image',0,ISTAT) CALL UHDAST(ID,STEM//'MASK',WEIGHT, : 'Drizzle, input weighting image', : 0,ISTAT) CALL UHDASR(ID,STEM//'WTSC',WTSCL, : 'Drizzle, weighting factor for input image',0,ISTAT) CALL UHDAST(ID,STEM//'KERN',KERNEL, : 'Drizzle, form of weight distribution kernel',0,ISTAT) CALL UHDASR(ID,STEM//'PIXF',PFRACT, : 'Drizzle, linear size of drop',0,ISTAT) CALL UHDAST(ID,STEM//'COEF',COEFFS, : 'Drizzle, coefficients file name ',0,ISTAT) CALL UHDASR(ID,STEM//'LAM',LAM, : 'Drizzle, wavelength applied for transformation (nm)', : 0,ISTAT) C Only put the next entries is we are NOT using WCS IF(.NOT.USEWCS) THEN CALL UHDASR(ID,STEM//'SCAL',SCALE, : 'Drizzle, scale (pixel size) of output image',0,ISTAT) C Convert the rotation angle back to degrees DROT=ROT/3.1415926535897*180.0 CALL UHDASR(ID,STEM//'ROT',DROT, : 'Drizzle, rotation angle, degrees anticlockwise',0,ISTAT) C Check the SCALE and units C The units are INPUT pixels on entry to this routine IF(SHFTUN.EQ.'output') THEN XSHI=XSH/SCALE YSHI=YSH/SCALE ELSE XSHI=XSH YSHI=YSH ENDIF CALL UHDASR(ID,STEM//'XSH',XSHI, : 'Drizzle, X shift applied',0,ISTAT) CALL UHDASR(ID,STEM//'YSH',YSHI, : 'Drizzle, Y shift applied',0,ISTAT) CALL UHDAST(ID,STEM//'SFTU',SHFTUN, : 'Drizzle, units used for shifts (output or input pixels)', : 0,ISTAT) CALL UHDAST(ID,STEM//'SFTF',SHFTFR, : 'Drizzle, frame in which shifts were applied', : 0,ISTAT) ENDIF CALL UHDAST(ID,STEM//'EXKY',EXPKEY, : 'Drizzle, exposure keyword name in input image',0,ISTAT) CALL UHDAST(ID,STEM//'INUN',INUN, : 'Drizzle, units of input image - counts or cps',0,ISTAT) CALL UHDAST(ID,STEM//'OUUN',OUTUN, : 'Drizzle, units of output image - counts or cps',0,ISTAT) CALL UHDAST(ID,STEM//'FVAL',FILSTR, : 'Drizzle, fill value for zero weight output pixels', : 0,ISTAT) IF(ALIGN.EQ.'corner') THEN OFF=0.5 ELSE OFF=1.0 ENDIF IF(.NOT.USEWCS) THEN CALL UHDASR(ID,STEM//'INXC',FLOAT(DNX/2)+OFF, : 'Drizzle, reference center of input image (X)',0,ISTAT) CALL UHDASR(ID,STEM//'INYC',FLOAT(DNY/2)+OFF, : 'Drizzle, reference center of input image (Y)',0,ISTAT) CALL UHDASR(ID,STEM//'OUXC',FLOAT(ONX/2)+OFF, : 'Drizzle, reference center of output image (X)',0,ISTAT) CALL UHDASR(ID,STEM//'OUYC',FLOAT(ONY/2)+OFF, : 'Drizzle, reference center of output image (Y)',0,ISTAT) ENDIF C If there are secondary parameters add these too IF(SECPAR) THEN CALL UHDASB(ID,STEM//'SECP',.TRUE., : 'Drizzle, there are secondary geometric parameters',0,ISTAT) CALL UHDASR(ID,STEM//'XSCL',XSCALE, : 'Drizzle, Secondary X scale applied',0,ISTAT) CALL UHDASR(ID,STEM//'YSCL',YSCALE, : 'Drizzle, Secondary Y scale applied',0,ISTAT) CALL UHDASR(ID,STEM//'XSH2',XSH2, : 'Drizzle, Secondary X shift applied',0,ISTAT) CALL UHDASR(ID,STEM//'YSH2',YSH2, : 'Drizzle, Secondary Y shift applied',0,ISTAT) C Convert the rotation angle back to degrees DROT=ROT2/3.1415926535897*180.0 CALL UHDASR(ID,STEM//'ROT2',DROT, : 'Drizzle, secondary rotation angle, degrees anticlockwise', : 0,ISTAT) CALL UHDAST(ID,STEM//'SFF2',SHFR2, : 'Drizzle, frame in which secondary shifts were applied', : 0,ISTAT) ELSE CALL UHDASB(ID,STEM//'SECP',.FALSE., : 'Drizzle, there are no secondary geometric parameters',0,ISTAT) ENDIF RETURN END SUBROUTINE DOBOX(DATA,WEI,NDAT,NCOU,NCON,DONE,DNX,DNY, : XMIN,XMAX,YMIN,YMAX,NOOVER, : KERNEL,XI,XO,YI,YO,XIB,XOB,YIB,YOB, : ONX,ONY,COTY,CONUM,XCO,YCO,WTSCL,ALIGN,INCPS,EXPIN, : PFRACT,SCALE,ROT,XSH,YSH,WCS,WCSOUT,ROTFIR, : SECPAR,XSH2,YSH2,ROT2,XSCALE,YSCALE,SHFR2,ROTF2, : CON,BITCON,INTAB,MAXIM,MAXEN,NEN,UNIQID, : UPDATE,IDD,USEWEI,USEWCS,IDW,ISTAT) C C This module does the actual mapping of input flux to output images using C "boxer", a code written by Bill Sparks for FOC geometric C distortion correction, rather than the "drizzling" approximation. C C This works by calculating the positions of the four corners of C a quadrilateral on the output grid corresponding to the corners C of the input pixel and then working out exactly how much of each C pixel in the output is covered, or not. C C In V1.6 this was simplified to use the DRIVAL routine and also C to include some limited multi-kernel support. C IMPLICIT NONE INTEGER DNX,DNY,ONX,ONY,IDD,IDW,NP INTEGER I,J,II,JJ,NMISS,ISTAT,NHIT,IS,NXI,NXA,NYI,NYA INTEGER XMIN,XMAX,YMIN,YMAX,LX,LY REAL DATA(DNX),WEI(DNX),W,EXPIN REAL NDAT(XMAX-XMIN+1,YMAX-YMIN+1) REAL NCOU(XMAX-XMIN+1,YMAX-YMIN+1) INTEGER COTY,CONUM REAL XCO(CONUM),YCO(CONUM) DOUBLE PRECISION WCS(8) DOUBLE PRECISION WCSOUT(8) REAL X,Y,XOUT(4),YOUT(4),VC,XSH,YSH,DOVER,DOW,DH,XX,YY REAL XIN(4),YIN(4) REAL D,ROT,XF,YF,SCALE,XXI,XXA,YYI,YYA REAL PFRACT,S2 REAL AC,WTSCL,JACO,DX,DY,PFO,PFO2,DD REAL XI(DNX,4),YI(DNX,4),XO(DNX,4),YO(DNX,4) REAL XIB(DNX),YIB(DNX),XOB(DNX),YOB(DNX) REAL R2,ES,EFAC,NSIG,XCEN,YCEN REAL OVER CHARACTER*80 CHARS CHARACTER*8 ALIGN,KERNEL LOGICAL ROTFIR,UPDATE,USEWEI,INCPS,USEWCS,NOOVER C Context related things LOGICAL CON LOGICAL BITCON INTEGER NCON(XMAX-XMIN+1,YMAX-YMIN+1) INTEGER DONE(XMAX-XMIN+1,YMAX-YMIN+1) INTEGER MAXEN,MAXIM,NEN,UNIQID INTEGER OLDCON,NEWCON,BV INTEGER INTAB(MAXIM,MAXEN) C Number of sigma to be included for gaussians PARAMETER (NSIG=2.5) C Secondary geometrical parameters, added in V1.5 LOGICAL SECPAR REAL XSH2,YSH2,ROT2,XSCALE,YSCALE CHARACTER*8 SHFR2 LOGICAL ROTF2 C-- C Some initial settings - note that the reference pixel C position is determined by the value of ALIGN NMISS=0 OLDCON=-1 C The bitmask - trimmed to the appropriate range NP=(UNIQID-1)/32+1 BV=2**(UNIQID-1-32*(NP-1)) C In the WCS case we can't use the scale to calculate the C Jacobian so we need to do it C C Note that we use the centre of the image rather than the C reference pixel as the reference here IF(USEWCS) THEN IF(ALIGN.EQ.'corner') THEN XCEN=FLOAT(DNX/2)+0.5 YCEN=FLOAT(DNY/2)+0.5 ELSE XCEN=FLOAT(DNX/2)+1.0 YCEN=FLOAT(DNY/2)+1.0 ENDIF XIN(1)=XCEN XIN(2)=XCEN XIN(3)=XCEN+1.0 XIN(4)=XCEN+1.0 YIN(1)=YCEN YIN(2)=YCEN+1.0 YIN(3)=YCEN+1.0 YIN(4)=YCEN CALL DRIVAL(XIN,YIN,4,DNX,DNY,ONX,ONY, : XSH,YSH,ROT,SCALE,ALIGN,ROTFIR, : SECPAR,XSH2,YSH2,ROT2,XSCALE,YSCALE,SHFR2,ROTF2, : USEWCS,WCS,WCSOUT, : COTY,CONUM,XCO,YCO,XOUT,YOUT) SCALE=SQRT(1.0/ABS(0.5*((XOUT(2)-XOUT(4))*(YOUT(1)-YOUT(3)) - : (XOUT(1)-XOUT(3))*(YOUT(2)-YOUT(4))))) ENDIF C Image subset size LX=XMAX-XMIN+1 LY=YMAX-YMIN+1 XF=1.0/SCALE YF=1.0/SCALE DH=0.5*PFRACT AC=1.0/(PFRACT*PFRACT) C Recalculate the area scaling factor S2=SCALE*SCALE C Offsets DX=REAL(XMIN-1) DY=REAL(YMIN-1) C Half pixfrac on output PFO=PFRACT/SCALE/2.0 PFO2=PFO**2 C Some Gaussian related numbers IF(KERNEL.EQ.'gaussian') THEN EFAC=2.3548**2*S2*AC/2.0 ES=EFAC/3.1415926 PFO=NSIG*PFRACT/2.3548/SCALE ENDIF C We skip all this if there is no overlap IF(.NOT.NOOVER) THEN C This is the outer loop over all the lines in the input image C Before we start we can fill the X arrays as they don't change C with Y X=0.0 IF(KERNEL.EQ.'square') THEN DO I=1,DNX X=X+1.0 XI(I,1)=X-DH XI(I,2)=X+DH XI(I,3)=X+DH XI(I,4)=X-DH ENDDO ELSE DO I=1,DNX X=X+1.0 XIB(I)=X ENDDO ENDIF C Loop over input lines Y=0.0 DO J=1,DNY Y=Y+1.0 C Read in one row of data CALL UIGL2R(IDD,J,DATA,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to read data image', : 1,0,ISTAT) RETURN ENDIF C If the input image is not in CPS we need to divide by C the exposure IF(.NOT.INCPS) CALL MULC(DATA,DNX,1,1.0/EXPIN) C If there is one, read in a line of the mask image IF(USEWEI) THEN CALL UIGL2R(IDW,J,WEI,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to read mask image', : 1,0,IS) RETURN ENDIF ENDIF C At this point we can handle the different kernels separately C First the cases where we just transform a single point rather C than four - every case except the "classic" square-pixel kernel IF(KERNEL.NE.'square') THEN C Fill up the position arrays ready to transform (only in Y, X is C already done) DO I=1,DNX YIB(I)=Y ENDDO C Transform onto the output grid CALL DRIVAL(XIB,YIB,DNX,DNX,DNY,ONX,ONY, : XSH,YSH,ROT,SCALE,ALIGN,ROTFIR, : SECPAR,XSH2,YSH2,ROT2,XSCALE,YSCALE,SHFR2,ROTF2, : USEWCS,WCS,WCSOUT, : COTY,CONUM,XCO,YCO,XOB,YOB) C Now we consider the different cases, first the "point" option C where a single pixel in the output image is affected IF(KERNEL.EQ.'point') THEN C Offset within the subset DO I=1,DNX II=NINT(XOB(I)-DX) JJ=NINT(YOB(I)-DY) C Check it is on the output image IF(II.GE.1 .AND. II.LE.LX .AND. : JJ.GE.1 .AND. JJ.LE.LY) THEN VC=NCOU(II,JJ) C Allow for stretching because of scale change D=DATA(I)*S2 C Scale the weighting mask by the scale factor C Note that we DON'T scale by the Jacobian as it hasn't been C calculated DOW=WEI(I)*WTSCL C If we are creating or modifying the context image we C do so here IF(CON .AND. DOW.GT.0.0) THEN IF(BITCON) THEN NCON(II,JJ)=IOR(NCON(II,JJ),BV) ELSE IF(DONE(II,JJ).EQ.0) THEN CALL UPCON( : NCON,II,JJ,OLDCON,NEWCON,DONE,LX,LY, : INTAB,NEN,MAXIM,MAXEN,UNIQID,ISTAT) ENDIF ENDIF ENDIF C Just a simple calculation without logical tests IF(VC.EQ.0.0) THEN NDAT(II,JJ)=D ELSE NDAT(II,JJ)= : (NDAT(II,JJ)*VC+DOW*D)/ : (VC+DOW) ENDIF NCOU(II,JJ)=VC+DOW ELSE NMISS=NMISS+1 ENDIF ENDDO C Next "tophat" - a circular kernel giving equal weight to C all points within a certain radius ELSE IF(KERNEL.EQ.'tophat') THEN DO I=1,DNX C Offset within the subset XX=XOB(I)-DX YY=YOB(I)-DY XXI=XOB(I)-DX-PFO XXA=XOB(I)-DX+PFO YYI=YOB(I)-DY-PFO YYA=YOB(I)-DY+PFO NXI=NINT(XXI) NXA=NINT(XXA) NYI=NINT(YYI) NYA=NINT(YYA) NHIT=0 C Allow for stretching because of scale change D=DATA(I)*S2 C Scale the weighting mask by the scale factor C and inversely by the Jacobian to ensure conservation C of weight in the output DOW=WEI(I)*WTSCL DD=DOW*D C Loop over output pixels which could be affected DO JJ=NYI,NYA DO II=NXI,NXA C Check it is on the output image IF(II.GE.1 .AND. II.LE.LX .AND. : JJ.GE.1 .AND. JJ.LE.LY) THEN C Radial distance R2=(XX-REAL(II))**2+(YY-REAL(JJ))**2 C Weight is one within the specified radius and zero outside C Note: weight isn't conserved in this case IF(R2.LE.PFO2) THEN C Count the hits NHIT=NHIT+1 VC=NCOU(II,JJ) C If we are creating or modifying the context image we C do so here IF(CON .AND. DOW.GT.0.0) THEN IF(BITCON) THEN NCON(II,JJ)=IOR(NCON(II,JJ),BV) ELSE IF(DONE(II,JJ).EQ.0) THEN CALL UPCON( : NCON,II,JJ,OLDCON,NEWCON,DONE,LX,LY, : INTAB,NEN,MAXIM,MAXEN,UNIQID,ISTAT) ENDIF ENDIF ENDIF C Just a simple calculation without logical tests IF(VC.EQ.0.0) THEN NDAT(II,JJ)=D ELSE NDAT(II,JJ)= : (NDAT(II,JJ)*VC+DD)/ : (VC+DOW) ENDIF NCOU(II,JJ)=VC+DOW ENDIF ENDIF ENDDO ENDDO C Count cases where the pixel is off the output image IF(NHIT.EQ.0) NMISS=NMISS+1 ENDDO C Next a gaussian weighting kernel with FWHM = pixfrac/scale ELSE IF(KERNEL.EQ.'gaussian') THEN DO I=1,DNX C Offset within the subset XX=XOB(I)-DX YY=YOB(I)-DY XXI=XOB(I)-DX-PFO XXA=XOB(I)-DX+PFO YYI=YOB(I)-DY-PFO YYA=YOB(I)-DY+PFO NXI=NINT(XXI) NXA=NINT(XXA) NYI=NINT(YYI) NYA=NINT(YYA) NHIT=0 C Allow for stretching because of scale change D=DATA(I)*S2 C Scale the weighting mask by the scale factor C and inversely by the Jacobian to ensure conservation C of weight in the output W=WEI(I)*WTSCL C Loop over output pixels which could be affected DO JJ=NYI,NYA DO II=NXI,NXA C Check it is on the output image IF(II.GE.1 .AND. II.LE.LX .AND. : JJ.GE.1 .AND. JJ.LE.LY) THEN C Radial distance R2=(XX-REAL(II))**2+(YY-REAL(JJ))**2 C Weight is a scaled gaussian function of radial distance DOVER=ES*EXP(-R2*EFAC) C Count the hits NHIT=NHIT+1 VC=NCOU(II,JJ) DOW=DOVER*W C If we are creating or modifying the context image we C do so here IF(CON .AND. DOW.GT.0.0) THEN IF(BITCON) THEN NCON(II,JJ)=IOR(NCON(II,JJ),BV) ELSE IF(DONE(II,JJ).EQ.0) THEN CALL UPCON( : NCON,II,JJ,OLDCON,NEWCON,DONE,LX,LY, : INTAB,NEN,MAXIM,MAXEN,UNIQID,ISTAT) ENDIF ENDIF ENDIF C Just a simple calculation without logical tests IF(VC.EQ.0.0) THEN NDAT(II,JJ)=D ELSE NDAT(II,JJ)= : (NDAT(II,JJ)*VC+DOW*D)/ : (VC+DOW) ENDIF NCOU(II,JJ)=VC+DOW ENDIF ENDDO ENDDO C Count cases where the pixel is off the output image IF(NHIT.EQ.0) NMISS=NMISS+1 ENDDO C The "turbo" option with a constant square pixfrac (as used in EIS Drizzle) ELSE IF(KERNEL.EQ.'turbo') THEN DO I=1,DNX C Offset within the subset XXI=XOB(I)-DX-PFO XXA=XOB(I)-DX+PFO YYI=YOB(I)-DY-PFO YYA=YOB(I)-DY+PFO NXI=NINT(XXI) NXA=NINT(XXA) NYI=NINT(YYI) NYA=NINT(YYA) NHIT=0 C Allow for stretching because of scale change D=DATA(I)*S2 C Scale the weighting mask by the scale factor C and inversely by the Jacobian to ensure conservation C of weight in the output W=WEI(I)*WTSCL C Loop over output pixels which could be affected DO JJ=NYI,NYA DO II=NXI,NXA C Check it is on the output image IF(II.GE.1 .AND. II.LE.LX .AND. : JJ.GE.1 .AND. JJ.LE.LY) THEN C Calculate the overlap using the simpler "aligned" box routine DOVER=OVER(II,JJ,XXI,XXA,YYI,YYA) IF(DOVER.GT.0.0) THEN C Correct for pixfrac area factor DOVER=DOVER*S2*AC C Count the hits NHIT=NHIT+1 VC=NCOU(II,JJ) DOW=DOVER*W C If we are creating or modifying the context image we C do so here IF(CON .AND. DOW.GT.0.0) THEN IF(BITCON) THEN NCON(II,JJ)=IOR(NCON(II,JJ),BV) ELSE IF(DONE(II,JJ).EQ.0) THEN CALL UPCON( : NCON,II,JJ,OLDCON,NEWCON,DONE,LX,LY, : INTAB,NEN,MAXIM,MAXEN,UNIQID,ISTAT) ENDIF ENDIF ENDIF C Just a simple calculation without logical tests IF(VC.EQ.0.0) THEN NDAT(II,JJ)=D ELSE NDAT(II,JJ)= : (NDAT(II,JJ)*VC+DOW*D)/ : (VC+DOW) ENDIF NCOU(II,JJ)=VC+DOW ENDIF ENDIF ENDDO ENDDO C Count cases where the pixel is off the output image IF(NHIT.EQ.0) NMISS=NMISS+1 ENDDO C End of the "single point" transform case statement ENDIF C Next the "classic" drizzle square kernel... C this is different because we have to transform all four corners C of the shrunken pixel ELSE C Fill up the position arrays ready to transform C Each is offset to the corner of the square - only in Y, C X is already done DO I=1,DNX YI(I,1)=Y+DH YI(I,2)=Y+DH YI(I,3)=Y-DH YI(I,4)=Y-DH ENDDO C Transform onto the output grid CALL DRIVAL(XI(1,1),YI(1,1),DNX,DNX,DNY,ONX,ONY, : XSH,YSH,ROT,SCALE,ALIGN,ROTFIR, : SECPAR,XSH2,YSH2,ROT2,XSCALE,YSCALE,SHFR2,ROTF2, : USEWCS,WCS,WCSOUT, : COTY,CONUM,XCO,YCO,XO(1,1),YO(1,1)) CALL DRIVAL(XI(1,2),YI(1,2),DNX,DNX,DNY,ONX,ONY, : XSH,YSH,ROT,SCALE,ALIGN,ROTFIR, : SECPAR,XSH2,YSH2,ROT2,XSCALE,YSCALE,SHFR2,ROTF2, : USEWCS,WCS,WCSOUT, : COTY,CONUM,XCO,YCO,XO(1,2),YO(1,2)) CALL DRIVAL(XI(1,3),YI(1,3),DNX,DNX,DNY,ONX,ONY, : XSH,YSH,ROT,SCALE,ALIGN,ROTFIR, : SECPAR,XSH2,YSH2,ROT2,XSCALE,YSCALE,SHFR2,ROTF2, : USEWCS,WCS,WCSOUT, : COTY,CONUM,XCO,YCO,XO(1,3),YO(1,3)) CALL DRIVAL(XI(1,4),YI(1,4),DNX,DNX,DNY,ONX,ONY, : XSH,YSH,ROT,SCALE,ALIGN,ROTFIR, : SECPAR,XSH2,YSH2,ROT2,XSCALE,YSCALE,SHFR2,ROTF2, : USEWCS,WCS,WCSOUT, : COTY,CONUM,XCO,YCO,XO(1,4),YO(1,4)) DO I=1,DNX C Offset within the subset XOUT(1)=XO(I,1)-DX XOUT(2)=XO(I,2)-DX XOUT(3)=XO(I,3)-DX XOUT(4)=XO(I,4)-DX YOUT(1)=YO(I,1)-DY YOUT(2)=YO(I,2)-DY YOUT(3)=YO(I,3)-DY YOUT(4)=YO(I,4)-DY C Work out the area of the quadrilateral on the output grid C Note that this expression expects the points to be in C clockwise order JACO=0.5*((XOUT(2)-XOUT(4))*(YOUT(1)-YOUT(3)) - : (XOUT(1)-XOUT(3))*(YOUT(2)-YOUT(4))) NHIT=0 C Allow for stretching because of scale change D=DATA(I)*S2 C Scale the weighting mask by the scale factor C and inversely by the Jacobian to ensure conservation C of weight in the output W=WEI(I)*WTSCL C Loop over output pixels which could be affected DO JJ=NINT(MIN(YOUT(1),YOUT(2),YOUT(3),YOUT(4))), : NINT(MAX(YOUT(1),YOUT(2),YOUT(3),YOUT(4))) DO II=NINT(MIN(XOUT(1),XOUT(2),XOUT(3),XOUT(4))), : NINT(MAX(XOUT(1),XOUT(2),XOUT(3),XOUT(4))) C Check it is on the output image IF(II.GE.1 .AND. II.LE.LX .AND. : JJ.GE.1 .AND. JJ.LE.LY) THEN C Call boxer to calculate overlap CALL BOXER(II,JJ,XOUT,YOUT,DOVER) IF(DOVER.GT.0.0) THEN C Re-normalise the area overlap using the Jacobian DOVER=DOVER/JACO C Count the hits NHIT=NHIT+1 VC=NCOU(II,JJ) DOW=DOVER*W C If we are creating or modifying the context image we C do so here IF(CON .AND. DOW.GT.0.0) THEN IF(BITCON) THEN NCON(II,JJ)=IOR(NCON(II,JJ),BV) ELSE IF(DONE(II,JJ).EQ.0) THEN CALL UPCON( : NCON,II,JJ,OLDCON,NEWCON,DONE,LX,LY, : INTAB,NEN,MAXIM,MAXEN,UNIQID,ISTAT) ENDIF ENDIF ENDIF C Just a simple calculation without logical tests IF(VC.EQ.0.0) THEN NDAT(II,JJ)=D ELSE NDAT(II,JJ)= : (NDAT(II,JJ)*VC+DOW*D)/ : (VC+DOW) ENDIF NCOU(II,JJ)=VC+DOW ENDIF ENDIF ENDDO ENDDO C Count cases where the pixel is off the output image IF(NHIT.EQ.0) NMISS=NMISS+1 ENDDO C End of the kernel "case" blocks ENDIF ENDDO ELSE NMISS=DNX*DNY ENDIF C Report number of points which were outside the output frame IF(NMISS.GT.0) THEN WRITE(CHARS,'(''! Warning, '',I7,'' points were '','// : '''outside the output image.'')') : NMISS CALL UMSPUT(CHARS,1,0,ISTAT) ENDIF C Set good status if we get this far ISTAT=0 RETURN END SUBROUTINE UPWCS(WCSIN,WCSOUT,DNX,DNY,ONX,ONY, : XSH,YSH,ROT,SCALE,ALIGN,ROTFIR, : SECPAR,XSH2,YSH2,ROT2,XSCALE,YSCALE,SHFR2,ROTF2, : USEWCS,COTY,CONUM,XCO,YCO) C C Update the WCS to include the drizzling transformations C C This is done by applying the transform to a unit square at C the reference pixel in the input whilst retaining the same C reference point on the sky. C C Richard Hook, ST-ECF, September 2000 C IMPLICIT NONE INTEGER DNX,DNY,ONX,ONY,ISTAT REAL XSH,YSH,ROT,SCALE CHARACTER*8 ALIGN,SHFR2 LOGICAL ROTFIR,SECPAR,ROTF2,USEWCS REAL XSH2,YSH2,ROT2,XSCALE,YSCALE INTEGER COTY,CONUM REAL XCO(CONUM),YCO(CONUM) DOUBLE PRECISION AM,BM,CM,DM,WCSIN(8),WCSOUT(8) DOUBLE PRECISION W5,W6,W7,W8,OFF REAL XIN(3),YIN(3),XOUT(3),YOUT(3) PARAMETER (OFF=0.1D0) C If we have the WCS already just return IF(USEWCS) RETURN C Set up a 1x1 box at the reference pixel (three sides only) XIN(1)=REAL(WCSIN(1)) YIN(1)=REAL(WCSIN(3)) XIN(2)=REAL(WCSIN(1)+OFF) YIN(2)=REAL(WCSIN(3)) XIN(3)=REAL(WCSIN(1)) YIN(3)=REAL(WCSIN(3)+OFF) C Transform CALL DRIVAL(XIN,YIN,3,DNX,DNY,ONX,ONY, : XSH,YSH,ROT,SCALE,ALIGN,ROTFIR, : SECPAR,XSH2,YSH2,ROT2,XSCALE,YSCALE,SHFR2,ROTF2, : USEWCS,WCSIN,WCSOUT, : COTY,CONUM,XCO,YCO,XOUT,YOUT) C We can immediately set the reference point on the sky WCSOUT(1)=DBLE(XOUT(1)) WCSOUT(3)=DBLE(YOUT(1)) WCSOUT(2)=WCSIN(2) WCSOUT(4)=WCSIN(4) C Now work out the effective CD matrix of the transformation AM=DBLE(XOUT(2)-XOUT(1))/OFF BM=DBLE(XOUT(3)-XOUT(1))/OFF CM=DBLE(YOUT(2)-YOUT(1))/OFF DM=DBLE(YOUT(3)-YOUT(1))/OFF C Check the determinant for singularity IF(AM*DM-BM*CM.EQ.0.0) THEN CALL UMSPUT('! Matrix is singular, cannot update WCS', : 1,0,ISTAT) ELSE C Invert the matrix CALL INMAT(AM,BM,CM,DM) C Correct the CD matrix W5=AM*WCSIN(5)+CM*WCSIN(7) W6=AM*WCSIN(6)+CM*WCSIN(8) W7=BM*WCSIN(5)+DM*WCSIN(7) W8=BM*WCSIN(6)+DM*WCSIN(8) WCSOUT(5)=W5 WCSOUT(6)=W6 WCSOUT(7)=W7 WCSOUT(8)=W8 ENDIF RETURN END SUBROUTINE BOXER (IS, JS, X, Y, DAREA) C C BOXER -- compute area of box overlap C C Calculate the area common to input clockwise polygon x(n), y(n) with C square (is, js) to (is+1, js+1). C This version is for a quadrilateral. C C W.B. Sparks STScI 2-June-1990. C Phil Hodge 20-Nov-1990 Change calling sequence; single precision. C Richard Hook ECF 24-Apr-1996 Change coordinate origin C so that centre of pixel has integer position C 03-Jan-2001 Removed accuracy check IMPLICIT NONE INTEGER IS, JS REAL X(*), Y(*) REAL DAREA C-- REAL PX(4), PY(4), SUM REAL SGAREA INTEGER I C Set up coords relative to unit square at origin C Note that the +0.5s were added when this code was C included in DRIZZLE DO I = 1, 4 PX(I) = X(I) - IS +0.5 PY(I) = Y(I) - JS +0.5 ENDDO C For each line in the polygon (or at this stage, input quadrilateral) C calculate the area common to the unit square (allow negative area for C subsequent `vector' addition of subareas). SUM = 0.0 DO I = 1, 3 SUM = SUM + SGAREA (PX(I), PY(I), PX(I+1), PY(I+1), IS, JS) ENDDO SUM = SUM + SGAREA (PX(4), PY(4), PX(1), PY(1), IS, JS) DAREA = SUM RETURN END REAL FUNCTION SGAREA (X1, Y1, X2, Y2, IS, JS) C C To calculate area under a line segment within unit square at origin. C This is used by BOXER C IMPLICIT NONE REAL X1, Y1, X2, Y2 INTEGER IS,JS REAL M, C, DX REAL XLO, XHI, YLO, YHI, XTOP LOGICAL NEGDX DX = X2 - X1 C Trap vertical line IF (DX .EQ. 0.0) THEN SGAREA = 0.0 GO TO 80 ENDIF C Order the two input points in x IF (X1 .LT. X2) THEN XLO = X1 XHI = X2 ELSE XLO = X2 XHI = X1 ENDIF C And determine the bounds ignoring y for now IF (XLO .GE. 1.0) THEN SGAREA = 0.0 GO TO 80 ENDIF IF (XHI .LE. 0.0) THEN SGAREA = 0.0 GO TO 80 ENDIF XLO = MAX (XLO, 0.0) XHI = MIN (XHI, 1.0) C Now look at y C basic info about the line y = mx + c NEGDX = (DX .LT. 0.0) M = (Y2 - Y1) / DX C = Y1 - M * X1 YLO = M * XLO + C YHI = M * XHI + C C Trap segment entirely below axis IF (YLO .LE. 0.0 .AND. YHI .LE. 0.0) THEN SGAREA = 0.0 GO TO 80 ENDIF C Adjust bounds if segment crosses axis (to exclude anything below axis) IF (YLO .LT. 0.0) THEN YLO = 0.0 XLO = -C/M ENDIF IF (YHI .LT. 0.0) THEN YHI = 0.0 XHI = -C/M ENDIF C There are four possibilities: both y below 1, both y above 1 C and one of each. IF (YLO .GE. 1.0 .AND. YHI .GE. 1.0) THEN C Line segment is entirely above square IF (NEGDX) THEN SGAREA = XLO - XHI ELSE SGAREA = XHI - XLO ENDIF GO TO 80 ENDIF IF (YLO .LE. 1.0 .AND. YHI .LE. 1.0) THEN C Segment is entirely within square IF (NEGDX) THEN SGAREA = 0.5 * (XLO-XHI) * (YHI+YLO) ELSE SGAREA = 0.5 * (XHI-XLO) * (YHI+YLO) END IF GO TO 80 ENDIF C otherwise it must cross the top of the square XTOP = (1.0 - C) / M IF (YLO .LT. 1.0) THEN IF (NEGDX) THEN SGAREA = -(0.5 * (XTOP-XLO) * (1.0+YLO) + XHI - XTOP) ELSE SGAREA = 0.5 * (XTOP-XLO) * (1.0+YLO) + XHI - XTOP END IF GO TO 80 ENDIF IF (NEGDX) THEN SGAREA = -(0.5 * (XHI-XTOP) * (1.0+YHI) + XTOP-XLO) ELSE SGAREA = 0.5 * (XHI-XTOP) * (1.0+YHI) + XTOP-XLO END IF 80 CONTINUE RETURN END SUBROUTINE GETPAR(GEOMOD,KERNEL,PFRACT,COEFFS,LAM, : SCALE,ROT,XSH,YSH,SHFTFR,SHFTUN,ALIGN, : XSCALE,YSCALE,XSH2,YSH2,ROT2,SHFR2, : EXPKEY,WTSTR,FILSTR,INUN,OUTUN, : DATA,WEIGHT,OUTDAT,OUTCOU,CONTIM,ONX,ONY, : OUTSCL,RACEN,DECCEN,ORIENT) C C Read all the parameter values needed. C C In this version they are obtained from the IRAF CL. C There is no error checking as these routines should C always succeed. C C This separate routine, October 2000 C C Modified for WCS option - STScI, November 2000 C IMPLICIT NONE C The parameters CHARACTER*8 KERNEL,GEOMOD,EXPKEY,SHFTFR,INUN,OUTUN,ALIGN CHARACTER*80 COEFFS,SHFTUN,WTSTR,FILSTR CHARACTER*80 DATA,WEIGHT,OUTDAT,OUTCOU,CONTIM REAL SCALE,ROT,PFRACT,XSH,YSH,LAM REAL XSH2,YSH2,ROT2,XSCALE,YSCALE DOUBLE PRECISION RACEN,DECCEN REAL ORIENT,OUTSCL CHARACTER*8 SHFR2 INTEGER ONX,ONY C Local variables INTEGER ISTAT C First the mode parameter - this is the WDRIZZLE change IF(GEOMOD.NE.'user') CALL UCLGST('geomode',GEOMOD,ISTAT) C Kernel shape and size CALL UCLGST('kernel',KERNEL,ISTAT) CALL UCLGSR('pixfrac',PFRACT,ISTAT) C The name of the geometric distortion coefficients file C and wavelength for the Trauger (WFPC2) case CALL UCLGST('coeffs',COEFFS,ISTAT) CALL UCLGSR('lambda',LAM,ISTAT) C Primary linear transformation parameters CALL UCLGSR('scale',SCALE,ISTAT) CALL UCLGSR('rot',ROT,ISTAT) CALL UCLGSR('xsh',XSH,ISTAT) CALL UCLGSR('ysh',YSH,ISTAT) CALL UCLGST('shft_fr',SHFTFR,ISTAT) CALL UCLGST('shft_un',SHFTUN,ISTAT) CALL UCLGST('align',ALIGN,ISTAT) C Secondary transformation parameters CALL UCLGSR('dr2gpar.xscale',XSCALE,ISTAT) CALL UCLGSR('dr2gpar.yscale',YSCALE,ISTAT) CALL UCLGSR('dr2gpar.xsh',XSH2,ISTAT) CALL UCLGSR('dr2gpar.ysh',YSH2,ISTAT) CALL UCLGSR('dr2gpar.rot',ROT2,ISTAT) CALL UCLGST('dr2gpar.shft_fr',SHFR2,ISTAT) C Exposure time, weighting and data scaling related keywords CALL UCLGST('expkey',EXPKEY,ISTAT) CALL UCLGST('wt_scl',WTSTR,ISTAT) CALL UCLGST('fillval',FILSTR,ISTAT) CALL UCLGST('in_un',INUN,ISTAT) CALL UCLGST('out_un',OUTUN,ISTAT) C Data array names CALL UCLGST('data',DATA,ISTAT) CALL UCLGST('in_mask',WEIGHT,ISTAT) CALL UCLGST('outdata',OUTDAT,ISTAT) CALL UCLGST('outweig',OUTCOU,ISTAT) CALL UCLGST('outcont',CONTIM,ISTAT) C Output array sizes CALL UCLGSI('outnx',ONX,ISTAT) CALL UCLGSI('outny',ONY,ISTAT) C WCS mode parameters of the output C - only in the case where GEOMOD is set appropriately IF(GEOMOD.EQ.'wcs') THEN CALL UCLGSR('outscl',OUTSCL,ISTAT) CALL UCLGSD('racen',RACEN,ISTAT) CALL UCLGSD('deccen',DECCEN,ISTAT) CALL UCLGSR('orient',ORIENT,ISTAT) ENDIF RETURN END SUBROUTINE ALLMEM(DNX,DNY,ONX,NSX,NSY,CON,BITCON, : PDATA,PWEI,PNDAT,PNCOU,PNCON,PCDON,PBUFF,PIBUFF, : PXI,PYI,PXO,PYO,PXIB,PYIB,PXOB,PYOB,ISTAT) C C Allocate the dynamic memory arrays needed for drizzle. C C This routine was added and the corresponding calls removed from C the main module, October 2000 C C Added BITCON flag for PCDON array - January 2001 C IMPLICIT NONE INTEGER DNX,DNY,ONX,NSX,NSY INTEGER PDATA,PWEI,PNDAT,PNCOU,PNCON,PCDON,PBUFF,PIBUFF INTEGER PXI,PYI,PXO,PYO INTEGER PXIB,PYIB,PXOB,PYOB,ISTAT LOGICAL CON LOGICAL BITCON C Set all the pointers to zero to start with - if they stay that C way we can assume the allocation has failed PDATA=0 PWEI=0 PNDAT=0 PNCOU=0 PNCON=0 PCDON=0 PBUFF=0 PIBUFF=0 PXI=0 PYI=0 PXO=0 PYO=0 PXIB=0 PYIB=0 PXOB=0 PYOB=0 C Allocate the memory arrays and return the pointers, check for C error status C Single line buffers for input data and weight images CALL UDMGET(DNX,6,PDATA,ISTAT) IF(ISTAT.NE.0) RETURN CALL UDMGET(DNX,6,PWEI,ISTAT) IF(ISTAT.NE.0) RETURN C Note that the next three, as they are subsets of the output C are of a different size CALL UDMGET(NSX*NSY,6,PNDAT,ISTAT) IF(ISTAT.NE.0) RETURN CALL UDMGET(NSX*NSY,6,PNCOU,ISTAT) IF(ISTAT.NE.0) RETURN C This one, the context image, is optional IF(CON) THEN CALL UDMGET(NSX*NSY,4,PNCON,ISTAT) IF(ISTAT.NE.0) RETURN IF(.NOT.BITCON) THEN CALL UDMGET(NSX*NSY,4,PCDON,ISTAT) IF(ISTAT.NE.0) RETURN ENDIF ENDIF C Scratch buffer space - the larger of the input and output IF(DNX.GT.ONX) THEN CALL UDMGET(DNX,6,PBUFF,ISTAT) CALL UDMGET(DNX,4,PIBUFF,ISTAT) ELSE CALL UDMGET(ONX,6,PBUFF,ISTAT) CALL UDMGET(ONX,4,PIBUFF,ISTAT) ENDIF IF(ISTAT.NE.0) RETURN CALL UDMGET(DNX*4,6,PXI,ISTAT) IF(ISTAT.NE.0) RETURN CALL UDMGET(DNX*4,6,PYI,ISTAT) IF(ISTAT.NE.0) RETURN CALL UDMGET(DNX*4,6,PXO,ISTAT) IF(ISTAT.NE.0) RETURN CALL UDMGET(DNX*4,6,PYO,ISTAT) IF(ISTAT.NE.0) RETURN CALL UDMGET(DNX,6,PXIB,ISTAT) IF(ISTAT.NE.0) RETURN CALL UDMGET(DNX,6,PYIB,ISTAT) IF(ISTAT.NE.0) RETURN CALL UDMGET(DNX,6,PXOB,ISTAT) IF(ISTAT.NE.0) RETURN CALL UDMGET(DNX,6,PYOB,ISTAT) IF(ISTAT.NE.0) RETURN C If we get here all is well ISTAT=0 RETURN END SUBROUTINE FREMEM(PDATA,PWEI,PNDAT,PNCOU,PNCON,PCDON,PBUFF,PIBUFF, : PXI,PYI,PXO,PYO,PXIB,PYIB,PXOB,PYOB) C C Free all allocated memory arrays. C C Successful allocations will have non-zero pointers. C C Added October 2000. C IMPLICIT NONE INTEGER PDATA,PWEI,PNDAT,PNCOU,PNCON,PCDON,PBUFF,PIBUFF INTEGER PXI,PYI,PXO,PYO INTEGER PXIB,PYIB,PXOB,PYOB,ISTAT IF(PDATA.NE.0) CALL UDMFRE(PDATA,6,ISTAT) IF(PWEI.NE.0) CALL UDMFRE(PWEI,6,ISTAT) IF(PNDAT.NE.0) CALL UDMFRE(PNDAT,6,ISTAT) IF(PNCOU.NE.0) CALL UDMFRE(PNCOU,6,ISTAT) IF(PNCON.NE.0) CALL UDMFRE(PNCON,4,ISTAT) IF(PCDON.NE.0) CALL UDMFRE(PCDON,4,ISTAT) IF(PIBUFF.NE.0) CALL UDMFRE(PIBUFF,4,ISTAT) IF(PXI.NE.0) CALL UDMFRE(PXI,6,ISTAT) IF(PYI.NE.0) CALL UDMFRE(PYI,6,ISTAT) IF(PXO.NE.0) CALL UDMFRE(PXO,6,ISTAT) IF(PYO.NE.0) CALL UDMFRE(PYO,6,ISTAT) IF(PXIB.NE.0) CALL UDMFRE(PXIB,6,ISTAT) IF(PYIB.NE.0) CALL UDMFRE(PYIB,6,ISTAT) IF(PXOB.NE.0) CALL UDMFRE(PXOB,6,ISTAT) IF(PYOB.NE.0) CALL UDMFRE(PYOB,6,ISTAT) RETURN END SUBROUTINE GETGEO(COEFFS,IDD,LAM,COTY,COMAX,CONUM, : XCO,YCO,ISTAT) C C Get the geometrical distortion information, either from C a file or the header. C C This new routine, October 2000 C Modified in drizzle 2.1 for more flexible coefficients C IMPLICIT NONE INTEGER IDD INTEGER COTY,COMAX,CONUM REAL LAM,XCO(COMAX),YCO(COMAX) CHARACTER*80 COEFFS INTEGER LUN,ISTAT,I LOGICAL NOCO C Interpret the coefficients file name and act appropriately NOCO=.FALSE. IF(COEFFS.EQ.' ') NOCO=.TRUE. C Now check for the special string "header" IF(COEFFS.EQ.'header') THEN C and try to read these from the input header CALL GTCOEF(IDD,COEFFS,LAM,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Cannot deduce image distortion information from header', : 1,0,ISTAT) ISTAT=1 RETURN ENDIF ENDIF C Set default values IF(NOCO) THEN DO I=1,COMAX XCO(I)=0.0 YCO(I)=0.0 ENDDO CONUM=1 COTY=0 ELSE C Open the coefficients file CALL UFOPEN(COEFFS,1,LUN,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Unable to open coefficients file',1,0, : ISTAT) ISTAT=1 RETURN ELSE CALL UMSPUT('-Opening coefficients file: '//COEFFS, : 1,0,ISTAT) ENDIF C Read in the coefficients from the file CALL GETCO(LUN,LAM,COTY,COMAX,CONUM,XCO,YCO,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Unable to read coefficients', : 1,0,ISTAT) CALL UFCLOS(LUN,ISTAT) ISTAT=1 RETURN ENDIF C Close the coefficients file CALL UFCLOS(LUN,ISTAT) ENDIF C Set a good status return ISTAT=0 RETURN END SUBROUTINE COPHED(IDI,IDO,COPALL,ISTAT) C C Copy over some important header items from the first C drizzled input image to the newly created output images. C C Richard Hook, March 2001 C IMPLICIT NONE INTEGER IDI,IDO,ISTAT LOGICAL COPALL CHARACTER*80 BUFFER C If we are to copy everything just do it directly IF(COPALL) THEN CALL UHDCPY(IDI,IDO,ISTAT) ELSE C Just copy a few things C First equinox and radesys - as strings CALL UHDGST(IDI,'EQUINOX',BUFFER,ISTAT) IF(ISTAT.EQ.0) : CALL UHDAST(IDO,'EQUINOX',BUFFER,' ',0,ISTAT) CALL UHDGST(IDI,'RADESYS',BUFFER,ISTAT) IF(ISTAT.EQ.0) : CALL UHDAST(IDO,'RADESYS',BUFFER,' ',0,ISTAT) C Now instrument related things CALL UHDGST(IDI,'INSTRUME',BUFFER,ISTAT) IF(ISTAT.EQ.0) : CALL UHDAST(IDO,'INSTRUME',BUFFER,' ',0,ISTAT) CALL UHDGST(IDI,'DETECTOR',BUFFER,ISTAT) IF(ISTAT.EQ.0) : CALL UHDAST(IDO,'DETECTOR',BUFFER,' ',0,ISTAT) CALL UHDGST(IDI,'CAMERA',BUFFER,ISTAT) IF(ISTAT.EQ.0) : CALL UHDAST(IDO,'CAMERA',BUFFER,' ',0,ISTAT) CALL UHDGST(IDI,'PHOTPLAM',BUFFER,ISTAT) IF(ISTAT.EQ.0) : CALL UHDAST(IDO,'PHOTPLAM',BUFFER,' ',0,ISTAT) ENDIF RETURN END