#ifndef STANDARD_MACROS 
#define STANDARD_MACROS
#include local.inc
macro Pi  { 3.14159265358979 }
macro Deg { 0.01745329251994 }
macro Deg_on_2 { 0.00872664625997 }
macro Rad { 57.2957795130823 }
macro Double_Quote { " }
macro IF { If( }
macro THEN { , }
macro ELSE { , }
macro ENDIF { ) }
macro Repeat(x, & n) 
	{              
      #if (n > 0) 
         x Repeat(x, Numeric(n-1)) 
      #endif 
   }         
macro Create_XDDs(n) 
   {
      move_to xdds
      move_to xdd_recs
      new(xdd, n)
   }
macro Numeric(n) 
   {
      #prm a = n; #out a
   }
macro Make_I_Local(c, v)
   {
      = c; c v min 0 
   }
macro Include_GUI_PDF 
	{
      #include Double_Quote##ROOT##GUI_PDF.inc##Double_Quote
	}
macro Include_Charge_Flipping 
   { 
      #include Double_Quote##ROOT##Charge_Flipping.inc##Double_Quote
   } 
macro Include_PDF_Generate 
   { 
      #include Double_Quote##ROOT##pdf-generate.inc##Double_Quote
   } 
macro Include_PDF_Generate_Common 
   { 
      #include Double_Quote##ROOT##pdf-generate-common.inc##Double_Quote
   } 
macro Not_Blank(s)
   {
      #m_ifarg s "" 
         0   
      #m_else
         1
      #m_endif
   }
macro & Odd(n)
   {
      Mod(n, 2)
   }
macro & Even(n)
   {
      Mod(n, 2) == 0
   }
macro Transition(t, l, h)
   {
      transition t
         use_layer l
         height = h;
   }
macro Transition(l, h)
   {
      transition l
         use_layer l
         height = h;
   }
macro SF_smooth(c, v, & s)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min 1e-6 max 20)
      hkl_I_d2Th_dL #m_unique idl
      gauss_fwhm = CeV(c, v) idl;
      th2_offset = Get(hkl_I_d2Th_dL, offset);
      peak_buffer_based_on = idl; peak_buffer_based_on_tol = Peak_Calculation_Step 0.25 s;
   }
macro SF_smooth(s)
   {
      SF_smooth(,1,s)
   }
macro Append_0(file)
   {
      out file
      #if (Run_Number > 0)
         append
      #endif
   }
macro Backup_INP
   {
      system_before_save_OUT  {
         copy INP_File##.inp INP_File##.backup
      }
   }
macro Save_Best
   {
      #if (Run_Number == 0)
         prm Best_Rwp_ = 1e99;		  
      #else
         prm Best_Rwp_ = #include Best_Rwp_.txt;
      #endif			
		out Best_Rwp_.txt Out(If(Get(r_wp) < Best_Rwp_, Get(r_wp), Best_Rwp_))
		out_file = If(Get(r_wp) < Best_Rwp_, Concat(String(INP_File), ".OUT"), "");
   }
macro File_Variable(c, x_start, dx)
   {
      #if (Run_Number == 0)
         #prm c = x_start;
      #else
         #prm c = #include c##.txt;
      #endif
      #prm c##_next = c + dx;
      out c##.txt Out(#out c##_next)
   }
macro File_Variables(a, a1, a2, da, b, b1, b2, db)
   {
      #if (Run_Number == 0)
         #prm a = a1;
         #prm b = b1;
      #else
         #prm a = #include a##.txt;
         #prm b = #include b##.txt;
      #endif
      #prm a##_next = If(b >= b2, a + da, a);
      #prm b##_next = If(b >= b2, b1, b + db);
      out a##.txt Out(#out a##_next)
      out b##.txt Out(#out b##_next)
   }
macro Randomize_File_Out_Normal(file)
   {
      xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f  " = X;
          " %11.6f\n  " = Max(Rand_Normal(Ycalc, Sqrt(Ycalc)), 0);
      }
   }
macro TOF_XYE(path, calc_step)
   {
      xdd path xye_format
         x_calculation_step calc_step
         neutron_data
         weighting = If(SigmaYobs < 1, 1, 1/SigmaYobs^2);
   }
macro TOF_GSAS(path, calc_step)
   {
      xdd path gsas_format
         neutron_data
         x_calculation_step calc_step
         weighting = If(SigmaYobs < 1, 1, 1/SigmaYobs^2);
   }
macro TOF_CS_L(c, v, & t1)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min 3 max 10000)
      lor_fwhm = Constant(t1 .1) D_spacing^2 / CeV(c, v);
   }
macro TOF_CS_G(c, v, & t1)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min 3 max 10000)
      gauss_fwhm = Constant(t1 .1) D_spacing^2 / CeV(c, v);
   }
macro TOF_Exponential(a0, a0v, a1, a1v, & wexp, & t1, lr)
   {
      #m_argu a0
      #m_argu a1
      If_Prm_Eqn_Rpt(a0, a0v, min = Max(Val .3, 1e-6); max = 2 Val + 1; )
      If_Prm_Eqn_Rpt(a1, a1v, min = Max(Val .3, 1e-6); max = 2 Val + 1; )
      exp_conv_const = lr Constant(t1) / (CeV(a0,a0v) + CeV(a1,a1v) / D_spacing^wexp);
   }
macro TOF_PV(fwhm, fwhmv, lor, lorv, & t1)
   {
      #m_argu fwhm
      If_Prm_Eqn_Rpt(fwhm, fwhmv, min .0001 max = 2 Val + 15 Yobs_dx_at(.5 (X1 + X2)); )
      peak_type pv
         pv_lor lor lorv
         pv_fwhm = CeV(fwhm, fwhmv) Constant(t1 0.00001) D_spacing;
   }
macro TOF_x_axis_calibration(t0, t0v, t1, t1v, t2, t2v)
   {
      #m_argu t0
      #m_argu t1
      #m_argu t2
      If_Prm_Eqn_Rpt(t0, t0v,
         min = Val - 15 Yobs_dx_at(.5(X1+X2));
         max = Val + 15 Yobs_dx_at(.5(X1+X2));
         del = .3 Peak_Calculation_Step; )
      If_Prm_Eqn_Rpt(t1, t1v,
         min = Val - 15 Yobs_dx_at(.5(X1+X2));
         max = Val + 15 Yobs_dx_at(.5(X1+X2));
         del .01)
      If_Prm_Eqn_Rpt(t2, t2v,
         min = Val - 15 Yobs_dx_at(.5(X1+X2));
         max = Val + 15 Yobs_dx_at(.5(X1+X2));
         del .005)
      pk_xo = CeV(t0,t0v) + CeV(t1,t1v) D_spacing + CeV(t2,t2v) D_spacing^2;
      gui_tof_t0 = CeV(t0,t0v);
      gui_tof_t1 = CeV(t1,t1v);
   }
macro TOF_LAM(w)
   {
      lam
         ymin_on_ymax w
         la 1 lo 0 lh 1
   }
macro Scherrer(& k, & fwhm, & s, csize)
   {
      local = If(Abs(fwhm) < Abs(s),
         -9999,
         .1 Rad Lam k / (Cos(Get(xo) Deg_on_2) Sqrt(fwhm^2-s^2))
         ); : csize
   }
macro e0_from_Strain(e0, sgc, sgv, slc, slv)
   {
      #m_argu sgc
      #m_argu slc
      local = Voigt_FWHM_GL(CeV_or_0(sgc, sgv), CeV_or_0(slc, slv)) .25 Pi/360; : e0
      #m_ifarg sgv "" #m_else
         Strain_G(sgc, sgv)
      #m_endif
      #m_ifarg slv "" #m_else
         Strain_L(slc, slv)
      #m_endif
   }
macro IB_from_CS(& csg, & csl)
   {
      #m_ifarg csg ""
         (1.570796326795/csl)
      #m_else
         #m_ifarg csl ""
            (1.064467019431/csg)
         #m_else
            Voigt_Integral_Breadth_GL(1/csg, 1/csl)
         #m_endif
      #m_endif
   }
macro Voigt_FWHM_from_CS(& csg, & csl)
   {
      #m_ifarg csg ""
         (1/csl)
      #m_else
         #m_ifarg csl ""
            (1/csg)
         #m_else
            Voigt_FWHM_GL(1/csg, 1/csl)
         #m_endif
      #m_endif
   }
macro LVol_FWHM_CS_G_L(& k, lvol, & kf, lvolf, csgc, csgv, cslc, cslv)
   {
      #m_argu csgc
      #m_argu cslc
      local = k / IB_from_CS(CV(csgc, csgv), CV(cslc, cslv)); : lvol
      local = kf / Voigt_FWHM_from_CS(CV(csgc, csgv), CV(cslc, cslv)); : lvolf
      #m_ifarg csgv "" #m_else
         CS_G(csgc, csgv)
      #m_endif
      #m_ifarg cslv "" #m_else
         CS_L(cslc, cslv)
      #m_endif
   }
macro PVII_m(v) {   PVII_m(,v) }
macro PVII_m(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min .55 max 30, del = 0.05 Val + 0.01; )
      m1 = CeV(c, v);
      m2 = Get(m1);
   }
macro PVII_h(v) {   PVII_h(,v) }
macro PVII_h(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min .0001 max 20)
      h1 = CeV(c, v) .5;
      h2 = Get(h1);
   }
macro SPVII_to_PVII_Peak_type { peak_type spvii }
macro CV(c, v)
   {
      #m_ifarg c ""
			#m_first_word v v
      #m_else
         c
      #m_endif
   }
macro & CeV(c, v)
   {
      CV(c, v)
   }
macro If_c_or_v_Ok(c, v, w)
   {
      #m_ifarg c ""
         #m_ifarg v ""
         #m_else
            w
         #m_endif
      #m_else
           w
      #m_endif
   }
macro If_Not_Blank(c, w)
   {
      #m_ifarg c ""
      #m_else
         w
      #m_endif
   }
macro Fix_If_Eqn(c, v)
   {
      #m_ifarg c #m_eqn
         = c; : v
      #m_else
         c v
      #m_endif
   }
macro CeV_or_0(c, v)
   {
      #m_ifarg v ""
         0
      #m_else
         CeV(c, v)
      #m_endif
   }
macro If_Prm_Eqn_Rpt(c, v, det1, det2)
   {
      #m_ifarg c #m_code
         prm c
         #m_ifarg c #m_code_refine
            #m_ifarg v #m_one_word
               v det1 det2
            #m_else
               v det2
            #m_endif
         #m_else
            #m_first_word v v
         #m_endif
      #m_else
         #m_ifarg c #m_eqn
            prm = c; : #m_first_word v v
         #m_endif
      #m_endif
   }
macro If_Prm_Eqn_Rpt(c, v, dets)
   {
      If_Prm_Eqn_Rpt(c, v, dets,)
   }
' ---------------------------------------
macro Angle_Distance_Restrain(& c, & t, t_calc, & tol, & wscale)
   {
      local = c; : t_calc
      penalty =
         wscale
         IF tol THEN
            IF c < (t - tol) THEN
               (c - (t - tol) )^2
            ELSE
               IF c > (t + tol) THEN
                  (c - (t + tol))^2
               ELSE
                  0
               ENDIF
            ENDIF
         ELSE
            (c - t)^2
         ENDIF
         ;
   }
macro Distance_Restrain(sites, c, t, t_calc, tol, wscale)
   {
      #m_ifarg c ""
         #m_unique_not_refine c
      #m_endif
      sites_distance c
         load site_to_restrain { sites }
      Angle_Distance_Restrain(c, t, t_calc, tol, wscale)
   }
macro Distance_Restrain(sites, t, t_calc, tol, wscale)
   {
      Distance_Restrain(sites,, t, t_calc, tol, wscale)
   }
macro Angle_Restrain(sites, c, t, t_calc, tol, wscale)
   {
      #m_ifarg c ""
         #m_unique_not_refine c
      #m_endif
      sites_angle c
         load site_to_restrain { sites }
      Angle_Distance_Restrain(c, t, t_calc, tol, wscale)
   }
macro Angle_Restrain(sites, t, t_calc, tol, wscale)
   {
      Angle_Restrain(sites,, t, t_calc, tol, wscale)
   }
macro Flatten(sites, c, t_calc, tol, & wscale)
   {
      #m_ifarg c ""
         #m_unique_not_refine c
      #m_endif
      sites_flatten c
         sites_flatten_tol = tol;
         load site_to_restrain { sites }

      penalty = c wscale; : t_calc
   }
macro Flatten(sites, t_calc, tol, wscale)
   {
      Flatten(sites,,t_calc, tol, wscale)
   }
