/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*
* Component: Magnetic_single_crystal
*
* %I
* Written by: Erik B Knudsen and Linda Udby
* Date: Jan 2020
* Origin: DTU Physics
* Release: McStas 2.6
*
* Mosaic magnetic single crystal with multiple scattering vectors.
*
* %D
* WARNING: This is an experimental component - no experimental validation has yet been done
*
* Single magnetic crystal with mosaic. Delta-D/D option for finite-size effects.
* Multiple scattering and secondary extinction included.
* The mosaic may EITHER be specified isotropic by setting the mosaic input
* parameter, OR anisotropic by setting the mosaic_h, mosaic_v, and mosaic_n
* parameters.
*
* The scattering is computed solely in an spin up-down configuration. That is the
* scattering is considered in relation to the externally defined vector (mx,my,mz), where it
* can be either SF or NSF.
* Simplifications and comments :
*  Lande splitting factor is assumed to be g=2
*  Magnetic form factors are set =1
*
* <b>Sample shape:</b>
* Sample shape may be a cylinder, a sphere, a box or any other shape
*   box/plate:       xwidth x yheight x zdepth
*   cylinder:        radius x yheight
*   sphere:          radius (yheight=0)
*   any shape:       geometry=OFF file
*
*   The complex geometry option handles any closed non-convex polyhedra.
*   It computes the intersection points of the neutron ray with the object  
*   transparently, so that it can be used like a regular sample object.
*   It supports the OFF and NOFF file format but not COFF (colored faces). 
*   Such files may be generated from XYZ data using qhull/powercrust, and 
*   viewed with geomview
*   The default size of the object depends of the OFF file data, but its 
*   bounding box may be resized using xwidth,yheight and zdepth.
*
* Also, always use a non-zero value of delta_d_d.
*
*
* %P
* INPUT PARAMETERS:
* radius:        [m] Outer radius of sample in (x,z) plane
* xwidth:        [m] Width of crystal
* yheight:       [m] Height of crystal
* zdepth:        [m] Depth of crystal (no extinction simulated)
* geometry:    [str] Name of an Object File Format (OFF) file for complex geometry. The OFF file may be generated from XYZ coordinates using qhull/powercrust
* delta_d_d:     [1] Lattice spacing variance, gaussian RMS
* mosaic:   [arcmin] Crystal mosaic (isotropic), gaussian RMS
* mosaic_h: [arcmin] Horizontal (rotation around Y) mosaic (anisotropic), gaussian RMS
* mosaic_v: [arcmin] Vertical (rotation around Z) mosaic (anisotropic), gaussian RMS
* mosaic_n: [arcmin] Out-of-plane (Rotation around X) mosaic (anisotropic), gaussian RMS
* recip_cell:    [1] Choice of direct/reciprocal (0/1) unit cell definition
* ax:  [AA or AA^-1]
* ay:  [AA or AA^-1] Coordinates of first (direct/recip) unit cell vector
* az:  [AA or AA^-1]
* bx:  [AA or AA^-1]
* by:  [AA or AA^-1] Coordinates of second (direct/recip) unit cell vector
* bz:  [AA or AA^-1]
* cx:  [AA or AA^-1]
* cy:  [AA or AA^-1] Coordinates of third (direct/recip) unit cell vector
* cz:  [AA or AA^-1]        
* mx:            [1]
* my:            [1] Coordinates of vector defining the SF/NSF direction
* mz:            [1]
* na:            [1]
* nb:            [1] Unit cell multipliers. The specified unit cell vectors are scaled by these factors. Note that the mulitpliers are applied directly to the raw input data. I.e. if recip. cell vectors  are given, multipliers should be <1 (= 1/n). F.i. used to specify a magnetic unit cell which is larger than the chemical unit cell. 
* nc:            [1]
* atom_sites:  [str] File name containing the atoms present in the unit cell. Use empty ("") or NULL for sigma_inc scattering only
* order:         [1] limit multiple scattering up to given order (0: all, 1: first, 2: second, ...)
* q_min:     [AA^-1] lower boundary of momentum transfer range to generate hkls in
* q_max:     [AA^-1] upper boundary of momentum transfer range to generate hkls in
*
* Optional input parameters:
*
* p_transmit:    [1] Monte Carlo probability for neutrons to be transmitted without any scattering. Used to improve statistics from weak reflections
* sigma_abs: [barns] absorption cross-section per unit cell at 2200 m/s
* sigma_inc: [barns] incoherent scattering cross-section per unit cell
*
* CALCULATED PARAMETERS:
*
* hkl_info: Internal [structure]
*
* Geomview and Object File Format (OFF) <http|://www.geomview.org>
* %E
****************************************************************************/

DEFINE COMPONENT Single_magnetic_crystal

SETTING PARAMETERS(string atom_sites=0, string geometry=NULL,
            xwidth=0, yheight=0, zdepth=0, radius=0, delta_d_d=1e-4,
            mosaic = -1, mosaic_h = -1, mosaic_v = -1, mosaic_n = -1,
            recip_cell=0, q_min=0,q_max=-1,mx=0,my=1,mz=0,
            na=1, nb=1, nc=1, ax = 0, ay = 0, az = 0, bx = 0, by = 0, bz = 0, cx = 0, cy = 0, cz = 0,
            p_transmit = -1, sigma_abs = 0, sigma_inc = 0,order=0)
NOACC

