acoustics-fit-toolkit.sh

Go to the documentation of this file.

#!/bin/zsh
script=${0##*/}

# -----------------------------------------------------------------------------------------------------
#
# Toolkit for the determination of the fixed parameters in the system, a.k.a. "pre-calibration".
# The list of parameters includes (x,y,z) positions of strings and tripods, stretching of strings,
# z-positions of the optical modules, orientations of the strings and possibly z-positions of anchors.
#
# The value of a given parameter is optimised by minimising the average chi2/NDF
# of a series of global fits applied to a predefined set of acoustic events.
# For each parameter, there exists a designated function.
#
# There are also steering functions to consistently determine the optimal values for a set of parameters.
# In these, the current accuracy of the parameter values is taken into account.
# These are therefore referred to as "stages".
# The pre-calibration thus consists of a sequence of stages.
# For each detector, a complete pre-calibration procedure is implemented in a designated script.
#
# Mandatory input variables:
#
# DETECTOR              # detector file
# TRIPOD                # tripod file
# INPUT_FILES           # list of input files (output of JAcousticsEventBuilder[.sh])
#
# Optional input variables:
#
# HYDROPHONE            # if defined, use hydrophones in global fit; else do not use hydrophones
#
# In the following:
#
# STRINGS corresponds to list of strings to be fitted; this can be changed with function fixStrings.
# TRIPODS corresponds to list of tripods to be fitted; this can be changed with function fixTripods.
#
# TRANSMITTERS corresponds to list of strings which have an emitter on the anchor (maybe empty).
# For these strings, the orientations and the z-positions of the anchors will separately be fitted.
# This list will also be changed with function fixStrings.
# -----------------------------------------------------------------------------------------------------

source $JPP_DIR/setenv.sh $JPP_DIR >& /dev/null

if do_usage $*; then
    usage "source $script"\
          "\nThis script contains auxiliary functions for the pre-calibration of specific detectors."\
          "\nIt should be sourced by the detector specific scripts."
fi

source JAcousticsToolkit.sh

set_variable  WORKDIR             `pwd`
set_variable  TMPDIR              $WORKDIR/.$$

typeset -A    TRIPODS             # tripods
typeset -a    CHI2                # chi2 values; index [1] corresponds to best value
set_variable  PRECISION  1.0E-5   # minimal change in chi2 to make step
set_variable  EPSILON    5.0E-4   # maximal distance to minimum of chi2
set_variable  RADIUS_M   1.0      # maximal horizontal distance between T-bar and emitter/hydrophone


# -----------------------------------------------------------------------------------------------------
# Initialise.
#
# Install input files in working directory and create temporary directory TMPDIR.
#
# -----------------------------------------------------------------------------------------------------
function initialise()
{
    eval `JPrintDetector -a $DETECTOR -O IDENTIFIER`
    eval `JPrintDetector -a $DETECTOR -O SUMMARY`
    eval `JPrintDetector -a $DETECTOR -O CAN`

    JAcoustics.sh  $DETECTOR_ID

    get_tripods  $TRIPOD  TRIPODS

    TRANSMITTERS=(`get_strings_with_transmitter $WORKDIR/transmitter.txt`)

    mkdir -p  $TMPDIR
}


# -----------------------------------------------------------------------------------------------------
# Remove temporary directory and restore input files in working directory.
#
# -----------------------------------------------------------------------------------------------------
function clean()
{
    rm -rf $TMPDIR >& /dev/null
    
    rm -f  {acoustics_fit_parameters,acoustics_trigger_parameters}.txt

    JAcoustics.sh  $DETECTOR_ID    
}


# -----------------------------------------------------------------------------------------------------
# Fix strings.
#
# \param  1-N     identifiers
# -----------------------------------------------------------------------------------------------------
function fixStrings()
{
    typeset -a BUFFER

    BUFFER=(`echo $*`)

    if (( ${#BUFFER} > 0 )); then
        STRINGS=(${STRINGS:|BUFFER})
        TRANSMITTERS=(${TRANSMITTERS:|BUFFER})
    fi
}


# -----------------------------------------------------------------------------------------------------
# Fix tripods.
#
# \param  1-N     identifiers
# -----------------------------------------------------------------------------------------------------
function fixTripods()
{
    typeset -a BUFFER

    BUFFER=(`echo $*`)
    
    if (( ${#BUFFER} > 0 )); then
        for ID in $BUFFER[*]; do
            unset "TRIPODS[${ID}]"
        done
    fi
}


# -----------------------------------------------------------------------------------------------------
# Get center-of-gravity of tripods.
#
# \return         x y z
# -----------------------------------------------------------------------------------------------------
function getCenter()
{
    typeset -A __TRIPODS__

    get_tripods $TRIPOD __TRIPODS__

    let "X = 0.0"
    let "Y = 0.0"
    let "Z = 0.0"

    if (( ${#__TRIPODS__} >= 1 )); then
        
        for _X _Y _Z in `echo ${(@v)__TRIPODS__}`; do
            let "X = $X + $_X"
            let "Y = $Y + $_Y"
            let "Z = $Z + $_Z"
        done

        let "X = $X / ${#__TRIPODS__}"
        let "Y = $Y / ${#__TRIPODS__}"
        let "Z = $Z / ${#__TRIPODS__}"
    fi

    printf "%15.5f %15.5f %15.5f"  $X  $Y  $Z   
}


# -----------------------------------------------------------------------------------------------------
# Set detector and tripods to given center-of-gravity.
#
# \param  1       X
# \param  2       Y
# \param  3       Z
# -----------------------------------------------------------------------------------------------------
function setCenter()
{
    getCenter | read X Y Z

    let "X = $X - $1"
    let "Y = $Y - $2"
    let "Z = $Z - $3"

    printf "Fix center %15.5f %15.5f %15.5f\n"  $X  $Y  $Z
    
    JEditTripod   -f $TRIPOD   -T "-1 sub $X $Y $Z" -o $TRIPOD   -d 0  >&  /dev/null
    JEditDetector -a $DETECTOR -S "-1 sub $X $Y $Z" -o $DETECTOR -d 0  >&  /dev/null
}


# -----------------------------------------------------------------------------------------------------
# Evaluate current chi2.
#
# \param  1       variable containing chi2 value on return
# -----------------------------------------------------------------------------------------------------
function getChi2()
{
    for (( m = 0; $m != 3; ++m )); do

        rm -f $TMPDIR/katoomba.{root,log} >& /dev/null

        date >& $TMPDIR/katoomba.log
        
        JKatoomba \
            -a $DETECTOR                              \
            -f "$INPUT_FILES[*]"                      \
            -o $TMPDIR/katoomba.root                  \
            -T $TRIPOD                                \
            -Y $WORKDIR/transmitter.txt               \
            -V $WORKDIR/sound_velocity.txt            \
            -M $WORKDIR/mechanics.txt                 \
            -@ $TMPDIR/acoustics_fit_parameters.txt   \
            -! $WORKDIR/disable.txt                   \
            ${HYDROPHONE:+-H $WORKDIR/hydrophone.txt} \
            -F 2                                      \
            -u                                        \
            -d 2 --!                                  >>& $TMPDIR/katoomba.log

        ERROR=$?

        if (( $ERROR == 0 )); then
        
            let "VALUE = $(JPrintResult -f $TMPDIR/katoomba.root:chi2 -F GetMean)" >>& $TMPDIR/katoomba.log

            ERROR=$?

            if (( $ERROR == 0 )); then

                eval $1=$VALUE

                return
            else
                printf "warning: exit code let %d\n" $ERROR
            fi
        else
            printf "warning: exit code JKatoomba %d\n" $ERROR
        fi

        # backup error status

        set_variable  UDIR  `mktemp -d $WORKDIR/XXXXXX`
        
        cp -p $TMPDIR/katoomba.log      $UDIR
        cp -p $TMPDIR/katoomba.root     $UDIR
        cp -p $DETECTOR                 $UDIR
        cp -p $TRIPOD                   $UDIR
        cp -p $WORKDIR/transmitter.txt  $UDIR

        printf "warning: backup %s\n" $UDIR
    done

    printf "warning: no valid chi2 -> exit\n"
    exit 1
}


# -----------------------------------------------------------------------------------------------------
# Fit given coordinate of string position.
#
# \param  1       string identifier
# \param  2       coordinate (X|Y|Z)
# \param  3       step
# -----------------------------------------------------------------------------------------------------
function fitPositionOfString()
{
    typeset -A BUFFER

    BUFFER=(X 0.0 Y 0.0 Z 0.0)
    
    BUFFER[$2]=$3

    if [[ "$2" == "Z" ]]; then
        OPTION=-s                  # optical modules only
    else
        OPTION=-S                  # all modules
    fi

    for (( ; $N != $NUMBER_OF_ITERATIONS; ++N )); do

        cp $DETECTOR $TMPDIR/detector.${DETECTOR:e}
        
        JEditDetector -a $DETECTOR $OPTION "$1 add ${(v)BUFFER}" -o $DETECTOR -q -d 0 >& /dev/null

        getChi2 CHI2\[2\]

        if (( $CHI2[2] >= $CHI2[1] - $PRECISION )); then

            mv $TMPDIR/detector.${DETECTOR:e} $DETECTOR

            break
        fi

        CHI2[1]=$CHI2[2]
    done

    for (( ; $N != $NUMBER_OF_ITERATIONS; ++N )); do
        
        cp $DETECTOR $TMPDIR/detector.${DETECTOR:e}

        JEditDetector -a $DETECTOR $OPTION "$1 sub ${(v)BUFFER}" -o $DETECTOR -q -d 0 >& /dev/null

        getChi2 CHI2\[2\]

        if (( $CHI2[2] >= $CHI2[1] - $PRECISION )); then

            mv $TMPDIR/detector.${DETECTOR:e} $DETECTOR

            break
        fi

        CHI2[1]=$CHI2[2]
    done

    printf "string     %04d  %s %6.4f        %6d %8.4f\n" $1 $2 $3 $N $CHI2[1]
}


# -----------------------------------------------------------------------------------------------------
# Fit stretching of string.
#
# \param  1       string identifier
# \param  2       multiplication factor
# -----------------------------------------------------------------------------------------------------
function fitStretchingOfString()
{
    let "Z = $CAN_ZMAX_M * $2"

    for (( ; $N != $NUMBER_OF_ITERATIONS; ++N )); do

        cp $DETECTOR $TMPDIR/detector.${DETECTOR:e}

        JEditDetector -a $DETECTOR -s "$1 mul $2" -s "$1 sub 0.0 0.0 $Z" -o $DETECTOR -q -d 0 >& /dev/null

        getChi2 CHI2\[2\]

        if (( $CHI2[2] >= $CHI2[1] - $PRECISION )); then

            mv $TMPDIR/detector.${DETECTOR:e} $DETECTOR
            
            break
        fi

        CHI2[1]=$CHI2[2]
    done

    for (( ; $N != $NUMBER_OF_ITERATIONS; ++N )); do
        
        cp $DETECTOR $TMPDIR/detector.${DETECTOR:e}

        JEditDetector -a $DETECTOR -s "$1 add 0.0 0.0 $Z" -s "$1 div $2" -o $DETECTOR -q -d 0 >& /dev/null

        getChi2 CHI2\[2\]

        if (( $CHI2[2] >= $CHI2[1] - $PRECISION )); then

            mv $TMPDIR/detector.${DETECTOR:e} $DETECTOR

            break
        fi

        CHI2[1]=$CHI2[2]
    done

    printf "string     %04d  M %6.4f        %6d %8.4f\n" $1 $2 $N $CHI2[1]
}


# -----------------------------------------------------------------------------------------------------
# Fit given coordinate of tripod position.
#
# \param  1       tripod identifier
# \param  2       coordinate (X|Y|Z)
# \param  3       step
# -----------------------------------------------------------------------------------------------------
function fitPositionOfTripod()
{
    typeset -A BUFFER

    BUFFER=(X 0.0 Y 0.0 Z 0.0)

    BUFFER[$2]=$3
    
    for (( ; $N != $NUMBER_OF_ITERATIONS; ++N )); do

        cp $TRIPOD $TMPDIR/tripod.txt
        
        JEditTripod -f $TRIPOD -T "$1 add ${(v)BUFFER}" -o $TRIPOD -q -d 0 >& /dev/null

        getChi2 CHI2\[2\]

        if (( $CHI2[2] >= $CHI2[1] - $PRECISION )); then

            mv $TMPDIR/tripod.txt $TRIPOD

            break
        fi

        CHI2[1]=$CHI2[2]
    done

    for (( ; $N != $NUMBER_OF_ITERATIONS; ++N )); do
        
        cp $TRIPOD $TMPDIR/tripod.txt

        JEditTripod -f $TRIPOD -T "$1 sub ${(v)BUFFER}" -o $TRIPOD -q -d 0 >& /dev/null

        getChi2 CHI2\[2\]

        if (( $CHI2[2] >= $CHI2[1] - $PRECISION )); then

            mv $TMPDIR/tripod.txt $TRIPOD

            break
        fi

        CHI2[1]=$CHI2[2]
    done

    printf "tripod       %2d  %s %6.4f        %6d %8.4f\n" $1 $2 $3 $N $CHI2[1]
}


# -----------------------------------------------------------------------------------------------------
# Fit module z-position.
#
# \param  1       string identifier
# \param  2       floor
# \param  3       step
# -----------------------------------------------------------------------------------------------------
function fitPositionOfModule()
{
    typeset -A BUFFER

    BUFFER=(X 0.0 Y 0.0 Z $3)

    set_variable  NUMBER_OF_STEPS   10       # limit number of steps
    set_variable  MODULE            `getModule -a $DETECTOR -L "$1 $2"`

    for (( i=0 ; $i != $NUMBER_OF_STEPS && $N != $NUMBER_OF_ITERATIONS; ++i, ++N )); do
        
        cp $DETECTOR $TMPDIR/detector.${DETECTOR:e}

        JEditDetector -a $DETECTOR -M "$MODULE add ${(v)BUFFER}" -o $DETECTOR -q -d 0 >& /dev/null

        getChi2 CHI2\[2\]

        if (( $CHI2[2] >= $CHI2[1] - $PRECISION )); then

            mv $TMPDIR/detector.${DETECTOR:e} $DETECTOR
            
            break
        fi

        CHI2[1]=$CHI2[2]
    done

    for (( i=0 ; $i != $NUMBER_OF_STEPS && $N != $NUMBER_OF_ITERATIONS; ++i, ++N )); do
        
        cp $DETECTOR $TMPDIR/detector.${DETECTOR:e}

        JEditDetector -a $DETECTOR -M "$MODULE sub ${(v)BUFFER}" -o $DETECTOR -q -d 0 >& /dev/null

        getChi2 CHI2\[2\]

        if (( $CHI2[2] >= $CHI2[1] - $PRECISION )); then

            mv $TMPDIR/detector.${DETECTOR:e} $DETECTOR

            break
        fi

        CHI2[1]=$CHI2[2]
    done

    printf "module     %04d.%02d %6.4f        %6d %8.4f\n" $1 $2 $3 $N $CHI2[1]
}


# -----------------------------------------------------------------------------------------------------
# Fit given coordinate of anchor position.
#
# Note that the position of the anchor corresponds to the position of the base module in the detector.
# The z-position can independently be changed whereas a change in the (x,y) position
# affects the complete string (i.e. including optical modules). 
#
# \param  1       string identifier
# \param  2       coordinate (X|Y|Z)
# \param  3       step
# -----------------------------------------------------------------------------------------------------
function fitPositionOfAnchor()
{
    typeset -A BUFFER

    BUFFER=(X 0.0 Y 0.0 Z 0.0)

    BUFFER[$2]=$3

    if [[ "$2" == "Z" ]]; then

        set_variable  MODULE  `getModule -a $DETECTOR -L "$1 0"`

        ADD="-M $MODULE add ${(v)BUFFER}"                  # base module only
        SUB="-M $MODULE sub ${(v)BUFFER}"                  # base module only
    else
        ADD="-S $1      add ${(v)BUFFER}"                  # all modules
        SUB="-S $1      sub ${(v)BUFFER}"                  # all modules
    fi

    for (( ; $N != $NUMBER_OF_ITERATIONS; ++N )); do

        cp $DETECTOR $TMPDIR/detector.${DETECTOR:e}

        JEditDetector -a $DETECTOR $ADD -o $DETECTOR -q -d 0 >& /dev/null

        getChi2 CHI2\[2\]

        if (( $CHI2[2] >= $CHI2[1] - $PRECISION )); then

            mv $TMPDIR/detector.${DETECTOR:e} $DETECTOR

            break
        fi

        CHI2[1]=$CHI2[2]
    done

    for (( ; $N != $NUMBER_OF_ITERATIONS; ++N )); do

        cp $DETECTOR $TMPDIR/detector.${DETECTOR:e}

        JEditDetector -a $DETECTOR $SUB -o $DETECTOR -q -d 0 >& /dev/null

        getChi2 CHI2\[2\]

        if (( $CHI2[2] >= $CHI2[1] - $PRECISION )); then

            mv $TMPDIR/detector.${DETECTOR:e} $DETECTOR

            break
        fi

        CHI2[1]=$CHI2[2]
    done

    printf "anchor     %04d  %s %6.4f        %6d %8.4f\n" $1 $2 $3 $N $CHI2[1]
}


# -----------------------------------------------------------------------------------------------------
# Fit rotation of anchor.
#
# Note that a rotation of the anchor corresponds to a change of the relative positions of
# the hydrophone and emitter with respect to the T-bar.
#
# \param  1       string identifier
# \param  2       rotation angle [rad]
# -----------------------------------------------------------------------------------------------------
function fitRotationOfAnchor()
{
    let "ROT = $2"

    for (( ; $N != $NUMBER_OF_ITERATIONS; ++N )); do

        cp $WORKDIR/hydrophone.txt  $TMPDIR/hydrophone.txt
        cp $WORKDIR/transmitter.txt $TMPDIR/transmitter.txt
        
        JEditHydrophone  -f $WORKDIR/hydrophone.txt  -S "$1 rot $ROT" -o $WORKDIR/hydrophone.txt  -q -d 0 >& /dev/null
        JEditTransmitter -f $WORKDIR/transmitter.txt -S "$1 rot $ROT" -o $WORKDIR/transmitter.txt -q -d 0 >& /dev/null

        getChi2 CHI2\[2\]

        if (( $CHI2[2] >= $CHI2[1] - $PRECISION )); then

            mv $TMPDIR/hydrophone.txt  $WORKDIR/hydrophone.txt
            mv $TMPDIR/transmitter.txt $WORKDIR/transmitter.txt

            break
        fi

        CHI2[1]=$CHI2[2]
    done

    let "ROT = -1.0 * $ROT"

    for (( ; $N != $NUMBER_OF_ITERATIONS; ++N )); do

        cp $WORKDIR/hydrophone.txt  $TMPDIR/hydrophone.txt
        cp $WORKDIR/transmitter.txt $TMPDIR/transmitter.txt
        
        JEditHydrophone  -f $WORKDIR/hydrophone.txt  -S "$1 rot $ROT" -o $WORKDIR/hydrophone.txt  -q -d 0 >& /dev/null
        JEditTransmitter -f $WORKDIR/transmitter.txt -S "$1 rot $ROT" -o $WORKDIR/transmitter.txt -q -d 0 >& /dev/null

        getChi2 CHI2\[2\]

        if (( $CHI2[2] >= $CHI2[1] - $PRECISION )); then

            mv $TMPDIR/hydrophone.txt  $WORKDIR/hydrophone.txt
            mv $TMPDIR/transmitter.txt $WORKDIR/transmitter.txt

            break
        fi

        CHI2[1]=$CHI2[2]
    done

    printf "string     %04d  R %6.4f        %6d %8.4f\n" $1 $2 $N $CHI2[1]
}


# -----------------------------------------------------------------------------------------------------
# Stage 0.
# This stage can be used to determine the (x,y,z) positions of tripods when a number of strings is fixed.
#
# -----------------------------------------------------------------------------------------------------
function stage_0()
{
    cat>$TMPDIR/acoustics_fit_parameters.txt<<EOF
`egrep Tmax_s $WORKDIR/acoustics_fit_parameters.txt`
Nmin        =   3;
sigma_s     = 100.0e-6;
stdev       =  10.0;
mestimator  =   0;
option      =   1;
EOF

    if (( ${#STRINGS} > 0 )); then
        JEditDetector  -a $DETECTOR  -o $TMPDIR/detector_0.datx  -k "$STRINGS[*]" -q -d 0 >& /dev/null
        JEditDetector  -a $DETECTOR  -o $DETECTOR                -r "$STRINGS[*]" -q -d 0 >& /dev/null
    fi
    
    getChi2 CHI2\[1\]

    set_variable  NUMBER_OF_ITERATIONS       1000

    for DX_M in 0.5 0.2 0.1; do              # determine (x,y,z) positions of tripods

        let "A_RAD = $DX_M / $RADIUS_M"      # use maximal horizontal distance between T-bar and emitter/hydrophone

        for (( N = 0; $N < $NUMBER_OF_ITERATIONS; ++N )); do

            CHI2[3]=$CHI2[1]
            
            for ID in ${(k)TRIPODS}; do
                fitPositionOfTripod $ID X $DX_M
                fitPositionOfTripod $ID Y $DX_M
                fitPositionOfTripod $ID Z $DX_M
            done

            for STRING in $TRANSMITTERS[*]; do  
                fitPositionOfAnchor $STRING X $DX_M
                fitPositionOfAnchor $STRING Y $DX_M
                fitPositionOfAnchor $STRING Z $DX_M
                fitRotationOfAnchor $STRING $A_RAD
            done

            if (( $CHI2[3] - $CHI2[1] < $EPSILON )); then
                break
            fi
        done

        if (( $N >= $NUMBER_OF_ITERATIONS )); then
            printf "warning: reached maximum number of iterations %d - convergence %7.3f\n" $N $(($CHI2[3] - $CHI2[1]))
        fi
    done

    if (( ${#STRINGS} > 0 )); then
        JMergeDetector  -a $DETECTOR  -a $TMPDIR/detector_0.datx  -o $DETECTOR >& /dev/null
    fi
}


# -----------------------------------------------------------------------------------------------------
# Stage 1A.
# This stage can be used to roughly determine the positions of the strings, tripods and modules.
# During this procedure, the strings are kept vertical in the global fits.
#
# -----------------------------------------------------------------------------------------------------
function stage_1A()
{
    cat>$TMPDIR/acoustics_fit_parameters.txt<<EOF
`egrep Tmax_s $WORKDIR/acoustics_fit_parameters.txt`
Nmin        =   3;
sigma_s     = 250.0e-6;
stdev       =  10.0;
mestimator  =   2;
option      =   0;
EOF

    getChi2 CHI2\[1\]

    set_variable  NUMBER_OF_ITERATIONS       1000

    for DX_M in 0.5 0.2; do                  # determine (x,y,z) positions of strings and tripods

        let "A_RAD = $DX_M / $RADIUS_M"      # use maximal horizontal distance between T-bar and emitter/hydrophone

        for (( N = 0; $N < $NUMBER_OF_ITERATIONS; ++N )); do

            CHI2[3]=$CHI2[1]
            
            for STRING in $STRINGS[*]; do       
                fitPositionOfString $STRING X $DX_M
                fitPositionOfString $STRING Y $DX_M
                fitPositionOfString $STRING Z $DX_M
            done

            for ID in ${(k)TRIPODS}; do
                fitPositionOfTripod $ID X $DX_M
                fitPositionOfTripod $ID Y $DX_M
                fitPositionOfTripod $ID Z $DX_M
            done

            for STRING in $TRANSMITTERS[*]; do  
                fitPositionOfAnchor $STRING X $DX_M
                fitPositionOfAnchor $STRING Y $DX_M
                fitPositionOfAnchor $STRING Z $DX_M
                fitRotationOfAnchor $STRING $A_RAD
            done

            if (( $CHI2[3] - $CHI2[1] < $EPSILON )); then
                break
            fi
        done

        if (( $N >= $NUMBER_OF_ITERATIONS )); then
            printf "warning: reached maximum number of iterations %d - convergence %7.3f\n" $N $(($CHI2[3] - $CHI2[1]))
        fi
    done

    getChi2 CHI2\[1\]

    set_variable  NUMBER_OF_ITERATIONS       2000

    for DX_M in 0.1; do                      # determine (z) positions of modules
        
        for (( N = 0; $N < $NUMBER_OF_ITERATIONS; ++N )); do
            
            CHI2[3]=$CHI2[1]

            for STRING in $STRINGS[*]; do
                for (( FLOOR = 1; $FLOOR <= 18; FLOOR += 1 )); do
                    fitPositionOfModule $STRING $FLOOR $DX_M
                done
            done

            if (( $CHI2[3] - $CHI2[1] < $EPSILON )); then
                break
            fi
        done

        if (( $N >= $NUMBER_OF_ITERATIONS )); then
            printf "warning: reached maximum number of iterations %d - convergence %7.3f\n" $N $(($CHI2[3] - $CHI2[1]))
        fi
    done
}


# -----------------------------------------------------------------------------------------------------
# Stage 1B.
# This stage can be used to determine the z-positions and stretching of individual strings.
#
# -----------------------------------------------------------------------------------------------------
function stage_1B()
{
    cat>$TMPDIR/acoustics_fit_parameters.txt<<EOF
`egrep Tmax_s $WORKDIR/acoustics_fit_parameters.txt`
Nmin        =   3;
sigma_s     = 100.0e-6;
stdev       =  10.0;
mestimator  =   0;
option      =   1;
EOF
 
    for STRING in $STRINGS[*]; do            # determine stretching and (z) position of given string          
        
        JEditDetector  -a $DETECTOR  -o $TMPDIR/detector_1.datx  -r "$STRING" -q -d 0 >& /dev/null
        JEditDetector  -a $DETECTOR  -o $DETECTOR                -k "$STRING" -q -d 0 >& /dev/null

        getChi2 CHI2\[1\]

        set_variable  NUMBER_OF_ITERATIONS       1000

        for MUL in 0.005 0.001; do

            DX_M=0.2

            for (( N = 0; $N < $NUMBER_OF_ITERATIONS; ++N )); do

                CHI2[3]=$CHI2[1]

                fitPositionOfString   $STRING Z  $DX_M
                fitStretchingOfString $STRING $MUL

                if (( $CHI2[3] - $CHI2[1] < $EPSILON )); then
                    break
                fi
            done

            if (( $N >= $NUMBER_OF_ITERATIONS )); then
                printf "warning: reached maximum number of iterations %d - convergence %7.3f\n" $N $(($CHI2[3] - $CHI2[1]))
            fi
        done

        JMergeDetector  -a $DETECTOR  -a $TMPDIR/detector_1.datx  -o $DETECTOR >& /dev/null
    done
}


# -----------------------------------------------------------------------------------------------------
# Stage 2.
# This stage can be used to improve the determination the positions of the strings, tripods and modules.
#
# -----------------------------------------------------------------------------------------------------
function stage_2()
{
    cat>$TMPDIR/acoustics_fit_parameters.txt<<EOF
`egrep Tmax_s $WORKDIR/acoustics_fit_parameters.txt`
Nmin        =   3;
sigma_s     = 100.0e-6;
stdev       =  10.0;
mestimator  =   0;
option      =   1;
EOF

    getChi2 CHI2\[1\]

    set_variable  NUMBER_OF_ITERATIONS       1000

    for DX_M in 0.2; do                      # determine (x,y,z) positions of strings and tripods

        let "A_RAD = $DX_M / $RADIUS_M"      # use maximal horizontal distance between T-bar and emitter/hydrophone
        
        for (( N = 0; $N < $NUMBER_OF_ITERATIONS; ++N )); do

            CHI2[3]=$CHI2[1]
            
            for STRING in $STRINGS[*]; do       
                fitPositionOfString $STRING X $DX_M
                fitPositionOfString $STRING Y $DX_M
                fitPositionOfString $STRING Z $DX_M
            done

            for ID in ${(k)TRIPODS}; do
                fitPositionOfTripod $ID X $DX_M
                fitPositionOfTripod $ID Y $DX_M
                fitPositionOfTripod $ID Z $DX_M
            done

            for STRING in $TRANSMITTERS[*]; do  
                fitPositionOfAnchor $STRING X $DX_M
                fitPositionOfAnchor $STRING Y $DX_M
                fitPositionOfAnchor $STRING Z $DX_M
                fitRotationOfAnchor $STRING $A_RAD
            done

            if (( $CHI2[3] - $CHI2[1] < $EPSILON )); then
                break
            fi
        done

        if (( $N >= $NUMBER_OF_ITERATIONS )); then
            warning "reached maximum number of iterations $N - convergence $(($CHI2[3] - $CHI2[1]))"
        fi
    done

    getChi2 CHI2\[1\]

    set_variable  NUMBER_OF_ITERATIONS       2000

    for DX_M in 0.1; do                      # determine (z) positions of modules and (x,y) positions of strings

        for (( N = 0; $N < $NUMBER_OF_ITERATIONS; ++N )); do

            CHI2[3]=$CHI2[1]
            
            for STRING in $STRINGS[*]; do       

                for (( FLOOR = 1; $FLOOR <= 18; FLOOR += 1 )); do
                    fitPositionOfModule $STRING $FLOOR $DX_M
                done

                fitPositionOfString $STRING X $DX_M
                fitPositionOfString $STRING Y $DX_M
            done

            if (( $CHI2[3] - $CHI2[1] < $EPSILON )); then
                break
            fi
        done

        if (( $N >= $NUMBER_OF_ITERATIONS )); then
            warning "reached maximum number of iterations $N - convergence $(($CHI2[3] - $CHI2[1]))"
        fi
    done
}

# -----------------------------------------------------------------------------------------------------
# Stage 3.
# This stage can be used to finalise the determination of the positions of the strings, tripods and modules.
#
# -----------------------------------------------------------------------------------------------------
function stage_3()
{
    cat>$TMPDIR/acoustics_fit_parameters.txt<<EOF
`egrep Tmax_s $WORKDIR/acoustics_fit_parameters.txt`
Nmin        =   3;
sigma_s     =  50.0e-6;
stdev       =  10.0;
mestimator  =   0;
option      =   1;
EOF

    getChi2 CHI2\[1\]

    set_variable  NUMBER_OF_ITERATIONS       1000

    for DX_M in 0.10 0.05; do                # determine (x,y,z) positions of strings and tripods

        let "A_RAD = $DX_M / $RADIUS_M"      # use maximal horizontal distance between T-bar and emitter/hydrophone
        
        for (( N = 0; $N < $NUMBER_OF_ITERATIONS; ++N )); do

            CHI2[3]=$CHI2[1]
            
            for STRING in $STRINGS[*]; do       
                fitPositionOfString $STRING X $DX_M
                fitPositionOfString $STRING Y $DX_M
                fitPositionOfString $STRING Z $DX_M
            done

            for ID in ${(k)TRIPODS}; do
                fitPositionOfTripod $ID X $DX_M
                fitPositionOfTripod $ID Y $DX_M
                fitPositionOfTripod $ID Z $DX_M
            done

            for STRING in $TRANSMITTERS[*]; do  
                fitPositionOfAnchor $STRING X $DX_M
                fitPositionOfAnchor $STRING Y $DX_M
                fitPositionOfAnchor $STRING Z $DX_M
                fitRotationOfAnchor $STRING $A_RAD
            done

            if (( $CHI2[3] - $CHI2[1] < $EPSILON )); then
                break
            fi
        done

        if (( $N >= $NUMBER_OF_ITERATIONS )); then
            warning "reached maximum number of iterations $N - convergence $(($CHI2[3] - $CHI2[1]))"
        fi
    done

    getChi2 CHI2\[1\]

    set_variable  NUMBER_OF_ITERATIONS       2000

    for DX_M in 0.05; do                     # determine (z) positions of modules and (x,y) positions of strings

        for (( N = 0; $N < $NUMBER_OF_ITERATIONS; ++N )); do

            CHI2[3]=$CHI2[1]
            
            for STRING in $STRINGS[*]; do       

                for (( FLOOR = 1; $FLOOR <= 18; FLOOR += 1 )); do
                    fitPositionOfModule $STRING $FLOOR $DX_M
                done

                fitPositionOfString $STRING X $DX_M
                fitPositionOfString $STRING Y $DX_M
            done

            if (( $CHI2[3] - $CHI2[1] < $EPSILON )); then
                break
            fi
        done

        if (( $N >= $NUMBER_OF_ITERATIONS )); then
            warning "reached maximum number of iterations $N - convergence $(($CHI2[3] - $CHI2[1]))"
        fi
    done
}


# -----------------------------------------------------------------------------------------------------
# Stage D.
# This stage can be used to determine the rotations of the strings and z-positions of the anchors.
# It is followed by a determination of the positions of the strings, tripods and modules.
#
# -----------------------------------------------------------------------------------------------------
function stage_D()
{
    cat>$TMPDIR/acoustics_fit_parameters.txt<<EOF
`egrep Tmax_s $WORKDIR/acoustics_fit_parameters.txt`
Nmin        =   3;
sigma_s     =  50.0e-6;
stdev       =  10.0;
mestimator  =   0;
option      =   1;
EOF

    getChi2 CHI2\[1\]

    set_variable  HYDROPHONE                 
    set_variable  NUMBER_OF_ITERATIONS       2000

    for DX_M in 0.10; do                     # determine rotation of string and z-position of anchors

        let "A_RAD = $DX_M / $RADIUS_M"      # use maximal horizontal distance between T-bar and emitter/hydrophone

        for (( N = 0; $N < $NUMBER_OF_ITERATIONS; ++N )); do

            CHI2[3]=$CHI2[1]

            for STRING in $STRINGS[*]; do
                fitRotationOfAnchor $STRING $A_RAD
                fitPositionOfAnchor $STRING Z $DX_M
            done

            if (( $CHI2[3] - $CHI2[1] < $EPSILON )); then
                break
            fi
        done

        if (( $N >= $NUMBER_OF_ITERATIONS )); then
            warning "reached maximum number of iterations $N - convergence $(($CHI2[3] - $CHI2[1]))"
        fi
    done

    getChi2 CHI2\[1\]

    set_variable  NUMBER_OF_ITERATIONS       2000

    for DX_M in 0.05; do                       # determine position of strings and tripods

        for (( N = 0; $N < $NUMBER_OF_ITERATIONS; ++N )); do

            CHI2[3]=$CHI2[1]

            for STRING in $STRINGS[*]; do
                fitPositionOfString $STRING X $DX_M
                fitPositionOfString $STRING Y $DX_M
                fitPositionOfString $STRING Z $DX_M
                fitPositionOfAnchor $STRING Z $DX_M
            done

            for ID in ${(k)TRIPODS}; do
                fitPositionOfTripod $ID X $DX_M
                fitPositionOfTripod $ID Y $DX_M
                fitPositionOfTripod $ID Z $DX_M
            done

            if (( $CHI2[3] - $CHI2[1] < $EPSILON )); then
                break
            fi
        done

        if (( $N >= $NUMBER_OF_ITERATIONS )); then
            warning "reached maximum number of iterations $N - convergence $(($CHI2[3] - $CHI2[1]))"
        fi
    done

    getChi2 CHI2\[1\]

    set_variable  NUMBER_OF_ITERATIONS       1000

    for DX_M in 0.05; do                     # determine (z) positions of modules and (x,y) positions of strings

        for (( N = 0; $N < $NUMBER_OF_ITERATIONS; ++N )); do

            CHI2[3]=$CHI2[1]

            for STRING in $STRINGS[*]; do

                for (( FLOOR = 1; $FLOOR <= 18; FLOOR += 1 )); do
                    fitPositionOfModule $STRING $FLOOR $DX_M
                done

                fitPositionOfString $STRING X $DX_M
                fitPositionOfString $STRING Y $DX_M
            done

            if (( $CHI2[3] - $CHI2[1] < $EPSILON )); then
                break
            fi
        done

        if (( $N >= $NUMBER_OF_ITERATIONS )); then
            warning "reached maximum number of iterations $N - convergence $(($CHI2[3] - $CHI2[1]))"
        fi
    done
}