' ---------------------------------------
macro Distance_Restrain_Keep_Within(sites, & r, & c, & wby, & num_cycle_iters)
   {
      #m_ifarg c ""
         #m_unique_not_refine c
      #m_endif
      sites_distance c
         load site_to_restrain { sites }
         penalty = wby
            IF And(Cycle_Iter < num_cycle_iters, c > r) THEN
               (c - r)^2
            ELSE
               0
            ENDIF
            ;
   }
macro Distance_Restrain_Keep_Within(sites, r, wby, num_cycle_iters)
   {
      Distance_Restrain_Keep_Within(sites, r,, wby, num_cycle_iters)
   }
macro Distance_Restrain_Keep_Out(sites, r, c, wby, num_cycle_iters)
   {
      #m_ifarg c ""
         #m_unique_not_refine c
      #m_endif
      sites_distance c
         load site_to_restrain { sites }
         penalty =
            IF And(Cycle_Iter < num_cycle_iters, c < r) THEN
               (c - r)^2
            ELSE
               0
            ENDIF
            ;
   }
macro Distance_Restrain_Keep_Out(sites, r, wby, num_cycle_iters)
   {
      Distance_Restrain_Keep_Out(sites, r,, wby, num_cycle_iters)
   }
' ---------------------------------------
macro AI_Anti_Bump(s1, s2, ro, c, wby, num_cycle_iters)
   {
      #m_ifarg c ""
         #m_unique_not_refine c
      #m_endif
      ai_anti_bump c
         ai_radius = ro;
         ai_sites_1 s1
         ai_sites_2 s2
         
         #m_ifarg num_cycle_iters ""
            penalty = (wby) c;
         #m_else
            penalty = (wby) If(Cycle_Iter < (num_cycle_iters), c, 0);
         #m_endif
   }
macro AI_Anti_Bump(s1, s2, ro, wby, num_cycle_iters)
   {
      AI_Anti_Bump(s1, s2, ro,, wby, num_cycle_iters)
   }
macro AI_Anti_Bump(s1, s2, ro, wby)
   {
      AI_Anti_Bump(s1, s2, ro,, wby,)
   }
macro Parabola_N(n1, n2, s1, s2, c, ro, wby)
   {
      #m_ifarg c ""
         #m_unique_not_refine c
      #m_endif
      box_interaction to_N n2 s1 s2 c =
         If(Or(Ri < n1, R < ro),
            (R-ro)^2,
            0
         );
      penalty = (wby) c;
   }
macro Parabola_N(n1, n2, s1, s2, ro, wby)
   {
      Parabola_N(n1, n2, s1, s2,, ro, wby)
   }
macro Rn_Rep(& q, & ro, & n)
   {
      (q / n) (ro^(n-1)) / R^n
   }
macro Box_Rn_Att_Rep(s1, s2, c, q, ro, n)
   {
      box_interaction s1 s2 c = q /R + Rn_Rep(Sqrt((q) (q)), ro, n);
   }
macro Grs_Interaction(s1, s2, & wqi, & wqj, c, ro, n, rsm)
   {
		grs_interaction qi = wqi; qj = wqj; s1 s2 c =
         If (R < rsm,
            (Abs( Get(qi) Get(qj) ) / n) (ro^(n-1)) ( (-n rsm^(-2 - n)/2) R^2 + rsm^(-n) + n/(2 rsm^n)),
            Rn_Rep(Abs( Get(qi) Get(qj) ), ro, n)
         );
   }		  
macro Grs_Interaction(s1, s2, wqi, wqj, c, ro, n)
   {
      Grs_Interaction(s1, s2, wqi, wqj, c, ro, n, 0.2)
   }
macro Grs_No_Repulsion(s1, s2, wqi, wqj, c)
   {
      grs_interaction qi = wqi; qj = wqj; s1 s2 c = 0;
   }
macro Box_Gauss_Rep_n(s1, s2, c, ro, y, n)
   {
      box_interaction s1 s2 c = Gauss_Rep_n(ro, y, n);
   }
macro Gauss_Rep_n(ro, y, n)
   {
      (Exp(R^n Ln(y) / ro^n))
   }
macro Exp_Rep(q, ro, b)
   {
      (Exp(- b R) Exp(b ro) Sqrt((q)^2) / (b ro^2))
   }
macro Exp_Rep_n(q, ro, b, n)
   {
      (Exp(b (ro - R)) n Sqrt((q)^2) / (b ro^(n+1)))
   }
macro Born_Mayer(q, ro, b)
   {
      (q / R + Exp_Rep(q, ro, b))
   }
macro Born_Mayer_n(q, ro, b, n)
   {
      (q / R^n + Exp_Rep_n(q, ro, b, n))
   }
macro Grs_BornMayer(s1, s2, wqi, wqj, c, ro, b)
   {
      Grs_BornMayer(s1, s2, wqi, wqj, c, ro, b, .2)
   }
macro Grs_BornMayer(s1, s2, wqi, wqj, c, ro, b, rsm)
   {
      grs_interaction qi wqi qj wqj s1 s2 c =
            Exp_Rep(Sqrt( ((wqi) (wqj))^2 ), ro, b);
   }
macro Grs_AB(s1, s2, wqi, wqj, c, & ro, & a)
   {
      grs_interaction qi wqi qj wqj s1 s2 c =
         If(R < ro,
            a (R-ro)^2,
            0
         );
   }
macro Temperature_Regime { load temperature }
macro Simulated_Annealing_1(v)
   {
      seed
      verbose -1
      continue_after_convergence
      chi2_convergence_criteria = (v);
      process_times
   }
macro TE_0(t, ts)
   {
      temperature = (t);
         move_to_the_next_temperature_regardless_of_the_change_in_rwp
         on_best_rewind
      load temperature { ts }
   }
macro TE_1(t, ts, x)
   {
      TE_0(t, ts)
		TE_0(t, ts)
      TE_0(t, ts)
      TE_0(t, ts)
      TE_0(t, ts)
      temperature 1 use_best_values
      load temperature { 1 1 1 1 1 1 1 1 1 1 }
		load temperature { 1 1 1 1 1 1 1 1 1 1 }
      temperature 1.0001
         on_best_rewind
      temperature = 100 (t);
         save_values_as_best_after_randomization
         move_to_the_next_temperature_regardless_of_the_change_in_rwp
   }
macro Auto_T(t)
   {
      Simulated_Annealing_1(0.001)
		quick_refine 0.01
      randomize_on_errors
      TE_1(t, 1 1 1, 1)
   }
macro Structure_Solution_Weighting
   {
      weighting = 1 / Sqrt(Max(Yobs, 1));
   }
macro Weighting_Individual_Atoms(fwhm)
   {
      weighting =
         Exp(- 4 Ln(2) ((Lam/(2 Sin(X Deg .5)))-1.3)^2 / (.25 Cycle_Iter + fwhm)^2) / Sqrt(Max(Yobs, 0.01 Get(max_yobs)));
         recal_weighting_on_iter
   }
macro Keep_Atom_Within_Box(& size)
   {
      move_to x min = Constant(Val) - size / Get(a); max = Constant(Val) + size / Get(a);
      move_to y min = Constant(Val) - size / Get(b); max = Constant(Val) + size / Get(b);
      move_to z min = Constant(Val) - size / Get(c); max = Constant(Val) + size / Get(c);
   }
' ---------------------------------------
macro XY(path, calc_step)
   {
      xdd path##.xy
      x_calculation_step calc_step
   }
