#! /bin/csh -f
#
#	Jiffy for transforming spots on ADSC images into a PDB file
#	of spots in reciprocal space
#
#	needs the DPS program dps_peaksearch in the path.
#
alias nawk awk

# peak-picking stuff
set dps_peaksearch = peaks
#set dps_peaksearch = /programs/dps_V1.34/bin/sgi/dps_peaksearch
set threshold = 2.2
set hires     = 3
set lores     = 10000
set zoom      = 1

# normally read from file, but can override here
set distance    = ""
set beam_center = ""
set wavelength  = ""
set omega       = ""
set spin        = ""

# filenames
set outfile  = peaks.pdb
set tempfile = adsc2pdb_temp


goto Setup
Help:
############################################################
cat << EOF

usage: adsc2pdb.com 1.0A 3-20A 100mm 94x94 frame_1_0*.img ...
or
usage: adsc2pdb.com $outfile 6-20A zoom 4

where:
threshold 2 - the peak-picking threshold (default: $threshold)
1.0A        - the x-ray wavelength (optional)
3-20A	    - resolution range to output (optional)
100mm       - the xtf distance (optional)
94x94	    - beam center (denzo/mosflm convention)
frame_1_0*.img - the images you want to transform
$outfile   - a pdb file generated by adsc2pdb.com that
		you want to re-size or cut-off with different
		resolution limits.

examples:
adsc2pdb.com frame_1_*.img theshold 2
adsc2pdb.com peaks.pdb 6-1000A zoom 0.4

EOF
exit 9
Return_from_Setup:
if($?pdbfile) goto zoom_pdb
if("$wavelength" == "") goto Help
if("$distance" == "") goto Help
if($frames == 0) goto Help
if(! $?pdbfile) set pdbfile = "$outfile"
############################################################


# now we just run the awk program deployed below on the peak-pick file
echo "transforming spots into $outfile"
cat ${tempfile}peaks |\
nawk -f ${tempfile}RS2xyz.awk |\
cat >! $outfile

zoom_pdb:
if("$outfile" == "$pdbfile") set outfile = "new_$pdbfile"

echo "applying $hires - $lores A cutoff to $pdbfile"
echo "and zooming by a factor of ${zoom}X..."
cat $pdbfile |\
nawk -v hires=$hires -v lores=$lores -v zoom=$zoom '{\
    occ=substr($0,55,6)+0; B=substr($0,61,6)+0; \
    X=substr($0, 31, 8); Y=substr($0, 39, 8);Z=substr($0, 47, 8);}\
    B>hires && B<lores {n=$2; \
    printf "ATOM  %5d  O   WAT  %4d    %8.3f%8.3f%8.3f %5.2f%6.2f\n",\
           n,n/10,zoom*X,zoom*Y,zoom*Z,occ,B}' |\
cat >! $outfile

# clean up
#rm -f ${tempfile}peaks ${tempfile}RS2xyz.awk >& /dev/null


cat << EOF

! suggested o macro:
sam_atom_in $outfile peaks
mol peaks
paint_ramp atom_wt ; red blue ;
obj peaks
zone ;
end

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

EOF

exit



