C   Adapted loosely from f2mtz.f
C   5/3/2004: Updated to show writing of crystals and datasets

C   This program writes a basic MTZ file, using information
C   in the DATA statements. Obviously, you will want to replace
C   this with your own data.

C   It should be run as "mtz_write HKLOUT foo.mtz"

      PROGRAM WRITEMTZ

C Max number of labels
      integer maxlab
      parameter (maxlab = 200)
C Max number of symmetry operations
      integer MAXSYM
      parameter (MAXSYM = 192)

      INTEGER I,J,MTZOUT,NLPRGO,NSYM,NSYMP,NSPGR,ISORT(5),NREFW
      REAL WAVELENGTH,CELL(6),RSYM(4,4,MAXSYM),ADATA(MAXLAB),
     +   COLUMNDATA1(MAXLAB),COLUMNDATA2(MAXLAB)
      CHARACTER TITLE*80,CRYSTALNAME*20,PROJECTNAME*20,
     +   DATASETNAME*20,
     +   LSPRGO(MAXLAB)*30,CTPRGO(MAXLAB)*1,
     +   XNAME(MAXLAB)*20,PNAME(MAXLAB)*20,DNAME(MAXLAB)*20,
     +   LTYPE*1,SPGRNAM*10,PGNAM*10

C Example data from distributed toxd example
      DATA TITLE /'From mtz_write'/
      DATA CELL /73.582,38.733,23.189,90.000,90.000,90.000/
      DATA WAVELENGTH /0.0/
      DATA PROJECTNAME /'TOXD'/
      DATA CRYSTALNAME /'NATIVE'/
      DATA DATASETNAME /'NATIVE'/
      DATA NLPRGO /5/
      DATA LSPRGO /'H','K','L','F','SIGF',195*' '/
      DATA CTPRGO /'H','H','H','F','Q',195*' '/
      DATA NSPGR /19/
      DATA LTYPE,SPGRNAM /'P','P212121'/
C Sort on H K L
      DATA ISORT /1,2,3,0,0/
C Reflection data
      DATA COLUMNDATA1 /0.0,0.0,2.0,626.00,112.00,195*0.0/
      DATA COLUMNDATA2 /0.0,0.0,4.0,2853.00,59.00,195*0.0/

C CCP4 initialisations
      call ccpfyp
C MTZ-specific initialisations
      call mtzini

C Open mtz file
      MTZOUT = 1
      write(6,*)'pjx opening file'
      call lwopen(MTZOUT,'HKLOUT')
      write(6,*)'pjx back from lwopen'

C Write title to MTZ header
      write(6,*)'pjx lwtitl'
      call lwtitl (MTZOUT,TITLE,0)
      write(6,*)'pjx back from lwtitl'

C Write cell to MTZ header
      write(6,*)'pjx lwcell'
      call lwcell(MTZOUT,CELL)
      write(6,*)'pjx back from lwcell'

C *** For CCP4 5.0 libraries
C Store the crystal, project and dataset names plus
C associated data in the mtz header:
      write(6,*)'pjx lwidx'
      call lwidx(MTZOUT,PROJECTNAME,CRYSTALNAME,DATASETNAME,
     +     CELL,WAVELENGTH)
      write(6,*)'pjx back from lwidx'
C Alternatively:
C *** For CCP4 4.2 (and earlier)
C Store the project and dataset names
C      call lwid(MTZOUT,PROJECTNAME,DATASETNAME)

C Store the column labels in the mtz header:
C Simpler alternative to LKYOUT/LWASSN
      write(6,*)'pjx lwclab'
      call lwclab(MTZOUT,LSPRGO,NLPRGO,CTPRGO,0)
      write(6,*)'pjx back from lwclab'

C *** For CCP4 5.0 libraries
C Associate columns with crystals and datasets
      do i=1,NLPRGO
        xname(i) = CRYSTALNAME
        dname(i) = DATASETNAME
      enddo
      write(6,*)'pjx lwidasx'
      call lwidasx(MTZOUT,NLPRGO,XNAME,DNAME,0)
      write(6,*)'pjx back from lwidasx'
C Alternatively:
C *** For CCP4 4.2 (and earlier)
C Associate columns with projects and datasets
C      do i=1,NLPRGO
C        pname(i) = PROJECTNAME
C        dname(i) = DATASETNAME
C      enddo
C      call lwidas(MTZOUT,NLPRGO,PNAME,DNAME,0)

C Get symmetry info from file SYMOP (default set by ccp4.setup)
C on channel 24
      CALL MSYMLB(24,NSPGR,SPGRNAM,PGNAM,NSYMP,NSYM,RSYM)
C Store the symmetry in the mtz header:
      call lwsymm (MTZOUT,NSYM,NSYMP,RSYM,LTYPE,NSPGR,SPGRNAM,PGNAM)

C Store program details in history header:
      call lwhstl (MTZOUT, ' ')

C  Store sort order to be used
      CALL LWSORT(MTZOUT,ISORT)

C   Loop over ouput reflections
C   Data for 2 reflections given here - normally calculated and
C   lots more!
      NREFW = 2
      DO I = 1,NREFW

C  Initialise ADATA with Missing Number Flags
        CALL EQUAL_MAGIC(MTZOUT,ADATA,NLPRGO)

C  Fill in data columns where possible (all in this simple e.g.)
        do 31 J = 1,NLPRGO
          if (i.eq.1) adata(J) = columndata1(J)
          if (i.eq.2) adata(J) = columndata2(J)
 31     continue

C  Write out reflection
        call lwrefl (MTZOUT, adata)

      ENDDO

C Close mtz file, print full header:
      call lwclos (MTZOUT, 3)

      call ccperr(0,'Normal termination')
      end