Page 69 out of 81 total pages , Page 12 out of 12 pages in this chapter
[7] J.E. Wertz and J.R. Bolton, Electron Spin Resonance. Elementary Theory and Practical Applications, McGraw-Hill Book Co., New York, New York, (1986), Chapman and Hall.
[8] Brink, D.M. and Satchler, G.R. (1962), Angular Momentum, Clarendon Press, Oxford.
/* IntGu_LC0.cc
**************************************************************-*-c++-*-
**
**
** Test Program for the GAMMA Library
**
** G Interaction Literature Comparison 0
**
**
**
** This program checks the G interaction class IntG in
**
** GAMMA. In particular it looks to see how well the class parallels
**
** the articles by Pascal P. Man
**
**
**
** "G Interactions", Encyclopedia of Magnetic Resonance,
**
** by Grant and Harris, Vol 6, Ped-Rel, page 3838-3869.
**
**
**
** and Alexander Vega
**
**
**
** "G Nuclei in Solids", Encyclopedia of Magnetic Resonance,
**
** by Grant and Harris, Vol 6, Ped-Rel, page 3869-3889.
**
**
**
** In particular, their PAS G Hamiltonians are generated and
**
** compared with the G interaction class Hamiltonians.
**
**
**
** Man's Hamiltonians will be generated from equations in (5) on page
**
** 3839 of the his article. Vega's Hamiltonians will be made from
**
** equations (28), (32) and (33) of his article. Note that his (32)
**
** is missing a factor of 1/3 on the <1|H|3> and <3|H|1> components.
**
**
**
** Author: S.A. Smith
**
** Date: 10/11/96
**
** Update: 10/11/96
**
** Version: 3.6
**
** Copyright: S. Smith. You can modify this program for personal use, but
**
** you must leave it intact if you re-distribute it.
**
**
**
**********************************************************************
*******************/
#include <gamma.h> // Include GAMMA
main (int argc, char* argv[])
{
// Set Up The G Interaction
int qn=1;
double I;
query_parameter(argc, argv, qn++, // Read in the coupling
"\n\tSpin Guantum Number? ", I);
double W; // G frequency
query_parameter(argc, argv, qn++, // Read in the coupling
"\n\tG Frequency(kHz)? ", W);
W *= 1.e3; // Put this in Hz
double Geta;
query_parameter(argc, argv, qn++, // Read in the coupling
"\n\tG Asymmetry [0, 1]? ", Geta);
// Construct GAMMA G Interaction
IntG G(I,wG2GCC(W,I),Geta,0.0,0.0);
// Here are the Operators To Build Man's Hamiltonians
int Ival = int(2.*I + 1); // For 1 spin SOp functions
matrix IE = Ie(Ival); // The operator 1
matrix IM = Im(Ival); // The operator I-
matrix IP = Ip(Ival); // The operator I+
matrix IZ = Iz(Ival); // The operator Iz
matrix IX = Ix(Ival); // The operator Ix
matrix IY = Iy(Ival); // The operator Iy
// Here's The H Accoring To Man's Equation (5a)
// (Note That His W is Half Of Our Definition)
matrix HMa = 3.0*IZ*IZ - (I*(I+1))*IE + Geta*((IX*IX)-(IY*IY));
HMa *= (W/6.0);
// Here's The H Accoring To Man'c Equation (5c)
// (Note That His W is Half Of Our Definition)
matrix HMb = 3.0*IZ*IZ - (I*(I+1))*IE + (Geta/2.)*((IP*IP)+(IM*IM));
HMb *= (W/6.0);
// Here's The H According To GAMMA
matrix HG = G.H();
// Here's The H Also According To GAMMA
matrix HGB = G.H(0.0, 0.0);
// Here Are Vegas V's According To Equations (22-27, 31)
// (Switches eta Sign To Account For Opposite PAS Definition)
double Eta = -G.eta();
double Vxx = 0.5*(- 1. - Eta);
double Vyy = 0.5*(-1. + Eta);
double Vzz = 1.0;
double Vxy = 0.0;
double Vxz = 0.0;
double Vyz = 0.0;
complex V1(-Vxz, -Vyz);
complex Vm1(Vxz, -Vyz);
complex V2(0.5*(Vxx-Vyy), Vxy);
complex Vm2(0.5*(Vxx-Vyy), -Vxy);
// Generate H According To Vega's Equations (32) Or (33)
matrix HVega;
if(I == 1)
{
HVega = matrix(3,3);
HVega.put(Vzz/6.0, 0, 0);
HVega.put(Vm1/sqrt(2.0), 0, 1);
HVega.put(Vm2/3., 0, 2); // Added 1/3 Factor!
HVega.put(-V1/sqrt(2.0), 1, 0);
HVega.put(-Vzz/3.0, 1, 1);
HVega.put(-Vm1/sqrt(2.0), 1, 2);
HVega.put(V2/3., 2, 0); // Added 1/3 Factor!
HVega.put(V1/sqrt(2.0), 2, 1);
HVega.put(Vzz/6.0, 2, 2);
HVega *= G.wG();
}
else if(I == 1.5)
{
HVega = matrix(4,4);
HVega.put(Vzz/2.0, 0, 0);
HVega.put(Vm1/sqrt(3.0), 0, 1);
HVega.put(Vm2/sqrt(3.0), 0, 2);
HVega.put(0.0, 0, 3);
HVega.put(-V1/sqrt(3.0), 1, 0);
HVega.put(-Vzz/2.0, 1, 1);
HVega.put(0.0, 1, 2);
HVega.put(Vm2/sqrt(3.0), 1, 3);
HVega.put(V2/sqrt(3.0), 2, 0);
HVega.put(0.0, 2, 1);
HVega.put(-Vzz/2.0, 2, 2);
HVega.put(-Vm1/sqrt(3.0), 2, 3);
HVega.put(0.0, 3, 0);
HVega.put(V2/sqrt(3.0), 3, 1);
HVega.put(V1/sqrt(3.0), 3, 2);
HVega.put(Vzz/2.0, 3, 3);
HVega *= G.wG();
}
// Generate H According To Vega's Equation (28)
matrix HV = Vzz*(3.*IZ*IZ-(I*(I+1.))*IE);
HV += (Vxx-Vyy)*(IX*IX-IY*IY);
HV += 2*Vxy*(IX*IY-IY*IX);
HV += 2*Vxz*(IX*IZ-IZ*IX);
HV += 2*Vyz*(IY*IZ-IZ*IY);
HV *= G.wG()/6.0;
// Output the Results for Visual Comparison
cout << "\n\t\t\t\tGAMMA's G H:\t" << HG;
cout << "\n\t\t\t\tGAMMA's Other G H:\t" << HGB;
cout << "\n\t\t\t\tMan's G H(a):\n\t" << HMa;
cout << "\n\t\t\t\tMan's G H(b):\n\t" << HMb;
if(I == 1.0 || I == 1.5)
cout << "\n\t\t\t\tVega's G H:\n\t" << HVega;
cout << "\n\t\t\t\tVega's Generic Guad H:\n\t" << HV;
}
/* IntGu_LC1.cc
**************************************************************-*-c++-*-
**
**
** Test Program for the GAMMA Library
**
** G Interaction Literature Comparison 1
**
**
**
** This program checks the G interaction class IntG in
**
** GAMMA. In particular it looks to see how well the class parallels
**
** the article by Alexander Vega -
**
**
**
** "G Nuclei in Solids", Encyclopedia of Magnetic Resonance,
**
** by Grant and Harris, Vol 6, Ped-Rel, page 3869-3889.
**
**
**
** Specifcally, herein we generate the spatial tensor components of
**
** an oriented G interaction and and compare the results to
**
** A.Vega's equations (22-27) and 31 on pages 3884-3885.
**
**
**
** Author: S.A. Smith
**
** Date: 10/11/96
**
** Update: 10/11/96
**
** Version: 3.6
**
** Copyright: S. Smith. You can modify this program as you see fit
**
** for personal use, but you must leave the program intact
**
** if you re-distribute it.
**
**
**
**********************************************************************
*******************/
#include <gamma.h> // Include GAMMA
main (int argc, char* argv[])
{
// Construct A G Interaction
int qn=1;
double W; // G frequency
query_parameter(argc, argv, qn++, // Read in the coupling
"\n\tG Frequency(kHz)? ", W);
W *= 1.e3; // Put this in Hz
double Geta;
query_parameter(argc, argv, qn++, // Read in the coupling
"\n\tG Asymmetry [0, 1]? ", Geta);
double Gtheta, Gphi;
query_parameter(argc, argv, qn++, // Read in the angle
"\n\tAngle down from z [0, 180]? ", Gtheta);
query_parameter(argc, argv, qn++, // Read in the angle
"\n\tAngle over from x [0, 360]? ", Gphi);
double I=1.0; // Use I=1, but this doesn't
double GCC = wG2GCC(W, I); // Heres quad. coupling
IntG G(I,GCC,Geta,Gtheta,Gphi); // matter for spatial parts
// Here Are Vegas V's According To Equations (22-27, 31)
// Note We Change Sign On ETA As He Using A Different PAS Definition
double Theta = G.theta()*DEG2RAD;
double Phi = G.phi()*DEG2RAD;
double Eta = -G.eta();
double Stheta = sin(Theta);
double Ctheta = cos(Theta);
double C2phi = cos(2.*Phi);
double S2phi = sin(2.*Phi);
double Vxx = 0.5*(3.*Stheta*Stheta - 1. - Eta*Ctheta*Ctheta*C2phi);
double Vxy = 0.5*Eta*Ctheta*S2phi;
double Vxz = -0.5*(Stheta*Ctheta*(3.0 + Eta*C2phi));
double Vyx = Vxy;
double Vyy = 0.5*(-1. + Eta*C2phi);
double Vyz = 0.5*Eta*Stheta*S2phi;
double Vzx = Vxz;
double Vzy = Vyz;
double Vzz = 0.5*(3.*Ctheta*Ctheta - 1. - Eta*Stheta*Stheta*C2phi);
complex V0(sqrt(1.5)*Vzz);
complex V1(-Vxz, -Vyz);
complex Vm1(Vxz, -Vyz);
complex V2(0.5*(Vxx-Vyy), Vxy);
complex Vm2(0.5*(Vxx-Vyy), -Vxy);
// Here Are The A's According To GAMMA G Interaction
// Need To Scale Our A's By (1/2)/sqrt[5/(24*PI)] To Get Vega's V's
double X = 0.5/RT5O24PI;
double Thetad = G.theta();
double Phid = G.phi();
double AGxx = X*G.Axx(Thetad, Phid);
double AGxy = X*G.Axy(Thetad, Phid);
double AGxz = X*G.Axz(Thetad, Phid);
double AGyy = X*G.Ayy(Thetad, Phid);
double AGyx = X*G.Ayx(Thetad, Phid);
double AGyz = X*G.Ayz(Thetad, Phid);
double AGzz = X*G.Azz(Thetad, Phid);
double AGzx = X*G.Azx(Thetad, Phid);
double AGzy = X*G.Azy(Thetad, Phid);
// Here Are The A's According To GAMMA G Interaction
// Need To Scale Our A's By (1/2)/sqrt[5/(24*PI)] To Get Vega's V's
double AG1xx = X*G.Axx();
double AG1xy = X*G.Axy();
double AG1xz = X*G.Axz();
double AG1yy = X*G.Ayy();
double AG1yx = X*G.Ayx();
double AG1yz = X*G.Ayz();
double AG1zz = X*G.Azz();
double AG1zx = X*G.Azx();
double AG1zy = X*G.Azy();
// Here Are The A's According To GAMMA G Interaction
// (Note That space_T Uses Azz>=Ayy>=Axx So ETA Opposite Vega's)
space_T Agen = A2(0.0, 1.0, Geta);
Agen = Agen.rotate(Phid, Thetad, 0.0);
Cartesian(Agen);
// Output Everyone For A Visual Comparison
cout << "\n " << " Vega" << " IntGA"
<< " IntGB" << " space_T";
cout << "\nVxx " << form("%8.3f", Vxx) << " " << form("%8.3f", AGxx)
<< " " << form("%8.3f", AG1xx) << " " << form("%8.3f", Agen.Ccomponent(0,0));
cout << "\nVxy " << form("%8.3f", Vxy) << " " << form("%8.3f", AGxy)
<< " " << form("%8.3f", AG1xy) << " " << form("%8.3f", Agen.Ccomponent(0,1));
cout << "\nVxz " << form("%8.3f", Vxz) << " " << form("%8.3f", AGxz)
<< " " << form("%8.3f", AG1xz) << " " << form("%8.3f", Agen.Ccomponent(0,2));
cout << "\nVyy " << form("%8.3f", Vyy) << " " << form("%8.3f", AGyy)
<< " " << form("%8.3f", AG1yy) << " " << form("%8.3f", Agen.Ccomponent(1,1));
cout << "\nVyx " << form("%8.3f", Vyx) << " " << form("%8.3f", AGyx)
<< " " << form("%8.3f", AG1yx) << " " << form("%8.3f", Agen.Ccomponent(1,0));
cout << "\nVyz " << form("%8.3f", Vyz) << " " << form("%8.3f", AGyz)
<< " " << form("%8.3f", AG1yz) << " " << form("%8.3f", Agen.Ccomponent(1,2));
cout << "\nVzz " << form("%8.3f", Vzz) << " " << form("%8.3f", AGzz)
<< " " << form("%8.3f", AG1zz) << " " << form("%8.3f", Agen.Ccomponent(2,2));
cout << "\nVzx " << form("%8.3f", Vzx) << " " << form("%8.3f", AGzx)
<< " " << form("%8.3f", AG1zx) << " " << form("%8.3f", Agen.Ccomponent(2,0));
cout << "\nVzy " << form("%8.3f", Vzy) << " " << form("%8.3f", AGzy)
<< " " << form("%8.3f", AG1zy) << " " << form("%8.3f", Agen.Ccomponent(2,1));
cout << "\nV0 " << V0 << " " << X*G.A0(Thetad, Phid)
<< " " << X*G.A0() << " " << Agen.component(2,0);
cout << "\nV1 " << V1 << " " << X*G.A1(Thetad, Phid)
<< " " << X*G.A1() << " " << Agen.component(2,1);
cout << "\nV-1" << Vm1 << " " << X*G.Am1(Thetad, Phid)
<< " " << X*G.Am1() << " " << Agen.component(2,-1);
cout << "\nV2 " << V2 << " " << X*G.A2(Thetad, Phid)
<< " " << X*G.A2() << " " << Agen.component(2,2);
cout << "\nV-2" << Vm2 << " " << X*G.Am2(Thetad, Phid)
<< " " << X*G.Am2() << " " << Agen.component(2,-2);
cout << "\n\n\n";
}
/* IntGu_PCT0.cc
************************************************************-*-c++-*-
**
**
** Example Program for the GAMMA Library
**
**
**
** This program calculates a powder average for a single spin which
**
** is associated with a G interaction. The high field
**
** approximation is invoked in that the G Hamiltonian is
**
** treated as a perturbation to the Zeeman Hamiltonian and taken to
**
** second order. Only G Hamiltonian terms which are
**
** rotationally invariant about the field axis (z) are maintained.
**
** Furthermore, only the central transtion will be considered.
**
**
**
** This will program is similar to IntGu_Pow2.cc but restricts the
**
** computation to only the central transition. In turn, that means
**
** only spins with I=m*1/2 where m is odd and larger than 1 are valid.
**
** Analog formula will be used to construct the spectrum.
**
**
**
** Later version of GAMMA will have the functions "scale" and "sum"
**
** in the library itself, so you will need to remove them from this
**
** program in that event.
**
**
**
** Author: S.A. Smith
**
** Date: 10/15/96
**
** Update: 10/15/96
**
** Version: 3.6
**
** Copyright: S. Smith. You can modify this program as you see fit
**
** for personal use, but you must leave the program intact
**
** if you re-distribute it.
**
**
**
**********************************************************************
*******************/
#include <gamma.h> // Include GAMMA
void addW(row_vector& vx, double Fst, double Ffi, double F, double I)
// Fst : Frequency of 1st point of vx (Hz)
// Ffi : Frequency of last point of vx (Hz)
// F : Transition frequency (Hz)
// Output void : The transition is added to the
{
if(F<Fst || F>Ffi) return; // Insure its in range
double Nm1 = double(vx.size()-1); // Freq. -> point conversion
double m = Nm1/(Ffi-Fst); // Slope Freq. -> pt
double dpt = m*(F-Fst); // Point index of F
int pt = int(dpt); // Main point for F
double drem = dpt - pt; // Part which isn't
if(!drem) vx.put(vx.get(pt)+I, pt); // Add if on a point
else if(drem > 0) // If in between points
{ // then just split it up
vx.put(vx.get(pt)+(1.0-drem)*I, pt); // between the two
vx.put(vx.get(pt+1)+drem*I, pt+1);
}
else
{
vx.put(vx.get(pt)+(1.0+drem)*I, pt);
vx.put(vx.get(pt-1)-drem*I, pt-1);
}
return;
}
main (int argc, char* argv[])
{
cout << "\n\t\tG Central Transition Powder Pattern";
cout << "\n\t\t (131Xe:3/2, 55Mn:5/2, 51V:7/2, ...)\n";
// First Make A G Interaction
String Iso; // Isotope of spin
int qn=1; // Guery index
query_parameter(argc, argv, qn++, // Get the isotope type
"\n\tIsotope Type [131Xe, 55Mn, 51V, ....]? ", Iso);
double wG; // Set Guad. frequency
query_parameter(argc, argv, qn++, // Get the G coupling
"\n\tG Frequency (kHz)? ", wG);
wG *= 1.e3; // Switch to Hz
double eta; // Set Guad. frequency
query_parameter(argc, argv, qn++, // Get the G coupling
"\n\tG eta Value [0, 1]? ", eta);
double Om;
query_parameter(argc, argv, qn++, // Get the field strength
"\n\tLarmor Frequency (MHz)? ", Om);
Om *= 1.e6; // Switch to MHz
Isotope S(Iso); // Make a spin isotope
double I = S.qn(); // This is isotope I value
IntG G(I,wG2GCC(wG, I), eta); // Set a Guad interaction
if(!int(2*I)%2)
{
cout << "\n\n\tSorry, I Must Be m*1/2, m Odd!\n\n";
exit(-1);
}
// Set Things Up For The Powder Average
int npts = 4096; // Block size
int Ntheta, Nphi=0; // Angle increment counts
query_parameter(argc, argv, qn++, // Get the theta increments
"\n\t# Theta (z down) Increments Spanning [0, 180]? ", Ntheta);
if(eta)
query_parameter(argc, argv, qn++, // Get the phi increments
"\n\t# Phi (x over) Increments Spanning [0, 360)? ", Nphi);
matrix ABC = G.wGcentral(Ntheta, Nphi); // Prep. for 2nd order shifts
// Powder Averaging
// Angle theta Is Down From The +z Axis, Angle phi Over From +x
// Note that since the 2nd order shift Wcentral(theta, phi) is symmetric with
// respect to both angles we need only average over parts of both angle ranges.
// For theta this means we sum the results from angles [0, 90) + half the result
// at 90. Twice that sum would produce the total theta average over [0, 180].
// For phi we usually average [0, 360) so this is reduced to a sum over 3/4 the
// result at 0 + the results from angles (0, 90) + 1/2 the result at 90. Four
// times that sum would produce the total phi average over [0, 360).
double dthe, Nm1o2 = double(Ntheta-1.0)/2.0; // For powder average
double dphi, Nm2o4 = double(Nphi)/4.0; // For powder average
int theta, phi; // Orientation angles
double W, WG = G.wG(); // Base Guad. frequency
double Ifact = I*(I+1) - 0.75; // Part of the prefactor
double prefact = -WG*WG*Ifact/Om; // Majority of the prefact
double Aaxis = (-1.0/9.0)*prefact; // For plot scaling
double Ctheta, Stheta, Cthetasq, Ctheta4; // We'll need these
double Fst = -2.5*Aaxis; // Starting plot limit
double Ffi = 1.5*Aaxis; // End plot limit
row_vector data(npts, complex0); // Array for spectrum
for(theta=0; theta<Ntheta; theta++) // Loop over theta angles
{
dthe = double(theta);
if(dthe <= Nm1o2) // Only look upper half
{ // of the sphere
Ctheta = ABC.getRe(0,theta); // Scale factor cos(theta)
Stheta = ABC.getRe(1,theta); // Scale factor sin(theta)
Cthetasq = Ctheta*Ctheta; // cosine(theta)^2
if(dthe == Nm1o2) Stheta *= 0.5; // Half scale if theta=90
if(!eta) // Without eta, no phi
{ // averaging is needed
W=(prefact/16.)*(1.-Cthetasq)*(9.*Cthetasq-1.); // Here is W adjustment
addW(data, Fst, Ffi, W, Stheta); // Add transition to spectrum
}
else
{
cout.flush();
Ctheta4 = Cthetasq*Cthetasq; // cosine(theta)^4
for(phi=0; phi<Nphi; phi++) // Loop over phi angles
{
dphi = double(phi); // Phi index as double
if(dphi <= Nm2o4) // Only sum 1st quarter
{
if(!phi) Stheta *= 0.75; // 3/4 scale if phi=0
else if(dphi == Nm2o4) Stheta *= 0.5; // 1/2 scale if phi=90
W = ABC.getRe(2,phi)*Ctheta4; // A part of W
W += ABC.getRe(3,phi)*Cthetasq; // B part of W
W += ABC.getRe(4,phi); // C part of W
W *= (prefact/6.); // Scale
addW(data, Fst, Ffi, W, Stheta); // Add transition to spectrum
}
}
}
}
}
double lb = 40.0; // Set a line broading factor
cout << "\n\n\tDone With Discrete Powder Average. Processing...";
cout.flush();
data = IFFT(data); // Put into time domain
exponential_multiply(data,-lb); // Apodize the "FID"
data = FFT(data); // Put back into frequency domain
GP_1D("spec.asc", data, 0, -2.5, 1.5); // Output the points in ASCII
GP_1Dplot("spec.gnu", "spec.asc"); // Call Gnuplot and make plot now
}
|
GAMMA Support Provided by the National High Magnetic Field Laboratory
© 1996 Scott A. Smith, The NHMFL, and The Florida State University. All Rights Reserved. |