SHARE
%{
  /* used for reading data table from file */
  %include "read_table-lib"
  %include "interoff-lib"
  %include "mccode-complex-lib"
  //  %include "columnfile"
  /* Declare structures and functions only once in each instrument. */
  #ifndef SINGLE_MAGNETIC_CRYSTAL_DECL
  #define SINGLE_MAGNETIC_CRYSTAL_DECL

  struct hkl_data {
    int h, k, l;                /* Indices for this reflection */
    double F2;                  /* Value of structure factor */
    cdouble f[4];               /* Structure factors (scattering amplitudes for different spin flips spin up->up, down->down, up->down, down-up */
    double tau_x, tau_y, tau_z; /* Coordinates in reciprocal space */
    double tau;                 /* Length of (tau_x, tau_y, tau_z) */
    double u1x, u1y, u1z;       /* First axis of local coordinate system */
    double u2x, u2y, u2z;       /* Second axis of local coordinate system */
    double u3x, u3y, u3z;       /* Third axis of local coordinate system */
    double sig1, sig2, sig3;    /* RMSs of Gauss axis */
    double sig123;              /* The product sig1*sig2*sig3 */
    double m1, m2, m3;          /* Diagonal matrix representation of Gauss */
    double cutoff;              /* Cutoff value for Gaussian tails */
  };

  struct hkl_info_struct {
    struct hkl_data* list;     /* Reflection array */
    int count;                 /* Number of reflections */
    struct tau_data* tau_list; /* Reflections close to Ewald Sphere */
    double m_delta_d_d;        /* Delta-d/d FWHM */
    double m_ax, m_ay, m_az;   /* First unit cell axis (direct space, AA) */
    double m_bx, m_by, m_bz;   /* Second unit cell axis */
    double m_cx, m_cy, m_cz;   /* Third unit cell axis */
    double asx, asy, asz;      /* First reciprocal lattice axis (1/AA) */
    double bsx, bsy, bsz;      /* Second reciprocal lattice axis */
    double csx, csy, csz;      /* Third reciprocal lattice axis */
    double m_a, m_b, m_c;      /* lattice parameter lengths */
    double m_aa, m_bb, m_cc;   /* lattice angles */
    double sigma_a, sigma_i;   /* abs and inc X sect */
    double rho;                /* density */
    double at_weight;          /* atomic weight */
    double at_nb;              /* nb of atoms in a cell */
    double V0;                 /* Unit cell volume (AA**3) */
    int column_order[5];       /* column signification [h,k,l,F,F2] */
    int recip;                 /* Flag to indicate if recip or direct cell axes given */
    int shape;                 /* 0:cylinder, 1:box, 2:sphere 3:any shape*/
    int flag_warning;          /* number of warnings */
    double tau_max;
    double tau_min;
    double refx, refy, refz; /*chosen polarisation reference direction*/
  };

  struct tau_data {
    int index; /* Index into reflection table */
    double refl;
    double xsect;
    double sigma_1, sigma_2;
    /* The following vectors are in local koordinates. */
    double kix, kiy, kiz;       /* Initial wave vector */
    double rho_x, rho_y, rho_z; /* The vector ki - tau */
    double rho;                 /* Length of rho vector */
    double ox, oy, oz;          /* Origin of Ewald sphere tangent plane */
    double nx, ny, nz;          /* Normal vector of Ewald sphere tangent */
    double b1x, b1y, b1z;       /* Spanning vectors of Ewald sphere tangent */
    double b2x, b2y, b2z;
    double l11, l12, l22; /* Cholesky decomposition L of 2D Gauss */
    double det_L;         /* Determinant of L */
    double y0x, y0y;      /* 2D Gauss center in tangent plane */
  };

  int
  read_hkl_data (char* atoms_file, struct hkl_info_struct* info, double SC_mosaic, double SC_mosaic_h, double SC_mosaic_v, double SC_mosaic_n) { /*{{{*/
    struct hkl_data* list = NULL;
    int size = 0;
    t_Table sTable; /* sample data table structure from atoms_file */
    int i = 0;
    double tmp_x, tmp_y, tmp_z;
    char** parsing;
    char flag = 0;

    double nb_atoms = 1;
    t_Table atoms;

    if (!atoms_file || !strlen (atoms_file) || !strcmp (atoms_file, "NULL") || !strcmp (atoms_file, "0")) {
      info->count = 0;
      flag = 1;
    }
    if (!flag) {
      int status;
      if ((status = Table_Read (&atoms, atoms_file, 0)) == -1) {
        fprintf (stderr, "Single_magnetic_crystal: Error reading atom list from file: %s\n", atoms_file);
        return (0);
      }

      printf ("index type x y z b\n");
      for (i = 0; i < atoms.rows; i++) {
        printf ("atom %d", (int)floor (Table_Index (atoms, i, 0)));
        int j;
        for (j = 1; j < atoms.columns; j++) {
          printf (" %g", Table_Index (atoms, i, j));
        }
        printf ("\n");
      }

      /* Compute reciprocal or direct lattice vectors. */
      if (!info->recip) { /*{{{*/
        vec_prod (tmp_x, tmp_y, tmp_z, info->m_bx, info->m_by, info->m_bz, info->m_cx, info->m_cy, info->m_cz);
        info->V0 = fabs (scalar_prod (info->m_ax, info->m_ay, info->m_az, tmp_x, tmp_y, tmp_z));
        printf ("V0=%g\n", info->V0);

        info->asx = 2 * PI / info->V0 * tmp_x;
        info->asy = 2 * PI / info->V0 * tmp_y;
        info->asz = 2 * PI / info->V0 * tmp_z;
        vec_prod (tmp_x, tmp_y, tmp_z, info->m_cx, info->m_cy, info->m_cz, info->m_ax, info->m_ay, info->m_az);
        info->bsx = 2 * PI / info->V0 * tmp_x;
        info->bsy = 2 * PI / info->V0 * tmp_y;
        info->bsz = 2 * PI / info->V0 * tmp_z;
        vec_prod (tmp_x, tmp_y, tmp_z, info->m_ax, info->m_ay, info->m_az, info->m_bx, info->m_by, info->m_bz);
        info->csx = 2 * PI / info->V0 * tmp_x;
        info->csy = 2 * PI / info->V0 * tmp_y;
        info->csz = 2 * PI / info->V0 * tmp_z;
      } else {
        info->asx = info->m_ax;
        info->asy = info->m_ay;
        info->asz = info->m_az;
        info->bsx = info->m_bx;
        info->bsy = info->m_by;
        info->bsz = info->m_bz;
        info->csx = info->m_cx;
        info->csy = info->m_cy;
        info->csz = info->m_cz;

        vec_prod (tmp_x, tmp_y, tmp_z, info->bsx / (2 * PI), info->bsy / (2 * PI), info->bsz / (2 * PI), info->csx / (2 * PI), info->csy / (2 * PI),
                  info->csz / (2 * PI));
        info->V0 = 1 / fabs (scalar_prod (info->asx / (2 * PI), info->asy / (2 * PI), info->asz / (2 * PI), tmp_x, tmp_y, tmp_z));
        printf ("V0=%g\n", info->V0);

        info->asx = 2 * PI / info->V0 * tmp_x;
        info->asy = 2 * PI / info->V0 * tmp_y;
        info->asz = 2 * PI / info->V0 * tmp_z;
        vec_prod (tmp_x, tmp_y, tmp_z, info->m_cx, info->m_cy, info->m_cz, info->m_ax, info->m_ay, info->m_az);
        info->bsx = 2 * PI / info->V0 * tmp_x;
        info->bsy = 2 * PI / info->V0 * tmp_y;
        info->bsz = 2 * PI / info->V0 * tmp_z;
        vec_prod (tmp_x, tmp_y, tmp_z, info->m_ax, info->m_ay, info->m_az, info->m_bx, info->m_by, info->m_bz);
        info->csx = 2 * PI / info->V0 * tmp_x;
        info->csy = 2 * PI / info->V0 * tmp_y;
        info->csz = 2 * PI / info->V0 * tmp_z;
      } /*}}}*/
      /*store lattice vector lengths for later reference*/
      info->m_a = sqrt (scalar_prod (info->m_ax, info->m_ay, info->m_az, info->m_ax, info->m_ay, info->m_az));
      info->m_b = sqrt (scalar_prod (info->m_bx, info->m_by, info->m_bz, info->m_bx, info->m_by, info->m_bz));
      info->m_c = sqrt (scalar_prod (info->m_cx, info->m_cy, info->m_cz, info->m_cx, info->m_cy, info->m_cz));

      /*give the atom list - calculate the chemical structure factors of a series of hkls that are within the q-range*/
      double as, bs, cs;
      as = sqrt (scalar_prod (info->asx, info->asy, info->asz, info->asx, info->asy, info->asz));
      bs = sqrt (scalar_prod (info->bsx, info->bsy, info->bsz, info->bsx, info->bsy, info->bsz));
      cs = sqrt (scalar_prod (info->csx, info->csy, info->csz, info->csx, info->csy, info->csz));

      cdouble f = cplx (0, 0);
      if (info->tau_max == -1)
        info->tau_max = 4 * as;

      double q;
      int h, k, l, m;

      /* allocate hkl_data array initially make room for 2048 reflections. will be expanded (realloc'd) if needed.*/
      size = 2048;
      list = (struct hkl_data*)malloc (size * sizeof (struct hkl_data));
      if (!list) {
        fprintf (stderr, "Single_magnetic_crystal: Error allocating reflection list\n");
        return (-1);
      }
      printf ("q=[%g %g]\n", info->tau_min, info->tau_max);
      printf ("as,bs,cs=(%g %g %g)\n", as, bs, cs);
      i = 0;
      for (h = -floor (info->tau_max / as); h < ceil (info->tau_max / as); h++) {
        for (k = -floor (info->tau_max / bs); k < ceil (info->tau_max / bs); k++) {
          for (l = -floor (info->tau_max / cs); l < ceil (info->tau_max / cs); l++) {
            if (h == 0 && k == 0 && l == 0)
              continue;
            double qx = (h * info->asx + k * info->bsx + l * info->csx);
            double qy = (h * info->asy + k * info->bsy + l * info->csy);
            double qz = (h * info->asz + k * info->bsz + l * info->csz);
            q = sqrt (qx * qx + qy * qy + qz * qz);
            if (q < info->tau_min || q > info->tau_max)
              continue;
            f = cplx (0, 0);
            cdouble f_tau = cplx (0, 0);
            cdouble* f_ = malloc (4 * sizeof (cdouble));
            if (!f_) {
              fprintf (stderr, "Single_magnetic_crystal: Error allocating cdouble f_[4]\n");
              return (-1);
            }
            f_[0] = f_[1] = f_[2] = f_[3] = cplx (0, 0);

            for (m = 0; m < atoms.rows; m++) {
              /*Tabulated values are usually in fm. Divide by 10 to end up with x-sections in barns*/
              double b = Table_Index (atoms, m, 5) / 10;
              f_tau = cexp (cplx (0, 2 * M_PI * (h * Table_Index (atoms, m, 2) + k * Table_Index (atoms, m, 3) + l * Table_Index (atoms, m, 4))));
              // printf("%g b=%g r=(%g %g %g) hkl=(%2d %2d %2d),
              // exp(-i*2pi*(H.r)=(%g%+gj)\n",atoms.data[m][0],b,atoms.data[m][2],atoms.data[m][3],atoms.data[m][4],h,k,l,creal(f_tau),cimag(f_tau));
              if (atoms.columns > 6) {
                /*the atom has a magnetic moment*/
                double S_x, S_y, S_z, L_x, L_y, L_z;
                double i_x = 0, i_y = 0, i_z = 0;
                /*G. Williams Definition of r0 / cm =-gamma*e^2/(m_e c^2)*/
                const double r0 = -0.5391e-12;
                double g, gs, gl;
                double M, beta = 0;
                gs = Table_Index (atoms, m, 6);
                /*S_a,S_b,S_c in the file are given in crystal coordinates - so convert that to the crystal cartesian coordinate system*/
                S_x = Table_Index (atoms, m, 7) * info->m_ax / info->m_a + Table_Index (atoms, m, 8) * info->m_bx / info->m_b
                      + Table_Index (atoms, m, 9) * info->m_cx / info->m_c;
                S_y = Table_Index (atoms, m, 7) * info->m_ay / info->m_a + Table_Index (atoms, m, 8) * info->m_by / info->m_b
                      + Table_Index (atoms, m, 9) * info->m_cy / info->m_c;
                S_z = Table_Index (atoms, m, 7) * info->m_az / info->m_a + Table_Index (atoms, m, 8) * info->m_bz / info->m_b
                      + Table_Index (atoms, m, 9) * info->m_cz / info->m_c;
                gl = Table_Index (atoms, m, 10);
                /*L_a,L_b,L_c in the file are given in crystal coordinates - so convert that to the crystal cartesian coordinate system*/
                L_x = Table_Index (atoms, m, 11) * info->m_ax / info->m_a + Table_Index (atoms, m, 12) * info->m_bx / info->m_b
                      + Table_Index (atoms, m, 13) * info->m_cx / info->m_c;
                L_y = Table_Index (atoms, m, 11) * info->m_ay / info->m_a + Table_Index (atoms, m, 12) * info->m_by / info->m_b
                      + Table_Index (atoms, m, 13) * info->m_cy / info->m_c;
                L_z = Table_Index (atoms, m, 11) * info->m_az / info->m_a + Table_Index (atoms, m, 12) * info->m_bz / info->m_b
                      + Table_Index (atoms, m, 13) * info->m_cz / info->m_c;
                g = gs + gl;

                M = fabs (r0) * g / 2 * 1e12; /*2nd factor to end up with x-sections in barns*/
                if ((S_x != 0 || S_y != 0 || S_z != 0)) {
                  /*S_|_=q_norm x (S x _norm) */
                  double S_orto_x, S_orto_y, S_orto_z;
                  vec_prod (S_orto_x, S_orto_y, S_orto_z, S_x, S_y, S_z, qx, qy, qz);
                  vec_prod (S_orto_x, S_orto_y, S_orto_z, qx, qy, qz, S_x, S_y, S_z);
                  S_orto_x /= (q * q);
                  S_orto_y /= (q * q);
                  S_orto_z /= (q * q);

                  /* Refer to coordinates where polarisation is along z*/
                  double mm1sq = scalar_prod (info->refx, info->refy, info->refz, info->refx, info->refy, info->refz);
                  /*gram-schmidt orthogonalization*/
                  double mm2x = 1, mm2y = 0, mm2z = 0, mm2sq = 1, mm3x = 0, mm3y = 0, mm3z = 1, mm3sq = 1;
                  if (info->refy == 0) {
                    if (info->refx == 0) {
                      mm2x = 1;
                      mm2y = 0;
                      mm2z = 0;
                      mm3x = 0;
                      mm3y = 1;
                      mm3z = 0;
                    } else {
                      mm2x = 0;
                      mm2y = 1;
                      mm2z = 0;
                      mm3x = 0;
                      mm3y = 0;
                      mm3z = 1;
                    }
                  }
                  double tmpx, tmpy, tmpz;
                  tmpx = mm2x - scalar_prod (mm2x, mm2y, mm2z, info->refx, info->refy, info->refz) / mm1sq * info->refx;
                  tmpy = mm2y - scalar_prod (mm2x, mm2y, mm2z, info->refx, info->refy, info->refz) / mm1sq * info->refy;
                  tmpz = mm2z - scalar_prod (mm2x, mm2y, mm2z, info->refx, info->refy, info->refz) / mm1sq * info->refz;
                  mm2x = tmpx;
                  mm2y = tmpy;
                  mm2z = tmpz;
                  mm2sq = scalar_prod (mm2x, mm2y, mm2z, mm2x, mm2y, mm2z);

                  tmpx = mm3x - scalar_prod (mm3x, mm3y, mm3z, info->refx, info->refy, info->refz) / mm1sq * info->refx;
                  tmpy = mm3y - scalar_prod (mm3x, mm3y, mm3z, info->refx, info->refy, info->refz) / mm1sq * info->refy;
                  tmpz = mm3z - scalar_prod (mm3x, mm3y, mm3z, info->refx, info->refy, info->refz) / mm1sq * info->refz;
                  mm3x = tmpx;
                  mm3y = tmpy;
                  mm3z = tmpz;
                  tmpx = mm3x - scalar_prod (mm3x, mm3y, mm3z, mm2x, mm2y, mm2z) / mm2sq * info->refx;
                  tmpy = mm3y - scalar_prod (mm3x, mm3y, mm3z, mm2x, mm2y, mm2z) / mm2sq * info->refy;
                  tmpz = mm3z - scalar_prod (mm3x, mm3y, mm3z, mm2x, mm2y, mm2z) / mm2sq * info->refz;
                  mm3x = tmpx;
                  mm3y = tmpy;
                  mm3z = tmpz;

                  S_x = scalar_prod (S_orto_x, S_orto_y, S_orto_z, mm2x, mm2y, mm2z);
                  S_y = scalar_prod (S_orto_x, S_orto_y, S_orto_z, mm3x, mm3y, mm3z);
                  S_z = scalar_prod (S_orto_x, S_orto_y, S_orto_z, info->refx, info->refy, info->refz);
                }
                f_[0] = cadd (f_[0], rmul (b - M * S_z + (beta / 2) * i_z, f_tau));
                f_[1] = cadd (f_[1], rmul (b + M * S_z - (beta / 2) * i_z, f_tau));
                f_[2] = cadd (f_[2], cmul (cadd (rmul (-M, cplx (S_x, S_y)), rmul (beta / 2, cplx (i_x, i_y))), f_tau));
                f_[3] = cadd (f_[3], cmul (cadd (rmul (-M, cplx (S_x, -S_y)), rmul (beta / 2, cplx (i_x, -i_y))), f_tau));
              } else {
                /*scattering is non-magnetic*/
                f_[0] = cadd (f_[0], rmul (b, f_tau));
                f_[1] = cadd (f_[1], rmul (b, f_tau));
                // f_[2]=f_[2]; // +0;
                // f_[3]=f_[3]; // +0;
              }
            }
            if (cabs (f_[0]) > FLT_EPSILON || cabs (f_[1]) > FLT_EPSILON || cabs (f_[2]) > FLT_EPSILON || cabs (f_[3]) > FLT_EPSILON) {
              list[i].h = h;
              list[i].k = k;
              list[i].l = l;
              for (m = 0; m < 4; m++) {
                list[i].f[m] = f_[m];
              }
              list[i].F2 = cabs (f_tau) * cabs (f_tau);
              if (++i == size) {
                size = 2 * size;
                list = (struct hkl_data*)realloc (list, size * sizeof (struct hkl_data));
                if (!list) {
                  fprintf (stderr, "Single_crystal: Error re-allocating reflection list\n");
                  return (-1);
                }
              }
              printf ("(hkl)=(%2d %2d %2d) ", h, k, l);
              printf (" |f++|^2=%g |f--|^2=%g |f+-|^2=%g |f-+|^2=%g", cabs (f_[0]) * cabs (f_[0]), cabs (f_[1]) * cabs (f_[1]), cabs (f_[2]) * cabs (f_[2]),
                      cabs (f_[3]) * cabs (f_[3]));
              printf (" f++=%g%+gj, f--=%g%+gj, f+-=%g%+gj, f-+=%g%+gj\n", creal (f_[0]), cimag (f_[0]), creal (f_[1]), cimag (f_[1]), creal (f_[2]),
                      cimag (f_[2]), creal (f_[3]), cimag (f_[3]));
            }
          }
        }
      }
      /*set size to actual number of reflections*/
      size = i;
    }

    if (flag)
      return (-1);

    /*re-loop over reflections and evaluate mosaicity coefficients etc.*/
    for (i = 0; i < size; i++) {
      double b1[3], b2[3];
      /* Precompute some values */
      list[i].tau_x = list[i].h * info->asx + list[i].k * info->bsx + list[i].l * info->csx;
      list[i].tau_y = list[i].h * info->asy + list[i].k * info->bsy + list[i].l * info->csy;
      list[i].tau_z = list[i].h * info->asz + list[i].k * info->bsz + list[i].l * info->csz;
      list[i].tau = sqrt (list[i].tau_x * list[i].tau_x + list[i].tau_y * list[i].tau_y + list[i].tau_z * list[i].tau_z);
      list[i].u1x = list[i].tau_x / list[i].tau;
      list[i].u1y = list[i].tau_y / list[i].tau;
      list[i].u1z = list[i].tau_z / list[i].tau;
      list[i].sig1 = FWHM2RMS * info->m_delta_d_d * list[i].tau;
      /* Find two arbitrary axes perpendicular to tau and each other. */
      normal_vec (&b1[0], &b1[1], &b1[2], list[i].u1x, list[i].u1y, list[i].u1z);
      vec_prod (b2[0], b2[1], b2[2], list[i].u1x, list[i].u1y, list[i].u1z, b1[0], b1[1], b1[2]);
      /* Find the two mosaic axes perpendicular to tau. */
      if (SC_mosaic > 0) {
        /* Use isotropic mosaic. */
        list[i].u2x = b1[0];
        list[i].u2y = b1[1];
        list[i].u2z = b1[2];
        list[i].sig2 = FWHM2RMS * list[i].tau * MIN2RAD * SC_mosaic;
        list[i].u3x = b2[0];
        list[i].u3y = b2[1];
        list[i].u3z = b2[2];
        list[i].sig3 = FWHM2RMS * list[i].tau * MIN2RAD * SC_mosaic;
      } else {
        /* Use anisotropic mosaic. */
        /*This is not implemeted fully yet (see todo below)- therefore exit with a warning*/
        fprintf (stderr, "Single_magnetic_crystal: Anisotropic mosaic not implemented yet - aborting.\n");
        return (0);
      }
      list[i].sig123 = list[i].sig1 * list[i].sig2 * list[i].sig3;
      list[i].m1 = 1 / (2 * list[i].sig1 * list[i].sig1);
      list[i].m2 = 1 / (2 * list[i].sig2 * list[i].sig2);
      list[i].m3 = 1 / (2 * list[i].sig3 * list[i].sig3);
      /* Set Gauss cutoff to 5 times the maximal sigma. */
      if (list[i].sig1 > list[i].sig2)
        if (list[i].sig1 > list[i].sig3)
          list[i].cutoff = 5 * list[i].sig1;
        else
          list[i].cutoff = 5 * list[i].sig3;
      else if (list[i].sig2 > list[i].sig3)
        list[i].cutoff = 5 * list[i].sig2;
      else
        list[i].cutoff = 5 * list[i].sig3;
    }
    info->list = list;
    info->count = i;
    info->tau_list = malloc (i * sizeof (*info->tau_list));
    if (!info->tau_list) {
      fprintf (stderr, "Single_crystal: Error: Out of memory!\n");
      return (0);
    }
    return (info->count = i);
  } /*}}}*/
  #endif /* !SINGLE_MAGNETIC_CRYSTAL_DECL */
%}

