#! /bin/tcsh -f
#
#	convert a CCP4 map into a PDB file with "rho" encoded as the B factor
#
#
set mapfile = "$1"

set tempfile = ${CCP4_SCR}/map2pdb$$


mapmask mapin $mapfile mapout ${tempfile}dumpme.map << EOF | tee ${tempfile}mapdump.log | grep dens
axis X Y Z
EOF
cat ${tempfile}mapdump.log |\
awk '/Grid sampling on x, y, z/{gx=$8;gy=$9;gz=$10}\
     /Maximum density /{max=$NF}\
     /Cell dimensions /{xs=$4/gx;ys=$5/gy;zs=$6/gz}\
     /Number of columns, rows, sections/{nc=$7;nr=$8;ns=$9}\
 END{print xs,ys,zs,nc,nr,ns,max}' >! ${tempfile}mapstuff.txt

awk '{print "grid points every:",$1,$2,$3}' ${tempfile}mapstuff.txt


set CELL = `awk '/Cell dimensions /{print $4,$5,$6,$7,$8,$9;exit}' ${tempfile}mapdump.log`
set SG   = `awk '/Space-group /{print $NF;exit}' ${tempfile}mapdump.log`


set size = `awk '{print 4*$4*$5*$6}' ${tempfile}mapstuff.txt`
set head = `ls -l ${tempfile}dumpme.map | awk -v size=$size '{print $5-size}'`
od -vf -w4 -j $head ${tempfile}dumpme.map | awk 'NF==2{print $2}'|\
cat ${tempfile}mapstuff.txt - |\
awk 'NR==1{xs=$1;ys=$2;zs=$3;nc=$4;nr=$5;ns=$6;max=$7;next} {print xs*c,ys*r,zs*s,$1,100*$1/max}\
  {++c} c>=nc{c=0;++r} r>=nr{r=0;++s}' |\
tee mapdump.txt |\
awk '{X=$1;Y=$2;Z=$3;rho=$4;nrho=$5;\
  printf("ATOM  %5d  E   RHO  %4d     %7.3f %7.3f %7.3f %5.2f%6.2f  \n",(++Atomnum)%100000,1,X,Y,Z,rho,nrho)}' |\
cat >! temp.pdb

pdbset xyzin temp.pdb xyzout mapdump.pdb << EOF >! ${tempfile}pdbset.log
CELL $CELL
SYMM $SG
EOF
if(! -s mapdump.pdb) exit 9

rm -f ${tempfile}dumpme.map ${tempfile}mapstuff.txt ${tempfile}mapdump.log ${tempfile}pdbset.log