subroutine images(filey)
C
C Write a raster*8 image file in HDF format, adding 
C subsequent images to the same file.  In 3D write
C 3 seperate files.  For non-cartesian grids, interpolate
C onto a cartesian grid of dimensions ig x jg
C
C##########################################################
C
      include 'global.h'
      include 'zone.h'
      character char*1, filex*8, filey*8, filez*8
      integer i, j, k, d8aimg, isdn
      real x(ig),y(jg)
      real dcon
C
      character*1 image(ii,jj), imagec(ig,jg)
c
      dcon = 253. / (cmax-cmin)
      k = 1
c
      if(ndim.eq.2) then
C
      if(ngeomy.lt.3) then
C
C just load the density into a character array
      do 10 j=1,jmax
      do 10 i=1,imax
      scrtch = (zrho(i,j,k) - cmin) * dcon
      ival = int(scrtch)
      image (i,jmax+1-j) = char(min(max(1,ival),253))
 10   continue
      isdn = d8aimg(filey,image,ii,jj,0)
C
      else
C          interpolate the density onto a cartesian grid of ig X jg points
c
c     construct grid coordinates
c
      xrght = zxa(imax)
      xlft  = 1.0e-10
      ytop  = xrght 
      ybot  = -ytop
c
      delx = (xrght - xlft) / ig
      dely = (ytop  - ybot) / jg
      do 24 i=1,ig
        x(i) = xlft + real(i-1) * delx
 24   continue
      do 25 j=1,jg
        y(j) = ybot + real(j-1) * dely
 25   continue
c
      do 28 j=1,jg
      do 28 i=1,ig
         radius = sqrt(y(j)**2 + x(i)**2)
         theta  = .5*pi - atan(y(j)/x(i))
         jpn = jpa(i,j)
         ipn = ipa(i,j)
         jmn = jpn - 1
         imn = ipn - 1
         if(ipn.eq.-1) then
            imagec(i,jg+1-j) = char(2)
         else
            c1 = (radius - zxa(ipn)) / (zxa(imn) - zxa(ipn))
            c2 = (radius - zxa(imn)) / (zxa(ipn) - zxa(imn))
            c3 = (theta  - zya(jpn)) / (zya(jmn) - zya(jpn))
            c4 = (theta  - zya(jmn)) / (zya(jpn) - zya(jmn))
            a1 = zrho(imn,jmn,k)*c1 + zrho(ipn,jmn,k)*c2
            a2 = zrho(imn,jpn,k)*c1 + zrho(ipn,jpn,k)*c2
            scrtch = ((a1*c3+a2*c4) - cmin) * dcon
            ival   = int(scrtch)
            imagec(i,jg+1-j) = char(min(max(1,ival),253))
         endif
 28   continue
      isdn = d8aimg(filey,imagec,ig,jg,0)
      endif
C
      else if (ndim .eq. 3) then
C Output 3 2d byte arrays of density for making a movie
C This 3D version only works for cubic cartesian grids (imax=jmax=kmax)
C
      filex = filey(1:7) // 'x'
      filez = filey(1:7) // 'z'
C
C Write out movie frame in XY plane at k = kmid
      kmid = kmax/2
      do 30 j=1,jmax
      do 30 i=1,imax
      scrtch = (zrho(i,j,kmid) - cmin) * dcon
      ival = int(scrtch)
      image (i,jmax+1-j) = char(min(max(1,ival),253))
 30   continue
      isdn = d8aimg(filez,image,ii,jj,0)
C
C Write out movie frame in YZ plane at i = imid
      imid = imax/2
      do 40 j=1,jmax
      do 40 k=1,kmax
      scrtch = (zrho(imid,j,k) - cmin) * dcon
      ival = int(scrtch)
      image (j,kmax+1-k) = char(min(max(1,ival),253))
 40   continue
      isdn = d8aimg(filex,image,ii,jj,0)
C
C Write out movie frame in ZX plane at j = jmid
      jmid = jmax/2
      do 50 k=1,kmax
      do 50 i=1,imax
      scrtch = (zrho(i,jmid,k) - cmin) * dcon
      ival = int(scrtch)
      image (j,kmax+1-k) = char(min(max(1,ival),253))
 50   continue
      isdn = d8aimg(filey,image,ii,jj,0)
C
      endif
C
Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
      return
      end