DECLARE
%{
  struct hkl_info_struct hkl_info;
  off_struct offdata;
%}

INITIALIZE
%{
  double as, bs, cs;

  /* transfer input parameters */
  hkl_info.m_delta_d_d = delta_d_d;
  hkl_info.m_ax = na * ax;
  hkl_info.m_ay = na * ay;
  hkl_info.m_az = na * az;
  hkl_info.m_bx = nb * bx;
  hkl_info.m_by = nb * by;
  hkl_info.m_bz = nb * bz;
  hkl_info.m_cx = nc * cx;
  hkl_info.m_cy = nc * cy;
  hkl_info.m_cz = nc * cz;
  hkl_info.sigma_a = sigma_abs;
  hkl_info.sigma_i = sigma_inc;
  hkl_info.recip = recip_cell;
  hkl_info.tau_min = q_min;
  hkl_info.tau_max = q_max;
  if (mx != 0 || my != 0 || mz != 0) {
    hkl_info.refx = mx;
    hkl_info.refy = my;
    hkl_info.refz = mz;
    NORM (hkl_info.refx, hkl_info.refy, hkl_info.refz);
  } else {
    hkl_info.refx = hkl_info.refz = 0;
    hkl_info.refy = 1;
  }
  /* Read in structure factors, and do some pre-calculations. */
  if (!read_hkl_data (atom_sites, &hkl_info, mosaic, mosaic_h, mosaic_v, mosaic_n))
    exit (-1);
  if (hkl_info.count)
    printf ("Single_crystal: %s: Read %d reflections from file '%s'\n", NAME_CURRENT_COMP, hkl_info.count, atom_sites);
  else
    printf ("Single_crystal: %s: Using incoherent elastic scattering only.\n", NAME_CURRENT_COMP, hkl_info.sigma_i);

  hkl_info.shape = -1; /* -1:no shape, 0:cyl, 1:box, 2:sphere, 3:any-shape  */
  if (geometry && strlen (geometry) && strcmp (geometry, "NULL") && strcmp (geometry, "0")) {
    if (off_init (geometry, xwidth, yheight, zdepth, 0, &offdata)) {
      hkl_info.shape = 3;
    }
  } else if (xwidth && yheight && zdepth)
    hkl_info.shape = 1; /* box */
  else if (radius > 0 && yheight)
    hkl_info.shape = 0; /* cylinder */
  else if (radius > 0 && !yheight)
    hkl_info.shape = 2; /* sphere */

  if (hkl_info.shape < 0)
    exit (fprintf (stderr,
                   "Single_magnetic_crystal: %s: sample has invalid dimensions.\n"
                   "ERROR           Please check parameter values (xwidth, yheight, zdepth, radius).\n",
                   NAME_CURRENT_COMP));

  printf ("Single_magnetic_crystal: %s: Vc=%g [Angs] sigma_abs=%g [barn] sigma_inc=%g [barn] reflections=%s\n", NAME_CURRENT_COMP, hkl_info.V0, hkl_info.sigma_a,
          hkl_info.sigma_i, atom_sites&& strlen (atom_sites) ? atom_sites : "NULL");

  printf ("WARNING; Single_magnetic_crystal has not yet been experimentally validated. Please use caution when intepreting results.\n");
%}

