C @(#)old_mcoadd.for	17.1.1.1 (ES0-DMD) 01/25/02 17:15:11
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
CI*I  subroutine acoadd   !I
      PROGRAM MCOADD      !M
C++
C
C MCOADD.F Version 1.1M (MIDAS version)
C
C Lucy co-addition, FFT version with acceleration option.
C
C This code is a preliminary and UNSUPPORTED implementation of
C an algorithm for combining images which have different PSFs.
C
C More details of the method are given in the ST-ECF Newsletter,
C 17, Feb 92, p10 and references given therein.
C
C Apart from the standard libraries the only other external
C dependencies are:
C
C fftetc.f - a set of FFT support routines including the Numerical
C            Recipes FFT code.
C
C midint.f - a small library of routines to emulate the F77/VOS
C            interface from IRAF in the MIDAS environment.
C
C Test version, Richard Hook, ST-ECF, Dec 1991
C               e-mail: rhook@eso.org (Internet)
C                       eso::rhook    (SPAN)
C                       rhook@dgaeso51.bitnet
C
C Acceleration added, Richard Hook, ST-ECF, Mar 1992
C Double precision (partly) version, Mar 1992
C Starting image option added, Mar 1992
C Fully double precision version, May 1992
C Weighting error corrected, May 1992
C Additional checks added, May 1992
C Added residual information, June 1992
C Version 1.0, 2nd June 1992
C One D support, 4th June 1992
C Sub-sampling support added, 29th July 1992
C MIDAS compatible version, 7th August 1992
C MIDAS version for release, 15th September 1992
C--
C
      IMPLICIT NONE
                   
      INTEGER MAXIMS,USEOF
      PARAMETER (MAXIMS=10)
      PARAMETER (USEOF=-2)
                          
C Iraf global storage
CI    double precision Memd(1) !I
CI    common /Mem/Memd         !I
                                 
C MIDAS global storage
      INTEGER MEMD(1)          !M
      COMMON /VMR/MEMD         !M
                                 
C Local variables
      INTEGER DATTYP(MAXIMS),NDIMS
      INTEGER DATPNT(MAXIMS),PSFPNT(MAXIMS)
      INTEGER ISTAT,NX,NY,N,M,NITER,IMSTAT,DIMS(2)
      INTEGER PSFFFT(MAXIMS),WORK,REST,W,DPS
      INTEGER CORRFAC,PH(MAXIMS)
      INTEGER IDD(MAXIMS),IDP(MAXIMS),IDDC,IDCOAD(MAXIMS),FID
      DOUBLE PRECISION WEIGHT(MAXIMS),T,FMAX,PT
      DOUBLE PRECISION RMSRES,RESMAX
      INTEGER IRMAX,JRMAX
      CHARACTER*80 DATA(MAXIMS),PSF(MAXIMS),FIM
      CHARACTER*80 COADDR,COADD(MAXIMS),DECON
      CHARACTER*80 CHARS
      CHARACTER*80 IMAGES,PSFS
      INTEGER NIM,NPSF
      INTEGER IMLD,PSFLD,IST,IEND
      INTEGER INR
      INTEGER NBX,NBY
      INTEGER I
      DOUBLE PRECISION DRV1,DRV2,DD1,DD2,HH,HHT
      DOUBLE PRECISION XLP,XL1,XL
      LOGICAL ACC
      LOGICAL VERBOSE
      LOGICAL FIRST
      LOGICAL POWER2
      LOGICAL SUBSAM
C++
   
C Start of code
               
      CALL STSPRO('MCOADD') !M
                              
C First get the boolean control logicals
      CALL UCLGSB('Verbose',VERBOSE,ISTAT)
      CALL UCLGSB('Accel',ACC,ISTAT)
                                    
C Announce the version
C
      IF(VERBOSE) CALL UMSPUT('+ MCOADD Version 1.1M (Sept. 1992)',
     :            1,0,ISTAT)
                            
C Initialise the template filename processing
      CALL UCLGST('Images',IMAGES,ISTAT)
      CALL TIMOTP(IMAGES,IMLD,ISTAT)
                                    
C Get the names of the input data images
      NIM=1
      IMSTAT=0
      DO I=1,MAXIMS+1
                               
C Check that there are not too many images
         IF(NIM.GT.MAXIMS) THEN
            CALL UMSPUT('Too many images specified',
     :                  1,0,ISTAT)
            GO TO 99
         ENDIF
              
         CALL TIMXTP(IMLD,DATA(NIM),IMSTAT)
         IF(IMSTAT.EQ.USEOF) GO TO 77
         CALL UIMOPN(DATA(NIM),1,IDD(NIM),ISTAT)
         IF(ISTAT.NE.0) THEN
               CALL UMSPUT('Unable to open data image',
     :                    1,0,ISTAT)
            GO TO 99
         ELSE
             IF(VERBOSE) CALL UMSPUT
     :         ('-Input Image: '//DATA(NIM),1,0,ISTAT)
         ENDIF
         NIM=NIM+1
      ENDDO
           
 77   CONTINUE
      NIM=NIM-1
               
      CALL TIMCTP(IMLD,ISTAT)
                             
C Open the list of PSF names
      CALL UCLGST('PSFs',PSFS,ISTAT)
      CALL TIMOTP(PSFS,PSFLD,ISTAT)
                                   
      NPSF=1
      IMSTAT=0
      DO I=1,MAXIMS+1
         CALL TIMXTP(PSFLD,PSF(NPSF),IMSTAT)
         IF(IMSTAT.EQ.USEOF) GO TO 78
         CALL UIMOPN(PSF(NPSF),1,IDP(NPSF),ISTAT)
         IF(ISTAT.NE.0) THEN
            CALL UMSPUT('Unable to open PSF image',
     :                    1,0,ISTAT)
            GO TO 99
         ELSE
             IF(VERBOSE) CALL UMSPUT('-Input PSF: '//PSF(NPSF),
     :                   1,0,ISTAT)
         ENDIF
         NPSF=NPSF+1
      ENDDO
           
 78   CONTINUE
      NPSF=NPSF-1
                 
      CALL TIMCTP(PSFLD,ISTAT)
                              
C Do some checks on PSFs size, number etc
      IF(NPSF.NE.NIM) THEN
         CALL UMSPUT('Different number of PSFs and Images',
     :                 1,0,ISTAT)
         GO TO 99
      ELSE
         WRITE(CHARS,'(''-There are '',i4,'' images/PSFs'')') NIM
         IF(VERBOSE) CALL UMSPUT(CHARS,1,0,ISTAT)
                                                 
C In either case we set the acceleration to 1.0 as default
         XL=1.0D0
      ENDIF
           
C Get the shapes and sizes of the images
      CALL UIMGID(IDD(1),DATTYP(1),NDIMS,DIMS,ISTAT)
                                                    
      IF(NDIMS.EQ.1) THEN
         DIMS(2)=1
      ENDIF
           
      NX=DIMS(1)
      NY=DIMS(2)
                
      DO N=1,NIM
        CALL UIMGID(IDD(N),DATTYP(N),NDIMS,DIMS,ISTAT)
                                                      
         IF(NDIMS.EQ.1) THEN
            DIMS(2)=1
         ENDIF
              
C Check they are OK
        IF(DIMS(1).NE.NX .OR.
     :     DIMS(2).NE.NY) THEN
           CALL UMSPUT('Data/PSF size/shape error',1,0,ISTAT)
           GO TO 99
        ENDIF
      ENDDO
           
C Check for powers of 2
      IF(.NOT.POWER2(DIMS(1)) .OR.
     :   .NOT.POWER2(DIMS(2))) THEN
         CALL UMSPUT('Image dimensions must be powers of 2',
     :               1,0,ISTAT)
         GO TO 99
      ENDIF
           
C Same for the PSFs
      DO N=1,NIM
        CALL UIMGID(IDP(N),DATTYP(N),NDIMS,DIMS,ISTAT)
                                                      
         IF(NDIMS.EQ.1) THEN
            DIMS(2)=1
         ENDIF
              
C Check they are OK
        IF(NX.NE.DIMS(1) .OR.
     :     NY.NE.DIMS(2)) THEN
           CALL UMSPUT('Data/PSF size/shape error',1,0,ISTAT)
           GO TO 99
        ENDIF
      ENDDO
           
C Allocate space for the data
      DO N=1,NIM
        CALL UDMGET(NX*NY,7,DATPNT(N),ISTAT)
        IF(ISTAT.NE.0) THEN
           CALL UMSPUT(
     :        'Unable to allocate memory for data array',
     :                 1,0,ISTAT)
           GO TO 99
        ENDIF
             
C Read the data into the memory
        CALL UIGS2D(IDD(N),1,NX,1,NY,MEMD(DATPNT(N)),ISTAT)
        IF(ISTAT.NE.0) THEN
           CALL UMSPUT(
     :        'Unable to read in data input image',
     :                 1,0,ISTAT)
           GO TO 99
        ENDIF
             
C Allocate space for the PSF
        CALL UDMGET(NX*NY,7,PSFPNT(N),ISTAT)
        IF(ISTAT.NE.0) THEN
           CALL UMSPUT(
     :        'Unable to allocate memory for PSF array',
     :                 1,0,ISTAT)
           GO TO 99
        ENDIF
             
C Read the PSF into the memory
        CALL UIGS2D(IDP(N),1,NX,1,NY,MEMD(PSFPNT(N)),ISTAT)
        IF(ISTAT.NE.0) THEN
           CALL UMSPUT(
     :        'Unable to read in PSF input image',
     :                 1,0,ISTAT)
           GO TO 99
        ENDIF
             
C Get working space arrays for the FFTs
        CALL UDMGET(NX*NY*2,7,PSFFFT(N),ISTAT)
        IF(ISTAT.NE.0) THEN
           CALL UMSPUT(
     :        'Unable to allocate memory for PSF FFT array',
     :                 1,0,ISTAT)
           GO TO 99
        ENDIF
      ENDDO
           
C Create the output data arrays
      CALL UCLGST('Decon',DECON,ISTAT)
      CALL UIMCRE(DECON,7,2,DIMS,IDDC,ISTAT)
      IF(ISTAT.NE.0) THEN
         CALL UMSPUT(
     :     'Unable to open output deconvolved image',
     :                 1,0,ISTAT)
         GO TO 99
      ELSE
         IF(VERBOSE) CALL UMSPUT(
     :     '-Created output (deconvolved) image: '//DECON,
     :     1,0,ISTAT)
      ENDIF
           
C Get the root name for the coadded output images
      CALL UCLGST('Coaddr',COADDR,ISTAT)
                                        
C Find the length of the string to avoid blanks in the names
      CALL LENSTR(COADDR,IST,IEND)
                                  
      DO N=1,NIM
         WRITE(COADD(N),'(a<iend-ist+1>,''_'',i1)')
     :      COADDR(IST:IEND),N
         CALL UIMCRE(COADD(N),7,2,DIMS,IDCOAD(N),ISTAT)
         IF(ISTAT.NE.0) THEN
            CALL UMSPUT('Unable to open output coadded image',
     :                    1,0,ISTAT)
            GO TO 99
         ELSE
             IF(VERBOSE) CALL UMSPUT(
     :        '-Created output (coadded) image: '//COADD(N),
     :        1,0,ISTAT)
         ENDIF
       ENDDO
            
89    CONTINUE
              
C Get the sub-sampling values in X and Y and check that
C the X and Y dimensions are multiples of them.
      CALL UCLGSI('Xsubsam',NBX,ISTAT)
      CALL UCLGSI('Ysubsam',NBY,ISTAT)
                                      
      IF(MOD(NX,NBX).NE.0 .OR.
     :   MOD(NY,NBY).NE.0) THEN
         CALL UMSPUT('Warning - arrays are not multiples'//
     :               ' of the sub-sampling',1,0,ISTAT)
      ENDIF
           
      IF(NBX.NE.1 .OR. NBY.NE.1) THEN
         SUBSAM=.TRUE.
      ELSE
         SUBSAM=.FALSE.
      ENDIF
           
C Prepare the arrays for the Numerical Recipes FFT routine
C and do the FFTs of the PSF and rotated PSF
C
C At this point we also check that the PSFs are normalised to
C a total of 1 - if they aren't we warn the user but continue
      DO N=1,NIM
         CALL TOTAL(MEMD(PSFPNT(N)),NX,NY,PT)
         IF(DABS(PT-1.0D0).GT.1.0E-5) THEN
            WRITE(CHARS,
     :       '(''Warning, PSF #'',i2,
     :       '' not normalised to a total of 1.0'',
     :       ''. sum is '',F15.13)') N,PT
            CALL UMSPUT(CHARS,1,0,ISTAT)
         ENDIF
              
         CALL DFILL(MEMD(PSFPNT(N)),NX,NY,MEMD(PSFFFT(N)))
         CALL DFOURN(MEMD(PSFFFT(N)),DIMS,2,+1)
                                               
C At this point we have finished with the PSFs so they can be closed
C and the memory freed.
         CALL UIMCLO(IDP(N),ISTAT)
         CALL UDMFRE(PSFPNT(N),7,ISTAT)
      ENDDO
           
C Create other arrays needed, workspace etc
        CALL UDMGET(NX*NY,7,W,ISTAT)
        IF(ISTAT.NE.0) THEN
           CALL UMSPUT(
     :        'Unable to allocate memory for working array',
     :                 1,0,ISTAT)
           GO TO 99
        ENDIF
             
        CALL UDMGET(NX*NY*2,7,WORK,ISTAT)
        IF(ISTAT.NE.0) THEN
           CALL UMSPUT(
     :        'Unable to allocate memory for workspace array',
     :                 1,0,ISTAT)
           GO TO 99
        ENDIF
             
        CALL UDMGET(NX*NY,7,REST,ISTAT)
        IF(ISTAT.NE.0) THEN
           CALL UMSPUT(
     :        'Unable to allocate memory for restored array',
     :                 1,0,ISTAT)
           GO TO 99
        ENDIF
             
C PHI arrays and DPS are only needed in the accelerated case
       IF(ACC)THEN
        DO N=1,NIM
           CALL UDMGET(NX*NY,7,PH(N),ISTAT)
           IF(ISTAT.NE.0) THEN
              CALL UMSPUT(
     :           'Unable to allocate memory for Phi array',
     :                    1,0,ISTAT)
              GO TO 99
           ENDIF
         ENDDO
              
         CALL UDMGET(NX*NY,7,DPS,ISTAT)
         IF(ISTAT.NE.0) THEN
           CALL UMSPUT(
     :        'Unable to allocate memory for delta Psi array',
     :                 1,0,ISTAT)
           GO TO 99
         ENDIF
       ENDIF
            
C Get space for and initialise the correction factor array
       CALL UDMGET(NX*NY,7,CORRFAC,ISTAT)
       IF(ISTAT.NE.0) THEN
          CALL UMSPUT('Unable to create corr fac image',
     :                    1,0,ISTAT)
          GO TO 99
       ENDIF
            
C Calculate the weights - just the normalised total counts in each image
       T=0.0D0
       DO N=1,NIM
          CALL TOTAL(MEMD(DATPNT(N)),NX,NY,WEIGHT(N))
          T=T+WEIGHT(N)
       ENDDO
            
       DO N=1,NIM
          WRITE(CHARS,'(''-Image '',i2,'' Total counts: '',
     :          F10.1,'' Weight: '',F5.3)') N,WEIGHT(N),WEIGHT(N)/T
          IF(VERBOSE) CALL UMSPUT(CHARS,1,0,ISTAT)
          WEIGHT(N)=WEIGHT(N)/T
       ENDDO
            
C Check whether we have a first estimate image, if so try to get it
      CALL UCLGSB('First',FIRST,ISTAT)
      IF(FIRST) THEN
         CALL UCLGST('FirstIm',FIM,ISTAT)
                                         
         CALL UIMOPN(FIM,1,FID,ISTAT)
         IF(ISTAT.NE.0) THEN
           CALL UMSPUT(
     :        'Unable to open first estimate image file',
     :                 1,0,ISTAT)
           GO TO 99
         ELSE
            IF(VERBOSE) CALL UMSPUT
     :         ('-First Estimate Image: '//FIM,1,0,ISTAT)
         ENDIF
              
         CALL UIGS2D(FID,1,NX,1,NY,MEMD(REST),ISTAT)
      ELSE
          
C Prepare the first estimate - just a flat image with the same total flux
          CALL FILCON(MEMD(REST),NX,NY,T/DFLOAT(NX*NY))
      ENDIF
           
C Now we can start the iterative procedure
      CALL UCLGSI('Niter',NITER,ISTAT)
                                      
      DO N=1,NITER
                  
         WRITE(CHARS,'(''--Starting iteration '',i5)') N
         IF(VERBOSE) CALL UMSPUT(CHARS,1,0,ISTAT)
                                                 
C The correction array has to be initialised each time
        CALL FILCON(MEMD(CORRFAC),NX,NY,0.0D0)
                                              
        DO M=1,NIM
                  
C We must first multiply the current estimate by the weight
C to get something with the same total flux as the data - this
C was omitted in earlier version and may have caused problems
C ___________________________________________________________
C
         CALL MULC(MEMD(REST),NX,NY,WEIGHT(M),MEMD(W))
                                                      
C Convolve the current estimate with the PSF to give Phi
         IF(ACC) THEN
          CALL DCONV(MEMD(W),NX,NY,MEMD(WORK),
     :              MEMD(PSFFFT(M)),MEMD(PH(M)),1)
         ELSE
          CALL DCONV(MEMD(W),NX,NY,MEMD(WORK),
     :              MEMD(PSFFFT(M)),MEMD(W),1)
         ENDIF
              
C Calculate and display the residual information to assist closeness
C of fit assessment
         IF(VERBOSE) THEN
          IF(ACC) THEN
             CALL RESINFO(MEMD(PH(M)),MEMD(DATPNT(M)),
     :                    NX,NY,RMSRES,RESMAX,IRMAX,JRMAX)
          ELSE
             CALL RESINFO(MEMD(W),MEMD(DATPNT(M)),
     :                    NX,NY,RMSRES,RESMAX,IRMAX,JRMAX)
          ENDIF
               
          WRITE(CHARS,'(''--Image #'',I2,'' RMS residual: '',D10.4,
     :                 '' Max Residual of '',D10.4,
     :                 '' at ('',I4'','',I4,'')'')')
     :             M,RMSRES,RESMAX,IRMAX,JRMAX
          CALL UMSPUT(CHARS,1,0,ISTAT)
         ENDIF
              
C If we are sub-sampling we must do so here.
      IF(SUBSAM) THEN
         IF(ACC) THEN
            CALL REBIN(MEMD(PH(M)),NX,NY,NBX,NBY)
         ELSE
            CALL REBIN(MEMD(W),NX,NY,NBX,NBY)
         ENDIF
      ENDIF
C Divide the data by this convolution
         IF(ACC) THEN
          CALL DIVIDE(MEMD(DATPNT(M)),MEMD(PH(M)),NX,NY,MEMD(W))
         ELSE
          CALL DIVIDE(MEMD(DATPNT(M)),MEMD(W),NX,NY,MEMD(W))
         ENDIF
              
C Convolve again, this time by the rotated PSF
         CALL DCONV(MEMD(W),NX,NY,MEMD(WORK),
     :             MEMD(PSFFFT(M)),MEMD(W),-1)
                                              
C Add the correction factor so deduced into the overall correction
C factor with the correct weighting
         CALL ADDW(MEMD(W),WEIGHT(M),NX,NY,MEMD(CORRFAC))
        ENDDO
             
C In the accelerated case calculate the maximum possible multiplication
C factor which is consistent with the non-negativity contraint
        IF(ACC) THEN
           CALL FIMAXF(MEMD(REST),MEMD(CORRFAC),NX,NY,FMAX)
           WRITE(CHARS,
     :   '(''--Max. speed-up possible'',
     :     '' within non-negativity (fmax): '',
     :      F10.5)') FMAX
           IF(VERBOSE) CALL UMSPUT(CHARS,1,0,ISTAT)
                                                   
C Work out the increment in Psi
           CALL PHINC(MEMD(REST),MEMD(CORRFAC),NX,NY,MEMD(DPS))
                                                               
C We now start the iteration loop for the Newton-Raphson search
                                                               
           IF(VERBOSE) THEN
             CALL UMSPUT('---Starting search for optimal speed'//
     :                 '-up factor...',1,0,ISTAT)
           ENDIF
                
           XL1=1.0D0
           INR=0
                
145        CONTINUE
           XLP=XL1
                  
           DRV1=0.0D0
           DRV2=0.0D0
           HHT=0.0D0
                    
           DO M=1,NIM
                     
C Convolve the increment in psi to get the increment in phi
              CALL DCONV(MEMD(DPS),NX,NY,MEMD(WORK),
     :                  MEMD(PSFFFT(M)),MEMD(W),1)
                                                  
C Multiply by the weight again to conserve the flux total
C (this is again a new addition)
              CALL MULC(MEMD(W),NX,NY,WEIGHT(M),MEMD(W))
                                                        
C If we are sub-sampling we must do so here again.
      IF(SUBSAM) THEN
         CALL REBIN(MEMD(W),NX,NY,NBX,NBY)
      ENDIF
           
C Calculate the derivatives needed for Newton-Raphson
              CALL DERIVS(MEMD(DATPNT(M)),MEMD(PH(M)),MEMD(W),
     :                    NX,NY,XL1,DD1,DD2,HH)
                                               
              HHT=HH+HHT
                        
C Add in the derivs with the correct weights
              DRV1=DRV1+DD1*WEIGHT(M)
              DRV2=DRV2+DD2*WEIGHT(M)
           ENDDO
                
           INR=INR+1
                    
C Calculate next estimate
           XL1=XL1-DRV1/DRV2
                            
           WRITE(CHARS,
     :       '(''---Iteration: '',I3,'' Speed up: '',F9.4,
     :         '' Log.Likelihood: '',
     :         D15.8)') INR,XL1,HHT
           IF(VERBOSE)  CALL UMSPUT(CHARS,1,0,ISTAT)
                                                    
C Check for negative likelihood gradient
           IF(INR.EQ.1 .AND.DRV1.LT.0.0D0) THEN
              CALL UMSPUT('Error - drv1<0 for inr=1',1,0,ISTAT)
              CALL UMSPUT('No acceleration this time round',1,0,
     :                    ISTAT)
              XL1=1.0D0
              GO TO 148
           ENDIF
                
C Check for speed-up too big for non-negativity
C In which case set the speed-up-factor at 70% of the max
C permitted factor
           IF(XL1.GT.FMAX) THEN
              XL1=1.0D0+0.7D0*(FMAX-1.0D0)
              IF(VERBOSE) CALL UMSPUT(
     :        '---Avoiding non-negativity'//
     :        ' (accel=1.0+0.7*(fmax-1))',
     :                    1,0,ISTAT)
              GO TO 148
           ENDIF
                
           IF(ABS(XL1-XLP).GT.0.02D0) GO TO 145
                                               
148        CONTINUE
           XL=XL1
           WRITE(CHARS,'(''---Speed up factor used: '',
     :           F10.5)') XL
           IF(VERBOSE) CALL UMSPUT(CHARS,1,0,ISTAT)
        ENDIF
             
C Multiply the last estimate by the total correction factor from all the
C images and hence obtain the next estimate
C We also remormalise at this stage to prevent accumulating
C normalisation errors
        CALL UPCORR(MEMD(REST),MEMD(CORRFAC),XL,NX,NY,
     :                MEMD(REST),T,VERBOSE)
      ENDDO
           
C Write out the results, deconvolved and then the coadded images
C which are the convolutions with the relavant PSFs
       CALL UIPS2D(IDDC,1,NX,1,NY,MEMD(REST),ISTAT)
       CALL UIMCLO(IDDC,ISTAT)
                              
       DO N=1,NIM
          CALL DCONV(MEMD(REST),NX,NY,MEMD(WORK),
     :              MEMD(PSFFFT(N)),MEMD(W),1)
                                              
C Multiply by the weight to get the correct scale
          CALL MULC(MEMD(W),NX,NY,WEIGHT(N),MEMD(W))
                                                    
          CALL UIPS2D(IDCOAD(N),1,NX,1,NY,MEMD(W),ISTAT)
          CALL UIMCLO(IDCOAD(N),ISTAT)
       ENDDO
            
99     CONTINUE
               
C Close all the images
      DO N=1,NIM
         CALL UIMCLO(IDD(N),ISTAT)
      ENDDO
           
C Free all the dynamic arrays
      DO N=1,NIM
         CALL UDMFRE(DATPNT(N),7,ISTAT)
         CALL UDMFRE(PSFFFT(N),7,ISTAT)
         IF(ACC) CALL UDMFRE(PH(N),7,ISTAT)
      ENDDO
           
      CALL UDMFRE(CORRFAC,7,ISTAT)
      CALL UDMFRE(WORK,7,ISTAT)
      CALL UDMFRE(REST,7,ISTAT)
      CALL UDMFRE(W,7,ISTAT)
      IF(ACC) CALL UDMFRE(DPS,7,ISTAT)
                                      
      CALL STSEPI    !M
                       
      END
         
      SUBROUTINE MULC(IN1,NX,NY,F,OUT)
C
C Just multiply one array by a constant
C
      IMPLICIT NONE
                   
      INTEGER NX,NY,I,J
      DOUBLE PRECISION IN1(NX,NY),OUT(NX,NY),F
                                              
      DO J=1,NY
         DO I=1,NX
            OUT(I,J)=IN1(I,J)*F
         ENDDO
      ENDDO
           
      RETURN
      END
         
      SUBROUTINE UPCORR(IN1,IN2,XL,NX,NY,OUT,T,VERBOSE)
C
C Apply the update, with a speed up factor
C
C Also renormalise to prevent accumulating errors
C
      IMPLICIT NONE
                   
      INTEGER NX,NY,I,J,ISTAT,NNEG
      DOUBLE PRECISION IN1(NX,NY),IN2(NX,NY)
      DOUBLE PRECISION T,V,OUT(NX,NY),XL,TOTAL
      CHARACTER*80 CHARS
      LOGICAL VERBOSE
                     
      TOTAL=0.0D0
      NNEG=0
            
C Work out the total counts allowing for negative number
C corrections
      DO J=1,NY
         DO I=1,NX
            V=IN1(I,J)*(1.0D0+XL*(IN2(I,J)-1.0D0))
            TOTAL=TOTAL+V
            IF(V.LT.0.0D0) THEN
               TOTAL=TOTAL-V
               OUT(I,J)=0.0D0
               NNEG=NNEG+1
            ELSE
               OUT(I,J)=V
            ENDIF
         ENDDO
      ENDDO
           
C Do the renormalisation
      DO J=1,NY
         DO I=1,NX
            OUT(I,J)=OUT(I,J)*T/TOTAL
         ENDDO
      ENDDO
           
      WRITE(CHARS,'(''--Renormalising restored image'',
     :'', (Factor: '',F20.18,'')'')')
     : T/TOTAL
      IF(VERBOSE) CALL UMSPUT(CHARS,1,0,ISTAT)
                                              
      IF(NNEG.GT.0) THEN
         WRITE(CHARS,'(''--A total of '',i8,
     :         '' negative points set to zero.'')') NNEG
         IF(VERBOSE) CALL UMSPUT(CHARS,1,0,ISTAT)
      ENDIF
           
      RETURN
      END
         
      SUBROUTINE DIVIDE(IN1,IN2,NX,NY,OUT)
C
C Just divide one array by another of the same size
C
      IMPLICIT NONE
                   
      INTEGER NX,NY,I,J
      DOUBLE PRECISION IN1(NX,NY),IN2(NX,NY),OUT(NX,NY)
                                                       
      DO J=1,NY
         DO I=1,NX
            IF(IN2(I,J).EQ.0.0D0) THEN
               OUT(I,J)=0.0D0
            ELSE
               OUT(I,J)=IN1(I,J)/IN2(I,J)
            ENDIF
         ENDDO
      ENDDO
           
      RETURN
      END
         
      SUBROUTINE FILCON(DATA,NX,NY,VAL)
C
C Fill up an array with a constant value
C
      IMPLICIT NONE
                   
      INTEGER NX,NY,I,J
      DOUBLE PRECISION DATA(NX,NY),VAL
                                      
      DO J=1,NY
         DO I=1,NX
            DATA(I,J)=VAL
          ENDDO
      ENDDO
           
      RETURN
      END
         
      SUBROUTINE TOTAL(DATA,NX,NY,TOT)
C
C Add up all elements in an image
C
      IMPLICIT NONE
                   
      INTEGER NX,NY,I,J
      DOUBLE PRECISION DATA(NX,NY),TOT
                                      
      TOT=0.0D0
      DO J=1,NY
         DO I=1,NX
             TOT=TOT+DATA(I,J)
         ENDDO
      ENDDO
           
      RETURN
      END
         
      SUBROUTINE ADDW(DATA,W,NX,NY,OUT)
C
C Add an image, weighted, to another
C
      IMPLICIT NONE
                   
      INTEGER NX,NY,I,J
      DOUBLE PRECISION DATA(NX,NY),W,OUT(NX,NY)
                                               
      DO J=1,NY
         DO I=1,NX
            OUT(I,J)=DATA(I,J)*W+OUT(I,J)
         ENDDO
      ENDDO
           
      RETURN
      END
         
      SUBROUTINE FIMAXF(PS,CF,NX,NY,FMAX)
C
C Find the maximum multiplying factor possible for acceleration
C which is consistent with the non-negativity contraint.
C
      IMPLICIT NONE
                   
      INTEGER NX,NY,I,J
      DOUBLE PRECISION PS(NX,NY),CF(NX,NY),FMAX,RT
                                                  
      FMAX=1.0D10
                 
      DO J=1,NY
         DO I=1,NX
             IF(PS(I,J).LT.1.0D-10 .OR.
     :          CF(I,J).EQ.1.0) GO TO 88
             RT=-1.0D0/(CF(I,J)-1.0D0)
             IF(RT.LT.0.0D0) GO TO 88
             IF(RT.LT.FMAX) FMAX=RT
 88          CONTINUE
          ENDDO
      ENDDO
           
      RETURN
      END
         
      SUBROUTINE DERIVS(PHT,PH,DPH,NX,NY,XL,DD1,DD2,HH)
C
C Calculate derivatives needed in Newton-Raphson search for
C optimal acceleration factor.
C
C Based on code by Leon Lucy
C
      IMPLICIT NONE
                   
      INTEGER NX,NY
      DOUBLE PRECISION PHT(NX,NY),PH(NX,NY),DPH(NX,NY)
      DOUBLE PRECISION XL
      DOUBLE PRECISION DD1,DD2,HH
                                 
      INTEGER I,J
      DOUBLE PRECISION DM,AA,DDPH,DPHT
                                      
      DD1=0.0D0
      DD2=0.0D0
      HH=0.0D0
              
C NB here we are ignoring negative values of dm - these
C may represent pathological situations which should be checked
C out.
      DO J=1,NY
         DO I=1,NX
            DDPH=DPH(I,J)
            DM=PH(I,J)+XL*DDPH
            IF(DM.GT.1.0D-20) THEN
               DPHT=PHT(I,J)
               AA=DDPH/DM
               DD1=DD1+DPHT*AA
               DD2=DD2-DPHT*AA*AA
               HH=HH+DPHT*DLOG(DM)
            ENDIF
         ENDDO
      ENDDO
           
      RETURN
      END
         
      SUBROUTINE PHINC(PS,CF,NX,NY,DPS)
C
C Calculate the increment in Psi, this is only needed in the
C accelerated case
C
      IMPLICIT NONE
                   
      INTEGER NX,NY
      DOUBLE PRECISION PS(NX,NY),CF(NX,NY),DPS(NX,NY)
      INTEGER I,J
                 
      DO J=1,NY
         DO I=1,NX
            DPS(I,J)=PS(I,J)*(CF(I,J)-1.0D0)
         ENDDO
      ENDDO
           
      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 RESINFO(IM1,IM2,NX,NY,RMSRES,RESMAX,
     :                   IRMAX,JRMAX)
C
C Calculate the RMS residual between two images and the
C largest residual.
C
      IMPLICIT NONE
                   
      INTEGER NX,NY,IRMAX,JRMAX
      DOUBLE PRECISION RMSRES,RESMAX,IM1(NX,NY),IM2(NX,NY)
                                                          
      INTEGER I,J
      DOUBLE PRECISION RES,T
                            
      T=0.0D0
      RESMAX=0.0D0
                  
      DO J=1,NY
         DO I=1,NX
            RES=IM1(I,J)-IM2(I,J)
            T=T+RES*RES
            IF(DABS(RES).GT.DABS(RESMAX)) THEN
               RESMAX=RES
               IRMAX=I
               JRMAX=J
            ENDIF
         ENDDO
      ENDDO
           
      RMSRES=DSQRT(T/DBLE(NX*NY))
                                 
      RETURN
      END
         
      SUBROUTINE REBIN(IN,NX,NY,NBX,NBY)
C
C Re-bin an array by replacing groups of pixels by their
C average. This is needed in the sub-sampling case.
C Flux is conserved and the operation is done 'in place'
C
      IMPLICIT NONE
                   
      INTEGER NX,NY,NBX,NBY
      DOUBLE PRECISION IN(NX,NY)
                                
      INTEGER I,J,K,L
      DOUBLE PRECISION T,DN
                           
      DN=DBLE(NBX*NBY)
                      
      DO J=1,NY,NBY
         DO I=1,NX,NBX
            T=0D0
            DO L=1,NBY
               DO K=1,NBX
                  T=T+IN(I+K-1,J+L-1)
               ENDDO
            ENDDO
            DO L=1,NBY
               DO K=1,NBX
                  IN(I+K-1,J+L-1)=T/DN
               ENDDO
            ENDDO
         ENDDO
      ENDDO
           
      RETURN
      END