#! /bin/csh -f
#
#	Phaser's fft script for making the phased map
#
#
##############################################################################
# set up awk
alias nawk /usr/bin/nawk
nawk 'BEGIN{print 1; exit}' >& /dev/null
if($status) alias nawk awk

# defaults
set MAPMAN   = /programs/bin/mapman
set BRIX     = /usr/bin/brix
set BONES    = /usr/bin/bones
set pdbfile  = ""
set mtzfile  = mtz/best_phased.mtz
set pdbfile  = pdb/build1.pdb
set tempfile = fft_temp

# oversampling makes maps look prettier
set hiRES    = ""
set F        = Fin
set SIGF     = SIGFin
set PHI      = PHIDM
set FOM      = FOMDM

# output files
set ccp4map  = ffted.map
set omap     = ffted.omap
set pickpdb  = fftpick.pdb
set omacro   = ffted.omacro
set bonefile = bones.o

set pick = pick

if("$1" == "") goto Help
goto Setup
#  Procedure (at bottom) to read command-line args
#  mtz, Fs, or resolution
Help:
cat << EOF

usage: $0 mtzfile.mtz [F] [PHI] [coverme.pdb]
   or: $0 mapfile.map [coverme.pdb]

where:
mtzfile.mtz  - a (phased) mtz file you want a map from
mapfile.map  - a pre-calculated map you want to convert to o format
coverme.pdb  - a pdb file you want the output map to cover

EOF
exit 9
ReturnFromSetup:

if($?user_mapfile) then
    # jump ahead if user specified a map
    cp $user_mapfile ${tempfile}.map
    goto normalize
endif


################################################################################
fft HKLIN $mtzfile MAPOUT ${tempfile}.map << EOF-fft
RESOLUTION 1000 $hiRES
#EXCLUDE SIG1 1
title  $hiRES A map of $FOM * $F @ $PHI
LABIN F1=$F $SIG1 PHI=$PHI $WFOM
END
EOF-fft
if($status) exit

normalize:
# normalize the map
echo "SCALE SIGMA" |\
mapmask mapin ${tempfile}.map mapout $ccp4map
rm -f ${tempfile}.map >& /dev/null

extend:
# extend the map

# quit if user wasn't interested in an O map
if("$omap" == "") exit

# get space group/CELL from map header
set CELL = `echo "go" | mapdump mapin $ccp4map | nawk '/Cell dimensions/{print $(NF-5), $(NF-4), $(NF-3), $(NF-2), $(NF-1), $NF; exit}'`
set SGnum = `echo "go" | mapdump mapin $ccp4map | nawk '/Space-group/{print $NF}'`
set SG = `nawk -v SGnum=$SGnum '$1 == SGnum {print $4;exit}' $CLIBD/symop.lib`

# decide on how to extend the map (cover or fill cell)
if(-e "$pdbfile") then
    set pickpdb = "$pdbfile"
    set pick = "build"
    set center = `echo "COM" | pdbset xyzin $pdbfile XYZOUT ${tempfile}.pdb | nawk '$1=="Center" && $3=="Mass:"{print $4, $5, $6}'`
    rm -f ${tempfile}.pdb >& /dev/null
else
    set center  = `echo $CELL | nawk '{print $1/2, $2/2, $3/2}'`
endif


# make an O macro for looking at the results
set temp = `dirname $pickpdb`
set temp = `cd $temp ; pwd`
set pickpdb = ${temp}/`basename $pickpdb`
cat << EOF >! $omacro

! read in a pdb file
sam_atom_in $pickpdb $pick
mol $pick
obj $pick
zone ;
end

sym_set ; ; $SG
sym_cell
! cen_xyz $center

! read in the sites
sam_atom_in `pwd`/o/sites.pdb sites
mol sites
obj sites
zone ;
end
sym_set ; ; $SG
sym_cell

! display them as big spheres
sketch_cpk sites
sym_sphere sites sym 30
sketch_cpk sym1
clear_flags
sketch_cpk sym2
clear_flags
sketch_cpk sym3
clear_flags
sketch_cpk sym4
clear_flags
sketch_cpk sym5
clear_flags
sketch_cpk sym6
clear_flags
sketch_cpk sym7
clear_flags
sketch_cpk sym8
clear_flags

read `pwd`/$bonefile
bone_setup skel bones 30 1 2 3 4 5
bone_draw

map_cache
map_active_center
map_file `pwd`/$omap
map_object map
!         dx dy dz sig color linestyle
map_param 25 25 25 1   white 0.5 0 1
map_draw

menu @$omacro on

EOF

# now actually do the map extension
if("$pick" == "build") then
    echo "border 10" |\
    mapmask mapin $ccp4map xyzin $pdbfile mapout ${tempfile}.map
else
    # just extend to unit cell (and then some)
    echo "XYZLIM -0.1 1.1   -0.1 1.1   -0.1 1.1 " |\
    mapmask mapin $ccp4map mapout ${tempfile}.map

    # pick peaks in it (just to have something to grab onto)
    peakmax MAPIN ${tempfile}.map XYZOUT $pickpdb << eof-pick
    THRESHOLD RMS 2
    OUTPUT BROOKHAVEN
    END
eof-pick
endif