macro BRML(w) { xdd w##.brml }
macro XYE(path_ext) {  xdd path_ext xye_format }
macro SCR(path) { xdd_scr path##.scr }
macro SHELX_HKL4(path) { xdd_scr path##.hkl }
macro DAT(path) { xdd path##.dat }
macro XDD(path) { xdd path##.xdd }
macro SST(path) { xdd path##.xdd }
macro RAW(path) { xdd path##.raw }
macro RAW(path, range_num) { xdd path##.raw range range_num }
macro STR(sg) { str space_group sg }
macro STR(sg, pn) { str space_group sg phase_name pn }
' ---------------------------------------
macro Cubic(cv)
   {
      a cv b = Get(a); c = Get(a);
   }
macro Tetragonal(a_cv, c_cv)
   {
      a a_cv b = Get(a); c c_cv
   }
macro Hexagonal(a_cv, c_cv)
   {
      a a_cv b = Get(a); c c_cv ga 120
   }
macro Rhombohedral(a_cv, al_cv)
   {
      a a_cv b = Get(a); c = Get(a);
      al al_cv be = Get(al); ga = Get(al);
   }
macro Trigonal { Hexagonal }
' ---------------------------------------
macro MM_CrystalAxis_Display(mxc, myc, mzc)
   {
     local = Get(mlx) Get(a); : mxc
     local = Get(mly) Get(b); : myc
     local = Get(mlz) Get(c); : mzc
   }
macro MM_CrystalAxis_Refine(mxc, mxv, myc, myv, mzc, mzv, mlx_v, mly_v, mlz_v)
   {
     #m_argu mxc
     #m_argu myc
     #m_argu mzc
     If_Prm_Eqn_Rpt(mxc, mxv, )
     If_Prm_Eqn_Rpt(myc, myv, )
     If_Prm_Eqn_Rpt(mzc, mzv, )

     mlx = CeV(mxc, mxv) / Get(a); : mlx_v
     mly = CeV(myc, myv) / Get(b); : mly_v
     mlz = CeV(mzc, mzv) / Get(c); : mlz_v
   }
macro MM_Cartesian_Display(mxc, myc, mzc)
   {
      local #m_unique !ax = Constant(Get(a));
      local #m_unique !bx = Constant(Get(b) Cos(Get(ga) Deg));
      local #m_unique !by = Constant(Get(b) Sin(Get(ga) Deg));
      local #m_unique !cx = Constant(Get(c) Cos(Get(be) Deg));
      local #m_unique !cy = Constant(Get(c) (Cos(Get(al) Deg) - Cos(Get(be) Deg) Cos(Get(ga) Deg)) / Sin(Get(ga) Deg));
      local #m_unique !cz = Constant(Sqrt(Get(c)^2 - cx^2 - cy^2));

      local = Get(mlx) ax + Get(mly) bx + Get(mlz) cx; : mxc
      local = Get(mly) by + Get(mlz) cy; : myc
      local = Get(mlz) cz; : mzc
   }
macro MM_Cartesian_Refine(mxc, mxv, myc, myv, mzc, mzv, mlx_v, mly_v, mlz_v)
   {
      #m_argu mxc
      #m_argu myc
      #m_argu mzc
      If_Prm_Eqn_Rpt(mxc, mxv, )
      If_Prm_Eqn_Rpt(myc, myv, )
      If_Prm_Eqn_Rpt(mzc, mzv, )

      local #m_unique !ax = Constant(Get(a));
      local #m_unique !bx = Constant(Get(b) Cos(Get(ga) Deg));
      local #m_unique !by = Constant(Get(b) Sin(Get(ga) Deg));
      local #m_unique !cx = Constant(Get(c) Cos(Get(be) Deg));
      local #m_unique !cy = Constant(Get(c) (Cos(Get(al) Deg) - Cos(Get(be) Deg) Cos(Get(ga) Deg)) / Sin(Get(ga) Deg));
      local #m_unique !cz = Constant(Sqrt(Get(c)^2 - cx^2 - cy^2));

      mlz =  CeV(mzc, mzv) / cz; : mlz_v
      mly = (CeV(myc, myv) - If(Abs(cy) > 1.0e-10, Get(mlz) cy, 0)) / by; : mly_v
      mlx = (CeV(mxc, mxv) - If(Abs(by) > 1.0e-10, Get(mly) by, 0) - If(Abs(cx) > 1.0e-10, Get(mlz) cx, 0)) / ax; : mlx_v
   }
' ---------------------------------------
macro Sites { load site x y z occ beq }
macro Lam_recs { load la lo lh }
macro Exclude { load exclude }
macro & Multiplicities_Sum( & sum_eqn)
   {
      Sum( sum_eqn, Mi = 0, Mi < M, Mi = Mi+1) / M
   }
' ---------------------------------------
macro Out_cf_hkl(file)
   {
      out file
      Out_String("loop_\n")
         Out_String("   _refln_index_h\n")
         Out_String("   _refln_index_k\n")
         Out_String("   _refln_index_l\n")
         Out_String("   _refln_F_squared_meas\n")
         Out_String("   _refln_F_squared_sigma\n")
      phase_out file append load out_record out_fmt out_eqn
      {
          " %4.0f" = H;
          " %4.0f" = K;
          " %4.0f" = L;
          " %9g" = I_no_scale_pks / M;
          " %1g\n" = 1;
      }
   }
macro Out_re_im_hkl(file)
   {
      phase_out file load out_record out_fmt out_eqn
      {
          " %4.0f" = H;
          " %4.0f" = K;
          " %4.0f" = L;
          " %11.5f" = A01;
          " %11.5f\n" = B01;
      }
   }
macro F_Real_positive { (A01-B11) }
macro F_Real_negative { (A01+B11) }
macro F_Imaginary_positive { (A11+B01) }
macro F_Imaginary_negative { (A11-B01) }
macro F2_positive { (F_Real_positive^2 + F_Imaginary_positive^2) }
macro F2_negative { (F_Real_negative^2 + F_Imaginary_negative^2) }
macro F2_Merged { (A01^2 + B01^2 + A11^2 + B11^2) }
macro Out_F2_Details(file)
   {
      out file
         Out_String("/")
         Out_String("*")
         Out_String("H --K --L --M -F_Real +ve -F_Real -ve -F_Imag +ve -F_Imag -ve -----F2 +ve -----F2 -ve --F2_Merged")
         Out_String(" *")
         Out_String("/\n")
      phase_out file append load out_record out_fmt out_eqn
      {
          "%3.0f" = H;
         " %3.0f" = K;
         " %3.0f" = L;
         " %3.0f" = M;
         " %11.5f" = F_Real_positive;
         " %11.5f" = F_Real_negative;
         " %11.5f" = F_Imaginary_positive;
         " %11.5f" = F_Imaginary_negative;
         " %11.5f" = F2_positive M .5;
         " %11.5f" = F2_negative M .5;
         " %11.5f\n" = F2_Merged M;
      }
   }
macro Out_Single_Crystal_Details(file)
   {
      phase_out file load out_record out_fmt out_eqn
      {
          " %4.0f" = H;
          " %4.0f" = K;
          " %4.0f" = L;
          " %4.0f" = Mobs;
          " %4.0f" = M;
          " %11.4f" = A01;
          " %11.4f" = A11;
          " %11.4f" = B01;
          " %11.4f" = B11;
          ' I_no_scale_pks
          '     = Get(scale) Mobs (A01-B11)^2 + (B01+A11)^2; when ignore_differences_in_Friedel_pairs is NOT defined
          '     = Get(scale) Mobs (A01^2 + B01^2 + A11^2 + B11^2); when ignore_differences_in_Friedel_pairs IS defined
          ' If there are no scale_pks then:
          '   I_no_scale_pks = I_after_scale_pks = Ycalc
          " %11.4f" = I_no_scale_pks;
          " %11.4f" = I_after_scale_pks;
          " %11.4f" = Ycalc;
          " %11.4f" = Yobs;
          " %11.4f\n" = SigmaYobs;
      }
   }
macro Out_A01_A11_B01_B11(file)
   {
      out file
         Out_String("/")
         Out_String("*")
         Out_String("H --K --L --M --------A01 --------A11 --------B01 --------B11")
         Out_String(" *")
         Out_String("/\n")
      phase_out file append load out_record out_fmt out_eqn
      {
          "%3.0f" = H;
         " %3.0f" = K;
         " %3.0f" = L;
         " %3.0f" = M;
         " %11.5f" = A01;
         " %11.5f" = A11;
         " %11.5f" = B01;
         " %11.5f\n" = B11;
      }
   }
macro Out_Yobs_Ycalc_and_Difference(file)
   {
      xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f  " = Yobs;
          " %11.6f  " = Ycalc;
          " %11.6f\n" = Yobs-Ycalc;
      }
   }
macro Out_Ycalc(file)
   {
      xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f\n" = Ycalc;
      }
   }
macro Out_X_Yobs(file)
   {
      xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f  " = X;
          " %11.6f\n" = Yobs;
      }
   }
macro Out_X_Y_(file, y)
   {
      xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f  " = X;
          " %11.6f\n" = y;
      }
   }   
macro Out_X_Ycalc(file)
   {
      xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f  " = X;
			 " %11.6f\n" = Ycalc;
      }
   }
macro Out_X_Yobs_Ycalc(file)
   {
      xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f  " = X;
          " %11.6f  " = Yobs;
          " %11.6f\n" = Ycalc;
      }
   }
macro Out_X_Difference(file)
   {
      xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f  " = X;
          " %11.6f\n" = Yobs-Ycalc;
      }
   }
macro Create_2Th_Ip_file(file, wI)
   {
      phase_out file load out_record out_fmt out_eqn
      {
         " %11.6f"   = 2 Rad Th;
         " %11.6f\n" = wI;
      }
   }
macro Create_d_Ip_file(file, wI)
   {
      phase_out file load out_record out_fmt out_eqn
      {
         " %11.6f"   = Lam / (2 Sin(Th));
         " %11.6f\n" = wI;
      }
   }
macro Create_hklm_d_Th2_Ip_file(file, wI)
   {
      phase_out file load out_record out_fmt out_eqn
      {
          "%3.0f" = H;
         " %3.0f" = K;
         " %3.0f" = L;
         " %3.0f" = M;
         " %11.5f"   = D_spacing;
         " %11.5f"   = 2 Rad Th;
         " %11.5f\n" = wI;
      }
   }
macro Create_2Th_Ip_file(file) {   Create_2Th_Ip_file(file, I_no_scale_pks) }
macro Create_d_Ip_file(file) { Create_d_Ip_file(file, I_no_scale_pks) }
macro Create_2Th_IScaled_file(file) {   Create_2Th_Ip_file(file, I_after_scale_pks) }
macro Create_d_IScaled_file(file) { Create_d_Ip_file(file, I_after_scale_pks) }
macro Create_hklm_d_Th2_Ip_file(file) { Create_hklm_d_Th2_Ip_file(file, I_no_scale_pks) }
macro Create_hklm_d_Th2_IScaled_file(file) { Create_hklm_d_Th2_Ip_file(file, I_after_scale_pks) }
' ---------------------------------------
macro Out_XDD_Start_Step(file)
   {
      out file         
         Out(X1) Out_String(" ") 
         Out(Yobs_dx_at(X1)) 
         Out_String("\n") 
   }
macro Out_XDD_xy(file, y)
   {
      xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f  " = X;
          " %11.6f\n" = y;
      }
   }             
macro Out_XDD_SST(file)
   {
      Out_XDD_Start_Step(file)
      xdd_out file append out_record out_fmt " %.9g\n" out_eqn 
   }
macro Out_String(fmt)
   {
      out_record
         out_fmt fmt
   }
macro Out(eqn, fmt)
   {
      out_record
         out_eqn = eqn; 
         out_fmt fmt
   }
macro Out(eqn)
   {
      out_record
         out_eqn = eqn; 
         out_fmt "%.9g"
   }
macro Out(eqn, fmt, fmt_err)
   {
      Out(eqn, fmt)
         out_fmt_err fmt_err
   }
macro Out_file(file, eqn, fmt)
   {
      out file append
      Out(eqn, fmt)
   }
macro Out_file(file, eqn, fmt, fmt_err)
   {
      out file append
      Out(eqn, fmt)
         out_fmt_err fmt_err
   }
