/********************************************************************************************************************************** * * McXtrace, X-ray tracing package * Copyright, All rights reserved * DTU Physics, Kgs. Lyngby, Denmark * Synchrotron SOLEIL, Saint-Aubin, France * * Component: Polycrystal * * %Identification * Written by: Martin Cramer Pedersen (mcpe@nbi.dk) * * Based on code by: Erik Knudsen, Kenichi Oikawa, and Alberto Cereser * Date: January 2015 * Origin: University of Copenhagen * Release: McXtrace 1.0 * * Polycrystal made from single crystal-like voxels * * %Description * The component creates a threedimensional grid of cubic instances of the Single_crystal-component,through * which the rays are propagated. The component relies on a list of possible orientations and stretches of the * initial unit cell and a list correlating each voxel of the polycrystal to an entry in the list of rotations * and stretches. * * Example: Polycrystal( MapFile= "polycrystal_1layer_2orts.map", OrientationsFile= "stretch_2orts.orts", * ReflectionsDatafile= "GeReduced.lau", xwidth= 200e-6, yheight= 200e-6, zdepth = 50e-6, * DeltadOverd = 0.001, Mosaicity = 1, SigmaAbsorbtion = 0.0, SigmaIncoherent = 0.0, * MaxNumberOfReflections = 1, ProbabilityOfTransmission = 0.5, * ax = 5.6579, ay = 0.0000, az = 0.0000, bx = 0.0000, by = 5.6579, bz = 0.0000, cx = 0.0000, cy = 0.0000, cz = 5.6579 ) * * %Parameters * MapFile: [str] File describing, which orientation is found in which voxel * OrientationsFile: [str] File describing the different orientations * ReflectionsDatafile: [str] File describing the reflections in the relevant lattice (usually in .lau-format) * MaterialDatafile: [str] File describing the scattering and absorbtion properties of the material * xwidth: [m] Width of the sample * yheight: [m] Height of the sample * zdepth: [m] Depth of the sample * ax: [AA] x-coordinate of the first internal unit cell vector in the sample * ay: [AA] y-coordinate of the first internal unit cell vector in the sample * az: [AA] z-coordinate of the first internal unit cell vector in the sample * bx: [AA] x-coordinate of the second internal unit cell vector in the sample * by: [AA] y-coordinate of the second internal unit cell vector in the sample * bz: [AA] z-coordinate of the second internal unit cell vector in the sample * cx: [AA] x-coordinate of the third internal unit cell vector in the sample * cy: [AA] y-coordinate of the third internal unit cell vector in the sample * cz: [AA] z-coordinate of the third internal unit cell vector in the sample * Mosaicity: [moa] Gaussian mosaicity * MosaicityA: [moa] Anisotropic mosaicity around the first unit cell vector * MosaicityB: [moa] Anisotropic mosaicity around the second unit cell vector * MosaicityC: [moa] Anisotropic mosaicity around the third unit cell vector * DeltadOverd: [1] Statistical description of the lattice spacing * ProbabilityOfTransmission: [0-1] Probability that a ray will not interact with the sample * SigmaAbsorbtion: [fm^2] Absorbtion crosssection of the sample * SigmaIncoherent: [fm^2] Incoherent crosssection of the sample * MaxNumberOfReflections: [1] Highest number of allowed scattering events in the entire crystal - to prevent computationally expensive high-order multiple scattering. If this parameter is set to 0, all possible orders of scattering are considered. * Reciprocal: [0/1] If this parameter is set to 0, then the lattice vectors should be given in real space. Anything else implies that the vectors are given in reciprocal space. * verbose: [0/1] If nonzero - output more info to the console. * %End ***********************************************************************************************************************************/ DEFINE COMPONENT Polycrystal SETTING PARAMETERS(string MapFile = "", string OrientationsFile = "", string ReflectionsDatafile = "Si.lau", string MaterialDatafile = "Si.txt", xwidth = 0.0, yheight = 0.0, zdepth = 0.0, ax = 5.43, ay = 0.00, az = 0.00, bx = 0.00, by = 5.43, bz = 0.00, cx = 0.00, cy = 0.00, cz = 5.43, Mosaicity = 3.0, MosaicityA = 0.0, MosaicityB = 0.0, MosaicityC = 0.0, DeltadOverd = 0.01, ProbabilityOfTransmission = 0.01, SigmaAbsorbtion = 0.0, SigmaIncoherent = 0.0, MaxNumberOfReflections = 1, Reciprocal = 0, verbose = 0) SHARE %{ // External libraries %include "read_table-lib" #ifndef MXPX_REFL_SLIST_SIZE #define MXPX_REFL_SLIST_SIZE 128 #endif // Structs for hkl-data struct hklDataStruct { // Miller indices int h; int k; int l; // Structure factor double F2; // Coordinates in reciprocal space double Tau[3]; double TauLength; // Local coordinate system double u1[3]; double u2[3]; double u3[3]; // Gaussians double Sigma[3]; double TotalSigma; double m[3]; double GaussianCutoff; }; struct TauDataStruct { // Scattering properties double Reflectivity; double Crosssection; // Index in reflection table int Index; // The following vectors are in local koordinates // Initial wave vector double ki[3]; // Rho = ki - tau double Rho[3]; double RhoLength; // Origin of the Ewald sphere tangent plane double Origin[3]; // Normal vector of the Ewald sphere tanget plane double NormalVector[3]; // Vectors spanning the Ewald sphere tangent plane double b1[3]; double b2[3]; // Cholesky decomposition of the matrix L double L11; double L12; double L22; double DeterminantL; // Center of Gaussian on tangent plane double GaussCenterx; double GaussCentery; }; struct hklStruct { // List of reflections struct hklDataStruct* hklData; int NumberOfReflections; // Format of columns in list: [h, k, l, F, F2] int ColumnOrder[5]; // Flags used to raise warnings int Recip; // Sigma of delta d / d double Deltad_d; // Coordinates of unit cell vectors double LatticeVectora[3]; double LatticeVectorb[3]; double LatticeVectorc[3]; double LatticeVectoraLength; double LatticeVectorbLength; double LatticeVectorcLength; double ReciprocalLatticeVectora[3]; double ReciprocalLatticeVectorb[3]; double ReciprocalLatticeVectorc[3]; // Lattice parameters and angles double LatticeAngleA; double LatticeAngleB; double LatticeAngleC; // Absorbtion crosssection double SigmaAbs; // Incoherent crosssection double SigmaInc; // Unit cell volume double VolumeOfUnitCell; }; // Struct for absorption-data struct AbsorbtionStruct { int muc; t_Table Table; }; // I/O Function for hkl-data int ReadReflectionData (char* ReflectionsDatafile, struct hklStruct* hklInfo, double Mosaicity, double MosaicityA, double MosaicityB, double MosaicityC, double* MosaicityAB) { // Declarations int i; int NumberOfRows; int PrintOrientations = 0; // Structs for storing read-in info struct hklDataStruct* hklData; t_Table sTable; // Temporary variables double Temp[3]; char** Parsing; int FlagMissingFile = 0; int NumberOfAtoms = 1; // hkl-info double asLength; double bsLength; double csLength; double b1[3]; double b2[3]; // Mosaicity double ConvertedMosaicityA; double ConvertedMosaicityB; double ConvertedMosaicityC; // Vectors used to calculate anisotropic mosaicity double XiA[3]; double XiALength; double XiB[3]; double XiBLength; double XiC[3]; double XiCLength; double Jacobiann[3]; double Jacobianv[3]; // Variables used to calculate in-plane anisotropic mosaicity double c1; double c2; double Determinant; double SigmaTauC; double em[3]; double TauA[3]; double TauB[3]; // Read file if (!ReflectionsDatafile || !strlen (ReflectionsDatafile) || !strcmp (ReflectionsDatafile, "NULL") || !strcmp (ReflectionsDatafile, "0")) { hklInfo->NumberOfReflections = 0; FlagMissingFile = 1; } if (!FlagMissingFile) { if (Table_Read (&sTable, ReflectionsDatafile, 1) <= 0 || sTable.columns < 4) { fprintf (stderr, "Error (%s): The number of columns in %s should be at least %d for [h, k, l, F2]\n", "ReadReflectionData", ReflectionsDatafile, 4); return (-1); } if (!sTable.rows) { fprintf (stderr, "Error (%s): The number of rows in %s should be at least %d\n", "ReadReflectionData", 1); return (-1); } else { NumberOfRows = sTable.rows; } // Parsing of header Parsing = Table_ParseHeader (sTable.header, "sigma_abs", "sigma_a", "sigma_inc", "sigma_i", "column_h", "column_k", "column_l", "column_F", "column_F2", "Delta_d/d", "lattice_a", "lattice_b", "lattice_c", "lattice_aa", "lattice_bb", "lattice_cc", "nb_atoms", "multiplicity", NULL); if (Parsing) { if (Parsing[0] && !hklInfo->SigmaAbs) hklInfo->SigmaAbs = atof (Parsing[0]); if (Parsing[1] && !hklInfo->SigmaAbs) hklInfo->SigmaAbs = atof (Parsing[1]); if (Parsing[2] && !hklInfo->SigmaInc) hklInfo->SigmaInc = atof (Parsing[2]); if (Parsing[3] && !hklInfo->SigmaInc) hklInfo->SigmaInc = atof (Parsing[3]); if (Parsing[4]) hklInfo->ColumnOrder[0] = atoi (Parsing[4]); if (Parsing[5]) hklInfo->ColumnOrder[1] = atoi (Parsing[5]); if (Parsing[6]) hklInfo->ColumnOrder[2] = atoi (Parsing[6]); if (Parsing[7]) hklInfo->ColumnOrder[3] = atoi (Parsing[7]); if (Parsing[8]) hklInfo->ColumnOrder[4] = atoi (Parsing[8]); if (Parsing[9] && hklInfo->Deltad_d < 0) hklInfo->Deltad_d = atof (Parsing[9]); if (Parsing[10] && !hklInfo->LatticeVectoraLength) hklInfo->LatticeVectoraLength = atof (Parsing[10]); if (Parsing[11] && !hklInfo->LatticeVectorbLength) hklInfo->LatticeVectorbLength = atof (Parsing[11]); if (Parsing[12] && !hklInfo->LatticeVectorcLength) hklInfo->LatticeVectorcLength = atof (Parsing[12]); if (Parsing[13] && !hklInfo->LatticeAngleA) hklInfo->LatticeAngleA = atof (Parsing[13]); if (Parsing[14] && !hklInfo->LatticeAngleB) hklInfo->LatticeAngleB = atof (Parsing[14]); if (Parsing[15] && !hklInfo->LatticeAngleC) hklInfo->LatticeAngleC = atof (Parsing[15]); if (Parsing[16]) NumberOfAtoms = atof (Parsing[16]); if (Parsing[17]) NumberOfAtoms = atof (Parsing[17]); for (i = 0; i < 18; ++i) { if (Parsing[i]) { free (Parsing[i]); } } free (Parsing); } } // Assign variables if (NumberOfAtoms > 1) { hklInfo->SigmaAbs *= NumberOfAtoms; hklInfo->SigmaInc *= NumberOfAtoms; } // Special cases for the structure definition if (hklInfo->LatticeVectora[0] || hklInfo->LatticeVectora[1] || hklInfo->LatticeVectora[2]) { hklInfo->LatticeVectoraLength = 0; } if (hklInfo->LatticeVectorb[0] || hklInfo->LatticeVectorb[1] || hklInfo->LatticeVectorb[2]) { hklInfo->LatticeVectorbLength = 0; } if (hklInfo->LatticeVectorc[0] || hklInfo->LatticeVectorc[1] || hklInfo->LatticeVectorc[2]) { hklInfo->LatticeVectorcLength = 0; } // Compute the norm from vector a if missing if (hklInfo->LatticeVectora[0] || hklInfo->LatticeVectora[1] || hklInfo->LatticeVectora[2]) { asLength = sqrt (hklInfo->LatticeVectora[0] * hklInfo->LatticeVectora[0] + hklInfo->LatticeVectora[1] * hklInfo->LatticeVectora[1] + hklInfo->LatticeVectora[2] * hklInfo->LatticeVectora[2]); if (!hklInfo->LatticeVectorb[0] && !hklInfo->LatticeVectorb[1] && !hklInfo->LatticeVectorb[2]) { hklInfo->LatticeVectoraLength = hklInfo->LatticeVectorbLength = asLength; } if (!hklInfo->LatticeVectorc[0] && !hklInfo->LatticeVectorc[1] && !hklInfo->LatticeVectorc[2]) { hklInfo->LatticeVectoraLength = hklInfo->LatticeVectorcLength = asLength; } } if (hklInfo->LatticeVectoraLength && !hklInfo->LatticeVectorbLength) { hklInfo->LatticeVectorbLength = hklInfo->LatticeVectoraLength; } if (hklInfo->LatticeVectorbLength && !hklInfo->LatticeVectorcLength) { hklInfo->LatticeVectorcLength = hklInfo->LatticeVectorbLength; } // Compute the lattive angles if not set from data file. Not used when in vector mode if (hklInfo->LatticeVectoraLength && !hklInfo->LatticeAngleA) { hklInfo->LatticeAngleA = 90; } if (hklInfo->LatticeAngleA && !hklInfo->LatticeAngleB) { hklInfo->LatticeAngleB = hklInfo->LatticeAngleA; } if (hklInfo->LatticeAngleB && !hklInfo->LatticeAngleC) { hklInfo->LatticeAngleC = hklInfo->LatticeAngleB; } // Parameters consistency checks if (!hklInfo->LatticeVectora[0] && !hklInfo->LatticeVectora[1] && !hklInfo->LatticeVectora[2] && !hklInfo->LatticeVectoraLength) { fprintf (stderr, "Error (%s): Wrong a lattice vector definition.\n", "ReadReflectionData"); return (0); } if (!hklInfo->LatticeVectorb[0] && !hklInfo->LatticeVectorb[1] && !hklInfo->LatticeVectorb[2] && !hklInfo->LatticeVectorbLength) { fprintf (stderr, "Error (%s): Wrong b lattice vector definition.\n", "ReadReflectionData"); return (0); } if (!hklInfo->LatticeVectorc[0] && !hklInfo->LatticeVectorc[1] && !hklInfo->LatticeVectorc[2] && !hklInfo->LatticeVectorcLength) { fprintf (stderr, "Error (%s): Wrong c lattice vector definition.\n", "ReadReflectionData"); return (0); } if (hklInfo->LatticeAngleA && hklInfo->LatticeAngleB && hklInfo->LatticeAngleC && hklInfo->Recip) { fprintf (stderr, "Error (%s): Selecting reciprocal cell and angles is unmeaningful.\n", "ReadReflectionData"); return (0); } // When lengths a, b, c and angles are given (instead of vectors a, b, c) if (hklInfo->LatticeAngleA && hklInfo->LatticeAngleB && hklInfo->LatticeAngleC) { if (hklInfo->LatticeVectoraLength) { asLength = hklInfo->LatticeVectoraLength; } else { asLength = sqrt (pow (hklInfo->LatticeVectora[0], 2) + pow (hklInfo->LatticeVectora[1], 2) + pow (hklInfo->LatticeVectora[2], 2)); } if (hklInfo->LatticeVectorbLength) { bsLength = hklInfo->LatticeVectorbLength; } else { bsLength = sqrt (pow (hklInfo->LatticeVectorb[0], 2) + pow (hklInfo->LatticeVectorb[1], 2) + pow (hklInfo->LatticeVectorb[2], 2)); } if (hklInfo->LatticeVectorcLength) { csLength = hklInfo->LatticeVectorcLength; } else { csLength = sqrt (pow (hklInfo->LatticeVectorc[0], 2) + pow (hklInfo->LatticeVectorc[1], 2) + pow (hklInfo->LatticeVectorc[2], 2)); } hklInfo->LatticeVectorb[2] = asLength; hklInfo->LatticeVectorb[1] = 0.0; hklInfo->LatticeVectorb[0] = 0.0; hklInfo->LatticeVectora[2] = bsLength * cos (hklInfo->LatticeAngleC * DEG2RAD); hklInfo->LatticeVectora[1] = bsLength * sin (hklInfo->LatticeAngleC * DEG2RAD); hklInfo->LatticeVectora[0] = 0.0; hklInfo->LatticeVectorc[2] = csLength * cos (hklInfo->LatticeAngleB * DEG2RAD); hklInfo->LatticeVectorc[1] = csLength * (cos (hklInfo->LatticeAngleA * DEG2RAD) - cos (hklInfo->LatticeAngleC * DEG2RAD) * cos (hklInfo->LatticeAngleB * DEG2RAD)) / sin (hklInfo->LatticeAngleC * DEG2RAD); hklInfo->LatticeVectorc[0] = sqrt (pow (csLength, 2) - pow (hklInfo->LatticeVectorc[2], 2) - pow (hklInfo->LatticeVectorc[1], 2)); if (PrintOrientations) { fprintf (stderr, "INFO (%s): \nStructure: \na = %g b = %g c = %g \naa = %g bb = %g cc = %g \n", "ReadReflectionData", asLength, bsLength, csLength, hklInfo->LatticeAngleA, hklInfo->LatticeAngleB, hklInfo->LatticeAngleC); } } else { if (!hklInfo->Recip) { if (PrintOrientations) { fprintf (stderr, "INFO (%s): \nStructure; \na = [%g, %g, %g] \nb = [%g, %g, %g] \nc = [%g, %g, %g] \n", "ReadReflectionData", hklInfo->LatticeVectora[0], hklInfo->LatticeVectora[1], hklInfo->LatticeVectora[2], hklInfo->LatticeVectorb[0], hklInfo->LatticeVectorb[1], hklInfo->LatticeVectorb[2], hklInfo->LatticeVectorc[0], hklInfo->LatticeVectorc[1], hklInfo->LatticeVectorc[2]); } } else { if (PrintOrientations) { fprintf (stderr, "INFO (%s): \nStructure: \na* = [%g, %g, %g] \nb* = [%g, %g, %g] \nc* = [%g, %g, %g] \n", "ReadReflectionData", hklInfo->LatticeVectora[0], hklInfo->LatticeVectora[1], hklInfo->LatticeVectora[2], hklInfo->LatticeVectorb[0], hklInfo->LatticeVectorb[1], hklInfo->LatticeVectorb[2], hklInfo->LatticeVectorc[0], hklInfo->LatticeVectorc[1], hklInfo->LatticeVectorc[2]); } } } // Compute reciprocal or direct lattice vectors if (!hklInfo->Recip) { vec_prod (Temp[0], Temp[1], Temp[2], hklInfo->LatticeVectorb[0], hklInfo->LatticeVectorb[1], hklInfo->LatticeVectorb[2], hklInfo->LatticeVectorc[0], hklInfo->LatticeVectorc[1], hklInfo->LatticeVectorc[2]); hklInfo->VolumeOfUnitCell = fabs (scalar_prod (hklInfo->LatticeVectora[0], hklInfo->LatticeVectora[1], hklInfo->LatticeVectora[2], Temp[0], Temp[1], Temp[2])); hklInfo->ReciprocalLatticeVectora[0] = 2.0 * M_PI / hklInfo->VolumeOfUnitCell * Temp[0]; hklInfo->ReciprocalLatticeVectora[1] = 2.0 * M_PI / hklInfo->VolumeOfUnitCell * Temp[1]; hklInfo->ReciprocalLatticeVectora[2] = 2.0 * M_PI / hklInfo->VolumeOfUnitCell * Temp[2]; vec_prod (Temp[0], Temp[1], Temp[2], hklInfo->LatticeVectorc[0], hklInfo->LatticeVectorc[1], hklInfo->LatticeVectorc[2], hklInfo->LatticeVectora[0], hklInfo->LatticeVectora[1], hklInfo->LatticeVectora[2]); hklInfo->ReciprocalLatticeVectorb[0] = 2.0 * M_PI / hklInfo->VolumeOfUnitCell * Temp[0]; hklInfo->ReciprocalLatticeVectorb[1] = 2.0 * M_PI / hklInfo->VolumeOfUnitCell * Temp[1]; hklInfo->ReciprocalLatticeVectorb[2] = 2.0 * M_PI / hklInfo->VolumeOfUnitCell * Temp[2]; vec_prod (Temp[0], Temp[1], Temp[2], hklInfo->LatticeVectora[0], hklInfo->LatticeVectora[1], hklInfo->LatticeVectora[2], hklInfo->LatticeVectorb[0], hklInfo->LatticeVectorb[1], hklInfo->LatticeVectorb[2]); hklInfo->ReciprocalLatticeVectorc[0] = 2.0 * M_PI / hklInfo->VolumeOfUnitCell * Temp[0]; hklInfo->ReciprocalLatticeVectorc[1] = 2.0 * M_PI / hklInfo->VolumeOfUnitCell * Temp[1]; hklInfo->ReciprocalLatticeVectorc[2] = 2.0 * M_PI / hklInfo->VolumeOfUnitCell * Temp[2]; } else { hklInfo->ReciprocalLatticeVectora[0] = hklInfo->LatticeVectora[0]; hklInfo->ReciprocalLatticeVectora[1] = hklInfo->LatticeVectora[1]; hklInfo->ReciprocalLatticeVectora[2] = hklInfo->LatticeVectora[2]; hklInfo->ReciprocalLatticeVectorb[0] = hklInfo->LatticeVectorb[0]; hklInfo->ReciprocalLatticeVectorb[1] = hklInfo->LatticeVectorb[1]; hklInfo->ReciprocalLatticeVectorb[2] = hklInfo->LatticeVectorb[2]; hklInfo->ReciprocalLatticeVectorc[0] = hklInfo->LatticeVectorc[0]; hklInfo->ReciprocalLatticeVectorc[1] = hklInfo->LatticeVectorc[1]; hklInfo->ReciprocalLatticeVectorc[2] = hklInfo->LatticeVectorc[2]; // Compute the direct cell parameters for completeness vec_prod (Temp[0], Temp[1], Temp[2], hklInfo->ReciprocalLatticeVectorb[0] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectorb[1] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectorb[2] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectorc[0] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectorc[1] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectorc[2] / (2 * M_PI)); hklInfo->VolumeOfUnitCell = 1.0 / fabs (scalar_prod (hklInfo->ReciprocalLatticeVectora[0] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectora[1] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectora[2] / (2.0 * M_PI), Temp[0], Temp[1], Temp[2])); hklInfo->LatticeVectora[0] = Temp[0] * hklInfo->VolumeOfUnitCell; hklInfo->LatticeVectora[1] = Temp[1] * hklInfo->VolumeOfUnitCell; hklInfo->LatticeVectora[2] = Temp[2] * hklInfo->VolumeOfUnitCell; vec_prod (Temp[0], Temp[1], Temp[2], hklInfo->ReciprocalLatticeVectorc[0] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectorc[1] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectorc[2] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectora[0] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectora[1] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectora[2] / (2 * M_PI)); hklInfo->LatticeVectorb[0] = Temp[0] * hklInfo->VolumeOfUnitCell; hklInfo->LatticeVectorb[1] = Temp[1] * hklInfo->VolumeOfUnitCell; hklInfo->LatticeVectorb[2] = Temp[2] * hklInfo->VolumeOfUnitCell; vec_prod (Temp[0], Temp[1], Temp[2], hklInfo->ReciprocalLatticeVectora[0] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectora[1] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectora[2] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectorb[0] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectorb[1] / (2.0 * M_PI), hklInfo->ReciprocalLatticeVectorb[2] / (2 * M_PI)); hklInfo->LatticeVectorc[0] = Temp[0] * hklInfo->VolumeOfUnitCell; hklInfo->LatticeVectorc[1] = Temp[1] * hklInfo->VolumeOfUnitCell; hklInfo->LatticeVectorc[2] = Temp[2] * hklInfo->VolumeOfUnitCell; } if (FlagMissingFile) { return (-1); } if (!hklInfo->ColumnOrder[0] || !hklInfo->ColumnOrder[1] || !hklInfo->ColumnOrder[2]) { fprintf (stderr, "Error (%s): Wrong h, k, l column definition\n", "ReadReflectionData"); return (0); } if (!hklInfo->ColumnOrder[3] && !hklInfo->ColumnOrder[4]) { fprintf (stderr, "Error (%s): Wrong F, F2 column definition\n", "ReadReflectionData"); return (0); } // Allocate hklDataStruct array hklData = (struct hklDataStruct*)malloc (NumberOfRows * sizeof (struct hklDataStruct)); for (i = 0; i < NumberOfRows; i++) { // Get data from table /*hklData[i].h = Table_Index(sTable, i, hklInfo->ColumnOrder[0] - 1);*/ hklData[i].h = Table_Index (sTable, i, 0); /*hklData[i].k = Table_Index(sTable, i, hklInfo->ColumnOrder[1] - 1);*/ hklData[i].k = Table_Index (sTable, i, 1); /*hklData[i].l = Table_Index(sTable, i, hklInfo->ColumnOrder[2] - 1);*/ hklData[i].l = Table_Index (sTable, i, 2); /*if (hklInfo->ColumnOrder[3]) { hklData[i].F2 = pow(Table_Index(sTable, i, hklInfo->ColumnOrder[3] - 1), 2); } else if (hklInfo->ColumnOrder[4]) { hklData[i].F2 = Table_Index(sTable, i, hklInfo->ColumnOrder[4] - 1); }*/ hklData[i].F2 = Table_Index (sTable, i, 6); // Precompute some values hklData[i].Tau[0] = hklData[i].h * hklInfo->ReciprocalLatticeVectora[0] + hklData[i].k * hklInfo->ReciprocalLatticeVectorb[0] + hklData[i].l * hklInfo->ReciprocalLatticeVectorc[0]; hklData[i].Tau[1] = hklData[i].h * hklInfo->ReciprocalLatticeVectora[1] + hklData[i].k * hklInfo->ReciprocalLatticeVectorb[1] + hklData[i].l * hklInfo->ReciprocalLatticeVectorc[1]; hklData[i].Tau[2] = hklData[i].h * hklInfo->ReciprocalLatticeVectora[2] + hklData[i].k * hklInfo->ReciprocalLatticeVectorb[2] + hklData[i].l * hklInfo->ReciprocalLatticeVectorc[2]; hklData[i].TauLength = sqrt (pow (hklData[i].Tau[0], 2) + pow (hklData[i].Tau[1], 2) + pow (hklData[i].Tau[2], 2)); hklData[i].u1[0] = hklData[i].Tau[0] / hklData[i].TauLength; hklData[i].u1[1] = hklData[i].Tau[1] / hklData[i].TauLength; hklData[i].u1[2] = hklData[i].Tau[2] / hklData[i].TauLength; hklData[i].Sigma[0] = FWHM2RMS * hklInfo->Deltad_d * hklData[i].TauLength; // Find two arbitrary axes perpendicular to tau and each other normal_vec (&(b1[0]), &(b1[1]), &(b1[2]), hklData[i].u1[0], hklData[i].u1[1], hklData[i].u1[2]); vec_prod (b2[0], b2[1], b2[2], hklData[i].u1[0], hklData[i].u1[1], hklData[i].u1[2], b1[0], b1[1], b1[2]); // Find the two mosaic axes perpendicular to tau if (Mosaicity > 0) { // Use isotropic mosaic hklData[i].u2[0] = b1[0]; hklData[i].u2[1] = b1[1]; hklData[i].u2[2] = b1[2]; hklData[i].Sigma[1] = FWHM2RMS * hklData[i].TauLength * MIN2RAD * Mosaicity; hklData[i].u3[0] = b2[0]; hklData[i].u3[1] = b2[1]; hklData[i].u3[2] = b2[2]; hklData[i].Sigma[2] = FWHM2RMS * hklData[i].TauLength * MIN2RAD * Mosaicity; } else if (MosaicityA > 0 && MosaicityB > 0 && MosaicityC > 0) { // Use anisotropic mosaic fprintf (stderr, "Warning (%s): You are using an experimental feature: anistropic mosaicity. Please examine your data carefully.\n", "ReadReflectionData"); // Compute the jacobian of (tau_v, tau_n) from rotations around the unit cell vectors struct hklDataStruct* CurrenthklData = &(hklData[i]); // Input parameters are in arc minutes ConvertedMosaicityA = MosaicityA * MIN2RAD; ConvertedMosaicityB = MosaicityB * MIN2RAD; ConvertedMosaicityC = MosaicityC * MIN2RAD; if (hklInfo->LatticeVectoraLength == 0) { hklInfo->LatticeVectoraLength = sqrt (pow (hklInfo->LatticeVectora[0], 2) + pow (hklInfo->LatticeVectora[1], 2) + pow (hklInfo->LatticeVectora[2], 2)); } if (hklInfo->LatticeVectorbLength == 0) { hklInfo->LatticeVectorbLength = sqrt (pow (hklInfo->LatticeVectorb[0], 2) + pow (hklInfo->LatticeVectorb[1], 2) + pow (hklInfo->LatticeVectorb[2], 2)); } if (hklInfo->LatticeVectorcLength == 0) { hklInfo->LatticeVectorcLength = sqrt (pow (hklInfo->LatticeVectorc[0], 2) + pow (hklInfo->LatticeVectorc[1], 2) + pow (hklInfo->LatticeVectorc[2], 2)); } CurrenthklData->u2[0] = b1[0]; CurrenthklData->u2[1] = b1[1]; CurrenthklData->u2[2] = b1[2]; CurrenthklData->u3[0] = b2[0]; CurrenthklData->u3[1] = b2[1]; CurrenthklData->u3[2] = b2[2]; XiA[0] = CurrenthklData->Tau[0] - (M_2_PI * hklData[i].h / hklInfo->LatticeVectoraLength) * hklInfo->ReciprocalLatticeVectora[0]; XiA[1] = CurrenthklData->Tau[1] - (M_2_PI * hklData[i].h / hklInfo->LatticeVectoraLength) * hklInfo->ReciprocalLatticeVectora[1]; XiA[2] = CurrenthklData->Tau[2] - (M_2_PI * hklData[i].h / hklInfo->LatticeVectoraLength) * hklInfo->ReciprocalLatticeVectora[2]; XiB[0] = CurrenthklData->Tau[0] - (M_2_PI * hklData[i].h / hklInfo->LatticeVectorbLength) * hklInfo->ReciprocalLatticeVectorb[0]; XiB[1] = CurrenthklData->Tau[1] - (M_2_PI * hklData[i].h / hklInfo->LatticeVectorbLength) * hklInfo->ReciprocalLatticeVectorb[1]; XiB[2] = CurrenthklData->Tau[2] - (M_2_PI * hklData[i].h / hklInfo->LatticeVectorbLength) * hklInfo->ReciprocalLatticeVectorb[2]; XiC[0] = CurrenthklData->Tau[0] - (M_2_PI * hklData[i].h / hklInfo->LatticeVectorcLength) * hklInfo->ReciprocalLatticeVectorc[0]; XiC[1] = CurrenthklData->Tau[1] - (M_2_PI * hklData[i].h / hklInfo->LatticeVectorcLength) * hklInfo->ReciprocalLatticeVectorc[1]; XiC[2] = CurrenthklData->Tau[2] - (M_2_PI * hklData[i].h / hklInfo->LatticeVectorcLength) * hklInfo->ReciprocalLatticeVectorc[2]; XiALength = sqrt (pow (XiA[0], 2) + pow (XiA[1], 2) + pow (XiA[2], 2)); XiBLength = sqrt (pow (XiB[0], 2) + pow (XiB[1], 2) + pow (XiB[2], 2)); XiCLength = sqrt (pow (XiC[0], 2) + pow (XiC[1], 2) + pow (XiC[2], 2)); vec_prod (Temp[0], Temp[1], Temp[2], CurrenthklData->Tau[0], CurrenthklData->Tau[1], CurrenthklData->Tau[2], CurrenthklData->u2[0], CurrenthklData->u2[1], CurrenthklData->u2[2]); Jacobiann[0] = XiALength / hklInfo->LatticeVectoraLength / CurrenthklData->TauLength * scalar_prod (hklInfo->ReciprocalLatticeVectora[0], hklInfo->ReciprocalLatticeVectora[1], hklInfo->ReciprocalLatticeVectora[2], Temp[0], Temp[1], Temp[2]); vec_prod (Temp[0], Temp[1], Temp[2], CurrenthklData->Tau[0], CurrenthklData->Tau[1], CurrenthklData->Tau[2], CurrenthklData->u2[0], CurrenthklData->u2[1], CurrenthklData->u2[2]); Jacobiann[1] = XiBLength / hklInfo->LatticeVectorbLength / CurrenthklData->TauLength * scalar_prod (hklInfo->ReciprocalLatticeVectorb[0], hklInfo->ReciprocalLatticeVectorb[1], hklInfo->ReciprocalLatticeVectorb[2], Temp[0], Temp[1], Temp[2]); vec_prod (Temp[0], Temp[1], Temp[2], CurrenthklData->Tau[0], CurrenthklData->Tau[1], CurrenthklData->Tau[2], CurrenthklData->u2[0], CurrenthklData->u2[1], CurrenthklData->u2[2]); Jacobiann[2] = XiCLength / hklInfo->LatticeVectorcLength / CurrenthklData->TauLength * scalar_prod (hklInfo->ReciprocalLatticeVectorc[0], hklInfo->ReciprocalLatticeVectorc[1], hklInfo->ReciprocalLatticeVectorc[2], Temp[0], Temp[1], Temp[2]); vec_prod (Temp[0], Temp[1], Temp[2], CurrenthklData->Tau[0], CurrenthklData->Tau[1], CurrenthklData->Tau[2], CurrenthklData->u3[0], CurrenthklData->u3[1], CurrenthklData->u3[2]); Jacobianv[0] = XiALength / hklInfo->LatticeVectoraLength / CurrenthklData->TauLength * scalar_prod (hklInfo->ReciprocalLatticeVectora[0], hklInfo->ReciprocalLatticeVectora[1], hklInfo->ReciprocalLatticeVectora[2], Temp[0], Temp[1], Temp[2]); vec_prod (Temp[0], Temp[1], Temp[2], CurrenthklData->Tau[0], CurrenthklData->Tau[1], CurrenthklData->Tau[2], CurrenthklData->u3[0], CurrenthklData->u3[1], CurrenthklData->u3[2]); Jacobianv[1] = XiBLength / hklInfo->LatticeVectorbLength / CurrenthklData->TauLength * scalar_prod (hklInfo->ReciprocalLatticeVectorb[0], hklInfo->ReciprocalLatticeVectorb[1], hklInfo->ReciprocalLatticeVectorb[2], Temp[0], Temp[1], Temp[2]); vec_prod (Temp[0], Temp[1], Temp[2], CurrenthklData->Tau[0], CurrenthklData->Tau[1], CurrenthklData->Tau[2], CurrenthklData->u3[0], CurrenthklData->u3[1], CurrenthklData->u3[2]); Jacobianv[2] = XiCLength / hklInfo->LatticeVectorcLength / CurrenthklData->TauLength * scalar_prod (hklInfo->ReciprocalLatticeVectorc[0], hklInfo->ReciprocalLatticeVectorc[1], hklInfo->ReciprocalLatticeVectorc[2], Temp[0], Temp[1], Temp[2]); // With the jacobian we can compute the sigmas in terms of the orthogonal vectors u2 and u3 CurrenthklData->Sigma[1] = ConvertedMosaicityA * fabs (Jacobianv[0]) + ConvertedMosaicityB * fabs (Jacobianv[1]) + ConvertedMosaicityC * fabs (Jacobianv[2]); CurrenthklData->Sigma[2] = ConvertedMosaicityA * fabs (Jacobiann[0]) + ConvertedMosaicityB * fabs (Jacobiann[1]) + ConvertedMosaicityC * fabs (Jacobiann[2]); } else if (MosaicityAB[0] != 0 && MosaicityAB[1] != 0) { if ((MosaicityAB[2] == 0 && MosaicityAB[3] == 0 && MosaicityAB[4] == 0) || (MosaicityAB[5] == 0 && MosaicityAB[6] == 0 && MosaicityAB[7] == 0)) { fprintf (stderr, "Error (%s): In-plane mosaics are specified but one (or both) in-plane reciprocal vector is the zero vector.\n", "ReadReflectionData"); return (-1); } fprintf (stderr, "Warning (%s): You are using an experimental feature: \"in-plane\" anistropic mosaicity. Please examine your data carefully.\n", "ReadReflectionData"); struct hklDataStruct* CurrenthklData = &(hklData[i]); // Convert Miller indices to taus if (hklInfo->LatticeVectoraLength == 0) { hklInfo->LatticeVectoraLength = sqrt (hklInfo->LatticeVectora[0] * hklInfo->LatticeVectora[0] + hklInfo->LatticeVectora[1] * hklInfo->LatticeVectora[1] + hklInfo->LatticeVectora[2] * hklInfo->LatticeVectora[2]); } if (hklInfo->LatticeVectorbLength == 0) { hklInfo->LatticeVectorbLength = sqrt (hklInfo->LatticeVectorb[0] * hklInfo->LatticeVectorb[0] + hklInfo->LatticeVectorb[1] * hklInfo->LatticeVectorb[1] + hklInfo->LatticeVectorb[2] * hklInfo->LatticeVectorb[2]); } if (hklInfo->LatticeVectorcLength == 0) { hklInfo->LatticeVectorcLength = sqrt (hklInfo->LatticeVectorc[0] * hklInfo->LatticeVectorc[0] + hklInfo->LatticeVectorc[1] * hklInfo->LatticeVectorc[1] + hklInfo->LatticeVectorc[2] * hklInfo->LatticeVectorc[2]); } TauA[0] = M_2_PI * ((MosaicityAB[2] / hklInfo->LatticeVectoraLength) * hklInfo->ReciprocalLatticeVectora[0] + (MosaicityAB[3] / hklInfo->LatticeVectorbLength) * hklInfo->ReciprocalLatticeVectorb[0] + (MosaicityAB[4] / hklInfo->LatticeVectorcLength) * hklInfo->ReciprocalLatticeVectorc[0]); TauA[1] = M_2_PI * ((MosaicityAB[2] / hklInfo->LatticeVectoraLength) * hklInfo->ReciprocalLatticeVectora[1] + (MosaicityAB[3] / hklInfo->LatticeVectorbLength) * hklInfo->ReciprocalLatticeVectorb[1] + (MosaicityAB[4] / hklInfo->LatticeVectorcLength) * hklInfo->ReciprocalLatticeVectorc[1]); TauA[2] = M_2_PI * ((MosaicityAB[2] / hklInfo->LatticeVectoraLength) * hklInfo->ReciprocalLatticeVectora[2] + (MosaicityAB[3] / hklInfo->LatticeVectorbLength) * hklInfo->ReciprocalLatticeVectorb[2] + (MosaicityAB[4] / hklInfo->LatticeVectorcLength) * hklInfo->ReciprocalLatticeVectorc[2]); TauB[0] = M_2_PI * ((MosaicityAB[5] / hklInfo->LatticeVectoraLength) * hklInfo->ReciprocalLatticeVectora[0] + (MosaicityAB[6] / hklInfo->LatticeVectorbLength) * hklInfo->ReciprocalLatticeVectorb[0] + (MosaicityAB[7] / hklInfo->LatticeVectorcLength) * hklInfo->ReciprocalLatticeVectorc[0]); TauB[1] = M_2_PI * ((MosaicityAB[5] / hklInfo->LatticeVectoraLength) * hklInfo->ReciprocalLatticeVectora[1] + (MosaicityAB[6] / hklInfo->LatticeVectorbLength) * hklInfo->ReciprocalLatticeVectorb[1] + (MosaicityAB[7] / hklInfo->LatticeVectorcLength) * hklInfo->ReciprocalLatticeVectorc[1]); TauB[2] = M_2_PI * ((MosaicityAB[5] / hklInfo->LatticeVectoraLength) * hklInfo->ReciprocalLatticeVectora[2] + (MosaicityAB[6] / hklInfo->LatticeVectorbLength) * hklInfo->ReciprocalLatticeVectorb[2] + (MosaicityAB[7] / hklInfo->LatticeVectorcLength) * hklInfo->ReciprocalLatticeVectorc[2]); // Check determinants to see how we should compute the linear combination of a and b (to match c) if ((Determinant = TauA[0] * TauB[1] - TauA[1] * TauB[0]) != 0) { c1 = (CurrenthklData->Tau[0] * TauB[1] - CurrenthklData->Tau[1] * TauB[0]) / Determinant; c2 = (TauA[0] * CurrenthklData->Tau[1] - TauA[1] * CurrenthklData->Tau[0]) / Determinant; } else if ((Determinant = TauA[1] * TauB[2] - TauA[2] * TauB[1]) != 0) { c1 = (CurrenthklData->Tau[1] * TauB[2] - CurrenthklData->Tau[2] * TauB[1]) / Determinant; c2 = (TauA[1] * CurrenthklData->Tau[2] - TauA[2] * CurrenthklData->Tau[1]) / Determinant; } else if ((Determinant = TauA[0] * TauB[2] - TauA[2] * TauB[0]) != 0) { c1 = (CurrenthklData->Tau[0] * TauB[2] - CurrenthklData->Tau[2] * TauB[0]) / Determinant; c2 = (TauA[0] * CurrenthklData->Tau[2] - TauA[2] * CurrenthklData->Tau[0]) / Determinant; } if ((c1 == 0) && (c2 == 0)) { fprintf (stderr, "WARNING (%s): Reflection tau[i] = (%g, %g, %g) has no component in defined mosaic plane.\n", "ReadReflectionData", CurrenthklData->Tau[0], CurrenthklData->Tau[1], CurrenthklData->Tau[2]); } // Compute linear combination => sig_tau_i = | c1*sig_tau_a + c2*sig_tau_b | - also add in the minute to radian scaling factor SigmaTauC = MIN2RAD * sqrt (pow (c1 * MosaicityAB[0], 2) + pow (c2 * MosaicityAB[1], 2)); CurrenthklData->u2[0] = b1[0]; CurrenthklData->u2[1] = b1[1]; CurrenthklData->u2[2] = b1[2]; CurrenthklData->u3[0] = b2[0]; CurrenthklData->u3[1] = b2[1]; CurrenthklData->u3[2] = b2[2]; // So now let's compute the rotation around planenormal TauA X TauB // g_bar (unit normal of rotation plane) = TauA X TauB / norm(TauA X TauB) vec_prod (Temp[0], Temp[1], Temp[2], TauA[0], TauA[1], TauA[2], TauB[0], TauB[1], TauB[2]); vec_prod (em[0], em[1], em[2], CurrenthklData->Tau[0], CurrenthklData->Tau[1], CurrenthklData->Tau[2], Temp[0], Temp[1], Temp[2]); NORM (em[0], em[1], em[2]); CurrenthklData->Sigma[1] = CurrenthklData->TauLength * SigmaTauC * fabs (scalar_prod (em[0], em[1], em[2], CurrenthklData->u2[0], CurrenthklData->u2[1], CurrenthklData->u2[2])); CurrenthklData->Sigma[2] = CurrenthklData->TauLength * SigmaTauC * fabs (scalar_prod (em[0], em[1], em[2], CurrenthklData->u3[0], CurrenthklData->u3[1], CurrenthklData->u3[2])); // Protect against collapsing gaussians - these seem to be sensible values if (CurrenthklData->Sigma[1] < 1e-5) { CurrenthklData->Sigma[1] = 1e-5; } if (CurrenthklData->Sigma[2] < 1e-5) { CurrenthklData->Sigma[2] = 1e-5; } } else { fprintf (stderr, "ERROR (%s): Either mosaic or (mosaic_a, mosaic_b, mosaic_c) must be given and be > 0.\n", "ReadReflectionData"); return (-1); } hklData[i].TotalSigma = hklData[i].Sigma[0] * hklData[i].Sigma[1] * hklData[i].Sigma[2]; hklData[i].m[0] = 1.0 / (2.0 * pow (hklData[i].Sigma[0], 2)); hklData[i].m[1] = 1.0 / (2.0 * pow (hklData[i].Sigma[1], 2)); hklData[i].m[2] = 1.0 / (2.0 * pow (hklData[i].Sigma[2], 2)); // Set Gauss cutoff to 5 times the maximal sigma if (hklData[i].Sigma[0] > hklData[i].Sigma[1]) { if (hklData[i].Sigma[0] > hklData[i].Sigma[2]) { hklData[i].GaussianCutoff = 5.0 * hklData[i].Sigma[0]; } else { hklData[i].GaussianCutoff = 5.0 * hklData[i].Sigma[2]; } } else { if (hklData[i].Sigma[1] > hklData[i].Sigma[2]) { hklData[i].GaussianCutoff = 5.0 * hklData[i].Sigma[1]; } else { hklData[i].GaussianCutoff = 5.0 * hklData[i].Sigma[2]; } } } Table_Free (&sTable); hklInfo->hklData = hklData; hklInfo->NumberOfReflections = i; return (i); } // Struct for description of absorption int ReadAbsorbtionData (char* AbsorptionFile, struct AbsorbtionStruct* AbsorbtionInfo) { // Declarations int j; int Status; char** Parsing; // Construct table if (AbsorptionFile && strlen (AbsorptionFile) && strcmp (AbsorptionFile, "NULL")) { if ((Status = Table_Read (&(AbsorbtionInfo->Table), AbsorptionFile, 0)) == -1) { fprintf (stderr, "Error (%s): Could not parse file %s\n", "ReadAbsorptionData", AbsorptionFile); exit (-1); } Parsing = Table_ParseHeader (AbsorbtionInfo->Table.header, "Z", "A[r]", "rho", NULL); if (AbsorbtionInfo->Table.columns == 3) { AbsorbtionInfo->muc = 1.0; } else { AbsorbtionInfo->muc = 5.0; } Table_Stat (&(AbsorbtionInfo->Table)); return 1; } else { // Create empty table Table_Init (&(AbsorbtionInfo->Table), 2, 2); AbsorbtionInfo->Table.data[0] = 0.0; AbsorbtionInfo->Table.data[1] = 0.0; AbsorbtionInfo->Table.data[2] = FLT_MAX; AbsorbtionInfo->Table.data[3] = 0.0; fprintf (stderr, "Warning (%s): MaterialDatafile file (%s) not found. Absorption set to 0.\n", "ReadAbsorptionData", AbsorptionFile); return 1; } } // Function used to get the dimensions of the voxels int GetVoxelDimensions (t_Table* TableOfMap, int* NumberOfVoxelsx, int* NumberOfVoxelsy, int* NumberOfVoxelsz) { // Declaration int xValue = 1; int yValue = 1; int zValue = 1; // Get dimensions of voxels in from table header char** parsing; parsing = Table_ParseHeader (TableOfMap->header, "xdim", "XDIM", "XDim", "ydim", "YDIM", "YDim", "zdim", "ZDIM", "ZDim", "Dimensions:", NULL); if (parsing[0]) xValue = strtol (parsing[0], NULL, 0); if (parsing[1]) xValue = strtol (parsing[1], NULL, 0); if (parsing[2]) xValue = strtol (parsing[2], NULL, 0); if (parsing[3]) yValue = strtol (parsing[3], NULL, 0); if (parsing[4]) yValue = strtol (parsing[4], NULL, 0); if (parsing[5]) yValue = strtol (parsing[5], NULL, 0); if (parsing[6]) zValue = strtol (parsing[6], NULL, 0); if (parsing[7]) zValue = strtol (parsing[7], NULL, 0); if (parsing[8]) zValue = strtol (parsing[8], NULL, 0); if (parsing[9]) sscanf (parsing[9], "%d x %d x %d", &xValue, &yValue, &zValue); if (xValue * yValue * zValue == TableOfMap->rows) { *NumberOfVoxelsx = xValue; *NumberOfVoxelsy = yValue; *NumberOfVoxelsz = zValue; return 0; } else { return -1; } } %} DECLARE %{ // Info on polycrystal struct hklStruct* PolyInfo; struct hklStruct hklInfo; struct AbsorbtionStruct AbsorbtionInfo; t_Table* TableOfMap; t_Table* TableOfOrientations; // Dimensions of voxels double Voxelxwidth; double Voxelyheight; double Voxelzdepth; int NumberOfVoxelsx; int NumberOfVoxelsy; int NumberOfVoxelsz; %} INITIALIZE %{ // Declarations int i; Coords aOriginal; Coords bOriginal; Coords cOriginal; Coords aCurrent; Coords bCurrent; Coords cCurrent; Rotation CurrentRotationMatrix; // In-plane anisotropy double* MosaicityAB; // Allocation of tables TableOfMap = malloc (sizeof (t_Table)); if ((i = Table_Read (TableOfMap, MapFile, 0)) == -1) { fprintf (stderr, "Error (%s): Unable to read map file \'%s\'. Aborting.\n", NAME_CURRENT_COMP, MapFile); exit (-1); } TableOfOrientations = malloc (sizeof (t_Table)); if ((i = Table_Read (TableOfOrientations, OrientationsFile, 0)) == -1) { fprintf (stderr, "Error (%s): Unable to read orientation file \'%s\'. Aborting.\n", NAME_CURRENT_COMP, OrientationsFile); exit (-1); } PolyInfo = calloc (TableOfOrientations->rows, sizeof (struct hklStruct)); // Compute dimensions of voxels if (GetVoxelDimensions (TableOfMap, &NumberOfVoxelsx, &NumberOfVoxelsy, &NumberOfVoxelsz) == -1) { fprintf (stderr, "Warning (%s): An error occured when assigning the dimensions of the polycrystal. Is the header appropriately formatted?\n", NAME_CURRENT_COMP); } Voxelxwidth = xwidth / NumberOfVoxelsx; Voxelyheight = yheight / NumberOfVoxelsy; Voxelzdepth = zdepth / NumberOfVoxelsz; // Coordinate calculations aOriginal = coords_set (ax, ay, az); bOriginal = coords_set (bx, by, bz); cOriginal = coords_set (cx, cy, cz); // Set up default hklInfo hklInfo.LatticeVectora[0] = ax; hklInfo.LatticeVectora[1] = ay; hklInfo.LatticeVectora[2] = az; hklInfo.LatticeVectorb[0] = bx; hklInfo.LatticeVectorb[1] = by; hklInfo.LatticeVectorb[2] = bz; hklInfo.LatticeVectorc[0] = cx; hklInfo.LatticeVectorc[1] = cy; hklInfo.LatticeVectorc[2] = cz; hklInfo.Deltad_d = DeltadOverd; hklInfo.SigmaAbs = SigmaAbsorbtion; hklInfo.SigmaInc = SigmaIncoherent; // Default format: h, k, l, F, F2 hklInfo.ColumnOrder[0] = 1; hklInfo.ColumnOrder[1] = 2; hklInfo.ColumnOrder[2] = 3; hklInfo.ColumnOrder[3] = 0; hklInfo.ColumnOrder[4] = 7; // Read in structure factors and absorption info ReadReflectionData (ReflectionsDatafile, &hklInfo, Mosaicity, MosaicityA, MosaicityB, MosaicityC, MosaicityAB); ReadAbsorbtionData (MaterialDatafile, &AbsorbtionInfo); // Print properties if (verbose) { fprintf (stderr, "INFO (%s): %d reflections found in %s.\n", NAME_CURRENT_COMP, hklInfo.NumberOfReflections, ReflectionsDatafile); fprintf (stderr, "INFO (%s): %d different orientations found in %s.\n", NAME_CURRENT_COMP, (int)TableOfOrientations->rows, OrientationsFile); fprintf (stderr, "INFO (%s): %d voxels organised in a %d by %d by %d-grid found in %s.\n", NAME_CURRENT_COMP, (int)TableOfMap->rows, NumberOfVoxelsx, NumberOfVoxelsy, NumberOfVoxelsz, MapFile); fprintf (stderr, "INFO (%s): The voxels of the polycrystals have dimensions of %g m by %g m by %g m.\n", NAME_CURRENT_COMP, Voxelxwidth, Voxelyheight, Voxelzdepth); fprintf (stderr, "INFO (%s): Volume of unit cell is %g AA^3.\n", NAME_CURRENT_COMP, hklInfo.VolumeOfUnitCell); } // Loop over all orientations for (i = 0; i < TableOfOrientations->rows; ++i) { PolyInfo[i] = hklInfo; PolyInfo[i].Recip = 0; // Apply rotation to crystal vectors memcpy (*CurrentRotationMatrix, &(TableOfOrientations->data[i * 9]), sizeof (CurrentRotationMatrix[0][0]) * 9); aCurrent = rot_apply (CurrentRotationMatrix, aOriginal); bCurrent = rot_apply (CurrentRotationMatrix, bOriginal); cCurrent = rot_apply (CurrentRotationMatrix, cOriginal); // Set the crystal parameters to the rotated ones coords_get (aCurrent, &(PolyInfo[i].LatticeVectora[0]), &(PolyInfo[i].LatticeVectora[1]), &(PolyInfo[i].LatticeVectora[2])); coords_get (bCurrent, &(PolyInfo[i].LatticeVectorb[0]), &(PolyInfo[i].LatticeVectorb[1]), &(PolyInfo[i].LatticeVectorb[2])); coords_get (cCurrent, &(PolyInfo[i].LatticeVectorc[0]), &(PolyInfo[i].LatticeVectorc[1]), &(PolyInfo[i].LatticeVectorc[2])); // Read reflections for the given orientation ReadReflectionData (ReflectionsDatafile, &(PolyInfo[i]), Mosaicity, MosaicityA, MosaicityB, MosaicityC, MosaicityAB); } %} TRACE %{ // Dummy integers int i; int j; // Sample propagation double l1 = 0.0; double l2 = 0.0; double l; int Intersect = 0; int IntersectVoxel = 1; int BackgroundVoxel; // Structs for hkl-data struct hklDataStruct* CurrenthklData; // List of reflections close to Ewald sphere struct TauDataStruct CurrentTauData[MXPX_REFL_SLIST_SIZE]; int TauCount; // Number of scattering events int NumberOfScatteringEvents; // Wavewectors double ki[3]; double kiLength; double kf[3]; double kfLength; // The vector ki - tau double Rho[3]; double RhoLength; // Properties of crystal double VolumeOfUnitCell; // Vectors spanning the Ewald sphere tangent plane double b1[3]; double b2[3]; // Matrix describing the 2D-Gaussian double N11; double N12; double N22; double DeterminantN; double NInverse11; double NInverse12; double NInverse22; // L = 1/2 * N^-1 (Cholesky decomposition) double L11; double L12; double L22; double DeterminantL; // Center of 2D-Gaussian double GaussCenterx; double GaussCentery; // Offset of 2D-Gaussian double Alpha; // Product of B * D * Origin double Productx; double Producty; // Coherent reflectivity double TotalReflectivity; // Dummy variables used to select the final wavevector from the 2D-Gaussian double z1; double z2; double y1; double y2; // Variables used in search for tau double r; double Sum; double TauMax; // Transmission, absorption and incoherence double AbsorbtionCrosssection; double AbsorbtionScatteringLength; double IncoherentCrosssection; double IncoherentScatteringLength; double CoherentCrosssection; double CoherentScatteringLength; double TotalCrosssection; double TotalScatteringLength; double CrosssectionFactor; double AbsorbtionMu; // Monte Carlo-properties double CalculatedProbabilityOfTransmission; double MCTransmission; double MCInteract; // Indices of current voxel int VoxelIndexx; int VoxelIndexy; int VoxelIndexz; int PreviousVoxelIndexx = -10000; int PreviousVoxelIndexy = -10000; int PreviousVoxelIndexz = -10000; double VoxelCenterx; double VoxelCentery; double VoxelCenterz; int OrientationIndex; // Compute intersection trajectory Intersect = box_intersect (&l1, &l2, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); if (l2 < 0) { Intersect = 0; } if (l1 > 0) { PROP_DL (l1); } NumberOfScatteringEvents = 0; // Begin computating if ray intersects crystal hull while (Intersect) { // Figure out which voxel the ray hit (the small added constant prevents errors on boundaries of the crystal hull VoxelIndexx = floor ((x + xwidth / 2.0) / xwidth * NumberOfVoxelsx + 1e-10); VoxelIndexy = floor ((y + yheight / 2.0) / yheight * NumberOfVoxelsy + 1e-10); VoxelIndexz = floor ((z + zdepth / 2.0) / zdepth * NumberOfVoxelsz + 1e-10); if (VoxelIndexx >= NumberOfVoxelsx || VoxelIndexx < 0 || VoxelIndexy >= NumberOfVoxelsy || VoxelIndexy < 0 || VoxelIndexz >= NumberOfVoxelsz || VoxelIndexz < 0) { Intersect = 0; break; } // Check, if ray is stuck if (PreviousVoxelIndexx == VoxelIndexx && PreviousVoxelIndexy == VoxelIndexy && PreviousVoxelIndexz == VoxelIndexz) { ABSORB; } // Determine coordinates of voxel center VoxelCenterx = (VoxelIndexx - (NumberOfVoxelsx - 1) / 2.0) * Voxelxwidth; VoxelCentery = (VoxelIndexy - (NumberOfVoxelsy - 1) / 2.0) * Voxelyheight; VoxelCenterz = (VoxelIndexz - (NumberOfVoxelsz - 1) / 2.0) * Voxelzdepth; // Get hklInfo struct from table - orientation index is in the 4th column of TableOfMap (-1 since array is zero-indexed) OrientationIndex = Table_Index (*TableOfMap, VoxelIndexx * (NumberOfVoxelsy * NumberOfVoxelsz) + VoxelIndexy * NumberOfVoxelsz + VoxelIndexz, 3) - 1; VolumeOfUnitCell = PolyInfo[OrientationIndex].VolumeOfUnitCell; // Check, if the ray is a void in the crystal if (OrientationIndex != -1) { BackgroundVoxel = 0; // Prepare voxel loop kiLength = sqrt (pow (kx, 2) + pow (ky, 2) + pow (kz, 2)); CrosssectionFactor = pow (2.0 * M_PI, 5.0 / 2.0) / (VolumeOfUnitCell * pow (kiLength, 2)); // Absorption cross-section AbsorbtionMu = Table_Value (AbsorbtionInfo.Table, kiLength * K2E, AbsorbtionInfo.muc); AbsorbtionCrosssection = /*PolyInfo[OrientationIndex].*/ SigmaAbsorbtion * AbsorbtionMu; AbsorbtionScatteringLength = AbsorbtionCrosssection / VolumeOfUnitCell; // Incoherent scattering IncoherentCrosssection = /*PolyInfo[OrientationIndex].*/ SigmaIncoherent; IncoherentScatteringLength = IncoherentCrosssection / VolumeOfUnitCell; // Get info from PolyInfo-struct CurrenthklData = PolyInfo[OrientationIndex].hklData; } else { BackgroundVoxel = 1; } // Begin loop over multiple scattering events IntersectVoxel = box_intersect (&l1, &l2, x - VoxelCenterx, y - VoxelCentery, z - VoxelCenterz, kx, ky, kz, Voxelxwidth, Voxelyheight, Voxelzdepth); while (IntersectVoxel) { // Remember current voxel PreviousVoxelIndexx = VoxelIndexx; PreviousVoxelIndexy = VoxelIndexy; PreviousVoxelIndexz = VoxelIndexz; // Check, if the x-ray is leaving the sample unexpectedly if (!IntersectVoxel || l2 < -1e-9 || l1 > 1e-9) { fprintf (stderr, "%s: Warning: Ray number %d has unexpectedly left the crystal.\n l1 = %g l2 = %g x = %g y = %g z = %g kx = %g ky = %g kz = %g.\n", NAME_CURRENT_COMP, (int)mcget_run_num (), l1, l2, x, y, z, kx, ky, kz); break; } if (BackgroundVoxel) { PROP_DL (l2); break; } // Copy incoming wave vector ki ki[0] = kx; ki[1] = ky; ki[2] = kz; // Calculate Intersection of Ewald sphere with reciprocal lattice points CoherentCrosssection = 0.0; TotalReflectivity = 0.0; TauMax = 2.0 * kiLength / (1.0 - 5.0 * DeltadOverd); j = 0; for (i = 0; i < PolyInfo[OrientationIndex].NumberOfReflections; ++i) { // End loop if the largest relevant length of tau has been reached if (TauMax < CurrenthklData[i].TauLength) { break; } // Check, if this reciprocal lattice point is close enough to the Ewald sphere to make scattering possible Rho[0] = ki[0] - CurrenthklData[i].Tau[0]; Rho[1] = ki[1] - CurrenthklData[i].Tau[1]; Rho[2] = ki[2] - CurrenthklData[i].Tau[2]; RhoLength = sqrt (pow (Rho[0], 2) + pow (Rho[1], 2) + pow (Rho[2], 2)); // Check if scattering is possible (cutoff of Gaussian tails) if (fabs (RhoLength - kiLength) < CurrenthklData[i].GaussianCutoff) { // Store reflection CurrentTauData[j].Index = i; // Get ki vector in local coordinates CurrentTauData[j].ki[0] = ki[0] * CurrenthklData[i].u1[0] + ki[1] * CurrenthklData[i].u1[1] + ki[2] * CurrenthklData[i].u1[2]; CurrentTauData[j].ki[1] = ki[0] * CurrenthklData[i].u2[0] + ki[1] * CurrenthklData[i].u2[1] + ki[2] * CurrenthklData[i].u2[2]; CurrentTauData[j].ki[2] = ki[0] * CurrenthklData[i].u3[0] + ki[1] * CurrenthklData[i].u3[1] + ki[2] * CurrenthklData[i].u3[2]; CurrentTauData[j].Rho[0] = CurrentTauData[j].ki[0] - CurrenthklData[i].TauLength; CurrentTauData[j].Rho[1] = CurrentTauData[j].ki[1]; CurrentTauData[j].Rho[2] = CurrentTauData[j].ki[2]; CurrentTauData[j].RhoLength = RhoLength; // Compute the tangent plane of the Ewald sphere CurrentTauData[j].NormalVector[0] = CurrentTauData[j].Rho[0] / CurrentTauData[j].RhoLength; CurrentTauData[j].NormalVector[1] = CurrentTauData[j].Rho[1] / CurrentTauData[j].RhoLength; CurrentTauData[j].NormalVector[2] = CurrentTauData[j].Rho[2] / CurrentTauData[j].RhoLength; CurrentTauData[j].Origin[0] = (kiLength - CurrentTauData[j].RhoLength) * CurrentTauData[j].NormalVector[0]; CurrentTauData[j].Origin[1] = (kiLength - CurrentTauData[j].RhoLength) * CurrentTauData[j].NormalVector[1]; CurrentTauData[j].Origin[2] = (kiLength - CurrentTauData[j].RhoLength) * CurrentTauData[j].NormalVector[2]; // Compute unit vectors b1 and b2 that span the tangent plane normal_vec (&(b1[0]), &(b1[1]), &(b1[2]), CurrentTauData[j].NormalVector[0], CurrentTauData[j].NormalVector[1], CurrentTauData[j].NormalVector[2]); vec_prod (b2[0], b2[1], b2[2], CurrentTauData[j].NormalVector[0], CurrentTauData[j].NormalVector[1], CurrentTauData[j].NormalVector[2], b1[0], b1[1], b1[2]); CurrentTauData[j].b1[0] = b1[0]; CurrentTauData[j].b1[1] = b1[1]; CurrentTauData[j].b1[2] = b1[2]; CurrentTauData[j].b2[0] = b2[0]; CurrentTauData[j].b2[1] = b2[1]; CurrentTauData[j].b2[2] = b2[2]; // Compute the 2D projection of the 3D Gauss of the reflection - the symmetric 2x2 matrix N describing the 2D gauss N11 = CurrenthklData[i].m[0] * pow (b1[0], 2) + CurrenthklData[i].m[1] * pow (b1[1], 2) + CurrenthklData[i].m[2] * pow (b1[2], 2); N12 = CurrenthklData[i].m[0] * b1[0] * b2[0] + CurrenthklData[i].m[1] * b1[1] * b2[1] + CurrenthklData[i].m[2] * b1[2] * b2[2]; N22 = CurrenthklData[i].m[0] * pow (b2[0], 2) + CurrenthklData[i].m[1] * pow (b2[1], 2) + CurrenthklData[i].m[2] * pow (b2[2], 2); // The (symmetric) inverse matrix of N DeterminantN = N11 * N22 - pow (N12, 2); NInverse11 = N22 / DeterminantN; NInverse12 = -N12 / DeterminantN; NInverse22 = N11 / DeterminantN; // The Cholesky decomposition of 1/2 * NInverse (lower triangular L) L11 = sqrt (NInverse11 / 2.0); L12 = NInverse12 / (2.0 * L11); L22 = sqrt (NInverse22 / 2.0 - L12 * L12); CurrentTauData[j].L11 = L11; CurrentTauData[j].L12 = L12; CurrentTauData[j].L22 = L22; DeterminantL = L11 * L22; CurrentTauData[j].DeterminantL = DeterminantL; // Compute product Productx = b1[0] * CurrenthklData[i].m[0] * CurrentTauData[j].Origin[0] + b1[1] * CurrenthklData[i].m[1] * CurrentTauData[j].Origin[1] + b1[2] * CurrenthklData[i].m[2] * CurrentTauData[j].Origin[2]; Producty = b2[0] * CurrenthklData[i].m[0] * CurrentTauData[j].Origin[0] + b2[1] * CurrenthklData[i].m[1] * CurrentTauData[j].Origin[1] + b2[2] * CurrenthklData[i].m[2] * CurrentTauData[j].Origin[2]; // Center of 2D Gauss in plane coordinates GaussCenterx = -Productx * NInverse11 - Producty * NInverse12; GaussCentery = -Productx * NInverse12 - Producty * NInverse22; CurrentTauData[j].GaussCenterx = GaussCenterx; CurrentTauData[j].GaussCentery = GaussCentery; // Factor Alpha for the distance of the 2D Gauss from the origin Alpha = CurrenthklData[i].m[0] * pow (CurrentTauData[j].Origin[0], 2) + CurrenthklData[i].m[1] * pow (CurrentTauData[j].Origin[1], 2) + CurrenthklData[i].m[2] * pow (CurrentTauData[j].Origin[2], 2) - (pow (GaussCenterx, 2) * N11 + pow (GaussCentery, 2) * N22 + 2.0 * GaussCenterx * GaussCentery * N12); CurrentTauData[j].Reflectivity = CrosssectionFactor * DeterminantL * exp (-Alpha) / CurrenthklData[i].TotalSigma; TotalReflectivity += CurrentTauData[j].Reflectivity; CurrentTauData[j].Crosssection = CurrentTauData[j].Reflectivity * CurrenthklData[i].F2; CoherentCrosssection += CurrentTauData[j].Crosssection; ++j; } } TauCount = j; if (TauCount == 0) { IntersectVoxel = 0; PROP_DL (l2); break; } // Probabilities of the different possible interactions TotalCrosssection = AbsorbtionCrosssection + IncoherentCrosssection + CoherentCrosssection; CoherentScatteringLength = CoherentCrosssection / VolumeOfUnitCell; TotalScatteringLength = TotalCrosssection / VolumeOfUnitCell; // Calculate and account for transmission CalculatedProbabilityOfTransmission = exp (-TotalScatteringLength * l2); if (!NumberOfScatteringEvents && ProbabilityOfTransmission >= 0.0 && ProbabilityOfTransmission <= 1.0) { MCTransmission = ProbabilityOfTransmission; } else { MCTransmission = CalculatedProbabilityOfTransmission; } MCInteract = 1.0 - MCTransmission; if (MCTransmission > 0.0 && (MCTransmission > 1.0 || rand01 () < MCTransmission)) { p *= CalculatedProbabilityOfTransmission / MCTransmission; SCATTER; IntersectVoxel = 0; PROP_DL (l2); break; } if (TotalScatteringLength < 0) { ABSORB; } if (MCInteract < 0) { IntersectVoxel = 0; PROP_DL (l2); break; } // Scatter the ray if (!NumberOfScatteringEvents) { p *= fabs (1.0 - CalculatedProbabilityOfTransmission) / MCInteract; } // Calculate, where the ray is scattered - use linear sampling for weak scatterers if (TotalScatteringLength * l2 < 1e-6) { l = rand0max (l2); } else { l = -log (1.0 - rand0max (1.0 - exp (-TotalScatteringLength * l2))) / TotalScatteringLength; } // Propagate the ray to this point PROP_DL (l); ++NumberOfScatteringEvents; // Account for absorption p *= (CoherentScatteringLength + IncoherentScatteringLength) / TotalScatteringLength; // Randomly select between coherent and incoherent scattering if (rand0max (CoherentScatteringLength + IncoherentScatteringLength) < IncoherentScatteringLength) { randvec_target_circle (&ki[0], &ki[1], &ki[2], NULL, kx, ky, kz, 0); kx = ki[0]; ky = ki[1]; kz = ki[2]; } else { if (TotalReflectivity < 0) { ABSORB; } // Decide, which reflection will interact with the ray r = rand0max (TotalReflectivity); Sum = 0; for (j = 0; j < TauCount; ++j) { Sum += CurrentTauData[j].Reflectivity; if (Sum > r) { break; } } // If no suitable reflection is found - use the last one in the list (this should not happen in practice) if (j >= TauCount) { fprintf (stderr, "%s: Error: Illegal tau search: r = %g, Sum = %g.\n", NAME_CURRENT_COMP, r, Sum); j = TauCount - 1; } i = CurrentTauData[j].Index; // Pick scattered wavevector kf from 2D Gauss distribution z1 = randnorm (); z2 = randnorm (); y1 = CurrentTauData[j].L11 * z1 + CurrentTauData[j].GaussCenterx; y2 = CurrentTauData[j].L12 * z1 + CurrentTauData[j].L22 * z2 + CurrentTauData[j].GaussCentery; kf[0] = CurrentTauData[j].Rho[0] + CurrentTauData[j].Origin[0] + CurrentTauData[j].b1[0] * y1 + CurrentTauData[j].b2[0] * y2; kf[1] = CurrentTauData[j].Rho[1] + CurrentTauData[j].Origin[1] + CurrentTauData[j].b1[1] * y1 + CurrentTauData[j].b2[1] * y2; kf[2] = CurrentTauData[j].Rho[2] + CurrentTauData[j].Origin[2] + CurrentTauData[j].b1[2] * y1 + CurrentTauData[j].b2[2] * y2; // Normalize kf to length of ki to account for plane approximation of the Ewald sphere kfLength = sqrt (kf[0] * kf[0] + kf[1] * kf[1] + kf[2] * kf[2]); kf[0] *= kiLength / kfLength; kf[1] *= kiLength / kfLength; kf[2] *= kiLength / kfLength; // Adjust weight p *= CurrentTauData[j].Crosssection * TotalReflectivity / (CoherentCrosssection * CurrentTauData[j].Reflectivity); // Adjust direction of ray kx = CurrenthklData[i].u1[0] * kf[0] + CurrenthklData[i].u2[0] * kf[1] + CurrenthklData[i].u3[0] * kf[2]; ky = CurrenthklData[i].u1[1] * kf[0] + CurrenthklData[i].u2[1] * kf[1] + CurrenthklData[i].u3[1] * kf[2]; kz = CurrenthklData[i].u1[2] * kf[0] + CurrenthklData[i].u2[2] * kf[1] + CurrenthklData[i].u3[2] * kf[2]; } SCATTER; /* Calculate new trajectory (the 1.0001 has been added to make sure that the ray leaves the current voxel*/ IntersectVoxel = box_intersect (&l1, &l2, x - VoxelCenterx, y - VoxelCentery, z - VoxelCenterz, kx, ky, kz, 1.0001 * Voxelxwidth, 1.0001 * Voxelyheight, 1.0001 * Voxelzdepth); /* End loop if maximum numer of scattering events has been reached */ if (MaxNumberOfReflections && NumberOfScatteringEvents >= MaxNumberOfReflections) { box_intersect (&l1, &l2, x - VoxelCenterx, y - VoxelCentery, z - VoxelCenterz, kx, ky, kz, xwidth, yheight, zdepth); IntersectVoxel = 0; Intersect = 0; PROP_DL (l2); break; } } } // Free appropriate pointers // free(CurrenthklData); // free(CurrentTauData); %} FINALLY %{ // Free remaining pointers // free(PolyInfo); Table_Free (TableOfMap); Table_Free (TableOfOrientations); %} MCDISPLAY %{ int i; int j; int k; double xmin; double xmax; double ymin; double ymax; double zmin; double zmax; for (i = 0; i < NumberOfVoxelsx; ++i) { xmin = -0.5 * xwidth + Voxelxwidth * i; xmax = -0.5 * xwidth + Voxelxwidth * (i + 1); for (j = 0; j < NumberOfVoxelsy; ++j) { ymin = -0.5 * yheight + Voxelyheight * j; ymax = -0.5 * yheight + Voxelyheight * (j + 1); for (k = 0; k < NumberOfVoxelsz; ++k) { zmin = -0.5 * zdepth + Voxelzdepth * k; zmax = -0.5 * zdepth + Voxelzdepth * (k + 1); line (xmin, ymin, zmin, xmax, ymin, zmin); line (xmin, ymin, zmin, xmin, ymax, zmin); line (xmax, ymin, zmin, xmax, ymax, zmin); line (xmin, ymax, zmin, xmax, ymax, zmin); line (xmin, ymin, zmin, xmin, ymin, zmax); line (xmax, ymin, zmin, xmax, ymin, zmax); line (xmin, ymax, zmin, xmin, ymax, zmax); line (xmax, ymax, zmin, xmax, ymax, zmax); line (xmin, ymin, zmax, xmax, ymin, zmax); line (xmin, ymin, zmax, xmin, ymax, zmax); line (xmax, ymin, zmax, xmax, ymax, zmax); line (xmin, ymax, zmax, xmax, ymax, zmax); } } } %} END