! Module file: printphasingpower
!
! CNS module
! **********
!
! Authors: Piet Gros and Axel T. Brunger
!
! copyright Yale University
!
! Function: 
! Print statistics on the phasing power <Fh>/<lack-of-closure>
!
! 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}

   &fh;                 {heavy atom structure factors}
   &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)
		


if (&epsilon="yes") then

   display                                     phasing power
   display                                     -------------
   display                                     sqrt<Fh^2>/sqrt<var_loc^2>  
   display                                     centrics  acentrics    all
   if (&text # "") then
      display &text
   end if
   
   bins ? evaluate ($old_bins=$result) 
   statistics 
      (sqrt(save(centr*abs(&fh)^2))/max($pp_tiny,sqrt(save(centr*&var))))
      (sqrt(save((1-centr)*abs(&fh)^2))/max($pp_tiny,sqrt(save((1-centr)*&var))))
      (sqrt(save(abs(&fh)^2))/max($pp_tiny,sqrt(save(&var))))
        sele=(&sel)
      output=&output
   end 
   bins=1
   statistics 
      (sqrt(save(centr*abs(&fh)^2))/max($pp_tiny,sqrt(save(centr*&var))))
      (sqrt(save((1-centr)*abs(&fh)^2))/max($pp_tiny,sqrt(save((1-centr)*&var))))
      (sqrt(save(abs(&fh)^2))/max($pp_tiny,sqrt(save(&var))))
      sele=(&sel)
      output=&output
   end 
   bins=$old_bins
   if (&table # "") then
      display TAB: &table ph(centric)= $expression1[F7.3] ph(acentric)= $expression2[F7.3] \
ph(all)= $expression3[F7.3]
   end if

else

   display                                     phasing power
   display                                     -------------
   display                                     sqrt<Fh^2>/sqrt<var_loc^2> 
   display                                     centrics  acentrics    all
   if (&text # "") then
      display &text
   end if
   
   bins ? evaluate ($old_bins=$result) 
   statistics 
      (sqrt(ave(centr*abs(&fh)^2))/max($pp_tiny,sqrt(ave(centr*&var))))
      (sqrt(ave((1-centr)*abs(&fh)^2))/max($pp_tiny,sqrt(ave((1-centr)*&var))))
      (sqrt(ave(abs(&fh)^2))/max($pp_tiny,sqrt(ave(&var))))
      sele=(&sel)
      output=&output
   end 
   bins=1
   statistics 
      (sqrt(ave(centr*abs(&fh)^2))/max($pp_tiny,sqrt(ave(centr*&var))))
      (sqrt(ave((1-centr)*abs(&fh)^2))/max($pp_tiny,sqrt(ave((1-centr)*&var))))
      (sqrt(ave(abs(&fh)^2))/max($pp_tiny,sqrt(ave(&var))))
      sele=(&sel)
      output=&output
   end 
   bins=$old_bins
   if (&table # "") then
      display TAB: &table ph(centric)= $expression1[F7.3] ph(acentric)= $expression2[F7.3] \
ph(all)= $expression3[F7.3]
   end if

end if

undeclare domain=reci name=centr end