Setup:
################################################################################
# go down command-line args
set frames = 0
set i = 0
while ( $i < $#argv )
    @ i = ( $i + 1 )
    @ j = ( $i + 1 )
    @ k = ( $i + 2 )
    if($j > $#argv) set j = $#argv
    if($k > $#argv) set k = $#argv
    
    set Arg = "$argv[$i]"
    set arg = `echo "$Arg" | nawk '{print tolower($0)}'`
    
    # any command-line options
    if(-e "$arg") continue
    
    # look for numbers
    set num = `echo "$arg" | nawk '/^[0-9\.-]/{print $1+0}'`
    if("$num" != "") then
	# we have a number
	if(("$Arg" =~ *A)||("$argv[$j]" == "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: wavelength or resolution?
		set test = `echo $arg | nawk '$1+0>2{print $1+0}'`
		if("$test" != "") then
		    # it's a resolution
		    set hires = "$test"
		    continue
		endif
		# its a wavelength
		set user_wave = `echo "$arg" | nawk '{print $1+0}'`
		continue
	    endif
	endif
	
	if(("$Arg" =~ *X)||("$argv[$j]" == "X")) then
	    # zoom is expressed in X distance units
	    set zoom = "$num"
	    continue
	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? 
	    endif
	endif
	
	if(("$arg" =~ *mm)||("$argv[$j]" == "mm")) then
	    # distance is expressed in mm
	    set distance = `echo "$arg" | nawk '{print $1+0}'`
	    continue
	endif
    endif
    
    # look for word defining spindle direction
    if(("$arg" =~ spin*)&&("$argv[$j]" =~ anti*))  set spin = "anti"
    if(("$arg" =~ spin*)&&("$argv[$j]" =~ clock*)) set spin = "clock"

    # look for "next" numbers
    set num = `echo "$argv[$j]" | nawk '/^[0-9\.-]/{print $1+0}'`
    if("$num" == "") continue
    
    if("$arg" =~ thresh*) set threshold  = $num
    if("$arg" =~ dist*)   set distance   = $num
    if("$arg" =~ wave*)   set wavelength = $num
    if("$arg" =~ omega*)  set omega      = $num
    if("$arg" == zoom)    set zoom       = $num
    if(("$arg" == "beam")&&( $k <= $#argv)) then
	set temp = `echo $argv[$j] $argv[$k] | nawk '$NF~/^[0-9\.-]/'`
	if($#temp == 2) set beam_center = ( $temp )
    endif
end
# update any user-specified parameters
set user_beam  = ( $beam_center )
set user_dist  = "$distance"
set user_wave  = "$wavelength"
set user_omega = "$omega"
set user_spin  = "$spin"

echo "threshold: $threshold"
echo "zoom:      $zoom"

# 2nd pass for image files
echo "" >! ${tempfile}peaks
set i = 0
while ( $i < $#argv )
    @ i = ( $i + 1 )
    set arg = "$argv[$i]"
    
    # only interested in files this time
    if(! -e "$arg") continue
    
    if("$arg" =~ *.pdb) then
	set pdbfile = "$arg"
	continue
    endif
    
    # get a safe image header to read:
    set safehead = `od -cv $arg | nawk '{for(i=2;i<=NF;++i) print $i}' | nawk '/\\n/{++n;c=0} {++c} c>1000{exit} END{print n+0}'`
    if($safehead == 0) continue 
    
    echo -n `basename $arg`
    head -$safehead $arg >! ${tempfile}header    
    
    # read the ADSC header
    set distance    = `grep "DISTANCE" ${tempfile}header | nawk -F "=" '{ print $2+0}'`
    set phi         = `grep "PHI" ${tempfile}header | nawk -F "=" '{ print $2+0}'`
    set osc         = `grep "OSC_RANGE" ${tempfile}header | nawk -F "=" '{ print $2+0}'`
    set wavelength  = `grep "WAVELENGTH" ${tempfile}header | nawk -F "=" '{ print $2+0}'`

    # site-dependent stuff?
    set SN = `nawk -F "=" '/DETECTOR_SN/{ print $2+0}' ${tempfile}header`
    set omega = 0
    set spin  = clock
    if("$SN" == 401) set omega =  90
    if("$SN" == 401) set spin  = anti
    if("$SN" == 402) set omega =  90
    if("$SN" == 402) set spin  = anti
    if("$SN" == 406) set omega =  90
    if("$SN" == 406) set spin  = clock
    if("$SN" == 411) set omega =  90
    if("$SN" == 411) set spin  = anti
    if("$SN" =~ 44?) set omega =  90
    if("$SN" =~ 44?) set spin  = anti

    # beam center
    set x_beam      = `nawk -F "=" '/BEAM_CENTER_X/{ print $2+0}' ${tempfile}header`
    set y_beam      = `nawk -F "=" '/BEAM_CENTER_Y/{ print $2+0}' ${tempfile}header`
	
    # convert beam center to appropriate coordinates
    set PIXEL       = `grep "PIXEL_SIZE" ${tempfile}header | nawk -F "=" '{ print $2+0}'`
    set Width       = `grep "SIZE2" ${tempfile}header | nawk -F "=" '{ print $2+0}' | tail -1`
    set beam_center = `echo "($Width * $PIXEL) - $y_beam" | bc`
    set beam_center = `echo $beam_center $x_beam`
    set beam_center = `echo "$Width $PIXEL $x_beam $y_beam" | nawk '{print $1*$2-$NF, $3}'` 

    # ALS puts "denzo" beam center in headers
    if("$SN" == 401) set beam_center = `echo $x_beam $y_beam`
    
    # USER override of various parameters
    if("$user_dist" != "") set distance = $user_dist
    if("$user_wave" != "") set wavelength = $user_wave
    if("$user_beam" != "") set beam_center = ( $user_beam )
    if("$user_omega" != "") set omega = $user_omega
    if("$user_spin" != "") set spin = $user_spin
    
    # create a temporary information file
    echo " ${wavelength}A ${distance}mm  beam $beam_center[1] $beam_center[2]  phi $phi  omega $omega  spindle $spin"
    cat << EOF >! ${tempfile}header
PIXEL $PIXEL $PIXEL
DIST $distance
WAVE $wavelength
BEAM $beam_center
PHI   $phi
OSC   $osc
OMEGA $omega
SPIN  $spin
EOF
    
    # do the peak-search
#    $dps_peaksearch $arg ${tempfile}.dpspeaks $threshold >& /dev/null
    $dps_peaksearch $arg > /dev/null
    if($status) then
	echo "WARNING: error reading $arg"
	continue
    endif
    cat peaks.list |\
    awk '{print $2,$1,$NF}' |\
    cat >! ${tempfile}.dpspeaks
    
    # count of frames sucessfully read
    @ frames = ( $frames + 1 )
    
    # concatenate peak files (for latter...)
    cat ${tempfile}header ${tempfile}.dpspeaks |\
    cat >> ${tempfile}peaks
end
rm -f ${tempfile}header ${tempfile}.dpspeaks >& /dev/null


################################################################################
# deploy the awk program for converting these
cat << EOF-awk >! ${tempfile}RS2xyz.awk
#! /usr/bin/awk -f
#
#   Awk program for converting an RS-type spot file
#   into reciprocal-space coordinates
#   phi, distance, etc. are indicated with keywords
#
#   the conventions of omega and spin are analogous to the 
#   DPS "rotation" and "yscale" parameters
#

BEGIN{
    Zd         = 150
    lambda     = 1.1
    Xbeam      = 94
    Ybeam      = 94
    OMEGA      = 180
    SPIN       = -1
    cutoff     = 0
    hires      = 3
    lores      = 10000
    zoom       = 1
    
    PI=2*atan2(1,0)
}

# "key-worded" input can be concatenated in front of the spot data
/^PIXEL/{Xpix=\$2;Ypix=\$3}
/^DIST/{Zd=\$2}
/^WAVE/{lambda=\$2}
/^BEAM/{Xbeam=\$2;Ybeam=\$3}
/^PHI/{PHI=\$2}
/^OSC/{OSC=\$2}
/^OMEG/{OMEGA=\$2*PI/180}
/^SIGN/{SPIN=\$2}
/^SPIN/ && /clock/{SPIN=1}
/^SPIN/ && /anti/{SPIN=-1}
/^ZOOM/{zoom=\$2}
/^RESO/{lores=\$2;hires=\$3; if(lores<hires){hires=\$2;lores=\$3}}

/^[A-Z]/{next}

NF>=2{
    # using mosflm standard detector coordinate system:
    #    image: Ys || fast axis     Xs || slow      origin: first pixel
    # detector: Zd || x-ray beam,   Yd || spindle   origin: direct beam
    #   camera:  x || x-ray beam,    z || spindle   origin: r.l. origin
    
    # read in a peak (Ys=fast, Xs=slow)
    Ys=\$1; Xs=\$2;
    
    # convert to mm coordinates, centered on direct beam
    Xmm=Xpix*Xs-Ybeam;
    Ymm=Ypix*Ys-Xbeam;
    
    # rotate the detector so that Yd is parallel to phi axis
    Xd  = Xmm*cos(OMEGA)-Ymm*sin(OMEGA);
    Yd  = Xmm*sin(OMEGA)+Ymm*cos(OMEGA);
    
    # now map each spot back onto the surface of the Ewald sphere
    # we are also converting from detector axes to camera axes.
    distortion = lambda*sqrt((Xd)^2+(Yd)^2+(Zd)^2)
    yo = Xd/distortion
    zo = Yd/distortion
    xo = Zd/distortion - 1/lambda
    
    # (xo, yo, zo) are now the dimensionless coordinates of the peak
    # in the Ewald construction at phi=PHI
    
    # reverse the phi to restore (xo,yo,zo) to their positions when phi==0.0
    # and convert phi to right-handed radians
    # (phi1=start;phi2=end)
    phi1 = -SPIN*PI/180*(PHI)
    phi2 = -SPIN*PI/180*(PHI+OSC)
    
    # apply the phi rotation about "z"
    x1 = yo*sin(phi1)+xo*cos(phi1)
    y1 = yo*cos(phi1)-xo*sin(phi1)
    z1 = zo
    x2 = yo*sin(phi2)+xo*cos(phi2)
    y2 = yo*cos(phi2)-xo*sin(phi2)
    z2 = zo

    # 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.1f%6.2f\\n",
           n,n/10,zoom*Zd*x1,zoom*Zd*y1,zoom*Zd*z1,PHI+OSC/2,d;
    printf "ATOM  %5d  O   WAT  %4d    %8.3f%8.3f%8.3f %5.1f%6.2f\\n",
           n,n/10,zoom*Zd*x2,zoom*Zd*y2,zoom*Zd*z2,PHI+OSC/2,d;
}
EOF-awk



goto Return_from_Setup