#! /bin/csh -f
#
#	reindex.com				-James Holton  2-22-05
#
#	simple interface to reindexing of merged or unmerged mtz data
#
#
# set this to wherever your awk program is
set nawk = nawk
$nawk 'BEGIN{print}' >& /dev/null
if($status) set nawk = awk
alias nawk $nawk

set mtzfile    = ""
set logfile    = "reindex.log"
set newmtzfile = "reindexed.mtz"
set tempfile   = ./reindex_temp

set message    = ""
set newSG      = ""

goto Setup
#
#	scan command line for mtz files, space groups, and preferences
#
Help:
cat << EOF-Help
usage: $0 mtzfile.mtz [SG] [flip] [newmtz.mtz]

where:

mtzfile.mtz	- the mtz data you want to re-index. (merged or unmerged)
SG		- the new space group (P43 for example)
flip		- flip ambiguous axes of P4x P3x and P6x cells
newmtz.mtz	- name you want for the reindexed mtz (default: reindexed.mtz)

Example:    $0 Fpp.mtz P43212 Fpp.P43212.mtz
Example:    $0 Fpp.mtz flip

Note: to put the screw axis of P2221 on the shortest cell axis, say P2122
The space group in the new file will be P2221 again, but the cell
will be flipped around appropriately.  Same goes for P2212, or
P21221 and P22121 with P21212

EOF-Help
exit
#
Return_from_Setup:

if(! -e "$mtzfile") goto Help

# okay to proceed
if("$message" != "") set message = " $message"
echo "reindexing $mtzfile to $newmtzfile (in ${newSG}${message})" | tee $logfile

if($?MERGED) then
    # apply any "flip" directives with "reindex"
    echo "$REINDEX" |\
    reindex HKLIN $mtzfile HKLOUT ${tempfile}reindexed.mtz >> $logfile
    if($status) set BAD
    
    # just use mtzutils to change SG record in header
    echo "SYMM $newSGnum" |\
    mtzutils HKLIN ${tempfile}reindexed.mtz HKLOUT ${tempfile}newSG.mtz >> $logfile
    if($status) set BAD
    
    rm -f ${tempfile}reindexed.mtz >& /dev/null
    
    # sort it again, just to be safe
    echo "H K L" |\
    sortmtz HKLIN ${tempfile}newSG.mtz HKLOUT $newmtzfile >> $logfile
    if($status) set BAD
    
    rm -f ${tempfile}newSG.mtz >& /dev/null
endif
    
if($?UNMERGED) then
    reindex HKLIN $mtzfile HKLOUT ${tempfile}reindexed.mtz << EOF-reindex >> $logfile
    SYMM $newSGnum
    $REINDEX
EOF-reindex
    if($status) set BAD

    # sort again
    echo "H K L M/ISYM BATCH I SIGI" |\
    sortmtz HKLIN ${tempfile}reindexed.mtz HKLOUT $newmtzfile >> $logfile
    if($status) set BAD
    
endif
# remove intermediate file
rm -f ${tempfile}reindexed.mtz >& /dev/null


if((-e "$newmtzfile")&&(! $?BAD)) then
    echo "$newmtzfile is $mtzfile in ${newSG}${message}"
    if("$axes" != "") then
	echo "head" | mtzdump hklin "$newmtzfile" |\
	nawk '/Cell Dimensions/{getline;getline;\
	print "new cell:", $1, $2, $3, $4, $5, $6}'
    endif
else
    # something went wrong
    echo "FAILED!  examine $logfile to see what went wrong."
endif

exit







Setup:
#############################################################
# check up on essentials
if(! $?CLIBD) then
    echo "Please set up CCP4, and then run $0 again"
    goto Help
endif
if(! -e $CLIBD/symop.lib) then
    echo "ERROR: no $CLIBD/symop.lib"
    echo "Please set up CCP4, and then run $0 again"
    goto Help
endif

set REINDEX
set axes
set FLIP = 0

