SUBROUTINE DRIZZE C++ C C DRIZZLE.F Version EIS V0.88 Variable Pixel Linear Reconstruction C C Also widely known as "drizzling" C C This task implements a simple method for transforming one image C onto the pixel grid of another. The method calculates where the C input pixels map on the output grid and performs an optimal weighted C combination of the input and output. Another image records the weighting C at a given pixel. C C It is based largely on ideas of Andy Fruchter and was coded by Richard Hook. C (Contacts: rhook@eso.org & fruchter@stsci.edu) C C This version is specially for the ESO Imaging Survey (EIS) C C This code uses the F77/VOS for IRAF i/o. C It may be compiled and linked using a command like: C C xc -p stsdas (-z) drizzle.f drutil.f eisut.f x_drizzle.x libiraf77.a C C Note that the optional -z switch (which tells IRAF to create a static C executable which doesn't use the IRAF sharable library) may be C necessary to avoid a virtual memory use limitation of around 167MB in C IRAF versions up to 2.10.4. C C History: C C EIS V0.1 created from drizzle.f V1.097, 9th April 1997 C C EIS V0.11 - first fully (?) working version, 16th April 1997 C C EIS V0.2 - removed many parameters which were irrelavant for EIS C as well as tidying up the output headers. 21st May 1997 C C EIS V0.21 - added WEIGHT keyword option and better WCS output. C also cleaned up output keywords. 1st July 1997 C C EIS V0.22 - simplified DOBOX and added new code to support much C more efficient code for pixfrac=1.0 to avoid repeat C calculations of the same numbers at pixel corners. 28th July 97 C C EIS V0.23 - trying to fix FITS output problems. 1st Aug 97 C C EIS V0.24 - changed to getting WGTSCALE from weights image, 4th Aug 97 C (this was subsequetly removed) C C EIS V0.25 - improved output header, 7th Aug 97 C C EIS V0.26 - corrected header update bug, 3rd September 1997 C C EIS V0.27 - corrected problem with negative jacobian, 7th October 1997 C C EIS V0.30 - introduced "context" image and index table, 24th October 1997 C C EIS V0.31 - new system for "context" image, 31st October 1997 C C EIS V0.32 - correct a few bugs and tidy up output, 4th November 1997 C C EIS V0.33 - put data file names into context image header and also C write the correct WCS to the context image. C C EIS V0.34 - try to speed up the context lookup by remembering the last C value. Also modified the reading of astrometric coefficients C so that WCS is used (CD matrix etc) if there are no LDAC C polynomial coefficients. C C EIS V0.35 - attempting a more general context implementation based C on a global context table and unique image identifiers. C C EIS V0.36 - support for more contexts (20000). This version doesn't C store context table in header, just in text file. 24/11/97 C C EIS V0.37 - added defaulting for unique image ID and also C wrote list of files to data and weight headers, not C just to context. 26/11/97 C C EIS V0.40 - added context counter to context table and prepared for C next "release" for use in EIS pipeline. 27/11/97 C C EIS V0.41 - added a new routine to check whether a line of input C overlaps the output. If it doesn't then the line is C ignored. This is an attempt to speed things up. 29/11/97 C C EIS V0.42 - support for 6 digit unique images ids. 4/12/97 C C EIS V0.43 - timing information and tests added. 8/12/97 C C EIS V0.44 - header improvements, better double precision value C reading from the CL and better information in the C context table. 17/12/97 C C EIS V0.45 - much reduced the header information to avoid very C lengthly files. 13/5/98 C C EIS V0.50 - modified drizzle box shape to use different code for C increased speed. 14/5/98 C C EIS V0.51 - another modification with an even simpler kernel C C EIS V0.52 - a slight change to the kernel again to ensure a correct C Jacobian at the centre of the field. C C EIS V0.53 - added a new array to record whether the context has been C updated at a given point. This was to avoid a major speed C problem. 16/7/98 C C EIS V0.54 - minor bug fix. C C EIS V0.55 - changed CDnnnnn keywords to DDnnnnn to avoid name clash C with CD matrix. 27/8/98 C C EIS V0.56 - put back the PFRACT != 1.0 option for Thomas Erben. This C should have no effect on normal EIS operations. C C EIS V0.57 - added scaling of flux according to the Jacobian, 14/12/98 C C EIS V0.58 - added defaulting of WGTSCALE to 1.0, 18/12/98 C C EIS V0.60 - added support for different astrometry headers styles C 15/1/99 C C EIS V0.61 - added support for multiple projections (TAN & COE so far) C 19/1/99 C C EIS V0.62 - small modification so that context images will not be C created if contab=' ' even if a name is specified for C the output context images. 19/2/99 C C EIS V0.70 - modified so that the context table is no longer held C in an external text file but instead as part of the C context image file. As a result the contab parameter C has gone. 11/3/99. C C EIS V0.71 - small bugfixs when testing for EIS/WFI. 29/6/99. C C EIS V0.80 - added code to handle the case where the input and output C pixel grids are the same except for an integer pixel offset C for "stacking" of de-dithered EIS frames. 23/12/99. C C EIS V0.82 - added WGTSCALE=1.0 for output images. 21/2/00. C C-----EIS V1.2 release C C EIS V0.83 - added facilities for handling COE format input C in a probably rather inefficient way. 2/3/00. C C EIS V0.83B - larger number of contexts (100000) for Thomas Erben C 20th July 2000, STScI C C EIS V0.85 - more header information in output. Put context C maximum to practical maximum (32767) C 10th August 2000, ST-ECF C C EIS V0.86 - modified code to continue when unable to open a weight C image, 25th October 2000 C C EIS V0.87 - modified use of TIME routines to work correctly on HPs C and also corrected the way the Jacobian was used. C Richard Hook, 10th September 2001 C C EIS V0.88 - increased number of images to maximum of 1000 C Richard Hook, 3rd May 2002 C C-- IMPLICIT NONE C Iraf global storage (real in this program) REAL MEMR(1) !I INTEGER*2 MEMS(1) !I DOUBLE PRECISION MEMD(1) !I COMMON /MEM/MEMR !I EQUIVALENCE (MEMR,MEMD) !I EQUIVALENCE (MEMR,MEMS) !I C Local variables INTEGER ISTAT,I,J INTEGER NDIMS,DNX,DNY,ONX,ONY,DDIMS(7),ODIMS(7),DATTYP INTEGER IDD,IDW,IDND,IDNC,IDCO,IC INTEGER PDATA,PWEI,PNDAT,PNCOU,PNCON,PCDON INTEGER IITIME,STIME REAL XSH,YSH,WTSCL REAL SCALE,ROT,PFRACT,FS REAL EXPIN,EXPOUT,EXPCUR,FILVAL DOUBLE PRECISION WCS(8) CHARACTER*80 DATA,WEIGHT,OUTDAT,OUTCOU,CHARS CHARACTER*80 WTSTR,FILSTR CHARACTER*45 VERS CHARACTER*8 KEY CHARACTER*8 CTYPE1,CTYPE2 CHARACTER*8 EXPKEY,OUTUN,ALIGN CHARACTER*24 ICTIME EXTERNAL IITIME,ICTIME CHARACTER*3 PROJ C EIS specific stuff INTEGER MPIX,NC DOUBLE PRECISION CRVAL1,CRPIX1,CRVAL2,CRPIX2 DOUBLE PRECISION CDELT1,CDELT2,EXCO(100),EYCO(100) DOUBLE PRECISION RAREF,DECREF,XPOFF,YPOFF,RPP DOUBLE PRECISION RARED,DECRED DOUBLE PRECISION RAOFF,DECOFF,FLXSCL,OUTSC INTEGER PXIN,PYIN,PXTRA,PYTRA,PBUFF,PXR1,PXR2,PYR1,PYR2 INTEGER XMIN,XMAX,YMIN,YMAX,PSBUFF CHARACTER*80 CONTIM,CONTAB INTEGER MAXEN,MAXIM,NEN,MAXHCN,I1,I2 INTEGER UNIQID PARAMETER (MAXEN=32768) PARAMETER (MAXIM=1000) PARAMETER (MAXHCN=20) INTEGER INTAB(MAXIM,MAXEN) LOGICAL EIS LOGICAL CON C Logical flags LOGICAL VERBOSE LOGICAL NOCO LOGICAL USEWEI LOGICAL UPDATE LOGICAL GOTWCS LOGICAL EXPUPD LOGICAL WRIWEI LOGICAL EXPSCL LOGICAL FILL LOGICAL CPS LOGICAL HCOEFF LOGICAL ADATA LOGICAL AWEI LOGICAL ANDAT LOGICAL ANCOU LOGICAL ANCON LOGICAL ODD LOGICAL ODW LOGICAL ODND LOGICAL ODNC LOGICAL ODCO LOGICAL SUB LOGICAL NOOVER LOGICAL EIDONE C Constant DOUBLE PRECISION PIBY PARAMETER (PIBY=3.141592653589793/180.0) C-- Start of executable code C Get the start time STIME=IITIME() VERS='DRIZZLE Version EIS V0.88 (3rd May 2002)' C First announce the version and the starting time CALL UMSPUT('+ '//VERS,1,0,ISTAT) CALL UMSPUT('+ Starting at '//ICTIME(STIME),1,0,ISTAT) VERBOSE=.TRUE. C Also get fraction of area of input pixels to be used CALL UCLGSR('pixfrac',PFRACT,ISTAT) C Check for small values IF(PFRACT.LT.0.001) PFRACT=0.001 C Get the filling value for undefined pixels CALL UCLGST('fillval',FILSTR,ISTAT) C Get the projection type (new in V0.61) CALL UCLGST('proj',PROJ,ISTAT) C Get the switch value for whether to write out in CPS or C counts (latter is the default) CALL UCLGST('out_un',OUTUN,ISTAT) IF(OUTUN.EQ.'cps') THEN CPS=.TRUE. ELSE CPS=.FALSE. ENDIF C Check for meaningful values IF(FILSTR.NE.'INDEF'.AND.FILSTR.NE.'indef') THEN READ(FILSTR,*,IOSTAT=ISTAT) FILVAL IF(ISTAT.NE.0) THEN CALL UMSPUT('! Invalid filling value specified', : 1,0,ISTAT) GO TO 99 ENDIF FILL=.TRUE. ELSE FILL=.FALSE. ENDIF C Set default values XSH=0.0 YSH=0.0 SCALE=1.0 ROT=0.0 ALIGN='CENTER' C Get the keyword used for exposure time CALL UCLGST('expkey',EXPKEY,ISTAT) C In v0.85 we always have EXPUPD=.FALSE. EXPUPD=.FALSE. C See whether we are going to scale the images by exposure C time before (and after) combination CALL UCLGSB('exp_sc',EXPSCL,ISTAT) C Check for compatible requests IF(.NOT.EXPUPD.AND.(EXPSCL.OR.CPS)) THEN CALL UMSPUT('! No exposure time keyword has been given', : 1,0,ISTAT) GO TO 99 ENDIF C Get the name of the input data image CALL UCLGST('data',DATA,ISTAT) C Try to access the input data image CALL UIMOPN(DATA,1,IDD,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Unable to open data image', : 1,0,ISTAT) ODD=.FALSE. GO TO 99 ELSE IF(VERBOSE) CALL UMSPUT( : '-Opening input data image: '//DATA, : 1,0,ISTAT) ODD=.TRUE. ENDIF C Get the array size and shape CALL UIMGID(IDD,DATTYP,NDIMS,DDIMS,ISTAT) IF(NDIMS.NE.2) THEN CALL UMSPUT('! Input image is not two dimensional', : 1,0,ISTAT) GO TO 99 ELSE DNX=DDIMS(1) DNY=DDIMS(2) ENDIF C This is the EIS only version NOCO=.FALSE. EIS=.TRUE. HCOEFF=.FALSE. C Read some parameters which are EIS specific CALL UCLGSD('raref',RARED,ISTAT) CALL UCLGSD('decref',DECRED,ISTAT) C Convert these angles to radians RAREF=RARED*PIBY DECREF=DECRED*PIBY CALL UCLGSD('xpoff',XPOFF,ISTAT) CALL UCLGSD('ypoff',YPOFF,ISTAT) CALL UCLGSD('outscl',OUTSC,ISTAT) C Convert from arcsecs/pix to radians/pixel RPP=OUTSC/(3600.0/PIBY) C Get all the numbers needed from the header of the input C EIS image CALL EIGTCO(IDD,CRVAL1,CRPIX1,CRVAL2,CRPIX2, : CDELT1,CDELT2,MPIX,FLXSCL, : EXCO,EYCO,RAOFF,DECOFF,NC,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('-Unable to get EIS header information', : 1,0,ISTAT) GO TO 99 ENDIF C Get the input projection CTYPE1=' ' CALL UHDGST(IDD,'CTYPE1',CTYPE1,ISTAT) IF(ISTAT.NE.0) CALL UMSPUT( : '! Warning, unable to read CTYPE1 from header',1,0,ISTAT) CTYPE2=' ' CALL UHDGST(IDD,'CTYPE2',CTYPE2,ISTAT) IF(ISTAT.NE.0) CALL UMSPUT( : '! Warning, unable to read CTYPE2 from header',1,0,ISTAT) C We can now decide whether the input image is actually part of C the grid of the output image - ie, whether it has the same scale, C projection and reference point. If so we set the MPIX=-2 code C which is then passed through to the appropriate, simple, code. IF(ABS(CRVAL1-RAREF).LT.1E-7 .AND. : ABS(CRVAL2-DECREF).LT.1E-7 .AND. : NC.EQ.2 .AND. : ABS(EXCO(2)/RPP+1.0).LT.1E-7 .AND. : ABS(EYCO(1)/RPP-1.0).LT.1E-7 .AND. : ABS(EXCO(1)).LT.1E-7 .AND. : ABS(EYCO(2)).LT.1E-7 .AND. : CTYPE1.EQ.'RA---'//PROJ .AND. : CTYPE2.EQ.'DEC--'//PROJ) THEN MPIX=-2 CALL UMSPUT('-Input image grid is part of output grid', : 1,0,ISTAT) ENDIF C Now we calculate the corners of the input image on the C output and hence decide which subset of the output is C affected CALL CORNER(DNX,DNY, : CRVAL1,CRPIX1,CRVAL2,CRPIX2, : CDELT1,CDELT2,MPIX, : EXCO,EYCO,RAOFF,DECOFF,NC, : RAREF,DECREF,XPOFF,YPOFF,RPP,PROJ, : XMIN,XMAX,YMIN,YMAX,ISTAT) C Get the input exposure time from the header, if possible IF(EXPKEY.NE.' ') THEN CALL UHDGSR(IDD,EXPKEY,EXPIN,ISTAT) IF(ISTAT.NE.0) THEN IF(EXPSCL.OR.CPS) THEN CALL UMSPUT( : '! Warning, failed to get exposure time keyword value', : 1,0,ISTAT) GO TO 99 ENDIF C If we can't get it, but don't need it... EXPUPD=.FALSE. ENDIF ENDIF C Try to also get a "weighting" image and a scale factor for it CALL UCLGST('in_mask',WEIGHT,ISTAT) CALL UCLGST('wt_scl',WTSTR,ISTAT) IF(WEIGHT.EQ.' ') THEN USEWEI=.FALSE. ODW=.FALSE. ELSE CALL UIMOPN(WEIGHT,1,IDW,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Unable to open input weighting image', : 1,0,ISTAT) CALL UMSPUT('! Giving all data points equal weight', : 1,0,ISTAT) ODW=.FALSE. USEWEI=.FALSE. ELSE IF(VERBOSE) CALL UMSPUT( : '-Opening input weighting image: '//WEIGHT, : 1,0,ISTAT) ODW=.TRUE. C Get the array size and shape of the mask image - must C be the same as the data CALL UIMGID(IDW,DATTYP,NDIMS,DDIMS,ISTAT) IF(NDIMS.NE.2 .OR. DDIMS(1).NE.DNX .OR. : DDIMS(2).NE.DNY) THEN CALL UMSPUT('! Weighting image is not the same shape '// : 'and size as data',1,0,ISTAT) GO TO 99 ENDIF USEWEI=.TRUE. ENDIF ENDIF C The scale is read as a string as it might be "exptime" C in which case we use the exposure time - if we have it IF(WTSTR.EQ.'exptime' .OR. WTSTR.EQ.'EXPTIME') THEN IF(EXPUPD) THEN WTSCL=EXPIN ELSE CALL UMSPUT('! No exposure time value is available', : 1,0,ISTAT) GO TO 99 ENDIF ELSE IF(WTSTR.EQ.'header' .OR. WTSTR.EQ.'HEADER') THEN C Note in V0.24 this was changed to use the weight image rather C than the data image - but it was changed back again and C currently uses the data image again CALL UHDGSR(IDD,'WGTSCALE',WTSCL,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Warning, failed to get weighting header keyword value', : 1,0,ISTAT) CALL UMSPUT( : '! using default value of 1.0', : 1,0,ISTAT) WTSCL = 1.0 ENDIF ELSE READ(WTSTR,*,IOSTAT=ISTAT) WTSCL IF(ISTAT.NE.0) THEN CALL UMSPUT('! Weighting scale is not a number', : 1,0,ISTAT) GO TO 99 ENDIF ENDIF C We get the name of the output images here CALL UCLGST('outdata',OUTDAT,ISTAT) CALL UCLGST('outweig',OUTCOU,ISTAT) CALL UCLGST('outcont',CONTIM,ISTAT) C If the output weighting image is null we just don't C create anything and don't write anything out IF(OUTCOU.EQ.' ') THEN WRIWEI=.FALSE. ELSE WRIWEI=.TRUE. ENDIF C Similarly for the context image IF(CONTIM.EQ.' ') THEN CON=.FALSE. ELSE C Prepare a context table file name, just in case we need it later CALL LENSTR(CONTIM,I1,I2) CONTAB=CONTIM(I1:I2)//'.CON' CON=.TRUE. ENDIF C It may be that the "new" image currently exists and we C want to update it by drizzling the new data on top IF(OUTDAT.EQ.' ') THEN CALL UMSPUT( : '! No output data array name specified', : 1,0,ISTAT) GO TO 99 ENDIF C We try to open it with update access CALL UIMOPN(OUTDAT,2,IDND,ISTAT) IF(ISTAT.NE.0) THEN ODND=.FALSE. UPDATE=.FALSE. ELSE CALL UMSPUT( : '-Opening current data image: '//OUTDAT, : 1,0,ISTAT) ODND=.TRUE. C Now we try to open the corresponding weighting image CALL UIMOPN(OUTCOU,2,IDNC,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to open current weighting image', : 1,0,ISTAT) ODNC=.FALSE. GO TO 99 ENDIF CALL UMSPUT( : '-Opening current weighting image: '//OUTCOU, : 1,0,ISTAT) ODNC=.TRUE. C and also the context image IF(CON) THEN CALL UIMOPN(CONTIM,2,IDCO,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to open current context image', : 1,0,ISTAT) ODCO=.FALSE. GO TO 99 ENDIF CALL UMSPUT( : '-Opening current context image: '//CONTIM, : 1,0,ISTAT) ODCO=.TRUE. ENDIF C We are doing well and probably can use the existing C images however we need to do a little more... C Get the array size and shape of old data frame CALL UIMGID(IDND,DATTYP,NDIMS,ODIMS,ISTAT) IF(NDIMS.NE.2) THEN CALL UMSPUT('! Current image is not two dimensional', : 1,0,ISTAT) GO TO 99 ELSE ONX=ODIMS(1) ONY=ODIMS(2) ENDIF C Also check that the current weighting image is the same size CALL UIMGID(IDNC,DATTYP,NDIMS,ODIMS,ISTAT) IF(NDIMS.NE.2) THEN CALL UMSPUT( : '! Current weighting image is not two dimensional', : 1,0,ISTAT) GO TO 99 ELSE IF(ODIMS(1).NE.ONX .OR. ODIMS(2).NE.ONY) THEN CALL UMSPUT('! Count image not the same size as data', : 1,0,ISTAT) GO TO 99 ENDIF C and the context image IF(CON) THEN CALL UIMGID(IDCO,DATTYP,NDIMS,ODIMS,ISTAT) IF(NDIMS.NE.2) THEN CALL UMSPUT( : '! Current context image is not two dimensional', : 1,0,ISTAT) GO TO 99 ELSE IF(ODIMS(1).NE.ONX .OR. ODIMS(2).NE.ONY) THEN CALL UMSPUT('! Context image not the same size as data', : 1,0,ISTAT) GO TO 99 ENDIF ENDIF C We can now set the update flag UPDATE=.TRUE. WRIWEI=.TRUE. C Get the current exposure time from the header, if possible IF(EXPKEY.NE.' ') THEN CALL UHDGSR(IDND,EXPKEY,EXPCUR,ISTAT) IF(ISTAT.NE.0) THEN IF(EXPSCL.OR.CPS) THEN CALL UMSPUT( : '! Warning, failed to get exposure time keyword value', : 1,0,ISTAT) GO TO 99 ENDIF C Else we can't update the exposure header... EXPUPD=.FALSE. ENDIF ENDIF C Before doing anything we find out if this image has already C been drizzled onto this output. If so we abort. IF(EIDONE(IDND,DATA)) THEN CALL UMSPUT( : '! This image has already been drizzled onto this output', : 1,0,ISTAT) GO TO 99 ENDIF C We now get some of the previous drizzle information from C the header of the current image and check for funny C combinations - in V0.45 only IC is returned CALL EGOLDH(IDND,IC,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Warning, unable to read old drizzle header keywords', : 1,0,ISTAT) ELSE WRITE(CHARS, : '(''-This is image # '',I3,'' to be drizzled onto this output'' : )') IC+1 CALL UMSPUT(CHARS,1,0,ISTAT) ENDIF ENDIF IF(.NOT.UPDATE) IC=0 C If we are using the context we need to get the unique image C identifier from the data file IF(CON) THEN CALL GTUNID(DATA,IDD,UNIQID,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Unable to get unique image identifier', : 1,0,ISTAT) CALL UMSPUT('! using running count',1,0,ISTAT) UNIQID=IC+1 ENDIF WRITE(CHARS,'(''-Unique data image id: '',I6)') UNIQID CALL UMSPUT(CHARS,1,0,ISTAT) ENDIF C If there is no extant output image C we will have to create output data and weighting images anyway C so we need to find out their dimensions IF(.NOT.UPDATE) THEN CALL UCLGSI('outnx',ONX,ISTAT) CALL UCLGSI('outny',ONY,ISTAT) ODIMS(1)=ONX ODIMS(2)=ONY C If we are not updating we need to create the output images CALL UIMCRE(OUTDAT,6,2,ODIMS,IDND,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to create output data image', : 1,0,ISTAT) ODND=.FALSE. GO TO 99 ELSE CALL UMSPUT( : '-Created new output data image: '//OUTDAT, : 1,0,ISTAT) ODND=.TRUE. ENDIF C If we are going to store the counts image we must create it IF(WRIWEI) THEN CALL UIMCRE(OUTCOU,6,2,ODIMS,IDNC,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to create output weighting image', : 1,0,ISTAT) ODNC=.FALSE. GO TO 99 ELSE CALL UMSPUT( : '-Created new output weighting image: '//OUTCOU, : 1,0,ISTAT) ODNC=.TRUE. ENDIF ELSE CALL UMSPUT( : '! Warning, output weighting image will not be saved', : 1,0,ISTAT) ENDIF C If we are going to store the context image we must create it IF(CON) THEN CALL UIMCRE(CONTIM,3,2,ODIMS,IDCO,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to create output context image', : 1,0,ISTAT) ODCO=.FALSE. GO TO 99 ELSE CALL UMSPUT( : '-Created new output context image: '//CONTIM, : 1,0,ISTAT) ODCO=.TRUE. ENDIF ELSE CALL UMSPUT( : '! Warning, output context image will not be saved', : 1,0,ISTAT) ENDIF ENDIF C Update the exposure times in the outputs IF(EXPKEY.NE.' ') THEN IF(CPS) THEN EXPOUT=1.0 ELSE IF(UPDATE) THEN EXPOUT=EXPIN+EXPCUR ELSE EXPOUT=EXPIN ENDIF ENDIF C Modified in V0.85 to ADD rather than PUT CALL UHDASR(IDND,EXPKEY,EXPOUT, : 'Total exposure time of all input images',0,ISTAT) IF(WRIWEI) CALL UHDASR(IDNC,EXPKEY,EXPOUT, : 'Total exposure time of all input images',0,ISTAT) IF(CON) CALL UHDASR(IDCO,EXPKEY,EXPOUT, : 'Total exposure time of all input images',0,ISTAT) ENDIF C Assume we are using the whole frame unless set otherwise SUB=.FALSE. NOOVER=.FALSE. C Check for no-overlap IF(XMAX.LT.1 .OR. YMAX.LT.1 .OR. : XMIN.GT.ONX .OR. YMIN.GT.ONY) THEN CALL UMSPUT( : '! Warning, there is no overlap with the output image', : 1,0,ISTAT) WRITE(CHARS, : '(''! Output pixel range would be: ['',I5,'':'',I5, : '','',I5,'':'',I5,'']'')') XMIN,XMAX,YMIN,YMAX CALL UMSPUT(CHARS,1,0,ISTAT) NOOVER=.TRUE. IF(UPDATE) THEN GO TO 99 ELSE XMIN=1 XMAX=1 YMIN=1 YMAX=1 ENDIF SUB=.TRUE. ELSE XMIN=MAX(1,XMIN) XMAX=MIN(ONX,XMAX) YMIN=MAX(1,YMIN) YMAX=MIN(ONY,YMAX) WRITE(CHARS, : '(''-Drizzling onto subset: ['',I5,'':'',I5, : '','',I5,'':'',I5,'']'')') XMIN,XMAX,YMIN,YMAX CALL UMSPUT(CHARS,1,0,ISTAT) SUB=.TRUE. ENDIF C Now we allocate space for all the arrays we will need, note that they C are all real except for the context array which is integer*2 (at C present) C C First one line of the data array CALL UDMGET(DNX,6,PDATA,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to allocate memory for data array line', : 1,0,ISTAT) ADATA=.FALSE. GO TO 99 ELSE ADATA=.TRUE. ENDIF C Mask array - just one line at a time CALL UDMGET(DNX,6,PWEI,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to allocate memory for weights array', : 1,0,ISTAT) AWEI=.FALSE. GO TO 99 ELSE AWEI=.TRUE. ENDIF C Output/current data array CALL UDMGET((XMAX-XMIN+1)*(YMAX-YMIN+1),6,PNDAT,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to allocate memory for output data array', : 1,0,ISTAT) ANDAT=.FALSE. GO TO 99 ELSE ANDAT=.TRUE. ENDIF C Context array ANCON=.FALSE. IF(CON) THEN CALL UDMGET((XMAX-XMIN+1)*(YMAX-YMIN+1),3,PNCON,ISTAT) CALL UDMGET((XMAX-XMIN+1)*(YMAX-YMIN+1),3,PCDON,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to allocate memory for context arrays', : 1,0,ISTAT) GO TO 99 ELSE ANCON=.TRUE. ENDIF ENDIF C Weighting output/current array CALL UDMGET((XMAX-XMIN+1)*(YMAX-YMIN+1),6,PNCOU,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to allocate memory for output weighting array', : 1,0,ISTAT) ANCOU=.FALSE. GO TO 99 ELSE ANCOU=.TRUE. ENDIF C Some more memory arrays for EIS CALL UDMGET((DNX+1)*4,7,PXIN,ISTAT) CALL UDMGET((DNX+1)*4,7,PYIN,ISTAT) CALL UDMGET((DNX+1)*4,7,PXR1,ISTAT) CALL UDMGET((DNX+1)*4,7,PYR1,ISTAT) CALL UDMGET((DNX+1)*4,7,PXR2,ISTAT) CALL UDMGET((DNX+1)*4,7,PYR2,ISTAT) CALL UDMGET(DNX*4,7,PXTRA,ISTAT) CALL UDMGET(DNX*4,7,PYTRA,ISTAT) CALL UDMGET(ONX,6,PBUFF,ISTAT) CALL UDMGET(ONX,3,PSBUFF,ISTAT) C In the EIS case we get the WCS in a different way so we don't C need the CD matrix etc. GOTWCS=.FALSE. C Here we set the WCS matrix to something sensible anyway DO I=1,8 WCS(I)=0.0D0 ENDDO C and if we are updating the EIS header we specify C the appropriate projection - modified in V0.61 WRITE(CTYPE1,'(''RA---'',A3)') PROJ WRITE(CTYPE2,'(''DEC--'',A3)') PROJ C Optional weights array IF(.NOT.USEWEI) CALL SETIM(MEMR(PWEI),DNX,1,1.0) C If we are updating we read in the current "new" images IF(UPDATE) THEN IF(SUB) THEN DO I=YMIN,YMAX CALL UIGL2R(IDND,I,MEMR(PBUFF),ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Failed to read extant data image', : 1,0,ISTAT) GO TO 99 ENDIF CALL COPY1D(MEMR(PBUFF+XMIN-1), : MEMR(PNDAT+(I-YMIN)*(XMAX-XMIN+1)), : XMAX-XMIN+1) ENDDO ELSE DO I=1,ONY CALL UIGL2R(IDND,I,MEMR(PNDAT+(I-1)*ONX),ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Failed to read extant data image', : 1,0,ISTAT) GO TO 99 ENDIF ENDDO ENDIF C Scale by exposure time IF(EXPSCL) THEN CALL MULC(MEMR(PNDAT),XMAX-XMIN+1,YMAX-YMIN+1, : 1.0/EXPCUR) ENDIF C Weighting image IF(SUB) THEN DO I=YMIN,YMAX CALL UIGL2R(IDNC,I,MEMR(PBUFF),ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Failed to read current weighting image', : 1,0,ISTAT) GO TO 99 ENDIF CALL COPY1D(MEMR(PBUFF+XMIN-1), : MEMR(PNCOU+(I-YMIN)*(XMAX-XMIN+1)), : XMAX-XMIN+1) ENDDO ELSE DO I=1,ONY CALL UIGL2R(IDNC,I,MEMR(PNCOU+(I-1)*ONX),ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Failed to read current weighting image', : 1,0,ISTAT) GO TO 99 ENDIF ENDDO ENDIF C Context image (short integers) IF(CON) THEN IF(SUB) THEN DO I=YMIN,YMAX CALL UIGL2S(IDCO,I,MEMS(PSBUFF),ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Failed to read extant context image', : 1,0,ISTAT) GO TO 99 ENDIF CALL COPY1S(MEMS(PSBUFF+XMIN-1), : MEMS(PNCON+(I-YMIN)*(XMAX-XMIN+1)), : XMAX-XMIN+1) ENDDO ELSE DO I=1,ONY CALL UIGL2S(IDCO,I,MEMS(PNCON+(I-1)*ONX),ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Failed to read context data image', : 1,0,ISTAT) GO TO 99 ENDIF ENDDO ENDIF ENDIF ELSE C The final case is where the output is empty. C Before we start drizzling we need to drain out the output arrays C if there are no old arrays CALL SETIM(MEMR(PNDAT),XMAX-XMIN+1,YMAX-YMIN+1,0.0) CALL SETIM(MEMR(PNCOU),XMAX-XMIN+1,YMAX-YMIN+1,0.0) IF(CON) THEN CALL SETIMS(MEMS(PNCON),XMAX-XMIN+1,YMAX-YMIN+1,0) ENDIF ENDIF C If we are using a context array we need to either read its C current values from the global context file, or, if the C context is new, initialise it with something sensible C C In version 0.7 this is read from the header of the context C image rather than from an external file IF(CON) THEN CALL SETIMS(MEMS(PCDON),XMAX-XMIN+1,YMAX-YMIN+1,0) IF(UPDATE) THEN CALL GTCOIN(IDCO,INTAB,MAXIM,MAXEN,NEN,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Note: Valid context table not found in context image header', : 1,0,ISTAT) NEN=0 ENDIF ELSE NEN=0 ENDIF ENDIF C Now do the actual combination using BOXER - only if there is an C overlap IF(.NOT.NOOVER) THEN CALL DOBOX(MEMR(PDATA),MEMR(PWEI), : MEMR(PNDAT),MEMR(PNCOU),MEMS(PNCON),MEMS(PCDON), : DNX,DNY, : XMIN,XMAX,YMIN,YMAX, : ONX,ONY,WTSCL,ALIGN, : PFRACT,SCALE,ROT,XSH,YSH,WCS, : UPDATE,EXPSCL,EXPIN,IDD,USEWEI,IDW, : EIS,CRVAL1,CRPIX1,CRVAL2,CRPIX2, : CDELT1,CDELT2,MPIX,FLXSCL, : EXCO,EYCO,RAOFF,DECOFF,NC, : RAREF,DECREF,XPOFF,YPOFF,RPP,PROJ, : MEMD(PXIN),MEMD(PYIN),MEMD(PXTRA),MEMD(PYTRA), : MEMD(PXR1),MEMD(PYR1),MEMD(PXR2),MEMD(PYR2), : CON,INTAB,MAXIM,MAXEN,NEN,UNIQID, : ISTAT) ENDIF C Check the exit status IF(ISTAT.NE.0) GO TO 99 C Scale the output if not CPS output FS=1.0 IF(EXPUPD)THEN IF(UPDATE.AND..NOT.CPS) THEN IF(EXPSCL) THEN FS=EXPOUT ELSE FS=EXPOUT/EXPCUR ENDIF ENDIF C If we have a single image and are converting to CPS C we need to divide by the exposure length IF(.NOT.UPDATE.AND.CPS) THEN FS=1.0/EXPIN ENDIF ENDIF C Scale the data values by exposure time IF(FS.NE.1.0) CALL MULC(MEMR(PNDAT),XMAX-XMIN+1, : YMAX-YMIN+1,FS) C Now set the scale factor to include the pre-scaling as well IF(UPDATE.AND.EXPSCL) FS=FS/EXPCUR C Put in the fill values (if defined) IF(FILL) CALL PUTFIL(MEMR(PNDAT),MEMR(PNCOU), : XMAX-XMIN+1,YMAX-YMIN+1,FILVAL) C Write out the two new images IF(VERBOSE) : CALL UMSPUT('-Writing output drizzled image: '// : OUTDAT,1,0,ISTAT) C If we are not updating (ie, this is the first drizzle C onto this output) then we modify the WCS and copy other stuff IF(.NOT.UPDATE) THEN CALL EUWCS(IDND,WCS,CTYPE1,CTYPE2,DECRED) CALL EUCPHD(IDD,IDND) ENDIF C Update the header - we do this first to make writing C out more efficient (?) CALL EUHEAD(IDND,VERS,DATA,FLXSCL,WTSCL) 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. (EXPUPD .AND. .NOT.CPS)) THEN IF(YMIN.GT.1) THEN DO I=1,YMIN-1 IF(UPDATE) THEN CALL UIGL2R(IDND,I,MEMR(PBUFF),ISTAT) IF(FS.NE.1.0) THEN CALL MULC(MEMR(PBUFF),ONX,1,FS) ENDIF ELSE CALL SETIM(MEMR(PBUFF),ONX,1,0.0) ENDIF CALL UIPL2R(IDND,I,MEMR(PBUFF),ISTAT) ENDDO ENDIF IF(YMAX.LT.ONY) THEN DO I=YMAX+1,ONY IF(UPDATE) THEN CALL UIGL2R(IDND,I,MEMR(PBUFF),ISTAT) IF(FS.NE.1.0) THEN CALL MULC(MEMR(PBUFF),ONX,1,FS) ENDIF ELSE CALL SETIM(MEMR(PBUFF),ONX,1,0.0) ENDIF CALL UIPL2R(IDND,I,MEMR(PBUFF),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(IDND,I,MEMR(PBUFF),ISTAT) IF(FS.NE.1.0) THEN CALL MULC(MEMR(PBUFF),ONX,1,FS) ENDIF ELSE CALL SETIM(MEMR(PBUFF),ONX,1,0.0) ENDIF CALL COPY1D(MEMR(PNDAT+(XMAX-XMIN+1)*(I-YMIN)), : MEMR(PBUFF+XMIN-1),XMAX-XMIN+1) CALL UIPL2R(IDND,I,MEMR(PBUFF),ISTAT) ENDDO IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to write output data image', : 1,0,ISTAT) GO TO 99 ENDIF ELSE DO I=1,ONY CALL UIPL2R(IDND,I,MEMR(PNDAT+(I-1)*ONX),ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to write output data image', : 1,0,ISTAT) GO TO 99 ENDIF ENDDO ENDIF C Now write out the weights array, if there is one IF(WRIWEI) THEN IF(VERBOSE) : CALL UMSPUT('-Writing output weighting image: '// : OUTCOU,1,0,ISTAT) C Update the header IF(.NOT.UPDATE) THEN CALL EUWCS(IDNC,WCS,CTYPE1,CTYPE2,DECRED) CALL EUCPHD(IDD,IDNC) ENDIF CALL EUHEAD(IDNC,VERS,DATA,FLXSCL,WTSCL) IF(SUB) THEN C First we read the ranges of Y for which the drizzle C has had no effect. IF(.NOT.UPDATE .OR. (EXPUPD .AND. .NOT.CPS)) THEN IF(YMIN.GT.1) THEN DO I=1,YMIN-1 IF(UPDATE) THEN CALL UIGL2R(IDNC,I,MEMR(PBUFF),ISTAT) ELSE CALL SETIM(MEMR(PBUFF),ONX,1,0.0) ENDIF CALL UIPL2R(IDNC,I,MEMR(PBUFF),ISTAT) ENDDO ENDIF IF(YMAX.LT.ONY) THEN DO I=YMAX+1,ONY IF(UPDATE) THEN CALL UIGL2R(IDNC,I,MEMR(PBUFF),ISTAT) ELSE CALL SETIM(MEMR(PBUFF),ONX,1,0.0) ENDIF CALL UIPL2R(IDNC,I,MEMR(PBUFF),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(IDNC,I,MEMR(PBUFF),ISTAT) ELSE CALL SETIM(MEMR(PBUFF),ONX,1,0.0) ENDIF CALL COPY1D(MEMR(PNCOU+(XMAX-XMIN+1)*(I-YMIN)), : MEMR(PBUFF+XMIN-1),XMAX-XMIN+1) CALL UIPL2R(IDNC,I,MEMR(PBUFF),ISTAT) ENDDO IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to write output weight image', : 1,0,ISTAT) GO TO 99 ENDIF ELSE DO I=1,ONY CALL UIPL2R(IDNC,I,MEMR(PNCOU+(I-1)*ONX),ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to write output weight image', : 1,0,ISTAT) GO TO 99 ENDIF ENDDO ENDIF ENDIF C Finally write the context array, which is INTEGER*2, if we have C one IF(CON) THEN IF(VERBOSE) : CALL UMSPUT('-Writing output context image: '// : CONTIM,1,0,ISTAT) C This was added in V0.33 IF(.NOT.UPDATE) THEN CALL EUWCS(IDCO,WCS,CTYPE1,CTYPE2,DECRED) CALL EUCPHD(IDD,IDCO) ENDIF IF(SUB) THEN C First we read the ranges of Y for which the drizzle C has had no effect. IF(.NOT.UPDATE) THEN IF(YMIN.GT.1) THEN DO I=1,YMIN-1 IF(UPDATE) THEN CALL UIGL2S(IDCO,I,MEMS(PSBUFF),ISTAT) ELSE CALL SETIMS(MEMS(PSBUFF),ONX,1,0) ENDIF CALL UIPL2S(IDCO,I,MEMS(PSBUFF),ISTAT) ENDDO ENDIF IF(YMAX.LT.ONY) THEN DO I=YMAX+1,ONY IF(UPDATE) THEN CALL UIGL2S(IDCO,I,MEMS(PSBUFF),ISTAT) ELSE CALL SETIMS(MEMS(PSBUFF),ONX,1,0) ENDIF CALL UIPL2S(IDCO,I,MEMS(PSBUFF),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(IDCO,I,MEMS(PSBUFF),ISTAT) ELSE CALL SETIMS(MEMS(PSBUFF),ONX,1,0) ENDIF CALL COPY1S(MEMS(PNCON+(XMAX-XMIN+1)*(I-YMIN)), : MEMS(PSBUFF+XMIN-1),XMAX-XMIN+1) CALL UIPL2S(IDCO,I,MEMS(PSBUFF),ISTAT) ENDDO IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to write output context image', : 1,0,ISTAT) GO TO 99 ENDIF ELSE DO I=1,ONY CALL UIPL2S(IDCO,I,MEMS(PNCON+(I-1)*ONX),ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to write output context image', : 1,0,ISTAT) GO TO 99 ENDIF ENDDO ENDIF ENDIF IF(CON) THEN C Write the name of the data file C Changed from CD to DD in version 0.55 WRITE(KEY,'(''DD'',I6)') UNIQID DO J=3,8 IF(KEY(J:J).EQ.' ') KEY(J:J)='0' ENDDO WRITE(CHARS,'(''Data image with unique id # '',I7)') UNIQID CALL UHDAST(IDND,KEY,DATA,CHARS,0,ISTAT) CALL UHDAST(IDNC,KEY,DATA,CHARS,0,ISTAT) IF(CON) CALL UHDAST(IDCO,KEY,DATA,CHARS,0,ISTAT) C Now update the global context table in header CALL PTCOIN(IDCO,CONTAB,MAXHCN,INTAB,MAXIM,MAXEN,NEN,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to update context in header', : 1,0,ISTAT) GO TO 99 ELSE CALL UMSPUT( : '-Context table written to context image header', : 1,0,ISTAT) ENDIF ENDIF C General abort point for this module 99 CONTINUE C Free up all allocated memory arrays IF(ADATA) CALL UDMFRE(PDATA,6,ISTAT) IF(AWEI) CALL UDMFRE(PWEI,6,ISTAT) IF(ANCOU) CALL UDMFRE(PNCOU,6,ISTAT) IF(ANCON) THEN CALL UDMFRE(PNCON,3,ISTAT) CALL UDMFRE(PCDON,3,ISTAT) ENDIF IF(ANDAT) CALL UDMFRE(PNDAT,6,ISTAT) C Close all the image files which may have been open IF(ODD) CALL UIMCLO(IDD,ISTAT) IF(ODW) CALL UIMCLO(IDW,ISTAT) IF(ODND) CALL UIMCLO(IDND,ISTAT) IF(ODNC) CALL UIMCLO(IDNC,ISTAT) IF(ODCO) CALL UIMCLO(IDCO,ISTAT) C End of main DRIZZLE module WRITE(CHARS, : '(''+ Total elapsed drizzle time: '',I5,'' (s)'')') : IITIME()-STIME CALL UMSPUT(CHARS,1,0,ISTAT) CALL UMSPUT('+ Drizzle ending at '//ICTIME(IITIME()),1,0,ISTAT) END SUBROUTINE EUHEAD(ID,VERS,DATA,FLXSCL,WTSCL) C C Update the header of the output image with all the C parameters of the current DRIZZLE run. EIS version. C C Much shortened in V0.45 C IMPLICIT NONE INTEGER ID,ISTAT,I,IITIME REAL WTSCL DOUBLE PRECISION FLXSCL CHARACTER*80 DATA,BUFFER CHARACTER*45 VERS CHARACTER*50 COMM CHARACTER*4 STEM CHARACTER*24 ICTIME 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 IF(STEM.EQ.'D001') THEN COMM='Drizzle, task version' ELSE COMM=' ' ENDIF CALL UHDAST(ID,STEM//'VER',VERS,COMM,0,ISTAT) C And the now the time of the run IF(STEM.EQ.'D001') THEN COMM='Drizzle, time of end of run' ELSE COMM=' ' ENDIF CALL UHDAST(ID,STEM//'TIME',ICTIME(IITIME()),COMM,0,ISTAT) C Now we continue to add the other items using the same C "stem" IF(STEM.EQ.'D001') THEN COMM='Drizzle, input data image' ELSE COMM=' ' ENDIF CALL UHDAST(ID,STEM//'DATA',DATA,COMM,0,ISTAT) IF(STEM.EQ.'D001') THEN COMM='Drizzle, flux scale applied to input' ELSE COMM=' ' ENDIF CALL UHDASD(ID,STEM//'FLSC',FLXSCL,COMM,0,ISTAT) IF(STEM.EQ.'D001') THEN COMM='Drizzle, weight scale applied to input' ELSE COMM=' ' ENDIF CALL UHDASR(ID,STEM//'WTSC',WTSCL,COMM,0,ISTAT) RETURN END SUBROUTINE DOBOX(DATA,WEI,NDAT,NCOU,NCON,DONE,DNX,DNY, : XMIN,XMAX,YMIN,YMAX, : ONX,ONY,WTSCL,ALIGN, : PFRACT,SCALE,ROT,XSH,YSH,WCS, : UPDATE,EXPSCL,EXPIN,IDD,USEWEI,IDW, : EIS,CRVAL1,CRPIX1,CRVAL2,CRPIX2, : CDELT1,CDELT2,MPIX,FLXSCL, : EXCO,EYCO,RAOFF,DECOFF,NC, : RAREF,DECREF,XPOFF,YPOFF,RPP,PROJ, : XIN,YIN,XTRA,YTRA, : XR1,YR1,XR2,YR2, : CON,INTAB,MAXIM,MAXEN,NEN,UNIQID, : 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 This version is for EIS C IMPLICIT NONE INTEGER DNX,DNY,ONX,ONY,IDD,IDW INTEGER I,J,II,JJ,NMISS,NSKIP,ISTAT,NHIT,IS,K INTEGER XMIN,XMAX,YMIN,YMAX REAL DATA(DNX),WEI(DNX),W,EXPIN REAL NDAT(XMAX-XMIN+1,YMAX-YMIN+1) REAL NCOU(XMAX-XMIN+1,YMAX-YMIN+1) INTEGER*2 NCON(XMAX-XMIN+1,YMAX-YMIN+1) INTEGER*2 DONE(XMAX-XMIN+1,YMAX-YMIN+1) REAL D,OVER DOUBLE PRECISION WCS(8),DXMIN,DYMIN REAL VC,XSH,YSH,DOVER,DH REAL XI,XA,YI,YA REAL XCEN,YCEN REAL ROT,XF,YF,SCALE REAL PFRACT,S2 REAL AC,WTSCL,JACO CHARACTER*80 CHARS CHARACTER*8 ALIGN CHARACTER*3 PROJ LOGICAL UPDATE,EXPSCL,USEWEI C EIS specific stuff INTEGER MPIX,NC,WC DOUBLE PRECISION CRVAL1,CRPIX1,CRVAL2,CRPIX2 DOUBLE PRECISION CDELT1,CDELT2,EXCO(100),EYCO(100) DOUBLE PRECISION RAREF,DECREF,XPOFF,YPOFF,RPP DOUBLE PRECISION RAOFF,DECOFF,FLXSCL REAL FLX,S2FLX,DOW LOGICAL EIS,PF1,CON,MATCH,DIDLAS INTEGER MAXEN,MAXIM,NEN,ICON,NN INTEGER UNIQID,INTAB(MAXIM,MAXEN) INTEGER*2 OLDCON,NEWCON INTEGER NEWMA(100) DOUBLE PRECISION DDH,DJ,XH,YH DOUBLE PRECISION XIN(DNX+1,4),YIN(DNX+1,4) DOUBLE PRECISION XTRA(DNX,4),YTRA(DNX,4) DOUBLE PRECISION XR1(DNX+1),YR1(DNX+1) DOUBLE PRECISION XR2(DNX+1),YR2(DNX+1) DOUBLE PRECISION XX(4),YY(4),XXOUT(4),YYOUT(4) INTEGER IXMIN,IXMAX,IYMIN,IYMAX INTEGER NXI,NXA,NYI,NYA C Constant DOUBLE PRECISION PIBY PARAMETER (PIBY=3.141592653589793/180.0) C-- C Some initial settings NMISS=0 NSKIP=0 FLX=SNGL(FLXSCL) OLDCON=-1 DIDLAS=.FALSE. C Calculate the position of the first whole pixel XF=1.0/SCALE YF=1.0/SCALE C Some more economies S2=1.0/(XF*YF) S2FLX=S2*FLX DH=0.5*PFRACT DDH=DBLE(DH) AC=1.0/(PFRACT*PFRACT) DXMIN=1.0-DBLE(XMIN) DYMIN=1.0-DBLE(YMIN) C We make a special case if PFRACT=1.0 so we check for this and C set a flag IF(ABS(PFRACT-1.0).LT.1.0E-4) THEN PF1=.TRUE. ELSE PF1=.FALSE. ENDIF C First we set the input X values, which are the same for all C rows in the image DO I=1,DNX+1 XIN(I,1)=DBLE(I) ENDDO C Before we start we work out the determinant of the Jacobian at the centre XX(1)=DBLE(DNX/2) XX(2)=DBLE(DNX/2+1) XX(3)=DBLE(DNX/2+1) XX(4)=DBLE(DNX/2) YY(1)=DBLE(DNY/2) YY(2)=DBLE(DNY/2) YY(3)=DBLE(DNY/2+1) YY(4)=DBLE(DNY/2+1) CALL EIEVAL(XX,YY,4, : CRVAL1,CRPIX1,CRVAL2,CRPIX2, : CDELT1,CDELT2,MPIX, : EXCO,EYCO,RAOFF,DECOFF,NC, : RAREF,DECREF,XPOFF,YPOFF,RPP,PROJ, : XXOUT,YYOUT,ISTAT) C Calculate the area of the quadrilateral and make sure C it is positive. Then calculate the size of a square of C the same area and the half-lengths of its sides JACO=0.5*((XXOUT(2)-XXOUT(4))*(YYOUT(1)-YYOUT(3)) - : (XXOUT(1)-XXOUT(3))*(YYOUT(2)-YYOUT(4))) IF(JACO.LT.0.0) JACO=JACO*-1.0 XH=0.5*SQRT(JACO*AC)*PFRACT YH=0.5*SQRT(JACO*AC)*PFRACT C Normalise flux to include the larger area if the Jacobian is C greater than 1.0 S2FLX=S2FLX/JACO C Outer loop over input image pixels (X,Y) DO J=1,DNY C Firstly check that this line overlaps the output image C at all. Only continue if it does XX(1)=1.0D0 XX(2)=DBLE(DNX) YY(1)=DBLE(J) YY(2)=YY(1) CALL EIEVAL(XX,YY,2, : CRVAL1,CRPIX1,CRVAL2,CRPIX2, : CDELT1,CDELT2,MPIX, : EXCO,EYCO,RAOFF,DECOFF,NC, : RAREF,DECREF,XPOFF,YPOFF,RPP,PROJ, : XXOUT,YYOUT,ISTAT) IXMIN=NINT(MIN(XXOUT(1),XXOUT(2))) IYMIN=NINT(MIN(YYOUT(1),YYOUT(2))) IXMAX=NINT(MAX(XXOUT(1),XXOUT(2))) IYMAX=NINT(MAX(YYOUT(1),YYOUT(2))) C Check for overlap - with a margin of 10 to be cautious IF(IXMAX.GT.-10 .AND. IXMIN.LT.ONX+10 .AND. : IYMAX.GT.-10 .AND. IYMIN.LT.ONY+10) THEN DJ=DBLE(J) DO I=1,DNX YIN(I,1)=DJ ENDDO CALL EIEVAL(XIN(1,1),YIN(1,1),DNX, : CRVAL1,CRPIX1,CRVAL2,CRPIX2, : CDELT1,CDELT2,MPIX, : EXCO,EYCO,RAOFF,DECOFF,NC, : RAREF,DECREF,XPOFF,YPOFF,RPP,PROJ, : XR1,YR1,ISTAT) 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 we are scaling the exposures we need to divide by C the exposure IF(UPDATE.AND.EXPSCL) 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 input weighting image', : 1,0,IS) RETURN ENDIF ENDIF DO I=1,DNX XCEN=XR1(I)+DXMIN YCEN=YR1(I)+DYMIN IF(MPIX.EQ.-2) THEN XI=XCEN XA=XI YI=YCEN YA=YI ELSE XI=XCEN-XH XA=XCEN+XH YI=YCEN-YH YA=YCEN+YH ENDIF NXI=NINT(XI) NXA=NINT(XA) NYI=NINT(YI) NYA=NINT(YA) C Initialise the pixels-affected counter NHIT=0 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.XMAX-XMIN+1 .AND. : JJ.GE.1 .AND. JJ.LE.YMAX-YMIN+1) THEN C Calculate the overlap using the simpler "aligned" box routine C or, in the case of a matched grid (MPIX=-2) then just set to C 1.0 IF(MPIX.EQ.-2) THEN DOVER=1.0 ELSE DOVER=OVER(II,JJ,XI,XA,YI,YA) ENDIF C Re-normalise the area overlap using the Jacobian C Note that in V0.27 and later this statement is done C before the overlap is checked for so that flipped images C are handled correctly C C The division was moved back after the test in V0.52 C when JACO is definitely positive IF(DOVER.GT.0.0) THEN DOVER=DOVER/JACO C Count the hits NHIT=NHIT+1 VC=NCOU(II,JJ) C Allow for stretching because of scale change D=DATA(I)*S2FLX 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 DOW=DOVER*W C Update the context image, if requested IF(CON .AND. W.GT.0.0 .AND. DONE(II,JJ).EQ.0) THEN 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 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 If using pixfrac=1 we now swap to the other way around IF(WC.EQ.1) THEN WC=2 ELSE WC=1 ENDIF ELSE C Count skipped lines NSKIP=NSKIP+1 IF(DIDLAS) DIDLAS=.FALSE. ENDIF ENDDO C Report number of points which were outside the output frame IF(NMISS.GT.0) THEN WRITE(CHARS,'(''! Note, '',I8,'' points were '','// : '''outside the output image.'')') : NMISS CALL UMSPUT(CHARS,1,0,ISTAT) ENDIF C and report the number of lines skipped IF(NSKIP.GT.0) THEN WRITE(CHARS,'(''! Note, '',I8,'// : ''' lines of the input image were skipped.'')') NSKIP CALL UMSPUT(CHARS,1,0,ISTAT) ENDIF C We now try to update the WCS C First CRPIX1,CRPIX2,CRVAL1,CRVAL2 WCS(1)=-XPOFF WCS(2)=RAREF/PIBY WCS(3)=-YPOFF WCS(4)=DECREF/PIBY WCS(5)=-RPP/PIBY WCS(6)=0.0 WCS(7)=0.0 WCS(8)=RPP/PIBY C Set good status if we get this far ISTAT=0 RETURN END 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 SUBROUTINE BOXER (IS, JS, X, Y, DAREA) 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 ISTAT,IS,JS REAL M, C, DX REAL XLO, XHI, YLO, YHI, XTOP CHARACTER*80 CHARS 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 C Check for a problem - note that this now has a small C error margin included to avoid errors coming from numerical C inaccuracy IF (XTOP .LT. XLO-1.0E-5 .OR. XTOP .GT. XHI+1.0E-5) THEN WRITE(CHARS,'(''! Warning, box overlap calculation'','// : ''' problem at ('',I5,'','',I5,'')'')') IS,JS CALL UMSPUT(CHARS,1,0,ISTAT) ENDIF 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 LOGICAL FUNCTION EIDONE(ID,DATA) C C This routine checks whether the input data array has already C been drizzled onto a given output image. If so it returns TRUE, C otherwise FALSE. C INTEGER ID CHARACTER*80 DATA CHARACTER*4 STEM CHARACTER*80 BUFFER INTEGER IC,ISTAT 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//'DATA',BUFFER,ISTAT) C Check for a status code indicating that this header C item doesn't exist IF(ISTAT.NE.0) THEN EIDONE=.FALSE. RETURN ENDIF C Check for a match IF(BUFFER.EQ.DATA) THEN EIDONE=.TRUE. RETURN ENDIF ENDDO EIDONE=.FALSE. RETURN END SUBROUTINE EGOLDH(ID,IC,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 is the EIS version, May 1997 C IMPLICIT NONE INTEGER IC,ISTAT,ID CHARACTER*4 STEM CHARACTER*80 BUFFER 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 C In Version 0.45 we just return here ISTAT=0 RETURN END SUBROUTINE EUWCS(ID,WCS,CTYPE1,CTYPE2,DECRED) C C Update the WCS header items using the 8 values C supplied in WCS and the values of CTYPE1 and CTYPE2. C C This is the EIS version which also updates C the PROJ keywords required by C the FITS definition of the COE projection. C C The image is assumed to be available with descriptor ID C and to already have all the parameters C C Modified so that the CD matrix is added - not modified - C Richard Hook, ST-ECF, July 1997 C C Modified so that full WCS is supplied and the items are C deleted before being added, Dec 97 C C After version 0.8 of EIS Drizzle CDELTs are no longer written C but EQUINOX is as well as RADESYS which are always set to C standard values (FK5 and 2000.0). C Richard Hook, ST-ECF, 25th January 2000 C IMPLICIT NONE INTEGER ID,ISTAT CHARACTER*8 CTYPE1,CTYPE2 DOUBLE PRECISION WCS(8),DECRED C Write a full WCS for the EIS COE projection CALL UHDDSP(ID,'CTYPE1',ISTAT) CALL UHDAST(ID,'CTYPE1',CTYPE1, : 'WCS projection type',0,ISTAT) CALL UHDDSP(ID,'CTYPE2',ISTAT) CALL UHDAST(ID,'CTYPE2',CTYPE2, : 'WCS projection type',0,ISTAT) CALL UHDDSP(ID,'CRPIX1',ISTAT) CALL UHDASD(ID,'CRPIX1',WCS(1), : 'Reference pixel (X)',0,ISTAT) CALL UHDDSP(ID,'CRVAL1',ISTAT) CALL UHDASD(ID,'CRVAL1',WCS(2), : 'Reference point on sky (RA,degrees)',0,ISTAT) CALL UHDDSP(ID,'CRPIX2',ISTAT) CALL UHDASD(ID,'CRPIX2',WCS(3), : 'Reference pixel (Y)',0,ISTAT) CALL UHDDSP(ID,'CRVAL2',ISTAT) CALL UHDASD(ID,'CRVAL2',WCS(4), : 'Reference point on sky (dec,degrees)',0,ISTAT) CALL UHDDSP(ID,'CD1_1',ISTAT) CALL UHDASD(ID,'CD1_1',WCS(5), : 'CD1_1 matrix element',0,ISTAT) CALL UHDDSP(ID,'CD2_1',ISTAT) CALL UHDASD(ID,'CD2_1',WCS(6), : 'CD2_1 matrix element',0,ISTAT) CALL UHDDSP(ID,'CD1_2',ISTAT) CALL UHDASD(ID,'CD1_2',WCS(7), : 'CD1_2 matrix element',0,ISTAT) CALL UHDDSP(ID,'CD2_2',ISTAT) CALL UHDASD(ID,'CD2_2',WCS(8), : 'CD2_2 matrix element',0,ISTAT) IF(CTYPE1.EQ.'RA---COE') THEN CALL UHDDSP(ID,'PROJP1',ISTAT) CALL UHDASD(ID,'PROJP1',DECRED, : 'COE projection reference declination (degrees)',0,ISTAT) CALL UHDDSP(ID,'PROJP2',ISTAT) CALL UHDASD(ID,'PROJP2',0.0D0, : 'COE projection - separation of reference decs', : 0,ISTAT) ENDIF CALL UHDDSP(ID,'EQUINOX',ISTAT) CALL UHDASR(ID,'EQUINOX',2000.0, : 'Equinox of positions (assumed)',0,ISTAT) CALL UHDDSP(ID,'RADESYS',ISTAT) CALL UHDAST(ID,'RADESYS','FK5', : 'Coordinate system (assumed)',0,ISTAT) CALL UHDDSP(ID,'FLXSCALE',ISTAT) CALL UHDASD(ID,'FLXSCALE',1.0D0,'Flux scale factor', : 0,ISTAT) CALL UHDDSP(ID,'WGTSCALE',ISTAT) CALL UHDASD(ID,'WGTSCALE',1.0D0,'Weight scale factor', : 0,ISTAT) RETURN END SUBROUTINE EUCPHD(IN,OUT) C C Copy some useful values from the input image to the output. C These are things which are not usually changed by the processing - C such as the instrument name, filter etc. C C If they are not present they are quietly ignored. C C Richard Hook, EIS/ST-ECF, August 2000 C IMPLICIT NONE INTEGER IN,OUT INTEGER KWL,ISTAT CHARACTER*8 VALUE,KEY C First the instrument name CALL UHDGST(IN,'INSTRUME',VALUE,ISTAT) IF(ISTAT.EQ.0) CALL UHDAST(OUT,'INSTRUME', : VALUE,'Instrument name from input frame',0,ISTAT) C Now copy all the '*FILT*' entries C Initialise header keyword matching... CALL UHDOKL(IN,'*FILT*',.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 Get the value CALL UHDGST(IN,KEY,VALUE,ISTAT) C Copy it over CALL UHDAST(OUT,KEY,VALUE, : 'Filter information from input header', : 0,ISTAT) ENDDO 77 CONTINUE C Close the keyword loop CALL UHDCKL(KWL,ISTAT) END