TOPAS Macro Browser

load u11 u22 u33 u12 u13 u23

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));

(q / R + Exp_Rep(q, ro, b))

(q / R^n + Exp_Rep_n(q, ro, b, n))

box_interaction s1 s2 c = Gauss_Rep_n(ro, y, n);

box_interaction s1 s2 c = q /R + Rn_Rep(Sqrt((q) (q)), ro, n);

(Exp(- b R) Exp(b ro) Sqrt((q)^2) / (b ro^2))

(Exp(b (ro - R)) n Sqrt((q)^2) / (b ro^(n+1)))

(Exp(R^n Ln(y) / ro^n))

grs_interaction qi wqi qj wqj s1 s2 c =
         If(R < ro,
            a (R-ro)^2,
            0
         );

Grs_BornMayer(s1, s2, wqi, wqj, c, ro, b, .2)

grs_interaction qi wqi qj wqj s1 s2 c =
            Exp_Rep(Sqrt( ((wqi) (wqj))^2 ), ro, b);

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)
         );

Grs_Interaction(s1, s2, wqi, wqj, c, ro, n, 0.2)

grs_interaction qi = wqi; qj = wqj; s1 s2 c = 0;

(q / n) (ro^(n-1)) / R^n

#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

AI_Anti_Bump(s1, s2, ro,, wby, num_cycle_iters)

AI_Anti_Bump(s1, s2, ro,, wby,)

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
         ;

#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)

Angle_Restrain(sites,, 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)

Distance_Restrain(sites,, t, t_calc, tol, wscale)

#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
            ;

Distance_Restrain_Keep_Out(sites, r,, 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
            ;

Distance_Restrain_Keep_Within(sites, r,, wby, num_cycle_iters)

#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

Flatten(sites,,t_calc, tol, wscale)

#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;

Parabola_N(n1, n2, s1, s2,, ro, wby)

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)

Set_Length(s0, s1, r, xc, yc, zc, cva1, cvb1)
      Set_Length(s0, s2, r, xc, yc, zc, cva2, cvb2)

Set_Lengths(s0, s1, s2, r, xc, yc, zc,
         cva1, cvb1,
         cva2, cvb2)
      Set_Length(s0, s3, r, xc, yc, zc, cva3, cvb3)

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)

0.5 (x1 + x2) + 0.5 (x2 - x1) Cos(2 Pi Cycle_Iter / n)

x1 + (x2-x1) Mod(Cycle, n) / (n-1)

Min(x2, Max(x, x1))

If(Abs(x) < lim, If(x < 0, -lim, lim), x)

I_parameter_names_have_hkl int
      out_A_matrix file
         A_matrix_prm_filter int*

pick_atoms_when = Mod(Cycle_Iter, n) == (n-1);

x1 + (x2-x1) Mod(Cycle_Iter, n) / (n-1)

x1 + (x2-x1) Mod(i, n2) / (n1-1)

If(Mod(Cycle_Iter, n2) >= n1,
         x2,
         Ramp_0(Cycle_Iter, x1, x2, n1, n2)
      )

If(Mod(Cycle_Iter, n) >= n,
         x2,
         Ramp_0(Cycle_Iter, x1, x2, n, n)
      )

x1 + (x2-x1) Mod(Run_Number, n) / (n-1)

0.5 (x1 + x2) + 0.5 (x2 - x1) Sin(2 Pi Cycle_Iter / n - Pi / 2)

0.5 (x1 + x2) + 0.5 (x2 - x1) Sin(2 Pi (Cycle_Iter - x0) / n - Pi / 2)

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;

Tangent(fraction_n, max_triplets, 35000)

If(x < X1, X1, If(x > (X2 - 1e-6), X2, x))

xdd_sum #m_unique pen = (Yobs - Get(bkg))^2 / Max(Get(bkg) Yobs, w_min^2);          
      penalty = pen c;

Deconvolution_Bkg_Penalty(c, 1)

process_times
      A0_matrix_is_constant
      penalties_weighting_K1 = c;
      save_best_chi2
      chi2_convergence_criteria 1e-5
      continue_after_convergence
      pen_weight 1

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

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))

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))) );

Load_Lam(coka3, yminymax)

Load_Lam(coka7_holzer, yminymax)

Load_Lam(cokb6_holzer, yminymax)

Load_Lam(crka7_holzer, yminymax)

Load_Lam(crkb5_holzer, yminymax)

#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
         )

Load_Lam(cuk1sharp, yminymax)

Load_Lam(cuka1, yminymax)

Load_Lam(cuka2, yminymax)

Load_Lam(cuka2_analyt, yminymax)

Load_Lam(cuka4_holzer, yminymax)

Load_Lam(cuka5, yminymax)

Load_Lam(cuka5_berger, yminymax)

Load_Lam(cukb4_holzer, yminymax)

Load_Lam(feka7_holzer, yminymax)

Load_Lam(fekb4_holzer, yminymax)

#include root##lam\##lam_file##.lam ymin_on_ymax yminymax

Load_Lam(ROOT, lam_file, yminymax)

Load_Lam(mnka7_holzer, yminymax)

