#!/usr/bin/tclsh
# by Jeff Comer jeffcomer at gmail dot com

if {[llength $argv] != 4} {
    error "Usage: $argv0 delTime histFile0 histFile1 outFile"
    exit
}
foreach {delTime histFile0 histFile1 outFile} $argv { break }

proc readHistData {fileName} {
    set in [open $fileName r]

    set entry -1
    set onData 0
    set data {}
    set retList {}
    while {[gets $in lin] >= 0} {
	set lin1 [string trim $lin]

	if {[string length $lin1] == 0} {
	    set onData 0
	} elseif {[string match "#*" $lin1]} {
	    if {$onData} {
		set onData 0
		# Store the current data and empty it.
		lappend retList $data
		set data {}
		incr entry
	    }
	    # If !$onData, this is just another comment. Ignore.
	} else {
	    if {![string is double [lindex $lin1 0]]} {
		error "ERROR! File `$fileName' has bad line `$lin'"
	    }

	    # This is data.
	    set onData 1
	    # Add this line to the list.
	    lappend data [concat $lin1]
	}
    }

    close $in
    return $retList
}

proc rmsd {l0 l1} {
    if {[llength $l0] != [llength $l1]} {
	error "ERROR! Lists must have the same length to calculate the RMSD. [llength $l0] != [llength $l1]"
    }

    set sum 0.0
    foreach item0 $l0 item1 $l1 {
	set dev [expr {[lindex $item0 1]-[lindex $item1 1]}]
	set sum [expr {$sum + $dev*$dev}]
    }

    return [expr {sqrt($sum/[llength $l0])}]
}

# Read the data from the history file.
set hist0 [readHistData $histFile0]
set dataN0 [llength $hist0]

puts "File `$histFile0' has $dataN0 frames."
set hist1 [readHistData $histFile1]
set dataN1 [llength $hist1]
puts "File `$histFile1' has $dataN1 frames."

# Choose the lowest number of frames.
set n [expr {($dataN0<=$dataN1)?$dataN0:$dataN1}]

# Write the output.
puts "Output will consider $n frames"
set out [open $outFile w]
for {set i 0} {$i < $n} {incr i} {
    set t [expr {$delTime*($i+1)}]
    set rmsd [rmsd [lindex $hist0 $i] [lindex $hist1 $i]]
    puts $out "$t $rmsd"
}
close $out