mapman:
# convert to O format (and do a bones trace)
setenv MAPSIZE `ls -ln ${tempfile}.map | nawk '{printf "%d", $5/3.9}'`
$MAPMAN << end-mapman |& grep -v Toodle
read map1 ${tempfile}.map CCP4
mappage map1 $omap
bone skel map1 1.5 0.5 100
bone conn $bonefile skel 5
quit
y
end-mapman
if(! $status) then
    # no need to do further conversions
    set BONES = ""
    set BRIX  = ""
endif

if(-e "$BRIX") then
    # use the brix program to make an o-readable file
    $BRIX ${tempfile}.map $omap
endif

if((-e "$BONES")&&(! -e "$pdbfile")) then
    $BONES ${tempfile}.map << EOF
    1.5 0.5
    5
$bonefile
skel
EOF
endif

rm -f ${tempfile}.map >& /dev/null



####################################################################
exit
####################################################################



Setup:
# scan the command line for files
foreach arg ( $* )
    if( "$arg" =~ *.mtz ) then
	if(! -e "$arg") then
	    echo "WARNING: $arg does not exist! "
	    continue
	endif
	set mtzfile  = "$arg"
	continue
    endif
    if(("$arg" =~ *.map)||("$arg" =~ *.ext)) then
	if(! -e "$arg") then
	    echo "WARNING: $arg does not exist! "
	    continue
	endif
	set user_mapfile  = "$arg"
	continue
    endif
    if(("$arg" =~ *.pdb)||("$arg" =~ *.brk)) then
	if(! -e "$arg") then
	    echo "WARNING: $arg does not exist! "
	    continue
	endif
	set pdbfile  = "$arg"
	continue
    endif
    if( "$arg" =~ [0-9]* ) then
	set temp = `echo "$arg" | nawk '$1+0>0.1{print $1+0}'`
	if("$temp" != "") set hiRES = "$temp"
    endif
end

# return early if a map was specified (not actually going to do an fft! )
if($?user_mapfile) goto ReturnFromSetup

#get variables from mtz file
echo "go" | mtzdump hklin $mtzfile |\
nawk '/OVERALL FILE STATISTICS/,/No. of reflections used/' |\
nawk 'NF>10 && $(NF-1) ~ /[FQPWADI]/' |\
cat >! ${tempfile}mtzdmp

# use completeness, or F/sigF to pick default F
cat ${tempfile}mtzdmp |\
nawk '$(NF-1) == "F"{F=$NF; meanF=$8; reso=$(NF-2); comp=substr($0,32)+0; \
      getline; S=$NF; if($8) meanF /= $8; print F, S, reso, comp, meanF;}' |\
sort +2n -3 +3nr -4 +4nr >! ${tempfile}F

# and extract all dataset types/labels
cat ${tempfile}mtzdmp |\
nawk 'NF>2{print $(NF-1), $NF, " "}' |\
cat >! ${tempfile}cards

#clean up
rm -f ${tempfile}mtzdmp

# pick F with best resolution, or <F>/<sigma(F)>
set F    = `head -1 ${tempfile}F`
if($#F > 2) then
    set SIGF = $F[2]
    set F    = $F[1]
endif

# pick most recent phase/FOM
grep "P $PHI" ${tempfile}cards >& /dev/null
if($status) then
    set temp = `nawk '/^P/{print $2}' ${tempfile}cards  | tail -1`
    if("$temp" != "") set PHI = "$temp"
endif
grep "W $FOM" ${tempfile}cards >& /dev/null
if($status) then
    set temp = `nawk '/^W/{print $2}' ${tempfile}cards | tail -1`
    if("$temp" != "") then
	set FOM = "$temp"
    else
	# there are no FOMs in this mtz file
	set FOM = ""
    endif
endif

# see if user specified an F, Phase, or FOM
foreach arg ( $* )
    set temp = `grep " $arg " ${tempfile}cards`
    if("$temp" =~ F*) then
	set F = "${arg}"
	set temp = `nawk -v arg="$arg" '$1==arg{print $2}' ${tempfile}F`
	if($#temp == 1) set SIGF = "$temp"
	continue
    endif
    if("$temp" =~ Q*) set SIGF = "${arg}"
    if("$temp" =~ P*) set PHI  = "${arg}"
    if("$temp" =~ W*) set FOM  = "${arg}"
    
    if($?NO && ("$arg" == FOM)) set FOM = ""
    
    unset NO
    if("$arg" == "no") set NO
end

# now check and see if the sigma is really there
grep "Q $SIGF" ${tempfile}cards >& /dev/null
if($status) then
    # no sigma availale (doesn't really matter anyway)
    set SIGF = ""
endif

rm -f ${tempfile}cards ${tempfile}F >& /dev/null


# assign the actual fft cards here (and blank them if they don't exist)
set SIG1 = "SIG1=$SIGF"
set WFOM = "W=$FOM"

if("$SIGF" == "") set SIG1 = ""
if("$FOM"  == "") then
    set FOM  = 1
    set WFOM = ""
endif

if(("$hiRES" == "")&&(-e "$mtzfile")&&(! $?user_mapfile)) then
    set hiRES = `echo "head" | mtzdump hklin $mtzfile | nawk '/Resolution Range/{getline;getline;print $6}' `
    # zero-fill for better-looking map
    set hiRES = `echo $hiRES |  nawk '$1>0{print 1/(2.0*((1/$1)^3))^(1/3)}'`
endif

goto ReturnFromSetup