# Draw displays for simulation time and pulling force 
#
# This file contains 3 routines (followed by other junk):
# disp_time (140ps display}
# disp_time3 (30ps display)
# disp_force (3000pN display}
# 
# Works best in orthographic projection mode!

proc disp_time {time} {

  # radius of display in Angstrom 
  set radd 6.0
  # radius of tube in Angstrom 
  set radt 0.5
  # resolution of tube, number of polygon edges   
  set poly 100
  # offset of display from center in Angstrom
  set xsep 0.0 
  set ysep 10.0 

  # initialize.
  set twopi [expr 4*asin(1.0)]
  set pi [expr 2*asin(1.0)]
  set delta [expr $twopi / $poly]

  # draw rim and face color
  draw color silver
  set center "$xsep $ysep 0.0"
  set angle1 [expr $twopi * 10.0 / 140.0]
  set sx $xsep
  set sy [expr $radd + $ysep]
  set old "$sx $sy 0.0"
  set angle $delta
  set i 0
  while {$angle <= $angle1} {
	set sx [expr $radd * sin($angle) + $xsep]
	set sy [expr $radd * cos($angle) + $ysep]
        set new "$sx $sy 0.0"
        draw color silver
 	draw cylinder $old $new radius $radt resolution 10
	set old $new
	set angle [expr $angle + $delta]
	incr i
  }	
  while {$i <= $poly} {
	set sx [expr $radd * sin($angle) + $xsep]
	set sy [expr $radd * cos($angle) + $ysep]
        set new "$sx $sy 0.0"
        draw color silver
	draw cylinder $old $new radius $radt resolution 10
        draw color blue
	draw materials off
	draw triangle $center $old $new
        draw materials on
	set old $new
	set angle [expr $angle + $delta]
	incr i
  }	
  set sx $xsep
  set sy [expr $radd + $ysep]
  set old "$sx $sy 0.0"
  set angle $delta
  set i 0
  while {$angle <= $angle1} {
	set sx [expr $radd * sin($angle) + $xsep]
	set sy [expr $radd * cos($angle) + $ysep]
        set new "$sx $sy 0.0"
        draw color red
	draw materials off
	draw triangle $center $old $new
        draw materials on
	set old $new
	set angle [expr $angle + $delta]
	incr i
  }	

  # draw hand
  draw color white
  set fsep [expr 2.0 * $radt]
  set angle [expr ($time / 140.0) * $twopi]
  set sx [expr 1.2 * $radd * sin($angle) + $xsep]
  set sy [expr 1.2 * $radd * cos($angle) + $ysep]
  set tx [expr -0.2 * $radd * sin($angle) + $xsep]
  set ty [expr -0.2 * $radd * cos($angle) + $ysep]
  draw cone "$tx $ty $fsep" "$sx $sy $fsep" radius $radt resolution 10
  set radk [expr 1.5 * $radt]
  draw sphere "$xsep $ysep $fsep" radius $radk resolution 10

  # draw face 
  set angle 0.0
  set angcount 0
  while {$angcount <= 70} {
     set sx [expr 1.1 * $radd * sin($angle) + $xsep]
     set sy [expr 1.1 * $radd * cos($angle) + $ysep]
     set tx [expr 0.7 * $radd * sin($angle) + $xsep]
     set ty [expr 0.7 * $radd * cos($angle) + $ysep]
     set radf [expr 0.6 * $radt]
     draw color silver
     draw cylinder "$tx $ty 0.0" "$sx $sy 0.0" radius $radf resolution 10 filled yes
     set ux [expr 1.3 * $radd * sin($angle) + $xsep]
     set uy [expr 1.3 * $radd * cos($angle) + $ysep]
     draw color white
     draw text "$ux $uy 0.0" "$angcount"
     set angle [expr $angle + $twopi * 10.0 / 140.0]
     set angcount [expr $angcount + 10]
  }
  while {$angle <= $twopi} {
     set sx [expr 1.1 * $radd * sin($angle) + $xsep]
     set sy [expr 1.1 * $radd * cos($angle) + $ysep]
     set tx [expr 0.7 * $radd * sin($angle) + $xsep]
     set ty [expr 0.7 * $radd * cos($angle) + $ysep]
     set radf [expr 0.6 * $radt]
     draw color silver
     draw cylinder "$tx $ty 0.0" "$sx $sy 0.0" radius $radf resolution 10 filled yes
     set ux [expr 1.3 * $radd * sin($angle) + $xsep - 1.5]
     set uy [expr 1.3 * $radd * cos($angle) + $ysep]
     draw color white
     draw text "$ux $uy 0.0" "$angcount"
     set angle [expr $angle + $twopi * 10.0 / 140.0]
     set angcount [expr $angcount + 10]
  }

  # draw caption
   set sx [expr $xsep - 2]
   set sy [expr -1.5 * $radd + $ysep]
   draw color white
   draw text "$sx $sy 0.0" "time (ps)"
}


