#
# Procedures for determining the distribution of local lipid orientations
#


namespace eval ::mbtools::analysis {}

namespace eval ::mbtools::analysis::localorient {
    variable av_localorient 0
    variable av_localorient_i 0
    variable verbose
    variable f_tvslocalorient
    variable nbins
    variable range

    namespace export setup_localorient
    namespace export analyze_localorient
    namespace export printav_localorient
    namespace export resetav_localorient
}

proc ::mbtools::analysis::localorient::resetav_localorient { } {
    # Do nothing
    }

#proc ::mbtools::analysis::localorients::printav_localorient {} {
#    variable av_localorient_i
#    variable av_localorient
#    variable f_localorient
#    global ::mbtools::analysis::outputdir
#    global ::mbtools::analysis::suffix
#    variable nbins
#    variable range

#    if { $av_localorient_i > 0 } {
#    set f_localorient [open "$outputdir/av_localo$suffix" "w" ]

#        # Write a header to the file
#        puts $f_localorient "\# local orientation distribution of lipids "
#        puts $f_localorient ""


#        set lobinwidth [expr $range/($nbins*1.0)]

#        # Write out the local orientation distribution to a file
#        for { set bin 0 } { $bin < [llength $av_localorient] } { incr bin } {
#            set currbin [lindex $av_localorient $bin]
#           puts $f_localorient "[expr $bin*$lobinwidth] [expr $currbin/(1.0*[llength $topology])]"
#        }
#        close $f_localorient
#    }
#
proc ::mbtools::analysis::localorient::printav_localorient { } {
    variable av_localorient
    variable av_localorient
    variable f_tvslocalorient
    global ::mbtools::analysis::time

    puts $f_tvslocalorient "$time [expr $av_localorient/(1.0*$av_localorient_i)]"
    flush $f_tvslocalorient
}



proc ::mbtools::analysis::localorient::setup_localorient { args } {
    variable verbose
    variable range
    variable nbins
    variable av_localorient

    set options {
        {verbose "print out lots of stuff" }
        {range.arg "1.0" "Range for localorients calculation" }
        {nbins.arg "100" "Number of bins used for analysis" }
    }
    set usage "Usage: setup_pressure verbose:range "
    array set params [::cmdline::getoptions args $options $usage]

    set verbose $params(verbose)
    set range $params(range)
    set nbins $params(nbins)
    set av_localorient 0
    unset av_localorient
    #Initialize av_localorient
    for { set bn 0 } { $bn < $nbins } { incr bn } {
        lappend av_localorient 0
    }
}


# ::mbtools::analysis::analyze_localorient --
#
# Calculates the projection of the lipid orientation vector onto the
# xy plane for each lipid and then bins the absolute values of these
# vectors.
#
proc ::mbtools::analysis::localorient::analyze_localorient { printflag } {

    variable av_localorient
    variable av_localorient_i
    variable nbins
    variable range
 


    ::mmsg::send [namespace current] "analyzing local orient"

    set nbins $nbins
    set binwidth [expr $range/($nbins*1.0)]

    # Internal copy of l_orients ( same molecule order as topology )
    set l_orient [lindex [ analyze get_lipid_orients ] 1]

    set localorient [lindex [analyze lipid_orient_order all] 1]

    foreach or $localorient {
        set x [lindex $or 0]
        set y [lindex $or 1]

        # First figure out the absolute value of x and y components together
        set xyproj [expr sqrt($x*$x + $y*$y)]

        # Work out the bin to which this belongs
        if { $xyproj < $range } {
            set thisbin [expr int(floor(($xyproj)/$binwidth))]
            # Increment the relevant bin
            lset av_localorient [namespace current]bin [expr [lindex $av_localorient $thisbin] + 1]
        }



    }

    incr av_localorient_i

}

proc ::mbtools::analysis::localorient::analyze_localorient {  } {
    ::mmsg::send [namespace current] "analyzing lipid orientation order "
    variable av_localorient
    variable av_localorient_i
    variable verbose

set tmp [analyze lipid_orient_order ]

    set av_localorient [expr $av_localorient + $tmp]
    incr av_localorient_i

    if { $verbose } {
        ::mmsg::send [namespace current] "s: [expr $av_localorient/(1.0*$av_localorient_i)]"
        flush stdout
    }

    ::mmsg::debug [namespace current] "done"
}