c
c ----------------------------------------------------------------------
c
c The routines CCP4MAPHEAD, CCP4MAPIN and CCP4MAPOUT were written
c by Kevin Cowtan. As of version 3.5 of CCP4 these routines will
c be placed in the library maplib.f, however they are also examples
c of using the maplib routines.
c
c This program will read a complete map into an array, automatically
c sorted as X,Y,Z. The map array is then manipulated in the subroutine
c 'mapmanipulate' (in this example the values in the map are inverted,
c but users can easily alter/add to the code to perform more complicated
c manipulations). The modified map is then written out to a second file.
c
c It should be run as:
c
c map_rw MAPIN foo1.map MAPOUT foo2.map
c
c Nb there are no keywords; see the parsing templates for examples
c of adding keywords.
c
c (Thanks to Ho-Leung Ng for fixing some problems with this template!)
c
      program map_rw
c
      implicit none
c
c     ..Parameters ..
      INTEGER MAXMAP
      PARAMETER (MAXMAP=1000000)
c
c     ..Variables to store the map info
      INTEGER nspgrp,nu,nv,nw,nu1,nv1,nw1,nu2,nv2,nw2,
     +        mapsize
      REAL cell(6)
      REAL my_map(maxmap)
      CHARACTER*120 title
c
c CCP4 initialisation
      call ccpfyp
c
c Get header info for map
c
      write(6,FMT=9100)
      write(6,fmt=9000)
      write(6,FMT=9100)
 9000 format('*****    RETRIEVING HEADER INFORMATION...   *****')
c
      call ccp4maphead('MAPIN',nspgrp,cell,nu,nv,nw,nu1,nv1,nw1,nu2,
     +                 nv2,nw2)
c
c Check that the map will fit into the array
c
      mapsize = (nu2-nu1)*(nv2-nv1)*(nw2-nw1)
c
      if (mapsize.gt.maxmap) then
        write(6,fmt='(a,i10,a)') 'Warning: increase parameter MAXMAP to'
     +                           ,mapsize,' and recompile'
        call ccperr (1, 'Map array too small')
      end if
c
c Write some information from the header to the standard output
c
      write(6,FMT=9100)
      write(6,FMT=9010)nspgrp
      write(6,FMT=9020)cell
      write(6,FMT=9030)nu1,nu2,nv1,nv2,nw1,nw2
c
 9010 format(/,'Spacegroup number ',I4)
 9020 format(/,'Cell dimensions:',3F8.3,' A',/,
     +         'Cell angles    :',3F8.1,' deg.')
 9030 format(/,'Size of map:     nu1  nu2  nv1  nv2  nw1  nw2',
     +       /,15X,6I5,/)
c
c Read in the map
c
      write(6,FMT=9100)
      write(6,FMT=9040)
      write(6,FMT=9100)
 9040 format('*****        READING IN THE MAP ...         *****')
c
      call ccp4mapin('MAPIN',title,my_map,nu1,nv1,nw1,nu2,nv2,nw2)