proc disp_time3 {time} {

  # radius of display in Angstrom 
  set radd 6.0
  # radius of tube in Angstrom 
  set radt 0.5
  # resolution of tube, number of polygon edges   
  set poly 100
  # offset of display from center in Angstrom
  set xsep 0.0 
  set ysep 10.0 

  # initialize.
  set twopi [expr 4*asin(1.0)]
  set pi [expr 2*asin(1.0)]
  set delta [expr $twopi / $poly]

  # draw rim and face color
  draw color silver
  set center "$xsep $ysep 0.0"
  set angle1 [expr $twopi * 10.0 / 30.0]
  set sx $xsep
  set sy [expr $radd + $ysep]
  set old "$sx $sy 0.0"
  set angle $delta
  set i 0
  while {$angle <= $angle1} {
	set sx [expr $radd * sin($angle) + $xsep]
	set sy [expr $radd * cos($angle) + $ysep]
        set new "$sx $sy 0.0"
        draw color silver
 	draw cylinder $old $new radius $radt resolution 10
	set old $new
	set angle [expr $angle + $delta]
	incr i
  }	
  while {$i <= $poly} {
	set sx [expr $radd * sin($angle) + $xsep]
	set sy [expr $radd * cos($angle) + $ysep]
        set new "$sx $sy 0.0"
        draw color silver
	draw cylinder $old $new radius $radt resolution 10
        draw color blue
	draw materials off
	draw triangle $center $old $new
        draw materials on
	set old $new
	set angle [expr $angle + $delta]
	incr i
  }	
  set sx $xsep
  set sy [expr $radd + $ysep]
  set old "$sx $sy 0.0"
  set angle $delta
  set i 0
  while {$angle <= $angle1} {
	set sx [expr $radd * sin($angle) + $xsep]
	set sy [expr $radd * cos($angle) + $ysep]
        set new "$sx $sy 0.0"
        draw color red
	draw materials off
	draw triangle $center $old $new
        draw materials on
	set old $new
	set angle [expr $angle + $delta]
	incr i
  }	

  # draw hand
  draw color white
  set fsep [expr 2.0 * $radt]
  set angle [expr ($time / 30.0) * $twopi]
  set sx [expr 1.2 * $radd * sin($angle) + $xsep]
  set sy [expr 1.2 * $radd * cos($angle) + $ysep]
  set tx [expr -0.2 * $radd * sin($angle) + $xsep]
  set ty [expr -0.2 * $radd * cos($angle) + $ysep]
  draw cone "$tx $ty $fsep" "$sx $sy $fsep" radius $radt resolution 10
  set radk [expr 1.5 * $radt]
  draw sphere "$xsep $ysep $fsep" radius $radk resolution 10

  # draw face 
  set angle 0.0
  set angcount 0
  while {$angcount <= 15} {
     set sx [expr 1.1 * $radd * sin($angle) + $xsep]
     set sy [expr 1.1 * $radd * cos($angle) + $ysep]
     set tx [expr 0.7 * $radd * sin($angle) + $xsep]
     set ty [expr 0.7 * $radd * cos($angle) + $ysep]
     set radf [expr 0.6 * $radt]
     draw color silver
     draw cylinder "$tx $ty 0.0" "$sx $sy 0.0" radius $radf resolution 10 filled yes
     set ux [expr 1.3 * $radd * sin($angle) + $xsep]
     set uy [expr 1.3 * $radd * cos($angle) + $ysep]
     draw color white
     draw text "$ux $uy 0.0" "$angcount"
     set angle [expr $angle + $twopi * 5.0 / 30.0]
     set angcount [expr $angcount + 5]
  }
  while {$angle <= $twopi} {
     set sx [expr 1.1 * $radd * sin($angle) + $xsep]
     set sy [expr 1.1 * $radd * cos($angle) + $ysep]
     set tx [expr 0.7 * $radd * sin($angle) + $xsep]
     set ty [expr 0.7 * $radd * cos($angle) + $ysep]
     set radf [expr 0.6 * $radt]
     draw color silver
     draw cylinder "$tx $ty 0.0" "$sx $sy 0.0" radius $radf resolution 10 filled yes
     set ux [expr 1.3 * $radd * sin($angle) + $xsep - 1.5]
     set uy [expr 1.3 * $radd * cos($angle) + $ysep]
     draw color white
     draw text "$ux $uy 0.0" "$angcount"
     set angle [expr $angle + $twopi * 5.0 / 30.0]
     set angcount [expr $angcount + 5]
  }

  # draw caption
   set sx [expr $xsep - 2]
   set sy [expr -1.5 * $radd + $ysep]
   draw color white
   draw text "$sx $sy 0.0" "time (ps)"
}