' ---------------------------------------
macro No_Th_Dependence { lam no_th_dependence   la 1 lo 0 lh 1 }
macro Load_Lam(root, lam_file, yminymax) { #include root##lam\##lam_file##.lam ymin_on_ymax yminymax }
macro Load_Lam(lam_file, yminymax) { Load_Lam(ROOT, lam_file, yminymax) }

macro CuKa1(yminymax) {   Load_Lam(cuka1, yminymax) }
macro CuKa2_analyt(yminymax) {   Load_Lam(cuka2_analyt, yminymax) }
macro CoKa3(yminymax) {   Load_Lam(coka3, yminymax) }
macro CoKa7_Holzer(yminymax) { Load_Lam(coka7_holzer, yminymax) }
macro CoKb6_Holzer(yminymax) { Load_Lam(cokb6_holzer, yminymax) }
macro CrKa7_Holzer(yminymax) { Load_Lam(crka7_holzer, yminymax) }
macro CrKb5_Holzer(yminymax) { Load_Lam(crkb5_holzer, yminymax) }
macro CuK1sharp(yminymax) { Load_Lam(cuk1sharp, yminymax) }
macro CuKa2(yminymax) {   Load_Lam(cuka2, yminymax) }
macro CuKa4_Holzer(yminymax) { Load_Lam(cuka4_holzer, yminymax) }
macro CuKa5(yminymax) {   Load_Lam(cuka5, yminymax) }
macro CuKa5_Berger(yminymax) { Load_Lam(cuka5_berger, yminymax) }
macro CuKb4_Holzer(yminymax) { Load_Lam(cukb4_holzer, yminymax) }
macro FeKa7_Holzer(yminymax) { Load_Lam(feka7_holzer, yminymax) }
macro FeKb4_Holzer(yminymax) { Load_Lam(fekb4_holzer, yminymax) }
macro MnKa7_Holzer(yminymax) { Load_Lam(mnka7_holzer, yminymax) }
macro MnKb5_Holzer(yminymax) { Load_Lam(mnkb5_holzer, yminymax) }
macro MoKa2(yminymax) {   Load_Lam(moka2, yminymax) }
macro NiKa5_Holzer(yminymax) { Load_Lam(nika5_holzer, yminymax) }
macro NiKb4_Holzer(yminymax) { Load_Lam(nikb4_holzer, yminymax) }
macro NoThDependence(yminymax) {   Load_Lam(nothdependence, yminymax) }
macro DeltaFunction
   {
      NoThDependence(0.0001)
      peak_type pv
         pv_lor 0
         pv_fwhm = 0.01 Peak_Calculation_Step;
   }
' ---------------------------------------
macro Point_for_site(site, x, y, z)
   {
      point_for_site site ux x uy y uz z
   }
macro Duplicate_Point(new_point, prev_point)
   {
      point_for_site new_point
         ux = Point(prev_point, rx);
         uy = Point(prev_point, ry);
         uz = Point(prev_point, rz);
   }
macro Duplicate_rotate_z(new_point, prev_point, angle_cv)
   {
      Duplicate_Point(new_point, prev_point)
      rotate angle_cv qz 1 operate_on_points new_point
   }
macro Triangle(s1, s2, s3, & r)
   {
        point_for_site s1 ux = -r Sqrt(0.75); uy = r .5;
        point_for_site s2 ux =  r Sqrt(0.75); uy = r .5;
        point_for_site s3 ux 0 uy = -r;
   }
macro Triangle(s0, s1, s2, s3, r)
   {
        point_for_site s0
        Triangle(s1, s2, s3, r)
   }
macro Translate(acv, bcv, ccv)
   {
      translate ta acv tb bcv tc ccv
   }
macro Translate(acv, bcv, ccv, ops)
   {
      Translate(acv, bcv, ccv)
         operate_on_points ops
   }
macro Translate_with_site_start_values(s0, xc, yc, zc)
   {
      Translate(xc 0, yc 0, zc 0)
         start_values_from_site s0
   }
macro Triangle(s0, s1, s2, s3, r, xc, yc, zc, cva, cvb, cvc)
   {
      rigid
         Triangle(s0, s1, s2, s3, r)
         Rotate_about_axies(cva, cvb, cvc)
         Translate_with_site_start_values(s0, xc, yc, zc)
   }
macro Tetrahedra(s0, s1, s2, s3, s4, r, xc, yc, zc, cva, cvb, cvc)
   {
      rigid
         point_for_site s0
         point_for_site s1 ux = -0.816496580927726 r; uy = 0.4714045207910317 r;  uz = -r/3;
         point_for_site s2 ux =  0.816496580927726 r; uy = 0.4714045207910317 r;  uz = -r/3;
         point_for_site s3 uy = -0.942809041582063 r; uz = -r/3;
         point_for_site s4 uz = r;
         Rotate_about_axies(cva, cvb, cvc)
         Translate_with_site_start_values(s0, xc, yc, zc)
   }
macro Rotate_about_axies(cva, cvb)
   {
      rotate cva qa 1
      rotate cvb qb 1
   }
macro Rotate_about_axies(cva, cvb, cvc)
   {
      Rotate_about_axies(cva, cvb)
      rotate cvc qc 1
   }
macro Octahedra(s0, s1, s2, s3, s4, s5, s6, r)
   {
      point_for_site s0 ux = 0; uy =  0; uz  =  0;
      point_for_site s1 ux = r; uy =  0; uz  =  0;
      point_for_site s2 ux =-r; uy =  0; uz  =  0;
      point_for_site s3 ux = 0; uy =  r; uz  =  0;
      point_for_site s4 ux = 0; uy = -r; uz  =  0;
      point_for_site s5 ux = 0; uy =  0; uz  =  r;
      point_for_site s6 ux = 0; uy =  0; uz  = -r;
   }
macro Octahedra(s0, s1, s2, s3, s4, s5, s6, r, xc, yc, zc, cva, cvb, cvc)
   {
      rigid
         Octahedra(s0, s1, s2, s3, s4, s5, s6, r)
         Rotate_about_axies(cva, cvb, cvc)
         Translate_with_site_start_values(s0, xc, yc, zc)
   }
macro Set_Length(s0, s1, r, xc, yc, zc, cva, cvb)
   {
      rigid
         Point_for_site(s0,  0,  0,  0)
         Point_for_site(s1,  0,  0,  r)
         Rotate_about_axies(cva, cvb)
         Translate_with_site_start_values(s0, xc, yc, zc)
   }
macro Set_Lengths(s0, s1, s2, r, xc, yc, zc,
   cva1, cvb1,
   cva2, cvb2)
   {
      Set_Length(s0, s1, r, xc, yc, zc, cva1, cvb1)
      Set_Length(s0, s2, r, xc, yc, zc, cva2, cvb2)
   }
macro Set_Lengths(s0, s1, s2, s3, r, xc, yc, zc,
   cva1, cvb1,
   cva2, cvb2,
   cva3, cvb3)
   {
      Set_Lengths(s0, s1, s2, r, xc, yc, zc,
         cva1, cvb1,
         cva2, cvb2)
      Set_Length(s0, s3, r, xc, yc, zc, cva3, cvb3)
   }
macro Set_Lengths(s0, s1, s2, s3, s4, r, xcv, ycv, zcv,
   cva1, cvb1,
   cva2, cvb2,
   cva3, cvb3,
   cva4, cvb4)
   {
      Set_Lengths(s0, s1, s2, s3, r, xcv, ycv, zcv,
      cva1, cvb1,
      cva2, cvb2,
      cva3, cvb3)
      Set_Length(s0, s4, r, xcv, ycv, zcv, cva4, cvb4)
   }
macro Point_ua(site_name) { Point_rx_to_ua(Point(site_name)) }         
macro Point_ub(site_name) { Point_ry_to_ub(Point(site_name)) }         
macro Point_uc(site_name) { Point_rz_to_uc(Point(site_name)) }         
macro Point(a) { Find_Child(Get(point_for_sites), a) }
macro Point(a, xyz) { Point_Current_Value(Get(Point(a), xyz)) }
macro Rotation_vector_from_points(a, b)
   {
      qx = Point(b, rx) - Point(a, rx);
      qy = Point(b, ry) - Point(a, ry);
      qz = Point(b, rz) - Point(a, rz);
   }
macro Translate_point_amount(a) { Translate_point_amount(a, ) }
macro Translate_point_amount(a, pm)
   {
      translate
         tx = pm Point(a, rx);
         ty = pm Point(a, ry);
         tz = pm Point(a, rz);
   }
macro Rotate_about_these_points(cv, a, b, ops)
   {
      Translate_point_amount(a, -) ops
      rotate cv Rotation_vector_from_points(a, b) ops
      Translate_point_amount(a) ops
   }
macro Rotate_about_points(cv, a, b)
   {
      Rotate_about_these_points(cv, a, b, operate_on_points Double_Quote !##a !##b * Double_Quote)
   }
macro Rotate_about_points(cv, a, b, pts)
   {
      Rotate_about_these_points(cv, a, b, operate_on_points pts)
   }
macro Hexagon_sitting_on_point_in_xy_plane(s1, s2, s3, s4, s5, s6, & a)
   {
      /*
         Centre at (0,0)

               s5
            s3    s1
            s4    s2
               s6
      */
      point_for_site s1 ux =  a Sqrt(3) .5; uy =  a .5;
      point_for_site s2 ux =  a Sqrt(3) .5; uy = -a .5;
      point_for_site s3 ux = -a Sqrt(3) .5; uy =  a .5;
      point_for_site s4 ux = -a Sqrt(3) .5; uy = -a .5;
      point_for_site s5 uy =  a;
      point_for_site s6 uy = -a;
   }
macro Hexagon_sitting_on_side_in_xy_plane(s1, s2, s3, s4, s5, s6, & a)
   {
      /*
         Centre at (0,0)

               s2 s1
             s6     s5
               s4 s3
      */
      point_for_site s1 ux =  a .5; uy =  a Sqrt(3) .5;
      point_for_site s2 ux = -a .5; uy =  a Sqrt(3) .5;
      point_for_site s3 ux =  a .5; uy = -a Sqrt(3) .5;
      point_for_site s4 ux = -a .5; uy = -a Sqrt(3) .5;
      point_for_site s5 ux =  a;
      point_for_site s6 ux = -a;
   }
macro Spherical_Coordinates_Theta_Phi(thc, thv, phic, phiv, ops)
   {
      #m_argu thc
      #m_argu phic
      rotate thc thv qc 1 ops
      rotate phic phiv qa = Sin(Deg CeV(thc, thv)); qb = -Cos(Deg CeV(thc, thv)); ops
   }
macro TI_(s, c)
   {
      If(Prm_There(s##_##c), s##_##c^2, 0)
   }
macro Texture_Index(s, t)
   {
      local = {
         return 1
         +1/5  (TI_(s,c20)+TI_(s,c21p)+TI_(s,c21m)+TI_(s,c22p)+TI_(s,c22m))
         +1/9  (TI_(s,c41)+TI_(s,c41p)+TI_(s,c41m)+TI_(s,c42p)+TI_(s,c42m)+TI_(s,c43p)+TI_(s,c43m)+TI_(s,c44p)+TI_(s,c44m))
         +1/13 (TI_(s,c60)+TI_(s,c61p)+TI_(s,c61m)+TI_(s,c62p)+TI_(s,c62m)+TI_(s,c63p)+TI_(s,c63m)+TI_(s,c64p)+TI_(s,c64m)+TI_(s,c65p)+TI_(s,c65m)+TI_(s,c66p)+TI_(s,c66m))
         +1/17 (TI_(s,c80)+TI_(s,c81p)+TI_(s,c81m)+TI_(s,c82p)+TI_(s,c82m)+TI_(s,c83p)+TI_(s,c83m)+TI_(s,c84p)+TI_(s,c84m)+TI_(s,c85p)+TI_(s,c85m)+TI_(s,c86p)+TI_(s,c86m)+TI_(s,c87p)+TI_(s,c87m)+TI_(s,c88p)+TI_(s,c88m));
      } : t
   }
' ---------------------------------------
macro & Gauss(& xo, & fwhm)
   {
      (2 Sqrt(Ln(2) / Pi) / fwhm) Exp(-4 Ln(2)((X - xo) / fwhm)^2)
   }
macro & Lorentzian(& xo, & fwhm)
   {
      (2 / (Pi fwhm)) / (1 + 4 ((X-xo) / fwhm)^2)
   }
macro PV(& a, & xo, & fwhm, & g)
   {
      fit_obj = a (1-g) Lorentzian(xo, fwhm) + a g Gauss(xo, fwhm);
   }
macro Plot_Fit_Obj(fit_obj_scale, name)
   {
      dummy_str
			space_group 1
         phase_name name
         scale = fit_obj_scale;
   }
macro Plot_Fit_Obj(name)
   {
      dummy_str
			space_group 1
         phase_name name
   }
macro PV(a, xo, fwhm, g, av, xov, fwhmv, gv)
   {
      #m_argu a
      #m_argu xo
      #m_argu g
      #m_argu fwhm
      PV(CeV(a,av), CeV(xo, xov), CeV(fwhm, fwhmv), CeV(g,gv))
         min_X = CeV(xo, xov) - 20 CeV(fwhm,fwhmv);
         max_X = CeV(xo, xov) + 20 CeV(fwhm,fwhmv);
      If_Prm_Eqn_Rpt(a, av,  min 0.00001)
      If_Prm_Eqn_Rpt(xo, xov, )
      If_Prm_Eqn_Rpt(fwhm, fwhmv,  min .00001)
      If_Prm_Eqn_Rpt(g, gv,  min 0 max 1)
   }
macro PV_Left_Right(& a, & xo, & fwhm1, & fwhm2, & g)
   {
      PV(2 a/(1+fwhm2/fwhm1), xo, fwhm1, g)
         min_X = xo - 20 fwhm1;
         max_X <= xo;
      PV(2 a/(1+fwhm1/fwhm2), xo, fwhm2, g)
         min_X > xo;
         max_X = xo + 20 fwhm2;
   }
macro PV_Left_Right(a, xo, fwhm1, fwhm2, g, av, xov, fwhm1v, fwhm2v, gv)
   {
      prm  a   av   min 0.00001
      prm  xo  xov
      prm  g   gv   min 0 max 1
      prm  fwhm1  fwhm1v  min 0.00001
      prm  fwhm2  fwhm2v  min 0.00001
      PV_Left_Right(a, xo, fwhm1, fwhm2, g)
   }
macro Bkg_Diffuse(b, bv, bb, bbv)
   {
      #m_argu b
      #m_argu bb
      If_Prm_Eqn_Rpt(b, bv, )
      If_Prm_Eqn_Rpt(bb, bbv, )
      fit_obj = b Sin(bb 4 Pi Sin(X Deg_on_2) / Lam) / (4 Pi Sin(X Deg_on_2) / Lam);
   }
' ---------------------------------------
macro Radius(r) { Radius(r,r) }
macro Radius(rp, rs) { Rp rp Rs rs }
macro Specimen_Tilt(v) { Specimen_Tilt(,v)  }
macro Specimen_Tilt(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min 0.0001 max 2)
      hat = Cos(Th) Rad CeV(c,v) / Rs;
   }
macro Source_Width(v) { Source_Width(,v) }
macro Source_Width(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min 0.00001 max 2)
      hat = Rad CeV(c,v) / Rs;
   }
macro Slit_Width { SW }
macro SW(v) { SW(,v) }
macro SW(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min 0.00001 max 2)
      hat = Rad CeV(c,v) / Rs;
   }
macro Sample_Thickness(v) { Sample_Thickness(,v) }
macro Sample_Thickness(dc, dv)
   {
      #m_argu dc
      If_Prm_Eqn_Rpt(dc, dv, min 0.00001 max 10,  del = CeV(dc,dv) .05+.001;) ' d in mm
      exp_limit = -CeV(dc,dv) Cos(Th) 2 Rad / Rs;
   }
macro Absorption { AB }
macro AB(v) { AB(,v) }
macro AB(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min 3 max 500)
      exp_conv_const = Ln(0.001) Sin(2 Th) / ( CeV(c,v) Rs 0.003490658 );
      ' Note: 0.003490658 = 2 Deg / 10
   }
macro Absorption_With_Sample_Thickness_mm_Shape(u, uv, d, dv)
   {
      #m_argu u
      #m_argu d
      #m_ifarg u ""
         #m_unique_not_refine u
      #m_endif
      AB(u, uv)
      Sample_Thickness(d, dv)
   }
macro Absorption_With_Sample_Thickness_mm_Intensity(u, uv, d, dv)
   {
      #m_argu u
      #m_argu d
      #m_ifarg u ""
         #m_unique_not_refine u
      #m_endif
      scale_pks = (1 - Exp(-.2 CeV(u, uv) CeV(d, dv) / Sin(Th)));
   }
macro Absorption_With_Sample_Thickness_mm_Shape_Intensity(u, uv, d, dv)
   {
      #m_argu u
      #m_argu d
      #m_ifarg u ""
         #m_unique_not_refine u
      #m_endif
      Absorption_With_Sample_Thickness_mm_Shape(u, uv, d, dv)
      scale_pks = (1 - Exp(-.2 CeV(u, uv) CeV(d, dv) / Sin(Th)));
   }
macro Cylindrical_I_Correction(uRv) { Cylindrical_I_Correction(,uRv) }
macro Cylindrical_I_Correction(uRc, uRv)
   {
      #m_argu uRc
      If_Prm_Eqn_Rpt(uRc, uRv, min 0.0001 max 12,)
      scale_pks = AL_Cyl_Corr(CeV(uRc, uRv)) Cos(Th)^2 + AB_Cyl_Corr(CeV(uRc, uRv)) Sin(Th)^2;
   }
macro Cylindrical_2Th_Correction(uRv) { Cylindrical_2Th_Correction(,uRv) }
macro Cylindrical_2Th_Correction(uRc, uRv)
   {
      #m_argu uRc
      If_Prm_Eqn_Rpt(uRc, uRv, min -12 max 12)
      th2_offset =
         0.000033 CeV(uRc, uRv) (Th Rad)^(1.168 - 0.22 CeV(uRc, uRv) + 0.0168 CeV(uRc, uRv)^2)
         (90 - Th Rad)^(1.155 + 0.2054 CeV(uRc, uRv) - 0.0224 CeV(uRc, uRv)^2);
   }
   
macro Capillary_Offset_Sin_2Th_mm(v) { Capillary_Offset_Sin_2Th_mm(, v) }
macro Capillary_Offset_Sin_2Th_mm(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min -10 max 10)
      th2_offset = CeV(c, v) Rad Sin(2 Th) / Rs;
   }
