! Module file: printanorcullis
!
! CNS module
! **********
!
! Authors: Piet Gros and Axel T. Brunger
!
! copyright Yale University
!
! Function: 
! Print statistics on the generalized anomalous Cullis R factor
!
! Requirements:
! 1. This module can only be called from xray
! 2. probability must be precomputed
!

module {printanorcullis}
(
   &text="";            {identifying text}
   &table="";           {identifying text for table}
   &epsilon="yes";      {"yes" | "no", use epsilon weighting in summing variances}

   &fph;               {derivative data}
   &var;                {lack-of-closure variance}
   &sel=(all);          {structure factor selection}
   &output=OUTPUT;      {output filename}
)

checkversion 1.3

set display=&output end

evaluate ($pp_tiny=1E-9)
evaluate ($pp_max=10)
		
declare domain=reci type=real name=selected end   {selected bijvoet pairs}


do (selected=0) ( all )      {select Bijvoet pairs}
do (selected=1) ( friedel_pair(&sel) and acentric )

if (&epsilon="yes") then

   display                          anomalous Cullis R-value
   display                          ------------------------
   display                          R-Cullis= sqrt<loc^2>/sqrt<(F+ - F-)^2>
   if (&text # "") then
      display &text
   end if
   
   bins ? evaluate ($old_bins=$result) 
   statistics 
      (min($pp_max,sqrt(save(&var))/max($pp_tiny,sqrt(save( ( abs(&fph)-abs(friedel(&fph)) )^2 )))))
      sele=(selected=1)
      output=&output
   end 
   bins=1
   statistics 
         (min($pp_max,sqrt(save(&var))/max($pp_tiny,sqrt(save( ( abs(&fph)-abs(friedel(&fph)) )^2 )))))
      sele=(selected=1)
      output=&output
   end 
   bins=$old_bins
   if (&table # "") then
      display TAB: &table anomalous R-Cullis= $expression1[F7.3]
   end if

else

   display                          anomalous Cullis R-value
   display                          ------------------------
   display                          R-Cullis= sqrt<loc^2>/sqrt<(F+ - F-)^2>
   if (&text # "") then
      display &text
   end if
   
   bins ? evaluate ($old_bins=$result) 
   statistics 
      (min($pp_max,sqrt(ave(&var))/max($pp_tiny,sqrt(ave( ( abs(&fph)-abs(friedel(&fph)) )^2 )))))
      sele=(selected=1)
      output=&output
   end 
   bins=1
   statistics 
         (min($pp_max,sqrt(ave(&var))/max($pp_tiny,sqrt(ave( ( abs(&fph)-abs(friedel(&fph)) )^2 )))))
      sele=(selected=1)
      output=&output
   end 
   bins=$old_bins
   if (&table # "") then
      display TAB: &table anomalous R-Cullis= $expression1[F7.3]
   end if

end if

undeclare domain=reci name=selected end   {selected bijvoet pairs}