proc disp_force {force} {

  # radius of display in Angstrom 
  set radd 6.0
  # radius of tube in Angstrom 
  set radt 0.5
  # resolution of tube, number of polygon edges   
  set poly 100
  # offset of display from center in Angstrom
  set xsep 0.0 
  set ysep -10.0 

  # initialize
  set twopi [expr 4*asin(1.0)]
  set pi [expr 2*asin(1.0)]
  set delta [expr $twopi / $poly]

  # draw rim and face color
  draw color silver
  set center "$xsep $ysep 0.0"
  set angle1 [expr $twopi * 50.0 / 360.0]
  set angle2 [expr $twopi * 150.0 / 360.0]
  set angle3 [expr $twopi * 210.0 / 360.0]
  set angle4 [expr $twopi * 310.0 / 360.0]
  set sx $xsep
  set sy [expr $radd  + $ysep]
  set old "$sx $sy 0.0"
  set angle $delta
  set i 0
  while {$angle <= $angle1} {
	set sx [expr $radd * sin($angle) + $xsep]
	set sy [expr $radd * cos($angle) + $ysep]
        set new "$sx $sy 0.0"
        draw color silver
 	draw cylinder $old $new radius $radt resolution 10
        draw color yellow
	draw materials off
	draw triangle $center $old $new
        draw materials on
	set old $new
	set angle [expr $angle + $delta]
	incr i
  }	
  while {$angle <= $angle2} {
	set sx [expr $radd * sin($angle) + $xsep]
	set sy [expr $radd * cos($angle) + $ysep]
        set new "$sx $sy 0.0"
        draw color silver
 	draw cylinder $old $new radius $radt resolution 10
        draw color red
	draw materials off
	draw triangle $center $old $new
        draw materials on
	set old $new
	set angle [expr $angle + $delta]
	incr i
  }	
  while {$angle <= $angle3} {
	set sx [expr $radd * sin($angle) + $xsep]
	set sy [expr $radd * cos($angle) + $ysep]
        set new "$sx $sy 0.0"
        draw color silver
 	draw cylinder $old $new radius $radt resolution 10
        draw color blue
	draw materials off
	draw triangle $center $old $new
        draw materials on
	set old $new
	set angle [expr $angle + $delta]
	incr i
  }	
  while {$angle <= $angle4} {
	set sx [expr $radd * sin($angle) + $xsep]
	set sy [expr $radd * cos($angle) + $ysep]
        set new "$sx $sy 0.0"
        draw color silver
 	draw cylinder $old $new radius $radt resolution 10
        draw color green
	draw materials off
	draw triangle $center $old $new
        draw materials on
	set old $new
	set angle [expr $angle + $delta]
	incr i
  }	
  while {$i <= $poly} {
	set sx [expr $radd * sin($angle) + $xsep]
	set sy [expr $radd * cos($angle) + $ysep]
        set new "$sx $sy 0.0"
        draw color silver
	draw cylinder $old $new radius $radt resolution 10
        draw color yellow
	draw materials off
	draw triangle $center $old $new
        draw materials on
	set old $new
	set angle [expr $angle + $delta]
	incr i
  }	

  # draw hand
  draw color white
  set fsep [expr 2.0 * $radt]
  set angle [expr $angle3 + ($force / 3000.0) * $twopi * 300.0 / 360.0]
  set sx [expr 1.2 * $radd * sin($angle) + $xsep]
  set sy [expr 1.2 * $radd * cos($angle) + $ysep]
  set tx [expr -0.2 * $radd * sin($angle) + $xsep]
  set ty [expr -0.2 * $radd * cos($angle) + $ysep]
  draw cone "$tx $ty $fsep" "$sx $sy $fsep" radius $radt resolution 10
  set radk [expr 1.5 * $radt]
  draw sphere "$xsep $ysep $fsep" radius $radk resolution 10

  # draw face 
  set angle 0.0
  set angcount 1500
  while {$angcount <= 3000} {
     set sx [expr 1.1 * $radd * sin($angle) + $xsep]
     set sy [expr 1.1 * $radd * cos($angle) + $ysep]
     set tx [expr 0.7 * $radd * sin($angle) + $xsep]
     set ty [expr 0.7 * $radd * cos($angle) + $ysep]
     set radf [expr 0.6 * $radt]
     draw color silver
     draw cylinder "$tx $ty 0.0" "$sx $sy 0.0" radius $radf resolution 10 filled yes
     set ux [expr 1.3 * $radd * sin($angle) + $xsep]
     set uy [expr 1.3 * $radd * cos($angle) + $ysep]
     draw color white
     draw text "$ux $uy 0.0" "$angcount"
     set angle [expr $angle + $twopi * 25.0 / 360.0]
     set angcount [expr $angcount + 250]
  }
  set angle [expr $twopi * 210.0 / 360.0]
  set angcount 0
  while {$angcount < 1500} {
     set sx [expr 1.1 * $radd * sin($angle) + $xsep]
     set sy [expr 1.1 * $radd * cos($angle) + $ysep]
     set tx [expr 0.7 * $radd * sin($angle) + $xsep]
     set ty [expr 0.7 * $radd * cos($angle) + $ysep]
     set radf [expr 0.6 * $radt]
     draw color silver
     draw cylinder "$tx $ty 0.0" "$sx $sy 0.0" radius $radf resolution 10 filled yes
     set ux [expr 1.3 * $radd * sin($angle) + $xsep - 2]
     set uy [expr 1.3 * $radd * cos($angle) + $ysep]
     draw color white
     draw text "$ux $uy 0.0" "$angcount"
     set angle [expr $angle + $twopi * 25.0 / 360.0]
     set angcount [expr $angcount + 250]
  }

  # draw caption
   set sx [expr $xsep - 2]
   set sy [expr -1.5 * $radd + $ysep]
   draw color white
   draw text "$sx $sy 0.0" "force (pN)"
}







#  loop through trajectory and animate 
#
# used for force-computing:
set kspring 10.0
set coeff 42.8

# production run
set offset 10.0         
set scale 0.1          
set num [molinfo top get numframes]
# loop through the frames
for {set i 0} {$i < $num} {incr i} {
    display update off
    draw delete all
    # go to the given frame
    animate goto $i
    set time [expr $scale * $i + $offset]
    set sel1 [atomselect top "segname PIH and name PG"]
    set sel2 [atomselect top "segname D6"]
    set coord1 [lindex [$sel1 get {x y z}] 0]
    set coord2 [lindex [$sel2 get {x y z}] 0]
    set restdist [veclength [vecsub $coord2 $coord1]]
    set force [expr $restdist * $coeff * $kspring] 
    disp_time $time 
    disp_force $force
    display update on
}

# equilibration / heatup
set offset 0.0         
set scale 0.1          
set num 100
# loop through the frames
for {set i 0} {$i < $num} {incr i} {
    display update off
    draw delete all
     # go to the given frame
    animate goto $i
    set time [expr $scale * $i + $offset]
    disp_time3 $time 
    disp_force 0.0
    display update on
}