macro Capillary_Offset_Cos_2Th_mm(v) { Capillary_Offset_Cos_2Th_mm(, v) }
macro Capillary_Offset_Cos_2Th_mm(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min -10 max 10)
      th2_offset = CeV(c, v) Rad Cos(2 Th) / Rs;
   }
macro Divergence(v) { Divergence(,v) }
macro Divergence(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min 0.0001 max 6)
      one_on_x_conv = -CeV(c,v)^2 Deg_on_2 / Tan(Th);
   }
macro Divergence_Sample_Length(c, v, slc, slv)
   {
      #m_argu c
      #m_argu slc
      If_Prm_Eqn_Rpt(c, v, min 0.0001 max = 2 Val + .1;, del = 0.01 Val + .001;)
      If_Prm_Eqn_Rpt(slc, slv, min 0.0001 max = 2 Val + 1;, del = 0.01 Val + .01;)
      one_on_x_conv = - Min(CeV(c,v), Rad Sin(Th) CeV(slc, slv)/Rs)^2  Deg_on_2 / Tan(Th);
      scale_pks = Min(CeV(c,v), Rad Sin(Th) CeV(slc, slv)/Rs) / CeV(c,v);
   }
macro FDS_BS_Shape(c, v, slc, slv)
   {
      #m_argu c
      #m_argu slc
      If_Prm_Eqn_Rpt(c, v, min 0.0001 max = 2 Val + .1;, del = 0.01 Val + .001;)
      If_Prm_Eqn_Rpt(slc, slv, min 0.0001 max = 2 Val + 1;, del = 0.01 Val + .01;)
      one_on_x_conv = - Min(CeV(c,v), Rad Sin(Th) CeV(slc, slv)/Rs)^2  Deg_on_2 / Tan(Th);
   }
macro Variable_Divergence_Shape(v) { Variable_Divergence_Shape(,v) }
macro Variable_Divergence_Shape(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min 0.0001 max 60)
      one_on_x_conv = -CeV(c,v)^2 Sin(2 Th) Rad / (4 Rs^2);
   }
macro Variable_Divergence_Intensity { scale_pks = Sin(Th); }
macro Variable_Divergence(v) { Variable_Divergence(,v) }
macro Variable_Divergence(c, v)
   {
      Variable_Divergence_Shape(c, v)
      Variable_Divergence_Intensity
   }
macro Simple_Axial_Model(v) { Simple_Axial_Model(,v) }
macro Simple_Axial_Model(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min 0.0001 max 50)
      circles_conv = -.5 Rad ( CeV(c,v) / Rs)^2 / Tan(2 Th);
   }
macro Full_Axial_Model(filament_cv, sample_cv, detector_cv)
   {
      axial_conv
         filament_length           filament_cv
         sample_length             sample_cv
         receiving_slit_length     detector_cv
   }
macro Full_Axial_Model(filament_cv, sample_cv, detector_cv, psol_cv, ssol_cv)
   {
      Full_Axial_Model(filament_cv, sample_cv, detector_cv)
      Sollers(psol_cv, ssol_cv)
   }
macro Sollers(cv_p, cv_s)
   {
      primary_soller_angle    cv_p
      secondary_soller_angle  cv_s
   }
macro Finger_et_al(s2, h2)
   {
      axial_conv
         filament_length            s2
         sample_length              = Get(source_len);
         receiving_slit_length      h2
         primary_soller_angle   0
   }
macro Tube_Tails(
   source_width_c, source_width_v, ' mm
   z1_c, z1_v, ' mm
   z2_c, z2_v, ' mm
   z1_z2_h_c, z1_z2_h_v ' fraction
   )
   {
      #m_argu source_width_c
      #m_argu z1_c
      #m_argu z2_c
      #m_argu z1_z2_h_c
      If_Prm_Eqn_Rpt(source_width_c, source_width_v, min 0.00001 max .5, del 0.005)
      If_Prm_Eqn_Rpt(z1_z2_h_c, z1_z2_h_v, min 0.00001 max .5, del .0001 )
      If_Prm_Eqn_Rpt(z1_c, z1_v, min -5 max -.00001, del -.005)
      If_Prm_Eqn_Rpt(z2_c, z2_v, min 0.00001 max 5, del .005 )
      stacked_hats_conv
         whole_hat = CeV(source_width_c, source_width_v) Rad /Rs;
         half_hat = CeV(z1_c, z1_v) Rad /Rs;
            hat_height = CeV(z1_z2_h_c, z1_z2_h_v);
         half_hat = CeV(z2_c, z2_v) Rad /Rs;
            hat_height = CeV(z1_z2_h_c, z1_z2_h_v);
   }           
' ---------------------------------------
macro I_Penalty(c, v)
   {
      ' http://topas.dur.ac.uk/unb/forum.php?req=thread&id=32
      #m_argu c
      c v
      = c^2;
   }
' ---------------------------------------
macro TCHZ_Peak_Type(u, v, w, z, x, y)
   {
      local #m_unique tch_p_l = x Tan(Th) + y / Cos(Th);
      local #m_unique tch_p_g = Sqrt( Abs( u Tan(Th)^2 + v Tan(Th) + w + z / Cos(Th)^2) );
      local #m_unique tch_p =
         (
                    tch_p_g^5 +
            2.69269 tch_p_g^4 tch_p_l +
            2.42843 tch_p_g^3 tch_p_l^2 +
            4.47163 tch_p_g^2 tch_p_l^3 +
            0.07842 tch_p_g   tch_p_l^4 +
                              tch_p_l^5
         )^0.2;
      local #m_unique tch_q = tch_p_l / tch_p;
      peak_type pv
         pv_lor = 1.36603 tch_q - 0.47719 tch_q^2 + 0.1116 tch_q^3;
         pv_fwhm = tch_p;
   }
macro TCHZ_Peak_Type(u, uv, v, vv, w, wv, z, zv, x, xv, y, yv)
   {
      #m_argu u
      #m_argu v
      #m_argu w
      #m_argu z
      #m_argu x
      #m_argu y
      If_Prm_Eqn_Rpt(u, uv, min = Max(-1, Val-.1); max = Min(2, Val+.1); del 1.0e-4)
      If_Prm_Eqn_Rpt(v, vv, min = Max(-1, Val-.1); max = Min(2, Val+.1); del 1.0e-4)
      If_Prm_Eqn_Rpt(w, wv, min = Max(-1, Val-.1); max = Min(2, Val+.1); del 1.0e-4)
      If_Prm_Eqn_Rpt(z, zv, min = Max(-1, Val-.1); max = Min(2, Val+.1); del 1.0e-4)
      If_Prm_Eqn_Rpt(x, xv, min = Max(0.0001, Val-.1); max = Min(2, Val+.1); del 1.0e-4 )
      If_Prm_Eqn_Rpt(y, yv, min = Max(0.0001, Val-.1); max = Min(2, Val+.1); del 1.0e-4 )
      TCHZ_Peak_Type(CeV(u, uv), CeV(v, vv), CeV(w, wv), CeV(z, zv), CeV(x, xv), CeV(y, yv))
   }
macro PV_Peak_Type(ha, hav, hb, hbv, hc, hcv, lora, lorav, lorb, lorbv, lorc, lorcv)
   {
      #m_argu ha
      #m_argu hb
      #m_argu hc
      #m_argu lora
      #m_argu lorb
      #m_argu lorc
      If_Prm_Eqn_Rpt(ha,   hav,   min 0.0001 max = 2 Val + .1;, del 0.001)
      If_Prm_Eqn_Rpt(hb,   hbv,   min 0.0001 max = 2 Val + .1;, del 0.001)
      If_Prm_Eqn_Rpt(hc,   hcv,   min 0.0001 max = 2 Val + .1;, del 0.001)
      If_Prm_Eqn_Rpt(lora, lorav, min 0.0001 max 1,            del 0.001)
      If_Prm_Eqn_Rpt(lorb, lorbv, min 0.0001 max = 2 Val + .1;, del 0.001)
      If_Prm_Eqn_Rpt(lorc, lorcv, min 0.0001 max = 2 Val + .1;, del 0.001)
      peak_type pv
         pv_lor = CeV(lora,lorav) + CeV(lorb,lorbv) Tan(Th) + CeV(lorc,lorcv) / Cos(Th);
         pv_fwhm = CeV(ha,hav) + CeV(hb,hbv) Tan(Th) + CeV(hc,hcv) / Cos(Th);
   }
macro PVII_Peak_Type(ha, hav, hb, hbv, hc, hcv, ma, mav, mb, mbv, mc, mcv)
   {
      #m_argu ha
      #m_argu hb
      #m_argu hc
      #m_argu ma
      #m_argu mb
      #m_argu mc
      If_Prm_Eqn_Rpt(ha, hav, min 0.0001 max 1)
      If_Prm_Eqn_Rpt(hb, hbv, min 0.0001 max 1)
      If_Prm_Eqn_Rpt(hc, hcv, min 0.0001 max 1)
      If_Prm_Eqn_Rpt(ma, mav, min 0.0001 max 20)
      If_Prm_Eqn_Rpt(mb, mbv, min 0.0001 max 5)
      If_Prm_Eqn_Rpt(mc, mcv, min 0.0001 max 5)
      peak_type spvii
         m1 = .6 + CeV(ma,mav) + CeV(mb,mbv) / Cos(Th) + CeV(mc,mcv) Tan(Th);
         m2 = Get(m1);
         h1 = CeV(ha,hav) + CeV(hb,hbv) / Cos(Th) + CeV(hc,hcv) Tan(Th);
         h2 = Get(h1);
   }
macro UVW(u, uv, v, vv, w, wv)
   {
      #m_argu u
      #m_argu v
      #m_argu w
      If_Prm_Eqn_Rpt(u, uv, min -1 max 2)
      If_Prm_Eqn_Rpt(v, vv, min -1 max 1)
      If_Prm_Eqn_Rpt(w, wv, min -1 max 2)
      gauss_fwhm = CeV(u, uv) Tan(Th)^2 + CeV(v, vv) Tan(Th) + CeV(w, wv);
   }
macro Sec_Th_Gaussian_Convolution { CS_G }
macro CS_G(v) { CS_G(,v) }
macro CS_G(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min .3 max = Min(Val 2 + .3, 10000);)
      gauss_fwhm = 0.1 Rad Lam / (Cos(Th) CeV(c,v));
   }
macro CS { CS_L }
macro Crystallite_Size { CS_L }
macro CS_L(v) { CS_L(,v) }
macro CS_L(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min .3 max = Min(Val 2 + .3, 10000); )
		lor_fwhm = 0.1 Rad Lam / (Cos(Th) CeV(c,v));
   }     
macro MS { Strain_L }
macro Microstrain { Strain_L }
macro Strain_L(v) { Strain_L(,v) }
macro Strain_L(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min .0001 max 5)
      lor_fwhm = CeV(c,v) Tan(Th);
   }
macro Strain_G(v) { Strain_G(,v) }
macro Strain_G(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min .0001 max 5)
      gauss_fwhm = CeV(c,v) Tan(Th);
   }
macro Strain(v) {   Strain(,v) }
macro Strain(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min .0001 max 2)
      gauss_fwhm = 2.69129 CeV(c,v) Tan(Th);
      ' Note: 2.69129 = Sqrt(4 Ln(2)/Pi) (5 / 100) (180 / Pi)
   }
macro User_Defined_Dependence_Convolution(convtype, x_dep, c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min .0001 max 100)
      #m_ifarg x_dep ""
         convtype = CeV(c,v);
      #m_else
         convtype = CeV(c,v) ( x_dep );
      #m_endif
   }
macro One_on_X(v) { One_on_X(,v) }
macro One_on_X(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min 0.001 )
      fit_obj = CeV(c,v) / X;
   }			  
macro LP_Factor_X(v) { LP_Factor_X(,v) }
macro LP_Factor_X(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, del 0.01 min 0 max 90)
      local #m_unique th = X Pi / 360;
      scale_phase_X = (1 + Cos(c Deg)^2 Cos(2 th)^2) / (Sin(th)^2 Cos(th));
   }
macro LP_Factor(v) { LP_Factor(,v) }
macro LP_Factor(c, v)
   {
      #m_argu c
		If_Prm_Eqn_Rpt(c, v, del 0.01 min 0 max 90)
      scale_pks = (1 + Cos(CeV(c,v) Deg)^2 Cos(2 Th)^2) / (Sin(Th)^2 Cos(Th));
   }     
macro Lorentz_Factor
   {
      scale_pks = 1 / (Sin(Th)^2 Cos(Th));
   }
macro LP_Factor_Synchrotron_Simple
   {
      Lorentz_Factor
   }