#first, get filenames from the command line
foreach arg ( $* )
    if("$arg" =~ *.mtz) then
	if("$mtzfile" == "") then
	    # havn't initialized this yet
	    if(-e "$arg") then
		# check to make sure it is readable
		echo "HEAD" | mtzdump HKLIN $arg >&! ${tempfile}mtzdump
		grep "Space group" ${tempfile}mtzdump >& /dev/null
		if($status) then
		    echo "WARNING: $arg is not an MTZ file! "
		    continue
		endif
		set mtzfile = "$arg"
	    else
		echo "WARNING: $arg does not exist! "
		continue
	    endif
	else
	    # already have an input mtz, so assume this is output
	    set newmtzfile = "$arg"
	    
	    # look for "hidden" space group in name? 
	    echo $newmtzfile |\
	     nawk '{for(i=1;i<=length($0);++i){c=substr($0,i,1);\
	       if(c ~ /[PpCcIiFfRrHh1-6]/){printf c}else{print ""}}}' |\
	     nawk '$1~/^[PpCcIiFfRrHh][1-6]/' >! ${tempfile}SGs
	    foreach sg ( `tail -10 ${tempfile}SGs` )
		set temp = `nawk -v SG=$sg '$4 == SG && $1 < 500 {print $4}' $CLIBD/symop.lib | head -1`
		if("$temp" != "") then
		    set newSG = "$temp"
		endif
	    end
	    rm -f ${tempfile}SGs >& /dev/null
	endif
    endif

    # check for new space group
    if("$arg" =~ [PpCcIiFfRrHh][1-6]*) then
	# check for SGs listed in library (but not the screwy ones)
	set temp = `echo $arg | nawk '{print toupper($1)}'`
	if($?CLIBD) then
	    set temp = `nawk -v SG=$temp '$4 == SG && $1 < 500 {print $4}' $CLIBD/symop.lib | head -1`
	endif
	if("$temp" != "") then
	    set newSG = "$temp"
	endif
	
	set temp = `echo $arg | nawk '{print toupper($1)}'`
	# check for orthorhombic "pseudo-spacegroup" language
	if("$temp" == P2221) then
	    # P2221 with screw along longest axis
	    set axes  = "a b c"
	    set newSG = "P2221"
	    continue
	endif
	if("$temp" == P2212) then
	    # P2221 with screw along mid-length axis
	    set axes  = "b c a"
	    set newSG = "P2221"
	    continue
	endif
	if("$temp" == P2122) then
	    # P2221 with screw along shortest axis
	    set axes  = "c a b"
	    set newSG = "P2221"
	    continue
	endif
	
	if("$temp" == P21212) then
	    # P21212 with non-screw along longest axis
	    set axes  = "a b c"
	    set newSG = "P21212"
	    continue
	endif
	if("$temp" == P21221) then
	    # P21212 with non-screw along mid-length axis
	    set axes  = "b c a"
	    set newSG = "P21212"
	    continue
	endif
	if("$temp" == P22121) then
	    # P21212 with non-screw along shotest axis
	    set axes  = "c a b"
	    set newSG = "P21212"
	    continue
	endif
	    
	if("$temp" == P112) then
	    # P2 with twofold along longest axis
	    #set axes  = "a b c"
	    set newSG = "P2"
	    continue
	endif
	if("$temp" == P121) then
	    # P2 with twofold along mid-length axis
	    set axes  = "b c a"
	    set newSG = "P2"
	    continue
	endif
	if("$temp" == P211) then
	    # P2 with twofold along shortest axis
	    set axes  = "c a b"
	    set newSG = "P2"
	    continue
	endif
	    
	if("$temp" == P1121) then
	    # P21 with twofold along longest axis
	    #set axes  = "a b c"
	    set newSG = "P21"
	    continue
	endif
	if("$temp" == P1211) then
	    # P21 with twofold along mid-length axis
	    set axes  = "b c a"
	    set newSG = "P21"
	    continue
	endif
	if("$temp" == P2111) then
	    # P21 with twofold along shortest axis
	    set axes  = "c a b"
	    set newSG = "P21"
	    continue
	endif
    endif
	
    # check for flipping of ambiguous axes
    if("$arg" == "flip") then
	# user requested "flip" of axes
	@ FLIP = ( $FLIP + 1 )
    endif
    
    # what about cubic? 
