sics/fourlib.h

/*
                         F O U R L I B

 This is a library of routines for doing transformations between the 
 various coordinate systems used on a four circle diffractometer as 
 used for neutron or X-ray diffraction. The coordinate systems used are
 described in  Busing, Levy, Acta Cryst (1967),22, 457 ff.

 Generally we have:
 
 Z = [OM][CHI][PHI][UB]h

 where: Z is  a vector in the diffractometer coordinate system
        OM CHI PHI are rotation matrices around the respective angles
        UB is the UB matrix
        h is the reciprocal lattice vector. 

 The vector Z cannot only be expressed in terms of the angles stt, om, chi,
 and phi put also by polar coordinates gamma and nu.

 This code is a reimplementation based on a F77 code from Gary McIntyre, ILL
 and code extracted from the ILL MAD control program. 

  Mark Koennecke, November - December 2001
*/
#ifndef FOURLIB
#define FOURLIB
#include "matrix/matrix.h"

/**---------------------- PSD specific code -------------------------------
  Some diffractometers are equipped with position sensitive detectors.  
  The scheme to handle them implemented here is to convert the detector
  coordinates to polar coordinates gamma and nu and use the normal beam
  functions for any further work. A lot of information about the detector 
  is required for this transformation which is held in the datastruture 
  as defined below
  ----------------------------------------------------------------------*/ 

typedef struct{
  double xScale;     /* scale factor pixel --> mm for x */
  double yScale;     /* scale factor pixel --> mm for y */
  double distance;   /* distance sample detector in mm */
  double gamma;      /* gamma == two theta position of the detector */
  double nu;         /* tilt angle of the detector */
  int xZero;         /* x pixel coordinate of the zero point of the detector */
  int yZero;         /* y pixel coordinate of the zero point of the detector */
} psdDescription;
/**
 * det2pol converts the pixel coordinates x and y on the detector described
 * by psd to polar coordinates gamma and nu.
*/

void det2pol(psdDescription *psd, int x, int y, double *gamma, double *nu);

/**
 * pol2det converts the polar coordinates gamma and nu to detector coordinates
 *  x and y on the detector described psd
*/
void pol2det(psdDescription *psd, double gamma, double nu, int *x, int *y);

/*------------------------------------------------------------------------
  calculation of four circle diffractometer angles stt, om, chi and phi
  in order to put a reflection onto the equatorial plane and into a
  diffraction condition. In order to cope with angular restrictions 
  imposed for instance by sample environment devices or the eulerian 
  cradle itself various variations are implemented.

  These routines start at or target the vector Z1:

    Z1 = UB *h 
 ------------------------------------------------------------------------*/
/**
 *  calculate stt, om, chi and phi in order to put z1 into the bissecting
 *  diffraction condition. Returns 1 on success and 0 if z1 and lambda
 *  were invalid. The m version acts upon a matrix.
*/
int z1ToBisecting(double lambda, double z1[3], double *stt, double *om,
                                                 double *chi, double *phi);

int z1mToBisecting(double lambda, MATRIX z1, double *stt, double *om,
                                                 double *chi, double *phi);
/**
 * calculates diffraction angles to put z1 into a detector with a given
 *  offset in omega against the bisecting position. Useful for tweaking
 *  if omega is restricted.
*/
int z1ToAnglesWithOffset(double lambda, MATRIX z1,double omOffset,
                         double *stt,  double *om, 
                         double *chi, double *phi);
 
/**
 * calculates a PSI for a given offset in omega from the bisecting
 *  position. This can be useful for tweaking reflections to be 
 *  measurable at omega restricted diffractometers.
*/ 
int psiForOmegaOffset(MATRIX z1, double omOffset, 
			       double chi, double phi, double *psi);

/**
 * calculate new setting angles for a psi rotation. Input are the angles
 * calculated for a bisecting diffraction position. Output angles represent
 * the new psi setting. This is the method described by Busing & Levy as the
 * version when om == theta is psi = 0. 
*/
void rotatePsi(double om, double chi, double phi, double psi,
               double *newom, double *newchi, double *newphi);

/**
 * calculate z1 from angles stt, om, chi and phi. The angles must not describe
 * a bissecting position, however it is required that the angles describe a
 * reflection measured in the horizontal plane.
*/
void z1FromAngles(double lambda, double stt, double om, 
		  double chi, double phi, double z1[3]);
/*-----------------------------------------------------------------------
  Normal beam calculations. Here the diffraction vector is expressed as
  polar angles gamma and nu and omega. 
  -----------------------------------------------------------------------*/
/**
 * calculate the normal beam angles omeganb, gamma and nu from the four
 * circle angles twotheta, omega, chi and phi. 
*/
int bisToNormalBeam(double twotheta, double omega, double chi, double phi,
		double *omeganb, double *gamma, double *nu);
/**
 * calculate normal beam angles from z1
 */
int z1mToNormalBeam(double lambda, MATRIX z1, double *gamma, double *om, double *nu);

/**
 * calculate the vector z1 from the normal beam angles omega, gamma and nu.
 * chi and phi either do not exist or are 0. 
*/
void z1FromNormalBeam(double lambda, double omega, double gamma, 
		      double nu, double z1[3]);

/**
 * calculate z1 from four circle angles plus gamma, nu
*/
void z1FromAllAngles(double lambda, double omega, double gamma, 
		     double nu, double chi, double phi, double z1[3]);
/*------------------------------------------------------------------------
              Utility
-------------------------------------------------------------------------*/
/**
 * return val put into the 0 - 360 degree range
*/
double circlify(double val);

/**
 * converts a vector to a Matrix type
*/
MATRIX vectorToMatrix(double z[3]);

/**
 * converts Matrix zm to the vector z
*/
void matrixToVector(MATRIX zm, double z[3]);

/**
 * chi rotation matrix for chi into chim
*/
void chimat(MATRIX chim, double chi);

/**
 * phi rotation matrix for phi into phim
*/
void phimat(MATRIX phim, double phi);
/**
 * calcTheta calculates theta for a z1
 * returns 1 on success, 0 else
 */
int calcTheta(double lambda, MATRIX z1, double *d, double *theta);
 
#endif