##
## Read an X-PLOR, ASCII format Electron Density Map file 
## into a volume set on the 'top' molecule in VMD. 
##
## $Id: readedm.tcl,v 1.5 2002/02/22 20:53:07 johns Exp johns $
## 
## Requires VMD 1.7.1 or newer
## 
## Authors: Deyu Lu and John Stone
## Theoretical Biophysics Group
## University of Illinois at Urbana-Champaign
##
## Example:
##   source readedm.tcl
##   mol load pdb mol_final.pdb
##     readedm mol_2fofc.map
##

## Remove leading spaces, and replace subsequent 
##   duplicate spaces with a single space, for use with "split".
proc eatspaces { string } {
  set newstr "";
  set newstr2 "";
  regsub -all -- {[\s]+} $string " " newstr
  regsub -- {[\s]+} $newstr "" newstr2
  return $newstr2
}

##
## Read an X-PLOR, ASCII format Electron Density Map file 
## into a volume set in VMD.
##
proc readedm { filename } {
  variable edmfile;

  puts "EDM Reader, opening $filename"
  set edmfile [open $filename r]

  ##
  ## Read the header information first 
  ##
  gets $edmfile;   # read first line, and ignore
  ## puts [lindex [split [eatspaces [gets $edmfile]]] 0]
  set ntitle [lindex [split [eatspaces [gets $edmfile]]] 0]
  ## puts "Number of title lines: $ntitle"

  ## read in title and comment lines
  for {set i 0} {$i < $ntitle} {incr i} {
    puts "   [gets $edmfile]" 
  }
  
  ## read in box dimensions and grid spacing deltas
  set diminfo [concat [split [eatspaces [gets $edmfile]]]] 
  set na    [lindex $diminfo 0]
  set amin  [lindex $diminfo 1]
  set amax  [lindex $diminfo 2]
  set nb    [lindex $diminfo 3]
  set bmin  [lindex $diminfo 4]
  set bmax  [lindex $diminfo 5]
  set nc    [lindex $diminfo 6]
  set cmin  [lindex $diminfo 7]
  set cmax  [lindex $diminfo 8]

  ## find number of samples in each dimension
  set xsize [expr $amax - $amin + 1]
  set ysize [expr $bmax - $bmin + 1]
  set zsize [expr $cmax - $cmin + 1]

  ## read in 6 values related to the box orientation
  set diminfo [concat [split [eatspaces [gets $edmfile]]]] 
  set a     [lindex $diminfo 0] 
  set b     [lindex $diminfo 1] 
  set c     [lindex $diminfo 2] 
  set alpha [lindex $diminfo 3] 
  set beta  [lindex $diminfo 4] 
  set gamma [lindex $diminfo 5] 

  ## find box coordinates
  set xdelta [expr $a / ($na * 1.0)]
  set ydelta [expr $b / ($nb * 1.0)]
  set zdelta [expr $c / ($nc * 1.0)]

  set xmin [expr $amin * $xdelta]
  set xmax [expr $amax * $xdelta]
  set ymin [expr $bmin * $ydelta]
  set ymax [expr $bmax * $ydelta]
  set zmin [expr $cmin * $zdelta]
  set zmax [expr $cmax * $zdelta]

  puts " Grid dimensions:  X: $xsize, Y: $ysize, Z: $zsize"
  puts "   Grid elements: [expr $xsize * $ysize * $zsize]" 
  puts "Grid orientation:      a: $a,    b: $b,     c: $c," 
  puts "                   alpha: $alpha, beta: $beta, gamma: $gamma" 

  ## set origin
  set xOrigin $xmin
  set yOrigin $ymin
  set zOrigin $zmin

  puts " xOrigin, yOrigin, zOrigin: $xOrigin, $yOrigin, $zOrigin"

  ##
  ## set angles to arc
  ##
  set pi [expr acos(-1)]
  set alpha1 [expr $alpha*$pi/180.0]
  set beta1 [expr $beta*$pi/180.0]
  set gamma1 [expr $gamma*$pi/180.0]

  ##
  ## Handle non-orthogonal unit cells
  ## set xVec yVec zVec cell side vectors
  ##
  set xVec [list [expr $xsize*$xdelta] 0 0 ]
  set yVec [list [expr cos($gamma1)*$xsize*$xdelta ] [expr sin($gamma1)*$ysize*$ydelta] 0]
  set z1 [expr cos($beta1)]
  set z2 [expr (cos($alpha1)-cos($beta1)*cos($gamma1))/sin($gamma1)]
  set z3 [expr sqrt(1.0-pow($z1,2)-pow($z2,2))]
  set zVec [list [expr $z1*$xsize*$xdelta] [expr $z2*$ysize*$ydelta] [expr $z3*$zsize*$zdelta] ]

  puts "xVec: $xVec"
  puts "yVec: $yVec"
  puts "zVec: $zVec" 

  ##   
  ## read in ZYX data ordering string 
  ##   
  puts " [gets $edmfile] "

  ##    
  ##  find the # of samples in 1 slide of the Z axis
  ##
  set n_gslide [expr $xsize*$ysize]
  puts "# of grids in 1 slide: $n_gslide"

  ## set the maximum number of data items per line in the file,
  ## and prepare to read in the actual density data values.
  set n_gline  6
  puts "# of grids in 1 line:  $n_gline"

  set n_line [expr floor($n_gslide/$n_gline)]
  if {[expr $n_line*$n_gline]<$n_gslide} {
    set n_line [expr $n_line+1]
  }
  puts "# of lines in 1 slide: $n_line"
 
  set zmap {};  # clear the density map to be emtpy. 

  ## read in all of the density values  
  for {set z 0} {$z<$zsize} {incr z} {
    gets $edmfile;  # read in the Z-plane index (and throw away)
    for {set line_id 1} {$line_id<$n_line } {incr line_id} {
      binary scan [gets $edmfile] a12a12a12a12a12a12 v1 v2 v3 v4 v5 v6
      lappend zmap $v1 $v2 $v3 $v4 $v5 $v6
    }
    # read the last line; binary scan returns how many items were read
    set n [binary scan [gets $edmfile] a12a12a12a12a12a12 v(1) v(2) v(3) v(4) v(5) v(6)]
    for { set i 1 } { $i <= $n } { incr i } {
      lappend zmap $v($i)
    }
  }
   
  puts [llength $zmap]
  puts [eatspaces [gets $edmfile]] 
  ## -9999 

  ## close the EDM file
  close $edmfile

  ## Add the complete volume set into the top molecule of VMD
  mol volume top "Electron Density" [list $xOrigin $yOrigin $zOrigin] $xVec $yVec $zVec $xsize $ysize $zsize $zmap

  ## At this point, VMD now has the data, and its up to the user to activate
  ## representations on it. 
}