As published in J. Appl. Cryst. 48(1) (2015), 189-194, Ectors et al. http://doi.org/10.1107/S1600576714026557 and J. Appl. Cryst. 48(6) (2015), 1998-2001, Ectors et al. http://doi.org/10.1107/S1600576715018488
Updated with an optional plotting macro PlotAnisoCS for the normals_plot window (Version 5+) July 2016
Fixed a bug concerning the automatic orientation tool (Abs(alr) instead of alr). Thank you Remo for the feedback. January 2018
For modifications to get error reporting see (updated) link https://topas.awh.durham.ac.uk/flarum/public/d/538-report-error-with-phase-out/. Thank you Adam for this. Feb 2024
Since my website is down heres a link to a quickstart guide. https://www.dropbox.com/scl/fi/ics0cza5smwv32c1he7g1/AnisoCS-Macro-collection-for-TOPAS.pdf?rlkey=fwmmw5u4x1lo351403zy3zndt&dl=0
Also an example including a .pro file compatible with TOPAS V7 https://www.dropbox.com/scl/fi/rgqisnkk4p8ho65wzydv1/PlotAnisoCS-Example.zip?rlkey=t52322nhzok234y2hmzxy823l&st=u4qwf812&dl=0
Application e.g. for a cylinder in a hexagonal/trigonal system:
str 'or hkl_Is ... AnisoCS( 2, 0, 0, 1, 1, 0, 0, rx, 5, rx, 5, rz, 10, !pc, 0, !tc, 0, !nc, 0, 0) AnisoCSg( tau, 1) AnisoCSgout(result.txt) PlotAnisoCS 'V5+ only ...
macro AnisoCS( mod, h11, k11, l11, h22, k22, l22, ac, av, bc, bv, cc, cv, pc, pv, tc, tv, nc, nv, sp)
{
' As published in Ectors et al. (2015) J. Appl. Cryst. 48, 189-194
' mod which model to use 1 for ellipsoid 2 for elliptic cylinder 3 for cuboid
' h11 k11 l11 vector defining the z-axis
' h22 k22 l22 vector defining the x-axis
' av bv cv rx ry rz of the model
' pv tv nv additional rotation parameters theta' theta' nu
' sp activate special method for triclinic systems 1/0 on/off
' --------- Code starts here ------
peak_buffer_step 0
'Initialize variables
#m_argu mod
#m_argu ac If_Prm_Eqn_Rpt(ac, av, min 0.0001 max 10000)
#m_argu bc If_Prm_Eqn_Rpt(bc, bv, min 0.0001 max 10000)
#m_argu cc If_Prm_Eqn_Rpt(cc, cv, min 0.0001 max 10000)
#m_argu pc If_Prm_Eqn_Rpt(pc, pv, min 0 max 180)
#m_argu tc If_Prm_Eqn_Rpt(tc, tv, min 0 max 180)
#m_argu nc If_Prm_Eqn_Rpt(nc, nv, min 0 max 360)
'G-Metrictensor just for output
local !m11 = Get(a)^2; local !m12 = Get(a)*Get(b)*Cos(Get(ga)*Deg); local !m13 = Get(a)*Get(c)*Cos(Get(be)*Deg);
local !m21 = m12; local !m22 = Get(b)^2; local !m23 = Get(b)*Get(c)*Cos(Get(al)*Deg);
local !m31 = m13; local !m32 = m23; local !m33 = Get(c)^2;
'G*-Reciprocal metrictensor and reciprocal lattice parameters
local !vol = Get(a)*Get(b)*Get(c)*Sqrt(1-((Cos(Get(al)*Deg))^2)-((Cos(Get(be)*Deg))^2)-((Cos(Get(ga)*Deg))^2)+2*(Cos(Get(al)*Deg)*Cos(Get(be)*Deg)*Cos(Get(ga)*Deg)));
local !ar = (Get(b)*Get(c)*Sin(Get(al)*Deg))/vol;
local !br = (Get(a)*Get(c)*Sin(Get(be)*Deg))/vol;
local !cr = (Get(b)*Get(a)*Sin(Get(ga)*Deg))/vol;
local !alr = (Cos(Get(be)*Deg)*Cos(Get(ga)*Deg)-Cos(Get(al)*Deg))/(Sin(Get(be)*Deg)*Sin(Get(ga)*Deg));
local !ber = (Cos(Get(al)*Deg)*Cos(Get(ga)*Deg)-Cos(Get(be)*Deg))/(Sin(Get(al)*Deg)*Sin(Get(ga)*Deg));
local !gar = (Cos(Get(be)*Deg)*Cos(Get(al)*Deg)-Cos(Get(ga)*Deg))/(Sin(Get(be)*Deg)*Sin(Get(al)*Deg));
local !r11 = ar^2; local !r12 = ar*br*gar; local !r13 = ar*cr*ber;
local !r21 = r12; local !r22 = br^2; local !r23 = br*cr*alr;
local !r31 = r13; local !r32 = r23; local !r33 = cr^2;
'Test if special method triclinic is activated
#m_ifarg sp 1
local !h1 = 0; local !h2 = If(Abs(alr) <= 0.00001, 0 , 1);
local !k1 = 0; local !k2 = If(Abs(alr) <= 0.00001, 1 , -(r31 / r32));
local !l1 = 1; local !l2 = 0;
#m_else
local !h1 = h11; local !h2 = h22;
local !k1 = k11; local !k2 = k22;
local !l1 = l11; local !l2 = l22;
#m_endif
'|hkl|Prim aka z-axis
local rowh1 = (h1 r11 + k1 r12 + l1 r13)*h1;
local rowk1 = (h1 r21 + k1 r22 + l1 r23)*k1;
local rowl1 = (h1 r31 + k1 r32 + l1 r33)*l1;
local dvec1 = Sqrt(rowh1+rowk1+rowl1);
'|hkl|Sec aka x-axis
local rowh2 = (h2 r11 + k2 r12 + l2 r13)*h2;
local rowk2 = (h2 r21 + k2 r22 + l2 r23)*k2;
local rowl2 = (h2 r31 + k2 r32 + l2 r33)*l2;
local dvec2 = Sqrt(rowh2+rowk2+rowl2);
'|hkl|Current
local rowH = (H r11 + K r12 + L r13)*H;
local rowK = (H r21 + K r22 + L r23)*K;
local rowL = (H r31 + K r32 + L r33)*L;
local dvec3 = Sqrt(rowH+rowK+rowL);
'Angle phi
local rowg1 = (H r11 + K r12 + L r13)*h1;
local rowg2 = (H r21 + K r22 + L r23)*k1;
local rowg3 = (H r31 + K r32 + L r33)*l1;
local sum1 = rowg1+rowg2+rowg3;
local cosphi = (sum1/(dvec1*dvec3));
local sinphi = Sqrt(Abs(1-(cosphi^2)));
'Angle theta
local rowp1 = (H r11 + K r12 + L r13)*h2;
local rowp2 = (H r21 + K r22 + L r23)*k2;
local rowp3 = (H r31 + K r32 + L r33)*l2;
local sum2 = rowp1+rowp2+rowp3;
local coseta =(sum2/(dvec2*dvec3));
local costheta = If(sinphi <= 0.00001, 1, coseta/sinphi);
local sintheta = Sqrt(Abs(1-(costheta^2)));
'Handedness
local p = h1*(k2*L - l2*K) + k1*(l2*H - h2*L) + l1*(h2*K - k2*H);
local pre = If (p <= 0 , -1 , 1);
'Rotation: Rotationmatrix and calculation of phirot and thetarot
local grot11 = (Cos(CeV(tc,tv) * Deg)*Cos(CeV(pc,pv) * Deg)*Cos(CeV(nc,nv) * Deg))-(Sin(CeV(tc,tv) * Deg)*Sin(CeV(nc,nv) * Deg));
local grot21 = -(Cos(CeV(tc,tv) * Deg)*Cos(CeV(pc,pv) * Deg)*Sin(CeV(nc,nv) * Deg))-(Sin(CeV(tc,tv) * Deg)*Cos(CeV(nc,nv) * Deg));
local grot31 = Cos(CeV(tc,tv) * Deg)*Sin(CeV(pc,pv) * Deg);
local grot12 = (Sin(CeV(tc,tv) * Deg)*Cos(CeV(pc,pv) * Deg)*Cos(CeV(nc,nv) * Deg))+(Cos(CeV(tc,tv) * Deg)*Sin(CeV(nc,nv) * Deg));
local grot22 = -(Sin(CeV(tc,tv) * Deg)*Cos(CeV(pc,pv) * Deg)*Sin(CeV(nc,nv) * Deg))+(Cos(CeV(tc,tv) * Deg)*Cos(CeV(nc,nv) * Deg));
local grot32 = Sin(CeV(tc,tv) * Deg)*Sin(CeV(pc,pv) * Deg);
local grot13 = -(Sin(CeV(pc,pv) * Deg)*Cos(CeV(nc,nv) * Deg));
local grot23 = Sin(CeV(pc,pv) * Deg)*Sin(CeV(nc,nv) * Deg);
local grot33 = Cos(CeV(pc,pv) * Deg);
local cosphirot = (costheta*sinphi*grot31)+(pre*sintheta*sinphi*grot32)+(cosphi*grot33);
local sinphirot = Sqrt(Abs(1-(cosphirot^2)));
local xfac = (costheta*sinphi*grot11)+(pre*sintheta*sinphi*grot12)+(cosphi*grot13);
local costhetarot = If(sinphirot <= 0.0001, 1, xfac / sinphirot);
local sinthetarot = Sqrt(Abs(1-(costhetarot^2)));
'Models:
'Ellipsoid
#m_ifarg mod 1
local volop = ((4/3) 3.14159265358979 CeV(ac,av) CeV(bc,bv) CeV(cc,cv));
local areaop = 3.14159265358979 Sqrt(((CeV(cc,cv) sinphirot)^2)(((CeV(ac,av) sinthetarot)^2)+(CeV(bc,bv) costhetarot)^2)+((CeV(ac,av) CeV(bc,bv) cosphirot)^2));
local sizeL = volop / areaop;
#m_endif
'Cylinder
#m_ifarg mod 2
local volop = (3.14159265358979 CeV(ac,av) CeV(bc,bv) (2 * CeV(cc,cv)));
local areaop = (Sqrt( ((CeV(ac,av) sinthetarot)^2) + ((CeV(bc,bv) costhetarot)^2)) 4 * CeV(cc,cv) * Abs(sinphirot)) + (3.14159265358979 CeV(ac,av) CeV(bc,bv) Abs(cosphirot));
local sizeL = volop / areaop;
#m_endif
'Cuboid
#m_ifarg mod 3
local volop = (8 CeV(ac,av) CeV(bc,bv) CeV(cc,cv));
local areaop = 4 CeV(ac,av) CeV(bc,bv) Abs(cosphirot) + 4 CeV(ac,av) CeV(cc,cv) Abs(sinphirot) Abs(sinthetarot) + 4 CeV(bc,bv) CeV(cc,cv) Abs(costhetarot) Abs(sinphirot) ;
local sizeL = volop / areaop;
#m_endif
'For Output purposes only
local ax = CeV(ac,av);
local bx = CeV(bc,bv);
local cx = CeV(cc,cv);
local px = CeV(pc,pv);
local tx = CeV(tc,tv);
local nx = CeV(nc,nv);
local modx = mod;
local vole = volop;
local larea = volop^(1/3);
local kk = larea / sizeL ;
local areaoe = areaop;
local phi = If(cosphi < -1, 180, If (cosphi > 1, 0, ArcCos(cosphi)*Rad));
local the = If(costheta < -1, 180, If (costheta > 1, 0, ArcCos(costheta)*Rad));
local phirot = If(cosphirot < -1, 180, If (cosphirot > 1, 0, ArcCos(cosphirot)*Rad));
local therot = If(costhetarot < -1, 180, If (costhetarot > 1, 0, ArcCos(costhetarot)*Rad));
'CS_Macro
'lor_fwhm = 0.1 180/pi^2 * Lambda / Cos(Th) * sizeL ... 0.1 for nm->Angstroem
lor_fwhm = 0.1 18.2378130556208 Lam / (Cos(Th) sizeL);
}
macro AnisoCSout(file) {
' As published in Ectors et al. (2015) J. Appl. Cryst. 48, 189-194
out file
Out_String("Analysis Report: ")
Out(Get(phase_name),"%s \n")
Out(Get (r_wp), " Rwp: %8.3f \n")
Out(Get (r_exp), " Rexp:%8.3f \n\n")
Out_String("G-Metrictensor \n\n")
Out(m11,"%8.3f") Out(m12,"%8.3f") Out(m13,"%8.3f\n")
Out(m21,"%8.3f") Out(m22,"%8.3f") Out(m23,"%8.3f\n")
Out(m31,"%8.3f") Out(m32,"%8.3f") Out(m33,"%8.3f\n\n")
Out_String("G*-Metrictensor \n\n")
Out(r11,"%8.3f") Out(r12,"%8.3f") Out(r13,"%8.3f\n")
Out(r21,"%8.3f") Out(r22,"%8.3f") Out(r23,"%8.3f\n")
Out(r31,"%8.3f") Out(r32,"%8.3f") Out(r33,"%8.3f\n\n")
Out_String("Rotationmatrix \n\n")
Out(grot11,"%8.3f") Out(grot12,"%8.3f") Out(grot13,"%8.3f\n")
Out(grot21,"%8.3f") Out(grot22,"%8.3f") Out(grot23,"%8.3f\n")
Out(grot31,"%8.3f") Out(grot32,"%8.3f") Out(grot33,"%8.3f\n\n")
Out_String("Model parameters \n")
Out_String(" (1) Ellipsoid / (2) Cylinder / (3) Cuboid \n")
Out(modx, " Model: %1.0f \n\n")
Out(h1," z-axis hkl: %8.3f") Out(k1," %8.3f") Out(l1," %8.3f \n")
Out(px, " Rotation phi: %8.3f") Out(tx, " theta: %8.3f \n\n")
Out(h2," x-axis hkl: %8.3f") Out(k2," %8.3f") Out(l2," %8.3f \n")
Out(nx, " Additional rotation nu: %8.3f \n\n")
Out(ax," rx-radius: %8.3f nm \n")
Out(bx," ry-radius: %8.3f nm \n")
Out(cx," rz-radius: %8.3f nm \n")
Out(vole," Volume: %8.3f nm^3 \n")
Out(larea," <True CS>: %8.3f nm \n\n")
Out_String(" H K L M 2Th D phi theta phirot thetarot <L_surf>\n")
phase_out file append load out_record out_fmt out_eqn
{
"%4.0f" = H;
"%4.0f" = K;
"%4.0f" = L;
"%4.0f" = M;
"%8.3f" = 2 Th Rad;
"%8.3f" = D_spacing;
"%8.3f" = phi;
"%8.3f" = the;
"%8.3f" = phirot;
"%8.3f" = therot;
"%8.3f\n" = sizeL;
}
}
macro AnisoCSg( tauc, tauv)
{
#m_argu tauc If_Prm_Eqn_Rpt(tauc, tauv, min 1 max 2)
'Parameters Eq.(5-7)
local para = Abs((costhetarot*sinphirot)/ax);
local parb = Abs((sinthetarot*sinphirot)/bx);
local parc = Abs((cosphirot)/cx);
'Parameters chapter 2.4
local pard = Sqrt(para^2+parb^2);
local pare = pard/parc;
local parf = Sqrt(1-pare^2);
'Here comes the tricky part
'Ellipsoid Eq.(4)
local ev = If(modx < 2,CeV(tauc,tauv) * (9/8) * sizeL,
'Cylinder
If(modx < 3,
'Eq. (11)
If(pard <= 0.0001,CeV(tauc, tauv) * sizeL,
'Eq. (12)
If(pard < parc, CeV(tauc, tauv) * ( (2/parc) - (8/(3.14159265358979*pard)) * (-(2/3)+(2/3)*parf+((pare*pare*parf)/3)+(pare*ArcSin(pare)))+
((2*parc)/(3.14159265358979*pard*pard))*(((pare*parf)/2)-((1/2)*(ArcSin(pare)))+(pare*pare*pare*parf)+(2*pare*pare*ArcSin(pare)))),
'Eq. (13)
CeV(tauc, tauv)*((16)/(3*3.14159265358979*pard))-((2*parc)/(4*pard*pard)))),
'Cuboid
'Eq. (8)
If(And(para>=parb, para>=parc),CeV(tauc, tauv)*(2/para)*(1-(parb/(3*para))-(parc/(3*para))+((parb*parc)/(6*para*para))),
'Eq. (9)
If(And(parb>=para, parb>=parc),CeV(tauc, tauv)*(2/parb)*(1-(para/(3*parb))-(parc/(3*parb))+((para*parc)/(6*parb*parb))),
'Eq. (10)
CeV(tauc, tauv)*(2/parc)*(1-(para/(3*parc))-(parb/(3*parc))+((para*parb)/(6*parc*parc))))))
);
'IB calculations Eq. (1-3)
local ibV = 1/ev;
local ibL = 1/(2*sizeL);
'ibG is 0 if ibV < ibL
local ibG = If (ibV > ibL, ((Sqrt((4*ibV^2)-(4.27679864*ibV*ibL)+(0.27679864*ibL^2)))/2), 0);
'CS_G macro
local sizeG = ibG/(Sqrt(3.14159265358979/(4*Ln(2))));
gauss_fwhm = 0.1 Rad Lam sizeG / Cos(Th);
'For Output purposes only
local taux = CeV(tauc, tauv);
}
macro AnisoCSgout(file) {
out file
Out_String("Analysis Report: ")
Out(Get(phase_name),"%s \n")
Out(Get (r_wp), " Rwp: %8.3f \n")
Out(Get (r_exp), " Rexp:%8.3f \n\n")
Out_String("G-Metrictensor \n\n")
Out(m11,"%8.3f") Out(m12,"%8.3f") Out(m13,"%8.3f\n")
Out(m21,"%8.3f") Out(m22,"%8.3f") Out(m23,"%8.3f\n")
Out(m31,"%8.3f") Out(m32,"%8.3f") Out(m33,"%8.3f\n\n")
Out_String("G*-Metrictensor \n\n")
Out(r11,"%8.3f") Out(r12,"%8.3f") Out(r13,"%8.3f\n")
Out(r21,"%8.3f") Out(r22,"%8.3f") Out(r23,"%8.3f\n")
Out(r31,"%8.3f") Out(r32,"%8.3f") Out(r33,"%8.3f\n\n")
Out_String("Rotationmatrix \n\n")
Out(grot11,"%8.3f") Out(grot12,"%8.3f") Out(grot13,"%8.3f\n")
Out(grot21,"%8.3f") Out(grot22,"%8.3f") Out(grot23,"%8.3f\n")
Out(grot31,"%8.3f") Out(grot32,"%8.3f") Out(grot33,"%8.3f\n\n")
Out_String("Model parameters \n")
Out_String(" (1) Ellipsoid / (2) Cylinder / (3) Cuboid \n")
Out(modx, " Model: %1.0f \n\n")
Out(h1," z-axis hkl: %8.3f") Out(k1," %8.3f") Out(l1," %8.3f \n")
Out(px, " Rotation phi: %8.3f") Out(tx, " theta: %8.3f \n\n")
Out(h2," x-axis hkl: %8.3f") Out(k2," %8.3f") Out(l2," %8.3f \n")
Out(nx, " Additional rotation nu: %8.3f \n\n")
Out(ax," rx-radius_surf: %8.3f nm ")
Out(ax*taux,"/ rx-radius_vol: %8.3f nm \n")
Out(bx," ry-radius_surf: %8.3f nm ")
Out(bx*taux,"/ ry-radius_vol: %8.3f nm \n")
Out(cx," rz-radius_surf: %8.3f nm ")
Out(cx*taux,"/ rz-radius_vol: %8.3f nm \n")
Out(vole," Volume_surf: %8.3f nm^3 ")
Out(vole*taux*taux*taux,"/ Volume_vol: %8.3f nm^3 \n")
Out(larea," <True CS>_surf: %8.3f nm ")
Out(larea*taux,"/ <True CS>_vol: %8.3f nm \n\n")
Out_String("Distribution parameters \n")
Out(taux," Tau: %8.4f \n")
Out_String(" If Lognormal of <True CS>\n")
Out(Ln(larea*taux)-(3.5*Ln(taux))," Logn. mean:%8.4f \n")
Out(Sqrt(Ln(taux))," Logn. var.: %8.4f \n\n")
Out_String(" H K L M 2Th D <L_surf><L_vol> <L_vol>/<L_surf> \n")
phase_out file append load out_record out_fmt out_eqn
{
"%4.0f" = H;
"%4.0f" = K;
"%4.0f" = L;
"%4.0f" = M;
"%8.3f" = 2 Th Rad;
"%8.3f" = D_spacing;
"%8.3f" = sizeL;
"%8.3f" = ev;
"%8.3f\n" = ev/sizeL;
}
}
macro PlotAnisoCS
{
'ellipsoid
normals_plot = If(modx < 2, sizeL*(3/4),
'cylinder
If(modx < 3,
1/( ( ((Abs((1/cx)*cosphirot))^(10))+((Abs((sinphirot))^(10))*( ((Abs((1/bx)*sinthetarot))^(2)) + ((Abs((1/ax)*costhetarot))^(2)) )^(5) ) )^(1/10) ),
'cuboid
1/( ( ((Abs((1/cx)*cosphirot))^(10))+((Abs((sinphirot))^(10))*( ((Abs((1/bx)*sinthetarot))^(10)) + ((Abs((1/ax)*costhetarot))^(10)) ) ) )^(1/10) )
)
);
}