#! /bin/csh -fe
#
#  Automatically generated scaling script for:
#
Help:
if(("" == "")||($?HELP)) then
    cat << EOF

usage: $0 raw.mtz [P21212] [1.7A] [aa]

where: P21212 	- is the space group you want to "reindex" to.
       1.7A	- is the outer resolution cutoff
       aa	- is the number of amino acids in your asymmetric unit

EOF
    if($?HELP) exit
endif

alias nawk /usr/bin/nawk
##############################################################################
set rawfile = raw.mtz
set outfile = merged.mtz
set tempfile = temp/Wedger.temp

##############################################################################
# this value was estimated for a Vm of 2.4 in P21212
set ASU  = 175
set SG   = ""
set CELL = ""
set RES  = 1.7

##############################################################################
goto Setup
# interpret command line
Return_From_Setup:
##############################################################################

echo "scaling, mergeing and truncating $rawfile into $outfile"
echo "to ${RES}A in $SG with ${ASU}aa per ASU"


##############################################################################
# make sure the output is sorted
#
##############################################################################
sortmtz hklout ${tempfile}.sorted.mtz  << end_sort
H K L M/ISYM BATCH I SIGI
$rawfile
end_sort
#
##############################################################################
# now run SCALA 
#
##############################################################################
scala hklin ${tempfile}.sorted.mtz hklout ${tempfile}.scaled.mtz  \
ROGUES rejected_spots.txt \
ANOMPLOT ${tempfile}.anomplot \
NORMPLOT ${tempfile}.normplot \
SCALES ${tempfile}.scales << EOF-scale
# save disk space?
#NODUMP
#ANALYSE NOPLOT

title "batch scaling of "

resolution $RES

scales batch bfactor on

# straight frame-by-frame scales (for discontinuities in dose per frame)
#scales batch bfactor on
# frame-by-frame scales with smooth B-factor (discontinuous dose, but no xtal decay)
#scales batch bfactor on brotation spacing 5 
# smooth scales with smooth B-factor (if you have no discontinuities in dose per frame)
#scales rotation spacing 5 bfactor on
#
# restrain scales
#tie rotation 0.1
#
# extend nearly-full partials (this is important if you have no fulls)
#intensities scale_partials 0.8
# or use summed partials as fulls:
#intensities partials
#
#########################################################
#  See SCALA documentation
SDCORR 1.3 0 0.03
#########################################################
# 
# converge to within 0.3 sigmas of deviation (maximum of 20 cycles)
cycles 20 converge 0.3

anomalous on

EOF-scale
#
##############################################################################
# run truncate to finish the job
#
##############################################################################
truncate:
truncate hklin ${tempfile}.scaled.mtz hklout ${tempfile}.merged.mtz << EOF-trunc
title "anomalous truncate of /fourier/lougheed/frames/miapt_2 frames 1 to 107"
# resolution used to determine B factor
#rscale 4 2.8
#ranges 0.01

nresidue $ASU
truncate yes
anomalous yes
LABOUT  F=F  SIGF=SIGF DANO=DANO SIGDANO=SIGDANO ISYM=ISYM
EOF-trunc

##############################################################################
# make sure completeness is correctly reported in mtz
#
##############################################################################
unique hklout ${tempfile}.unique.mtz << EOF-unique
SYMM $SGnum
CELL $CELL
RESOLUTION $RES
LABOUT  F=F  SIGF=SIGF
EOF-unique
mtzutils hklin1 ${tempfile}.merged.mtz \
         hklin2 ${tempfile}.unique.mtz \
         hklout $outfile << EOF-util
UNIQUE
RUN
EOF-util

# clean up
rm -f ${tempfile}.merged.mtz >& /dev/null
rm -f ${tempfile}.unique.mtz >& /dev/null
rm -f ${tempfile}.scaled.mtz >& /dev/null
rm -f ${tempfile}.sorted.mtz >& /dev/null
rm -f ${tempfile}.anomplot   >& /dev/null
rm -f ${tempfile}.normplot   >& /dev/null
rm -f ${tempfile}.${SG}.mtz  >& /dev/null
exit