end

if(! -e "$mtzfile") then
    goto Help
endif

# get info from the MTZ (should be left over from command-line check)
set CELL  = `nawk '/Cell Dimensions/{getline;getline;print}' ${tempfile}mtzdump`
set SGnum = ` nawk '/Space group/{print $NF+0}' ${tempfile}mtzdump `
cat $CLIBD/symop.lib |\
nawk -v SGnum=$SGnum '$1==SGnum{type=substr(tolower($6),1,1); key=substr($6,4,1);\
    if((type=="t")&&(key=="G")) type="h"; if((type=="t")&&(key=="C")) type="a";\
    print $1, $5, type substr($4,1,1), $6,$4}' |\
head -1 >! ${tempfile}sgdata
set SGnum = `nawk '{print $1}' ${tempfile}sgdata`
set SG    = `nawk '{print $5}' ${tempfile}sgdata`
set PG    = `nawk '{print $2}' ${tempfile}sgdata`
set BRAV  = `nawk '{print $3}' ${tempfile}sgdata`
set LATT  = `nawk '{print $4}' ${tempfile}sgdata`

if("$newSG" == "") set newSG = "$SG"
cat $CLIBD/symop.lib |\
nawk -v SG=$newSG '$4==SG{type=substr(tolower($6),1,1); key=substr($6,4,1);\
    if((type=="t")&&(key=="G")) type="h"; if((type=="t")&&(key=="C")) type="a";\
    print $1, $5, type substr($4,1,1), $6;}' |\
head -1 >! ${tempfile}sgdata
set newSGnum = `nawk '{print $1}' ${tempfile}sgdata`
set newPG    = `nawk '{print $2}' ${tempfile}sgdata`
set newBRAV  = `nawk '{print $3}' ${tempfile}sgdata`
set newLATT  = `nawk '{print $4}' ${tempfile}sgdata`
rm -f ${tempfile}sgdata >& /dev/null

# decide how/if to flip
if($FLIP) then
    if(("$SG" =~ [IP][46]*)||("$SG" =~ P3*12)||("$SG" =~ [PIF]2*3)) then
	set REINDEX = "reindex k, h, -l"
	set message = "with a and b axes flipped"
    endif
    if(("$SG" =~ [RH]32)||("$SG" =~ P3*21)) then
	set REINDEX = "reindex -h, -k, l"
	set message = "with a and b axes inverted"
    endif
    if(("$SG" =~ [PRH]3)||("$SG" =~ P3[12])) then
	# four possibilities here...
	if($FLIP == 1) then
	    set REINDEX = "reindex k, h, -l"
	    set message = "with a and b axes flipped"
	endif
	if($FLIP == 2) then
	    set REINDEX = "reindex -h, -k, l"
	    set message = "with a and b axes inverted"
	endif
	if($FLIP >  2) then
	    set REINDEX = "reindex -k, -h, -l"
	    set message = "with a and b axes flipped and inverted"
	endif
    endif
    if("$SG" == P1) then
	# permute the axes (from their cannonical order)
	if($FLIP == 1) set axes = "a b c"
	if($FLIP == 2) set axes = "b c a"
	if($FLIP == 3) set axes = "c a b"
    endif
    if("$PG" == PG2) then
	# flip the a and c axes?
	set REINDEX = "reindex l, -k, h"
	set message = "with a and c axes flipped"
    endif
    # for other space groups, "flip" is simply ignored
endif

if(("$newSG" == "")&&(! $FLIP)) then
    # WTF?
    echo "nothing to do! "
    rm -f ${tempfile}mtzdump >& /dev/null
    goto Help
endif

# now that whole command-line has been read, make some holistic decisions


# check to see if it's merged or unmerged
grep "Number of Batches" ${tempfile}mtzdump >& /dev/null
if($status) then
    set MERGED
else
    set UNMERGED
endif
rm -f ${tempfile}mtzdump


