{+ file: mad_bijvoet_ave.inp +}
{+ directory: xtal_phase +}
{+ description: Average amplitudes at different wavelengths, keeping 
                Friedel mates separate. +}
{+ comment:
            The anomalous diffraction ratios are used for relative
            weighting between wavelengths. +}
{+ authors: Axel T. Brunger +}
{+ copyright: Yale University +}

{- Guidelines for using this file:
   - all strings must be quoted by double-quotes
   - logical variables (true/false) are not quoted
   - do not remove any evaluate statements from the file -}

{- begin block parameter definition -} define(

{====================== crystallographic data ========================}

{* space group *}
{* use International Table conventions with subscripts substituted
   by parenthesis *}
{===>} sg="P2(1)2(1)2(1)";

{* unit cell parameters in Angstroms and degrees *}
{+ table: rows=1 "cell" cols=6 "a" "b" "c" "alpha" "beta" "gamma" +}
{===>} a=65.508;
{===>} b=72.216;
{===>} c=45.035;
{===>} alpha=90;
{===>} beta=90;
{===>} gamma=90;

{* input filename with merged diffraction data *}
{===>} reflection_infile="mbp_scale.hkl";

{* use anomalous diffraction ratio to weigh data *}
{* if false weights are set to 1 *}
{+ choice: true false +}
{===>} weight=true;

{* number of bins for empirical averaging weight calculation *}
{===>} bins=10;

{* number wavelengths to be averaged *}
{===>} n_wave=2;

{* amplitude of 1st wavelength to be averaged *}
{===>} f_w_1="f_w2";
{* sigma of 1st wavelength to be averaged *}
{===>} s_w_1="s_w2";

{* amplitude of 2nd wavelength to be averaged *}
{===>} f_w_2="f_w3";
{* sigma of 2nd wavelength to be averaged *}
{===>} s_w_2="s_w3";

{* amplitude of 3rd wavelength to be averaged *}
{===>} f_w_3="";
{* sigma of 3rd wavelength to be averaged *}
{===>} s_w_3="";

{* amplitude of 4th wavelength to be averaged *}
{===>} f_w_4="";
{* sigma of 4th wavelength to be averaged *}
{===>} s_w_4="";

{* amplitude of 5th wavelength to be averaged *}
{===>} f_w_5="";
{* sigma of 5th wavelength to be averaged *}
{===>} s_w_5="";

{=========================== output files ============================}

{* name for average amplitudes *}
{===>} f_ave="f_ave";
{* name for corresponding sigmas *}
{===>} s_ave="s_ave";

{* output reflection filename which contains averaged amplitudes and sigmas *}
{===>} reflection_outfile="mad_bijvoet_ave.hkl";

{===========================================================================}
{         things below this line do not normally need to be changed         }
{===========================================================================}

 ) {- end block parameter definition -}


 checkversion 1.3

 evaluate ($log_level=quiet)


xray                                    

   @CNS_XTALLIB:spacegroup.lib (sg=&sg;sgparam=$sgparam;) 

   a=&a b=&b c=&c  alpha=&alpha beta=&beta gamma=&gamma

   reflection @@&reflection_infile end

   set echo=off end
   evaluate ($i_wave=1)
   while ($i_wave <= &n_wave) loop ccc
      if ( &BLANK%f_w_$i_wave = false ) then
         query name=&STRIP%f_w_$i_wave domain=reciprocal end
         if ( $object_exist = false ) then
            display 
            display ****************************************************************
            display wavelength $i_wave : amplitude array &f_w_$i_wave does not exist
            display ****************************************************************
            display 
            abort
         end if
         {- this array can be of any type -}
         
         if ( &BLANK%s_w_$i_wave = true ) then
            display 
            display **********************************************
            display wavelength $i_wave : sigma array not specified
            display **********************************************
            display 
            abort
         else
            query name=&STRIP%s_w_$i_wave domain=reciprocal end
            if ( $object_exist = false ) then
               display 
               display ************************************************************
               display wavelength $i_wave : sigma array &s_w_$i_wave does not exist
               display ************************************************************
               display 
               abort
            end if
            if ( $object_type # "REAL" ) then
               display 
               display *********************************************************************
               display wavelength $i_wave : sigma array &s_w_$i_wave has the wrong data type
               display *********************************************************************
               display
               abort
            end if
         end if
       end if
      evaluate ($i_wave=$i_wave+1) 
   end loop ccc 
   
   if ( &BLANK%f_ave = true ) then
      display 
      display ********************************************
      display output average amplitude array not specified
      display ********************************************
      display 
      abort
   elseif ( &BLANK%s_ave = true ) then
      display 
      display ****************************************
      display output average sigma array not specified
      display ****************************************
      display 
      abort
   end if
   
   set echo=on end

   declare name=&STRIP%f_ave domain=reciprocal type=real end
   declare name=&STRIP%s_ave domain=reciprocal type=real end
   declare name=count domain=reciprocal type=integer end
   declare name=www   domain=reciprocal type=real end
   
   do (&STRIP%f_ave=0) ( all )
   do (&STRIP%s_ave=0) ( all ) 
   do (count=0) ( all ) 

   show max (d) (all )
   evaluate ($madbave_low_res=$result+0.001) 
   show min (d) (all )
   evaluate ($madbave_high_res=$result-0.001) 
   binresolution $madbave_low_res $madbave_high_res
  
   evaluate ($ii=1)
   
   while ($ii <= &n_wave) loop main

      if (&weight=true) then
         do (www=0) ( all ) 
         do (www=(sqrt(save(abs(amplitude(&STRIP%f_w_$ii)-amplitude(friedel(&STRIP%f_w_$ii)))^2 ))/
                     sqrt(save( ((amplitude(&STRIP%f_w_$ii)+amplitude(friedel(&STRIP%f_w_$ii)) )/2)^2)) ) )                                                
                    ( friedel_pair( abs(&STRIP%f_w_$ii)>0 ) )
         do (www=1/www) ( abs(www)>0) 
      else
         do (www=1.) ( all )
      end if

      do (&STRIP%f_ave=&STRIP%f_ave + www*&STRIP%f_w_$ii)      ( abs(&STRIP%f_w_$ii)>0 )
      do (&STRIP%s_ave=&STRIP%s_ave + www*&STRIP%s_w_$ii ^ 2 ) ( abs(&STRIP%f_w_$ii)>0 )
      do (count=count+www) ( abs(&STRIP%f_w_$ii)>0 )
   
      statistics
         ( www ) 
         selection= ( abs(www)>0 )  
         output=OUTPUT
      end
      
      eval ($ii=$ii+1)       
   end loop main
  
   do (&STRIP%f_ave=&STRIP%f_ave / count) ( count > 0 ) 
   do (&STRIP%s_ave=sqrt ( &STRIP%s_ave / count ) ) ( count > 0 ) 
   
end

set display=&reflection_outfile end

@CNS_XTALMODULE:write_hkl_header (sg=&STRIP%sg;
                                  sgparam=$sgparam;)

xray
   
   write reflection output= &reflection_outfile &STRIP%f_ave &STRIP%s_ave end 
   
end
  
set display=OUTPUT end
 
stop