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

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

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

checkversion 1.3

set display=&output end


declare domain=reci type=real name=centr end 

do (centr=0) (acentric)
do (centr=1) (centric)

evaluate ($pp_tiny=1E-9)
evaluate ($pp_max=10)
		{*print data*}

if (&epsilon="yes") then

   display                                      generalized Cullis R-value
   display                                      --------------------------
   display                                      sqrt<loc^2>/sqrt<delta F^2>
   display                                      centrics  acentrics    all
   if (&text # "") then
      display &text
   end if
   bins ? evaluate ($old_bins=$result) 
   statistics 
      (min($pp_max,sqrt(save(centr* &var ))/max($pp_tiny,sqrt(save(centr* (abs(&fp)-abs(&fph))^2)))))
      (min($pp_max,sqrt(save((1-centr)* &var ))/max($pp_tiny,sqrt(save((1-centr)* (abs(&fp)-abs(&fph))^2 )))))
      (min($pp_max,sqrt(save(&var))/max($pp_tiny,sqrt(save( (abs(&fp)-abs(&fph))^2 )))))
      sele=(&sel)
      output=&output
   end 
   bins=1
   statistics 
      (min($pp_max,sqrt(save(centr* &var ))/max($pp_tiny,sqrt(save(centr* (abs(&fp)-abs(&fph))^2 )))))
      (min($pp_max,sqrt(save((1-centr)* &var ))/max($pp_tiny,sqrt(save((1-centr)* (abs(&fp)-abs(&fph))^2 )))))
      (min($pp_max,sqrt(save(&var))/max($pp_tiny,sqrt(save( (abs(&fp)-abs(&fph))^2 )))))
       sele=(&sel)
      output=&output
   end 
   bins=$old_bins
   if (&table # "") then
     display TAB: &table R-Cullis(centric)= $expression1[F7.3] R-Cullis(acentric)= $expression2[F7.3] \
R-Cullis(all)= $expression3[F7.3]
   end if

else

   display                                      generalized Cullis R-value
   display                                      --------------------------
   display                                      sqrt<loc^2>/sqrt<delta F^2>
   display                                      centrics  acentrics    all
   if (&text # "") then
      display &text
   end if
   bins ? evaluate ($old_bins=$result) 
   statistics 
      (min($pp_max,sqrt(ave(centr* &var ))/max($pp_tiny,sqrt(ave(centr* (abs(&fp)-abs(&fph))^2 )))))
      (min($pp_max,sqrt(ave((1-centr)* &var ))/max($pp_tiny,sqrt(ave((1-centr)* (abs(&fp)-abs(&fph))^2 )))))
      (min($pp_max,sqrt(ave(&var))/max($pp_tiny,sqrt(ave( (abs(&fp)-abs(&fph))^2 )))))
       sele=(&sel)
      output=&output
   end 
   bins=1
   statistics 
      (min($pp_max,sqrt(ave(centr* &var ))/max($pp_tiny,sqrt(ave(centr* (abs(&fp)-abs(&fph))^2 )))))
      (min($pp_max,sqrt(ave((1-centr)* &var ))/max($pp_tiny,sqrt(ave((1-centr)* (abs(&fp)-abs(&fph))^2 )))))
      (min($pp_max,sqrt(ave(&var))/max($pp_tiny,sqrt(ave( (abs(&fp)-abs(&fph))^2 )))))
       sele=(&sel)
      output=&output
   end 
   bins=$old_bins   
   if (&table # "") then
      display TAB: &table R-Cullis(centric)= $expression1[F7.3] R-Cullis(acentric)= $expression2[F7.3] \
R-Cullis(all)= $expression3[F7.3]
   end if

end if

undeclare domain=reci name=centr end