package require La
set outfile [open axis.dat w];                                                           
set sel [atomselect top "not water"]               
set n [molinfo top get numframes]                                                              
for { set i 0 } { $i <= $n } { incr i } {                                                                                                                                          
    set weights [ $sel get mass ]                                                
    set x [ $sel get x ]                                                              
    set y [ $sel get y ]                                                              
    set z [ $sel get z ]                                                              
    set m $weights                                                                    
                                                                                           
    set comx 0                                                                        
    set comy 0                                                                        
    set comz 0                                                                        
    set totalm 0 
    set COM [list $comx $comy $comz]                                                                     
    foreach xx $x yy $y zz $z mm $m {                                                 
        # use the abs of the weights                                                  
        set mm [expr abs($mm)]                                                        
	set comx [ expr "$comx + $xx*$mm" ]                                           
	set comy [ expr "$comy + $yy*$mm" ]                                           
	set comz [ expr "$comz + $zz*$mm" ]                                           
	set totalm [ expr "$totalm + $mm" ]                                           
    }                                                                                 
    set comx [ expr "$comx / $totalm" ]                                               
    set comy [ expr "$comy / $totalm" ]                                               
    set comz [ expr "$comz / $totalm" ]                                                                                           
    set x [ $sel get x ]                                                              
    set y [ $sel get y ]                                                              
    set z [ $sel get z ]                                                              
    set m $weights                                                                    
                                                                                      
    # compute I                                  hah                                     
    set Ixx 0                                                                         
    set Ixy 0                                                                         
    set Ixz 0                                                                         
    set Iyy 0                                                                         
    set Iyz 0                                                                         
    set Izz 0                                                                         
    foreach xx $x yy $y zz $z mm $m {                                                 
        # use the abs of the weights                                                  
       set mm [expr abs($mm)]                                                        
                                                                                      
        # subtract the COM                                                            
        set xx [expr $xx - [lindex $COM 0]]                                           
        set yy [expr $yy - [lindex $COM 1]]                                           
        set zz [expr $zz - [lindex $COM 2]]                                           
                                                                                      
        set rr [expr $xx + $yy + $zz]                                                 
                                                                                      
        set Ixx [expr $Ixx + $mm*($yy*$yy+$zz*$zz)]                                   
        set Ixy [expr $Ixy - $mm*($xx*$yy)]                                           
        set Ixz [expr $Ixz - $mm*($xx*$zz)]                                           
        set Iyy [expr $Iyy + $mm*($xx*$xx+$zz*$zz)]                                   
        set Iyz [expr $Iyz - $mm*($yy*$zz)]                                           
        set Izz [expr $Izz + $mm*($xx*$xx+$yy*$yy)]  
                                          
     set   I  [list 2 3 3 $Ixx $Ixy $Ixz $Ixy $Iyy $Iyz $Ixz $Iyz $Izz]
    
     set a1   " [lindex $I 3] [lindex $I 6] [lindex $I 9]" 
     set a2   " [lindex $I 4] [lindex $I 7] [lindex $I 10]"  
     set a3   " [lindex $I 5] [lindex $I 8] [lindex $I 11]" 
     La::mevsvd_br I evals  
    set vector1  " [lindex $I 3] [lindex $I 6] [lindex $I 9]"  
    set a        [vecnorm $vector1]
    set b    {0 0 1}
    
    set vec1 [vecnorm $vector1]                
                                               
    set b         { 0  0  1}                                                       
    set amag     [veclength $a]                                                    
    set bmag      [veclength $b ]                                                  
    set dotprod  [vecdot $a $b]                                                    
    set angle    [expr 57.2958 * acos($dotprod / ($amag * $bmag))]                
                                                    
                                                           
    }    
    
  $sel frame $i                                                                                                      
  $sel update        
  puts   $outfile " [expr {$i+1}]            $angle            {$a}     {$I}    "  
 }                                             
                                                     
  close $outfile