#include "copyright.h"

! Routines for writing volumetric data to file according to the DX file format.
! A DX file format description can be found, for example, here:
! http://www.poissonboltzmann.org/file-formats/mesh-and-data-formats/
! opendx-scalar-data

! `voldata` is assumed to be a grid of values with shape (1:xm, 1:ym, 1:zm).
! `xm, ym, zm` define the number of nodes along x, y, and z direction.
! `gox, goy, goz` define the grid origin in x, y, and z directions.
! `h` defines the grid spacing in all dimensions.

! Originally written by .... whom?
! Consolidated by Jan-Philip Gehrcke, TU Dresden, Germany (2012)


#define _REAL_ double precision


subroutine gen_dx_file(&
        xm, ym, zm, h, gox, goy, goz, voldata, dxfilename, dxfileid, dataname)
    integer, intent(in) :: xm, ym, zm, dxfileid
    _REAL_, intent(in) ::  h, gox, goy, goz
    _REAL_, intent(in) ::  voldata(1:xm, 1:ym, 1:zm)
    character(len=*), intent(in) :: dxfilename
    character(len=*), intent(in) :: dataname
    integer :: cnt, i, j, k

    ! Unit number must be unique within the program.
    ! Trailing spaces in a file name can be annoying. Trim them.
    open(unit=dxfileid, file=trim(dxfilename))
    write(dxfileid, 100) xm, ym, zm
    write(dxfileid, 110) gox+h, goy+h, goz+h
    write(dxfileid, 120) sngl(h), 0, 0
    write(dxfileid, 121) 0, sngl(h), 0
    write(dxfileid, 122) 0, 0, sngl(h)
    write(dxfileid, 130) xm, ym, zm
    write(dxfileid, 140) xm*ym*zm

    cnt = 0
    do i = 1, xm
        do j = 1, ym
            do k = 1, zm
                cnt = cnt+1
                write(dxfileid, 150, advance='no') voldata(i,j,k)
                if (MOD(cnt, 3) == 0) write(dxfileid, *) ''
            end do
        end do
    end do
    if (.NOT. MOD(cnt, 3) == 0) write(dxfileid, *) ''

    write(dxfileid, *) 'attribute "dep" string "positions"'
    write(dxfileid, *) 'object "', dataname, '" class field'
    write(dxfileid, *) 'component "positions" value 1'
    write(dxfileid, *) 'component "connections" value 2'
    write(dxfileid, *) 'component "data" value 3'
    close(dxfileid)

    100 format('object 1 class gridpositions counts', 3I5)
    110 format('origin', 3F9.4)
    120 format('delta', F10.7, 2I2)
    121 format('delta', I2, F10.7, I2)
    122 format('delta', 2I2, F10.7)
    130 format('object 2 class gridconnections counts', 3I5)
    140 format('object 3 class array type double rank 0 items', I10,&
        ' data follows')
    150 format(ES19.10, ' ')
end subroutine gen_dx_file


! subroutine printphidx(phigrid, xm, ym, zm, gox, goy, goz, h)
!     use poisson_boltzmann, only : frcfac, eps0
!     !Passed Variables
!     integer xm, ym, zm
!     _REAL_ phigrid(1:xm, 1:ym, 1:zm)
!     _REAL_ gox, goy, goz, h
!     !Internal Variables
!     _REAL_ phitemp(xm, ym, zm)
!     character(len=6) :: phi_filename
!     character(len=7) :: phi_setname
!     integer :: unitid = 8888
!
!     phi_filename = "dxfile"
!     phi_setname = "dataset"
!
!     phitemp = phigrid / eps0 !*frcfac
!     call gen_dx_file(&
!         xm, ym, zm, h, &
!         gox, goy, goz, &
!         phitemp, unitid, phi_filename, phi_setname)
! end subroutine printphidx


! subroutine gen_integer_dx_file(&
!         xm,ym,zm,h,gox,goy,goz,voldata,filename,fn,dataname)
!     !This subroutine is not actually in use and has a newline bug.
!     implicit none
!     !passed variables
!     integer xm,ym,zm,fn
!     _REAL_  h, gox, goy, goz
!     integer  voldata(1:xm,1:ym,1:zm)
!     character*(*), intent(in) :: filename, dataname
!     !internal variables
!     integer cnt, i, j, k
!
!     !write(6,*) 'writing potential map in dx format'
!     open(fn,file=filename)
!     write(fn,200) xm,ym,zm
!     write(fn,210) gox+h,goy+h,goz+h
!     write(fn,220) sngl(h),0,0
!     write(fn,221) 0,sngl(h),0
!     write(fn,222) 0,0,sngl(h)
!     write(fn,230) xm,ym,zm
!     write(fn,240) xm*ym*zm
!     cnt=0
!     do i=1,xm; do j=1,ym;
!         do k=1,zm
!             cnt = cnt+1
!             write(fn,250,advance='no') voldata(i,j,k)
!             if (cnt >= 3) then
!                 cnt = 0
!                 write(fn,*)''
!             end if
!         end do
!     end do; end do
!     write(fn,*) 'attribute "dep" string "positions"'
!     write(fn,*) 'object "',dataname,'" class field'
!     write(fn,*) 'component "positions" value 1'
!     write(fn,*) 'component "connections" value 2'
!     write(fn,*) 'component "data" value 3'
!     close(fn)
!     200 format('object 1 class gridpositions counts',3I5)
!     210 format('origin',3F9.4)
!     220 format('delta',F10.7,2I2)
!     221 format('delta',I2,F10.7,I2)
!     222 format('delta',2I2,F10.7)
!     230 format('object 2 class gridconnections counts',3I5)
!     240 format('object 3 class array type double rank 0 items',I10,' data follows')
!     250 format(I12,' ')
! end subroutine gen_integer_dx_file