Load_Lam(mnkb5_holzer, yminymax)

Load_Lam(moka2, yminymax)

Load_Lam(nika5_holzer, yminymax)

Load_Lam(nikb4_holzer, yminymax)

lam no_th_dependence   la 1 lo 0 lh 1

Load_Lam(nothdependence, yminymax)

#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);

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)

Bkg_Straight_Line_xy(n, x, y)

prm !bkg_x_##n = x;
      prm !bkg_y_##n = y;

(2 Sqrt(Ln(2) / Pi) / fwhm) Exp(-4 Ln(2)((X - xo) / fwhm)^2)

(2 / (Pi fwhm)) / (1 + 4 ((X-xo) / fwhm)^2)

One_on_X(,v)

#m_argu c
      If_Prm_Eqn_Rpt(c, v, min 0.001 )
      fit_obj = CeV(c,v) / X;

dummy_str
			space_group 1
         phase_name name
         scale = fit_obj_scale;

dummy_str
			space_group 1
         phase_name name

fit_obj = a (1-g) Lorentzian(xo, fwhm) + a g Gauss(xo, fwhm);

#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)

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;

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)

#include Double_Quote##ROOT##Charge_Flipping.inc##Double_Quote

#include Double_Quote##ROOT##GUI_PDF.inc##Double_Quote

#include Double_Quote##ROOT##pdf-generate.inc##Double_Quote

#include Double_Quote##ROOT##pdf-generate-common.inc##Double_Quote

try_space_groups "196 197 195"

try_space_groups "5 3"

try_space_groups "22 23 21 16"

try_space_groups "79 75"

try_space_groups "2"

try_space_groups "146 143"

try_space_groups "196"

try_space_groups "197"

try_space_groups "195"

index_x0 = 1 / ( (worder) d)^2;

Index_x0_from_d(d, 1)

index_x0 = ( 2 Sin(th2 Pi/360) / (Get(index_lam) (worder))  )^2;

Index_x0_from_th2(th2, 1)

load dummy ndx_sg ndx_uni ndx_vol ndx_gof ndx_alp ndx_blp ndx_clp ndx_allp ndx_belp ndx_galp

load dummy ndx_sg ndx_status ndx_uni ndx_vol ndx_gof ndx_alp ndx_blp ndx_clp ndx_allp ndx_belp ndx_galp

load dummy ndx_sg ndx_uni ndx_vol ndx_gof ndx_ze ndx_alp ndx_blp ndx_clp ndx_allp ndx_belp ndx_galp

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

try_space_groups "5"

try_space_groups "3"

try_space_groups "21"

try_space_groups "22"

try_space_groups "23"

try_space_groups "16"

try_space_groups "79"

try_space_groups "75"

try_space_groups "2"

try_space_groups "143"

try_space_groups "146"

try_space_groups "228 219 203 210 196 230 220 206 214 197 222 218 201 205 212 198 195"

try_space_groups "9 5 14 7 4 3"

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"

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"

try_space_groups "2"

try_space_groups "161 146 184 159 158 169 144 173 143"

AB(,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

AB

#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)));

#m_argu u
      #m_argu d
      #m_ifarg u ""
         #m_unique_not_refine u
      #m_endif
      AB(u, uv)
      Sample_Thickness(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)));

Capillary_Offset_Cos_2Th_mm(, v)

#m_argu c
      If_Prm_Eqn_Rpt(c, v, min -10 max 10)
      th2_offset = CeV(c, v) Rad Cos(2 Th) / Rs;

Capillary_Offset_Sin_2Th_mm(, v)

#m_argu c
      If_Prm_Eqn_Rpt(c, v, min -10 max 10)
      th2_offset = CeV(c, v) Rad Sin(2 Th) / Rs;

Cylindrical_2Th_Correction(,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);

Cylindrical_I_Correction(,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;

Divergence(,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);

#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);

#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);

axial_conv
         filament_length            s2
         sample_length              = Get(source_len);
         receiving_slit_length      h2
         primary_soller_angle   0

axial_conv
         filament_length           filament_cv
         sample_length             sample_cv
         receiving_slit_length     detector_cv

Full_Axial_Model(filament_cv, sample_cv, detector_cv)
      Sollers(psol_cv, ssol_cv)

Radius(r,r)

Rp rp Rs rs

Sample_Thickness(,v)

#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;

Simple_Axial_Model(,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);

#m_argu k
      If_Prm_Eqn_Rpt(k, kv, min 0.001 max 100)
      Sinc_Convolution(CV(k, kv), x_range)

user_defined_convolution = If(Abs(X) < 1e-5, 1, (Sin(k X) /(k X))^2);
         min = -x_range; max = x_range;

SW

primary_soller_angle    cv_p
      secondary_soller_angle  cv_s

Source_Width(,v)

#m_argu c
      If_Prm_Eqn_Rpt(c, v, min 0.00001 max 2)
      hat = Rad CeV(c,v) / Rs;

SW(,v)

#m_argu c
      If_Prm_Eqn_Rpt(c, v, min 0.00001 max 2)
      hat = Rad CeV(c,v) / Rs;

#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);

Variable_Divergence(,v)

Variable_Divergence_Shape(c, v)
      Variable_Divergence_Intensity

scale_pks = Sin(Th);

Variable_Divergence_Shape(,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);

scale_pks = 1 / (Sin(Th)^2 Cos(Th));

Lorentz_Factor

LP_Factor_X(,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));

