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.0 and BLOT 0.6, October 2000 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 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 of the form given in the C WFPC2 instrument handbook C IMPLICIT NONE REAL X,Y,CO(10) 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 SUBROUTINE GETCO(LUN,LAM,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 IMPLICIT NONE INTEGER LUN,ISTAT,I,J CHARACTER*80 LINE REAL LAM,MGF2,N REAL A,B,C REAL XCO(10),YCO(10) 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 ISTAT=0 GO TO 99 C Next "cubic" style - 20 coefficients, read in rows C of five ELSE IF(LINE(1:5).EQ.'cubic') THEN J=1 DO WHILE(J.LT.20) CALL UFGLIN(LUN,LINE,ISTAT) IF(LINE.NE.' ' .AND. LINE(1:1).NE.'#') THEN IF(J.LE.10) THEN READ(LINE,*,IOSTAT=ISTAT) (XCO(I),I=J,J+4) IF(ISTAT.NE.0) GO TO 99 ELSE READ(LINE,*,IOSTAT=ISTAT) (YCO(I-10),I=J,J+4) IF(ISTAT.NE.0) GO TO 99 ENDIF J=J+5 ENDIF ENDDO 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 UHDPST(ID,'CTYPE1',CTYPE1,ISTAT) CALL UHDPST(ID,'CTYPE2',CTYPE2,ISTAT) CALL UHDPSD(ID,'CRPIX1',WCS(1),ISTAT) CALL UHDPSD(ID,'CRVAL1',WCS(2),ISTAT) CALL UHDPSD(ID,'CRPIX2',WCS(3),ISTAT) CALL UHDPSD(ID,'CRVAL2',WCS(4),ISTAT) CALL UHDPSD(ID,'CD1_1',WCS(5),ISTAT) CALL UHDPSD(ID,'CD2_1',WCS(6),ISTAT) CALL UHDPSD(ID,'CD1_2',WCS(7),ISTAT) CALL UHDPSD(ID,'CD2_2',WCS(8),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,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 IMPLICIT NONE INTEGER IC,ISTAT,ID 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 UHDGSR(ID,STEM//'SCAL',LASTSC,ISTAT) IF(ISTAT.NE.0) RETURN CALL UHDGST(ID,STEM//'OUUN',LASTUN,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 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,NONLIN, : XCEN,YCEN,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 NONLIN - logical, whether there are non-linear distortion coefficients C C XCEN,YCEN - reals, the centre of the chip, for reference C C XCO,YCO - double precision 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 IMPLICIT NONE C Input parameters DOUBLE PRECISION WCSIN(8),WCSOUT(8) REAL XCO(*),YCO(*),XCEN,YCEN LOGICAL NONLIN 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,XX,YY INTEGER I,ISTAT C Set up a central square for reference XIN(1)=DBLE(XCEN) XIN(2)=DBLE(XCEN) XIN(3)=DBLE(XCEN+1.0) XIN(4)=DBLE(XCEN+1.0) YIN(1)=DBLE(YCEN) YIN(2)=DBLE(YCEN+1.0) YIN(3)=DBLE(YCEN+1.0) YIN(4)=DBLE(YCEN) C Transform these points into the WCS 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(NONLIN) THEN XIN(I)=EVALCU(XX-XCEN,YY-YCEN,XCO) YIN(I)=EVALCU(XX-XCEN,YY-YCEN,YCO) ELSE XIN(I)=XX-XCEN YIN(I)=YY-YCEN ENDIF ENDDO C Now we have the inputs and outputs and can derive the linear C transform between them C First we get the linear shift from the reference pixel XT=SNGL(XOUT(1))-XIN(1) YT=SNGL(YOUT(1))-YIN(1) C And now the scales (ie, the CD matrix) XC=SNGL(0.5*((XOUT(4)-XOUT(1))/(XIN(4)-XIN(1)) + : (XOUT(3)-XOUT(2))/(XIN(3)-XIN(2)))) YS=SNGL(-0.5*((XOUT(2)-XOUT(1))/(YIN(2)-YIN(1)) + : (XOUT(3)-XOUT(4))/(YIN(3)-YIN(4)))) XS=SNGL(0.5*((YOUT(4)-YOUT(1))/(XIN(4)-XIN(1)) + : (YOUT(3)-YOUT(2))/(XIN(3)-XIN(2)))) YC=SNGL(0.5*((YOUT(2)-YOUT(1))/(YIN(2)-YIN(1)) + : (YOUT(3)-YOUT(4))/(YIN(3)-YIN(4)))) 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, : 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 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) REAL XCO(10),YCO(10) 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, : 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, : 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, : 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 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,XCO(*),YCO(*) REAL XS,XC,YS,YC REAL EVALCU,SINTH,COSTH,XOFF,YOFF,XT,YT,XP,YP LOGICAL ROTFIR,USEWCS,NONLIN,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 If all the non-linear coefficients are zero we can C ignore that part of things NONLIN=.FALSE. DO I=1,10 IF(XCO(I).NE.0.0 .OR. YCO(I).NE.0.0) NONLIN=.TRUE. ENDDO 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,NONLIN, : XCEN,YCEN,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(NONLIN) THEN XCORN=EVALCU(XSCALE*(XIN(I)-XCEN), : YSCALE*(YIN(I)-YCEN),XCO) YCORN=EVALCU(XSCALE*(XIN(I)-XCEN), : YSCALE*(YIN(I)-YCEN),YCO) 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-XP 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(NONLIN) THEN XCORN=EVALCU(XIN(I)-XCEN,YIN(I)-YCEN,XCO) YCORN=EVALCU(XIN(I)-XCEN,YIN(I)-YCEN,YCO) 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 COPY1S(IN,OUT,N) C C Copy a 1d short integer array from one place to another. C IMPLICIT NONE INTEGER N INTEGER*2 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 IMPLICIT NONE INTEGER ISTAT,K INTEGER MAXEN,MAXIM,NEN,ICON,NN,LX,LY INTEGER UNIQID,INTAB(MAXIM,MAXEN) INTEGER*2 OLDCON,NEWCON INTEGER NEWMA(100),II,JJ INTEGER*2 NCON(LX,LY) INTEGER*2 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 Note that we have been here DONE(II,JJ)=1 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 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 GETIMS(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*2 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*2 P(XMAX-XMIN+1,YMAX-YMIN+1) INTEGER*2 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 UIGL2S(ID,I,BUFF,ISTAT) IF(ISTAT.NE.0) RETURN C Copy the part of the line into the subset image CALL COPY1S(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 UIGL2S(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 PUTIMS(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*2 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*2 P(XMAX-XMIN+1,YMAX-YMIN+1) INTEGER*2 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 UIGL2S(ID,I,BUFF,ISTAT) ELSE CALL SETIMS(BUFF,NX,1,0.0) ENDIF CALL UIPL2S(ID,I,BUFF,ISTAT) ENDDO ENDIF IF(YMAX.LT.NY) THEN DO I=YMAX+1,NY IF(UPDATE) THEN CALL UIGL2S(ID,I,BUFF,ISTAT) ELSE CALL SETIMS(BUFF,NX,1,0.0) ENDIF CALL UIPL2S(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 UIGL2S(ID,I,BUFF,ISTAT) ELSE CALL SETIMS(BUFF,NX,1,0.0) ENDIF CALL COPY1S(P(1,I-YMIN+1),BUFF(XMIN),XMAX-XMIN+1) CALL UIPL2S(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 UIPL2S(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 SETIMS(A,NX,NY,V) C C Set a 2d integer*2 image to a constant value C INTEGER NX,NY INTEGER*2 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*(*) 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*(*) PNAME CHARACTER*80 CHARS DOUBLE PRECISION VAL INTEGER ISTAT CALL UCLGST(PNAME,CHARS,ISTAT) READ(CHARS,*,IOSTAT=ISTAT) VAL RETURN END