- Many fixes to the triple axis stuff

* update after a1-a6 drive
  * intrduction of targets
- POLDI writing
- Moved HKL calculation 4 TRICS to fourlib
This commit is contained in:
cvs
2002-01-25 14:48:50 +00:00
parent 8c043c8cd1
commit 1e60f3be82
39 changed files with 3513 additions and 1160 deletions

807
hkl.c
View File

@ -10,6 +10,9 @@
Mark Koennecke, February 1998
Updated to use fourlib.
Mark Koennecke, December 2001
-----------------------------------------------------------------------------*/
#include <math.h>
#include <ctype.h>
@ -18,9 +21,15 @@
#include "motor.h"
#include "selector.h"
#include "selvar.h"
#include "fourlib.h"
#include "matrix/matrix.h"
#include "hkl.h"
#include "hkl.i"
/*
the space we leave in omega in order to allow for a scan to be done
*/
#define SCANBORDER 3.
/*-------------------------------------------------------------------------*/
static int HKLSave(void *pData, char *name, FILE *fd)
{
@ -359,318 +368,371 @@
self->iNOR = iNOB;
return 1;
}
/*-------------------------------------------------------------------------*/
static int hamth[8] = {0,0,1,1,0,0,1,1};
static int hamom[8] = {0,0,1,1,0,0,1,1};
static int hamch[8] = {0,1,1,1,1,1,1,0};
static int hamph[8] = {0,1,0,1,1,0,1,0};
static int hamhkl[8] = {0,1,0,1,0,1,0,1};
static const float RD = 57.2957795, pi = 3.1415926;
#define ABS(x) (x < 0 ? -(x) : (x))
/*-----------------------------------------------------------------------*/
static int checkBisecting(pHKL self,
double stt, double om, double chi, double phi)
{
int iTest;
float fHard, fLimit;
char pError[132];
/*-------------------------------------------------------------------------*/
static void rtan(double y, double x, double *rt)
{
if( (x == 0.) && (y == 0.) )
{
*rt = 0.;
return;
}
if( x == 0.)
{
*rt = pi/2.;
if(y < 0.)
{
*rt = -*rt;
}
return;
}
if(ABS(y) < ABS(x))
{
*rt = atan(ABS(y/x));
if(x < 0.)
{
*rt = pi - *rt;
}
if(y < 0.)
{
*rt = - *rt;
}
return;
}
else
{
*rt = pi/2 - atan(ABS(x/y));
}
if(x < 0.)
{
*rt = pi - *rt;
}
if( y < 0.)
{
*rt = - *rt;
}
return;
/* check two theta */
iTest = MotorCheckBoundary(self->pTheta,(float)stt, &fHard,pError,131);
if(!iTest)
{
return -1;
}
/* for omega check against the limits +- SCANBORDER in order to allow for
a omega scan
*/
MotorGetPar(self->pOmega,"softlowerlim",&fLimit);
if((float)om < fLimit + SCANBORDER){
iTest = 0;
} else {
iTest = 1;
MotorGetPar(self->pOmega,"softupperlim",&fLimit);
if((float)om > fLimit - SCANBORDER){
iTest = 0;
} else {
iTest = 1;
}
}
/* check chi and phi*/
iTest += MotorCheckBoundary(self->pChi,(float)chi, &fHard,pError,131);
iTest += MotorCheckBoundary(self->pPhi,(float)phi, &fHard,pError,131);
if(iTest == 3) /* none of them burns */
{
return 1;
}
return 0;
}
/*-----------------------------------------------------------------------*/
static int checkNormalBeam(double om, double gamma, double nu,
float fSet[4], SConnection *pCon, pHKL self)
{
int iTest;
char pError[132];
float fHard;
/* check omega, gamma and nu */
iTest = MotorCheckBoundary(self->pOmega,(float)om, &fHard,pError,131);
iTest += MotorCheckBoundary(self->pTheta,(float)gamma, &fHard,pError,131);
iTest += MotorCheckBoundary(self->pNu,(float)nu, &fHard,pError,131);
if(iTest == 3) /* none of them burns */
{
fSet[0] = (float)gamma;
fSet[1] = (float)om;
fSet[2] = (float)nu;
return 1;
}
return 0;
}
/*-----------------------------------------------------------------------
tryOmegaTweak tries to calculate a psi angle in order to put an
offending omega back into range.
-----------------------------------------------------------------------*/
static int tryOmegaTweak(pHKL self, MATRIX z1, double stt, double *om,
double *chi, double *phi){
int status;
float fLower, fUpper, omTarget, omOffset;
double dumstt, offom, offchi, offphi;
status = checkBisecting(self,stt,*om,*chi,*phi);
if(status < 0){
return 0; /* stt is burning */
} else if(status == 1){
return 1;
}
/*------------------------------------------------------------------------*/
static double fsign(double a, double a2)
/*
Is omega really the problem?
*/
omTarget = -9999;
MotorGetPar(self->pOmega,"softlowerlim",&fLower);
MotorGetPar(self->pOmega,"softupperlim",&fUpper);
if(*om < fLower + SCANBORDER) {
omTarget = fLower + SCANBORDER + .5;
}
if(*om > fUpper - SCANBORDER){
omTarget = fUpper - SCANBORDER - .5;
}
if(omTarget < -7000){
return 0;
}
/*
calculate omega offset
*/
omOffset = *om - omTarget;
omOffset = -omOffset;
/*
calculate angles with omega offset
*/
status = z1ToAnglesWithOffset(self->fLambda,z1, omOffset, &dumstt,
&offom, &offchi, &offphi);
if(!status){
return 0;
}
if(checkBisecting(self,dumstt,offom,offchi,offphi)){
*om = offom;
*chi = offchi;
*phi = offphi;
return 1;
}
return 0;
}
/*-----------------------------------------------------------------------*/
static MATRIX calculateScatteringVector(pHKL self, float fHKL[3])
{
MATRIX z1, hkl, ubm;
int i;
hkl = mat_creat(3,1,ZERO_MATRIX);
ubm = mat_creat(3,3,ZERO_MATRIX);
for(i = 0; i < 3; i++)
{
hkl[i][0] = fHKL[i];
ubm[0][i] = self->fUB[i];
ubm[1][i] = self->fUB[i+3];
ubm[2][i] = self->fUB[i+6];
}
z1 = mat_mul(ubm,hkl);
mat_free(ubm);
mat_free(hkl);
return z1;
}
/*---------------------------------------------------------------------*/
static int calculateBisecting(MATRIX z1, pHKL self, SConnection *pCon,
float fSet[4], double myPsi, int iRetry)
{
double stt, om, chi, phi, psi, ompsi, chipsi, phipsi;
int i, test;
/*
just the plain angle calculation
*/
if(!z1mToBisecting(self->fLambda,z1,&stt,&om,&chi,&phi))
{
if(a2 >= 0.)
return 0;
}
/*
check if angles in limits. If omega problem: try to tweak
*/
chi = circlify(chi);
phi = circlify(phi);
if(iRetry > 1)
{
if(tryOmegaTweak(self,z1,stt,&om,&chi,&phi)){
fSet[0] = (float)stt;
fSet[1] = (float)om;
fSet[2] = (float)circlify(chi);
fSet[3]= (float)circlify(phi);
return 1;
}
}
/*
if this does not work, try rotating through psi in order to
find a useful setting
*/
for(i = 0; i < iRetry; i++)
{
if(iRetry > 1)
{
return ABS(a);
psi = i*10.;
}
else
{
return -ABS(a);
psi = myPsi;
}
}
/*------------------------------------------------------------------------*/
static void sincos(double *sn, double *cs, double *loc)
{
if(ABS(*sn) >= 1.)
rotatePsi(om,chi,phi,psi,&ompsi,&chipsi,&phipsi);
chipsi = circlify(chipsi);
phipsi = circlify(phipsi);
test = checkBisecting(self,stt,ompsi,chipsi,phipsi);
if(test == 1)
{
*sn = fsign(1.0,*sn);
*cs = 0.;
return;
}
else
{
*cs = sqrt(1. - *sn * *sn);
return;
fSet[0] = (float)stt;
fSet[1] = (float)ompsi;
fSet[2] = (float)chipsi;
fSet[3]= (float)phipsi;
return 1;
}
}
/*-------------------------------------------------------------------------*/
static int ICAL(pHKL self, float fHKL[3], float fPsi, int iHamil,
int iQuad, float fSet[4], float fDelom)
/*
giving up!
*/
for(i = 0; i < 4; i++)
{
double z[3],tang[2],r,r1,phi,delomr,t,dstar,v,theta,omega,chi,cosp;
double sinp,d,q,zer,d2,all,den1,den2,cosps,a,sinc,cosc,b,c,chibi;
double sinps,sinchi,extrap, psi, loc;
int i,ip,ivariz,ipara,nh,nhamil,isign,icry,ii;
delomr = fDelom/(RD * -1);
iHamil -= 1; /* to get iHamil from 1 -- 8 & negative if not requested*/
psi = fPsi;
/* HKl --> Lab Axis through UB */
for(i = 0; i < 3; i++)
{
ii = 3*i;
z[i] = self->fUB[ii]*fHKL[0] + self->fUB[ii+1]*fHKL[1] +
self->fUB[ii+2]*fHKL[2];
}
fSet[i] = .0;
}
return 0;
}
/*----------------------------------------------------------------------*/
#define ABS(x) (x < 0 ? -(x) : (x))
/*-----------------------------------------------------------------------*/
static int calculateNormalBeam(MATRIX z1, pHKL self, SConnection *pCon,
float fSet[4], double myPsi, int iRetry)
{
int i, iTest;
double stt, om, chi, phi, gamma, nu, psi;
float currentPhi, currentChi;
double ompsi, chipsi, phipsi;
MATRIX chim, phim, z4, z3;
/* bissecting calculation! */
if(!self->iNOR)
{
/* four circle calculation */
t = z[0]*z[0] + z[1]*z[1] + z[2]*z[2];
dstar = sqrt(t);
for(i = 0; i < 3; i++)
{
z[i] /= dstar;
}
v = self->fLambda *.5 * dstar;
if(v > 1.)
{
return -1;
}
theta = 2 * asin(v);
omega = (.5 * theta) - delomr;
sinchi = -z[2]/cos(delomr);
if(ABS(sinchi) > 1.)
{
return -1;
}
chi = asin(sinchi);
if(iQuad == 1)
{
chi = pi - chi;
if(chi > pi)
{
chi = chi -(2.*pi);
}
}
d = delomr;
q = sin(d);
r = cos(d) * cos(chi);
cosp = (q*z[0] + r*z[1])/(q*q + r*r);
sinp = (z[0] - q*cosp)/r;
rtan(cosp,sinp,&phi);
chi = -chi;
if(fDelom != 0.)
{
/* recalculate psi */
r1 = -z[2];
if(r1 == 1.)
{
psi = phi;
}
else
{
chibi = asin(r1);
if(iQuad == 1)
{
chibi = pi - chibi;
if(chibi > pi)
{
chibi = chibi -(2*pi);
}
}
sinps = tan(chibi)*tan(delomr);
cosps = (cos(chibi) -sinps*sin(delomr)*sinchi)/cos(chi);
rtan(sinps,cosps,&psi);
psi = -psi;
}
psi = RD*psi;
goto hamil;
}
if(psi != 0.)
{
/* non zero psi */
a = psi/RD;
sinps = sin(a);
cosps = cos(a);
sinc = sin(chi);
cosc = cos(chi);
sinp = sin(phi);
cosp = cos(phi);
a = sinps*sinp - sinc*cosp*cosps;
b = -sinps*cosp - sinc*sinp*cosps;
/* do chi */
c = sqrt(a*a + b*b);
d = cosps*cosc;
rtan(c,d,&chi);
/* do phi */
rtan(-b,-a,&phi);
/* do omega */
a = -sinps *cosc;
rtan(a,sinc,&omega);
omega = omega + theta *.5;
}
hamil:
theta = theta *RD;
omega = omega *RD;
chi = chi*RD;
phi = phi*RD;
if(iHamil >= 0)
{
/* hamilton logic */
if(hamth[iHamil]) theta = theta * -1;
if(hamom[iHamil]) omega = omega -ABS(theta);
if(hamph[iHamil])
{
phi += 180.;
if(phi > 180.)
{
phi -= 360;
}
}
if(hamch[iHamil])
{
if( (iHamil == 1) || (iHamil == 6) )
{
chi *= -1;
}else if((iHamil == 1) || (iHamil == 5))
{
chi -= 180.;
if(chi < -180.)
{
chi += 360.;
}
}
else
{
chi = chi* -1 - 180;
if(chi > 180.)
{
chi -= 360;
}
}
}
}
/* hamhkl left out */
fSet[0] = theta;
fSet[1] = omega;
fSet[2] = chi;
fSet[3] = phi;
return 1;
/*
The usual condition for normal beam calculations is that both chi
and phi are 0. This is not the case at TRICS. Therefore we have to
multiply the scattering vector first with the chi and phi rotations
before we start.
*/
iTest = MotorGetSoftPosition(self->pChi,pCon,&currentChi);
if(iTest != 1)
{
return 0;
}
iTest = MotorGetSoftPosition(self->pPhi,pCon,&currentPhi);
if(iTest != 1)
{
return 0;
}
phim = mat_creat(3,3,ZERO_MATRIX);
phimat(phim,(double)currentPhi);
z4 = mat_mul(phim,z1);
chim = mat_creat(3,3,ZERO_MATRIX);
chimat(chim,(double)currentChi);
z3 = mat_mul(chim,z4);
mat_free(phim);
mat_free(chim);
mat_free(z4);
} /* end four circle */
else
{
/* start normal bissecting calculation */
/* ignore psi values */
d2 = 0.;
for(i = 0; i < 3; i++)
{
d2 += z[i]*z[i];
}
d = sqrt(d2);
d2 *= self->fLambda*self->fLambda;
tang[0] = z[2]*self->fLambda;
all = tang[0]*tang[0];
den1 = sqrt(d2-all);
if(all > 1.)
{
return -1;
}
den2 = sqrt(1. - all);
if(!((den1+den2 > 1.) && (ABS(den2-den1) < 1.)))
{
return -1;
}
isign = 1;
fSet[3] = 0.;
sincos(&tang[0],&tang[1],&loc);
fSet[2] = RD*atan2(tang[0],tang[1]);
tang[1] = (2-d2)/(2*den2);
sincos(&tang[1],&tang[0],&loc);
fSet[0] = RD*atan2(tang[0],tang[1])*isign;
tang[1] = d2/(2*den1);
sincos(&tang[1],&tang[0],&loc);
all = atan2(tang[0],tang[1]);
all -= isign*atan2(z[1],z[0]);
fSet[1] = 180. - RD*all -90.;
if(fSet[1] < 0.)
{
fSet[1] = 360. + fSet[1];
}
/*
do the bisecting angles first
*/
if(!z1mToBisecting(self->fLambda,z3,&stt,&om,&chi,&phi))
{
return 0;
}
return 1;
}
}
if(ABS(chi -90.) < .001 && ABS(phi-180.) < .001)
{
chi = .0;
phi = .0;
}
/*
in order to cope with all those limitations: rotate through psi
*/
for(i = 0; i < iRetry; i++)
{
if(iRetry > 1)
{
psi = 10. *i;
}
else
{
psi = myPsi;
}
rotatePsi(om,chi,phi,psi,&ompsi,&chipsi,&phipsi);
if(bisToNormalBeam(stt,ompsi,chipsi,phipsi,
&om, &gamma, &nu))
{
if(checkNormalBeam(om, gamma, nu,fSet,pCon,self))
{
return 1;
}
}
}
return 0;
}
/*---------------------------------------------------------------------*/
static int calculateNormalBeamOmega(MATRIX z1, pHKL self,
SConnection *pCon,
float fSet[4], double omOffset)
{
int iTest;
double stt, om, chi, phi, gamma, nu, psi;
float currentPhi, currentChi;
double ompsi, chipsi, phipsi;
MATRIX chim, phim, z4, z3;
/*
The usual condition for normal beam calculations is that both chi
and phi are 0. This is not the case at TRICS. Therefore we have to
multiply the scattering vector first with the chi and phi rotations
before we start.
*/
iTest = MotorGetSoftPosition(self->pChi,pCon,&currentChi);
if(iTest != 1)
{
return 0;
}
iTest = MotorGetSoftPosition(self->pPhi,pCon,&currentPhi);
if(iTest != 1)
{
return 0;
}
phim = mat_creat(3,3,ZERO_MATRIX);
phimat(phim,(double)currentPhi);
z4 = mat_mul(phim,z1);
chim = mat_creat(3,3,ZERO_MATRIX);
chimat(chim,(double)currentChi);
z3 = mat_mul(chim,z4);
mat_free(phim);
mat_free(chim);
mat_free(z4);
/*
do the bisecting angles first
*/
if(!z1ToAnglesWithOffset(self->fLambda,z3, omOffset, &stt,
&ompsi, &chi, &phi))
{
return 0;
}
if(ABS(chi -90.) < .001 && ABS(phi-180.) < .001)
{
chi = .0;
phi = .0;
}
if(bisToNormalBeam(stt,ompsi,chi,phi,
&om, &gamma, &nu))
{
if(checkNormalBeam(om, gamma, nu,fSet,pCon,self))
{
return 1;
}
}
return 0;
}
/*-------------------------------------------------------------------------
calculates the four circle settings. If the position can not be reached
because of a limit violation, then psi is rotated in 10 degree steps
until either the loop ends or we finally succeed. If there is a omega
violation we first try to calculate a delta omega which puts omega
into the right range. This is a fix because the omega movement is quite
often restricted due to the crygenic garbage around the sample.
often restricted due to the cryogenic garbage around the sample.
*/
int CalculateSettings(pHKL self, float fHKL[3], float fPsi, int iHamil,
float fSet[4], SConnection *pCon)
{
int iRet,iRetry, i;
int iQuad = 0;
int iTest;
float fDelom = 0.;
char pError[132];
char pBueffel[512];
float fHard, fVal, fUpper, fLower;
float myPsi = fPsi;
double myPsi = fPsi;
MATRIX z1;
double stt, om, chi, phi, ompsi, chipsi, phipsi;
int i,iRetry, status;
/* catch shitty input */
if( (fHKL[0] == 0.) && (fHKL[1] == 0.) && (fHKL[2] == 0.))
@ -681,19 +743,12 @@
}
/* some people are stupid.......... */
if(myPsi > 360.)
{
myPsi -= 360;
}
if(myPsi < .0)
{
myPsi += 360.;
}
/* no retries if normal beam calculation
or specific Hamilton or specific
psi requested */
if( (self->iNOR) || (iHamil != 0) || (myPsi > 0.1) )
myPsi = circlify(myPsi);
/*
no retries if specific psi requested.
*/
if((myPsi > 0.1) )
{
iRetry = 1;
}
@ -701,101 +756,35 @@
{
iRetry = 35;
}
/* loop till success */
for(i = 0, myPsi = fPsi; i < iRetry; i++, myPsi += 10.)
z1 = calculateScatteringVector(self,fHKL);
if(self->iNOR == 0)
{
/* just try it*/
iRet = ICAL(self,fHKL, myPsi, iHamil, self->iQuad,fSet,fDelom);
if(iRet < 0 ) /* could not do it */
{
sprintf(pBueffel,"ERROR: cannot calculate %4.1f %4.1f %4.1f",
status = calculateBisecting(z1,self,pCon,fSet, myPsi, iRetry);
}
else if(self->iNOR == 1)
{
status = calculateNormalBeam(z1,self,pCon,fSet, myPsi, iRetry);
}
else
{
myPsi = fPsi;
status = calculateNormalBeamOmega(z1,self,pCon,fSet, myPsi);
}
if(!status)
{
sprintf(pBueffel,"ERROR: cannot calculate %4.1f %4.1f %4.1f",
fHKL[0], fHKL[1], fHKL[2]);
SCWrite(pCon,pBueffel,eError);
return 0;
}
/* check two theta */
iTest = MotorCheckBoundary(self->pTheta,fSet[0], &fHard,pError,131);
if(!iTest)
{
/* this cannot be fixed */
sprintf(pBueffel,
"ERROR: %4.1f, %4.1f, %4.1f, %5.2f violates two theta limits",
fSet[0], fHKL[0], fHKL[1],fHKL[2]);
SCWrite(pCon,pBueffel,eError);
return 0;
}
/* check nu and omega if normal beam */
if(self->iNOR)
{
/* check omega */
iTest = MotorCheckBoundary(self->pOmega,fSet[1], &fHard,pError,131);
iTest += MotorCheckBoundary(self->pNu,fSet[2], &fHard,pError,131);
if(iTest != 2)
{
sprintf(pBueffel,
"ERROR: %4.1f, %4.1f, %4.1f, %5.2f %5.2f %5.2f violates nu limits",
fHKL[0], fHKL[1],fHKL[2],fSet[0], fSet[1],fSet[2]);
SCWrite(pCon,pBueffel,eError);
return 0;
}
else
{
return 1;
}
}
/* check chi and phi and omega, but first put into 0-360 degrees */
if(fSet[2] < 0.0)
{
fSet[2] = 360.0 + fSet[2];
}
/*
if(fSet[3] < 0.0)
{
fSet[3] = 360.0 + fSet[3];
}
*/
iTest = MotorCheckBoundary(self->pOmega,fSet[1], &fHard,pError,131);
if(iTest == 0 && i == 0)
{
/*
calculate a delta omega to put omega into center of range
*/
MotorGetPar(self->pOmega,"softupperlim",&fUpper);
MotorGetPar(self->pOmega,"softlowerlim",&fLower);
if(fSet[1] > fUpper)
{
fVal = fUpper - 5.; /* 5 degree safety for scanning */
}
else if (fSet[1] < fLower)
{
fVal = fLower + 5.; /* same */
}
fDelom = fSet[1] - fVal;
fDelom = -fDelom;
}
iTest += MotorCheckBoundary(self->pChi,fSet[2], &fHard,pError,131);
iTest += MotorCheckBoundary(self->pPhi,fSet[3], &fHard,pError,131);
if(iTest == 3) /* none of them burns */
{
sprintf(pBueffel,"Solution found at psi = %4.1f",myPsi);
SCWrite(pCon,pBueffel,eWarning);
return 1;
}
}
if(iRetry != 1)
{
sprintf(pBueffel,
"ERROR: failed to find a possible setting for %4.1f %4.1f %4.1f %s",
fHKL[0], fHKL[1], fHKL[2], "\n Even tried 36 psi settings");
SCWrite(pCon,pBueffel,eError);
return 0;
SCWrite(pCon,pBueffel,eError);
}
mat_free(z1);
return status;
return 0;
}
/*-------------------------------------------------------------------------*/
@ -1098,32 +1087,19 @@ ente:
static int angle2HKL(pHKL self ,double tth, double om,
double chi, double phi, float fHKL[3])
{
double dTh, dOm, dChi, dPhi, dHM;
MATRIX res, rez;
MATRIX rez, z1m;
double z1[3];
int i;
/* conversion to radians */
dTh = tth/(2*RD);
dOm = (om - tth/2.)/RD;
dChi = chi/RD;
dPhi = phi/RD;
dHM = 2* sin(dTh)/self->fLambda;
/* angles to XYZ, stolen from DIFRAC code in PRPXYZ */
res = mat_creat(3,1,ZERO_MATRIX);
res[0][0] = dHM*(cos(dChi)*cos(dPhi)*cos(dOm) -
sin(dPhi)*sin(dOm));
res[1][0] = dHM*(cos(dChi)*sin(dPhi)*cos(dOm) +
cos(dPhi)*sin(dOm));
res[2][0] = dHM*sin(dChi)*cos(dOm);
z1FromAngles(self->fLambda,tth,om,chi,phi,z1);
z1m = vectorToMatrix(z1);
/* multiply with UBinv in order to yield HKL */
rez = mat_mul(self->UBinv,res);
rez = mat_mul(self->UBinv,z1m);
for(i = 0; i < 3; i++)
{
fHKL[i] = rez[i][0];
}
mat_free(res);
fHKL[i] = (float)rez[i][0];
}
mat_free(z1m);
mat_free(rez);
return 1;
@ -1448,9 +1424,18 @@ ente:
return 0;
}
iRet = CalculateSettings(self,fHKL, fPsi, iHamil, fSet,pCon);
sprintf(pBueffel," theta = %f, omega = %f, chi = %f, phi = %f",
if(self->iNOR)
{
sprintf(pBueffel," gamma = %f, omega = %f, nu = %f",
fSet[0], fSet[1], fSet[2]);
SCWrite(pCon,pBueffel,eValue);
}
else
{
sprintf(pBueffel," theta = %f, omega = %f, chi = %f, phi = %f",
fSet[0], fSet[1], fSet[2],fSet[3]);
SCWrite(pCon,pBueffel,eValue);
SCWrite(pCon,pBueffel,eValue);
}
if(!iRet)
{
SCWrite(pCon,