#! /bin/csh -f
#
#	script for converting mosflm spot output with into the
#	reciprocal-space x,y,z coordinates used by DPS
#
########
# 
# use the most recently-created .spt file by default
set spotfile  = `ls -1rt *.spt |& tail -1`
set tempfile  = spt2pdb_temp

# camera settings (not available in *.spt file)
set wavelength  = ""
set distance    = ""
set twotheta    = 0
set beam_center = ""
set hiRES       = 1
set loRES       = 1000

goto Setup
Help:
############################################################
cat << EOF
usage: spt2xyz.com mosflm.spt 1.0A 100mm P212121
  or
usage: spt2xyz.com mosflm.spt mosflm.com P212121

where:
mosflm.spt  - the "spot" file (picked with mosflm)
1.0A        - the x-ray wavelength
100mm       - the xtf distance
mosflm.com  - a mosflm script containing CELL, DIST, WAVE, etc., 
3-20A	    - resolution range to use (optional)
50sig	    - I/sigma cutoff to use (default: 0)

EOF
exit 9
Return_from_Setup:

if("$wavelength" == "") goto Help
if("$distance" == "") goto Help

############################################################
# retrieve useful parameters from the mosflm spot file
echo "$0 $*"
echo "using ${spotfile} ($frames images)"
if($#beam_center != 2) set beam_center = `nawk 'NR==3 && NF==2{print;exit}' $spotfile`

set rotation    = `nawk 'NR==1{rot=$5} NR==2{print (360-rot+180*$1)%360;exit}' $spotfile`
set ysign       = `nawk 'NR==2{print 1-($1*2);exit}' $spotfile`
set raster      = `nawk '{print $3;exit}' $spotfile`
set yscale      = `nawk '{print $4;exit}' $spotfile`
set beamstop    = `echo $beam_center $raster $yscale | nawk '{printf "r %d s %d", $1/$3,$2*$4/$3}'`
set sscale      = `echo $ysign $yscale | nawk '{print $1/$2}'`
set mode        = automatic

#set ysign = -1
#set rotation = 180

# report values input into spt2xyz
echo "wavelength:  $wavelength"
echo "beam center: $beam_center"
echo "beam center: $beamstop"
echo "distance:    $distance"
echo "2theta:      $twotheta"
echo "sscale:      $sscale"
echo "rotation:    $rotation"
echo "resolution:  $loRES $hiRES"
echo "sigma cut:   $cutoff"
echo ""

cat << EOF >! ${tempfile}info
DIST $distance
WAVE $wavelength
BEAM $beam_center
OMEG $rotation
SIGN $ysign
RESO $loRES $hiRES
EOF

# create a temporary awk program for doing the conversion...
cat << EOF-awk >! ${tempfile}spt2xyz.awk
#! /usr/bin/awk -f
#
#   Awk program for converting a mosflm-formatted *.spt file
#   into reciprocal-space coordinates
#
#   the conventions of omega and sign are analogous to the 
#   DPS "rotation" and "yscale" parameters
#

BEGIN{
    D          = $distance
    lambda     = $wavelength
    Xbeam      = $beam_center[1]
    Ybeam      = $beam_center[2]
    OMEGA      = $rotation
    SIGN       = $ysign
    cutoff     = $cutoff
    hires      = $hiRES
    lores      = $loRES
    
    PI=2*atan2(1,0)
}

# "key-worded" input can be concatenated in front of the *.spt file
/^DIST/{D=\$2}
/^WAVE/{lambda=\$2}
/^BEAM/{Xbeam=\$2;Ybeam=\$3}
/^OMEG/{OMEGA=\$2*PI/180}
/^SIGN/{SIGN=\$2}
/^SIGMA/{cutoff=\$2}
/^RESO/{lores=\$2;hires=\$3; if(lores<hires){hires=\$2;lores=\$3}}

NF>5 && \$1!="-99.00" && \$1!="-999.00" && \$6+0>0{

    # read in a spot from the mosflm *.spt format
    Xo=\$1-Xbeam;
    Yo=\$2-Ybeam;
    phi1=\$4-\$3;
    phi2=\$4+\$3;
    height=\$5;
    sigma=\$6;
    
    # now transform into some kind of standard system:
    # Z == x-ray beam  Y == spindle
    
    # make the phi spin clockwise (phi1=start;phi2=end)
    phi1 = SIGN*PI/180*phi1;
    phi2 = SIGN*PI/180*phi2;
    
    # transform detector coordinates so that y is phi axis
    X  = Xo*cos(omega)-Yo*sin(omega);
    Y  = Xo*sin(omega)+Yo*cos(omega);
    
    # now map each spot back onto the surface of the Ewald sphere
    # at the moment, the center of the Ewald sphere is the origin
    distortion = lambda*sqrt((X)^2+(Y)^2+(D)^2)
    xo = X/distortion
    yo = Y/distortion
    zo = D/distortion
    
    # shift the spots back onto the true origin of reciprocal space
    # (the surface of the Ewald sphere)
    zo -= 1/lambda;
    
    # apply the phi rotation about "Y"
    x1 = xo*cos(phi1)-zo*sin(phi1)
    y1 = yo
    z1 = xo*sin(phi1)+zo*cos(phi1)
    x2 = xo*cos(phi2)-zo*sin(phi2)
    y2 = yo
    z2 = xo*sin(phi2)+zo*cos(phi2)

    # apply sigma cutoff
    if(height/sigma<cutoff) next
    # apply resolution cutoff
    d = 1/sqrt(x1*x1+y1*y1+z1*z1)
    if(d<hires||d>lores)next

    # count up number of spots...
    ++n;
    printf "ATOM  %5d  O   WAT  %4d    %8.3f%8.3f%8.3f %5.2f%6.2f\\n",
           n,n/10,x1*D,y1*D,z1*D,d,phi1;
    printf "ATOM  %5d  O   WAT  %4d    %8.3f%8.3f%8.3f %5.2f%6.2f\\n",
           n,n/10,x2*D,y2*D,z2*D,d,phi1;
}
EOF-awk


# convert the mosflm spot file (X,Y,phi) into reciprocal (x,y,z) coordinates
cat ${tempfile}info $spotfile |\
nawk -f ${tempfile}spt2xyz.awk |\
nawk 'length($0)==66' |\
cat >! peaks.pdb

cat << EOF
sam_atom_in peaks.pdb peaks
mol peaks
paint_ramp atom_wt ; red blue ;
paint_ramp atom_b ; red blue ;
obj peaks
zone ;
end

sym_set ; 1000 1000 1000 90 90 90 p-1
sym_cell
sym_obj peaks sym 10
EOF


# clean up
rm -f ${tempfile}info >& /dev/null
rm -f ${tempfile}spt2xyz.awk >& /dev/null
exit
############################################################


Setup:

nawk 'BEGIN{exit}' >& /dev/null
if($status) alias nawk awk
set cutoff = 0

# scan the command line for files
set i = 0
while( $i < $#argv )
    @ i = ( $i + 1 )
    @ nexti = ( $i + 1 )
    @ lasti = ( $i - 1 )
    if($nexti > $#argv) set nexti = $#argv
    if($lasti < 1) set lasti = 1
    set arg = "$argv[$i]"
    set ARG = `echo "$arg" | nawk '{print toupper($1)}'`
    
    if("$arg" =~ *.spt) then
	if(! -e "$arg") then
	    echo "WARNING: $arg does not exist! "
	    continue
	endif
	set spotfile = "$arg"
	
	continue
    endif
    
    if(-e "$arg") then
	# assume this is the mosflm.com file
	set temp = `nawk '/^WAVE/ && $2+0>0.1 && $2+0<10{print $2}' $arg`
	if("$temp" != "") set wavelength = "$temp"
	set temp = `nawk '/^DIST/ && $2+0>10 && $2+0<1000{print $2}' $arg`
	if("$temp" != "") set distance = "$temp"
	set temp = `nawk '/^TWOTHETA/ && $2+0>-90 && $2+0<90{print $2}' $arg`
	if("$temp" != "") set twotheta = "$temp"

	# get upper resolution limit
	set temp = `nawk '$1 ~ /^RESO/{for(i=2;i<=NF;++i){print $i}}' $arg | nawk '$1+0>0.3 && $1+0 < 15' | tail -1`
	if($#temp == 1) then
	    if("$temp" =~ [0-9]*) set hiRES  = "$temp"
	endif
	continue
    endif
    
    # now read things that aren't files

    if("$arg" =~ [0-9]*) then
	# we have a number
	if(("$arg" =~ *A)||("$argv[$nexti]" == "A")) then
	    # user-preferred resolution limits
	    set temp = `echo "$arg" | nawk 'BEGIN{FS="-"} $1+0 > 0.1{print $1+0} $2+0 > 0.1{print $2+0}'`
	    if($#temp != 1) then
		set temp = `echo $temp | nawk '$1>$2{print $1, $2} $2>$1{print $2, $1}'`
		if($#temp == 2) then
		    # a range 3-10A is a resolution
		    set hiRES = `echo "$temp[2]" | nawk '{print $1+0}'`
		    set loRES = `echo "$temp[1]" | nawk '{print $1+0}'`
		    continue
		endif
	    else
		# single 1.1A is a wavelength
		set wavelength = `echo "$arg" | nawk '{print $1+0}'`
		continue
	    endif
	endif
	
	if("$arg" =~ *x*) then
	    # user-preferred beam center
	    set temp = `echo "$arg" | nawk 'BEGIN{FS="x"} {print $1+0, $2+0}'`
	    if($#temp == 2) then
		# a cross 95x94mm is a beam center
		set beam_center = ( $temp )
		continue
	    else
		# means nothing? 
		continue
	    endif
	endif
	
	if(("$arg" =~ *mm)||("$argv[$nexti]" == "mm")) then
	    # distance is expressed in mm
	    set distance = `echo "$arg" | nawk '{print $1+0}'`
	    continue
	endif
	
	if(("$arg" =~ *sig)||("$argv[$nexti]" == "sig")) then
	    # sigma cutoff for spot use
	    set cutoff = `echo "$arg" | nawk '{print $1+0}'`
	    continue
	endif
    endif
    
    
    unset NO
    if("$ARG" == NO) set NO
    if("$ARG" == NOT) set NO
end

# get the number of images that were peak-picked into this spotfile
set frames = `nawk '$NF=="-99.0" || $NF=="-999.0"' $spotfile | wc -l`
set spots = `nawk 'NF>5 && $1!="-99.00" && $1!="-999.00"' $spotfile | wc -l`

goto Return_from_Setup