Merge pull request #13 from csiro/fits-regrid-maxdim-fixes

Some small fixes to fits and regrid.
Generalise maxdim parameter handling - move all to CMakeLists.txt

Author: mwiering

CMakeLists.txt

@@ -29,17 +29,20 @@ if(LARGEN)
   set(MAXCHAN 4096)
   set(MAXAVER 20000000)
   set(MAXWIN 16)
-  set(MAXWIDE 18)
   set(MAXSRC 36)
 else()
   set(MAXBUF 1024*1024*32)
   set(MAXANT 50)
   set(MAXCHAN 70000)
   set(MAXAVER 2000000)
   set(MAXWIN 60)
-  set(MAXWIDE 18)
   set(MAXSRC 1024)
 endif()
+set(MAXDIM 65536)
+set(MAXWIDE 18)
+set(MAXFBIN 16)
+set(MAXPNT 20000)
+
 set(MIRROOT ${CMAKE_INSTALL_PREFIX})
 
 # Preprocessing flags

inc/maxdim.h.in

@@ -1,29 +1,18 @@
 C=======================================================================
 C-----------------------------------------------------------------------
 C     Size of an INTEGER array used to implement a memory heap.  This
-C     array is the sole variable in blank COMMON in Miriad.  Trial-and-
-C     error compilations on an x86-64 system with gcc/g77 show that the
-C     limit on MAXBUF for which most tasks build successfully is
-C     1073741823 (2**30 - 1) which corresponds to 4GiB.  With MAXDIM
-C     less than 32768 the limit for which all tasks build successfully
-C     is about 260000000; unsuccessful links produce messages about
-C     truncated relocations, imom being the worst offender.
+C     array is the sole variable in blank COMMON in Miriad.
 C     The default value allocates 128MiB (for normal 4-byte INTEGERs).
       INTEGER   MAXBUF
       PARAMETER(MAXBUF = @MAXBUF@)
 
 C     Maximum image axis length.  Array dimensions are typically a few
 C     times MAXDIM (never MAXDIM**2) so MAXDIM is associated with a much
-C     smaller allocation of static memory than MAXBUF.  Thus the default
-C     value of MAXDIM is quite generous.
-C     The current value of MAXDIM is the largest value that keeps the
-C     memory allocation to store a single image plane below the current
-C     limit of 8GB. Note that invert is limited to 32768, the highest
-C     power of 2 below the limit.
+C     smaller allocation of static memory than MAXBUF.
       INTEGER   MAXDIM
-      PARAMETER(MAXDIM = 32768)
+      PARAMETER(MAXDIM = @MAXDIM@)
 
-C     Maximum number of antennas (ATA=64).
+C     Maximum number of antennas.
       INTEGER   MAXANT
       PARAMETER(MAXANT = @MAXANT@)
 
@@ -45,9 +34,9 @@ C     Maximum number of wideband channels.
 
 C     Maximum number of frequency bins in gains table
       INTEGER   MAXFBIN
-      PARAMETER(MAXFBIN = 16)
+      PARAMETER(MAXFBIN = @MAXFBIN@)
 
 C     Maximum number of mosaic pointings
       INTEGER   MAXPNT
-      PARAMETER(MAXPNT = 20000)
+      PARAMETER(MAXPNT = @MAXPNT@)
 C=======================================================================

inc/maxdimc.h.in

@@ -25,7 +25,7 @@
 #define MAXBUF  16777216
 
 /* Maximum number of pixels on an image axis. */
-#define MAXDIM  32768
+#define MAXDIM  @MAXDIM@
 
 /* Maximum number of image axes. */
 #define MAXNAX  7
@@ -46,4 +46,4 @@
 #define MAXWIDE @MAXWIDE@
 
 /* Maximum number of mosaic pointings. */
-#define MAXPNT 20000
+#define MAXPNT @MAXPNT@

src/prog/fits.for

@@ -1652,8 +1652,6 @@ c
         call output('  Using antenna table information')
         nconfig = nconfig + 1
         call ftabInfo(lu,'STABXYZ',type,units,n,nxyz)
-        if (nconfig.gt.MAXCONFG)
-     *    call bug('f','Too many array configurations')
         if (nxyz.ne.3 .or. n.le.0 .or. type.ne.'D')
      *    call bug('f','Something is screwy with the antenna table')
         if (n.gt.MAXANT) call bug('f','Too many antennas for me')
@@ -1728,6 +1726,10 @@ c
         emok = emok .and. .not.badmnt
 
         call ftabNxt(lu,'AIPS AN',found)
+        if (nconfig.ge.MAXCONFG) then
+          call bug('w','Too many array configurations')
+          found = .false.
+        endif
       enddo
 c
 c  If no antenna table was found, try for an OB table!

src/prog/regrid.for

@@ -25,7 +25,7 @@ c       (WCS) standard as defined in "WCS Paper II",
 c       Calabretta & Greisen (2002), A&A 395, 1077.
 c@ in
 c       The input image name.  In the first instance, coordinate
-c       descriptors for the output image are copied from the input.  
+c       descriptors for the output image are copied from the input.
 c       They may then be overridden by other parameters as described
 c       below.  No default.
 c@ out