TRACE
%{
  double t1, t2 = 0;          /* Entry and exit times in sample */
  struct hkl_data* L;         /* Structure factor list */
  int i;                      /* Index into structure factor list */
  struct tau_data* T;         /* List of reflections close to Ewald sphere */
  int j;                      /* Index into reflection list */
  int event_counter;          /* scattering event counter */
  double kix, kiy, kiz, ki;   /* Initial wave vector [1/AA] */
  double kfx, kfy, kfz;       /* Final wave vector */
  double v;                   /* Neutron velocity */
  double tau_max;             /* Max tau allowing reflection at this ki */
  double rho_x, rho_y, rho_z; /* the vector ki - tau */
  double rho;
  double diff;                      /* Deviation from Bragg condition */
  double ox, oy, oz;                /* Origin of Ewald sphere tangent plane */
  double b1x, b1y, b1z;             /* First vector spanning tangent plane */
  double b2x, b2y, b2z;             /* Second vector spanning tangent plane */
  double n11, n12, n22;             /* 2D Gauss description matrix N */
  double det_N;                     /* Determinant of N */
  double inv_n11, inv_n12, inv_n22; /* Inverse of N */
  double l11, l12, l22;             /* Cholesky decomposition L of 1/2*inv(N) */
  double det_L;                     /* Determinant of L */
  double Bt_D_O_x, Bt_D_O_y;        /* Temporaries */
  double y0x, y0y;                  /* Center of 2D Gauss in plane coordinates */
  double alpha;                     /* Offset of 2D Gauss center from 3D center */
  int tau_count;                    /* Number of reflections within cutoff */
  double V0;                        /* Volume of unit cell */
  double l_full;                    /* Neutron path length for transmission */
  double l;                         /* Path length to scattering event */
  double abs_xsect, abs_xlen;       /* sigma_abs cross section and length */
  double inc_xsect, inc_xlen;       /* sigma_inc cross section and length */
  double coh_xsect, coh_xlen;       /* Coherent cross section and length */
  double tot_xsect, tot_xlen;       /* Total cross section and length */
  double z1, z2, y1, y2;            /* Temporaries to choose kf from 2D Gauss */
  double adjust, coh_refl;          /* Temporaries */
  double r, sum;                    /* Temporaries */
  double xsect_factor;              /* Common factor in coherent cross-section */
  double p_trans;                   /* Transmission probability */
  double mc_trans, mc_interact;     /* Transmission, interaction MC choices */
  int intersect = 0;

  /* Intersection neutron trajectory / sample (sample surface) */
  if (hkl_info.shape == 0)
    intersect = cylinder_intersect (&t1, &t2, x, y, z, vx, vy, vz, radius, yheight);
  else if (hkl_info.shape == 1)
    intersect = box_intersect (&t1, &t2, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);
  else if (hkl_info.shape == 2)
    intersect = sphere_intersect (&t1, &t2, x, y, z, vx, vy, vz, radius);
  else if (hkl_info.shape == 3)
    intersect = off_intersect (&t1, &t2, NULL, NULL, x, y, z, vx, vy, vz, 0, 0, 0, offdata);

  if (t2 < 0)
    intersect = 0; /* we passed sample volume already */

  if (intersect) { /* Neutron intersects crystal */
    if (t1 > 0)
      PROP_DT (t1); /* Move to crystal surface if not inside */
    v = sqrt (vx * vx + vy * vy + vz * vz);
    ki = V2K * v;
    event_counter = 0;
    abs_xsect = hkl_info.sigma_a * 2200 / v;
    inc_xsect = hkl_info.sigma_i;
    V0 = hkl_info.V0;
    abs_xlen = 1e2 * abs_xsect / V0;
    inc_xlen = 1e2 * inc_xsect / V0;
    L = hkl_info.list;
    T = hkl_info.tau_list;

    do { /* Loop over multiple scattering events */

      if (hkl_info.shape == 0)
        intersect = cylinder_intersect (&t1, &t2, x, y, z, vx, vy, vz, radius, yheight);
      else if (hkl_info.shape == 1)
        intersect = box_intersect (&t1, &t2, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);
      else if (hkl_info.shape == 2)
        intersect = sphere_intersect (&t1, &t2, x, y, z, vx, vy, vz, radius);
      else if (hkl_info.shape == 3)
        intersect = off_intersect (&t1, &t2, NULL, NULL, x, y, z, vx, vy, vz, 0, 0, 0, offdata);
      if (!intersect || t2 * v < -1e-9 || t1 * v > 1e-9) {
        /* neutron is leaving the sample */
        if (hkl_info.flag_warning < 100)
          fprintf (stderr,
                   "Single_crystal: %s: Warning: neutron has unexpectedly left the crystal!\n"
                   "                t1=%g t2=%g x=%g y=%g z=%g vx=%g vy=%g vz=%g\n",
                   NAME_CURRENT_COMP, t1, t2, x, y, z, vx, vy, vz);
        hkl_info.flag_warning++;
        break;
      }

      l_full = t2 * v;

      /* (1). Compute incoming wave vector ki */
      kix = V2K * vx;
      kiy = V2K * vy;
      kiz = V2K * vz;

      /* (2). Intersection of Ewald sphere with reciprocal lattice points */

      /* Max possible tau with 5*sigma delta-d/d cutoff. */
      tau_max = 2 * ki / (1 - 5 * hkl_info.m_delta_d_d);

      coh_xsect = 0;
      coh_refl = 0;
      xsect_factor = pow (2 * PI, 5.0 / 2.0) / (V0 * ki * ki);
      for (i = j = 0; i < hkl_info.count; i++) {
        /* Assuming reflections are sorted, stop search when max tau exceeded. */
        if (L[i].tau > tau_max)
          break;
        /* Check if this reciprocal lattice point is close enough to the
           Ewald sphere to make scattering possible. */
        rho_x = kix - L[i].tau_x;
        rho_y = kiy - L[i].tau_y;
        rho_z = kiz - L[i].tau_z;
        rho = sqrt (rho_x * rho_x + rho_y * rho_y + rho_z * rho_z);
        diff = fabs (rho - ki);

        /* Check if scattering is possible (cutoff of Gaussian tails). */
        if (diff <= L[i].cutoff) {
          /* Store reflection. */
          T[j].index = i;
          /* Get ki vector in local coordinates. */
          T[j].kix = kix * L[i].u1x + kiy * L[i].u1y + kiz * L[i].u1z;
          T[j].kiy = kix * L[i].u2x + kiy * L[i].u2y + kiz * L[i].u2z;
          T[j].kiz = kix * L[i].u3x + kiy * L[i].u3y + kiz * L[i].u3z;
          T[j].rho_x = T[j].kix - L[i].tau;
          T[j].rho_y = T[j].kiy;
          T[j].rho_z = T[j].kiz;
          T[j].rho = rho;
          /* Compute the tangent plane of the Ewald sphere. */
          T[j].nx = T[j].rho_x / T[j].rho;
          T[j].ny = T[j].rho_y / T[j].rho;
          T[j].nz = T[j].rho_z / T[j].rho;
          ox = (ki - T[j].rho) * T[j].nx;
          oy = (ki - T[j].rho) * T[j].ny;
          oz = (ki - T[j].rho) * T[j].nz;
          T[j].ox = ox;
          T[j].oy = oy;
          T[j].oz = oz;
          /* Compute unit vectors b1 and b2 that span the tangent plane. */
          normal_vec (&b1x, &b1y, &b1z, T[j].nx, T[j].ny, T[j].nz);
          vec_prod (b2x, b2y, b2z, T[j].nx, T[j].ny, T[j].nz, b1x, b1y, b1z);
          T[j].b1x = b1x;
          T[j].b1y = b1y;
          T[j].b1z = b1z;
          T[j].b2x = b2x;
          T[j].b2y = b2y;
          T[j].b2z = b2z;
          /* Compute the 2D projection of the 3D Gauss of the reflection. */
          /* The symmetric 2x2 matrix N describing the 2D gauss. */
          n11 = L[i].m1 * b1x * b1x + L[i].m2 * b1y * b1y + L[i].m3 * b1z * b1z;
          n12 = L[i].m1 * b1x * b2x + L[i].m2 * b1y * b2y + L[i].m3 * b1z * b2z;
          n22 = L[i].m1 * b2x * b2x + L[i].m2 * b2y * b2y + L[i].m3 * b2z * b2z;
          /* The (symmetric) inverse matrix of N. */
          det_N = n11 * n22 - n12 * n12;
          inv_n11 = n22 / det_N;
          inv_n12 = -n12 / det_N;
          inv_n22 = n11 / det_N;
          /* The Cholesky decomposition of 1/2*inv_n (lower triangular L). */
          l11 = sqrt (inv_n11 / 2);
          l12 = inv_n12 / (2 * l11);
          l22 = sqrt (inv_n22 / 2 - l12 * l12);
          T[j].l11 = l11;
          T[j].l12 = l12;
          T[j].l22 = l22;
          det_L = l11 * l22;
          T[j].det_L = det_L;
          /* The product B^T D o. */
          Bt_D_O_x = b1x * L[i].m1 * ox + b1y * L[i].m2 * oy + b1z * L[i].m3 * oz;
          Bt_D_O_y = b2x * L[i].m1 * ox + b2y * L[i].m2 * oy + b2z * L[i].m3 * oz;
          /* Center of 2D Gauss in plane coordinates. */
          y0x = -(Bt_D_O_x * inv_n11 + Bt_D_O_y * inv_n12);
          y0y = -(Bt_D_O_x * inv_n12 + Bt_D_O_y * inv_n22);
          T[j].y0x = y0x;
          T[j].y0y = y0y;
          /* Factor alpha for the distance of the 2D Gauss from the origin. */
          alpha = L[i].m1 * ox * ox + L[i].m2 * oy * oy + L[i].m3 * oz * oz - (y0x * y0x * n11 + y0y * y0y * n22 + 2 * y0x * y0y * n12);
          T[j].refl = xsect_factor * det_L * exp (-alpha) / L[i].sig123; /* intensity of that Bragg */
          coh_refl += T[j].refl;                                         /* total scatterable intensity */
          /*some logic here to figure out cross-sections for NSF and SF scattering*/
          cdouble F;
          double F2, spin_up, spin_down;
          /*fractions of spin-up and spin-down incoming neutrons*/
          spin_up = (1 + scalar_prod (sx, sy, sz, mx, my, mz)) / 2;
          spin_down = (1 - scalar_prod (sx, sy, sz, mx, my, mz)) / 2;
          F2 = spin_up * (cabs (L[i].f[0]) * cabs (L[i].f[0]) + cabs (L[i].f[2]) * cabs (L[i].f[2]))
               + spin_down * (cabs (L[i].f[1]) * cabs (L[i].f[1]) + cabs (L[i].f[3]) * cabs (L[i].f[3]));
          // F2=cabs(F)*cabs(F);
          T[j].xsect = T[j].refl * F2;
          coh_xsect += T[j].xsect;
          j++;
        }
      } /* end for */
      tau_count = j;

      /* (3). Probabilities of the different possible interactions. */
      tot_xsect = abs_xsect + inc_xsect + coh_xsect;
      /* Cross-sections are in barns = 10**-28 m**2, and unit cell volumes are
         in AA**3 = 10**-30 m**2. Hence a factor of 100 is used to convert
         scattering lengths to m**-1 */
      coh_xlen = 1e2 * coh_xsect / V0;
      tot_xlen = 1e2 * tot_xsect / V0;

      /* (5). Transmission */
      p_trans = exp (-tot_xlen * l_full);
      if (!event_counter && p_transmit >= 0 && p_transmit <= 1) {
        mc_trans = p_transmit; /* first event */
      } else {
        mc_trans = p_trans;
      }
      mc_interact = 1 - mc_trans;
      if (mc_trans > 0 && (mc_trans >= 1 || rand01 () < mc_trans)) /* Transmit */
      {
        p *= p_trans / mc_trans;
        intersect = 0;
        break;
      }
      if (tot_xlen <= 0)
        ABSORB;

      if (mc_interact <= 0) /* Protect against rounding errors */
      {
        intersect = 0;
        break;
      }
      if (!event_counter)
        p *= fabs (1 - p_trans) / mc_interact;
      /* Select a point at which to scatter the neutron, taking
         secondary extinction into account. */
      /* dP(l) = exp(-tot_xlen*l)dl
         P(l<l_0) = [-1/tot_xlen*exp(-tot_xlen*l)]_0^l_0
                  = (1 - exp(-tot_xlen*l0))/tot_xlen
         l = -log(1 - tot_xlen*rand0max(P(l<l_full)))/tot_xlen
       */
      if (tot_xlen * l_full < 1e-6)
        /* For very weak scattering, use simple uniform sampling of scattering
           point to avoid rounding errors. */
        l = rand0max (l_full);
      else
        l = -log (1 - rand0max ((1 - exp (-tot_xlen * l_full)))) / tot_xlen;
      PROP_DT (l / v);
      event_counter++;

      /* (4). Account for the probability of sigma_abs */
      p *= (coh_xlen + inc_xlen) / tot_xlen;
      /* Choose between coherent and sigma_inc scattering */
      if (coh_xlen == 0 || rand0max (coh_xlen + inc_xlen) <= inc_xlen) {
        /* (6). sigma_inc scattering */
        randvec_target_circle (&kix, &kiy, &kiz, NULL, vx, vy, vz, 0);
        vx = kix; /* ki vector is used as tmp var with norm v */
        vy = kiy;
        vz = kiz; /* Go for next scattering event */
      } else {
        /* 7. Coherent scattering. Select reciprocal lattice point. */
        if (coh_refl <= 0)
          ABSORB;
        r = rand0max (coh_refl);
        sum = 0;
        for (j = 0; j < tau_count; j++) {
          sum += T[j].refl;
          if (sum > r)
            break;
        }
        if (j >= tau_count) {
          fprintf (stderr,
                   "Single_crystal: Error: Illegal tau search "
                   "(r = %g, sum = %g).\n",
                   r, sum);
          j = tau_count - 1;
        }
        i = T[j].index;
        /* (8). Pick scattered wavevector kf from 2D Gauss distribution. */
        z1 = randnorm ();
        z2 = randnorm ();
        y1 = T[j].l11 * z1 + T[j].y0x;
        y2 = T[j].l12 * z1 + T[j].l22 * z2 + T[j].y0y;
        kfx = T[j].rho_x + T[j].ox + T[j].b1x * y1 + T[j].b2x * y2;
        kfy = T[j].rho_y + T[j].oy + T[j].b1y * y1 + T[j].b2y * y2;
        kfz = T[j].rho_z + T[j].oz + T[j].b1z * y1 + T[j].b2z * y2;
        /* Normalize kf to length of ki, to account for planer
          approximation of the Ewald sphere. */
        adjust = ki / sqrt (kfx * kfx + kfy * kfy + kfz * kfz);
        kfx *= adjust;
        kfy *= adjust;
        kfz *= adjust;
        /* Adjust neutron weight (see manual for explanation). */
        p *= T[j].xsect * coh_refl / (coh_xsect * T[j].refl);

        vx = K2V * (L[i].u1x * kfx + L[i].u2x * kfy + L[i].u3x * kfz);
        vy = K2V * (L[i].u1y * kfx + L[i].u2y * kfy + L[i].u3y * kfz);
        vz = K2V * (L[i].u1z * kfx + L[i].u2z * kfy + L[i].u3z * kfz);

        /*adjust neutron pol vector*/
        double spin_up, spin_down, P_i, P_f;
        /*fractions of spin-up an0d spin-down for incoming neutrons*/
        P_i = scalar_prod (sx, sy, sz, mx, my, mz);
        spin_up = (1 + P_i) / 2;
        spin_down = (1 - P_i) / 2;
        /*from P-vec get fractions of up and down. New fractions of up and down will be weighted by up->up + down->up
          etc. Then we construct a new P-vector which obeys this*/
        double sigma_tot, polar;
        sigma_tot = spin_up * (cabs (L[i].f[0]) * cabs (L[i].f[0]) + cabs (L[i].f[2]) * cabs (L[i].f[2]))
                    + spin_down * (cabs (L[i].f[1]) * cabs (L[i].f[1]) + cabs (L[i].f[3]) * cabs (L[i].f[3]));
        polar = (spin_up * (cabs (L[i].f[0]) * cabs (L[i].f[0])) + spin_down * (cabs (L[i].f[3]) * cabs (L[i].f[3]))
                 - spin_down * (cabs (L[i].f[1]) * cabs (L[i].f[1])) - spin_up * (cabs (L[i].f[2]) * cabs (L[i].f[2])))
                / sigma_tot;
        if (polar < -1 || polar > 1) {
          fprintf (stderr, "Single_magnetic_crystal: inconsistent polarisation (%g %g %g %g)\n", spin_up, spin_down, sigma_tot, polar);
        }
        double ss = sqrt (scalar_prod (sx, sy, sz, sx, sy, sz));
        if (fabs (P_i) < FLT_EPSILON) {
          sx = mx * polar;
          sy = my * polar;
          sz = mz * polar;
        } else {
          sx *= (polar / P_i);
          sy *= (polar / P_i);
          sz *= (polar / P_i);
        }
      }
      SCATTER;
      /* exit if multiple scattering order has been reached */
      if (order && event_counter >= order) {
        intersect = 0;
        break;
      }
      /* Repeat loop for next scattering event. */
    } while (intersect); /* end do (intersect) (multiple scattering loop) */
  } /* if intersect */
%}

