#! /bin/csh -f
#
#
#      Calculate an O map of a given PDB rendered at the provided resolution
#
#
#
set pdbin       = RNase_H.pdb
set SG		= P212121
set MAPMAN      = /programs/o/rave/rave_irix6/6d_mapman
set center      = center.pdb  	# center of view in o
set sfall_reso  = 0.6		# resolution to use in sfall (keep low to prevent aliasing)
set grid_reso   = 0.8		# resolution to use for grid spacing
set reso        = 1.0		
set WilsonB     = 5
set logfile     = rolloff.log
set outfile     = this.omap

# scan command line for user-specified pdb and resolution
foreach arg ( $* )
    if("$arg" =~ *.pdb) set pdbin = $arg
    if("$arg" =~ [0-9]*) set reso = $arg
end
if(! -e "$center") set center = "$pdbin"

# use the B-factor to roll-off resolution
set B = `echo $reso | nawk '{print 78.956835*($1/3)^2}'`
set WilsonB = `nawk '/^ATOM/{B += substr($0, 61, 6); ++n} END{if(n) print B/n}' $pdbin`
set Badd = `echo $B $WilsonB | nawk '{print $1 - $2}'`

# calculate Fs to $render_reso A
echo "calculating Fs to ${render_reso}A"
sfall xyzin $pdbin HKLOUT temp.mtz << EOF >> $logfile
MODE SFCALC XYZIN
SYMM $SG
RESOLUTION 1000 $render_reso
EOF

# make the map (at constant grid spacing)
echo "calculating map with Badd = $Badd on ${grid_reso}A grid"
fft hklin temp.mtz mapout fft.map << EOF >! $logfile
RESOLUTION $grid_reso
LABIN F1=FC PHI=PHIC
SCALE F1 1 $Badd
EOF
# normalize it
mapmask mapin fft.map mapout norm.map << EOF >> $logfile
SCALE SIGMA
EOF
# cover the central point-of-interest
mapmask mapin norm.map xyzin $center mapout ccp4.map << EOF >> $logfile
BORDER 12
EOF

# convert to O format
rm -f $outfile >& /dev/null
setenv MAPSIZE `ls -l ccp4.map | nawk '{print $5/3.9}'`
$MAPMAN << EOF   >& $logfile
read map ccp4.map CCP4
mappage map $outfile
quit
y
EOF
if(-e $outfile) echo "done"
if(! -e $outfile) echo "crap! "

# clean up
rm -f temp.pdb temp.mtz
rm -f fft.map norm.map ccp4.map