@@ -130,7 +130,7 @@ c           AIR  Airy (1,D)
 c         Cylindricals:
 c           CYP  Cylindrical perspective (2,G|D)
 c           CEA  Cylindrical equal area (1,E,G)
-c           CAR  Plate carrée (aka Cartesian) (0,G) - please note that
+c           CAR  Plate Carrée (aka Cartesian) (0,G) - please note that
 c                this is NOT the same as the simple linear system used
 c                previously unless the reference coordinates (CRVAL)
 c                are (0,0)
@@ -200,7 +200,7 @@ c@ options
 c       Extra processing options that alter the axis description defined
 c       by the template image, axis descriptors, or input image.
 c       Several can be given, separated by commas, with minimum-match.
-c         altprj    Interpret a CAR (plate carée) projection in the
+c         altprj    Interpret a CAR (Plate Carrée) projection in the
 c                   input ot template image as a simple linear
 c                   coordinate system with an additional 1/cos(lat0)
 c                   scaling factor applied when computing the longitude,
@@ -232,6 +232,11 @@ c@ tol
 c       Interpolation tolerance.  Tolerate an error of the specified
 c       amount in converting pixel locations in the input to the output.
 c       Must be less that 0.5.  The default is 0.05.
+c@ div
+c       Number of divisions of the coordinate axes to use with offset
+c       option. Some projections, like TAN, require finer division near
+c       discontinuities. Default is 10 on first 3 axes. Suggested
+c       maximum: 200. Specify up to 3 values.
 c
 c--
 c  History:
@@ -255,7 +260,7 @@ c-----------------------------------------------------------------------
      *          maxc(3), maxv(3), minc(3), minv(3), n, naxes,
      *          nAxIn(MAXNAX), nAxOut(MAXNAX), nAxTem(MAXNAX), nblank,
      *          nBuf(3), ndesc, nIAxes, npv, nTAxes, nxy, off(3),
-     *          offset, order(3), xyzero
+     *          offset, order(3), xyzero, div(3), ndiv
       ptrdiff   xv, yv, zv, valid, rBuf, lBuf
       real      tol, fblank
       double precision cdelt, crpix, crval, desc(4,MAXNAX), latpol,
@@ -319,6 +324,15 @@ c     Get the input parameters.
       if (tol.lt.0.0 .or. tol.ge.0.5)
      *  call bug('f','Invalid value for the tol parameter')
 
+      if (doOff) then
+        call mkeyi('div',div,3,ndiv)
+        do i=ndiv+1,3
+          div(i) = 10
+        enddo
+        if (div(1).le.0.or.div(2).le.0.or.div(3).le.0)
+     *    call bug('f',"Invalid div values")
+      endif
+
       call keyfin
 
 c     Open the input dataset.
@@ -515,7 +529,7 @@ c     Set up output celestial coordinates.
       call setCel(lIn, cOut, doDesc, doEqEq, doGalEq)
 
 c     Set up offset coordinates.
-      if (doOff) call doOffset(lIn,nAxIn,cOut,nAxOut)
+      if (doOff) call doOffset(lIn,nAxIn,cOut,nAxOut,div)
 
 c     Initialise coordinate conversions.
       call pCvtInit(cOut,lIn)
@@ -848,29 +862,27 @@ c       Put coordinates in the normal orientation, i.e. with north up.
 
 c***********************************************************************
 
-      subroutine doOffset(lIn,nAxIn,lOut,nAxOut)
+      subroutine doOffset(lIn,nAxIn,lOut,nAxOut,div)
 
-      integer lIn,lOut,nAxIn(3),nAxOut(3)
+      integer lIn,lOut,nAxIn(3),nAxOut(3),div(3)
 c-----------------------------------------------------------------------
       include 'maxdim.h'
 
-      integer NV
       double precision TOL
-      parameter (NV = 10, TOL = 0.49d0)
+      parameter (TOL = 0.49d0)
 
-      logical first, valid(NV,NV,NV), warned, weird(3)
+      logical first, valid(div(1),div(2),div(3)), warned, weird(3)
       integer i, j, k, l, maxv(3), minv(3), nv1, nv2, nv3
-      double precision crpix, in(3), out(3,NV,NV,NV)
+      double precision crpix, in(3), out(3,div(1),div(2),div(3))
 
       external  itoaf
       character itoaf*2
 c-----------------------------------------------------------------------
       call pCvtInit(lIn,lOut)
 
-      nv1 = min(max(3,nAxIn(1)),NV)
-      nv2 = min(max(3,nAxIn(2)),NV)
-      nv3 = min(max(3,nAxIn(3)),NV)
-
+      nv1 = min(max(3,nAxIn(1)),div(1))
+      nv2 = min(max(3,nAxIn(2)),div(2))
+      nv3 = min(max(3,nAxIn(3)),div(3))
       first = .true.
       do k = 1, nv3
         do j = 1, nv2
@@ -1161,7 +1173,7 @@ c-----------------------------------------------------------------------
               zv(i,j) = real(out(3))
             else
               zv(i,j) = real(in(3))
-            endif            
+            endif
           endif
         enddo
       enddo