macro LP_Factor_Synchrotron(pp, ppv, mono, monov)
   /*
        By Ian Madsen, CSIRO Minerals, Australia
        pp is the polarisation in the plane of the synchrotron
        pp = 1.0 for circularly polarised X-rays (i.e. laboratory X-ray tubes)
        pp = 0.0 for fully polarised X-rays (ideal synchrotron source)
        expect pp ~ 0.05 for 'real' synchrotron source
        mono = the 2theta diffraction angle of a crystal monochromator
        NOTE - pp and mono will be ~100% correlated - do not attempt to refine both together!
        Last modification - 12/02/2008
   */
   {
      #m_argu pp
      #m_argu mono
      If_Prm_Eqn_Rpt(pp, ppv, min 0.0000001 max 1.0)
      If_Prm_Eqn_Rpt(mono, monov, min 0.0000001 max 90.0)
      scale_pks = (1/( Sin(Th)^2 Cos(Th))) * ((1 + CeV(pp,ppv) * Cos(CeV(mono,monov) Deg)^2 Cos(2 Th)^2)/(1 + CeV(pp,ppv) * Cos(CeV(mono,monov) Deg)^2 ));
   }
macro Surface_Roughness_Pitschke_et_al(a1c, a1v, a2c, a2v)
   {
      #m_argu a1c
      #m_argu a2c
      If_Prm_Eqn_Rpt(a1c, a1v, min .0001 max 1)
      If_Prm_Eqn_Rpt(a2c, a2v, min .0001 max 1)
      scale_pks = (1 - CeV(a1c, a1v) (1 / Sin(Th) - CeV(a2c, a2v) / Sin(Th)^2)) / (1 - CeV(a1c, a1v) + CeV(a1c, a1v) CeV(a2c, a2v));
   }
macro Surface_Roughness_Suortti(a1c, a1v, a2c, a2v)
   {
      #m_argu a1c
      #m_argu a2c
      If_Prm_Eqn_Rpt(a1c, a1v, min .0001 max 2)
      If_Prm_Eqn_Rpt(a2c, a2v, min .0001 max 2)
      scale_pks = (CeV(a1c, a1v) + (1 - CeV(a1c, a1v)) Exp( -CeV(a2c, a2v) /Sin(Th))) / (CeV(a1c, a1v) + (1 - CeV(a1c, a1v)) Exp( -CeV(a2c, a2v)));
   }
macro Zero_Error_No_Reset(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v,
         min = Max(Val - 20 Yobs_dx_at(X1), -100 Yobs_dx_at(X1));
         max = Min(Val + 20 Yobs_dx_at(X2),  100 Yobs_dx_at(X2));
         del = .01 Yobs_dx_at(X1);
         )
      th2_offset = CeV(c,v);
   }
macro Zero_Error { ZE }
macro ZE(v) { ZE(,v) }
macro ZE(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v,
         min = Max(Val - 20 Yobs_dx_at(X1), -100 Yobs_dx_at(X1));
         max = Min(Val + 20 Yobs_dx_at(X2),  100 Yobs_dx_at(X2));
         del = .01 Yobs_dx_at(X1);
         _v 0
         )
      th2_offset = CeV(c,v);
   }
macro Specimen_Displacement { SD }
macro SD(v)  { SD(, v) }
macro SD(c, v)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min = Val-.1; max = Val+.1; del 0.001)
      th2_offset = -2 Rad CeV(c,v) Cos(Th) / Rs;
   }
macro & PO_eqn(& c, & ang)
   {
      (c^2 Cos(ang)^2 + Sin(ang)^2 / c) ^ (-1.5)
   }
macro Preferred_Orientation { PO }
macro PO(c, v, ang, hkl)
   {
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min 0.0001 max = 2 Val + .5;)
      PO(CeV(c, v), ang, hkl)
   }
macro PO(& c, md, hkl)
   {
      #m_ifarg md ""
         #m_unique_not_refine md
      #m_endif
      march_dollase md hkl = c;
      scale_pks = md;
   }
macro PO_Two_Directions(
   c1, v1, ang1, hkl1,
   c2, v2, ang2, hkl2,
   w1c, w1v
   )
   {
      #m_ifarg ang1 ""
         #m_unique_not_refine ang1
      #m_endif
      #m_ifarg ang2 ""
         #m_unique_not_refine ang2
      #m_endif
      #m_argu c1
      #m_argu c2
      #m_argu w1c
      If_Prm_Eqn_Rpt(c1, v1, min 0.0001 max = 2 Val + .5;)
      If_Prm_Eqn_Rpt(c2, v2, min 0.0001 max = 2 Val + .5;)
      If_Prm_Eqn_Rpt(w1c, w1v, min 0.0001 max 1)
      march_dollase ang1 hkl1 = CeV(c1,v1);
      march_dollase ang2 hkl2 = CeV(c2,v2);
      scale_pks = CeV(w1c, w1v) ang1 + (1 - CeV(w1c, w1v)) ang2;
   }
macro PO_Spherical_Harmonics(sh, order)
   {
      #m_ifarg sh ""
         #m_unique_not_refine sh
      #m_endif
      spherical_harmonics_hkl sh
         sh_order order
      scale_pks = sh;
   }
macro ADPs { load u11 u22 u33 u12 u13 u23 }
macro ADPs_Keep_PD
   {
      move_to u12 min 0 max = Sqrt(Get(u11) Get(u22));
      move_to u13 min 0 max = Sqrt(Get(u11) Get(u33));
      move_to u23 min 0 max = Sqrt(Get(u22) Get(u33));
   }
' ---------------------------------------
macro Remove_Phase(& wt, & fr)
   {
      remove_phase = And(Get(weight_percent) < wt, Error(Get(weight_percent)) > fr Get(weight_percent));
   }
' ---------------------------------------
macro Robust_Refinement
   {
      /*
      ' Rescale peaks according to robust refinement algorithm
      ' K. H. Stone, S. H. Lapidus and P. W. Stephens; J. Appl. Cryst. (2009). 42, 385-391
      */
      weighting =
         {
            def test = Get(r_exp);
            def N = 1 / test^2;
            def p0 = 0.40007404;
            def p1 = -2.5949286;
            def p2 = 4.3513542;
            def p3 = -1.7400101;
            def p4 = 3.6140845e-1;
            def p5 = -4.45247609e-2;
            def p6 = 3.5986364e-3;
            def p7 = -1.8328008e-4;
            def p8 = 5.7937184e-6;
            def p9 = -1.035303e-7;
            def p10 = 7.9903166e-10;
            def t = (Yobs - Ycalc) / SigmaYobs;
            return
               If(t < 0.8,
                  N / Max(SigmaYobs^2, 1),
                  If(t < 21,
                     N ((((((((((p10 t + p9) t + p8) t + p7) t + p6) t + p5) t + p4) t + p3) t + p2) t + p1) t + p0) / (Yobs - Ycalc)^2,
                     N (2.0131 Ln(t) + 3.9183) / (Yobs - Ycalc)^2
                  )
               );
         }
		recal_weighting_on_iter
   }
' ---------------------------------------
macro Sinc_Convolution(k, kv, x_range)
   {
      #m_argu k
      If_Prm_Eqn_Rpt(k, kv, min 0.001 max 100)
      Sinc_Convolution(CV(k, kv), x_range)
   }
macro Sinc_Convolution(& k, & x_range)
   {
      user_defined_convolution = If(Abs(X) < 1e-5, 1, (Sin(k X) /(k X))^2);
         min = -x_range; max = x_range;
   }
' ---------------------------------------
macro Known_Weight_Percent(& w)
   {
      scale = (w/(100-w)) Get(sum_smvs_minus_this) / (Get(cell_mass) Get(cell_volume));
   }     
macro MVW(m_v, v_v, w_v)
   {
      cell_mass       m_v
      cell_volume     v_v
      weight_percent  w_v
   }
macro Apply_Brindley_Spherical_R_PD(& R, & PD)
   {
      brindley_spherical_r_cm = R PD;
   }
macro Mixture_LAC_1_on_cm(mlac)
   {
      local = Get(mixture_MAC) Get(mixture_density_g_on_cm3); : mlac
   }
macro Phase_Density_g_on_cm3(pd)
   {
      local = Phase_Density_Eqn_g_on_cm3; : pd
   }
macro Phase_LAC_1_on_cm(u)
   {
      local = Phase_LAC_Eqn_1_on_cm; : u
   }
macro Phase_LAC_Eqn_1_on_cm
   {
      (Get(phase_MAC) Phase_Density_Eqn_g_on_cm3)
   }
macro & Phase_Density_Eqn_g_on_cm3
   {
      1.6605402 Get(cell_mass) / Get(cell_volume)
   }
' ---------------------------------------
macro K_Factor_MAC_K(mac, k, amor)
   {
      move_to xdd
      local !k_factor_mac_local_ mac
      local !k_factor_k_local_ k
      local !k_factor_sum_wps_  100 : amor
   }
macro K_Factor_WP(result)
   {
      local k_factor_wp_ = 1.6605402 Get(smv) k_factor_mac_local_  / k_factor_k_local_; : result
      if Prm_There(k_factor_sum_wps_) {
         existing_prm k_factor_sum_wps_ -= k_factor_wp_;
      }
   }
' ----------------w   -----------------------
macro Out_FCF(file)
   {
      out file
         Out_String("\ndata_")
         Out(Get(a), "\n_cell_length_a %V")
         Out(Get(b), "\n_cell_length_b %V")
         Out(Get(c), "\n_cell_length_c %V")
         Out(Get(al), "\n_cell_angle_alpha %V")
         Out(Get(be), "\n_cell_angle_beta  %V")
         Out(Get(ga), "\n_cell_angle_gamma %V")
         Out_String("\n_shelx_F_squared_multiplier 1")

         Out_String("\nloop_\n_symmetry_equiv_pos_as_xyz")
         Out(Get(sp_xyzs_txt),  "%s")

         Out_String("\nloop_")
         Out_String("\n_refln_index_h")
         Out_String("\n_refln_index_k")
         Out_String("\n_refln_index_l")
         Out_String("\n_refln_F_squared_calc")
         Out_String("\n_refln_F_squared_meas")
         Out_String("\n_refln_F_squared_sigma")
         Out_String("\n_refln_observed_status")
         phase_out file append
            load out_record out_fmt out_eqn
            {
               "\n%4.0f" = H;
               " %4.0f" = K;
               " %4.0f" = L;
               " %12.2f" = I_no_scale_pks;
               " %12.2f" = Iobs_no_scale_pks;
               " %10.2f o" = Iobs_no_scale_pks_err;
            }
   }
macro Out_CIF_STR(file)
   {
      Out_CIF_STR(file,)
   }
macro Out_CIF_STR(file, with_id)
   {
      ' https://www.iucr.org/__data/iucr/cifdic_html/1/cif_core.dic/Iatom_site_symmetry_multiplicity.html 
      out file
         Out_String("\ndata_")
         Out(Get(phase_name), "\n_chemical_name_mineral ?%s?")
         Out(Get(a), "\n_cell_length_a  %V")
         Out(Get(b), "\n_cell_length_b  %V")
         Out(Get(c), "\n_cell_length_c  %V")
         Out(Get(al), "\n_cell_angle_alpha %V")
         Out(Get(be), "\n_cell_angle_beta  %V")
         Out(Get(ga), "\n_cell_angle_gamma %V")
         Out(Get(cell_volume), "\n_cell_volume %V")
         Out(Get(sp_grp_char), "\n_symmetry_space_group_name_H-M %s")
         Out_String("\nloop_")
         #m_ifarg with_id ""
            Out_String("\n\t_symmetry_equiv_pos_as_xyz")
            Out(Get(sp_xyzs_txt),  "%s")
         #m_else
            Out_String("\n\t_symmetry_equiv_pos_site_id\n\t_symmetry_equiv_pos_as_xyz")
            Out(Get(sp_xyzs_txt_with_id),  "%s")
         #m_endif
         Out_String("\nloop_")
            Out_String("\n_atom_site_label")
            Out_String("\n_atom_site_type_symbol")
            Out_String("\n_atom_site_site_symmetry_multiplicity")
            Out_String("\n_atom_site_fract_x")
            Out_String("\n_atom_site_fract_y")
            Out_String("\n_atom_site_fract_z")
            Out_String("\n_atom_site_occupancy")
            Out_String("\n_atom_site_B_iso_or_equiv")
            atom_out file append
               load out_record out_fmt out_eqn
               {
                  "\n%s" = Get_From_String(Get(current_atom), "site");
                  " %s" = Get_From_String(Get(current_atom), "atom");
                  " %3.0f" = Get_From_String(Get(current_atom), "num_posns");
                  " %V" = Get_From_String(Get(current_atom), "x");
                  " %V" = Get_From_String(Get(current_atom), "y");
                  " %V" = Get_From_String(Get(current_atom), "z");
                  " %V" = Get_From_String(Get(current_atom), "occ");
						" %V" = Get_From_String(Get(current_atom), "beq");
               }
   }
