#! /bin/csh -f
#
#	./scripts/maps.com 	- make refmac's maps
#
#
#	use this script to get the "current" map out of
#	refmacout.mtz, even if refmac is still running
#
alias nawk /usr/bin/nawk

set refmacmtz  = ./refmacout.mtz
set pdbin      = ./pdb/starthere.pdb	# PDB to "cover"
set xyzlim     = ""			# important for ARP maps

# need this to make O maps
set MAPMAN = /programs/o/rave/rave_irix6/6d_mapman

set tempfile   = ${CCP4_SCR}ffttemp

# user override defaults
foreach arg ( $* )
    if("$arg" =~ *.mtz) set refmacmtz = "$arg"
    if("$arg" =~ *.pdb) set pdbin     = "$arg"
    if("$arg" =~ *mapman*) set MAPMAN     = "$arg"
end

if(! -e "$refmacmtz") then
    echo "ERROR: no $refmacmtz ! "
    exit 9
endif

FFT:
###############################################################################
#
# make maps.
#
#  Sigmaa style 2mfo-dfc map with FREE R removed (for ARP)
#
echo "making arp's maps/2fofc_free.map"
fft hklin ./refmacout.mtz mapout ${tempfile}.map <<eof | tee ${tempfile}fft.log
title Sigmaa style 2mfo-dfc map calculated with refmac coefficients
labi F1=2FOFCWT PHI=PH2FOFCWT  FREE = FreeR_flag
end
eof
# take a stab at predicting arp's acceptable map
if("$xyzlim" == "") then
    # add one extra grid point for ARP
    set xyzlim = `nawk '/Map limits in grid points on xyz/{print "xyzlim", $(NF-5),$(NF-4)+1,$(NF-3),$(NF-2)+1,$(NF-1),$NF+1}' ${tempfile}fft.log | tail -1`
endif
rm -f ${tempfile}fft.log >& /dev/null
mapmask MAPIN ${tempfile}.map MAPOUT maps/2fofc_free.map << EOF
$xyzlim
END
EOF
#
#
#  Sigmaa style 2mfo-dfc map _without_ FREE R removed (for O)
#
echo "making maps/2fofc.map"
fft hklin $refmacmtz mapout maps/2fofc.map <<eof 
title Sigmaa style 2mfo-dfc map calculated with refmac coefficients
labi F1=2FOFCWT PHI=PH2FOFCWT
end
eof
#
#   Sigmaa style mfo-dfc map (FREE R already removed)
#
echo "making maps/fofc.map"
fft hklin $refmacmtz mapout ${tempfile}.map \
<<eof 
title Sigmaa style mfo-dfc map calculated with refmac coefficients
labi F1=FOFCWT PHI=PHFOFCWT  
end 
eof
mapmask MAPIN ${tempfile}.map MAPOUT maps/fofc.map << EOF 
$xyzlim
END
EOF
if(! $?debug) rm -f ${tempfile}.map



###############################################################################
# convert maps to O format
if(! -e ./o) goto done
if(! -e $MAPMAN) then
    echo "MAPMAN unavailable! "
    echo "cannot make O maps. "
    goto pick_peaks
endif

echo "making ./o/2fofc.omap covering $pdbin"
mapmask MAPIN maps/2fofc.map \
    XYZIN $pdbin \
    MAPOUT ${tempfile}.map << eof-extend 
SCALE SIGMA
BORDER 10
eof-extend
set mapsize = `ls -l ${tempfile}.map | awk '{print $5/3.9}'`
$MAPMAN -b mapsize $mapsize << eof-mapman |& cat
read map1 ${tempfile}.map CCP4
mappage map1 ./o/2fofc.omap
quit
eof-mapman

echo "making ./o/fofc.omap covering $pdbin"
mapmask MAPIN maps/fofc.map \
        XYZIN $pdbin \
        MAPOUT ${tempfile}.map << eof-extend 
SCALE SIGMA
BORDER 10
eof-extend
set mapsize = `ls -l ${tempfile}.map | awk '{print $5/3.9}'`
$MAPMAN -b mapsize $mapsize << eof-mapman |& cat
read map1 ${tempfile}.map CCP4
mappage map1 ./o/fofc.omap
quit
eof-mapman
rm -f ${tempfile}.map >& /dev/null
###############################################################################
pick_peaks:
# list difference peaks
echo "picking peaks in maps/fofc.map"
peakmax MAPIN maps/fofc.map XYZOUT ${tempfile}pick.pdb << eof-pick 
THRESHOLD RMS 4 NEGATIVES
OUTPUT BROOKHAVEN
END
eof-pick

# convert to an O macro
cat ${tempfile}pick.pdb |\
nawk '/^ATOM/{print substr($0, 61, 6)+0, substr($0, 31, 8), substr($0, 39, 8), substr($0, 47, 8)}' |\
sort -nr |\
nawk '{print "! ", $1, "sigma peak"; print "centre_xyz", $2,$3,$4; \
print "sym_sphere ; ; ;";print "@map";}' |\
cat >! ${tempfile}peaklist
rm -f ${tempfile}pick.pdb >& /dev/null

if($?NO_MACROS) goto done
# make some nice O macros for visiting large difference features
echo "generating o macros in ./o/"
cat ${tempfile}peaklist >! ./o/visit_peaks.omac

set temp = `awk '/^ATOM/{++n;X+=substr($0, 31, 8);Y+=substr($0, 39, 8);Z+=substr($0, 47, 8)} END{if(n) print "cen_x", X/n, Y/n, Z/n}' $pdbin`
cat << EOF-latest >! ./o/latest.omac
sam_atom_in `pwd`/$pdbin latest
mol latest
obj latest
zone ;
end
$temp

paint_ramp atom_b ; blue red
obj Bfac
zone ;
end

sym_set ; ; P212121 ; 
sym_cell

menu @visit_peaks.omac on
menu @map on
EOF-latest

if(! -e ./o/map) then
    cat << EOF-map >! ./o/map
map_cache
map_active_centre
map_file `pwd`/o/2fofc.omap
map_object sigmaa
map_parameter 15 15 15 1 tan 0.5 0 1
map_draw

map_cache
map_active_centre
map_file `pwd`/o/fofc.omap
map_object good
map_parameter 15 15 15 3 green 0.5 0 1
map_draw

map_cache
map_active_centre
map_file `pwd`/o/fofc.omap
map_object bad
map_parameter 15 15 15 -3 red 0.5 0 1
map_draw

EOF-map
endif

done:
# list the difference peaks
cat ${tempfile}peaklist |\
nawk '/sigma/{printf "%s | ", $0} /cen/{print}'
rm -f ${tempfile}peaklist >& /dev/null
echo "maps ready."
exit