FINALLY
%{
  if (hkl_info.flag_warning)
    fprintf (stderr, "Single_crystal: %s: Error message was repeated %i times with absorbed neutrons.\n", NAME_CURRENT_COMP, hkl_info.flag_warning);
%}

MCDISPLAY
%{
  magnify ("xyz");
  if (hkl_info.shape == 0) { /* cylinder */
    circle ("xz", 0, yheight / 2.0, 0, radius);
    circle ("xz", 0, -yheight / 2.0, 0, radius);
    line (-radius, -yheight / 2.0, 0, -radius, +yheight / 2.0, 0);
    line (+radius, -yheight / 2.0, 0, +radius, +yheight / 2.0, 0);
    line (0, -yheight / 2.0, -radius, 0, +yheight / 2.0, -radius);
    line (0, -yheight / 2.0, +radius, 0, +yheight / 2.0, +radius);
  } else if (hkl_info.shape == 1) { /* box */
    double xmin = -0.5 * xwidth;
    double xmax = 0.5 * xwidth;
    double ymin = -0.5 * yheight;
    double ymax = 0.5 * yheight;
    double zmin = -0.5 * zdepth;
    double zmax = 0.5 * zdepth;
    multiline (5, xmin, ymin, zmin, xmax, ymin, zmin, xmax, ymax, zmin, xmin, ymax, zmin, xmin, ymin, zmin);
    multiline (5, xmin, ymin, zmax, xmax, ymin, zmax, xmax, ymax, zmax, xmin, ymax, zmax, xmin, ymin, zmax);
    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);
  } else if (hkl_info.shape == 2) { /* sphere */
    circle ("xy", 0, 0.0, 0, radius);
    circle ("xz", 0, 0.0, 0, radius);
    circle ("yz", 0, 0.0, 0, radius);
  } else if (hkl_info.shape == 3) { /* OFF file */
    off_display (offdata);
  }
%}
END