c
c Write more information to standard output
c
      write(6,FMT=9100)
      write(6,FMT=9050)title
 9050 format(/,'Map title:',//,A120)
c
c ********************************************************
c **** Put code in MAPMANIPULATE to use/alter the map ****
c ********************************************************
c
      call mapmanipulate(my_map,nu1,nv1,nw1,nu2,nv2,nw2)
c
c Write out the modified map
c
      write(6,FMT=9100)
      write(6,FMT=9060)
      write(6,FMT=9100)
 9060 format('*****      WRITING OUT THE NEW MAP ...      *****')
 9100 format('*************************************************')
c
      call ccp4mapout('MAPOUT',title,my_map,nspgrp,cell,
     +                nu,nv,nw,nu1,nv1,nw1,nu2,nv2,nw2)
c
c End the program cleanly
c
      call ccperr(0,'Normal termination')
      end
c
c ----------------------------------------------------------------------
c
c_BEGIN_MAPMANIPULATE
c
      subroutine mapmanipulate(map,nu1,nv1,nw1,nu2,nv2,nw2)
c     =====================================================
c
c---- The map held in map(nu1:nu2,nv1:nv2,nw1:nw2) can be manipulated
c     within this subroutine. In this example the values in the map
c     are simply scaled by a factor of -1.
c
c_END_MAPMANIPULATE
c
      implicit none
c
c   ..Subroutine arguments
      integer nu1,nv1,nw1,nu2,nv2,nw2
      real map(nu1:nu2,nv1:nv2,nw1:nw2)
c
c   ..Local variables
      integer ic,jc,kc
c
c *********************************************
c **** Put code here to manipulate the map ****
c *********************************************
c
c In this example the map is inverted, i.e. scaled by a factor of -1.
c
      do 10 ic = nu1,nu2
        do 20 jc = nv1,nv2
          do 30 kc = nw1,nw2
c
            map(ic,jc,kc) = -map(ic,jc,kc)
c
 30       continue
 20     continue
 10   continue
c
      return
c
      end
c
c

c ----------------------------------------------------------------------
c
c_BEGIN_CCP4MAPHEAD
c
      subroutine ccp4maphead(name,nspgrp,cell,nu,nv,nw,nu1,nv1,nw1,nu2,
     +                       nv2,nw2)
c     ================================================================
c
c CCP4MAPhead - get limits of mask from file
c
c---- Function: read header information from a map file. It is used to
c     get the map limits before calling ccp4mapin.
c
c  Call: CALL ccp4maphead(name,nspgrp,cell,nu,nv,nw,nu1,nv1,nw1,nu2,
c       +                 nv2,nw2)
c
c---- Note that This differs from MRDHDS (which it calls) in that the
c     file is not left open but is closed (by a call to MRCLOS) before
c     the subroutine terminates. Note also that the map limits are
c     returned in x,y,z order rather than in fast, medium, slow order.
c
c
c---- Arguments:
c     =========
c
c  name    (I)  Logical file name (type character*8) e.g.'MAPIN'
c  nspgrp  (O)  Space group number (integer)
c  cell    (O)  6 word array for cell dimensions in Angstroms and
c               degrees (real)
c  nu      (O)  Sampling interval along whole cell on X (integer)
c  nv      (O)  Sampling interval along whole cell on Y (integer)
c  nw      (O)  Sampling interval along whole cell on Z (integer)
c  nu1     (O)  Start of map on X axis, in grid units (integer)
c  nv1     (O)  Start of map on Y axis, in grid units (integer)
c  nw1     (O)  Start of map on Z axis, in grid units (integer)
c  nu2     (O)  End of map on X axis (integer)
c  nv2     (O)  End of map on Y axis (integer)
c  nw2     (O)  End of map on Z axis (integer)
c
c_END_CCP4MAPHEAD
c
      implicit none
c
      character name*(*)
      integer nspgrp,nu,nv,nw,nu1,nv1,nw1,nu2,nv2,nw2
      real cell(6)
c
      integer jfms(3),juvw(3),mxyz(3),m1(3),m2(3),nsec,mode
      integer i,ifail,iprint
      real rmin,rmax,rmean,rrms
      character title*80
c
c read the file headers
      ifail=0
      iprint=1
      call mrdhds(7,name,title,nsec,jfms,mxyz,m1(3),m1(1),m2(1),
     +  m1(2),m2(2),cell,nspgrp,mode,rmin,rmax,rmean,rrms,ifail,iprint)
      call mrclos(7)
      m2(3)=m1(3)+nsec-1
      do 110 i=1,3
       juvw(jfms(i))=i
 110  continue
      nu1=m1(juvw(1))
      nu2=m2(juvw(1))
      nv1=m1(juvw(2))
      nv2=m2(juvw(2))
      nw1=m1(juvw(3))
      nw2=m2(juvw(3))
      nu=mxyz(1)
      nv=mxyz(2)
      nw=mxyz(3)
c
      return
      end
c
c ----------------------------------------------------------------------
c
c_BEGIN_CCP4MAPIN
c
      subroutine ccp4mapin(name,title,map,nu1,nv1,nw1,nu2,nv2,nw2)
c     ============================================================
c
c CCP4MAPIN - read a ccp4 logical map and store in xyz order
c
c---- Function: read whole map into an array and store in x,y,z order.
c     The map limits are required as input to dimension the array holding
c     the map, and can be obtained with a call to the subroutine
c     ccp4maphead.
c
c  Call: CALL ccp4mapin (name,title,map,nu1,nv1,nw1,nu2,nv2,nw2)
c
c---- ccp4mapin is a "wrapper" subroutine which utilises calls to the
c     following maplib routines: MRDHDS, MGULPR, MRCLOS.
c
c
c---- Arguments:
c     ==========
c
c  name    (O)  logical file name (type character) e.g. 'MAPIN'
c  title   (O)  Map title (type character)
c  map     (O)  Real array of dimension (nu1:nu2,nv1:nv2,nw1:nw2)
c               which stores the map which is read in
c  nu1     (I)  Start of map on X axis, in grid units (integer)
c  nv1     (I)  Start of map on Y axis, in grid units (integer)
c  nw1     (I)  Start of map on Z axis, in grid units (integer)
c  nu2     (I)  End of map on X axis (integer)
c  nv2     (I)  End of map on Y axis (integer)
c  nw2     (I)  End of map on Z axis (integer)
c
c_END_CCP4MAPIN
c
      implicit none
c
      integer nu1,nv1,nw1,nu2,nv2,nw2
      real map(nu1:nu2,nv1:nv2,nw1:nw2)
      character name*(*),title*(*)
c
      real lsec(0:100000)
c
      integer ifast,imedm,islow,lfast,lmedm,lslow,ierr,ifail,iprint
      integer m1(3),m2(3),ifms(3),jfms(3),juvw(3),mxyz(3)
      integer i,iu,iv,iw
      integer nsec,mode,mspgrp
      real cell(6),rmin,rmax,rmean,rrms
c
c
c Note: juvw convert from fast/med/slow to u/v/w
c       jfms convert from u/v/w to fast/med/slow
c
c now open map header and read map, re-ordering it as necessary
c
      ifail=0
      iprint=0
      call mrdhds(7,name,title,nsec,jfms,mxyz,m1(3),m1(1),m2(1),
     +  m1(2),m2(2),cell,mspgrp,mode,rmin,rmax,rmean,rrms,ifail,iprint)
      m2(3)=m1(3)+nsec-1
      do 110 i=1,3
       juvw(jfms(i))=i
 110  continue
c check we got the right header info:
      if (nu1.ne.m1(juvw(1)).or.nu2.ne.m2(juvw(1)).or.
     +    nv1.ne.m1(juvw(2)).or.nv2.ne.m2(juvw(2)).or.
     +    nw1.ne.m1(juvw(3)).or.nw2.ne.m2(juvw(3)))
     +  call ccperr(1,'ccp4mapin - mask grid sizes are inconsistent')
c find out the map grid dimensions
      lfast=m2(1)-m1(1)+1
      lmedm=m2(2)-m1(2)+1
      lslow=m2(3)-m1(3)+1
c
      if (lfast*lmedm.gt.100000)
     +  call ccperr(1,' ccp4mapin - Map section > lsec: recompile')
c now read the map in the order it is on file
      do 200 islow=0,lslow-1
c get a section
       ierr=0
       call mgulpr(7,lsec,ierr)
       if (ierr.ne.0) call ccperr(1,' ccp4mapin - ccp4 read error')
c now sort into uvw map
       ifms(3)=islow+m1(3)
       do 190 imedm=0,lmedm-1
       ifms(2)=imedm+m1(2)
       do 190 ifast=0,lfast-1
       ifms(1)=ifast+m1(1)
        iu=ifms(juvw(1))
        iv=ifms(juvw(2))
        iw=ifms(juvw(3))
        map(iu,iv,iw)=lsec(ifast+lfast*imedm)
 190   continue
 200  continue
c close the map file
      call mrclos(7)
c
      write (*,910)
 910  format (/' MAP/MASK READ SUCCESSFUL'//)
c
      return
      end
c
c ----------------------------------------------------------------------
c
c_BEGIN_CCP4MAPOUT
c
      subroutine ccp4mapout(name,title,map,nspgrp,cell,nu,nv,nw,nu1,
     +                     nv1,nw1,nu2,nv2,nw2)
c     ==============================================================
c
c CCP4MAPOUT - write a ccp4 logical map in xyz order
c
c---- Function: Write out a whole map in x,y,z order.
c
c  Call: CALL ccpmap4out(name,title,map,nspgrp,cell,nu,nv,nw,nu1,
c       +               nv1,nw1,nu2,nv2,nw2)
c
c---- ccpmap4out is a "wrapper" routine which utilises calls to the
c     following maplib subroutines: MWRHDL, MSYWRT, MSPEW, MWCLOSE. There
c     is also a call to the symlib routine MSYMLB.
c
c
c---- Arguments:
c     ==========
c
c  name    (I)  Logical file name (type character) e.g.'MAPIN'
c  title   (I)  Map title (type character)
c  map     (I)  Real array of dimension (nu1:nu2,nv1:nv2,nw1:nw2)
c               which stores the map being written out
c  nspgrp  (I)  Space group number (integer)
c  cell    (I)  6 word array for cell dimensions in Angstroms and degrees (real)
c  nu      (I)  Sampling interval along whole cell on X (integer)
c  nv      (I)  Sampling interval along whole cell on Y (integer)
c  nw      (I)  Sampling interval along whole cell on Z (integer)
c  nu1     (I)  Start of map on X axis, in grid units (integer)
c  nv1     (I)  Start of map on Y axis, in grid units (integer)
c  nw1     (I)  Start of map on Z axis, in grid units (integer)
c  nu2     (I)  End of map on X axis (integer)
c  nv2     (I)  End of map on Y axis (integer)
c  nw2     (I)  End of map on Z axis (integer)
c
c_END_CCP4MAPOUT
c
      implicit none
c
      integer nspgrp,nu,nv,nw,nu1,nv1,nw1,nu2,nv2,nw2
      real cell(6),map(nu1:nu2,nv1:nv2,nw1:nw2)
      character name*(*),title*(*)
c
      real lsec(0:100000)
c
      integer ifast,imedm,islow,lfast,lmedm,lslow
      integer m1(3),m2(3),ifms(3),jfms(3),juvw(3),mxyz(3)
      integer i,iu,iv,iw
      integer nsec,mode,nsym,nsymp
      real rsym(4,4,192)
      character*10 namspg,nampg
c
c
c spacegroups with yxz=1 and zxy=2 axis ordering
      integer axis(230)
      data axis/2,2,2,2,1,1,1,1,1,2,1,1,1,1,1,2,2,2,1,2,2,1,2,207*1/
c
c Note: juvw convert from fast/med/slow to u/v/w
c       jfms convert from u/v/w to fast/med/slow
c
c now open map header and read map, re-ordering it as necessary
c
      if (axis(mod(nspgrp,1000)).eq.1) then
       jfms(1)=2
       jfms(2)=1
       jfms(3)=3
      else
       jfms(1)=3
       jfms(2)=1
       jfms(3)=2
      endif
c
      do 110 i=1,3
       juvw(jfms(i))=i
 110  continue
c
      mxyz(1)=nu
      mxyz(2)=nv
      mxyz(3)=nw
      m1(juvw(1))=nu1
      m2(juvw(1))=nu2
      m1(juvw(2))=nv1
      m2(juvw(2))=nv2
      m1(juvw(3))=nw1
      m2(juvw(3))=nw2
      nsec=m2(3)-m1(3)+1
      mode=2
c
      call msymlb(7,nspgrp,namspg,nampg,nsym,nsymp,rsym)
c
      call mwrhdl(7,name,title,nsec,jfms,mxyz,
     +  m1(3),m1(1),m2(1),m1(2),m2(2),cell,nspgrp,mode)
      call msywrt(7,nsym,rsym)
c
c find out the mask grid dimensions
      lfast=m2(1)-m1(1)+1
      lmedm=m2(2)-m1(2)+1
      lslow=m2(3)-m1(3)+1
c
      if (lfast*lmedm.gt.100000)
     +  call ccperr(1,' ccp4mapout - Mask section > lsec: recompile')
c now write the map in the new order
      do 200 islow=0,lslow-1
c sort onto section
       ifms(3)=islow+m1(3)
       do 190 imedm=0,lmedm-1
       ifms(2)=imedm+m1(2)
       do 190 ifast=0,lfast-1
       ifms(1)=ifast+m1(1)
        iu=ifms(juvw(1))
        iv=ifms(juvw(2))
        iw=ifms(juvw(3))
        lsec(ifast+lfast*imedm)=map(iu,iv,iw)
 190   continue
c write the section
       call mspew(7,lsec)
 200  continue
c close the map file
      call mwclose(7)
c
      return
      end
c
c ----------------------------------------------------------------------
c