#m_argu c
      If_Prm_Eqn_Rpt(c, v, min 0.0001 max = 2 Val + .5;)
      PO(CeV(c, v), ang, hkl)

#m_ifarg md ""
         #m_unique_not_refine md
      #m_endif
      march_dollase md hkl = c;
      scale_pks = md;

(c^2 Cos(ang)^2 + Sin(ang)^2 / c) ^ (-1.5)

#m_ifarg sh ""
         #m_unique_not_refine sh
      #m_endif
      spherical_harmonics_hkl sh
         sh_order order
      scale_pks = sh;

#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;

PO

#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));

#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)));

a cv b = Get(a); c = Get(a);

a a_cv b = Get(a); c c_cv ga 120

a a_cv b = Get(a); c = Get(a);
      al al_cv be = Get(al); ga = Get(al);

a a_cv b = Get(a); c c_cv

Hexagonal

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

#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

local = Get(mlx) Get(a); : mxc
     local = Get(mly) Get(b); : myc
     local = Get(mlz) Get(c); : mzc

#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

load exclude

#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));

#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;

' http://topas.dur.ac.uk/unb/forum.php?req=thread&id=32
      #m_argu c
      c v
      = c^2;

#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

load la lo lh

Sum( sum_eqn, Mi = 0, Mi < M, Mi = Mi+1) / M

#prm a = n; #out a

remove_phase = And(Get(weight_percent) < wt, Error(Get(weight_percent)) > fr Get(weight_percent));

#if (n > 0) 
         x Repeat(x, Numeric(n-1)) 
      #endif

load site x y z occ beq

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

#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

#m_argu c
      If_Prm_Eqn_Rpt(c, v, min 3 max 10000)
      gauss_fwhm = Constant(t1 .1) D_spacing^2 / CeV(c, v);

#m_argu c
      If_Prm_Eqn_Rpt(c, v, min 3 max 10000)
      lor_fwhm = Constant(t1 .1) D_spacing^2 / CeV(c, v);

#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);

xdd path gsas_format
         neutron_data
         x_calculation_step calc_step
         weighting = If(SigmaYobs < 1, 1, 1/SigmaYobs^2);

lam
         ymin_on_ymax w
         la 1 lo 0 lh 1

#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;

#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);

xdd path xye_format
         x_calculation_step calc_step
         neutron_data
         weighting = If(SigmaYobs < 1, 1, 1/SigmaYobs^2);

out_record
         out_eqn = eqn; 
         out_fmt fmt

out_record
         out_eqn = eqn; 
         out_fmt "%.9g"

Out(eqn, fmt)
         out_fmt_err fmt_err

out file append
      Out(eqn, fmt)

out file append
      Out(eqn, fmt)
         out_fmt_err fmt_err

out_record
         out_fmt fmt

xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f  " = X;
          " %11.6f\n" = Yobs-Ycalc;
      }

xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f  " = X;
          " %11.6f\n" = y;
      }

xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f  " = X;
			 " %11.6f\n" = Ycalc;
      }

xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f  " = X;
          " %11.6f\n" = Yobs;
      }

xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f  " = X;
          " %11.6f  " = Yobs;
          " %11.6f\n" = Ycalc;
      }

Out_XDD_Start_Step(file)
      xdd_out file append out_record out_fmt " %.9g\n" out_eqn

out file         
         Out(X1) Out_String(" ") 
         Out(Yobs_dx_at(X1)) 
         Out_String("\n")

xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f  " = X;
          " %11.6f\n" = y;
      }

xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f\n" = Ycalc;
      }

xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f  " = Yobs;
          " %11.6f  " = Ycalc;
          " %11.6f\n" = Yobs-Ycalc;
      }

xdd_out file load out_record out_fmt out_eqn
      {
          " %11.6f  " = X;
          " %11.6f\n  " = Max(Rand_Normal(Ycalc, Sqrt(Ycalc)), 0);
      }

phase_out file load out_record out_fmt out_eqn
      {
         " %11.6f"   = 2 Rad Th;
         " %11.6f\n" = wI;
      }

Create_2Th_Ip_file(file, I_no_scale_pks)

Create_2Th_Ip_file(file, I_after_scale_pks)

phase_out file load out_record out_fmt out_eqn
      {
         " %11.6f"   = Lam / (2 Sin(Th));
         " %11.6f\n" = wI;
      }

Create_d_Ip_file(file, I_no_scale_pks)

Create_d_Ip_file(file, I_after_scale_pks)

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;
      }

Create_hklm_d_Th2_Ip_file(file, I_no_scale_pks)

Create_hklm_d_Th2_Ip_file(file, I_after_scale_pks)

#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);

#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);

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;

#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))

#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);

brindley_spherical_r_cm = R PD;