macro Out_CIF_ADPs(file)
   {
      out file append
         Out_String("\nloop_")
         Out_String("\n_atom_site_aniso_label")
         Out_String("\n_atom_site_type_symbol")
         Out_String("\n_atom_site_aniso_U_11")
         Out_String("\n_atom_site_aniso_U_22")
         Out_String("\n_atom_site_aniso_U_33")
         Out_String("\n_atom_site_aniso_U_12")
         Out_String("\n_atom_site_aniso_U_13")
         Out_String("\n_atom_site_aniso_U_23")
         atom_out file append adps
            load out_record out_fmt out_eqn
            {
               "\n%s" = Get_From_String(Get(current_atom), "site");
               " %s" = Get_From_String(Get(current_atom), "atom");
               " %V" = Get_From_String(Get(current_atom), "u11");
               " %V" = Get_From_String(Get(current_atom), "u22");
               " %V" = Get_From_String(Get(current_atom), "u33");
               " %V" = Get_From_String(Get(current_atom), "u12");
               " %V" = Get_From_String(Get(current_atom), "u13");
               " %V" = Get_From_String(Get(current_atom), "u23");
             }
   }
macro Out_CIF_Bonds_Angles(file)
   {
      out file append
         Out(Get(cif_bonds_angles), "%s")
   }
' ---------------------------------------
macro Unique_Cubic_sgs { try_space_groups "228 219 203 210 196 230 220 206 214 197 222 218 201 205 212 198 195" }
macro Unique_Trigonal_Hexagonal_sgs { try_space_groups "161 146 184 159 158 169 144 173 143" }
macro Unique_Tetragonal_sgs { try_space_groups "142 110 141 109 108 88 80 79 130 126 133 103 104 106 137 138 134 125 114 105 102 101 100 86 85 92 94 76 77 90 75" }
macro Unique_Orthorhombic_sgs { try_space_groups "70 43 22 68 73 37 45 41 46 36 39 20 23 21 52 56 60 61 48 54 50 33 34 32 30 29 27 31 26 19 18 17 16" }
macro Unique_Monoclinic_sgs { try_space_groups "9 5 14 7 4 3" }
macro Unique_Triclinic_sgs { try_space_groups "2" }

macro Bravais_Cubic_sgs { try_space_groups "196 197 195" }
macro Bravais_Trigonal_Hexagonal_sgs { try_space_groups "146 143" }
macro Bravais_Tetragonal_sgs { try_space_groups "79 75" }
macro Bravais_Orthorhombic_sgs { try_space_groups "22 23 21 16" }
macro Bravais_Monoclinic_sgs { try_space_groups "5 3" }
macro Bravais_Triclinic_sgs { try_space_groups "2" }

macro Indexing_Solutions_With_Zero_Error { load dummy ndx_sg ndx_uni ndx_vol ndx_gof ndx_ze ndx_alp ndx_blp ndx_clp ndx_allp ndx_belp ndx_galp }
macro Indexing_Solutions { load dummy ndx_sg ndx_uni ndx_vol ndx_gof ndx_alp ndx_blp ndx_clp ndx_allp ndx_belp ndx_galp }

macro Indexing_Solutions_With_Zero_Error_2 { load dummy ndx_sg ndx_status ndx_uni ndx_vol ndx_gof ndx_ze ndx_alp ndx_blp ndx_clp ndx_allp ndx_belp ndx_galp }
macro Indexing_Solutions_2 { load dummy ndx_sg ndx_status ndx_uni ndx_vol ndx_gof ndx_alp ndx_blp ndx_clp ndx_allp ndx_belp ndx_galp }

macro Index_x0_from_d(d, worder) {   index_x0 = 1 / ( (worder) d)^2; }
macro Index_x0_from_d(d) {   Index_x0_from_d(d, 1) }
macro Index_x0_from_th2(th2, worder) { index_x0 = ( 2 Sin(th2 Pi/360) / (Get(index_lam) (worder))  )^2; }
macro Index_x0_from_th2(th2) { Index_x0_from_th2(th2, 1) }

macro Cubic_F { try_space_groups "196" }
macro Cubic_I { try_space_groups "197" }
macro Cubic_P { try_space_groups "195" }
macro Trigonal_Hexagonal_R { try_space_groups "146" }
macro Trigonal_Hexagonal_P { try_space_groups "143" }
macro Tetragonal_I { try_space_groups "79" }
macro Tetragonal_P { try_space_groups "75" }
macro Orthorhombic_F { try_space_groups "22" }
macro Orthorhombic_I { try_space_groups "23" }
macro Orthorhombic_C { try_space_groups "21" }
macro Orthorhombic_P { try_space_groups "16" }
macro Monoclinic_C { try_space_groups "5" }
macro Monoclinic_P { try_space_groups "3" }
macro Triclinic_P { try_space_groups "2" }
' ---------------------------------------
macro Decompose(diff_toll)
   {
      x_calculation_step = diff_toll;
      yobs_to_xo_posn_yobs = Get(x_calculation_step);
   }
' ---------------------------------------
macro Out_for_cf(file)
   {
      I_parameter_names_have_hkl int
      out_A_matrix file
         A_matrix_prm_filter int*
   }
macro Tangent(& fraction_n, max_triplets, & max_num_e_values)
   {
      tangent_num_h_read = Get(num_observed_reflections_above_d_min);
      tangent_num_h_keep = fraction_n Get(tangent_num_h_read);
      tangent_num_k_read = Get(tangent_num_h_read);
      tangent_min_triplets_per_h = 1;
      tangent_max_triplets_per_h = max_triplets;
   }
macro Tangent(fraction_n, max_triplets)
   {
      Tangent(fraction_n, max_triplets, 35000)
   }
macro Restart_CF
   {
      use_Fc
      randomize_initial_phases_by 0
   }
macro & Sine_Wave(& x1, & x2, & n)
   {
      0.5 (x1 + x2) + 0.5 (x2 - x1) Sin(2 Pi Cycle_Iter / n - Pi / 2)
   }
macro & Sine_Wave(& x0, & x1, & x2, & n)
   {
      0.5 (x1 + x2) + 0.5 (x2 - x1) Sin(2 Pi (Cycle_Iter - x0) / n - Pi / 2)
   }
macro & Cosine_Wave(& x1, & x2, & n)
   {
      0.5 (x1 + x2) + 0.5 (x2 - x1) Cos(2 Pi Cycle_Iter / n)
   }
macro & Ramp_Run_Number(& x1, & x2, & n)
   {
      x1 + (x2-x1) Mod(Run_Number, n) / (n-1)
   }
macro & Ramp_0(& i, & x1, & x2, & n1, & n2)
   {
      x1 + (x2-x1) Mod(i, n2) / (n1-1)
   }
macro & Ramp(& x1, & x2, & n)
   {
      x1 + (x2-x1) Mod(Cycle_Iter, n) / (n-1)
   }
macro & Limit(& x, & x1, & x2)
   {
      Min(x2, Max(x, x1))
   }
macro Limit_with_sign(& x, & lim) 
	{
		If(Abs(x) < lim, If(x < 0, -lim, lim), x) 
	}
macro & Ramp_Clamp(& x1, & x2, & n1, & n2)
   {
      If(Mod(Cycle_Iter, n2) >= n1,
         x2,
         Ramp_0(Cycle_Iter, x1, x2, n1, n2)
      )
   }
macro & Ramp_Clamp(& x1, & x2, & n)
   {
      If(Mod(Cycle_Iter, n) >= n,
         x2,
         Ramp_0(Cycle_Iter, x1, x2, n, n)
      )
   }
macro & Cycle_Ramp(& x1, & x2, & n)
   {
      x1 + (x2-x1) Mod(Cycle, n) / (n-1)
   }
macro Pick(& n)
   {
      pick_atoms_when = Mod(Cycle_Iter, n) == (n-1);
   }
' ---------------------------------------
macro Bkg_Straight_Line_xy(n, x, y)
   {
      prm !bkg_x_##n = x;
      prm !bkg_y_##n = y;
   }
macro Bkg_Straight_Line_First(n, x, y)
   {
      Bkg_Straight_Line_xy(n, x, y)
   }
