#!/bin/bash
# -*-shell-script-*-

# Extracts various data from Starlab log files.
# Copyright (C) 2005, 2007 Ernest N. Mamikonyan
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.


action=com_position
force=false
output=/dev/stdout
quiet=false


# enable colors but only when output is a tty
if [[ -t 2 ]]; then
    ec='\e[31m' # error color (red)
    wc='\e[33m' # warning color (yellow)
    sc='\e[32m' # success color (green)
    fc='\e[35m' # file color (magenta)
    nc='\e[39m' # normal color
fi


error() {
    builtin echo -n "${0##*/}: $ec$1" >&2
    [[ "$2" ]] && builtin echo " - exiting...$nc" >&2 && exit $2
    builtin echo "$nc" >&2
}


warn() { echo "${0##*/}: $wc$@$nc" >&2; }


usage() {
    echo "\
Usage: ${0##*/} [OPTIONS] [LOG FILE] ...
Extracts various data from Starlab log files.
Defaults are shown in parentheses.

  -a	center of density position, X-Y projection
  -b	center of density galacto-centric distance
  -c	center of mass position, X-Y projection (default)
  -d	center of mass galacto-centric distance
  -D DF look in the associated dyn file, DF, for scale factors
  -f	overwrite existing files silently ($force)
  -h	display this help message and exit
  -M	mass loss due to dynamical evaporation and stellar evolution;
	output will be 3 columns: time, dm_escapers, dm_massloss
  -m	mass (in stars) bound to the cluster
  -n	number of stars bound to the cluster
  -o F	output to file F ($output)
  -r	cluster half-mass radius
  -q	suppress all diagnostic output
  -S SF assume the scale factors SF given in the order: mass, length, time,
	i.e., -S '91902.6 0.039 0.000377257' (search input)
  -V	display version information and exit
" >&2
    exit $1
}


to_clobber() {
    [[ -e "$1" ]] || return 0
    if $force; then warn "$fc$1$nc exists - ${wc}overwriting$nc..."; return 0
    else
        read -n1 -p"$(builtin echo -n $fc$1$nc exists - overwrite\?) " </dev/tty
        if [[ "$REPLY" != [Yy] ]]; then
            if [[ "$REPLY" ]]; then builtin echo "\b${wc}No$nc" >&2
	    else builtin echo "${wc}No$nc" >&2
	    fi
            return 1
        fi
        builtin echo "\b${ec}Yes$nc" >&2
        return 0
    fi
}


shopt -s expand_aliases nullglob xpg_echo

while getopts :abcD:dfhMmno:qrS:V option; do
    case $option in
	a) action=cod_position;;
	b) action=cod_distance;;
	c) action=com_position;;
	d) action=com_distance;;
	D) dyn="$OPTARG";;
	f) force=true;;
	h) usage 0;;
	M) action=mass_loss;;
	m) action=mass;;
	n) action=star_number;;
	o) output="$OPTARG";;
	q) quiet=true;;
	r) action=rhalf;;
	S) sf=($OPTARG);;
	V) echo $'\
log_extract version 0.6, Copyright  2005 Ernest N. Mamikonyan.
log_extract comes with ABSOLUTELY NO WARRANTY.
This is free software, and you are welcome to redistribute it
under certain conditions; see the source for details.
'
	    exit 0
	    ;;
	:) echo "The \`-$OPTARG' option requires an argument!" >&2; usage 1;;
	*) echo "\`-$OPTARG' is not a valid option!" >&2; usage 1;;
    esac
done

shift $((OPTIND-1))


$quiet && echo() { :; } || echo() { builtin echo "$@" >&2; }


# scale factors in this order: mass, length, time
# e.g., mass: ${sf[0]}, length: ${sf[1]}, time: ${sf[2]}
if [[ "$sf" ]]; then
    echo "\
accepted the following scale factors,
mass: ${sf[0]}, length: ${sf[1]}, time: ${sf[2]}"
elif [[ "$dyn" ]]; then
    echo -n "searching for scale factors in $fc$dyn$nc..."
    sf=($(awk '
BEGIN {
  while (i < 3 && getline) {
    if (/^[[:blank:]]*mass_scale[[:blank:]]*=/) { ++i; print 1/$3 }
    else if (/^[[:blank:]]*size_scale[[:blank:]]*=/) { ++i; print 2.255e-8/$3 }
    else if (/^[[:blank:]]*time_scale[[:blank:]]*=/) { ++i; print 1/$3 }
  }
}' "$dyn"))
else
    echo -n "searching for scale factors in input..."
    sf=($(awk '
BEGIN {
  while (i < 3 && getline) {
    if (/^[[:blank:]]*\[m\]:/) { ++i; print $2 }
    else if (/^[[:blank:]]*\[R\]:/) { ++i; print $2 }
    else if (/^[[:blank:]]*\[T\]:/) { ++i; print $2 }
  }
}' "$@"))
fi

if [[ "$sf" ]]; then
    echo "${sc}done$nc\nmass: ${sf[0]}, length: ${sf[1]}, time: ${sf[2]}"
else
    echo "${ec}failed$nc"
    warn "can't find scale factors - continuing"
    sf=(1 1 1)
fi


[[ "$output" == /dev/stdout ]] || to_clobber "$output" || exit


case $action in
    cod_distance)
	echo -n computing the center of density galacto-centric distance...
	awk "
/^Time = / { t = ${sf[2]}*\$3 }
/density_center = / {
  print t, sqrt(\$3^2+\$4^2+\$5^2)*${sf[1]}
}" "$@" > "$output"
	echo done
	;;
    cod_position)
	echo -n extracting the center of density position \(X-Y plane\)...
	awk "
/density_center = / {
  print \$3*${sf[1]}, \$4*${sf[1]}
}" "$@" > "$output"
	echo done
	;;
    com_distance)
	echo -n computing the center of mass galacto-centric distance...
	awk "
/^Time = / { t = ${sf[2]}*\$3 }
/center of mass position = / {
  print t, sqrt(\$6^2+\$7^2+\$8^2)*${sf[1]}
}" "$@" > "$output"
	echo done
	;;
    com_position)
	echo -n extracting the center of mass position \(X-Y plane\)...
	awk "
/center of mass position = / {
  print \$6*${sf[1]}, \$7*${sf[1]}
}" "$@" > "$output"
	echo done
	;;
    mass_loss)
	echo -n extracting mass loss due to evaporation and stellar evolution...
	awk "
/^Time = / { t = ${sf[2]}*\$3 }
/dm_escapers = / { me = ${sf[0]}*\$3 }
/dm_massloss = / {
  print t, me, ${sf[0]}*\$3
}" "$@" > "$output"
	echo done
	;;
    mass)
	echo -n computing the bound mass...
	awk "
/^Time = / {
  print \$3*${sf[2]}, substr(\$12,2)*${sf[0]}
}" "$@" > "$output"
	echo done
	;;
    rhalf)
	echo -n extracting the half-mass radius...
	awk "
/^Time = / { t = ${sf[2]}*\$3 }
/rhalf = / {
  print t, \$3*${sf[1]}
}" "$@" > "$output"
	echo done
	;;
    star_number)
	echo -n extracting the number of bound stars...
	awk "
/^Time = / {
  print \$3*${sf[2]}, substr(\$7,2)
}" "$@" > "$output"
	echo done
	;;
esac
