program imrot
c-----------------------------------------------------------------------
c= IMROT - Derotate Q,U images given a specified Rotation Measure
c& dpr
c: image analysis
c+
c	IMROT forms new Q,U images corrected for a user-specified Rotation
c       measure. The rotation measure is applied across the entire map, and
c       affects all planes. Position angle is positive N -> E.
c@ in
c	The input Q, U images, respectively.
c	Wild card expansion is possibly supported. 
c@ out  
c	The input Q, U images, respectively.
c@ rm
c	Rotation measure in rad/m/m. Sense?
c--
c  History:
c    dpr  20aug01 Original version from impol.
c------------------------------------------------------------------------
      implicit none
c
      include 'maxdim.h'
      include 'maxnax.h'
      character version*(*)
      parameter (version = 'ImRot: version 20-Aug-2001')
cc
      real qline(maxdim), uline(maxdim), qoline(maxdim), 
     +  uoline(maxdim),  
     +  qepoch, uepoch,
     +  rm
      double precision qcdelt(maxnax), ucdelt(maxnax), 
     +  qcrval(maxnax), ucrval(maxnax),
     +  qcrpix(maxnax), ucrpix(maxnax)
c
      integer lq, lu, lqout(2), luout(2), 
     +  qsize(maxnax), usize(maxnax), qnaxis, unaxis,
     +  qstkax, ustkax, nin, nout
c
      character ins(2)*64, outs(2)*64,
     +  qout*64, uout*64, qin*64, uin*64,
     +  qctype(maxnax)*9, uctype(maxnax)*9, 
     +  bflag
c
      logical qflags(maxdim), 
     +  uflags(maxdim), qoflags(maxdim), uoflags(maxdim), 
     +  relax

c
      integer len1
      integer nkeys
      parameter (nkeys = 25)
      character keyw(nkeys)*8
c
      data keyw/     'obstime ','epoch   ','history ','instrume',
     +    'niters  ','object  ','restfreq','telescop','vobs    ',
     +    'obsra   ','obsdec  ','observer','cellscal',
     +    'bmaj    ','bmin    ','bpa     ','pbfwhm  ','lstart  ',
     +    'lstep   ','ltype   ','lwidth  ','vobs    ','pbtype  ',
     +    'btype   ','bunit   '/
      data lqout, luout  /2*0, 2*0/
c-------------------------------------------------------------------------
      call output (version)
      call output (' ')
c
c Get the inputs
c
      call keyini
      call mkeyf ('in', ins, 2, nin)
      call mkeya ('out', outs, 2, nout)
      call keyr ('rm', rm, 0.0)
      call keyfin
c
c Process the inputs
c
      if (nin.ne.2) 
     +  call bug ('f', 'You must give two input images')
      qin = ins(1)
      uin = ins(2)
      if (nout.ne.2)
     +  call bug ('f', 'You must give two output images')
      qout = outs(1)
      uout = outs(2)
c
c
      bflag = 'f'
      if (relax) bflag = 'w'
c
c Open the input images 
c