# decide on new axis ordering (for asymmetric orthorhombics)
if("$newSG" == "P222")    set axes = "a b c"
if("$newSG" == "P212121") set axes = "a b c"
if("$axes" != "") then
    # get current axis ordering
    # find out what the cannonical one would be
    # then decide how to go from current ordering to the desired one
    echo "$CELL" | nawk '{\
	# print out current axis order \
	print $1, "h"; print $2, "k"; print $3, "l"}' |\
    sort -n |\
    nawk '\
	# add cannonical axis names\
	NR==1{print $0, "a"} NR==2{print $0, "b"} NR==3{print $0, "c"}' |\
    nawk -v axes="$axes" 'BEGIN{split(axes, abc)} {\
	# write desired axis ordering in front of cannonical one \
	print abc[NR], $0}' |\
    sort |\
    nawk '# print out new hkl order \
          {printf "%s ", $3} END{print ""}' |\
    nawk '$1 $2 $3 $1 $2 $3 !~ /hkl/{$3 = "-" $3} {print $1, $2, $3}' |\
    cat >! ${tempfile}order    
    set REINDEX = "reindex "`nawk '{print $1 ",", $2 ",", $3}'  ${tempfile}order`

    # this should give us a mapping between any two orthorhombics
    
    set temp = `nawk '{print $NF}' ${tempfile}order`
    rm -f ${tempfile}order
    if("$newSG" == "P2221") then
	set message = "with screw along $temp axis"
    endif
    if("$newSG" == "P21212") then
	set message = "with non-screw along $temp axis"
    endif
    if("$newSG" == "P2") then
	set message = "with twofold along $temp axis"
    endif
    if("$newSG" == "P21") then
	set message = "with twofold screw along $temp axis"
    endif
    if("$newSG" == "P1") then
	set message = "with axes in $axes order."
    endif
endif


# known ways to convert between bravais lattices
if(("$newBRAV" == "oC")&&("$BRAV" =~ h[PR])) then
    # this lattice happens to fit, even if it's merged
    set REINDEX = "reindex h+k, k-h, l"
    set message = "with old hkl -> h+k,k-h,l"
    set OKAY
endif
set change_type = `echo $SG $newSG | nawk 'substr($1,1,1) != substr($2,1,1)'`
# cubic to just about anything should be okay
if(("$change_type" == "")&&("$BRAV" =~ c*)&&("$newBRAV" !~ h*)) set OKAY
# tetratonal to anything below it should be okay
if(("$change_type" == "")&&("$BRAV" =~ t*)&&("$newBRAV" !~ h*)&&("$newBRAV" !~ c*)) set OKAY
# orthorhombic to anything below it is okay

# anything can be converted to P1
if("$newBRAV" == aP) set OKAY

# now check for unhandlable situations


if(($?MERGED)&&("$newPG" != "$PG")&&("$PG" != "")&&(! $?OKAY)) then
    # can't do this
    echo "WARNING: You Shouldn't reindex $SG to $newSG for a merged mtz! "
    echo ""
    if($SGnum < $newSGnum) then
	cat << EOF
Since $newSG has higher symmetry than $SG, we would have
to merge some "unique" spots in $mtzfile together.  Although
this is technically possible, it is definitely not advisable.
If your space group really is $newSG, you will probably get
much better scaling and mergeing performance with the new
symmetry imposed during data processing, and probably in
data reduction as well.

EOF
    else
	cat << EOF
Since $SG has higher symmetry than $newSG, we would have
to "unmerge" the unique HKLs in $mtzfile to form multiple, 
yet completely identical Fs.  So, all this will do is
make your mtz bigger, so we will just change the spacegroup
name.

If you think you "overmerged" your data in $SG, then you
should go back to your mergeing step, and change $SG to
$newSG there.

EOF
    endif
#    exit
endif

if($?MERGED) goto Return_from_Setup
# now we are dealing with unmerged data

# what do we do about THIS! 
set temp = `echo "$newSG $SG" | nawk 'substr($1,1,1) != substr($2,1,1)'`
if("$temp" != "") then
    # crystal system has changed!
    if($newSGnum < $SGnum) then
	# we are going to loose some data
	
    else
	
    endif
    echo "WARNING: changing the lattice symmetry "
endif

# known ways to drop down to lower-symmetry bravais lattice

goto Return_from_Setup