##############################################################################
##############################################################################
##############################################################################
##############################################################################
Setup:
# command-line arguments
foreach arg ( $* )
    if("$arg" =~ *.mtz) then
	if($?usermtz) then
	    set rawfile = "$arg"
	else
	    set rawfile = "$rawfile $arg"
	endif
	set usermtz
    endif
    if(("$arg" =~ [0-9]*)&&("$arg" =~ *A)) then
        # a number ending with "A" must mean resolution
        set RES = `echo "$arg" | nawk '{print $1+0}'`
	continue
    endif
    if(("$arg" =~ [0-9]*)&&("$arg" !~ [a-z,A-Z]*)) then
	# a pure number is on the command line
	if("$arg" =~ *.*) then
	    # take it as resolution
	    set RES = `echo "$arg" | nawk '{print $1+0}'`
	else
	    # take it as ASU size
	    set ASU = `echo "$arg" | nawk '{print $1+0}'`
	endif
    endif
    
    # new space group
    if("$arg" =~ [PpCcIiFfRr][1-6]*) then
	set temp = `nawk -v SG=$arg '$4 == toupper(SG) && $1<500{print $4}' $CLIBD/symop.lib | head -1`
	if("$temp" =~ [PpCcIiFfRr][1-6]*) then
	    # known space group
	    set new_SG = "$temp"
	else
	    # check for "pseudo-spacegroup" language
	    if("$arg" =~ [Pp]2212) then
		# P2221 with screw along current "b"
		set new_SG = "P2221"
		set REINDEX = "reindex k, l, h"
	    endif
	    if("$arg" =~ [Pp]2122) then
		# P2221 with screw along current "a"
		set new_SG = "P2221"
		set REINDEX = "reindex l, h, k"
	    endif
	    if("$arg" =~ [Pp]21221) then
		# P21212 with non-screw along current "b"
		set new_SG = "P21212"
		set REINDEX = "reindex k, l, h"
	    endif
	    if("$arg" =~ [Pp]22121) then
		# P21212 with non-screw along current "a"
	    	set new_SG = "P21212"
		set REINDEX = "reindex l, h, k"
	    endif
	    if("$arg" =~ [Cc]2212) then
		# C2221 with screw along current "b"
		set new_SG = "C2221"
		set REINDEX = "reindex k, l, h"
	    endif
	    if("$arg" =~ [Cc]2122) then
		# C2221 with screw along current "a"
		set new_SG = "C2221"
		set REINDEX = "reindex l, h, k"
	    endif
	endif
    endif
end


##############################################################################
# check input MTZ file
if(! -e "$rawfile") then
    set HELP
    goto Help
endif

##############################################################################
# update cell and SG variables from the raw MTZ header
#
##############################################################################
echo "HEADER" | mtzdump hklin $rawfile >! $tempfile
set SG =   `nawk '/Space group/{gsub("\047","");print $5}' $tempfile`
set SGnum = ` nawk '/Space group/{print $NF+0}' ${tempfile} `
set SG = ` nawk -F "[\047]" '/Space group/{print $2}' ${tempfile} `
set SG = ` nawk -v num=$SGnum '$1==num && NF>5{print $4}' ${CLIBD}/symop.lib `
set CELL = `nawk '/Cell Dimensions/{getline;getline;print}' $tempfile`
set mtzRES = `nawk '/Resolution Range/{getline;getline;print $6}' $tempfile`
rm -f $tempfile

set temp = `echo "$RES $mtzRES" | nawk '$1>=$2{print $1} $1<$2{print $2}'`
if("$temp" != "$RES") then
    echo "$rawfile resolution is only ${temp}A"
    set RES = "$temp"
endif

##############################################################################
# change of space group? 
#
##############################################################################
if(("$?new_SG")||($?REINDEX)) then
    # need a new space group
    if(! $?REINDEX) set REINDEX = ""
    
    # change the space group
    echo "reindexing $rawfile to $new_SG $REINDEX"
    reindex HKLIN $rawfile HKLOUT ${tempfile}.${new_SG}.mtz << end_reindex
SYMMETRY $new_SG
$REINDEX
END
end_reindex
    set SG = $new_SG
    set rawfile = ${tempfile}.${SG}.mtz
endif

goto Return_From_Setup