move_to xdd
      local !k_factor_mac_local_ mac
      local !k_factor_k_local_ k
      local !k_factor_sum_wps_  100 : amor

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_;
      }

scale = (w/(100-w)) Get(sum_smvs_minus_this) / (Get(cell_mass) Get(cell_volume));

local = Get(mixture_MAC) Get(mixture_density_g_on_cm3); : mlac

cell_mass       m_v
      cell_volume     v_v
      weight_percent  w_v

1.6605402 Get(cell_mass) / Get(cell_volume)

local = Phase_Density_Eqn_g_on_cm3; : pd

local = Phase_LAC_Eqn_1_on_cm; : u

(Get(phase_MAC) Phase_Density_Eqn_g_on_cm3)

point_for_site new_point
         ux = Point(prev_point, rx);
         uy = Point(prev_point, ry);
         uz = Point(prev_point, rz);

Duplicate_Point(new_point, prev_point)
      rotate angle_cv qz 1 operate_on_points new_point

/*
         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;

/*
         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;

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;

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)

Find_Child(Get(point_for_sites), a)

Point_Current_Value(Get(Point(a), xyz))

point_for_site site ux x uy y uz z

Point_rx_to_ua(Point(site_name))

Point_ry_to_ub(Point(site_name))

Point_rz_to_uc(Point(site_name))

rotate cva qa 1
      rotate cvb qb 1

Rotate_about_axies(cva, cvb)
      rotate cvc qc 1

Rotate_about_these_points(cv, a, b, operate_on_points Double_Quote !##a !##b * Double_Quote)

Rotate_about_these_points(cv, a, b, operate_on_points pts)

Translate_point_amount(a, -) ops
      rotate cv Rotation_vector_from_points(a, b) ops
      Translate_point_amount(a) ops

qx = Point(b, rx) - Point(a, rx);
      qy = Point(b, ry) - Point(a, ry);
      qz = Point(b, rz) - Point(a, rz);

#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

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)

translate ta acv tb bcv tc ccv

Translate(acv, bcv, ccv)
         operate_on_points ops

Translate_point_amount(a, )

translate
         tx = pm Point(a, rx);
         ty = pm Point(a, ry);
         tz = pm Point(a, rz);

Translate(xc 0, yc 0, zc 0)
         start_values_from_site s0

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;

point_for_site s0
        Triangle(s1, s2, s3, r)

rigid
         Triangle(s0, s1, s2, s3, r)
         Rotate_about_axies(cva, cvb, cvc)
         Translate_with_site_start_values(s0, xc, yc, zc)

CS_L

CS_G(,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));

CS_L(,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));

#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

#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

Strain_L

Strain_L

local = If(Abs(fwhm) < Abs(s),
         -9999,
         .1 Rad Lam k / (Cos(Get(xo) Deg_on_2) Sqrt(fwhm^2-s^2))
         ); : csize

CS_G

Specimen_Tilt(,v)

#m_argu c
      If_Prm_Eqn_Rpt(c, v, min 0.0001 max 2)
      hat = Cos(Th) Rad CeV(c,v) / Rs;

#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)
         )

#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)
         )

#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);

#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
         )

#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
         )

#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)
         )

#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)
         )

#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
         )

#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)
         )

#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)
         )

#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)
         )

Strain(,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)

Strain_G(,v)

#m_argu c
      If_Prm_Eqn_Rpt(c, v, min .0001 max 5)
      gauss_fwhm = CeV(c,v) Tan(Th);

Strain_L(,v)

#m_argu c
      If_Prm_Eqn_Rpt(c, v, min .0001 max 5)
      lor_fwhm = CeV(c,v) Tan(Th);

#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

#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;

SF_smooth(,1,s)

transition t
         use_layer l
         height = h;

transition l
         use_layer l
         height = h;

(A01^2 + B01^2 + A11^2 + B11^2)

(F_Real_negative^2 + F_Imaginary_negative^2)

(F_Real_positive^2 + F_Imaginary_positive^2)

(A11-B01)

(A11+B01)

(A01+B11)

(A01-B11)

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;
      }

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;
      }

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");
             }

out file append
         Out(Get(cif_bonds_angles), "%s")

Out_CIF_STR(file,)

' 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");
               }

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;
      }

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;
            }

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;
      }

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;
      }

Simulated_Annealing_1(0.001)
		quick_refine 0.01
      randomize_on_errors
      TE_1(t, 1 1 1, 1)

x_calculation_step = diff_toll;
      yobs_to_xo_posn_yobs = Get(x_calculation_step);

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);

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)
         )

seed
      verbose -1
      continue_after_convergence
      chi2_convergence_criteria = (v);
      process_times

weighting = 1 / Sqrt(Max(Yobs, 1));

temperature = (t);
         move_to_the_next_temperature_regardless_of_the_change_in_rwp
         on_best_rewind
      load temperature { ts }

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

load temperature

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

xdd w##.brml

move_to xdds
      move_to xdd_recs
      new(xdd, n)

xdd path##.dat

xdd path##.raw

xdd path##.raw range range_num

xdd_scr path##.scr

xdd_scr path##.hkl

xdd path##.xdd

str space_group sg

str space_group sg phase_name pn

xdd path##.xdd

xdd path##.xy
      x_calculation_step calc_step

xdd path_ext xye_format

SD(, 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;

SD

ZE(,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);

ZE

#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);

#if (Not_Blank(capillary_scan))
			if (Prm_There(p0)) {
				user_y cap capillary_scan
					#if (capillary_rebin)			  
						user_y_rebin_with_dx_of = capillary_rebin;
					#endif
				fit_obj = (p0 + p1 x0) cap;
				existing_prm bkg_tot_ += (p0 + p1 x0) cap;
			}
		#endif

process_times
		A0_matrix_is_constant
		penalties_weighting_K1 1
		save_best_chi2
		chi2_convergence_criteria 1e-5
		continue_after_convergence
		pen_weight 1

prm !da = Ramp_Clamp(100, 1, 0, niters);
		prm #m_unique gap = yobs - bkg_0_;

		xdd_sum #m_unique bkg_1_ =
			If(yobs < bkg_0_, 1e4, 1)
			If(Or(yobs < bkg_0_, gap < (da Sqrt(yobs))),
				gap,
				da Sqrt(yobs)
			)^2 / Max(yobs, 1);
		penalty = bkg_1_ damp;

prm #m_unique gap = yobs - bkg_0_;
		xdd_sum #m_unique bkg_1_ =
			If(yobs < bkg_0_, 1e4, 1)
			gap^2 / Max(yobs, 1);
		penalty = bkg_1_ damp;

(4 Pi / Lam) Sin(th2 Pi / 360)

(4 Pi / Lam) Sin(th2 0.5)

Th_Rad_from_Q(q) 360 / Pi

2 Th_Rad_from_Q(q)

Th_Rad_from_Q(q) 180 / Pi

If(q Lam > (4 Pi), Throw_Error("q Lam > (4 Pi)"), 1)  ArcSin(q Lam / (4 Pi))

out file
      #if (Run_Number > 0)
         append
      #endif

system_before_save_OUT  {
         copy INP_File##.inp INP_File##.backup
      }

CV(c, v)

#m_ifarg v ""
         0
      #m_else
         CeV(c, v)
      #m_endif

CS_L

#m_ifarg c ""
			#m_first_word v v
      #m_else
         c
      #m_endif

0.01745329251994

0.00872664625997

NoThDependence(0.0001)
      peak_type pv
         pv_lor 0
         pv_fwhm = 0.01 Peak_Calculation_Step;

"

,

)

Mod(n, 2) == 0

#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)

#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)

#m_ifarg c #m_eqn
         = c; : v
      #m_else
         c v
      #m_endif

If(

#m_ifarg c ""
         #m_ifarg v ""
         #m_else
            w
         #m_endif
      #m_else
           w
      #m_endif

#m_ifarg c ""
      #m_else
         w
      #m_endif

#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

If_Prm_Eqn_Rpt(c, v, dets,)

LP_Factor(,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));

= c; c v min 0

#m_ifarg s "" 
         0   
      #m_else
         1
      #m_endif

Mod(n, 2)

3.14159265358979

PVII_h(,v)

#m_argu c
      If_Prm_Eqn_Rpt(c, v, min .0001 max 20)
      h1 = CeV(c, v) .5;
      h2 = Get(h1);

PVII_m(,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);

57.2957795130823

use_Fc
      randomize_initial_phases_by 0

/*
      ' 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

#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"), "");

peak_type spvii

,

If(Prm_There(s##_##c), s##_##c^2, 0)

#m_argu c1
		#m_argu c2
		#m_argu c3
		#m_argu c4
		If_Prm_Eqn_Rpt(c1, v1, min 1e-6 max 10 val_on_continue = Rand(.01, 2); )
		If_Prm_Eqn_Rpt(c2, v2, min 1e-6 max 10 val_on_continue = Rand(.01, 2); )
		If_Prm_Eqn_Rpt(c3, v3, min 1e-6 max 10 val_on_continue = Rand(.01, 2); )
		If_Prm_Eqn_Rpt(c4, v4, min 1e-6 max 500 val_on_continue = Rand(.01, 2); )
		beq = CeV(c1, v1) Sqrt(Abs(1 - CeV(c2, v2) /  X - CeV(c3, v3) /  X^2 + 1e-5 CeV(c4, v4)  X^2));

#m_argu v1c
	#m_argu v2c
	#m_argu v3c
	#m_argu v4c
	If_Prm_Eqn_Rpt(v1c, v1v, min 0.000001 max 10, del = 0.01 Val;)
	If_Prm_Eqn_Rpt(v2c, v2v, min 0.000001 max 10, del = 0.01 Val;)
	If_Prm_Eqn_Rpt(v3c, v3v, min 0.000001 max 10, del = 0.01 Val;)
	If_Prm_Eqn_Rpt(v4c, v4v, min 0.000001 max 10, del = 0.01 Val;)
	beq = CeV(v1c, v1v) + CeV(v2c, v2v) X + CeV(v3c, v3v) / X + CeV(v4c, v4v) / X^2;

#m_argu uiso
	#m_argu rcut
	#m_argu sratio
	#m_argu delta1
	#m_argu delta2
	#m_argu qbroad
	If_Prm_Eqn_Rpt(uiso, uisov,  min 0.000001 max 1, del = 0.01 Val;)
	If_Prm_Eqn_Rpt(rcut,  rcutv,   min 0.000001 max 10, del = 0.01 Val;)
	If_Prm_Eqn_Rpt(sratio,sratiov, min 0.000001 max 10, del = 0.01 Val;)
	If_Prm_Eqn_Rpt(delta1,delta1v, min 0.000001 max 1, del = 0.01 Val;)
	If_Prm_Eqn_Rpt(delta2,delta2v, min 0.000001 max 1, del = 0.01 Val;)
	If_Prm_Eqn_Rpt(qbroad,qbroadv, min 0.000001 max 1, del = 0.01 Val;)
	beq = If(X < CeV(rcut, rcutv), 0.5 CeV(sratio, sratiov), 1) CeV(uiso, uisov) 8 Pi^2 (1 - Min(((CeV(delta1, delta1v)/X) + (CeV(delta2, delta2v)/X^2) - (CeV(qbroad, qbroadv)^2 X^2)),1));

#m_argu beqc
	#m_argu r
	#m_argu r2
	If_Prm_Eqn_Rpt(beqc, beqv, min 0.001 max 10, del = 0.01 Val;)
	If_Prm_Eqn_Rpt(r, rv, min 0.001 max 10, del = 0.01 Val;)
	If_Prm_Eqn_Rpt(r2, r2v, min 0.00001 max 10, del = 0.01 Val;)
	beq = CeV(beqc, beqv) + CeV(r, rv) X + CeV(r2, r2v) X^2;

#m_argu beqlo
	#m_argu beqhi
	#m_argu rcut
	If_Prm_Eqn_Rpt(beqlo, beqlov, min 0.001 max 10, del = 0.05 Val;)
	If_Prm_Eqn_Rpt(beqhi, beqhiv, min 0.001 max 10, del = 0.05 Val;)
	beq = IF X < rcut THEN
				CeV(beqlo, beqlov)
			ELSE
				CeV(beqhi, beqhiv)
			ENDIF;

#m_argu rcut
	#m_argu beqcut
	#m_argu rlo
	#m_argu beqlo
	#m_argu beqhi
	#m_argu radius
	If_Prm_Eqn_Rpt(rcut, rcutv, min 0.1 max 10, del = 0.05 Val;)
	If_Prm_Eqn_Rpt(beqcut, beqcutv, min 0.001 max 10, del = 0.05 Val;)
	If_Prm_Eqn_Rpt(rlo, rlov, min 0.1 max 100, del = 0.05 Val;)
	If_Prm_Eqn_Rpt(beqlo, beqlov, min 0.001 max 10, del = 0.05 Val;)
	If_Prm_Eqn_Rpt(beqhi, beqhiv, min 0.001 max 100, del = 0.05 Val;)
	If_Prm_Eqn_Rpt(radius, radiusv, min 1 max 1000, del = 0.05 Val;)
	beq = IF X < CeV(rcut, rcutv) THEN
				CeV(beqcut, beqcutv)
			ELSE
				IF X < CeV(rlo, rlo) THEN
					CeV(beqlo, beqlov)
				ELSE
					IF X < (2 CeV(radius, radiusv)) + CeV(rlo, rlov) THEN
						CeV(beqlo, beqlov) + (CeV(beqhi, beqhiv)-CeV(beqlo, beqlov)) (1-((Pi (X+CeV(rlo, rlov))^2 ((0.25 ((X+CeV(rlo, rlov))/(CeV(radius, radiusv)+CeV(rlo, rlov)))^3)-(3 (X+CeV(rlo, rlov))/(CeV(radius, radiusv)+CeV(rlo, rlov)))+4)) / (4 Pi (X+CeV(rlo, rlov))^2)))
					ELSE
						CeV(beqhi, beqhiv)
					ENDIF
				ENDIF
			ENDIF;

beq_rcut_rlo_spherical(rcut,rcutv,beqcut,beqcutv,!rlo,0,beqlo,beqlov,beqhi,beqhiv,radius,radiusv)

beq_rcut_rlo_spherical(!rcut,0,!beqcut,0,rlo,rlov,beqlo,beqlov,beqhi,beqhiv,radius,radiusv)

beq_rcut_rlo_spherical(!rcut,0,!beqcut,0,!rlo,0,beqlo,beqlov,beqhi,beqhiv,radius,radiusv)

#m_argu c
      If_Prm_Eqn_Rpt(c, v, min 1e-6 max 5 val_on_continue = Val Rand(.1, 2); )
      beq = CeV(c, v) Erf_Approx(erf_a X);         'error function approximation

#m_argu c
      If_Prm_Eqn_Rpt(c, v, min 1e-6 max 5 val_on_continue = Val Rand(.1, 2); )
      beq = CeV(c, v) Sqrt(Abs(1 - delt1 /(1e-10+X) - delt2 /(1e-10+X)^2 + Qb^2 X^2));

#m_argu c
		If_Prm_Eqn_Rpt(c, v, min 1e-6 max 5 val_on_continue = Val Rand(.1, 2); )
		beq = CeV(c, v) Sqrt(1- Min((delt1/X + delt2/X^2 - Qb^2 X^2),1));	'8 Pi^2

'Convolution to account for a Q dependent broadening term in the S(Q) data
   #m_argu alpha
   If_Prm_Eqn_Rpt(alpha, alphav, min 0.00001 max 0.1, del 0.001)
	local #m_unique fwhm = Max(CeV(alpha,alphav) Xo, 0.01);
   pdf_convolute = Gauss(0, fwhm);
   	min_X = Min( Ln(0.001) / fwhm,0) ; 
   	max_X = Max(-Ln(0.001) / fwhm,0) ;
   	convolute_X_recal = If(Xo,1,1) 1;

'Convolution to account for the use of a Lorch function prior to Fourier transform of the S(Q)
   #m_argu Qmax
   If_Prm_Eqn_Rpt(Qmax, Qmaxv, min 5 max 60, del 0.1)
   local #m_unique fwhm = (Pi / CeV(Qmax,Qmaxv))/Ln(2);
   pdf_convolute = Gauss(0, fwhm);
   	min_X = Ln(0.001) / fwhm;
   	max_X =-Ln(0.001) / fwhm;
   	convolute_X_recal = 1;

'Convolution to include Sinc function due to termeination of S(Q) at finite Qmax during Fourier transform
	'Use in combination with start_X 1 to avoid very low-r data 
	'Ref: J. S. Chung and M. F. Thorpe, Phys. Rev. B, 55(3), 1997, 1545
	#m_argu Qmax
	If_Prm_Eqn_Rpt(Qmax, Qmaxv, min 1 max 100, del 0.01)
	local #m_unique conv_max = (5 CeV(Qmax,Qmaxv) - Mod(5 CeV(Qmax,Qmaxv),1))/5 2 Pi / CeV(Qmax,Qmaxv);
	pdf_convolute = If(Abs(X) > Yobs_dx_at(Xo),(Sin(CeV(Qmax,Qmaxv) X)/(X)),CeV(Qmax,Qmaxv));
		min_X = Min(-Xo,-conv_max) ;
		max_X = Max( Xo, conv_max) ;

'Covolution to account for the use of the Soper-Lorch function prior to Fourier transform of the S(Q)
	'd_zero is the "Width of broadening in r space", as used in GudrunN
	'Ref: A. K. Soper and E. R. Barney, J. Appl. Cryst. (2012). 45, 1314�1317
   #m_argu d_zero
   #m_argu d_one
   If_Prm_Eqn_Rpt(d_zero, d_zerov, min 0.001 max 1.0, del 0.01)
   If_Prm_Eqn_Rpt(d_one, d_onev, min 0 max 0.1, del 0.01)
   local #m_unique fwhm = 2.3548 (CeV(d_zero,d_zerov)/Ln(2) Ln(2 Cosh( (CeV(d_one,d_onev)/CeV(d_zero,d_zerov)) Ln(2) Xo)));
   pdf_convolute = Gauss(0, fwhm);
   	min_X = Ln(0.001) / fwhm;
   	max_X =-Ln(0.001) / fwhm;
   	convolute_X_recal = 1;

convolute_SoperLorch(d_zero,d_zerov,!d_one,0)

'Models dampening of the PDF intensity as a function of r due to a constant width reciprocal space peak shape
	'Ref: Acta Cryst. (1992). A48, 336-346
	#m_argu dQ 'Instrumental FWHM of S(Q) data
	If_Prm_Eqn_Rpt(dQ, dQv, min 0.001 max 0.2, del 0.0001)
	scale_phase_X = Exp(-0.5  X^2 (CeV(dQ,dQv)/2.35482)^2);

'A modification of dQ_dampening to correct for a Lorentzian component to the reciprocal space peak shape
	#m_argu dQ 'Instrumental FWHM of S(Q) data
	#m_argu lor 'Lorentzian contribute to peak shape (between 0 and 1)
	If_Prm_Eqn_Rpt(dQ, dQv, min 0.001 max 0.2, del 0.0001)
	If_Prm_Eqn_Rpt(lor, lorv, min 0 max 1)
	scale_phase_X = (1-CeV(lor,lorv)) Exp(-0.5  X^2 (CeV(dQ,dQv)/2.35482)^2) + CeV(lor,lorv) Exp(-0.5 CeV(dQ,dQv) X);

2.3548 ((beq1 + beq2) / (8 Pi^2))^0.5

'Rescaling to G'(r) according to Keen, J. Appl. Cryst. (2001). 34, 172-177
	'num_density is the number density in atoms per cubic Angstrom for the phase (within str) or sample (within xdd) 
	#m_argu num_density
	If_Prm_Eqn_Rpt(num_density, num_densityv, min 0 max 1, del 0.001)
	scale_phase_X = If(X>0.01,1/(4 Pi CeV(num_density,num_densityv) X),0);
   fit_obj = If(X>0.01,1,0);

'Rescaling to G(r) according to Keen, J. Appl. Cryst. (2001). 34, 172-177
	'num_density is the number density in atoms per cubic Angstrom for the phase (within str) or sample (within xdd)
	'weighted_av is the weighted average structure factor, <F>^2
	#m_argu num_density
	#m_argu weighted_av
	If_Prm_Eqn_Rpt(num_density, num_densityv, min 0 max 1, del 0.001)
	If_Prm_Eqn_Rpt(weighted_av, weighted_avv, min 0 max 2, del 0.001)
	scale_phase_X = If(X>0.01,CeV(weighted_av, weighted_avv) 1/(4 Pi CeV(num_density,num_densityv) X),0);
   fit_obj = If(X>0.01,0,-CeV(weighted_av, weighted_avv));

spherical_damping(r,rv)

#m_argu r
	If_Prm_Eqn_Rpt(r, rv, min 1 max 1000, del = 0.01 Val;)
	scale_phase_X = 
	IF X > 2 CeV(r,rv) THEN
		1
	ELSE
		If(X>0.01,1-((Pi X^2 ((0.25 (X/CeV(r,rv))^3)-(3 X/CeV(r,rv))+4)) / (4 Pi X^2 1)),0)
	ENDIF;

xdd_out file load out_record out_fmt out_eqn
   {
      " %11.6f  " = X; 
   	" %11.6f  " = Yobs;
      " %11.6f  " = Ycalc;
      " %11.6f\n" = Yobs-Ycalc;
   }

#m_argu CN
   #m_argu xo
   #m_argu beq
   If_Prm_Eqn_Rpt(CN, CNv,  min 0.00001)
   If_Prm_Eqn_Rpt(xo, xov, )
   If_Prm_Eqn_Rpt(beq, beqv,  min 0.00001 del 0.001)   
   fit_obj = (CeV(CN,CNv) / CeV(xo,xov)) Gauss(CeV(xo,xov), 2.3548 (2 CeV(beq,beqv))^0.5); ''2 (2 Ln(2))^0.5 = 2.3548

#m_argu num_density
	If_Prm_Eqn_Rpt(num_density, num_densityv, min 0 max 1, del 0.001)
	fit_obj = -4 Pi X CeV(num_density,num_densityv);

pdf_gauss_fwhm = 2.3548 ((beq1 + beq2) / (8 Pi^2))^0.5;

#m_argu xo
   #m_argu a
   #m_argu fwhm
   If_Prm_Eqn_Rpt(xo, xov, )
   If_Prm_Eqn_Rpt(a, av,  min 1e-5)
   If_Prm_Eqn_Rpt(fwhm, fwhmv,  min 0.001 del 0.001)   
   fit_obj = CeV(a,av) Gauss(CeV(xo,xov), CeV(fwhm,fwhmv));

#m_argu xo
   #m_argu a
   #m_argu fwhm
   If_Prm_Eqn_Rpt(xo, xov, )
   If_Prm_Eqn_Rpt(a, av, )
   If_Prm_Eqn_Rpt(fwhm, fwhmv,  min 0.001 del 0.001)   
   fit_obj = CeV(a,av) Gauss(CeV(xo,xov), CeV(fwhm,fwhmv)); ''2 (2 Ln(2))^0.5 = 2.3548

#m_argu c ' size of whole hat
	  #m_argu num_hats
      If_Prm_Eqn_Rpt(c, v, min 1e-6)
		pdf_convolute = num_hats; min_X = -CeV(c, v) 0.5; max_X = CeV(c, v) 0.5;

pkshape_dQ_alpha(dQ,dQv,!alpha,0,lor,lorv)

'Bragg peak shape of the form (dQ^2 + alpha^2 Q^2)^0.5 for 2theta data
	'	Q = (2 Pi/D_spacing)
	'	deltaQ/Q = (dQ^2 + alpha^2 Q^2)^0.5 / Q
	'	delta2Th = 2 deltaQ/Q Tan(Th)
	#m_argu dQ 		'Fixed peak width, in Q
	#m_argu alpha 	'Q dependent peak width term
	#m_argu lor		'Lorch fraction
	If_Prm_Eqn_Rpt(dQ, dQv, min 0.001 max 0.1, del 0.0001)
	If_Prm_Eqn_Rpt(alpha, alphav, min 0.00001 max 0.01, del 0.001)
	If_Prm_Eqn_Rpt(lor, lorv, min 0 max 1)
	peak_type pv
		pv_lor = CeV(lor,lorv);		
		pv_fwhm = 2 (CeV(dQ,dQv)^2 + CeV(alpha,alphav)^2 (2 Pi/D_spacing)^2)^0.5 / (2 Pi/D_spacing)  Tan(Th) Rad;

#m_argu r
	If_Prm_Eqn_Rpt(r, rv, min 1 max 1000, del = 0.01 Val;)
	scale_phase_X = 
	IF X > 2 CeV(r,rv) THEN
		0
	ELSE
		If(X>0.01,(Pi X^2 ((0.25 (X/CeV(r,rv))^3)-(3 X/CeV(r,rv))+4)) / (4 Pi X^2 1),1)
	ENDIF;

0 0 0

1 0 0