macro Bkg_Straight_Line(n1, n2, x2, y2)
   {
      fit_obj = (X - bkg_x_##n1) ((y2)-bkg_y_##n1)/((x2)-bkg_x_##n1) + bkg_y_##n1;
         #m_ifarg n1 "0"
            min_X >= bkg_x_##n1 - 1.0e-9;
         #m_else
            min_X >= bkg_x_##n1 + 1.0e-9;
         #m_endif
         max_X <= x2 + 1.0e-10;
      Bkg_Straight_Line_xy(n2, x2, y2)
   }
' ---------------------------------------
macro Stephens_lor_gauss(etac, etav, mhkl)
   {
      #m_argu etac  If_Prm_Eqn_Rpt(etac, etav, min 0 max 1)
      local #m_unique pp_ = D_spacing^2 * Sqrt(Max(mhkl,0)) Tan(Th) 0.0018/3.1415927;
		gauss_fwhm =  pp_ (1-CeV(etac, etav));
      lor_fwhm   =  pp_ CeV(etac, etav);
   }	  
macro Stephens_triclinic
   (
      etac, etav,
      s400c,s400v,
      s040c,s040v,
      s004c,s004v,
      s220c,s220v,
      s202c,s202v,
      s022c,s022v,
      s310c,s310v,
      s103c,s103v,
      s031c,s031v,
      s130c,s130v,
      s301c,s301v,
      s013c,s013v,
      s211c,s211v,
      s121c,s121v,
      s112c,s112v
   )
   {
      #m_argu s400c If_Prm_Eqn_Rpt(s400c,s400v,)
      #m_argu s040c If_Prm_Eqn_Rpt(s040c,s040v,)
      #m_argu s004c If_Prm_Eqn_Rpt(s004c,s004v,)
      #m_argu s220c If_Prm_Eqn_Rpt(s220c,s220v,)
      #m_argu s202c If_Prm_Eqn_Rpt(s202c,s202v,)
      #m_argu s022c If_Prm_Eqn_Rpt(s022c,s022v,)
      #m_argu s310c If_Prm_Eqn_Rpt(s310c,s310v,)
      #m_argu s103c If_Prm_Eqn_Rpt(s103c,s103v,)
      #m_argu s031c If_Prm_Eqn_Rpt(s031c,s031v,)
      #m_argu s130c If_Prm_Eqn_Rpt(s130c,s130v,)
      #m_argu s301c If_Prm_Eqn_Rpt(s301c,s301v,)
      #m_argu s013c If_Prm_Eqn_Rpt(s013c,s013v,)
      #m_argu s211c If_Prm_Eqn_Rpt(s211c,s211v,)
      #m_argu s121c If_Prm_Eqn_Rpt(s121c,s121v,)
      #m_argu s112c If_Prm_Eqn_Rpt(s112c,s112v,)
      Stephens_lor_gauss(
         etac, etav,
         CeV(s400c,s400v) H^4 +
         CeV(s040c,s040v) K^4 +
         CeV(s004c,s004v) L^4 +
         CeV(s220c,s220v) H^2 K^2 +
         CeV(s202c,s202v) H^2 L^2 +
         CeV(s022c,s022v) K^2 L^2 +
         CeV(s310c,s310v) H^3 K +
         CeV(s103c,s103v) H L^3 +
         CeV(s031c,s031v) K^3 L +
         CeV(s130c,s130v) H K^3 +
         CeV(s301c,s301v) H^3 L +
         CeV(s013c,s013v) K L^3 +
         CeV(s211c,s211v) H^2 K L +
         CeV(s121c,s121v) H K^2 L +
         CeV(s112c,s112v) H K L^2
         )
   }
macro Stephens_monoclinic
   (
      etac, etav,
      s400c,s400v,
      s040c,s040v,
      s004c,s004v,
      s220c,s220v,
      s202c,s202v,
      s022c,s022v,
      s121c,s121v,
      s103c,s103v,
      s301c,s301v
   )
   {
      #m_argu s400c If_Prm_Eqn_Rpt(s400c, s400v,)
      #m_argu s040c If_Prm_Eqn_Rpt(s040c, s040v,)
      #m_argu s004c If_Prm_Eqn_Rpt(s004c, s004v,)
      #m_argu s220c If_Prm_Eqn_Rpt(s220c, s220v,)
      #m_argu s202c If_Prm_Eqn_Rpt(s202c, s202v,)
      #m_argu s022c If_Prm_Eqn_Rpt(s022c, s022v,)
      #m_argu s121c If_Prm_Eqn_Rpt(s121c, s121v,)
      #m_argu s103c If_Prm_Eqn_Rpt(s103c, s103v,)
      #m_argu s301c If_Prm_Eqn_Rpt(s301c, s301v,)
      Stephens_lor_gauss(
         etac, etav,
         CeV(s400c,s400v) H^4 +
         CeV(s040c,s040v) K^4 +
         CeV(s004c,s004v) L^4 +
         CeV(s220c,s220v) H^2 K^2 +
         CeV(s202c,s202v) H^2 L^2 +
         CeV(s022c,s022v) K^2 L^2 +
         CeV(s121c,s121v) H K^2 L +
         CeV(s103c,s103v) H L^3 +
         CeV(s301c,s301v) H^3 L
         )
   }
macro Stephens_orthorhombic
   (
      etac, etav,
      s400c,s400v,
      s040c,s040v,
      s004c,s004v,
      s220c,s220v,
      s202c,s202v,
      s022c,s022v
   )
   {
      #m_argu s400c If_Prm_Eqn_Rpt(s400c, s400v,)
      #m_argu s040c If_Prm_Eqn_Rpt(s040c, s040v,)
      #m_argu s004c If_Prm_Eqn_Rpt(s004c, s004v,)
      #m_argu s220c If_Prm_Eqn_Rpt(s220c, s220v,)
      #m_argu s202c If_Prm_Eqn_Rpt(s202c, s202v,)
      #m_argu s022c If_Prm_Eqn_Rpt(s022c, s022v,)
      Stephens_lor_gauss(
         etac, etav,
         CeV(s400c,s400v) H^4 +
         CeV(s040c,s040v) K^4 +
         CeV(s004c,s004v) L^4 +
         CeV(s220c,s220v) H^2 K^2 +
         CeV(s202c,s202v) H^2 L^2 +
         CeV(s022c,s022v) K^2 L^2
         )
   }
macro Stephens_tetragonal_low
   (
      etac, etav,
      s400c,s400v,
      s004c,s004v,
      s220c,s220v,
      s202c,s202v,
      s310c,s310v
   )
   {
      #m_argu s400c If_Prm_Eqn_Rpt(s400c, s400v,)
      #m_argu s004c If_Prm_Eqn_Rpt(s004c, s004v,)
      #m_argu s220c If_Prm_Eqn_Rpt(s220c, s220v,)
      #m_argu s202c If_Prm_Eqn_Rpt(s202c, s202v,)
      #m_argu s310c If_Prm_Eqn_Rpt(s310c, s310v,)
      Stephens_lor_gauss(
         etac, etav,
         CeV(s400c,s400v) (H^4 + K^4) +
         CeV(s004c,s004v) L^4 +
         CeV(s220c,s220v) H^2 K^2 +
         CeV(s202c,s202v) (H^2 L^2 + K^2 L^2) +
         CeV(s310c,s310v) (H^3 K - H K^3)
         )
   }
macro Stephens_tetragonal_high
   (
      etac, etav,
      s400c,s400v,
      s004c,s004v,
      s220c,s220v,
      s202c,s202v
   )
   {
      #m_argu s400c If_Prm_Eqn_Rpt(s400c, s400v,)
      #m_argu s004c If_Prm_Eqn_Rpt(s004c, s004v,)
      #m_argu s220c If_Prm_Eqn_Rpt(s220c, s220v,)
      #m_argu s202c If_Prm_Eqn_Rpt(s202c, s202v,)
      Stephens_lor_gauss(
         etac, etav,
         CeV(s400c,s400v) (H^4 + K^4) +
         CeV(s004c,s004v) L^4 +
         CeV(s220c,s220v) H^2 K^2 +
         CeV(s202c,s202v) (H^2 L^2 + K^2 L^2)
         )
   }
macro Stephens_trigonal_low
   (
      etac, etav,
      s400c,s400v,
      s004c,s004v,
      s220c,s220v,
      s202c,s202v,
      s310c,s310v
   )
   {
      #m_argu s400c If_Prm_Eqn_Rpt(s400c, s400v,)
      #m_argu s004c If_Prm_Eqn_Rpt(s004c, s004v,)
      #m_argu s220c If_Prm_Eqn_Rpt(s220c, s220v,)
      #m_argu s202c If_Prm_Eqn_Rpt(s202c, s202v,)
      #m_argu s310c If_Prm_Eqn_Rpt(s310c, s310v,)
      Stephens_lor_gauss(
         etac, etav,
         CeV(s400c,s400v) (H^4 + 2 H^3 K + 3 H^2 K^2 + 2 H K^3 + K^4)+
         CeV(s004c,s004v) L^4 +
         CeV(s220c,s220v) H^2 K^2 +
         CeV(s202c,s202v) (H^2 L^2 + K^2 L^2) +
         CeV(s310c,s310v) (H^3 K - H K^3)
         )
   }
macro Stephens_trigonal_high
   (
      etac, etav,
      s400c,s400v,
      s004c,s004v,
      s202c,s202v,
      s211c,s211v
   )
   {
      #m_argu s400c If_Prm_Eqn_Rpt(s400c, s400v,)
      #m_argu s004c If_Prm_Eqn_Rpt(s004c, s004v,)
      #m_argu s202c If_Prm_Eqn_Rpt(s202c, s202v,)
      #m_argu s211c If_Prm_Eqn_Rpt(s211c, s211v,)
      Stephens_lor_gauss(
         etac, etav,
         CeV(s400c,s400v) (H^4 + 2 (H^3 K + H K^3) + 3 H^2 K^2 + K^4) +
         CeV(s004c,s004v) L^4 +
         CeV(s202c,s202v) (H^2 L^2 + H K L^2 + K^2 L^2) +
         CeV(s211c,s211v) (H^2 K L - H K^2 L + (H^3 L - K^3 L) 2 / 3)
         )
   }
macro Stephens_trigonal_high_2
   (
      etac, etav,
      s400c,s400v,
      s004c,s004v,
      s202c,s202v,
      s211c,s211v
   )
   {
      #m_argu s400c If_Prm_Eqn_Rpt(s400c, s400v,)
      #m_argu s004c If_Prm_Eqn_Rpt(s004c, s004v,)
      #m_argu s202c If_Prm_Eqn_Rpt(s202c, s202v,)
      #m_argu s211c If_Prm_Eqn_Rpt(s211c, s211v,)
      Stephens_lor_gauss(
         etac, etav,
         CeV(s400c,s400v) (H^4 + 2 H^3 K + 3 H^2 K^2 + 2 H K^3 + K^4) +
         CeV(s004c,s004v) L^4 +
         CeV(s202c,s202v) (H^2 L^2 + H K L^2 + K^2 L^2) +
         CeV(s211c,s211v) (H^2 K L - H K^2 L)
         )
   }
macro Stephens_hexagonal
   (
      etac, etav,
      s400c,s400v,
      s004c,s004v,
      s202c,s202v
   )
   {
      #m_argu s400c If_Prm_Eqn_Rpt(s400c, s400v,)
      #m_argu s004c If_Prm_Eqn_Rpt(s004c, s004v,)
      #m_argu s202c If_Prm_Eqn_Rpt(s202c, s202v,)
      Stephens_lor_gauss(
         etac, etav,
         CeV(s400c,s400v) (H^4+ 2 H^3 K + 3 H^2 K^2 + 2 H K^3 + K^4) +
         CeV(s004c,s004v) L^4 +
         CeV(s202c,s202v) (H^2 L^2 + H K L^2 + K^2 L^2)
         )
   }
macro Stephens_cubic
   (
      etac, etav,
      s400c,s400v,
      s220c,s220v
   )
   {
      #m_argu s400c If_Prm_Eqn_Rpt(s400c, s400v,)
      #m_argu s220c If_Prm_Eqn_Rpt(s220c, s220v,)
      Stephens_lor_gauss(
         etac, etav,
         CeV(s400c,s400v) (H^4 + K^4 + L^4) +
         CeV(s220c,s220v) (H^2 K^2 + H^2 L^2 + K^2 L^2)
         )
   }
' ---------------------------------------
macro Absorption_Edge_Correction_Eqn(& edge, & a_white, & b_white, & a_erf, & edge_extra)
   {
      modify_peak_eqn =
         (
            Get(current_peak) +
            a_white Exp(- b_white (Get(current_peak_x) - edge)^2) / Tan(Th)
            ' Not necessary
            ' If(Get(current_peak_x) < edge, 1/Get(current_peak_x)^3, 1 )
         )
         (edge_extra + 0.5 (1 + Erf_Approx( a_erf (Get(current_peak_x) - edge))) );
   }
macro Absorption_Edge_Correction(
   & max_lam,
   cedge,       vedge,
   ca_white,    va_white,
   cb_white,    vb_white,
   ca_erf,      va_erf,
   cedge_extra, vedge_extra)                               
   {
      current_peak_max_x = max_lam;
      #m_argu ca_erf
      #m_argu ca_white
      #m_argu cb_white
      #m_argu cedge
      #m_argu cedge_extra
      If_Prm_Eqn_Rpt(cedge, vedge, min 1e-6)
      If_Prm_Eqn_Rpt(ca_white, va_white, min 1e-6)
      If_Prm_Eqn_Rpt(cb_white, vb_white, min 1e-6)
      If_Prm_Eqn_Rpt(ca_erf, va_erf, min 60 max 3000)
      If_Prm_Eqn_Rpt(cedge_extra, vedge_extra, min 1e-6)
      Absorption_Edge_Correction_Eqn(
         CeV(cedge, vedge),
         CeV(ca_white, va_white),
         CeV(cb_white, vb_white),
         CeV(ca_erf, va_erf),
         CeV(cedge_extra, vedge_extra))
   }
macro Cu6_Ni_Edge
   (
      vymin_on_ymax,
      ccuKb, vcuKb,
      max_lam,
      cedge,       vedge,
      ca_white,    va_white,
      cb_white,    vb_white,
      ca_erf,      va_erf,
      cedge_extra, vedge_extra
   )
   {
      #m_argu ccuKb
      If_Prm_Eqn_Rpt(ccuKb, vcuKb,)
      lam
         ymin_on_ymax vymin_on_ymax
         la  0.0159 lo  1.534753 lh  3.6854
         la  0.5791 lo  1.540596 lh  0.437
            lo_ref
         la  0.0762 lo  1.541058 lh  0.6
         la  0.2417 lo  1.54441  lh  0.52
         la  0.0871 lo  1.544721 lh  0.62

         la  = CeV(ccuKb, vcuKb); lo  1.3922160 lh  0.5502872

         Absorption_Edge_Correction
         (
            max_lam,
            cedge,       vedge,
            ca_white,    va_white,
            cb_white,    vb_white,
            ca_erf,      va_erf,
            cedge_extra, vedge_extra
         )
   }
' ---------------------------------------
macro Pen_Wt
   {
      Write_Pen_Wt(
         Max(Abs(Get(pen_Aii)),1e-15),
         Max(Abs(Get(pen_Ai )),1e-15),
         Max(Abs(Get(pen_Pii)),1e-15),
         Max(Abs(Get(pen_Pi )),1e-15)
         )
   }
' ---------------------------------------
macro Deconvolution_Init(c)
   {
      process_times
      A0_matrix_is_constant
      penalties_weighting_K1 = c;
      save_best_chi2
      chi2_convergence_criteria 1e-5
      continue_after_convergence
      pen_weight 1
   }
macro Clamp_X(& x)
   {
      If(x < X1, X1, If(x > (X2 - 1e-6), X2, x))
   }
macro Deconvolution_Intensity_Penalty(i_name, fn_name)
   {
      fn fn_name (x, a0, a1) = (a0 - a1)^2 / Max(a0 + a1, 1e-6);
      default_I_attributes 1e-6 min 0 _v = Val Rand(0.99, 1.01); 
      create_pks_fn fn_name
      create_pks_name % i_name
   }
macro Deconvolution_Bkg_Penalty(& c, & w_min)
   {
      xdd_sum #m_unique pen = (Yobs - Get(bkg))^2 / Max(Get(bkg) Yobs, w_min^2);          
      penalty = pen c;   
   }
macro Deconvolution_Bkg_Penalty(& c)
   {
      Deconvolution_Bkg_Penalty(c, 1)
   }        
' ---------------------------------------
#include pdf.inc
macro Gr_bkg(ac, av, bc, bv, cc, cv)		 
	{
      #m_argu ac
      #m_argu bc
      #m_argu cc
      If_Prm_Eqn_Rpt(ac, av,)
      If_Prm_Eqn_Rpt(bc, bv,)
      If_Prm_Eqn_Rpt(cc, cv,)
		fit_obj = CeV(ac, av) Cos(CeV(bc, bv) X + CeV(cc, cv));
	}			 
macro Gr_hat(c, v, numhats)
	{
      #m_argu c
      If_Prm_Eqn_Rpt(c, v, min 1e-5 max = 2 Val + 0.1;)
      pdf_convolute = numhats; min_X = -CeV(c, v) 0.5; max_X = CeV(c, v) 0.5;
	}			 		 
' ---------------------------------------
#endif								  
                                  