c
c Q
c
        call openin (bflag, maxdim, maxnax, qin, lq, qnaxis, qsize, 
     +     qepoch, qcrpix, qcdelt, qcrval, qctype, qstkax)
        if (qstkax.ne.0 .and. qcrval(qstkax).ne.2) call bug (bflag, 
     +     qin(1:len1(qin))//' does not appear to be a Q image')
c
c U
c
        call openin (bflag, maxdim, maxnax, uin, lu, unaxis, usize,
     +      uepoch, ucrpix, ucdelt, ucrval, uctype, ustkax)
        if (ustkax.ne.0 .and. ucrval(ustkax).ne.3) call bug (bflag, 
     +     uin(1:len1(uin))//' does not appear to be a U image')
c
c Compare images for consistency
c
        call chkdes (bflag, qin, uin, qnaxis, unaxis, qsize, usize, 
     +     qcrpix, ucrpix, qcdelt, ucdelt, qcrval, ucrval, qepoch,
     +     uepoch, qctype, uctype, qstkax, ustkax)
c
c Open output images.
c
        call openout (lq, qnaxis, qsize, nkeys, keyw,
     +     qcrval, qcrpix, qcdelt, qctype, qout, 
     +     'intensity', version, lqout(1))
        call openout (lu, unaxis, usize, nkeys, keyw,
     +     ucrval, ucrpix, ucdelt, uctype, uout, 
     +     'intensity', version, luout(1))
c
c Now compute and write out the output image(s)
c
        call polout (lq, lu, lqout, luout, qnaxis, qsize, 
     +   qcrpix, qcrval, qcdelt, qctype, rm,  qline, uline,
     +   qoline, uoline, qflags, 
     +   uflags, qoflags, uoflags)
c
c Close up
c
        call xyclose (lq)
        call xyclose (lu)
        call xyclose (lqout(1))
        call xyclose (luout(1))


c
      end
c
c
      subroutine getopt (debias, radians, relax, zero)
c----------------------------------------------------------------------
c     Decode options array into named variables.
c
c   Output:
c     debias    Debias the polarized intensity image
c     radians   Output position nagle in radians
c     relax     Warnings only for axis descriptor mismatches
c     zero      Output zeros rather than setting flagging mask
c-----------------------------------------------------------------------
      implicit none
c
      logical debias, radians, relax, zero
cc
      integer maxopt
      parameter (maxopt = 4)
c
      character opshuns(maxopt)*7
      logical present(maxopt)
      data opshuns /'bias', 'radians', 'relax', 'zero'/
c-----------------------------------------------------------------------
      call options ('options', opshuns, present, maxopt)
c
      debias  = .not.present(1)
      radians = present(2)
      relax   = present(3)
      zero    = present(4)
c
      end
c
c
      subroutine openin (bflag, maxdim, maxnax, in, lun, naxis, size,
     +                   epoch, crpix, cdelt, crval, ctype, stkax)
c-----------------------------------------------------------------------
c     Open an image and return some information about it
c
c  Input
c    bflag      Bug flag
c    maxdim     Maximum size a row can be
c    maxnax     Maximum number of axes image can have
c    in         Image name
c  Output
c    lun        Handle
c    naxis      Number of axes
c    size       Size of each axis
c    epoch      EPoch of image
c    crpix      Refernce pixels
c    cdelt      Increments
c    crval      Reference values
c    ctype      Axis types
c    stkax      Stokes axis
c-----------------------------------------------------------------------
      implicit none
c
      integer maxdim, maxnax, lun, naxis, size(maxnax), stkax
      double precision cdelt(maxnax), crval(maxnax), crpix(maxnax)
      real epoch
      character*(*) ctype(maxnax), in, bflag*1
cc
      integer len1, i
      character*80 aline
c-----------------------------------------------------------------------     
      call xyopen (lun, in, 'old', maxnax, size)
      call rdhdi (lun, 'naxis', naxis, 0)
      if (naxis.eq.0) then
        aline = in(1:len1(in))//' has zero dimensions !!'
        call bug ('f', aline)
      end if
c
      if (size(1).gt.maxdim) then
        aline = 'First dimension of '//in(1:len1(in))//
     +             ' too large for storage'
        call bug ('f', aline)
      end if
      call hedinf (lun, naxis, size, epoch, crpix, cdelt, crval, ctype)
c
      stkax = 0
      do i = 1, naxis
        if (ctype(i).eq.'STOKES') stkax = i
      end do
      if (stkax.eq.0) then
        aline = 'Could not find Stokes axis in '//in(1:len1(in))
        call bug (bflag, aline)
      end if
c
      end
c
c
      subroutine hedinf (lun, naxis, size, epoch, crpix, cdelt,
     +                   crval, ctype)
c------------------------------------------------------------------------
c     Get some header keywords from the image associated with LUN
c 
c     Input
c       lun      Handle of image
c       naxis    Number of dimensions in image
c       size     Size of each axis
c     Output
c       epoch    Epoch of image
c       crpix    Array of image reference pixels
c       cdelt    Array of image increments (natural inits; rad)
c       crval    Array of image reference values (natural units)
c       ctype    Array of image axis types
c
c------------------------------------------------------------------------
      implicit none
c
      integer lun, naxis, size(naxis)
      real epoch
      double precision cdelt(naxis), crval(naxis), crpix(naxis)
      character*(*) ctype(naxis)
cc
      integer i
      character str*1, itoaf*1
c---------------------------------------------------------------------
      do i = 1, naxis
        str = itoaf(i)
c
        call rdhdd (lun, 'crpix'//str, crpix(i), dble(size(i))/2.0d0)
        call rdhdd (lun, 'cdelt'//str, cdelt(i), 1.0d0)
        call rdhda (lun, 'ctype'//str, ctype(i), ' ')
        call rdhdd (lun, 'crval'//str, crval(i), 0.0d0)
      end do
      call rdhdr (lun, 'epoch', epoch, 0.0)
c
      end 
c
c
      subroutine chkdes (bflag, im1, im2, naxis1, naxis2, size1, size2,
     +   crpix1, crpix2, cdelt1, cdelt2, crval1, crval2, epoch1, 
     +   epoch2, ctype1, ctype2, stkax1, stkax2)
c-----------------------------------------------------------------------
c     Compare axis descriptors 
c
c  Input:
c   im1,2        Images
c   naxis1,2     Number of axes
c   size1,2      Sizes of each dimension
c   crpix1,2     Reference pixels
c   cdelt1,2     Increments
c   crval1,2     Refernce values
c   ctype1,2     types of axes
c   epoch1,2     Epochs
c   stkax1,2     Stokes axis
c  Output
c   stkax        Stokes axis
c-----------------------------------------------------------------------
      implicit none
c
      integer naxis1, naxis2, size1(*), size2(*), stkax1, stkax2
      character*(*) im1, im2, ctype1(*), ctype2(*), bflag
      real epoch1, epoch2
      double precision crval1(*), crval2(*), cdelt1(*), cdelt2(*),
     +  crpix1(*), crpix2(*)
cc
      integer k, l1, l2, len1
      character line*130
c-----------------------------------------------------------------------
      l1 = len1(im1)
      l2 = len1(im2)
c
      if (epoch1.ne.epoch2) then
        line = 'Unequal epochs for images '//im1(1:l1)//' & '//im2(1:l2)
        call bug (bflag, line)
      end if
c
      if (naxis1.ne.naxis2) then
        line = 'Unequal number dimensions for images '//
     +         im1(1:l1)//' & '//im2(1:l2)
        call bug (bflag, line)
      end if
c
      do k = 1, min(naxis1,naxis2)
        if (size1(k).ne.size2(k)) then
          write (line, 10) im1(1:l1), im2(1:l2), k
10        format ('Unequal sizes for images ', a, ' & ', a, 
     +            ' on axis ', i1)
          call bug (bflag, line)
        end if
c
        if (ctype1(k).ne.ctype2(k)) then
          write (line, 20) im1(1:l1), im2(1:l2), k
20        format ('Unequal ctype for images ', a, ' & ', a, 
     +            ' on axis ', i1)
          call bug (bflag, line)
        end if
c
        call chkds2 (bflag, 'crpix', k, im1(1:l1), im2(1:l2), 
     +               crpix1(k), crpix2(k))
        call chkds2 (bflag, 'cdelt', k, im1(1:l1), im2(1:l2), 
     +               cdelt1(k), cdelt2(k))
        if (k.ne.stkax1 .or. k.ne.stkax2)
     +    call chkds2 (bflag, 'crval', k, im1(1:l1), im2(1:l2), 
     +                 crval1(k), crval2(k))
      end do
c
      end
c
c
      subroutine chkds2 (bflag, type, iaxis, im1, im2, des1, des2)
c-----------------------------------------------------------------------
c     Compare an axis descriptor from two images
c
c  Input:
c    type    Type fo descriptor
c    iaxis   Axis number
c    im1,2   Images
c    des1,2  Descriptors
c
c-----------------------------------------------------------------------
      implicit none
c
      character*(*) type, im1, im2, bflag
      integer iaxis
      double precision des1, des2
cc
      character line*130
c-----------------------------------------------------------------------
      if (abs(des1-des2).gt.0.01*max(abs(des1),abs(des2))) then
        write (line, 10) type, im1, im2, iaxis
10      format ('Unequal ', a, ' for images ', a, ' & ', a, 
     +          ' on axis ', i1)
        call bug (bflag, line)
      end if
c
      end
c
c
      subroutine openout (lin, naxis, size, nkeys, keyw, crval, crpix,
     +                    cdelt, ctype, out, btype, version, lout)
c-----------------------------------------------------------------------
c     Open an output image, copy header keywords across and write
c     the history
c
c  Input
c   lin    Image to copy keuwrods from
c   naxis  Number of axes
c   size   Size of axes
c   nkeys  Number of header keywords to copy
c   keyw   Keywords
c   crval  Refernce values
c   crpix  Reference pixels
c   cdelt  Increments
c   ctype  Axis types
c   out    Name of output image
c   btype  'intensity'
c   versionVersion of this program
c  Output
c   lout   Handle for output image
c
c-----------------------------------------------------------------------
      implicit none
c
      integer lin, lout, naxis, size(naxis), nkeys
      double precision crval(naxis), cdelt(naxis), crpix(naxis)
      character*(*) keyw(nkeys), out, version, ctype(naxis), btype*(*)
cc
      integer i
      character itoaf*1, istr*1, aline*80
c-----------------------------------------------------------------------
      call xyopen (lout, out, 'new', naxis, size)
      do i = 1, nkeys
        call hdcopy (lin, lout, keyw(i))
      end do
c
c Do these separately because we had to strip the Stokes
c axis from the input image
c
      do i = 1, naxis
        istr = itoaf(i)
        call wrhdd (lout, 'crval'//istr, crval(i))
        call wrhdd (lout, 'crpix'//istr, crpix(i))
        call wrhdd (lout, 'cdelt'//istr, cdelt(i))
        call wrhda (lout, 'ctype'//istr, ctype(i))
      end do

c
      call hisopen  (lout, 'append')
      aline = 'IMROT: Miriad '//version
      call hiswrite (lout, aline)
      call hisinput (lout, 'IMROT')
      call hisclose (lout)
c
      end
c
      subroutine polout (lq, lu, lqout, luout, naxis,
     *   size, crpix, crval, cdelt, ctype,  rm,
     *   qline, uline, qoline, uoline,  qflags,
     *   uflags, qoflags,uoflags)
c-----------------------------------------------------------------------
c     Compute the new Q and U images
c
c-----------------------------------------------------------------------
      implicit none
c
      integer lq, lu, lqout(2), luout(2), naxis, 
     *  size(naxis)
      real qline(*), uline(*), qoline(*), uoline(*), rm
      double precision crval(naxis), cdelt(naxis), crpix(naxis)
      logical qflags(*), uflags(*), qoflags(*), uoflags(*)
      character*(*) ctype(naxis)
cc
      include 'mirconst.h'
      double precision r2d
      parameter (r2d = 180.0d0/dpi)
c
      integer i, j, k, frqax
      double precision fac
      real freq, parot
      character ustr*8, aline*80
c-----------------------------------------------------------------------
      fac = 1.0
      ustr = ' radians'

c
c Find frequency axis if rotating position angles back
c
      if (rm.ne.0.0) then
        frqax = 0
        do i = 1, naxis
          if (index(ctype(i),'FREQ').ne.0) frqax = i
        end do
c
        if (frqax.eq.0) call bug ('f', 
     +    'Could not find frequency axis with which to apply RM')
        if (frqax.le.2) call bug ('f',
     +    'Frequency axis is either 1 or 2.  These should be spatial')
c
        if (frqax.gt.3 .or. size(frqax).eq.1) then
c
c Find frequency of pixel one
c
          freq = (1.0 - crpix(frqax))*cdelt(frqax) + crval(frqax)
          freq = freq * 1.0e9
          parot = rm * (dcmks / freq)**2
          write (aline, 10) fac*parot, ustr
10        format ('Rotating position angles back by ', 1pe11.4, a)
          call output (aline)
        end if
      else
        parot = 0.0
      end if
c
c Loop over planes
c
      do k = 1, size(3)
        if (rm.ne.0.0 .and. frqax.eq.3 .and. size(frqax).gt.1) then
          freq = 1.0e9 * ((real(k)-crpix(3))*cdelt(3) + crval(3))
          parot = rm * (dcmks / freq)**2 
        end if
c
c Set planes
c
        call xysetpl (lq, 1, k)
        call xysetpl (lu, 1, k)
        call xysetpl (lqout(1), 1, k)
        call xysetpl (luout(1), 1, k)
c
c Read lines of data
c
        do j = 1, size(2)
          call xyread  (lq, j, qline)
          call xyflgrd (lq, j, qflags)
          call xyread  (lu, j, uline)
          call xyflgrd (lu, j, uflags)
c
c Work out everything possible from this row. By default, snclip(2)=0
c so ISNR=1 means no I based blanking by default
c
          do i = 1, size(1)
c
c Init all output arrays with zeros and bad flags
c
            call allblnk (qoline(i), qoflags(i), uoline(i), uoflags(i))
c
c See what we can validly work out
c
            if ( qflags(i) .and. uflags(i) ) then

c Work out all the output quantities here. 

                qoline(i) = qline(i)*cos(parot*2.0) + 
     *                        uline(i)*sin(parot*2.0)
                uoline(i) = -1.0*qline(i)*sin(parot*2.0) + 
     *                        uline(i)*cos(parot*2.0)
                qoflags(i) = .true.
                uoflags(i) = .true.
            end if
          end do

c
c Write out images
c
          if (lqout(1).ne.0) then
            call xywrite (lqout(1), j, qoline)
            call xyflgwr (lqout(1), j, qoflags)
          end if
c
          if (luout(1).ne.0) then
            call xywrite (luout(1), j, uoline)
            call xyflgwr (luout(1), j, uoflags)
          end if
        end do
      end do
      end
c
c
      subroutine allblnk (qp, qf, up, uf)
      implicit none
      real qp, up
      logical qf,uf
c
      qp = 0.0
      qf = .false.
      up = 0.0
      uf = .false.
c
      end
c
c
c