diff --git a/Makefile b/Makefile
index 1fd458eb..ce812bc2 100644
--- a/Makefile
+++ b/Makefile
@@ -48,7 +48,8 @@ SOBJ = network.o ifile.o conman.o SCinter.o splitter.o passwd.o \
circular.o el755driv.o maximize.o sicscron.o tecsdriv.o sanscook.o \
tasinit.o tasutil.o t_rlp.o t_conv.o d_sign.o d_mod.o \
tasdrive.o tasscan.o synchronize.o definealias.o swmotor.o t_update.o \
- hmcontrol.o userscan.o slsmagnet.o rs232controller.o lomax.o
+ hmcontrol.o userscan.o slsmagnet.o rs232controller.o lomax.o \
+ polterwrite.o fourlib.o
MOTOROBJ = motor.o el734driv.o simdriv.o el734dc.o pipiezo.o pimotor.o
COUNTEROBJ = countdriv.o simcter.o counter.o
@@ -61,16 +62,17 @@ VELOOBJ = velo.o velosim.o velodorn.o velodornier.o
#----- comment or uncomment the following according to operating system
#------------- for Digital Unix
-#HDFROOT=/data/koenneck
-#CC=cc
-#EXTRA=
-#CFLAGS = -I$(HDFROOT)/include -Ihardsup -DHDF4 -DHDF5 -I. -std1 \
-# -g -warnprotos -c
+BINTARGET = bin
+HDFROOT=/data/lnslib
+CC=cc
+EXTRA=
+CFLAGS = -I$(HDFROOT)/include -Ihardsup -DHDF4 -DHDF5 -I. -std1 \
+ -g -warnprotos -c
#CFLAGS = -I$(HDFROOT)/include -DFORTIFY -DHDF4 -DHDF5 -Ihardsup -g \
# -std1 -warnprotos -c
-#LIBS = -L$(HDFROOT)/lib -Lhardsup -lhlib -Lmatrix -lmatrix -Ltecs \
-# -ltecsl -ltcl8.0 -lfor $(HDFROOT)/lib/libhdf5.a \
-# -lmfhdf -ldf $(HDFROOT)/lib/libjpeg.a -lz -lm -ll -lc
+LIBS = -L$(HDFROOT)/lib -Lhardsup -lhlib -Lmatrix -lmatrix -Ltecs \
+ -ltecsl -ltcl8.0 -lfor $(HDFROOT)/lib/libhdf5.a \
+ -lmfhdf -ldf $(HDFROOT)/lib/libjpeg.a -lz -lm -ll -lc
#------- for cygnus
#HDFROOT=../HDF411
@@ -81,16 +83,16 @@ VELOOBJ = velo.o velosim.o velodorn.o velodornier.o
# -lmfhdf -ldf -ljpeg -lz -lm
#---------- for linux
-BINTARGET=../bin
-HDFROOT=/usr/local
-CC=gcc
-CFLAGS = -I$(HDFROOT)/include -DHDF4 -Ihardsup -fwritable-strings \
- -DCYGNUS -DNONINTF -g -c
+#BINTARGET=../bin
+#HDFROOT=/usr/local
+#CC=gcc
+#CFLAGS = -I$(HDFROOT)/include -DHDF4 -Ihardsup -fwritable-strings \
+# -DCYGNUS -DNONINTF -g -c
#CFLAGS = -I$(HDFROOT)/include -Ihardsup -fwritable-strings -DFORTIFY \
# -DCYGNUS -DNONINTF -g -c
-LIBS= -L$(HDFROOT)/lib -Lhardsup -Ltecs -ltecsl -Lmatrix -lmatrix -lhlib \
- -ltcl8.0 -lmfhdf -ldf -ljpeg -lz -lm -lg2c -ldl
-EXTRA=nintf.o
+#LIBS= -L$(HDFROOT)/lib -Lhardsup -Ltecs -ltecsl -Lmatrix -lmatrix -lhlib \
+# -ltcl8.0 -lmfhdf -ldf -ljpeg -lz -lm -lg2c -ldl
+#EXTRA=nintf.o
#---------------------------------
.c.o:
diff --git a/choco.c b/choco.c
index fdf02d52..2279446c 100644
--- a/choco.c
+++ b/choco.c
@@ -174,7 +174,8 @@
*/
extern pCodri MakeSimChopper(void);
-extern pCodri MakeDoChoDriver(char *pHost, int iPort, int iChannel);
+extern pCodri MakeDoChoDriver(char *pHost, int iPort, int iChannel,
+ int iSingle);
extern pCodri MakeCookerDriver(char *pHost, int iPort, int iChannel);
/*-----------------------------------------------------------------------*/
int ChocoFactory(SConnection *pCon, SicsInterp *pSics, void *pData,
@@ -185,6 +186,7 @@ extern pCodri MakeCookerDriver(char *pHost, int iPort, int iChannel);
pObjectDescriptor pDes = NULL;
char pBueffel[132];
int iRet, iPort, iChannel;
+ int iSingle = 0;
if(argc < 3)
{
@@ -229,7 +231,19 @@ extern pCodri MakeCookerDriver(char *pHost, int iPort, int iChannel);
SCWrite(pCon,pBueffel,eError);
return 0;
}
- pDriv = MakeDoChoDriver(argv[3],iPort,iChannel);
+ if(argc > 7)
+ {
+ iRet = Tcl_GetInt(pSics->pTcl,argv[6],&iSingle);
+ if(iRet != TCL_OK)
+ {
+ sprintf(pBueffel,
+ "ERROR: expected integer as single flag, got %s",
+ argv[6]);
+ SCWrite(pCon,pBueffel,eError);
+ return 0;
+ }
+ }
+ pDriv = MakeDoChoDriver(argv[3],iPort,iChannel,iSingle);
}
else if(strcmp(argv[2],"sanscook") == 0)
{
diff --git a/crysconv.c b/crysconv.c
index e6e14dbf..64def425 100644
--- a/crysconv.c
+++ b/crysconv.c
@@ -534,3 +534,5 @@ L999:
return 0;
} /* erreso_ */
+
+
diff --git a/danu.dat b/danu.dat
index 12e2756e..3cecb0a9 100644
--- a/danu.dat
+++ b/danu.dat
@@ -1,3 +1,3 @@
- 123
+ 174
NEVER, EVER modify or delete this file
You'll risk eternal damnation and a reincarnation as a cockroach!|n
\ No newline at end of file
diff --git a/doc/user/tricsman b/doc/user/tricsman
index f45d0005..f3672c6c 100644
--- a/doc/user/tricsman
+++ b/doc/user/tricsman
@@ -344,6 +344,7 @@ H H L
%html tricspsd.htm 1
%html histogram.htm 2
%html nextrics.htm 2
+%html peaksearch.htm 2
%html trscan.htm 2
%html psddata.htm 1
diff --git a/doc/user/tricspsd.htm b/doc/user/tricspsd.htm
index db9eb037..c83b920a 100644
--- a/doc/user/tricspsd.htm
+++ b/doc/user/tricspsd.htm
@@ -5,12 +5,12 @@
TRICS with Position Sensitive Detectors
-There are no PSD's available yet, but the software is already ready. TRICS
- with a PSD requires the following special features.
+TRICS with a PSD requires the following special features.
diff --git a/docho.c b/docho.c
index aadddc32..f513f5da 100644
--- a/docho.c
+++ b/docho.c
@@ -13,6 +13,10 @@
Mark Koennecke, January 1999
+
+ Modified to support a single chopper only,
+
+ Uwe Filges, Mark Koennecke; November 2001
--------------------------------------------------------------------------*/
#include
#include
@@ -41,6 +45,8 @@
int iError;
int iBusy;
float fRatio;
+ int iSingle;
+ char pError[80];
} DoCho, *pDoCho;
/*
pHost, iPort and iChannel combined are the adress of the chopper
@@ -69,6 +75,9 @@
other parameters, its target value cannot be extracted from the chopper
status message.
+ iSingle is a flag which is true if only a single chopper is controlled
+ through this driver. This supports the POLDI single choper case.
+
*/
/*----------------------------------------------------------------------
ERROR CODES:
@@ -78,13 +87,64 @@
#define PARERROR -8004
#define BADSYNC -8005
#define BADSTOP -8006
+#define CHOPERROR -8007
+
+extern char *trim(char *pTrim); /* trim.c */
+
+/*----------------------------------------------------------------------*/
+static void SplitChopperReply(pCodri self, char *prefix, char *pBueffel)
+{
+ char pToken[30], pValue[20];
+ char *pPtr, *pTok, *pVal;
+ int iCount, iRet;
+ pDoCho pPriv = NULL;
+
+ pPriv = (pDoCho)self->pPrivate;
+
+ /* decompose pBueffel and store into string dictionary */
+ pPtr = strtok(pBueffel,";");
+ while(pPtr != NULL)
+ {
+ iCount = sscanf(pPtr,"%s %s",pToken,pValue);
+ if(iCount == 2)
+ {
+ pTok = trim(pToken);
+ pVal = trim(pValue);
+ sprintf(pToken,"%s.%s",prefix,pTok);
+ iRet = StringDictUpdate(pPriv->pPar,pToken,pVal);
+ if(!iRet)
+ {
+ StringDictAddPair(pPriv->pPar,pToken,pVal);
+ strcat(self->pParList,pToken);
+ strcat(self->pParList,",");
+ }
+ }
+ else
+ {
+ /* this fixes a bug with oversized messages in dphas */
+ if(strstr(pPtr,"dphas") != NULL)
+ {
+ sprintf(pToken,"%s.dphas",prefix);
+ iRet = StringDictUpdate(pPriv->pPar,
+ pToken,pPtr+5);
+ if(!iRet)
+ {
+ StringDictAddPair(pPriv->pPar,pToken,
+ pPtr+5);
+ strcat(self->pParList,pToken);
+ strcat(self->pParList,",");
+ }
+ }
+ }
+ pPtr = strtok(NULL,";");
+ }
+}
/*-------------------------------------------------------------------------
Well, DoChoStatus sends a status request to the Dornier chopper control
system. There is a gotcha, you need three reads to get the full information.
Then the answer is parsed and decomposed into parameter content for the
string dictionary. The single status components are separated by ;.
-------------------------------------------------------------------------*/
-extern char *trim(char *pTrim); /* trim.c */
static int DoChoStatus(pCodri self)
{
@@ -115,86 +175,21 @@ extern char *trim(char *pTrim); /* trim.c */
pPriv->iError = iRet;
return 0;
}
- /* decompose pBueffel and store into string dictionary */
- pPtr = strtok(pBueffel,";");
- while(pPtr != NULL)
- {
- iCount = sscanf(pPtr,"%s %s",pToken,pValue);
- if(iCount == 2)
- {
- pTok = trim(pToken);
- pVal = trim(pValue);
- sprintf(pToken,"chopper1.%s",pTok);
- iRet = StringDictUpdate(pPriv->pPar,pToken,pVal);
- if(!iRet)
- {
- StringDictAddPair(pPriv->pPar,pToken,pVal);
- strcat(self->pParList,pToken);
- strcat(self->pParList,",");
- }
- }
- else
- {
- /* this fixes a bug with oversized messages in dphas */
- if(strstr(pPtr,"dphas") != NULL)
- {
- iRet = StringDictUpdate(pPriv->pPar,
- "chopper1.dphas",pPtr+5);
- if(!iRet)
- {
- StringDictAddPair(pPriv->pPar,"chopper1.dphas",
- pPtr+5);
- strcat(self->pParList,"chopper1.dphas");
- strcat(self->pParList,",");
- }
- }
- }
- pPtr = strtok(NULL,";");
- }
+ SplitChopperReply(self,"chopper1",pBueffel);
- /* second send: get next second chopper line */
- iRet = SerialWriteRead(&(pPriv->pData),"",pBueffel,1023);
- if(iRet < 0)
+ if(!pPriv->iSingle)
{
- pPriv->iError = iRet;
- return 0;
- }
- /* decompose pBueffel and store into string dictionary */
- pPtr = strtok(pBueffel,";");
- while(pPtr != NULL)
- {
- iCount = sscanf(pPtr,"%s %s",pToken,pValue);
- if(iCount == 2)
+ /* second send: get next second chopper line */
+ iRet = SerialWriteRead(&(pPriv->pData),"",pBueffel,1023);
+ if(iRet < 0)
{
- pTok = trim(pToken);
- pVal = trim(pValue);
- sprintf(pToken,"chopper2.%s",pTok);
- iRet = StringDictUpdate(pPriv->pPar,pToken,pVal);
- if(!iRet)
- {
- StringDictAddPair(pPriv->pPar,pToken,pVal);
- strcat(self->pParList,pToken);
- strcat(self->pParList,",");
- }
+ pPriv->iError = iRet;
+ return 0;
}
- else
- {
- /* this fixes a bug with oversized messages in dphas */
- if(strstr(pPtr,"dphas") != NULL)
- {
- iRet = StringDictUpdate(pPriv->pPar,
- "chopper2.dphas",pPtr+5);
- if(!iRet)
- {
- StringDictAddPair(pPriv->pPar,"chopper2.dphas",
- pPtr+5);
- strcat(self->pParList,"chopper2.dphas");
- strcat(self->pParList,",");
- }
- }
- }
- pPtr = strtok(NULL,";");
+ SplitChopperReply(self,"chopper2",pBueffel);
}
+
+
return 1;
}
/*-------------------------------------------------------------------------*/
@@ -393,7 +388,16 @@ extern char *trim(char *pTrim); /* trim.c */
pPriv->iError = iRet;
return 0;
}
- pPriv->iError = 0;
+ if(strstr(pReply,"error") != NULL)
+ {
+ pPriv->iError = CHOPERROR;
+ strncpy(pPriv->pError,pReply,79);
+ return 0;
+ }
+ else
+ {
+ pPriv->iError = 0;
+ }
pPriv->iBusy = 1;
return 1;
}
@@ -610,6 +614,9 @@ extern char *trim(char *pTrim); /* trim.c */
case BADSYNC:
strncpy(pError,"Cannot drive slave chopper",iLen);
break;
+ case CHOPERROR:
+ strncpy(pError,pPriv->pError,iLen);
+ break;
case BADSTOP:
strncpy(pError,
"User called STOP. WARNING: chopper is still untamed!",
@@ -680,7 +687,7 @@ extern char *trim(char *pTrim); /* trim.c */
return CHFAIL;
}
/*-------------------------------------------------------------------------*/
- pCodri MakeDoChoDriver(char *pHost, int iPort, int iChannel)
+ pCodri MakeDoChoDriver(char *pHost, int iPort, int iChannel, int iSingle)
{
pCodri pNew = NULL;
pDoCho pPriv = NULL;
@@ -706,6 +713,7 @@ extern char *trim(char *pTrim); /* trim.c */
pPriv->iRefreshIntervall = 60;
pPriv->pPar = CreateStringDict();
pPriv->tRefresh = time(NULL);
+ pPriv->iSingle = iSingle;
if(!pPriv->pPar)
{
free(pText);
diff --git a/fourlib.c b/fourlib.c
new file mode 100644
index 00000000..63809e98
--- /dev/null
+++ b/fourlib.c
@@ -0,0 +1,791 @@
+/*
+ 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
+*/
+#include
+#include
+#include
+#include "fortify.h"
+#include "fourlib.h"
+
+#define PI 3.141592653589793
+#define RD 57.30
+#define ABS(x) (x < 0 ? -(x) : (x))
+
+/*------------------------------------------------------------------------
+ a safe routine for calculating the atan of y/x
+ ------------------------------------------------------------------------*/
+static double myatan(double y, double x){
+ if(ABS(x) < 0.0001){
+ if(y > .0){
+ return PI/2;
+ }else {
+ return -PI/2;
+ }
+ }
+
+ if(x > .0){
+ return atan(y/x);
+ } else {
+ if(y > .0){
+ return PI + atan(y/x);
+ } else {
+ return -PI + atan(y/x);
+ }
+ }
+}
+/*-------------------------------------------------------------------------*/
+static double rtan(double y, double x){
+ double val;
+
+ if( (x == 0.) && (y == 0.) ) {
+ return .0;
+ }
+ if( x == 0.) {
+ if(y < 0.){
+ return -PI/2.;
+ } else {
+ return PI/2.;
+ }
+ }
+ if(ABS(y) < ABS(x)) {
+ val = atan(ABS(y/x));
+ if(x < 0.) {
+ val = PI - val;
+ }
+ if(y < 0.){
+ val = -val;
+ }
+ return val;
+ } else {
+ val = PI/2. - atan(ABS(x/y));
+ if(x < 0.) {
+ val = PI - val;
+ }
+ if( y < 0.) {
+ val = - val;
+ }
+ }
+ return val;
+}
+/*----------------------------------------------------------------------
+ clear3x3 sets a 3x3 matrix to 0
+ -----------------------------------------------------------------------*/
+static void clear3x3(MATRIX target){
+ int i, j;
+ for(i = 0; i < 3; i++){
+ for(j = 0; j < 3; j++){
+ target[i][j] = .0;
+ }
+ }
+}
+/*---------------------------------------------------------------------*/
+MATRIX vectorToMatrix(double z[3]){
+ int i;
+ MATRIX res;
+
+ res = mat_creat(3,1,ZERO_MATRIX);
+ for(i = 0; i < 3; i++){
+ res[i][0] = z[i];
+ }
+ return res;
+}
+/*---------------------------------------------------------------------*/
+void matrixToVector(MATRIX zm, double z[3]){
+ int i;
+ for(i = 0; i < 3; i++){
+ z[i] = zm[i][0];
+ }
+}
+/*----------------------------------------------------------------------
+ A Busing & levy PSI matrix
+ ----------------------------------------------------------------------*/
+static void psimat(MATRIX target, double psi){
+ double psird = psi/RD;
+
+ clear3x3(target);
+
+ target[0][0] = 1.;
+ target[1][1] = cos(psird);
+ target[1][2] = sin(psird);
+ target[2][1] = -target[1][2];
+ target[2][2] = target[1][1];
+}
+/*----------------------------------------------------------------------
+ A Busing & levy CHI matrix
+ ----------------------------------------------------------------------*/
+void chimat(MATRIX target, double chi){
+ double chird = chi/RD;
+
+ clear3x3(target);
+
+ target[0][0] = cos(chird);
+ target[0][2] = sin(chird);
+ target[1][1] = 1.;
+ target[2][0] = -target[0][2];
+ target[2][2] = target[0][0];
+}
+/*----------------------------------------------------------------------
+ A Busing & levy PHI matrix
+ ----------------------------------------------------------------------*/
+void phimat(MATRIX target, double phi){
+ double phird = phi/RD;
+
+ clear3x3(target);
+
+ target[0][0] = cos(phird);
+ target[0][1] = sin(phird);
+ target[1][0] = -target[0][1];
+ target[1][1] = target[0][0];
+ target[2][2] = 1.;
+}
+/*------------------- PSD specific code ----------------------------*/
+void det2pol(psdDescription *psd, int x, int y, double *gamma, double *nu) {
+ double xobs, yobs, b, z, d;
+
+ assert(ABS(psd->distance ) > .001);
+
+ xobs = (x - psd->xZero)*psd->xScale;
+ yobs = (y - psd->yZero)*psd->yScale;
+
+ b = psd->distance*cos(yobs/psd->distance);
+ z = psd->distance*sin(yobs/psd->distance);
+ d = sqrt(xobs*xobs + b*b);
+ *gamma = psd->gamma + myatan(xobs,b)*RD;
+ *nu = psd->nu + myatan(z,d)*RD;
+}
+/*----------------------------------------------------------------------*/
+void pol2det(psdDescription *psd, double gamma, double nu, int *x, int *y){
+ double delga, delnu, td, tn, e, f, g, zobs, xobs;
+
+ delga = gamma - psd->gamma;
+ delnu = nu - psd->nu;
+ td = tan(delga/RD);
+ tn = tan(delnu/RD);
+ e = sqrt(1. + td*td);
+ f = tn*e;
+ g = psd->distance*(1. + tn*tn +tn*tn*td*td);
+ zobs = psd->distance*(atan(tn*e) + asin(f/g));
+ xobs = td*(psd->distance*cos(zobs/psd->distance));
+ *y = (int)nint(psd->yZero + zobs/psd->yScale);
+ *x = (int)nint(psd->xZero + xobs/psd->xScale);
+}
+
+/*--------------- bisecting geometry code -----------------------------*/
+
+/*
+ turn chi and phi in order to get Z1 into the equatorial plane
+*/
+static void turnEquatorial(MATRIX z1, double *chi, double *phi){
+ if(ABS(z1[0][0]) < .0001 || ABS(z1[1][0]) < .0001){
+ *phi = .0;
+ *chi = 90.;
+ if(z1[2][0] < .0){
+ *chi = - *chi;
+ }
+ } else {
+ *phi = myatan(z1[1][0],z1[0][0])*RD;
+ *chi = myatan(z1[2][0],
+ sqrt(z1[0][0]*z1[0][0] + z1[1][0]*z1[1][0]))*RD;
+ }
+}
+/*----------------------------------------------------------------------
+calculate d-spacing and theta from z1
+-----------------------------------------------------------------------*/
+static int calcTheta(double lambda, MATRIX z1, double *d, double *theta){
+ double dstar, sintheta;
+
+ dstar = sqrt(z1[0][0]*z1[0][0] + z1[1][0]*z1[1][0] + z1[2][0]*z1[2][0]);
+ if(dstar < .0001){
+ *d = .0;
+ *theta = 0;
+ return 0;
+ }
+
+ *d = 1./dstar;
+ sintheta = lambda * dstar/2.;
+ if(ABS(sintheta) > 1.0){
+ *d = .0;
+ *theta = .0;
+ return 0;
+ }
+ *theta = asin(sintheta)*RD;
+ return 1;
+}
+/*-------------------------------------------------------------------*/
+int z1ToBisecting(double lambda, double z1[3], double *stt, double *om,
+ double *chi, double *phi) {
+ MATRIX zz1;
+ int status;
+
+ zz1 = vectorToMatrix(z1);
+ status = z1mToBisecting(lambda,zz1,stt,om,chi,phi);
+ mat_free(zz1);
+ return status;
+}
+/*------------------------------------------------------------------------*/
+int z1mToBisecting(double lambda, MATRIX z1, double *stt, double *om,
+ double *chi, double *phi) {
+ double d;
+ int status;
+
+ status = calcTheta(lambda,z1,&d,om);
+ if(!status){
+ *stt = .0;
+ *om = .0;
+ *chi = .0;
+ *phi = .0;
+ return status;
+ }
+ turnEquatorial(z1,chi,phi);
+ *stt = *om*2.;
+ *chi = 180. - *chi;
+ *phi = 180. + *phi;
+
+ return 1;
+}
+/*---------------------------------------------------------------------*/
+int z1ToAnglesWithOffset(double lambda, MATRIX z1m,double omOffset,
+ double *stt, double *om,
+ double *chi, double *phi){
+ int status, i;
+ double d, delOm, sinchi, q, cosp, sinp, dstar;
+ MATRIX z1;
+
+ status = calcTheta(lambda,z1m,&d,om);
+ if(!status){
+ *stt = .0;
+ *om = .0;
+ *chi = .0;
+ *phi = .0;
+ return status;
+ }
+ *stt = 2. * *om;
+ *om += omOffset;
+
+ z1 = mat_copy(z1m); /* routine messes z1m up! */
+ dstar = sqrt(z1[0][0]*z1[0][0] + z1[1][0]*z1[1][0] + z1[2][0]*z1[2][0]);
+ for(i = 0; i < 3; i++){
+ z1[i][0] /= dstar;
+ }
+
+ delOm = -omOffset/RD;
+ sinchi = z1[2][0]/cos(delOm);
+ if(ABS(sinchi) > 1.){
+ mat_free(z1);
+ return 0;
+ }
+ *chi = asin(sinchi);
+ *chi = PI - *chi;
+ if(*chi > PI){
+ *chi = *chi - 2* PI;
+ }
+
+
+ q = cos(delOm)*cos(*chi);
+ cosp = (sin(delOm)*z1[0][0] + q*z1[1][0])/
+ (sin(delOm)*sin(delOm) + q*q);
+ sinp = (z1[0][0] - sin(delOm)*cosp)/ q;
+ *phi = rtan(cosp,sinp);
+
+ *phi *= RD;
+ *chi *= RD;
+
+ mat_free(z1);
+ return 1;
+}
+/*---------------------------------------------------------------------*/
+int psiForOmegaOffset(MATRIX z1m, double omOffset,
+ double chi, double phi, double *psi){
+ double chibi, sinps, cosps, sinchi, dstar;
+ MATRIX z1;
+
+ z1 = mat_copy(z1m); /* routine messes z1m up! */
+ dstar = sqrt(z1[0][0]*z1[0][0] + z1[1][0]*z1[1][0] + z1[2][0]*z1[2][0]);
+ z1[2][0] /= dstar;
+
+ omOffset /= RD;
+ if(-z1[2][0] == 1.){
+ *psi = phi;
+ } else {
+ chibi = PI - asin(-z1[2][0]);
+ if(chibi > PI){
+ chibi -= 2*PI;
+ }
+ sinps = tan(chibi)*tan(-omOffset);
+ sinchi = -z1[2][0]/cos(-omOffset);
+ if(ABS(sinchi) > 1.){
+ mat_free(z1);
+ return 0;
+ }
+ cosps = (cos(chibi) - sinps*sin(-omOffset)*sinchi)/cos(chi/RD);
+ *psi = -rtan(sinps,cosps)*RD;
+ }
+ mat_free(z1);
+ return 1;
+
+}
+
+/*----------------------------------------------------------------------*/
+void rotatePsi(double om, double chi, double phi, double psi,
+ double *newom, double *newchi, double *newphi){
+ MATRIX chim, phim, psim, r0, r0psi;
+
+ chim = mat_creat(3,3,ZERO_MATRIX);
+ chimat(chim,chi);
+ phim = mat_creat(3,3,ZERO_MATRIX);
+ phimat(phim,phi);
+ r0 = mat_mul(chim,phim);
+ psim = mat_creat(3,3,ZERO_MATRIX);
+ psimat(psim,psi);
+ r0psi = mat_mul(psim,r0);
+
+ *newchi = rtan(sqrt(r0psi[2][0]*r0psi[2][0] + r0psi[2][1]*r0psi[2][1]),
+ r0psi[2][2])*RD;
+ *newphi = rtan(-r0psi[2][1],-r0psi[2][0])*RD;
+ *newom = om + rtan(-r0psi[1][2],r0psi[0][2])*RD;
+ mat_free(chim);
+ mat_free(phim);
+ mat_free(psim);
+ mat_free(r0);
+ mat_free(r0psi);
+}
+/*--------------------------------------------------------------------
+ calculate Z1 = [PHI]TZ3
+ The T means transposed matrixes for the return calculation!
+ ---------------------------------------------------------------------*/
+static void z1fromz2(double z1[3], MATRIX z2,
+ double phi){
+ MATRIX phim, phimt, zz1;
+
+ phim = mat_creat(3,3,ZERO_MATRIX);
+ phimat(phim, phi);
+ phimt = mat_tran(phim);
+ zz1 = mat_mul(phimt,z2);
+ matrixToVector(zz1,z1);
+ mat_free(phim);
+ mat_free(phimt);
+ mat_free(zz1);
+}
+/*--------------------------------------------------------------------
+ calculate Z1 = [PHI]T*[CHI]T*Z3
+ The T means transposed matrixes for the return calculation!
+ ---------------------------------------------------------------------*/
+static void z1fromz3(double z1[3], MATRIX z3, double chi,
+ double phi){
+ MATRIX chim, chimt, z2;
+
+ chim = mat_creat(3,3,ZERO_MATRIX);
+ chimat(chim, chi);
+ chimt = mat_tran(chim);
+ z2 = mat_mul(chimt,z3);
+ z1fromz2(z1,z2,phi);
+ mat_free(chim);
+ mat_free(chimt);
+ mat_free(z2);
+}
+/*--------------------------------------------------------------------
+ calculate Z1 = [PHI]T*[CHI]T*[OM]T*Z4
+ The T means transposed matrixes for the return calculation!
+ ---------------------------------------------------------------------*/
+static void z1fromz4(double z1[3], MATRIX z4, double om, double chi,
+ double phi){
+ MATRIX oma, oma2, z3;
+
+ oma = mat_creat(3,3,ZERO_MATRIX);
+ phimat(oma, om);
+ oma2 = mat_tran(oma);
+ z3 = mat_mul(oma2,z4);
+ z1fromz3(z1,z3,chi,phi);
+ mat_free(oma);
+ mat_free(oma2);
+ mat_free(z3);
+}
+/*---------------------------------------------------------------------*/
+void z1FromAngles(double lambda, double stt, double om,
+ double chi, double phi, double z1[3]){
+ MATRIX z4;
+ double th;
+
+ z4 = mat_creat(3,1,ZERO_MATRIX);
+ th = (stt/2.)/RD;
+ z4[0][0] = (2. * sin(th)*cos(th))/lambda;
+ z4[1][0] = (-2. *sin(th)*sin(th))/lambda;
+ z4[2][0] = .0;
+ z1fromz4(z1,z4,om,chi,phi);
+ mat_free(z4);
+}
+/*----------------------------------------------------------------------
+ Normal Beam geometry specific calculations. This means the reflection is
+ expressed through the three angles omega, gamma (two theta detector) and
+ nu (detector tilt).
+
+ bisToNormalBeam was lifted from HKLGEN.
+ -----------------------------------------------------------------------*/
+double sign(double a, double b){
+ if(b >= .0){
+ return ABS(a);
+ } else {
+ -ABS(a);
+ }
+}
+/*----------------------------------------------------------------------*/
+int bisToNormalBeam(double twotheta, double omega, double chi, double phi,
+ double *omeganb, double *gamma, double *nu){
+ double tth, fom, fchi, fphi, sint, sinnu, del;
+
+ tth = twotheta/RD;
+ fom = omega/RD;
+ fchi = chi/RD;
+ fphi = phi/RD;
+
+
+ /* nu calculation */
+ sint = sin(tth/2.);
+ *nu = 2. * sin(fchi)*sint;
+ if(*nu > 1.){
+ return 0;
+ }
+ *nu = asin(*nu);
+ *nu = ABS(*nu)*sign(1.,fchi);
+ sinnu = sin(*nu);
+
+ /* gamma calculation */
+ *gamma = cos(tth)/cos(*nu);
+ if(ABS(*gamma) > 1.0){
+ return 0;
+ }
+ *gamma = acos(*gamma);
+ *gamma = *gamma * sign(1.,tth);
+
+ /* omega normal beam */
+ del = asin(sint/cos(fchi));
+ *omeganb = del + fphi;
+
+ *omeganb *= RD;
+ *gamma *= RD;
+ *nu *= RD;
+
+ if(ABS(*omeganb) > 180.){
+ if(*omeganb < -180.){
+ *omeganb += 360.;
+ }
+ if(*omeganb > 180.){
+ *omeganb -= 360.;
+ }
+ }
+ return 1;
+}
+/* --------------------------------------------------------------------*/
+static void z4FromNormalBeam(MATRIX z4, double lambda, double gamma,
+ double nu){
+
+ double gar, nur;
+
+ gar = gamma/RD;
+ nur = nu/RD;
+
+ /*
+ this is basically a conversion from polar coordinates to
+ x,y,z coordinates. However, sin and cos are exchanged in
+ comparison to textbook formulas for this conversion. But this
+ is only a shift of origin
+ */
+ z4[0][0] = (sin(gar) * cos(nur))/lambda;
+ z4[1][0] = (cos (gar) * cos(nur) -1.)/lambda;
+ z4[2][0] = (sin(nur))/lambda;
+}
+/*----------------------------------------------------------------------*/
+void z1FromNormalBeam(double lambda, double omega, double gamma,
+ double nu, double z1[3]){
+ MATRIX z4, omm, omt, z1m;
+
+ z4 = mat_creat(3,1,ZERO_MATRIX);
+ z4FromNormalBeam(z4,lambda, gamma,nu);
+ omm = mat_creat(3,3,ZERO_MATRIX);
+ phimat(omm,omega);
+ omt = mat_tran(omm);
+ z1m = mat_mul(omt,z4);
+ matrixToVector(z1m,z1);
+ mat_free(z4);
+ mat_free(omm);
+ mat_free(omt);
+ mat_free(z1m);
+}
+/*--------------------------------------------------------------------*/
+double circlify(double val){
+ while(val > 360.){
+ val -= 360.;
+ }
+ while(val < 0.){
+ val += 360.;
+ }
+ return val;
+}
+/*----------------------------------------------------------------------*/
+void z1FromAllAngles(double lambda, double omega , double gamma,
+ double nu, double chi, double phi, double z1[3]){
+ MATRIX z3;
+
+ z1FromNormalBeam(lambda, omega, gamma, nu, z1);
+ z3 = vectorToMatrix(z1);
+ z1fromz3(z1,z3,chi,phi);
+ mat_free(z3);
+}
+/*------------------- a test program ------------------------------------*/
+#ifdef TESTCODE
+int main(int argc, char *argv[]){
+ psdDescription psd;
+ double gamma,nu, lambda, stt, om, chi, phi, z1[3], psi;
+ double psiom, psichi, psiphi;
+ double sttex, omex, chiex, phiex;
+ int x, y, status, i;
+ MATRIX UB, UBinv, zz1,zz1b, hkl, hklback;
+ double gaex, omnbex, nuex, omeganb;
+
+ /* initializations */
+ UB = mat_creat(3,3,ZERO_MATRIX);
+ UB[0][0] = 0.016169;
+ UB[0][1] = 0.011969;
+ UB[0][2] = 0.063195;
+ UB[1][0] = -.000545;
+ UB[1][1] = 0.083377;
+ UB[1][2] = -0.009117;
+ UB[2][0] = -0.162051;
+ UB[2][1] = 0.000945;
+ UB[2][2] = 0.006321;
+ UBinv = mat_inv(UB);
+ lambda = 1.179;
+
+
+ /* test PSD code */
+ psd.xScale = .7;
+ psd.yScale = 1.4;
+ psd.xZero = 128;
+ psd.yZero = 64;
+ psd.distance = 450.;
+ psd.gamma = 43.;
+ psd.nu = 12.3;
+
+ det2pol(&psd,200, 111,&gamma, &nu);
+ pol2det(&psd,gamma,nu, &x, &y);
+ if(x == 200 && y == 111){
+ printf("PSD test: OK\n");
+ }
+
+ /* test bisecting calculations. The tests were taken from a system known
+ to work
+ */
+ status = 1;
+ sttex = 30.58704;
+ omex = 15.293520;
+ chiex = 129.053604;
+ phiex = 134.191132;
+ hkl = mat_creat(3,1,ZERO_MATRIX);
+ hkl[0][0] = -2.;
+ hkl[1][0] = -2.;
+ hkl[2][0] = 4.;
+ zz1 = mat_mul(UB,hkl);
+ z1mToBisecting(lambda,zz1,&stt,&om,&chi,&phi);
+ if(ABS(stt - sttex) > .01 ||
+ ABS(om - omex) > .01 ||
+ ABS(chi - chiex) > .01 ||
+ ABS(phi - phiex) > .01){
+ printf("Bisecting Angles calculation FAILED!!!!!!!!!\n");
+ printf("Expected: %12.4f %12.4f %12.4f %12.4f\n", sttex, omex, chiex,phiex);
+ printf("Got: %12.4f %12.4f %12.4f %12.4f\n", stt,om,chi,phi);
+ status = 0;
+ }
+ z1FromAngles(lambda,sttex,omex,chiex,phiex,z1);
+ for(i = 0; i < 3; i++){
+ if(ABS(zz1[i][0] - z1[i]) > .0001){
+ printf("Calculation of Z1 from Angles FAILED!!!!\n");
+ printf("Expected: %8.4f %8.4f %8.4f\n", zz1[0][0], zz1[1][0],
+ zz1[2][0]);
+ printf("GOT: %8.4f %8.4f %8.4f\n",z1[0],z1[1],z1[2]);
+ status = 0;
+ }
+ }
+ hklback = mat_mul(UBinv,zz1);
+ for(i = 0; i < 3; i++){
+ if(ABS(hklback[i][0] - hkl[i][0]) > .0001){
+ printf("Back Calculation of HKL FAILED!!!!!!!!!\n");
+ printf("Expected: %8.4f %8.4f %8.4f\n", hkl[0][0], hkl[1][0],
+ hkl[2][0]);
+ printf("GOT: %8.4f %8.4f %8.4f\n",hklback[0][0],
+ hklback[1][0],hklback[2][0]);
+ status = 0;
+ }
+ }
+ if(status == 1){
+ printf("Bisecting test: OK\n");
+ }
+
+
+
+ /* try turning in psi */
+ rotatePsi(omex,chiex,phiex,30.,&psiom,&psichi,&psiphi);
+ omex = 37.374298;
+ chiex = 123.068192;
+ phiex = 170.8209099;
+ if(ABS(psiom - omex) > .01 ||
+ ABS(psichi - chiex) > .01 ||
+ ABS(psiphi - phiex) > .01){
+ printf("Psi Rotation test FAILED!!!!!!!!!\n");
+ printf("Expected: %12.4f %12.4f %12.4f\n", omex, chiex,phiex);
+ printf("Got: %12.4f %12.4f %12.4f\n", psiom,psichi,psiphi);
+ status = 0;
+ } else {
+ printf("Psi rotation test: OK\n");
+ }
+
+ /*
+ try to calculate a psi for a given omega offset
+ */
+ omex = 15.293520;
+ chiex = 129.053604;
+ phiex = 134.191132;
+ if(psiForOmegaOffset(zz1,-5.,chiex, phiex, &psi) != 1){
+ printf("Failed to calculate PSI for a Omega Offset\n");
+ } else {
+ rotatePsi(omex,chiex,phiex,psi,&psiom,&psichi,&psiphi);
+ if(ABS(psiom +5. - omex) > .5){
+ printf("Failed to reach desired omega with calculated psi\n");
+ printf("Expected: %12.4f Got: %12.4f PSI = %12.4f\n", omex -10.,
+ psiom, psi);
+ } else {
+ printf("Psi for Omega Offset: OK\n");
+ }
+ }
+
+ /*
+ try to calculate angles with a omega offset
+ */
+ omex = 10.29;
+ chiex = 128.78;
+ phiex = 126.24;
+ sttex = 30.59;
+ z1ToAnglesWithOffset(lambda, zz1, -5., &stt, &om, &chi, &phi);
+ if(ABS(stt - sttex) > .01 ||
+ ABS(om - omex) > .01 ||
+ ABS(chi - chiex) > .01 ||
+ ABS(phi - phiex) > .01){
+ printf("Angles calculation with Omega Offset FAILED!!!!!!!!!\n");
+ printf("Expected: %12.4f %12.4f %12.4f %12.4f\n", sttex, omex, chiex,phiex);
+ printf("Got: %12.4f %12.4f %12.4f %12.4f\n", stt,om,chi,phi);
+ status = 0;
+ } else {
+ printf("Angles with Omega Offset: OK\n");
+ }
+
+
+ /*
+ test normal beam codes
+ */
+ omex = 37.374298;
+ gaex = 19.322;
+ omnbex = -21.057;
+ nuex = 24.186;
+ z1mToBisecting(lambda,zz1,&stt,&om,&chi,&phi);
+ chi = 50.949;
+ phi = -45.8088;
+ if(!bisToNormalBeam(stt,om,chi,phi,
+ &omeganb, &gamma, &nu)){
+ printf("Error on normal beam calculation!!\n");
+ }
+ if(ABS(omeganb - omnbex) > .01 ||
+ ABS(gamma - gaex) > .01 ||
+ ABS(nu - nuex) > .01){
+ printf("Normal Beam Calculation test FAILED!!!!!!!!!\n");
+ printf("Expected: %12.4f %12.4f %12.4f\n", omnbex, gaex,nuex);
+ printf("Got: %12.4f %12.4f %12.4f\n", omeganb,gamma,nu);
+ } else {
+ printf("Normal Beam Angle Calculation: OK\n");
+ }
+
+ z1FromNormalBeam(lambda, omnbex, gaex, nuex,z1);
+ status = 1;
+ for(i = 0; i < 3; i++){
+ if(ABS(zz1[i][0] - z1[i]) > .0001){
+ printf("Calculation of Z1 from Normal Beam Angles FAILED!!!!\n");
+ printf("Expected: %8.4f %8.4f %8.4f\n", zz1[0][0], zz1[1][0],
+ zz1[2][0]);
+ printf("GOT: %8.4f %8.4f %8.4f\n",z1[0],z1[1],z1[2]);
+ break;
+ status = 0;
+ }
+ }
+ if(status){
+ printf("Normal Beam Backwards Calculation: OK\n");
+ }
+
+ /*
+ for interest: try to see how normal beam angles go if we rotate psi
+ */
+ /*
+ sttex = 30.58704;
+ om = 15.293520;
+ chi = 50.949;
+ phi = -45.8088;
+ printf(" PSI Omega Gamma Nu\n");
+ for(i = 0; i < 36; i++){
+ psi = i*10.;
+ rotatePsi(om,chi,phi,psi,&psiom,&psichi,&psiphi);
+ status = bisToNormalBeam(sttex,psiom,psichi,psiphi,&omeganb, &gamma, &nu);
+ if(status){
+ printf("%12.4f%12.4f%12.4f%12.4f\n",psi,omeganb,gamma,nu);
+ }
+ }
+ */
+
+}
+
+#endif
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/fourlib.h b/fourlib.h
new file mode 100644
index 00000000..55cc7e5a
--- /dev/null
+++ b/fourlib.h
@@ -0,0 +1,173 @@
+/*
+ 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 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);
+#endif
+
+
+
+
+
+
+
+
diff --git a/hardsup/Makefile b/hardsup/Makefile
index a93165ef..e3d1f462 100644
--- a/hardsup/Makefile
+++ b/hardsup/Makefile
@@ -11,11 +11,11 @@ OBJ= el734_utility.o asynsrv_utility.o stredit.o \
makeprint.o StrMatch.o
#---------- for Redhat linux
-CC= gcc
-CFLAGS= -I/usr/local/include -I. -I../ -DLINUX -g -c
+#CC= gcc
+#CFLAGS= -I/usr/local/include -I. -I../ -DLINUX -g -c
#------------ for DigitalUnix
-#CC=cc
-#CFLAGS= -I/data/koenneck/include -I. -I../ -std1 -g -c
+CC=cc
+CFLAGS= -I/data/koenneck/include -I. -I../ -std1 -g -c
#------------ for DigitalUnix with Fortify
## CC=cc
## CFLAGS= -DFORTIFY -I. -I../ -std1 -g -c
diff --git a/hkl.c b/hkl.c
index 4a5158fa..89b19ad4 100644
--- a/hkl.c
+++ b/hkl.c
@@ -10,6 +10,9 @@
Mark Koennecke, February 1998
+ Updated to use fourlib.
+
+ Mark Koennecke, December 2001
-----------------------------------------------------------------------------*/
#include
#include
@@ -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,¤tChi);
+ if(iTest != 1)
+ {
+ return 0;
+ }
+ iTest = MotorGetSoftPosition(self->pPhi,pCon,¤tPhi);
+ 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,¤tChi);
+ if(iTest != 1)
+ {
+ return 0;
+ }
+ iTest = MotorGetSoftPosition(self->pPhi,pCon,¤tPhi);
+ 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,
diff --git a/lomax.c b/lomax.c
index 99612754..59fe81d8 100644
--- a/lomax.c
+++ b/lomax.c
@@ -103,6 +103,7 @@ static int testMaximum(int *iData, int xsize, int ysize,
{
int testValue, x, y, half;
int *iPtr;
+ int equalCount = 0;
testValue = iData[j * xsize + i];
half = window/2;
@@ -114,8 +115,16 @@ static int testMaximum(int *iData, int xsize, int ysize,
{
if(iPtr[x] > testValue)
return 0;
+ if(iPtr[x] == testValue)
+ equalCount++;
}
}
+ /*
+ if(equalCount > 3)
+ {
+ return 0;
+ }
+ */
return 1;
}
/*--------------------------------------------------------------------*/
@@ -144,7 +153,7 @@ int testLocalMaximum(int *iData, int xsize, int ysize,
}
/*-------------------------------------------------------------------*/
int calculateCOG(int *iData, int xsize, int ysize,
- int *i, int *j, int *intensity,
+ int *i, int *j, int *intensity,int *nCount,
int cogWindow,
float contour)
{
@@ -153,6 +162,9 @@ int calculateCOG(int *iData, int xsize, int ysize,
float cogTotal, cogX, cogY;
int *iPtr;
+ if(!testBoundaries(xsize,ysize,cogWindow,*i,*j))
+ return 0;
+
/*
preparations
*/
@@ -175,6 +187,7 @@ int calculateCOG(int *iData, int xsize, int ysize,
/*
build the sums
*/
+ *nCount = 0;
cogTotal = cogY = cogX = .0;
for(y = yLow; y < yMax; y++)
{
@@ -183,13 +196,14 @@ int calculateCOG(int *iData, int xsize, int ysize,
{
if(iPtr[x] > threshold)
{
+ *nCount++;
cogTotal += iPtr[x];
cogY += y * iPtr[x];
cogX += x * iPtr[x];
}
}
}
- if(cogTotal < .0)
+ if(cogTotal <= .0)
{
return 0;
}
@@ -200,6 +214,82 @@ int calculateCOG(int *iData, int xsize, int ysize,
return 1;
}
+/*-------------------------------------------------------------------*/
+void calculateStatistics(int *iData, int xsize, int ysize,
+ float *average, float *maximum)
+{
+ int i, iLength;
+ int max = -9999999999999;
+ long sum = 0;
+
+ iLength = xsize*ysize;
+ for(i = 0; i < iLength; i++)
+ {
+ sum += iData[i];
+ if(iData[i] > max)
+ max = iData[i];
+ }
+ *average = (float)sum/(float)iLength;
+ *maximum = (float)max;
+}
+/*-------------------------------------------------------------------*/
+int wellFormed(int *iData, int xsize, int ysize,
+ int i, int j, int window, float contour,
+ int maxBad)
+{
+ int testValue, x, y, half;
+ int *iPtr;
+ int badCount = 0;
+
+
+ testValue = (int)((float)iData[j * xsize + i]*contour);
+ half = window/2;
+
+ /*
+ test upper row
+ */
+ iPtr = iData + (j - half) * xsize + i - half;
+ for(x = 0; x < window; x++)
+ {
+ if(iPtr[x] > testValue)
+ badCount++;
+ }
+
+ /*
+ test lower row
+ */
+ iPtr = iData + (j + half) * xsize + i - half;
+ for(x = 0; x < window; x++)
+ {
+ if(iPtr[x] > testValue)
+ badCount++;
+ }
+
+ /*
+ test columns
+ */
+ for(y = j - half; y < j + half; y++)
+ {
+ /*
+ left
+ */
+ if(iData[y*xsize + i - half] > testValue)
+ badCount++;
+ /*
+ right
+ */
+ if(iData[y*xsize + i + half] > testValue)
+ badCount++;
+ }
+
+ if(badCount > maxBad)
+ {
+ return 0;
+ }
+
+ return 1;
+}
+
/*-------------------------------------------------------------------*/
static int checkHM(pHistMem *pHM, SicsInterp *pSics, SConnection *pCon,
char *name, int *iDim)
@@ -247,7 +337,7 @@ int LoMaxAction(SConnection *pCon, SicsInterp *pSics,
{
pLoMax self = NULL;
char pBueffel[256], pNum[20];
- int iRet, i, j, intensity;
+ int iRet, i, j, intensity, count, badMax;
int iDim[10], nDim;
pHistMem pHM = NULL;
CommandList *pCom = NULL;
@@ -256,6 +346,7 @@ int LoMaxAction(SConnection *pCon, SicsInterp *pSics,
int *iData;
double dVal;
ObPar *ob = NULL;
+ float average, maximum;
self = (pLoMax)pData;
assert(pCon);
@@ -291,9 +382,9 @@ int LoMaxAction(SConnection *pCon, SicsInterp *pSics,
return 0;
}
Tcl_DStringInit(&result);
- Tcl_DStringStartSublist(&result);
iData = GetHistogramPointer(pHM,pCon);
window = (int)ObVal(self->pParam,WINDOW);
+ count = 0;
for(i = 0 + window/2; i < iDim[0] - window/2; i++)
{
for(j = 0 + window/2; j < iDim[1] - window/2; j++)
@@ -305,18 +396,20 @@ int LoMaxAction(SConnection *pCon, SicsInterp *pSics,
(int)ObVal(self->pParam,THRESHOLD),
&intensity))
{
- Tcl_DStringStartSublist(&result);
- sprintf(pNum,"%d", i);
- Tcl_DStringAppendElement(&result,pNum);
- sprintf(pNum,"%d", j);
- Tcl_DStringAppendElement(&result,pNum);
+ if(count != 0)
+ {
+ Tcl_DStringAppend(&result,"@",strlen("@"));
+ }
+ sprintf(pNum,"%d ", i);
+ Tcl_DStringAppend(&result,pNum,strlen(pNum));
+ sprintf(pNum,"%d ", j);
+ Tcl_DStringAppend(&result,pNum,strlen(pNum));
sprintf(pNum,"%d", intensity);
- Tcl_DStringAppendElement(&result,pNum);
- Tcl_DStringEndSublist(&result);
+ Tcl_DStringAppend(&result,pNum,strlen(pNum));
+ count++;
}
}
}
- Tcl_DStringEndSublist(&result);
SCWrite(pCon,Tcl_DStringValue(&result),eValue);
Tcl_DStringFree(&result);
return 1;
@@ -347,27 +440,87 @@ int LoMaxAction(SConnection *pCon, SicsInterp *pSics,
return 0;
}
Tcl_DStringInit(&result);
- Tcl_DStringStartSublist(&result);
iData = GetHistogramPointer(pHM,pCon);
window = (int)ObVal(self->pParam,COGWINDOW);
- iRet = calculateCOG(iData,iDim[0], iDim[1], &i, &j, &intensity,
+ iRet = calculateCOG(iData,iDim[0], iDim[1], &i, &j, &intensity,&count,
window, ObVal(self->pParam,COGCONTOUR));
if(!iRet)
{
SCWrite(pCon,"ERROR: no intensity in data",eError);
return 0;
}
- sprintf(pNum,"%d", i);
- Tcl_DStringAppendElement(&result,pNum);
- sprintf(pNum,"%d", j);
- Tcl_DStringAppendElement(&result,pNum);
- sprintf(pNum,"%d", intensity);
- Tcl_DStringAppendElement(&result,pNum);
- Tcl_DStringEndSublist(&result);
+ sprintf(pNum,"%d ", i);
+ Tcl_DStringAppend(&result,pNum,strlen(pNum));
+ sprintf(pNum,"%d ", j);
+ Tcl_DStringAppend(&result,pNum,strlen(pNum));
+ sprintf(pNum,"%d ", intensity);
+ Tcl_DStringAppend(&result,pNum,strlen(pNum));
+ sprintf(pNum,"%d ", count);
+ Tcl_DStringAppend(&result,pNum,strlen(pNum));
SCWrite(pCon,Tcl_DStringValue(&result),eValue);
Tcl_DStringFree(&result);
return 1;
}
+ else if(strcmp(argv[1],"wellformed") == 0) /* test for wellformedness */
+ {
+ if(argc < 6)
+ {
+ sprintf(pBueffel,
+ "ERROR: insufficient number of arguments to %s.wellformed",
+ argv[0]);
+ SCWrite(pCon,pBueffel,eError);
+ return 0;
+ }
+ if(!checkHM(&pHM, pSics,pCon, argv[2],iDim))
+ {
+ return 0;
+ }
+ if(Tcl_GetInt(pSics->pTcl,argv[3],&i)!= TCL_OK)
+ {
+ sprintf(pBueffel,"ERROR: expected number, got %s",argv[2]);
+ SCWrite(pCon,pBueffel,eError);
+ return 0;
+ }
+ if(Tcl_GetInt(pSics->pTcl,argv[4],&j)!= TCL_OK)
+ {
+ sprintf(pBueffel,"ERROR: expected number, got %s",argv[2]);
+ SCWrite(pCon,pBueffel,eError);
+ return 0;
+ }
+ if(Tcl_GetInt(pSics->pTcl,argv[5],&badMax)!= TCL_OK)
+ {
+ sprintf(pBueffel,"ERROR: expected number, got %s",argv[2]);
+ SCWrite(pCon,pBueffel,eError);
+ return 0;
+ }
+ iData = GetHistogramPointer(pHM,pCon);
+ window = (int)ObVal(self->pParam,COGWINDOW);
+ iRet = wellFormed(iData,iDim[0], iDim[1], i, j,
+ window, ObVal(self->pParam,COGCONTOUR),
+ badMax);
+ sprintf(pBueffel,"%5d", iRet);
+ SCWrite(pCon,pBueffel,eValue);
+ return 1;
+ }
+ else if(strcmp(argv[1],"stat") == 0)
+ {
+ if(argc < 3)
+ {
+ sprintf(pBueffel,"ERROR: insufficient number of arguments to %s.search",
+ argv[0]);
+ SCWrite(pCon,pBueffel,eError);
+ return 0;
+ }
+ if(!checkHM(&pHM, pSics,pCon, argv[2],iDim))
+ {
+ return 0;
+ }
+ iData = GetHistogramPointer(pHM,pCon);
+ calculateStatistics(iData,iDim[0],iDim[1],&average,&maximum);
+ sprintf(pBueffel," %f %f", average, maximum);
+ SCWrite(pCon,pBueffel,eValue);
+ return 1;
+ }
else
{
/* we are handling one of the parameter commands
diff --git a/lomax.h b/lomax.h
index b27c338a..823a1b3a 100644
--- a/lomax.h
+++ b/lomax.h
@@ -27,10 +27,15 @@
int window, int threshold, int steepness,
int *intensity);
int calculateCOG(int *iData, int xsize, int ysize,
- int *i, int *j, int *intensity,
+ int *i, int *j, int *intensity, int *count,
int cogWindow,
float contour);
-
+ void calculateStatistics(int *iData, int xsize, int ysize,
+ float *average, float *maximum);
+ int wellFormed(int *iData, int xsize, int ysize,
+ int x, int y, int window, float contour,
+ int maxBad);
+
diff --git a/lomax.w b/lomax.w
index d6684646..f3fb2595 100644
--- a/lomax.w
+++ b/lomax.w
@@ -46,10 +46,15 @@ ever needed.
int window, int threshold, int steepness,
int *intensity);
int calculateCOG(int *iData, int xsize, int ysize,
- int *i, int *j, int *intensity,
+ int *i, int *j, int *intensity, int *count,
int cogWindow,
float contour);
-
+ void calculateStatistics(int *iData, int xsize, int ysize,
+ float *average, float *maximum);
+ int wellFormed(int *iData, int xsize, int ysize,
+ int x, int y, int window, float contour,
+ int maxBad);
+
@}
testLocalMaxima checks if point i, j is a valid local maximum. It
@@ -63,6 +68,15 @@ calculateCOG calculates the center of gravity for point i, j,. Points within a
calculation. The position of the maximum (i,j) and its intensity are
updated by calculateCOG.
+calculateStatistics finds the average and the maximum value of the data
+in iData.
+
+wellFormed checks a candidate peak position at x, y for well
+formedness. Basically it checks if there are points above
+contour * maximum at the borders of the COG window. If this count is
+larger then maxBad 0 is returned else 1. This should catch powder
+lines and guard against overlapped peaks.
+
@o lomax.h @{
/*-------------------------------------------------------------------------
diff --git a/napi.h b/napi.h
index 84d734c0..20d6244c 100644
--- a/napi.h
+++ b/napi.h
@@ -34,7 +34,6 @@
#ifndef NEXUSAPI
#define NEXUSAPI
-
/* NeXus HDF45 */
#define NEXUS_VERSION "2.1.0." /* major.minor.patch */
@@ -162,6 +161,7 @@ typedef struct {
# define NXgetgroupinfo MANGLE(nxigetgroupinfo)
# define NXinitgroupdir MANGLE(nxiinitgroupdir)
# define NXinitattrdir MANGLE(nxiinitattrdir)
+# define NXsetcache MANGLE(nxisetcache)
/* FORTRAN helpers - for NeXus internal use only */
# define NXfopen MANGLE(nxifopen)
# define NXfclose MANGLE(nxifclose)
@@ -171,7 +171,6 @@ typedef struct {
# define NXfcompress MANGLE(nxifcompress)
# define NXfputattr MANGLE(nxifputattr)
-
/*
* Standard interface
*/
@@ -220,8 +219,10 @@ NX_EXTERNAL NXstatus CALLING_STYLE NXfree(void** data);
*/
NX_EXTERNAL void CALLING_STYLE NXMSetError(void *pData, void (*ErrFunc)(void *pD, char *text));
+/*
+ another special function for setting the default cache size for HDF-5
+*/
+NX_EXTERNAL NXstatus CALLING_STYLE NXsetcache(long newVal);
#endif /*NEXUSAPI*/
-
-
diff --git a/nextrics.c b/nextrics.c
index 1f8d225d..e41236b8 100644
--- a/nextrics.c
+++ b/nextrics.c
@@ -307,7 +307,7 @@ name of hkl object holding crystallographic information
char pBueffel[512];
CounterMode eMode;
HistInt lData[DET1X*DET1Y], i, lVal;
- int32 iVal;
+ int32 iVal,lBeam;
float fVal, fTTheta;
pMotor pMot;
@@ -347,24 +347,28 @@ name of hkl object holding crystallographic information
eMode = GetHistCountMode(self->pHistogram1);
fVal = GetHistPreset(self->pHistogram1);
lVal = GetHistMonitor(self->pHistogram1,1,pCon);
+ lBeam = GetHistMonitor(self->pHistogram1,4,pCon);
}
if(self->pHistogram2 != NULL)
{
eMode = GetHistCountMode(self->pHistogram2);
fVal = GetHistPreset(self->pHistogram2);
lVal = GetHistMonitor(self->pHistogram2,1,pCon);
+ lBeam = GetHistMonitor(self->pHistogram2,4,pCon);
}
if(self->pHistogram3 != NULL)
{
eMode = GetHistCountMode(self->pHistogram3);
fVal = GetHistPreset(self->pHistogram3);
lVal = GetHistMonitor(self->pHistogram3,1,pCon);
+ lBeam = GetHistMonitor(self->pHistogram2,4,pCon);
}
if(self->pCount != NULL)
{
eMode = GetCounterMode(self->pCount);
fVal = GetCounterPreset(self->pCount);
lVal = GetMonitor(self->pCount,1,pCon);
+ lBeam = GetMonitor(self->pCount,4,pCon);
}
if(eMode == eTimer)
{
@@ -378,6 +382,8 @@ name of hkl object holding crystallographic information
NXDputalias(hfil,self->pDict,"framepreset",&fVal);
iVal = (int32)lVal;
NXDputalias(hfil,self->pDict,"framemonitor",&iVal);
+ iVal = (int32)lBeam;
+ NXDputalias(hfil,self->pDict,"sinqmonitor",&iVal);
/* write detector1 histogram */
if(self->pHistogram1 != NULL)
@@ -676,7 +682,7 @@ name of hkl object holding crystallographic information
SCWrite(pCon,"ERROR: failed to write detecor y zero point",eError);
}
pVar = NULL;
- pVar = FindVariable(pServ->pSics,"det1dist");
+ pVar = FindVariable(pServ->pSics,"detdist1");
if(pVar)
{
fVal = pVar->fVal;
@@ -772,7 +778,7 @@ name of hkl object holding crystallographic information
eError);
}
pVar = NULL;
- pVar = FindVariable(pServ->pSics,"det2dist");
+ pVar = FindVariable(pServ->pSics,"detdist2");
if(pVar)
{
fVal = pVar->fVal;
@@ -799,7 +805,7 @@ name of hkl object holding crystallographic information
/*
third detector, but only if present
*/
- if(self->pHistogram2 != NULL)
+ if(self->pHistogram3 != NULL)
{
strcpy(pBueffel,"detector3");
NXDupdate(self->pDict,"dnumber",pBueffel);
@@ -811,7 +817,7 @@ name of hkl object holding crystallographic information
iRet = NXDputalias(hfil,self->pDict,"ddescription",pBueffel);
if(iRet != NX_OK)
{
- SCWrite(pCon,"ERROR: failed to write detector2 description",eError);
+ SCWrite(pCon,"ERROR: failed to write detector3 description",eError);
}
for(i = 0; i < DET3X; i++)
{
@@ -870,7 +876,7 @@ name of hkl object holding crystallographic information
eError);
}
pVar = NULL;
- pVar = FindVariable(pServ->pSics,"det3dist");
+ pVar = FindVariable(pServ->pSics,"detdist3");
if(pVar)
{
fVal = pVar->fVal;
diff --git a/nxamor.c b/nxamor.c
index d9e68c24..e69e99b2 100644
--- a/nxamor.c
+++ b/nxamor.c
@@ -76,6 +76,7 @@ pAmor2T pAmor = NULL;
float fVal;
/* open files */
+ NXsetcache(TOFBLOCK);
iRet = NXopen(file,NXACC_CREATE5,&hfil);
if(iRet != NX_OK)
{
diff --git a/ofac.c b/ofac.c
index a60e87b9..5079f931 100644
--- a/ofac.c
+++ b/ofac.c
@@ -106,6 +106,7 @@
#include "hmcontrol.h"
#include "rs232controller.h"
#include "lomax.h"
+#include "polterwrite.h"
/*----------------------- Server options creation -------------------------*/
static int IFServerOption(SConnection *pCon, SicsInterp *pSics, void *pData,
int argc, char *argv[])
@@ -287,6 +288,7 @@
AddCommand(pInter,"MakeHMControl",MakeHMControl,NULL,NULL);
AddCommand(pInter,"MakeRS232Controller",RS232Factory,NULL,NULL);
AddCommand(pInter,"MakeMaxDetector",LoMaxFactory,NULL,NULL);
+ AddCommand(pInter,"PolterInstall",PolterInstall,NULL,NULL);
}
/*---------------------------------------------------------------------------*/
static void KillIniCommands(SicsInterp *pSics)
@@ -345,6 +347,7 @@
RemoveCommand(pSics,"MakeHMControl");
RemoveCommand(pSics,"MakeRS232Controller");
RemoveCommand(pSics,"MakeMaxDetector");
+ RemoveCommand(pSics,"PolterInstall");
}
diff --git a/peaksearch.tcl b/peaksearch.tcl
index c5b6423b..d33ae393 100644
--- a/peaksearch.tcl
+++ b/peaksearch.tcl
@@ -5,6 +5,8 @@
#---------------------------------------------------------------------------
proc initPeakSearch {} {
+
+ #------- phi Range
VarMake ps.phiStart Float User
ps.phiStart 0
VarMake ps.phiEnd Float User
@@ -12,6 +14,7 @@ proc initPeakSearch {} {
VarMake ps.phiStep Float User
ps.phiStep 3.
+ #-------- chi range
VarMake ps.chiStart Float User
ps.chiStart 0
VarMake ps.chiEnd Float User
@@ -19,6 +22,7 @@ proc initPeakSearch {} {
VarMake ps.chiStep Float User
ps.chiStep 12.
+ #-------- omega range
VarMake ps.omStart Float User
ps.omStart 0
VarMake ps.omEnd Float User
@@ -26,6 +30,7 @@ proc initPeakSearch {} {
VarMake ps.omStep Float User
ps.omStep 3.
+ #------- two theta range
VarMake ps.sttStart Float User
ps.sttStart 5
VarMake ps.sttEnd Float User
@@ -33,6 +38,7 @@ proc initPeakSearch {} {
VarMake ps.sttStep Float User
ps.sttStep 3.
+ #------- maximum finding parameters
VarMake ps.threshold Int User
ps.threshold 30
VarMake ps.steepness Int User
@@ -44,21 +50,51 @@ proc initPeakSearch {} {
VarMake ps.cogcontour Float User
ps.cogcontour .2
+ #-------- counting parameters
VarMake ps.countmode Text User
ps.countmode monitor
-
VarMake ps.preset Float User
ps.preset 1000
+ #-------- final scan counting parameters
+ VarMake ps.scanpreset Float User
+ ps.scanpreset 1000000
+ VarMake ps.scansteps Int User
+ ps.scansteps 24
+
+ #--------- file to which to write the results
+ VarMake ps.listfile Text User
+ ps.listfile peaksearch.dat
+
+ #--------- conversion factors from Pixel to mm
+ VarMake xfactor Float Mugger
+ xfactor 0.715
+ VarMake yfactor Float Mugger
+ yfactor 1.42
+
+ #--------- published functions
Publish ps.phirange User
Publish ps.chirange User
Publish ps.omrange User
Publish ps.sttrange User
Publish ps.countpar User
+ Publish ps.scanpar User
Publish ps.maxpar User
+ Publish ps.list User
+ Publish ps.listpeaks User
+ Publish ps.run User
+ Publish ps.continue User
+ Publish ps.scanlist User
#------- these are for debugging only!
Publish checkomega User
+ Publish optimizedetector User
Publish optimizepeak User
+ Publish printpeak User
+ Publish initsearch User
+ Publish catchdrive User
+ Publish catchdriveval User
+ Publish scandetectorsD User
+ Publish scandetectors User
}
#--------------------------------------------------------------------------
proc ps.phirange args {
@@ -122,6 +158,16 @@ proc ps.countpar args {
[SplitReply [ps.countmode]] [SplitReply [ps.preset]]]
}
#-------------------------------------------------------------------------
+proc ps.scanpar args {
+ if { [llength $args] >= 2 } {
+ ps.scanpreset [lindex $args 0]
+ ps.scansteps [lindex $args 1]
+ }
+ clientput "Peak Search Scan Parameters:"
+ return [format " Count Preset = %12.2f, No. Steps %4d" \
+ [SplitReply [ps.scanpreset]] [SplitReply [ps.scansteps]]]
+}
+#-------------------------------------------------------------------------
proc ps.maxpar args {
if { [llength $args] >= 5 } {
ps.window [lindex $args 0]
@@ -131,28 +177,43 @@ proc ps.maxpar args {
ps.cogcontour [lindex $args 4]
}
clientput "Peak Search Maximum Detection Parameters:"
- set t1 [format " Window = %d, Threshold = %d, Steepness = %d" \
+ set t1 [format " Window = %d, Threshold = %d * average, Steepness = %d" \
[SplitReply [ps.window]] [SplitReply [ps.threshold]] \
[SplitReply [ps.steepness] ]]
set t2 [format " COGWindow = %d, COGcontour = %f " \
[SplitReply [ps.cogwindow]] [SplitReply [ps.cogcontour]]]
return [format "%s\n%s" $t1 $t2]
}
+#-----------------------------------------------------------------------
+proc ps.list {} {
+ clientput [ps.sttrange]
+ clientput [ps.omrange]
+ clientput [ps.chirange]
+ clientput [ps.phirange]
+ clientput [ps.countpar]
+ clientput [ps.scanpar]
+ clientput [ps.maxpar]
+}
#------------------------------------------------------------------------
-proc checknewomega {omega maxIntensity} {
- if {[catch {drive $omega} msg] != 0} {
+proc string2list {txt} {
+ return [split [string trim $txt \{\}]]
+}
+#------------------------------------------------------------------------
+proc checknewomega {hm x y omega maxIntensity} {
+ if {[catch {drive om $omega} msg] != 0} {
error $msg
}
- if {[catch {hmc start [SplitReply [ps.preset]] [SplitReply \
- [ps.countmode]]} msg] != 0} {
+ if {[catch {hmc start [SplitReply [ps.preset]] [string trim [SplitReply \
+ [ps.countmode]]]} msg] != 0} {
error $msg
}
if {[catch {Success} msg] != 0} {
error $msg
}
if { [catch {lomax cog $hm $x $y} result] != 0} {
- error "Failed to calculate COG"
+ error "Failed to calculate COG: $result"
}
+ set result [split $result " "]
if {[lindex $result 2] > $maxIntensity } {
return $result
} else {
@@ -160,19 +221,45 @@ proc checknewomega {omega maxIntensity} {
}
}
#------------------------------------------------------------------------
+proc optimizedetector {hm} {
+ if { [catch {lomax stat $hm} result] != 0} {
+#--- This can be due to the fact that the detector is missing. Sigh ....
+ return
+ }
+ set l2 [split [string trim $result]]
+ lomax threshold [expr [lindex $l2 0] * [SplitReply [ps.threshold]]]
+ set result [lomax search $hm]
+ set oldom [SplitReply [om]]
+ set result [split $result @]
+ for {set i 0} { $i < [llength $result]} {incr i} {
+ if { [catch {drive om $oldom} msg] != 0} {
+ error $msg
+ }
+ set piecks [split [lindex $result $i] " "]
+ set x [lindex $piecks 0]
+ set y [lindex $piecks 1]
+ if { [catch {optimizepeak $hm $x $y} msg] != 0} {
+ clientput [format "Aborted peak at %d %d with %s" $x $y $msg]
+ }
+ }
+}
+#------------------------------------------------------------------------
proc optimizepeak {hm x y} {
if { [catch {lomax cog $hm $x $y} result] != 0} {
- error "Failed to calculate COG"
+ error "Failed to calculate COG: $result"
}
+ set result [split $result " "]
set xMax [lindex $result 0]
- set ymax [lindex $result 1]
+ set yMax [lindex $result 1]
set maxIntensity [lindex $result 2]
set maxOmega [SplitReply [om]]
set startOmega $maxOmega
+ set omSearchStep .1
#--------- move to positive omega until maximum found
while {1} {
- set newOm [expr [SplitReply [om]] + [SplitReply [ps.omStep]]]
- if {[catch {checknewomega $newOm $maxIntensity} result] != 0} {
+ set newOm [expr [SplitReply [om]] + $omSearchStep]
+ if {[catch {checknewomega $hm $xMax $yMax $newOm $maxIntensity} \
+ result] != 0} {
error $result
}
if {$result != 0} {
@@ -189,8 +276,9 @@ proc optimizepeak {hm x y} {
# negative direction
if {$maxOmega == $startOmega} {
while {1} {
- set newOm [expr [SplitReply [om]] - [SplitReply [ps.omStep]]]
- if {[catch {checknewomega $newOm $maxIntensity} result] != 0} {
+ set newOm [expr [SplitReply [om]] - $omSearchStep]
+ if {[catch {checknewomega $hm $xMax $yMax $newOm $maxIntensity} \
+ result] != 0} {
error $result
}
if {$result != 0} {
@@ -204,5 +292,208 @@ proc optimizepeak {hm x y} {
}
}
#----------- print the results we have found
+ printpeak $hm $xMax $yMax $maxOmega $maxIntensity
+#------------ scan the peak for Oksana
+# set scanStart [expr $maxOmega - 0.1*([SplitReply [ps.scansteps]]/2)]
+# if { [catch {tricsscan $scanStart .1 [SplitReply [ps.scansteps]] \
+# [SplitReply [ps.countmode]] [SplitReply [ps.scanpreset]]} msg] } {
+# error $msg
+# }
}
+#---------------------------------------------------------------------------
+proc printpeak {hm x y om intensity} {
+ set phval [SplitReply [phi]]
+ set chval [SplitReply [chi]]
+ set gamOffset .0
+ switch $hm {
+ hm1 {
+ set tilt [SplitReply [dg1]]
+ set gamOffset .0
+ }
+ hm2 {
+ set tilt [SplitReply [dg2]]
+ set gamOffset 0.
+ }
+ hm3 {
+ set tilt [SplitReply [dg3]]
+ set gamOffset 45.
+ }
+ default {error "Invalid hm requested in printpeak"}
+ }
+ set sttval [expr [SplitReply [stt]] + $gamOffset]
+ set zero [SplitReply [$hm configure dim0]]
+ set xval [expr $x * [SplitReply [xfactor]]]
+ set zero [SplitReply [$hm configure dim1]]
+ set yval [expr $y * [SplitReply [yfactor]]]
+ set line [format "%7.2f%7.2f%7.2f%7.2f%7.2f%7.2f%7.2f%10d" \
+ $xval $yval $sttval $om $chval $phval $tilt $intensity]
+ clientput "Found Peak at:"
+ clientput $line
+ set f [open [string trim [SplitReply [ps.listfile]]] a+]
+ puts $f $line
+ close $f
+}
+#--------------------------------------------------------------------------
+proc ps.listpeaks {} {
+ clientput "Peakse found so far: "
+ clientput " X Y STT OM CHI PHI TILT INTENSITY"
+ set f [open [string trim [SplitReply [ps.listfile]]] r]
+ while {[gets $f line] > 0} {
+ clientput [format "%s" $line]
+ }
+ close $f
+}
+#-------------------------------------------------------------------------
+proc initsearch {filename} {
+#----- stow away filename and empty it
+ ps.listfile $filename
+ set f [open $filename w]
+ close $f
+#----- tell lomax its parameters
+ lomax threshold [SplitReply [ps.threshold]]
+ lomax steepness [SplitReply [ps.steepness]]
+ lomax window [SplitReply [ps.window]]
+ lomax cogwindow [SplitReply [ps.cogwindow]]
+ lomax cogcontour [SplitReply [ps.cogcontour]]
+#----- drive to start
+ if { [catch {drive stt [SplitReply [ps.sttStart]] \
+ om [SplitReply [ps.omStart]] \
+ chi [SplitReply [ps.chiStart]] \
+ phi [SplitReply [ps.phiStart]] } msg] != 0} {
+ error $msg
+ }
+}
+#------------------------------------------------------------------------
+# This code means: ignore any errors except interrupts when searching
+proc scandetectors {} {
+# set names [list hm1 hm2 hm3]
+ set names [list hm2 hm3]
+ if {[catch {hmc start [SplitReply [ps.preset]] [string trim [SplitReply \
+ [ps.countmode]]]} msg] != 0} {
+ if{[string compare [getint] continue] != 0} {
+ error $msg
+ } else {
+ clientput [format "Ignoring: %s" $msg]
+ }
+ }
+ if {[catch {Success} msg] != 0} {
+ if{[string compare [getint] continue] != 0} {
+ error $msg
+ } else {
+ clientput [format "Ignoring: %s" $msg]
+ }
+ }
+ for {set i 0} { $i < [llength $names]} {incr i} {
+ set ret [catch {optimizedetector [lindex $names $i]} msg]
+ if { $ret != 0} {
+ if {[string compare [getint] continue] != 0} {
+ error $msg
+ } else {
+ clientput [format "Ignoring problem: %s" $msg]
+ }
+ }
+ }
+}
+#------------------------------------------------------------------------
+# loop debugging
+
+proc scandetectorsD {} {
+ clientput [format "stt = %6.2f, om = %6.2f, chi = %6.2f, phi = %6.2f" \
+ [SplitReply [stt]] [SplitReply [om]] \
+ [SplitReply [chi]] [SplitReply [phi]]]
+ wait 1
+}
+#-----------------------------------------------------------------------
+proc catchdrive { mot step} {
+ set ret [catch {drive $mot [expr [SplitReply [$mot]] + $step]} msg]
+ if {$ret != 0} {
+ if {[string compare [getint] continue] != 0} {
+ error $msg
+ } else {
+ clientput [format "Ignoring: %s" $msg]
+ }
+ }
+}
+#-----------------------------------------------------------------------
+proc catchdriveval { mot val} {
+ set ret [catch {drive $mot $val} msg]
+ if {$ret != 0} {
+ if {[string compare [getint] continue] != 0} {
+ error $msg
+ } else {
+ clientput [format "Ignoring: %s" $msg]
+ }
+ }
+}
+#------------------------------------------------------------------------
+# The actual loop. It is written in a way which allows for the continuation
+# of a search
+
+proc searchloop { } {
+ set sttStep [SplitReply [ps.sttStep]]
+ set sttEnd [SplitReply [ps.sttEnd]]
+ set chiStep [SplitReply [ps.chiStep]]
+ set chiEnd [SplitReply [ps.chiEnd]]
+ set phiStep [SplitReply [ps.phiStep]]
+ set phiEnd [SplitReply [ps.phiEnd]]
+ set omStep [SplitReply [ps.omStep]]
+ set omEnd [SplitReply [ps.omEnd]]
+ while {[SplitReply [stt]] + $sttStep <= $sttEnd} {
+ while {[SplitReply [chi]] + $chiStep <= $chiEnd} {
+ while {[SplitReply [om]] + $omStep <= $omEnd} {
+ while {[SplitReply [phi]] + $phiStep <= $phiEnd} {
+ scandetectors
+ catchdrive phi $phiStep
+ }
+ catchdrive om $omStep
+ catchdriveval phi [SplitReply [ps.phiStart]]
+ }
+ catchdrive chi $chiStep
+ catchdriveval om [SplitReply [ps.omStart]]
+ }
+ catchdrive stt $sttStep
+ catchdriveval chi [SplitReply [ps.chiStart]]
+ }
+ return "Peak Search finished normally"
+}
+#---------------------------------------------------------------------------
+proc ps.run {filename} {
+ initsearch $filename
+ searchloop
+}
+#-------------------------------------------------------------------------
+proc ps.continue {} {
+ searchloop
+}
+#------------------------------------------------------------------------
+proc ps.scanlist {} {
+ if { [catch {set f [open [string trim [SplitReply [ps.listfile]]] "r"]} \
+ msg ] != 0} {
+ error $msg
+ }
+ while { [gets $f line] > 0} {
+ set n [stscan $line "%f %f %f %f %f %f" x y stt om chi phi]
+ if {$n < 6} {
+ clientput [format "Skipping invalid line: %s" line]
+ continue
+ }
+ if { [catch {drive stt $stt om $om chi $chi phi $phi} msg] != 0 } {
+ clientput $msg
+ if {[string compare [getint] continue] != 0} {
+ error "ps.scanlist interupted"
+ }
+ }
+ set scanStart [expr $om - 0.1*([SplitReply [ps.scansteps]]/2)]
+ if { [catch {tricsscan $scanStart .1 [SplitReply [ps.scansteps]] \
+ [SplitReply [ps.countmode]] [SplitReply [ps.scanpreset]]} msg] \
+ != 0 } {
+ clientput $msg
+ if {[string compare [getint] continue] != 0} {
+ error "ps.scanlist interupted"
+ }
+ }
+ }
+ close $f
+ return "Scanning list finished"
+}
\ No newline at end of file
diff --git a/polterwrite.c b/polterwrite.c
new file mode 100644
index 00000000..d4861859
--- /dev/null
+++ b/polterwrite.c
@@ -0,0 +1,375 @@
+/*--------------------------------------------------------------------------
+ P O L T E R W R I T E
+
+ fowrite is an object for writing POLTERDI data files.
+
+ copyright: see copyright.h
+
+ Uwe Filges, November 2001
+----------------------------------------------------------------------------*/
+#include
+#include
+#include "fortify.h"
+#include "sics.h"
+#include "counter.h"
+#include "HistMem.h"
+#include "nxdict.h"
+#include "nxutil.h"
+#include "motor.h"
+#include "sicsvar.h"
+#include "polterwrite.h"
+
+/*
+ diaphragm1 - chopper
+*/
+#define DIA1DIST 3000
+/*
+ diaphragm2 - diaphragm1
+*/
+#define DIA2DIST 500
+
+/*
+ diaphragm2 - sample position
+*/
+#define SADIST 500
+
+/*
+ histogram memory name
+*/
+#define HM "hm"
+/*
+ detector distance - sample
+*/
+#define DETDIST 700
+/*
+ no detectors, change in poldi.dic as well
+*/
+#define NODET 400
+/*
+ theta spacing
+*/
+#define THSPACE .3
+
+
+
+typedef struct {
+ pObjectDescriptor pDes;
+ char *dictfile;
+}Polterdi, *pPolterdi;
+
+/*----------------------------------------------------------------------*/
+static void KillPolterdi(void *pData){
+ pPolterdi self = (pPolterdi)pData;
+ if(!self)
+ return;
+
+ if(self->pDes){
+ DeleteDescriptor(self->pDes);
+ }
+ if(self->dictfile){
+ free(self->dictfile);
+ }
+ free(self);
+}
+
+/*-----------------------------------------------------------------------*/
+int PolterInstall(SConnection *pCon, SicsInterp *pSics,
+ void *pData, int argc, char *argv[])
+{
+ pPolterdi pNew = NULL;
+
+ if(argc < 2){
+ SCWrite(pCon,"ERROR: insufficient number of arguments to PolterInstall",
+ eError);
+ return 0;
+ }
+
+ pNew = (pPolterdi)malloc(sizeof(Polterdi));
+ if(!pNew){
+ SCWrite(pCon,"ERROR: out of memory in PolterInstall",
+ eError);
+ return 0;
+ }
+ memset(pNew,0,sizeof(Polterdi));
+ pNew->pDes = CreateDescriptor("PolterdiWrite");
+ pNew->dictfile = strdup(argv[1]);
+ if(!pNew->pDes || !pNew->dictfile){
+ SCWrite(pCon,"ERROR: out of memory in PolterInstall",
+ eError);
+ return 0;
+ }
+ return AddCommand(pSics,"storedata",PolterAction,KillPolterdi,pNew);
+}
+/*-------------------------------------------------------------------*/
+static void writePolterdiGlobal(NXhandle hfil, NXdict hdict, SConnection *pCon,
+ SicsInterp *pSics){
+ char pBueffel[512];
+
+ SNXSPutVariable(pSics,pCon,hfil,hdict,"etitle","title");
+ SNXFormatTime(pBueffel,511);
+
+ NXDputalias(hfil,hdict,"estart",pBueffel);
+ SNXSPutVariable(pSics,pCon,hfil,hdict,"iname","instrument");
+ sprintf(pBueffel,"%d",strlen("SINQ, PSI, Switzerland"));
+ NXDupdate(hdict,"strdim",pBueffel);
+ NXDputalias(hfil,hdict,"sname","SINQ, PSI, Switzerland");
+ sprintf(pBueffel,"%d",strlen("continous spallation source"));
+ NXDupdate(hdict,"strdim",pBueffel);
+ NXDputalias(hfil,hdict,"stype","continous spallation source");
+ NXDupdate(hdict,"strdim","132");
+}
+/*-------------------------------------------------------------------*/
+static void writeChopper(NXhandle hfil, NXdict hdict, SConnection *pCon,
+ SicsInterp *pSics){
+ SNXSPutVariable(pSics,pCon,hfil,hdict,"cname","choppername");
+ SNXSPutDrivable(pSics,pCon,hfil,hdict,"chopperspeed","crot");
+}
+/*-------------------------------------------------------------------*/
+static void writeDiaphragm1(NXhandle hfil, NXdict hdict, SConnection *pCon,
+ SicsInterp *pSics){
+ float dist;
+
+ SNXSPutMotor(pSics,pCon,hfil,hdict,"dia1x","d1m");
+ SNXSPutMotorNull(pSics,pCon,hfil,hdict,"dia1x0","d1m");
+ SNXSPutMotor(pSics,pCon,hfil,hdict,"dia1y","d1o");
+ SNXSPutMotorNull(pSics,pCon,hfil,hdict,"dia1y0","d1o");
+ dist = (float)DIA1DIST;
+ NXDputalias(hfil,hdict,"dia1dist",&dist);
+}
+
+/*------------------------------------------------------------------*/
+
+static void writeDiaphragm2(NXhandle hfil, NXdict hdict, SConnection *pCon,
+ SicsInterp *pSics){
+ float dist2;
+
+ SNXSPutMotor(pSics, pCon, hfil, hdict, "dia2x_plus", "d2m");
+ SNXSPutMotorNull(pSics, pCon, hfil, hdict, "dia2xplus0", "d2m");
+ SNXSPutMotor(pSics, pCon, hfil, hdict, "dia2x_minus", "d2u");
+ SNXSPutMotorNull(pSics, pCon, hfil, hdict, "dia2xminus0", "d2u");
+ SNXSPutMotor(pSics, pCon, hfil, hdict, "dia2z_plus", "d2o");
+ SNXSPutMotorNull(pSics, pCon, hfil, hdict, "dia2zplus0", "d2o");
+ SNXSPutMotor(pSics, pCon, hfil, hdict, "dia2z_minus", "d2l");
+ SNXSPutMotorNull(pSics, pCon, hfil, hdict, "dia2zminus0", "d2l");
+
+ dist2 = (float) DIA2DIST;
+
+}
+/*--------------------------------------------------------------------*/
+static void writeSample(NXhandle hfil, NXdict hdict, SConnection *pCon,
+ SicsInterp *pSics)
+{
+ float sadist = (float)SADIST;
+
+ NXDputalias(hfil, hdict, "sdist", &sadist);
+ SNXSPutVariable(pSics,pCon,hfil,hdict,"saname","sample");
+ SNXSPutVariable(pSics,pCon,hfil,hdict,"senvir","environment");
+ SNXSPutEVVar(hfil,hdict,"temperature",pCon,"stemp","stddev");
+
+ SNXSPutMotor(pSics, pCon, hfil, hdict, "srsu", "rsu");
+ SNXSPutMotorNull(pSics, pCon, hfil, hdict, "srsu0", "rsu");
+ SNXSPutMotor(pSics, pCon, hfil, hdict, "srsl", "rsl");
+ SNXSPutMotorNull(pSics, pCon, hfil, hdict, "srsl0", "rsl");
+ SNXSPutMotor(pSics, pCon, hfil, hdict, "srsa", "rsa");
+ SNXSPutMotorNull(pSics, pCon, hfil, hdict, "srsa0", "rsa");
+ SNXSPutMotor(pSics, pCon, hfil, hdict, "sshu", "shu");
+ SNXSPutMotorNull(pSics, pCon, hfil, hdict, "sshu0", "shu");
+ SNXSPutMotor(pSics, pCon, hfil, hdict, "sshl", "shl");
+ SNXSPutMotorNull(pSics, pCon, hfil, hdict, "sshl0", "shl");
+ SNXSPutMotor(pSics, pCon, hfil, hdict, "ssgu", "sgu");
+ SNXSPutMotorNull(pSics, pCon, hfil, hdict, "ssgu0", "sgu");
+ SNXSPutMotor(pSics, pCon, hfil, hdict, "ssgl", "sgl");
+ SNXSPutMotorNull(pSics, pCon, hfil, hdict, "ssgl0", "sgl");
+ SNXSPutMotor(pSics, pCon, hfil, hdict, "ssv", "sv");
+ SNXSPutMotorNull(pSics, pCon, hfil, hdict, "ssv0", "sv");
+}
+/*-----------------------------------------------------------------------*/
+static void writeDetector(NXhandle hfil, NXdict hdict, SConnection *pCon,
+ SicsInterp *pSics)
+{
+ const float *fTime = NULL;
+ CounterMode eMode;
+ pHistMem pHist= NULL;
+ float *fTime2 = NULL, fTheta[NODET], fVal;
+ int iLength, i;
+ char pBueffel[80];
+ pMotor sa;
+ HistInt *lData = NULL;
+ long lVal;
+
+ pHist = (pHistMem)FindCommandData(pSics,HM,"HistMem");
+ if(!pHist)
+ {
+ SCWrite(pCon,"ERROR: Histogram memory NOT found",eError);
+ return;
+ }
+
+ /*
+ write time binning
+ */
+ fTime = GetHistTimeBin(pHist,&iLength);
+ fTime2 = (float *)malloc(iLength*sizeof(float));
+ if(fTime2)
+ {
+ for(i = 0;i < iLength; i++)
+ {
+ fTime2[i] = fTime[i]/10.;
+ }
+ sprintf(pBueffel,"%d",iLength);
+ NXDupdate(hdict,"timebin",pBueffel);
+ NXDputalias(hfil,hdict,"dtime",fTime2);
+ free(fTime2);
+ }
+ else
+ {
+ SCWrite(pCon,"ERROR: out of memory while writing time binning",eError);
+ }
+
+ /*
+ do theta
+ */
+ fVal = -999.;
+ sa = FindMotor(pSics,"sa");
+ if(sa)
+ {
+ MotorGetSoftPosition(sa,pCon,&fVal);
+ }
+ else
+ {
+ SCWrite(pCon,"ERROR: failed to locate two theta motor",eError);
+ }
+ for(i = 0; i < NODET; i++)
+ {
+ fTheta[i] = fVal + i * THSPACE;
+ }
+ NXDputalias(hfil,hdict,"dtheta",fTheta);
+
+ fVal = (float)DETDIST;
+ NXDputalias(hfil,hdict,"ddist",&fVal);
+
+ /*
+ write Histogram
+ */
+ lData = GetHistogramPointer(pHist,pCon);
+ if(lData)
+ {
+ NXDputalias(hfil,hdict,"dcounts",lData);
+ }
+ else
+ {
+ SCWrite(pCon,"ERROR: failed to get histogram data",eError);
+ }
+
+ /*
+ write counting data
+ */
+ fVal = GetHistCountTime(pHist,pCon);
+ NXDputalias(hfil,hdict,"cntime",&fVal);
+ lVal = GetHistMonitor(pHist,1,pCon);
+ NXDputalias(hfil,hdict,"cnmon1",&lVal);
+ eMode = GetHistCountMode(pHist);
+ if(eMode == eTimer)
+ {
+ strcpy(pBueffel,"timer");
+ }
+ else
+ {
+ strcpy(pBueffel,"monitor");
+ }
+ NXDputalias(hfil,hdict,"cnmode",pBueffel);
+ fVal = GetHistPreset(pHist);
+ NXDputalias(hfil,hdict,"cnpreset",&fVal);
+}
+/*---------------------------------------------------------------------*/
+static void makeLinks(NXhandle hfil, NXdict hdict, SConnection *pCon,
+ SicsInterp *pSics)
+{
+ NXDaliaslink(hfil,hdict,"dana","dcounts");
+ NXDaliaslink(hfil,hdict,"dana","dtheta");
+ NXDaliaslink(hfil,hdict,"dana","dtime");
+}
+/*--------------------------------------------------------------------*/
+static int StoreData(SConnection *pCon, SicsInterp *pSics, pPolterdi self)
+{
+ char pBueffel[256], *pFile = NULL;
+ NXhandle hfil = NULL;
+ NXdict hdict= NULL;
+ int status;
+
+ /* create filename */
+ pFile = SNXMakeFileName(pSics,pCon);
+ if(!pFile)
+ {
+ SCWrite(pCon,"ERROR: Extra severe: failed to create data file name",
+ eError);
+ return 0;
+ }
+
+ /* create a Nexus file */
+ NXopen(pFile,NXACC_CREATE,&hfil);
+ if(!hfil)
+ {
+ SCWrite(pCon,"ERROR: cannot create data file ",eError);
+ return 0;
+ }
+
+ /* tell Uwe User what we are doing */
+ sprintf(pBueffel,"Writing %s ......",pFile);
+ SCWrite(pCon,pBueffel,eWarning);
+
+ /* write globals */
+ SNXSPutGlobals(hfil,pFile,"POLTERDI",pCon);
+ free(pFile);
+ pFile = NULL;
+
+ /* open nxdict */
+ status = NXDinitfromfile(self->dictfile,&hdict);
+ if(status != NX_OK)
+ {
+ sprintf(pBueffel,"ERROR: failed to open dictionary file %s",
+ self->dictfile);
+ SCWrite(pCon,pBueffel,eError);
+ SCWrite(pCon,"ERROR: Aborting data file writing",eError);
+ SCWrite(pCon,"ERROR: This is a SERIOUS problem!",eError);
+ SCWrite(pCon,"ERROR: DATA NOT WRITTEN",eError);
+ NXclose(&hfil);
+ return 0;
+ }
+
+
+ writePolterdiGlobal(hfil,hdict,pCon,pSics);
+
+ writeChopper(hfil,hdict,pCon,pSics);
+
+ writeDiaphragm1(hfil,hdict,pCon,pSics);
+
+ writeDiaphragm2(hfil,hdict,pCon,pSics);
+
+ writeSample(hfil,hdict,pCon,pSics);
+
+ writeDetector(hfil,hdict,pCon,pSics);
+
+ makeLinks(hfil,hdict,pCon,pSics);
+
+ /* close everything */
+ NXclose(&hfil);
+ NXDclose(hdict,NULL);
+
+ SCSendOK(pCon);
+ return 1;
+}
+/*----------------------------------------------------------------------*/
+int PolterAction(SConnection *pCon, SicsInterp *pSics,
+ void *pData, int argc, char *argv[])
+{
+ pPolterdi self = (pPolterdi)pData;
+
+ assert(self);
+ assert(pCon);
+ assert(pSics);
+
+ return StoreData(pCon,pSics,self);
+}
+
+
+
diff --git a/polterwrite.h b/polterwrite.h
new file mode 100644
index 00000000..68db1bb1
--- /dev/null
+++ b/polterwrite.h
@@ -0,0 +1,21 @@
+/*--------------------------------------------------------------------------
+ P O L T E R W R I T E
+
+ fowrite is an object for writing POLTERDI data files.
+
+ copyright: see copyright.h
+
+ Uwe Filges, November 2001
+----------------------------------------------------------------------------*/
+#ifndef POLTERDIWRITE
+#define POLTERDIWRITE
+
+ int PolterInstall(SConnection *pCon, SicsInterp *pSics,
+ void *pData, int argc, char *argv[]);
+
+
+ int PolterAction(SConnection *pCon, SicsInterp *pSics,
+ void *pData, int argc, char *argv[]);
+
+#endif
+
diff --git a/sanscook.c b/sanscook.c
index d4ab4ea2..ba13d1a9 100644
--- a/sanscook.c
+++ b/sanscook.c
@@ -4,7 +4,7 @@
This is a controller driver for a heaiting device developed especially
for use with SANS at SINQ by somebody at the Some God Forsaken University
- (SGFU). As this device somes with two motors in addition to the heater
+ (SGFU). As this device comes with two motors in addition to the heater
the general controller mechanism as described in choco.tex is used.
copyright: see copyright.h
diff --git a/sicsstat.tcl b/sicsstat.tcl
index 499edb7f..29a305d4 100644
--- a/sicsstat.tcl
+++ b/sicsstat.tcl
@@ -1,343 +1,8 @@
-scaninfo 7,en,-0.300000,0.100000
-scaninfo setAccess 0
-sicsdatapath /data/koenneck/src/sics/tmp/
-sicsdatapath setAccess 1
-arx2 4.290000
-arx2 setAccess 1
-arx1 0.150000
-arx1 setAccess 1
-mrx2 10.420000
-mrx2 setAccess 1
-mrx1 0.280000
-mrx1 setAccess 1
-lpa 0
-lpa setAccess 2
-dt 0.000000
-dt setAccess 2
-dqm 0.000000
-dqm setAccess 2
-qm 0.000000
-qm setAccess 2
-etaa 0.000000
-etaa setAccess 2
-etas 0.000000
-etas setAccess 2
-etam 0.000000
-etam setAccess 2
-wav 0.000000
-wav setAccess 2
-den 0.100000
-den setAccess 2
-dql 0.000000
-dql setAccess 2
-dqk 0.000000
-dqk setAccess 2
-dqh 0.000000
-dqh setAccess 2
-dkf 0.000000
-dkf setAccess 2
-def 0.500000
-def setAccess 2
-dki 0.000000
-dki setAccess 2
-dei 0.500000
-dei setAccess 2
-dagl 0.000000
-dagl setAccess 2
-dsgu 0.000000
-dsgu setAccess 2
-dsgl 0.000000
-dsgl setAccess 2
-dmgl 0.000000
-dmgl setAccess 2
-datu 0.000000
-datu setAccess 2
-datl 0.000000
-datl setAccess 2
-dstu 0.000000
-dstu setAccess 2
-dstl 0.000000
-dstl setAccess 2
-dmtu 0.000000
-dmtu setAccess 2
-dmtl 0.000000
-dmtl setAccess 2
-dach 0.000000
-dach setAccess 2
-dsro 0.000000
-dsro setAccess 2
-dmcv 0.000000
-dmcv setAccess 2
-da6 0.000000
-da6 setAccess 2
-da5 0.000000
-da5 setAccess 2
-da4 0.100000
-da4 setAccess 2
-da3 0.000000
-da3 setAccess 2
-da2 0.000000
-da2 setAccess 2
-da1 0.000000
-da1 setAccess 2
-bet4 0.000000
-bet4 setAccess 2
-bet3 0.000000
-bet3 setAccess 2
-bet2 0.000000
-bet2 setAccess 2
-bet1 0.000000
-bet1 setAccess 2
-alf4 0.000000
-alf4 setAccess 2
-alf3 0.000000
-alf3 setAccess 2
-alf2 3.000000
-alf2 setAccess 2
-alf1 11.000000
-alf1 setAccess 2
-local Mordahl Schlawadini
-local setAccess 2
-output a4
-output setAccess 2
-lastcommand sc en 0 den .1 np 7 ti 2
-lastcommand setAccess 2
-user Billy Looser
-user setAccess 2
-title SimSulfid Test
-title setAccess 2
-f2 0
-f2 setAccess 2
-f1 0
-f1 setAccess 2
-swunit 0
-swunit setAccess 2
-hz 0.000000
-hz setAccess 2
-hy 0.000000
-hy setAccess 2
-hx 0.000000
-hx setAccess 2
-helm 0.000000
-helm setAccess 2
-if2h 0.000000
-if2h setAccess 2
-if1h 0.000000
-if1h setAccess 2
-if2v 0.000000
-if2v setAccess 2
-if1v 0.000000
-if1v setAccess 2
-mn 700
-mn setAccess 2
-ti 2.000000
-ti setAccess 2
-np 7
-np setAccess 2
-fx 2
-fx setAccess 2
-sa -1
-sa setAccess 2
-ss 1
-ss setAccess 2
-sm 1
-sm setAccess 2
-da 3.354000
-da setAccess 1
-dm 3.354000
-dm setAccess 1
-en 0.300000
-en setAccess 2
-ql 0.000000
-ql setAccess 2
-qk 0.000000
-qk setAccess 2
-qh 2.000000
-qh setAccess 2
-kf 1.964944
-kf setAccess 2
-ef 8.000000
-ef setAccess 2
-ki 2.001447
-ki setAccess 2
-ei 8.300000
-ei setAccess 2
-bz 1.000000
-bz setAccess 2
-by 0.000000
-by setAccess 2
-bx 0.000000
-bx setAccess 2
-az 0.000000
-az setAccess 2
-ay 0.000000
-ay setAccess 2
-ax 1.000000
-ax setAccess 2
-cc 90.000000
-cc setAccess 2
-bb 90.000000
-bb setAccess 2
-aa 90.000000
-aa setAccess 2
-cs 5.000000
-cs setAccess 2
-bs 5.000000
-bs setAccess 2
-as 5.000000
-as setAccess 2
-# Counter counter
-counter SetPreset 2.000000
-counter SetMode Timer
-# Motor agl
-agl SoftZero -0.490000
-agl SoftLowerLim -9.510000
-agl SoftUpperLim 10.490000
-agl Fixed -1.000000
-agl sign 1.000000
-agl InterruptMode 0.000000
-agl AccessCode 2.000000
-# Motor sgu
-sgu SoftZero 0.000000
-sgu SoftLowerLim -16.000000
-sgu SoftUpperLim 16.000000
-sgu Fixed -1.000000
-sgu sign 1.000000
-sgu InterruptMode 0.000000
-sgu AccessCode 2.000000
-# Motor sgl
-sgl SoftZero 1.550000
-sgl SoftLowerLim -17.549999
-sgl SoftUpperLim 14.450000
-sgl Fixed -1.000000
-sgl sign 1.000000
-sgl InterruptMode 0.000000
-sgl AccessCode 2.000000
-# Motor mgl
-mgl SoftZero 1.300000
-mgl SoftLowerLim -11.300000
-mgl SoftUpperLim 8.700000
-mgl Fixed -1.000000
-mgl sign 1.000000
-mgl InterruptMode 0.000000
-mgl AccessCode 2.000000
-# Motor atu
-atu SoftZero -0.880000
-atu SoftLowerLim -16.120001
-atu SoftUpperLim 17.759998
-atu Fixed -1.000000
-atu sign 1.000000
-atu InterruptMode 0.000000
-atu AccessCode 2.000000
-# Motor atl
-atl SoftZero 0.000000
-atl SoftLowerLim -17.000000
-atl SoftUpperLim 17.000000
-atl Fixed -1.000000
-atl sign 1.000000
-atl InterruptMode 0.000000
-atl AccessCode 2.000000
-# Motor stu
-stu SoftZero 0.000000
-stu SoftLowerLim -30.000000
-stu SoftUpperLim 30.000000
-stu Fixed -1.000000
-stu sign 1.000000
-stu InterruptMode 0.000000
-stu AccessCode 2.000000
-# Motor stl
-stl SoftZero 0.000000
-stl SoftLowerLim -30.000000
-stl SoftUpperLim 30.000000
-stl Fixed -1.000000
-stl sign 1.000000
-stl InterruptMode 0.000000
-stl AccessCode 2.000000
-# Motor mtu
-mtu SoftZero 2.850000
-mtu SoftLowerLim -19.850000
-mtu SoftUpperLim 14.150000
-mtu Fixed -1.000000
-mtu sign 1.000000
-mtu InterruptMode 0.000000
-mtu AccessCode 2.000000
-# Motor mtl
-mtl SoftZero 0.000000
-mtl SoftLowerLim -17.000000
-mtl SoftUpperLim 17.000000
-mtl Fixed -1.000000
-mtl sign 1.000000
-mtl InterruptMode 0.000000
-mtl AccessCode 2.000000
-# Motor ach
-ach SoftZero 0.000000
-ach SoftLowerLim -0.500000
-ach SoftUpperLim 11.500000
-ach Fixed -1.000000
-ach sign 1.000000
-ach InterruptMode 0.000000
-ach AccessCode 2.000000
-# Motor sro
-sro SoftZero 0.000000
-sro SoftLowerLim 0.000000
-sro SoftUpperLim 351.000000
-sro Fixed -1.000000
-sro sign 1.000000
-sro InterruptMode 0.000000
-sro AccessCode 2.000000
-# Motor mcv
-mcv SoftZero 0.000000
-mcv SoftLowerLim -9.000000
-mcv SoftUpperLim 124.000000
-mcv Fixed -1.000000
-mcv sign 1.000000
-mcv InterruptMode 0.000000
-mcv AccessCode 2.000000
-# Motor a6
-a6 SoftZero 0.040000
-a6 SoftLowerLim -116.040001
-a6 SoftUpperLim 165.960007
-a6 Fixed -1.000000
-a6 sign 1.000000
-a6 InterruptMode 0.000000
-a6 AccessCode 2.000000
-# Motor a5
-a5 SoftZero 176.479996
-a5 SoftLowerLim -200.000000
-a5 SoftUpperLim 380.000000
-a5 Fixed -1.000000
-a5 sign 1.000000
-a5 InterruptMode 0.000000
-a5 AccessCode 2.000000
-# Motor a4
-a4 SoftZero -0.710000
-a4 SoftLowerLim -134.389999
-a4 SoftUpperLim 124.110001
-a4 Fixed -1.000000
-a4 sign 1.000000
-a4 InterruptMode 0.000000
-a4 AccessCode 2.000000
-# Motor a3
-a3 SoftZero 11.540000
-a3 SoftLowerLim -188.839996
-a3 SoftUpperLim 165.760010
-a3 Fixed -1.000000
-a3 sign 1.000000
-a3 InterruptMode 0.000000
-a3 AccessCode 2.000000
-# Motor a2
-a2 SoftZero -0.010000
-a2 SoftLowerLim 33.109997
-a2 SoftUpperLim 120.010002
-a2 Fixed -1.000000
-a2 sign 1.000000
-a2 InterruptMode 0.000000
-a2 AccessCode 2.000000
-# Motor a1
-a1 SoftZero 0.590000
-a1 SoftLowerLim -0.590000
-a1 SoftUpperLim 110.410004
-a1 Fixed -1.000000
-a1 sign 1.000000
-a1 InterruptMode 0.000000
-a1 AccessCode 2.000000
+# Motor gurke
+gurke SoftZero 0.000000
+gurke SoftLowerLim -20.000000
+gurke SoftUpperLim 20.000000
+gurke Fixed -1.000000
+gurke sign 1.000000
+gurke InterruptMode 0.000000
+gurke AccessCode 2.000000
diff --git a/sicsstatus.tcl b/sicsstatus.tcl
index 8d6f7ef1..4cf57677 100644
--- a/sicsstatus.tcl
+++ b/sicsstatus.tcl
@@ -1,50 +1,33 @@
-ps.preset 20.000000
-ps.preset setAccess 2
-ps.countmode timer
-ps.countmode setAccess 2
-ps.cogcontour 0.200000
-ps.cogcontour setAccess 2
-ps.cogwindow 60
-ps.cogwindow setAccess 2
-ps.window 7
-ps.window setAccess 2
-ps.steepness 2
-ps.steepness setAccess 2
-ps.threshold 50
-ps.threshold setAccess 2
-ps.sttstep 5.000000
-ps.sttstep setAccess 2
-ps.sttend 50.000000
-ps.sttend setAccess 2
-ps.sttstart 10.000000
-ps.sttstart setAccess 2
-ps.omstep 5.000000
-ps.omstep setAccess 2
-ps.omend 20.000000
-ps.omend setAccess 2
-ps.omstart 10.000000
-ps.omstart setAccess 2
-ps.chistep 10.000000
-ps.chistep setAccess 2
-ps.chiend 150.000000
-ps.chiend setAccess 2
-ps.chistart 10.000000
-ps.chistart setAccess 2
-ps.phistep 5.000000
-ps.phistep setAccess 2
-ps.phiend 120.000000
-ps.phiend setAccess 2
-ps.phistart 10.000000
-ps.phistart setAccess 2
-hm3 CountMode timer
-hm3 preset 10.000000
-hm2 CountMode timer
-hm2 preset 5.000000
-hm1 CountMode timer
-hm1 preset 5.000000
+a5l.length 80.000000
+flightpathlength 0.000000
+flightpathlength setAccess 1
+flightpath 0.000000
+flightpath setAccess 1
+delay 2500.000000
+delay setAccess 1
+hm CountMode timer
+hm preset 100.000000
+hm genbin 120.000000 35.000000 512
+hm init
+datafile focus-1001848.hdf
+datafile setAccess 3
+hm2 CountMode monitor
+hm2 preset 2.000000
+hm1 CountMode monitor
+hm1 preset 2.000000
+dbfile UNKNOWN
+dbfile setAccess 2
+# Motor th
+th SoftZero 0.000000
+th SoftLowerLim 4.000000
+th SoftUpperLim 113.000000
+th Fixed -1.000000
+th sign 1.000000
+th InterruptMode 0.000000
+th AccessCode 2.000000
#Crystallographic Settings
hkl lambda 1.179000
-hkl setub 0.016169 0.011969 0.063195 -0.000545 0.083377 -0.009117 -0.162051 0.000945 0.006312
+hkl setub 0.004076 -0.080526 -0.018163 -0.008113 -0.023908 0.061299 -0.161515 -0.000831 -0.003537
det3dist 300.000000
det3dist setAccess 1
det3zeroy 128.000000
@@ -70,25 +53,25 @@ monodescription setAccess 1
# Motor om
om SoftZero 0.000000
om SoftLowerLim -73.000000
-om SoftUpperLim 39.000000
+om SoftUpperLim 134.000000
om Fixed -1.000000
om sign 1.000000
om InterruptMode 0.000000
om AccessCode 2.000000
# Motor stt
stt SoftZero 0.000000
-stt SoftLowerLim -120.000000
-stt SoftUpperLim 120.000000
+stt SoftLowerLim 4.000000
+stt SoftUpperLim 113.000000
stt Fixed -1.000000
stt sign 1.000000
stt InterruptMode 0.000000
stt AccessCode 2.000000
# Motor ch
ch SoftZero 0.000000
-ch SoftLowerLim 0.000000
-ch SoftUpperLim 360.000000
+ch SoftLowerLim 80.000000
+ch SoftUpperLim 212.000000
ch Fixed -1.000000
-ch sign 0.500000
+ch sign 1.000000
ch InterruptMode 0.000000
ch AccessCode 1.000000
# Motor ph
@@ -96,7 +79,7 @@ ph SoftZero 0.000000
ph SoftLowerLim -360.000000
ph SoftUpperLim 360.000000
ph Fixed -1.000000
-ph sign 1.000000
+ph sign -1.000000
ph InterruptMode 0.000000
ph AccessCode 2.000000
# Motor dg3
@@ -136,35 +119,37 @@ phi SoftZero 0.000000
phi SoftLowerLim -360.000000
phi SoftUpperLim 360.000000
phi Fixed -1.000000
-phi sign 1.000000
+phi sign -1.000000
phi InterruptMode 0.000000
phi AccessCode 2.000000
# Motor chi
chi SoftZero 0.000000
-chi SoftLowerLim 0.000000
-chi SoftUpperLim 360.000000
+chi SoftLowerLim 80.000000
+chi SoftUpperLim 212.000000
chi Fixed -1.000000
-chi sign 0.500000
+chi sign 1.000000
chi InterruptMode 0.000000
chi AccessCode 1.000000
# Motor omega
omega SoftZero 0.000000
omega SoftLowerLim -73.000000
-omega SoftUpperLim 39.000000
+omega SoftUpperLim 134.000000
omega Fixed -1.000000
omega sign 1.000000
omega InterruptMode 0.000000
omega AccessCode 2.000000
# Motor twotheta
twotheta SoftZero 0.000000
-twotheta SoftLowerLim -120.000000
-twotheta SoftUpperLim 120.000000
+twotheta SoftLowerLim 4.000000
+twotheta SoftUpperLim 113.000000
twotheta Fixed -1.000000
twotheta sign 1.000000
twotheta InterruptMode 0.000000
twotheta AccessCode 2.000000
lastscancommand cscan a4 10. .1 10 5
lastscancommand setAccess 2
+banana CountMode timer
+banana preset 100.000000
sample_mur 0.000000
sample_mur setAccess 2
email UNKNOWN
@@ -176,7 +161,7 @@ phone setAccess 2
adress UNKNOWN
adress setAccess 2
# Counter counter
-counter SetPreset 5.000000
+counter SetPreset 1.000000
counter SetMode Timer
# Motor som
som SoftZero 0.000000
@@ -444,9 +429,9 @@ a1 InterruptMode 0.000000
a1 AccessCode 2.000000
user Uwe Filges
user setAccess 2
-sample test
+sample D20 30K SNP Okt 2001 GS
sample setAccess 2
-title uwe_test1
+title snp gs apd 30K
title setAccess 2
-starttime 2001-09-21 16:55:53
+starttime 2001-11-16 08:55:46
starttime setAccess 2
diff --git a/sinqhm/SinqHM_bootParamsConfig.c b/sinqhm/SinqHM_bootParamsConfig.c
index 595fbd6a..943d0ddd 100755
--- a/sinqhm/SinqHM_bootParamsConfig.c
+++ b/sinqhm/SinqHM_bootParamsConfig.c
@@ -735,6 +735,11 @@
p_arg[1], p_arg[2], p_arg[3], p_arg[4], p_arg[5],
p_arg[6], p_arg[7], p_arg[8], p_arg[9]);
if (!status) return status;
+ }else if (strcmp (p_arg[0], "trans") == 0) {
+ status = sbpc_config (SQHM__TRANS,
+ p_arg[1], p_arg[2], p_arg[3], p_arg[4], p_arg[5],
+ p_arg[6], p_arg[7], p_arg[8], p_arg[9]);
+ if (!status) return status;
}else {
sprintf (errmsg, "\"%s\" is an unrecognised thing to configure!\n",
p_arg[0]);
diff --git a/sinqhm/SinqHM_srv_filler.c b/sinqhm/SinqHM_srv_filler.c
index 8b7fd25c..1ed029c3 100755
--- a/sinqhm/SinqHM_srv_filler.c
+++ b/sinqhm/SinqHM_srv_filler.c
@@ -1103,12 +1103,16 @@
if (lwl_hdr.ui4 == LWL_FIFO_EMPTY) {
taskDelay (0); /* If FIFO is empty, we can take a breather! */
}else {
+#ifdef TRANNY
VmioBase[VMIO_PORT_A] = 0xff; /* Set timer level (if present) */
*my_char_nxt = lwl_Packet_Read (lwl_hdr.ui4, &my_char_nxt[1]);
Bytes_free = Bytes_free - *my_char_nxt - 1;
my_char_nxt = my_char_nxt + *my_char_nxt + 1;
if (Bytes_free < 24) my_char_nxt = selectNewBuffer (my_char_nxt);
VmioBase[VMIO_PORT_A] = 0x00; /* Reset timer level (if present) */
+#else
+ printf("%#8.8x\n",lwl_hdr.ui4);
+#endif
}
}
/* Our flag has been set. There should be ..
diff --git a/t_conv.c b/t_conv.c
index f68a636f..e0de4c1a 100644
--- a/t_conv.c
+++ b/t_conv.c
@@ -1,3 +1,4 @@
+#line 1 "t_conv.f"
/* t_conv.f -- translated by f2c (version 20000817).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
@@ -5,6 +6,7 @@
#include "f2c.h"
+#line 1 "t_conv.f"
/* Common Block Declarations */
struct {
@@ -28,6 +30,7 @@ static doublereal c_b12 = 360.;
/* Mark Koennecke, November 2000 */
/* ------------------------------------------------------------------------- */
+/*< >*/
/* Subroutine */ int inicurve_(integer *midx, real *mrx1, real *mrx2, integer
*aidx, real *arx1, real *arx2, real *mmin, real *mmax, real *amin,
real *amax)
@@ -36,23 +39,60 @@ static doublereal c_b12 = 360.;
/* Initializes a common with the curvature parameters. */
/* In: monochrmoator curvatuure motor index and parameters */
/* In: analyzer curvature motor index + parameters */
+/*< INTEGER MIDX, AIDX >*/
+/*< REAL*4 MRX1, MRX2, ARX1, ARX2, MMIN, MMAX, AMIN, AMAX >*/
+/*< REAL*4 CM1RX, CM2RX, CA1RX, CA2RX, RMMIN, RMMAX >*/
+/*< REAL*4 RAMIN, RAMAX >*/
+/*< INTEGER ICRM, ICRA, ICLM, INX >*/
+/*< LOGICAL AT_SINQ >*/
+/*< >*/
+/*< ICRM = MIDX >*/
+#line 21 "t_conv.f"
curve_1.icrm = *midx;
+/*< ICRA = AIDX >*/
+#line 22 "t_conv.f"
curve_1.icra = *aidx;
+/*< CM1RX = MRX1 >*/
+#line 23 "t_conv.f"
curve_1.cm1rx = *mrx1;
+/*< CM2RX = MRX2 >*/
+#line 24 "t_conv.f"
curve_1.cm2rx = *mrx2;
+/*< CA1RX = ARX1 >*/
+#line 25 "t_conv.f"
curve_1.ca1rx = *arx1;
+/*< CA2RX = ARX2 >*/
+#line 26 "t_conv.f"
curve_1.ca2rx = *arx2;
+/*< RMMIN = MMIN >*/
+#line 27 "t_conv.f"
curve_1.rmmin = *mmin;
+/*< RMMAX = MMAX >*/
+#line 28 "t_conv.f"
curve_1.rmmax = *mmax;
+/*< RAMIN = AMIN >*/
+#line 29 "t_conv.f"
curve_1.ramin = *amin;
+/*< RAMAX = AMAX >*/
+#line 30 "t_conv.f"
curve_1.ramax = *amax;
+/*< INX = 0 >*/
+#line 31 "t_conv.f"
curve_1.inx = 0;
+/*< AT_SINQ = .TRUE. >*/
+#line 32 "t_conv.f"
curve_1.at_sinq__ = TRUE_;
+/*< ICLM = 0 >*/
+#line 33 "t_conv.f"
curve_1.iclm = 0;
+/*< RETURN >*/
+#line 34 "t_conv.f"
return 0;
+/*< END >*/
} /* inicurve_ */
/* -------------------------------------------------------------------------- */
+/*< >*/
/* Subroutine */ int t_conv__(real *ei, real *aki, real *ef, real *akf, real *
qhkl, real *en, real *hx, real *hy, real *hz, integer *if1, integer *
if2, logical *ldk, logical *ldh, logical *ldf, logical *lpa, real *dm,
@@ -134,17 +174,44 @@ static doublereal c_b12 = 360.;
/* TARGET OF EI(EF) IS UPDATED IS KI(KF) IS DRIVEN */
/* TARGET OF VARIABLE ENERGY IS UPDATED IF EN IS DRIVEN */
/* ----------------------------------------------------------------------- */
-/* File [MAD.SRC]T_CONV.FOR */
+/*< implicit none >*/
+/*< REAL*4 EPS1, EPS4 >*/
+/*< parameter (EPS1 = 1.D-1) >*/
+/*< parameter (EPS4 = 1.D-4) >*/
+/*< INTEGER ICRM, ICRA, ICLM, INX >*/
+/*< LOGICAL AT_SINQ >*/
+/*< REAL*4 CM1RX, CM2RX, CA1RX, CA2RX, RMMIN, RMMAX >*/
+/*< REAL*4 RAMIN, RAMAX >*/
+/*< >*/
/* include 'curve.inc' */
/* ----------------------------------------------------------------------- */
/* Define the dummy arguments */
+/*< real*4 ei, aki, ef, akf, qhkl(3), en >*/
+/*< real*4 hx, hy, hz >*/
+/*< integer*4 if1, if2 >*/
+/*< logical*4 ldk(8), ldh, ldf, lpa >*/
+/*< real*4 dm, da >*/
+/*< real*4 helm, f1h, f1v, f2h, f2v, f >*/
+/*< integer*4 ifx, iss, ism, isa >*/
+/*< real*4 t_a(6), t_rm, t_alm, t_ra, qm >*/
+/*< logical*4 ldra(6), ldr_rm, ldr_alm, ldr_ra >*/
+/*< real*4 p_ih(8), c_ih(4) >*/
+/*< integer*4 ier >*/
/* ----------------------------------------------------------------------- */
/* LOCAL VARIABLES */
+/*< integer*4 i, id, imod, iq >*/
+/*< logical*4 lmoan(2), lqhkle >*/
+/*< double precision a1, a2, a3, a4, a5, a6 >*/
+/*< double precision ala, alm, dakf, daki, dqm, dqs >*/
+/*< double precision def, dei >*/
+/*< double precision ra, rm >*/
+/*< double precision edef(2), akdef(2), dqhkl(3), dqt(3) >*/
+/*< double precision ddm, dda >*/
/* ----------------------------------------------------------------------- */
/* SET UP */
/* IMOD INDEX FOR ERROR TREATMENAT BY ERRESO */
@@ -154,66 +221,164 @@ static doublereal c_b12 = 360.;
/* AND VARIABLE ENERGY IN EDEF(2) */
/* IF ISA IS NUL SET IFX TO 1 AND PUT EF,KF, EQUAL TO EI,KI */
+/*< IMOD = 3 >*/
+#line 135 "t_conv.f"
/* Parameter adjustments */
+#line 135 "t_conv.f"
--c_ih__;
+#line 135 "t_conv.f"
--p_ih__;
+#line 135 "t_conv.f"
--ldra;
+#line 135 "t_conv.f"
--t_a__;
+#line 135 "t_conv.f"
--ldk;
+#line 135 "t_conv.f"
--qhkl;
+#line 135 "t_conv.f"
+#line 135 "t_conv.f"
/* Function Body */
+#line 135 "t_conv.f"
imod = 3;
+/*< DDM = DM >*/
+#line 136 "t_conv.f"
ddm = *dm;
+/*< DDA = DA >*/
+#line 137 "t_conv.f"
dda = *da;
+/*< DO I = 1,2 >*/
+#line 138 "t_conv.f"
for (i__ = 1; i__ <= 2; ++i__) {
+/*< LMOAN(I) = .FALSE. >*/
+#line 139 "t_conv.f"
lmoan[i__ - 1] = FALSE_;
+/*< ENDDO >*/
+#line 140 "t_conv.f"
}
+/*< LQHKLE = .FALSE. >*/
+#line 141 "t_conv.f"
lqhkle = FALSE_;
+/*< DO IQ = 5,8 >*/
+#line 142 "t_conv.f"
for (iq = 5; iq <= 8; ++iq) {
+/*< LQHKLE = (LQHKLE .OR. LDK(IQ)) >*/
+#line 143 "t_conv.f"
lqhkle = lqhkle || ldk[iq];
+/*< ENDDO >*/
+#line 144 "t_conv.f"
}
+/*< DAKI = AKI >*/
+#line 145 "t_conv.f"
daki = *aki;
+/*< DAKF = AKF >*/
+#line 146 "t_conv.f"
dakf = *akf;
+/*< IF (ISA .EQ. 0) IFX = 1 >*/
+#line 147 "t_conv.f"
if (*isa == 0) {
+#line 147 "t_conv.f"
*ifx = 1;
+#line 147 "t_conv.f"
}
+/*< EDEF(IFX) = EI >*/
+#line 148 "t_conv.f"
edef[*ifx - 1] = *ei;
+/*< AKDEF(IFX) = AKI >*/
+#line 149 "t_conv.f"
akdef[*ifx - 1] = *aki;
+/*< EDEF(3-IFX) = EF >*/
+#line 150 "t_conv.f"
edef[3 - *ifx - 1] = *ef;
+/*< AKDEF(3-IFX) = AKF >*/
+#line 151 "t_conv.f"
akdef[3 - *ifx - 1] = *akf;
+/*< IF( ISA .EQ. 0) THEN >*/
+#line 152 "t_conv.f"
if (*isa == 0) {
+/*< EDEF(2) = EDEF(1) >*/
+#line 153 "t_conv.f"
edef[1] = edef[0];
+/*< AKDEF(2) = AKDEF(1) >*/
+#line 154 "t_conv.f"
akdef[1] = akdef[0];
+/*< LDK(3) = .TRUE. >*/
+#line 155 "t_conv.f"
ldk[3] = TRUE_;
+/*< LDK(4) = .TRUE. >*/
+#line 156 "t_conv.f"
ldk[4] = TRUE_;
+/*< T_A(5) = 0. >*/
+#line 157 "t_conv.f"
t_a__[5] = 0.f;
+/*< T_A(6) = 0. >*/
+#line 158 "t_conv.f"
t_a__[6] = 0.f;
+/*< LDRA(5) = .TRUE. >*/
+#line 159 "t_conv.f"
ldra[5] = TRUE_;
+/*< LDRA(6) = .TRUE. >*/
+#line 160 "t_conv.f"
ldra[6] = TRUE_;
+/*< ENDIF >*/
+#line 161 "t_conv.f"
}
/* ----------------------------------------------------------------------- */
/* FIRST TAKE IN ACCOUNT THE FIXED ENERGY PB */
+/*< IF (LDK(2*IFX-1) .OR. LDK(2*IFX)) THEN >*/
+#line 165 "t_conv.f"
if (ldk[(*ifx << 1) - 1] || ldk[*ifx * 2]) {
+/*< LMOAN(IFX) = .TRUE. >*/
+#line 166 "t_conv.f"
lmoan[*ifx - 1] = TRUE_;
+/*< IF (LDK(2*IFX-1)) THEN >*/
+#line 167 "t_conv.f"
if (ldk[(*ifx << 1) - 1]) {
- *ier = 1;
+/*< IER = 1 + 8 >*/
+#line 168 "t_conv.f"
+ *ier = 9;
+/*< IF(EDEF(1) .LT. EPS1) GOTO 999 >*/
+#line 169 "t_conv.f"
if (edef[0] < .1f) {
+#line 169 "t_conv.f"
goto L999;
+#line 169 "t_conv.f"
}
+/*< IER = 0 >*/
+#line 170 "t_conv.f"
*ier = 0;
+/*< AKDEF(1) = SQRT(EDEF(1)/F) >*/
+#line 171 "t_conv.f"
akdef[0] = sqrt(edef[0] / *f);
+/*< ELSE >*/
+#line 172 "t_conv.f"
} else {
- *ier = 1;
+/*< IER = 1 + 8 >*/
+#line 173 "t_conv.f"
+ *ier = 9;
+/*< IF(AKDEF(1) .LT. EPS1) GOTO 999 >*/
+#line 174 "t_conv.f"
if (akdef[0] < .1f) {
+#line 174 "t_conv.f"
goto L999;
+#line 174 "t_conv.f"
}
+/*< IER = 0 >*/
+#line 175 "t_conv.f"
*ier = 0;
+/*< EDEF(1) = F*AKDEF(1)**2 >*/
/* Computing 2nd power */
+#line 176 "t_conv.f"
d__1 = akdef[0];
+#line 176 "t_conv.f"
edef[0] = *f * (d__1 * d__1);
+/*< ENDIF >*/
+#line 177 "t_conv.f"
}
+/*< ENDIF >*/
+#line 178 "t_conv.f"
}
/* ----------------------------------------------------------------------- */
/* NOW TAKE IN ACCOUNT THE VARIABLE ENERGY PB */
@@ -229,127 +394,331 @@ static doublereal c_b12 = 360.;
/* IF KF IS CONSTANT USE ALWAYS THE VALUE IN TARGET AND */
/* DO A DRIVE OF KF TO KEEP A5/A6 IN RIGHT POSITION */
+/*< IF (LDK(5-2*IFX) .OR. LDK(6-2*IFX)) THEN >*/
+#line 193 "t_conv.f"
if (ldk[5 - (*ifx << 1)] || ldk[6 - (*ifx << 1)]) {
+/*< LMOAN(3-IFX) = .TRUE. >*/
+#line 194 "t_conv.f"
lmoan[3 - *ifx - 1] = TRUE_;
+/*< IF (LDK(5-2*IFX)) THEN >*/
+#line 195 "t_conv.f"
if (ldk[5 - (*ifx << 1)]) {
- *ier = 1;
+/*< IER = 1 + 8 >*/
+#line 196 "t_conv.f"
+ *ier = 9;
+/*< IF(EDEF(2) .LT. EPS4) GOTO 999 >*/
+#line 197 "t_conv.f"
if (edef[1] < 1e-4f) {
+#line 197 "t_conv.f"
goto L999;
+#line 197 "t_conv.f"
}
+/*< IER = 0 >*/
+#line 198 "t_conv.f"
*ier = 0;
+/*< AKDEF(2) = SQRT(EDEF(2)/F) >*/
+#line 199 "t_conv.f"
akdef[1] = sqrt(edef[1] / *f);
+/*< ELSE >*/
+#line 200 "t_conv.f"
} else {
- *ier = 1;
+/*< IER = 1 + 8 >*/
+#line 201 "t_conv.f"
+ *ier = 9;
+/*< IF(AKDEF(2) .LT. EPS4) GOTO 999 >*/
+#line 202 "t_conv.f"
if (akdef[1] < 1e-4f) {
+#line 202 "t_conv.f"
goto L999;
+#line 202 "t_conv.f"
}
+/*< IER = 0 >*/
+#line 203 "t_conv.f"
*ier = 0;
+/*< EDEF(2) = F*AKDEF(2)**2 >*/
/* Computing 2nd power */
+#line 204 "t_conv.f"
d__1 = akdef[1];
+#line 204 "t_conv.f"
edef[1] = *f * (d__1 * d__1);
+/*< ENDIF >*/
+#line 205 "t_conv.f"
}
+/*< EN = (3-2*IFX)*(EDEF(1)-EDEF(2)) >*/
+#line 206 "t_conv.f"
*en = (3 - (*ifx << 1)) * (edef[0] - edef[1]);
+/*< ELSEIF (LQHKLE) THEN >*/
+#line 207 "t_conv.f"
} else if (lqhkle) {
+/*< LMOAN(3-IFX) = .TRUE. >*/
+#line 208 "t_conv.f"
lmoan[3 - *ifx - 1] = TRUE_;
+/*< EDEF(2) = EDEF(1)+(2*IFX-3)*EN >*/
+#line 209 "t_conv.f"
edef[1] = edef[0] + ((*ifx << 1) - 3) * *en;
- *ier = 1;
+/*< IER = 1 + 8 >*/
+#line 210 "t_conv.f"
+ *ier = 9;
+/*< IF(EDEF(2) .LT. EPS4) GOTO 999 >*/
+#line 211 "t_conv.f"
if (edef[1] < 1e-4f) {
+#line 211 "t_conv.f"
goto L999;
+#line 211 "t_conv.f"
}
+/*< IER = 0 >*/
+#line 212 "t_conv.f"
*ier = 0;
+/*< AKDEF(2) = SQRT(EDEF(2)/F) >*/
+#line 213 "t_conv.f"
akdef[1] = sqrt(edef[1] / *f);
+/*< ENDIF >*/
+#line 214 "t_conv.f"
}
/* ----------------------------------------------------------------------- */
/* CALCULATE MONOCHROMATOR AND ANALYSER ANGLES */
+/*< IF(LMOAN(1)) THEN >*/
+#line 218 "t_conv.f"
if (lmoan[0]) {
+/*< DEI = EDEF(IFX) >*/
+#line 219 "t_conv.f"
dei = edef[*ifx - 1];
+/*< DAKI = AKDEF(IFX) >*/
+#line 220 "t_conv.f"
daki = akdef[*ifx - 1];
+/*< CALL EX_CASE(DDM,ISM,DAKI,A1,A2,RM,ALM,0,IER) >*/
+#line 221 "t_conv.f"
ex_case__(&ddm, ism, &daki, &a1, &a2, &rm, &alm, &c__0, ier);
+/*< IF (IER .EQ. 0) THEN >*/
+#line 222 "t_conv.f"
if (*ier == 0) {
+/*< AKI = DAKI >*/
+#line 223 "t_conv.f"
*aki = daki;
+/*< EI = DEI >*/
+#line 224 "t_conv.f"
*ei = dei;
+/*< T_A(1) = A1 >*/
+#line 225 "t_conv.f"
t_a__[1] = a1;
+/*< T_A(2) = A2 >*/
+#line 226 "t_conv.f"
t_a__[2] = a2;
+/*< LDRA(1) = .TRUE. >*/
+#line 227 "t_conv.f"
ldra[1] = TRUE_;
+/*< LDRA(2) = .TRUE. >*/
+#line 228 "t_conv.f"
ldra[2] = TRUE_;
+/*< if (icrm .ne. 0) then >*/
+#line 229 "t_conv.f"
if (curve_1.icrm != 0) {
+/*< T_RM = RM >*/
+#line 230 "t_conv.f"
*t_rm__ = rm;
+/*< LDR_RM = .TRUE. >*/
+#line 231 "t_conv.f"
*ldr_rm__ = TRUE_;
+/*< endif >*/
+#line 232 "t_conv.f"
}
+/*< if ((iclm .ne. 0) .and. (inx .ne. 0)) then >*/
+#line 233 "t_conv.f"
if (curve_1.iclm != 0 && curve_1.inx != 0) {
+/*< T_ALM = ALM >*/
+#line 234 "t_conv.f"
*t_alm__ = alm;
+/*< LDR_ALM = .TRUE. >*/
+#line 235 "t_conv.f"
*ldr_alm__ = TRUE_;
+/*< endif >*/
+#line 236 "t_conv.f"
}
+/*< ELSE >*/
+#line 237 "t_conv.f"
} else {
+/*< IER = IER + 8 >*/
+#line 238 "t_conv.f"
+ *ier += 8;
+/*< GOTO 999 >*/
+#line 239 "t_conv.f"
goto L999;
+/*< ENDIF >*/
+#line 240 "t_conv.f"
}
+/*< ENDIF >*/
+#line 241 "t_conv.f"
}
+/*< IF(LMOAN(2)) THEN >*/
+#line 242 "t_conv.f"
if (lmoan[1]) {
+/*< DEF = EDEF(3-IFX) >*/
+#line 243 "t_conv.f"
def = edef[3 - *ifx - 1];
+/*< DAKF = AKDEF(3-IFX) >*/
+#line 244 "t_conv.f"
dakf = akdef[3 - *ifx - 1];
+/*< CALL EX_CASE(DDA,ISA,DAKF,A5,A6,RA,ALA,1,IER) >*/
+#line 245 "t_conv.f"
ex_case__(&dda, isa, &dakf, &a5, &a6, &ra, &ala, &c__1, ier);
+/*< IF (IER .EQ. 0) THEN >*/
+#line 246 "t_conv.f"
if (*ier == 0) {
+/*< AKF = DAKF >*/
+#line 247 "t_conv.f"
*akf = dakf;
+/*< EF = DEF >*/
+#line 248 "t_conv.f"
*ef = def;
+/*< T_A(5) = A5 >*/
+#line 249 "t_conv.f"
t_a__[5] = a5;
+/*< T_A(6) = A6 >*/
+#line 250 "t_conv.f"
t_a__[6] = a6;
+/*< LDRA(5) = .TRUE. >*/
+#line 251 "t_conv.f"
ldra[5] = TRUE_;
+/*< LDRA(6) = .TRUE. >*/
+#line 252 "t_conv.f"
ldra[6] = TRUE_;
+/*< if (icra .ne. 0) then >*/
+#line 253 "t_conv.f"
if (curve_1.icra != 0) {
+/*< T_RA = RA >*/
+#line 254 "t_conv.f"
*t_ra__ = ra;
+/*< LDR_RA = .TRUE. >*/
+#line 255 "t_conv.f"
*ldr_ra__ = TRUE_;
+/*< endif >*/
+#line 256 "t_conv.f"
}
+/*< ELSE >*/
+#line 257 "t_conv.f"
} else {
+/*< IER = IER + 8 >*/
+#line 258 "t_conv.f"
+ *ier += 8;
+/*< GOTO 999 >*/
+#line 259 "t_conv.f"
goto L999;
+/*< ENDIF >*/
+#line 260 "t_conv.f"
}
+/*< ENDIF >*/
+#line 261 "t_conv.f"
}
/* ----------------------------------------------------------------------- */
/* USE (QH,QK,QL) TO CALCULATE A3 AND A4 */
/* CALCULATE Q1 AND Q2 IN SCATTERING PLANE */
+/*< IMOD = 2 >*/
+#line 266 "t_conv.f"
imod = 2;
+/*< IF (LQHKLE) THEN >*/
+#line 267 "t_conv.f"
if (lqhkle) {
+/*< DO ID = 1,3 >*/
+#line 268 "t_conv.f"
for (id = 1; id <= 3; ++id) {
+/*< DQHKL(ID) = QHKL(ID) >*/
+#line 269 "t_conv.f"
dqhkl[id - 1] = qhkl[id];
+/*< ENDDO >*/
+#line 270 "t_conv.f"
}
+/*< CALL RL2SPV(DQHKL,DQT,DQM,DQS,IER) >*/
+#line 271 "t_conv.f"
rl2spv_(dqhkl, dqt, &dqm, &dqs, ier);
+/*< IF (IER .NE. 0) GOTO 999 >*/
+#line 272 "t_conv.f"
if (*ier != 0) {
+#line 272 "t_conv.f"
goto L999;
+#line 272 "t_conv.f"
}
+/*< CALL SAM_CASE(DQT,DQM,DQS,DAKI,DAKF,A3,A4,ISS,IER) >*/
+#line 273 "t_conv.f"
sam_case__(dqt, &dqm, &dqs, &daki, &dakf, &a3, &a4, iss, ier);
+/*< IF (IER .EQ. 0) THEN >*/
+#line 274 "t_conv.f"
if (*ier == 0) {
+/*< T_A(3) = A3 >*/
+#line 275 "t_conv.f"
t_a__[3] = a3;
+/*< T_A(4) = A4 >*/
+#line 276 "t_conv.f"
t_a__[4] = a4;
+/*< LDRA(3) = .TRUE. >*/
+#line 277 "t_conv.f"
ldra[3] = TRUE_;
+/*< LDRA(4) = .TRUE. >*/
+#line 278 "t_conv.f"
ldra[4] = TRUE_;
+/*< QM = DQM >*/
+#line 279 "t_conv.f"
*qm = dqm;
+/*< ELSE >*/
+#line 280 "t_conv.f"
} else {
+/*< IER = IER + 4 >*/
+#line 281 "t_conv.f"
+ *ier += 4;
+/*< GOTO 999 >*/
+#line 282 "t_conv.f"
goto L999;
+/*< ENDIF >*/
+#line 283 "t_conv.f"
}
+/*< ENDIF >*/
+#line 284 "t_conv.f"
}
/* ----------------------------------------------------------------------- */
/* DEAL WITH FLIPPERS AND HELMOTZ COILS IF LPA */
+/*< IF (LPA .AND. (LMOAN(1) .OR. LMOAN(2))) THEN >*/
+#line 288 "t_conv.f"
if (*lpa && (lmoan[0] || lmoan[1])) {
+/*< IF >*/
+#line 289 "t_conv.f"
if (*ldf) {
+#line 289 "t_conv.f"
flip_case__(if1, if2, &p_ih__[1], f1v, f1h, f2v, f2h, aki, akf,
ier);
+#line 289 "t_conv.f"
}
+/*< IF >*/
+#line 291 "t_conv.f"
if (*ldh) {
+#line 291 "t_conv.f"
helm_case__(hx, hy, hz, &p_ih__[1], &c_ih__[1], aki, akf, &a4, qm,
helm, ier);
+#line 291 "t_conv.f"
}
+/*< endif >*/
+#line 293 "t_conv.f"
}
/* ----------------------------------------------------------------------- */
+/*< 999 CONTINUE >*/
+#line 295 "t_conv.f"
L999:
+/*< IF (IER .NE. 0) CALL ERRESO(IMOD,IER) >*/
+#line 296 "t_conv.f"
if (*ier != 0) {
+#line 296 "t_conv.f"
erreso_(&imod, ier);
+#line 296 "t_conv.f"
}
+/*< RETURN >*/
+#line 297 "t_conv.f"
return 0;
+/*< END >*/
} /* t_conv__ */
+/*< SUBRO >*/
/* Subroutine */ int ex_case__(doublereal *dx, integer *isx, doublereal *akx,
doublereal *ax1, doublereal *ax2, doublereal *rx, doublereal *alx,
integer *mon_or_anal__, integer *ier)
@@ -392,61 +761,128 @@ L999:
/* 2 'KI OR KF TOO SMALL', */
/* 3 'KI OR KF TOO BIG', */
/* ----------------------------------------------------------------------- */
-/* Part of T_CONV.FOR */
+/*< implicit none >*/
+/*< double precision PI, RD, EPS1 >*/
+/*< PARAMETER (PI = 3.14159265358979323846264338327950D0) >*/
+/*< PARAMETER (RD = 57.29577951308232087679815481410517D0) >*/
+/*< PARAMETER (EPS1 = 1.D-1) >*/
/* ----------------------------------------------------------------------- */
/* Define the dummy arguments */
+/*< double precision dx >*/
+/*< integer*4 isx >*/
+/*< double precision akx, ax1, ax2, rx, alx >*/
+/*< integer*4 mon_or_anal, ier >*/
/* ----------------------------------------------------------------------- */
/* include 'curve.inc' */
/* include 'motdef.inc' */
/* include 'iolsddef.inc' */
+/*< INTEGER ICRM, ICRA, ICLM, INX >*/
+/*< LOGICAL AT_SINQ >*/
+/*< REAL*4 CM1RX, CM2RX, CA1RX, CA2RX, RMMIN, RMMAX >*/
+/*< REAL*4 RAMIN, RAMAX >*/
+/*< >*/
/* real*4 tbut(5,NBMOT) */
/* equivalence (rbut, tbut) */
/* ----------------------------------------------------------------------- */
/* LOCAL VAR */
+/*< double precision arg, dc1rx, dc2rx, drmin, drmax, dcl1r, my_rx >*/
+/*< integer*4 ios, indx >*/
/* ----------------------------------------------------------------------- */
/* INIT AND TEST */
+/*< ier = 0 >*/
+#line 363 "t_conv.f"
*ier = 0;
+/*< ax1 = 0.0 >*/
+#line 364 "t_conv.f"
*ax1 = 0.f;
+/*< ax2 = 0.0 >*/
+#line 365 "t_conv.f"
*ax2 = 0.f;
+/*< rx = 0.0 >*/
+#line 366 "t_conv.f"
*rx = 0.f;
+/*< alx = 0.0 >*/
+#line 367 "t_conv.f"
*alx = 0.f;
/* ---------------------------------------------------------------------- */
/* Check validity of inputs. */
+/*< if (dx .lt. EPS1) ier = 1 >*/
+#line 370 "t_conv.f"
if (*dx < .1) {
+#line 370 "t_conv.f"
*ier = 1;
+#line 370 "t_conv.f"
}
+/*< if (akx .lt. EPS1) ier = 2 >*/
+#line 371 "t_conv.f"
if (*akx < .1) {
+#line 371 "t_conv.f"
*ier = 2;
+#line 371 "t_conv.f"
}
+/*< arg = PI/(dx * akx) >*/
+#line 372 "t_conv.f"
arg = 3.1415926535897932384626433832795 / (*dx * *akx);
+/*< if (abs(arg) .gt. 1.0) ier = 3 >*/
+#line 373 "t_conv.f"
if (abs(arg) > 1.f) {
+#line 373 "t_conv.f"
*ier = 3;
+#line 373 "t_conv.f"
}
+/*< if (ier .ne. 0) goto 999 >*/
+#line 374 "t_conv.f"
if (*ier != 0) {
+#line 374 "t_conv.f"
goto L999;
+#line 374 "t_conv.f"
}
/* ---------------------------------------------------------------------- */
+/*< if (mon_or_anal .eq. 0) then ! Use monochr or anal params? >*/
+#line 376 "t_conv.f"
if (*mon_or_anal__ == 0) {
-/* Use monochr or anal params? */
+/*< indx = icrm ! Monochr, so set up params. >*/
+#line 377 "t_conv.f"
indx = curve_1.icrm;
-/* Monochr, so set up params. */
+/*< dc1rx = cm1rx >*/
+#line 378 "t_conv.f"
dc1rx = curve_1.cm1rx;
+/*< dc2rx = cm2rx >*/
+#line 379 "t_conv.f"
dc2rx = curve_1.cm2rx;
+/*< dcl1r = ICLM >*/
+#line 380 "t_conv.f"
dcl1r = (doublereal) curve_1.iclm;
+/*< drmin = rmmin >*/
+#line 381 "t_conv.f"
drmin = curve_1.rmmin;
+/*< drmax = rmmax >*/
+#line 382 "t_conv.f"
drmax = curve_1.rmmax;
+/*< else >*/
+#line 383 "t_conv.f"
} else {
+/*< indx = icra ! Analyser, so set up params. >*/
+#line 384 "t_conv.f"
indx = curve_1.icra;
-/* Analyser, so set up params. */
+/*< dc1rx = ca1rx ! There is no ALX in this case. >*/
+#line 385 "t_conv.f"
dc1rx = curve_1.ca1rx;
-/* There is no ALX in this case. */
+/*< dc2rx = ca2rx >*/
+#line 386 "t_conv.f"
dc2rx = curve_1.ca2rx;
+/*< drmin = ramin >*/
+#line 387 "t_conv.f"
drmin = curve_1.ramin;
+/*< drmax = ramax >*/
+#line 388 "t_conv.f"
drmax = curve_1.ramax;
+/*< endif >*/
+#line 389 "t_conv.f"
}
/* if (indx .ne. 0) then ! Include zero-offset in min/max */
@@ -458,17 +894,33 @@ L999:
/* ----------------------------------------------------------------------- */
/* Calculation of the two angles */
+/*< if (isx .eq. 0) then ! "Straight-through" mode? >*/
+#line 400 "t_conv.f"
if (*isx == 0) {
-/* "Straight-through" mode? */
+/*< ax1 = 0.0 ! Yes. >*/
+#line 401 "t_conv.f"
*ax1 = 0.f;
-/* Yes. */
+/*< ax2 = 0.0 >*/
+#line 402 "t_conv.f"
*ax2 = 0.f;
+/*< rx = drmin >*/
+#line 403 "t_conv.f"
*rx = drmin;
+/*< alx = 0.0 >*/
+#line 404 "t_conv.f"
*alx = 0.f;
+/*< return >*/
+#line 405 "t_conv.f"
return 0;
+/*< endif >*/
+#line 406 "t_conv.f"
}
+/*< ax1 = asin (arg) * isx * rd >*/
+#line 408 "t_conv.f"
*ax1 = asin(arg) * *isx * 57.29577951308232087679815481410517;
+/*< ax2 = 2.0d0 * ax1 >*/
+#line 409 "t_conv.f"
*ax2 = *ax1 * 2.;
/* ----------------------------------------------------------------------- */
/* Calculation of mono curvature RM or analyser curvature RA */
@@ -487,40 +939,64 @@ L999:
/* rmmin = 0. */
/* rmmax = 20. */
/* ----------------------------------------------------------------------- */
+/*< if (mon_or_anal .eq. 0) then ! Monochr or analyser? >*/
+#line 427 "t_conv.f"
if (*mon_or_anal__ == 0) {
-/* Monochr or analyser? */
+/*< if (inx .ne. 0) then ! Monochr. Is there a translation? >*/
+#line 428 "t_conv.f"
if (curve_1.inx != 0) {
-/* Monochr. Is there a translation? */
+/*< if (iclm .ne. 0) then ! Yes, IN8 case. If there's a .. >*/
+#line 429 "t_conv.f"
if (curve_1.iclm != 0) {
-/* Yes, IN8 case. If there's a .. */
+/*< alx = (dcl1r/sin(ax2/rd)) * cos(ax2/rd) ! .. motor, do the .. >*/
+#line 430 "t_conv.f"
*alx = dcl1r / sin(*ax2 / 57.29577951308232087679815481410517)
* cos(*ax2 / 57.29577951308232087679815481410517);
-/* .. motor, do the .. */
+/*< rx = dc2rx * sqrt(sin(abs(ax2)/rd)) - dc1rx ! .. calculation. >*/
+#line 431 "t_conv.f"
*rx = dc2rx * sqrt(sin(abs(*ax2) /
57.29577951308232087679815481410517)) - dc1rx;
-/* .. calculation. */
+/*< rx = dmin1 (dmax1 (rx, drmin), drmax) >*/
/* Computing MIN */
+#line 432 "t_conv.f"
d__1 = max(*rx,drmin);
+#line 432 "t_conv.f"
*rx = min(d__1,drmax);
+/*< return >*/
+#line 433 "t_conv.f"
return 0;
+/*< endif >*/
+#line 434 "t_conv.f"
}
+/*< else ! Not IN8 case so, .. >*/
+#line 435 "t_conv.f"
} else {
-/* Not IN8 case so, .. */
+/*< my_rx = dc1rx + dc2rx/sin(abs(ax1)/rd) ! .. simply calculate. >*/
+#line 436 "t_conv.f"
my_rx__ = dc1rx + dc2rx / sin(abs(*ax1) /
57.29577951308232087679815481410517);
-/* .. simply calculate. */
+/*< endif >*/
+#line 437 "t_conv.f"
}
+/*< else ! Analyser. >*/
+#line 438 "t_conv.f"
} else {
-/* Analyser. */
+/*< my_rx = dc1rx + dc2rx * sin(abs(ax1)/rd) ! Simply calculate. >*/
+#line 439 "t_conv.f"
my_rx__ = dc1rx + dc2rx * sin(abs(*ax1) /
57.29577951308232087679815481410517);
-/* Simply calculate. */
+/*< endif >*/
+#line 440 "t_conv.f"
}
+/*< if (indx .ne. 0) then ! If there's a motor, return the curvature. >*/
+#line 442 "t_conv.f"
if (indx != 0) {
-/* If there's a motor, return the curvature. */
+/*< rx = dmin1 (dmax1 (my_rx, drmin), drmax) >*/
/* Computing MIN */
+#line 443 "t_conv.f"
d__1 = max(my_rx__,drmin);
+#line 443 "t_conv.f"
*rx = min(d__1,drmax);
/* if (rx .ne. my_rx) then */
/* write (iolun, 101, iostat=ios) motnam(indx), my_rx */
@@ -528,13 +1004,21 @@ L999:
/* + 'or high limits!'/ */
/* + ' Calculated curvature was', f9.2) */
/* endif */
+/*< endif >*/
+#line 450 "t_conv.f"
}
/* ----------------------------------------------------------------------- */
+/*< 999 continue >*/
+#line 452 "t_conv.f"
L999:
+/*< return >*/
+#line 453 "t_conv.f"
return 0;
+/*< end >*/
} /* ex_case__ */
+/*< SUBRO >*/
/* Subroutine */ int sam_case__(doublereal *qt, doublereal *qm, doublereal *
qs, doublereal *aki, doublereal *akf, doublereal *ax3, doublereal *
ax4, integer *iss, integer *ier)
@@ -568,47 +1052,94 @@ L999:
/* 3 'Q MODULUS TOO SMALL', */
/* 4 'Q MODULUS TOO BIG', */
/* ----------------------------------------------------------------------- */
-/* Part of T_CONV.FOR */
+/*< implicit none >*/
+/*< double precision RD, EPS3, EPS6 >*/
+/*< PARAMETER (RD = 57.29577951308232087679815481410517D0) >*/
+/*< PARAMETER (EPS3 = 1.D-3) >*/
+/*< PARAMETER (EPS6 = 1.D-6) >*/
/* ----------------------------------------------------------------------- */
/* Define the dummy arguments */
+/*< double precision qt(3) >*/
+/*< double precision qm, qs, aki, akf, ax3, ax4 >*/
+/*< integer*4 iss, ier >*/
/* ----------------------------------------------------------------------- */
/* Local variables */
+/*< double precision arg, sax3 >*/
/* ----------------------------------------------------------------------- */
/* INIT AND TEST */
+/*< IER = 0 >*/
+#line 494 "t_conv.f"
/* Parameter adjustments */
+#line 494 "t_conv.f"
--qt;
+#line 494 "t_conv.f"
+#line 494 "t_conv.f"
/* Function Body */
+#line 494 "t_conv.f"
*ier = 0;
+/*< IF ((ABS(QS) .LT. EPS6) .OR. (ABS(QM) .LT. EPS3)) THEN >*/
+#line 495 "t_conv.f"
if (abs(*qs) < 1e-6 || abs(*qm) < .001) {
+/*< IER = 3 >*/
+#line 496 "t_conv.f"
*ier = 3;
+/*< GOTO 999 >*/
+#line 497 "t_conv.f"
goto L999;
+/*< ENDIF >*/
+#line 498 "t_conv.f"
}
/* ----------------------------------------------------------------------- */
/* CALCULATE A3 AND MOVE IT INTHE -180 ,+180 INTERVAL */
+/*< ARG = (AKI**2 + AKF**2 - QS)/(2.D0*AKI*AKF) >*/
/* Computing 2nd power */
+#line 502 "t_conv.f"
d__1 = *aki;
/* Computing 2nd power */
+#line 502 "t_conv.f"
d__2 = *akf;
+#line 502 "t_conv.f"
arg = (d__1 * d__1 + d__2 * d__2 - *qs) / (*aki * 2. * *akf);
+/*< IF(ABS(ARG) .GT. 1.D0)THEN >*/
+#line 503 "t_conv.f"
if (abs(arg) > 1.) {
+/*< IER = 4 >*/
+#line 504 "t_conv.f"
*ier = 4;
+/*< GOTO 999 >*/
+#line 505 "t_conv.f"
goto L999;
+/*< ELSE >*/
+#line 506 "t_conv.f"
} else {
+/*< AX4 = ACOS(ARG)*ISS*RD >*/
+#line 507 "t_conv.f"
*ax4 = acos(arg) * *iss * 57.29577951308232087679815481410517;
+/*< ENDIF >*/
+#line 508 "t_conv.f"
}
+/*< >*/
/* Computing 2nd power */
+#line 509 "t_conv.f"
d__1 = *akf;
/* Computing 2nd power */
+#line 509 "t_conv.f"
d__2 = *aki;
+#line 509 "t_conv.f"
*ax3 = (-atan2(qt[2], qt[1]) - acos((d__1 * d__1 - *qs - d__2 * d__2) / (*
qm * -2. * *aki)) * d_sign(&c_b10, ax4)) *
57.29577951308232087679815481410517;
+/*< sax3 = Dsign(1.D0,ax3) >*/
+#line 511 "t_conv.f"
sax3 = d_sign(&c_b10, ax3);
+/*< AX3 = DMOD(AX3+sax3*180.D0,360.D0)-sax3*180.D0 >*/
+#line 512 "t_conv.f"
d__1 = *ax3 + sax3 * 180.;
+#line 512 "t_conv.f"
*ax3 = d_mod(&d__1, &c_b12) - sax3 * 180.;
/* IF(LPLATE) AX3 = -ATAN(SIN(AX4/RD)/(LSA*TAN(AX5/RD)/(ALMS*C */
@@ -617,11 +1148,17 @@ L999:
/* IF( A3 .LT. -180.D0) AX3 = AX3+360.D0 */
/* IF(LPLATE .AND. (A3 .GT. 0.0)) AX3 = AX3-180 */
/* C----------------------------------------------------------------------- */
+/*< 999 CONTINUE >*/
+#line 520 "t_conv.f"
L999:
+/*< RETURN >*/
+#line 521 "t_conv.f"
return 0;
+/*< END >*/
} /* sam_case__ */
+/*< SUBRO >*/
/* Subroutine */ int helm_case__(real *hx, real *hy, real *hz, real *t_ih__,
real *c_ih__, real *aki, real *akf, doublereal *a4, real *qm, real *
helm, integer *ier)
@@ -659,41 +1196,90 @@ L999:
/* There is a 3rd coil for Hz. */
/* ----------------------------------------------------------------------- */
-/* Part of T_CONV.FOR */
+/*< implicit none >*/
/* include 'common_sinq.inc' */
+/*< double precision PI, RD, EPS3, EPS4 >*/
+/*< PARAMETER (PI = 3.14159265358979323846264338327950D0) >*/
+/*< PARAMETER (RD = 57.29577951308232087679815481410517D0) >*/
+/*< PARAMETER (EPS3 = 1.0D-3) >*/
+/*< PARAMETER (EPS4 = 1.0D-4) >*/
/* ----------------------------------------------------------------------- */
/* Define the dummy arguments */
+/*< real*4 hx, hy, hz >*/
+/*< real*4 t_ih(8) >*/
+/*< real*4 c_ih(4) >*/
+/*< real*4 aki, akf >*/
+/*< double precision a4 >*/
+/*< real*4 qm, helm >*/
+/*< integer*4 ier >*/
+/*< LOGICAL AT_SINQ >*/
/* ----------------------------------------------------------------------- */
/* Local variables */
+/*< integer*4 ic, ncoef >*/
+/*< double precision hdir, hdir2, hrad, phi, qpar, qperp >*/
/* ----------------------------------------------------------------------- */
/* INIT AND TEST */
+/*< AT_SINQ = .TRUE. >*/
+#line 569 "t_conv.f"
/* Parameter adjustments */
+#line 569 "t_conv.f"
--c_ih__;
+#line 569 "t_conv.f"
--t_ih__;
+#line 569 "t_conv.f"
+#line 569 "t_conv.f"
/* Function Body */
+#line 569 "t_conv.f"
at_sinq__ = TRUE_;
+/*< ncoef = 4 >*/
+#line 570 "t_conv.f"
ncoef = 4;
+/*< if (at_sinq) ncoef = 3 >*/
+#line 571 "t_conv.f"
if (at_sinq__) {
+#line 571 "t_conv.f"
ncoef = 3;
+#line 571 "t_conv.f"
}
+/*< IER = 1 >*/
+#line 573 "t_conv.f"
*ier = 1;
+/*< IF (ABS(QM) .LT. EPS4) goto 999 >*/
+#line 574 "t_conv.f"
if (dabs(*qm) < 1e-4) {
+#line 574 "t_conv.f"
goto L999;
+#line 574 "t_conv.f"
}
+/*< IER = 0 >*/
+#line 575 "t_conv.f"
*ier = 0;
+/*< DO IC = 1,ncoef >*/
+#line 576 "t_conv.f"
i__1 = ncoef;
+#line 576 "t_conv.f"
for (ic = 1; ic <= i__1; ++ic) {
+/*< IF (C_IH(IC) .LT. EPS4) IER = 2 >*/
+#line 577 "t_conv.f"
if (c_ih__[ic] < 1e-4) {
+#line 577 "t_conv.f"
*ier = 2;
+#line 577 "t_conv.f"
}
+/*< ENDDO >*/
+#line 578 "t_conv.f"
}
+/*< IF (IER .NE. 0) GOTO 999 >*/
+#line 579 "t_conv.f"
if (*ier != 0) {
+#line 579 "t_conv.f"
goto L999;
+#line 579 "t_conv.f"
}
/* ----------------------------------------------------------------------- */
/* CALCULATE MODULE AND ANGLES OF IN PLANE FIELD H */
@@ -702,92 +1288,217 @@ L999:
/* HDIR ! Direction of H relative to PHI (in radians) */
/* HDIR2 ! Angle between field and axis of Coil 1 (in radians) */
+/*< qpar = aki - akf * cos(a4/RD) >*/
+#line 587 "t_conv.f"
qpar = *aki - *akf * cos(*a4 / 57.29577951308232087679815481410517);
+/*< qperp = -akf * sin(a4/RD) >*/
+#line 588 "t_conv.f"
qperp = -(*akf) * sin(*a4 / 57.29577951308232087679815481410517);
+/*< if (abs(qpar) .gt. EPS3 .and. abs(qperp) .gt. EPS3) then >*/
+#line 589 "t_conv.f"
if (abs(qpar) > .001 && abs(qperp) > .001) {
+/*< phi = atan2 (abs(qperp), abs(qpar)) >*/
+#line 590 "t_conv.f"
phi = atan2((abs(qperp)), (abs(qpar)));
+/*< if (qpar .gt. 0 .and. qperp .lt. 0) then >*/
+#line 591 "t_conv.f"
if (qpar > 0. && qperp < 0.) {
+/*< phi = -phi >*/
+#line 592 "t_conv.f"
phi = -phi;
+/*< elseif (qpar .lt. 0 .and. qperp .gt. 0) then >*/
+#line 593 "t_conv.f"
} else if (qpar < 0. && qperp > 0.) {
+/*< phi = PI - phi >*/
+#line 594 "t_conv.f"
phi = 3.1415926535897932384626433832795 - phi;
+/*< elseif (qpar .lt. 0 .and. qperp .lt. 0) then >*/
+#line 595 "t_conv.f"
} else if (qpar < 0. && qperp < 0.) {
+/*< phi = phi - PI >*/
+#line 596 "t_conv.f"
phi += -3.1415926535897932384626433832795;
+/*< endif >*/
+#line 597 "t_conv.f"
}
+/*< elseif (abs(qpar) .gt. EPS3) then >*/
+#line 598 "t_conv.f"
} else if (abs(qpar) > .001) {
+/*< if (qpar .ge. 0.0) phi = 0.0 >*/
+#line 599 "t_conv.f"
if (qpar >= 0.f) {
+#line 599 "t_conv.f"
phi = 0.f;
+#line 599 "t_conv.f"
}
+/*< if (qpar .lt. 0.0) phi = PI >*/
+#line 600 "t_conv.f"
if (qpar < 0.f) {
+#line 600 "t_conv.f"
phi = 3.1415926535897932384626433832795;
+#line 600 "t_conv.f"
}
+/*< elseif (abs(qperp) .gt. EPS3) then >*/
+#line 601 "t_conv.f"
} else if (abs(qperp) > .001) {
+/*< if (qperp .ge. 0.0) phi = 0.5 * PI >*/
+#line 602 "t_conv.f"
if (qperp >= 0.f) {
+#line 602 "t_conv.f"
phi = 1.5707963267948966;
+#line 602 "t_conv.f"
}
+/*< if (qperp .lt. 0.0) phi = -0.5 * PI >*/
+#line 603 "t_conv.f"
if (qperp < 0.f) {
+#line 603 "t_conv.f"
phi = -1.5707963267948966;
+#line 603 "t_conv.f"
}
+/*< else >*/
+#line 604 "t_conv.f"
} else {
+/*< phi = 0.0 >*/
+#line 605 "t_conv.f"
phi = 0.f;
+/*< endif >*/
+#line 606 "t_conv.f"
}
+/*< hrad = sqrt (hx**2 + hy**2) >*/
/* Computing 2nd power */
+#line 608 "t_conv.f"
r__1 = *hx;
/* Computing 2nd power */
+#line 608 "t_conv.f"
r__2 = *hy;
+#line 608 "t_conv.f"
hrad = sqrt(r__1 * r__1 + r__2 * r__2);
+/*< if (abs(hx) .gt. EPS3 .and. abs(hy) .gt. EPS3) then >*/
+#line 609 "t_conv.f"
if (dabs(*hx) > .001 && dabs(*hy) > .001) {
+/*< hdir = atan2 (abs(hy), abs(hx)) >*/
+#line 610 "t_conv.f"
hdir = atan2((dabs(*hy)), (dabs(*hx)));
+/*< if (hx .gt. 0 .and. hy .lt. 0) then >*/
+#line 611 "t_conv.f"
if (*hx > 0.f && *hy < 0.f) {
+/*< hdir = -hdir >*/
+#line 612 "t_conv.f"
hdir = -hdir;
+/*< elseif (hx .lt. 0 .and. hy .gt. 0) then >*/
+#line 613 "t_conv.f"
} else if (*hx < 0.f && *hy > 0.f) {
+/*< hdir = PI - hdir >*/
+#line 614 "t_conv.f"
hdir = 3.1415926535897932384626433832795 - hdir;
+/*< elseif (hx .lt. 0 .and. hy .lt. 0) then >*/
+#line 615 "t_conv.f"
} else if (*hx < 0.f && *hy < 0.f) {
+/*< hdir = hdir - PI >*/
+#line 616 "t_conv.f"
hdir += -3.1415926535897932384626433832795;
+/*< endif >*/
+#line 617 "t_conv.f"
}
+/*< elseif (abs(hx) .gt. EPS3) then >*/
+#line 618 "t_conv.f"
} else if (dabs(*hx) > .001) {
+/*< if (hx .ge. 0.0) hdir = 0.0 >*/
+#line 619 "t_conv.f"
if (*hx >= 0.f) {
+#line 619 "t_conv.f"
hdir = 0.f;
+#line 619 "t_conv.f"
}
+/*< if (hx .lt. 0.0) hdir = PI >*/
+#line 620 "t_conv.f"
if (*hx < 0.f) {
+#line 620 "t_conv.f"
hdir = 3.1415926535897932384626433832795;
+#line 620 "t_conv.f"
}
+/*< elseif (abs(hy) .gt. EPS3) then >*/
+#line 621 "t_conv.f"
} else if (dabs(*hy) > .001) {
+/*< if (hy .ge. 0.0) hdir = 0.5 * PI >*/
+#line 622 "t_conv.f"
if (*hy >= 0.f) {
+#line 622 "t_conv.f"
hdir = 1.5707963267948966;
+#line 622 "t_conv.f"
}
+/*< if (hy .lt. 0.0) hdir = -0.5 * PI >*/
+#line 623 "t_conv.f"
if (*hy < 0.f) {
+#line 623 "t_conv.f"
hdir = -1.5707963267948966;
+#line 623 "t_conv.f"
}
+/*< else >*/
+#line 624 "t_conv.f"
} else {
+/*< hdir = 0.0 >*/
+#line 625 "t_conv.f"
hdir = 0.f;
+/*< endif >*/
+#line 626 "t_conv.f"
}
+/*< hdir2 = hdir + phi - (helm/RD) >*/
+#line 628 "t_conv.f"
hdir2 = hdir + phi - *helm / 57.29577951308232087679815481410517;
/* ----------------------------------------------------------------------- */
/* !CALC CURRENTS */
/* !POSITION OF PSP FOR COIL I */
+/*< if (.not. at_sinq) then >*/
+#line 633 "t_conv.f"
if (! at_sinq__) {
+/*< hdir2 = hdir2 + 0.5 * PI ! ??? >*/
+#line 634 "t_conv.f"
hdir2 += 1.5707963267948966;
-/* ??? */
+/*< do ic = 1,3 >*/
+#line 635 "t_conv.f"
for (ic = 1; ic <= 3; ++ic) {
+/*< t_ih(ic+4) = cos(hdir2+(ic-1)*2.*PI/3.)*hrad/c_ih(ic)/1.5 >*/
+#line 636 "t_conv.f"
t_ih__[ic + 4] = cos(hdir2 + (ic - 1) * 2.f *
3.1415926535897932384626433832795 / 3.f) * hrad / c_ih__[
ic] / 1.5f;
+/*< enddo >*/
+#line 637 "t_conv.f"
}
+/*< t_ih(8) = hz/c_ih(4) >*/
+#line 638 "t_conv.f"
t_ih__[8] = *hz / c_ih__[4];
+/*< else >*/
+#line 639 "t_conv.f"
} else {
+/*< t_ih(5) = cos(hdir2) * hrad/c_ih(1) >*/
+#line 640 "t_conv.f"
t_ih__[5] = cos(hdir2) * hrad / c_ih__[1];
+/*< t_ih(6) = sin(hdir2) * hrad/c_ih(2) >*/
+#line 641 "t_conv.f"
t_ih__[6] = sin(hdir2) * hrad / c_ih__[2];
+/*< t_ih(7) = hz/c_ih(3) >*/
+#line 642 "t_conv.f"
t_ih__[7] = *hz / c_ih__[3];
+/*< endif >*/
+#line 643 "t_conv.f"
}
/* ----------------------------------------------------------------------- */
+/*< 999 CONTINUE >*/
+#line 645 "t_conv.f"
L999:
+/*< RETURN >*/
+#line 646 "t_conv.f"
return 0;
+/*< END >*/
} /* helm_case__ */
+/*< SUBRO >*/
/* Subroutine */ int flip_case__(integer *if1, integer *if2, real *t_ih__,
real *f1v, real *f1h, real *f2v, real *f2h, real *aki, real *akf,
integer *ier)
@@ -798,33 +1509,76 @@ L999:
/* CALCULATE P_IF CURRENTS FOR THE TWO FLIPPERS */
/* ----------------------------------------------------------------------- */
/* Define the dummy arguments */
-/* Part of T_CONV.FOR */
+/*< integer*4 if1, if2 >*/
+/*< real*4 t_ih(8) >*/
+/*< real*4 f1v, f1h, f2v, f2h >*/
+/*< real*4 aki, akf >*/
+/*< integer*4 ier >*/
/* ----------------------------------------------------------------------- */
/* INIT AND TEST */
+/*< IER = 0 >*/
+#line 665 "t_conv.f"
/* Parameter adjustments */
+#line 665 "t_conv.f"
--t_ih__;
+#line 665 "t_conv.f"
+#line 665 "t_conv.f"
/* Function Body */
+#line 665 "t_conv.f"
*ier = 0;
/* ----------------------------------------------------------------------- */
+/*< IF (IF1 .EQ. 1) THEN >*/
+#line 668 "t_conv.f"
if (*if1 == 1) {
+/*< T_IH(1) = F1V >*/
+#line 669 "t_conv.f"
t_ih__[1] = *f1v;
+/*< T_IH(2) = AKI*F1H >*/
+#line 670 "t_conv.f"
t_ih__[2] = *aki * *f1h;
+/*< ELSE >*/
+#line 671 "t_conv.f"
} else {
+/*< T_IH(1) = 0. >*/
+#line 672 "t_conv.f"
t_ih__[1] = 0.f;
+/*< T_IH(2) = 0. >*/
+#line 673 "t_conv.f"
t_ih__[2] = 0.f;
+/*< ENDIF >*/
+#line 674 "t_conv.f"
}
+/*< IF (IF2 .EQ. 1) THEN >*/
+#line 675 "t_conv.f"
if (*if2 == 1) {
+/*< T_IH(3) = F2V >*/
+#line 676 "t_conv.f"
t_ih__[3] = *f2v;
+/*< T_IH(4) = AKF*F2H >*/
+#line 677 "t_conv.f"
t_ih__[4] = *akf * *f2h;
+/*< ELSE >*/
+#line 678 "t_conv.f"
} else {
+/*< T_IH(3) = 0. >*/
+#line 679 "t_conv.f"
t_ih__[3] = 0.f;
+/*< T_IH(4) = 0. >*/
+#line 680 "t_conv.f"
t_ih__[4] = 0.f;
+/*< ENDIF >*/
+#line 681 "t_conv.f"
}
/* ----------------------------------------------------------------------- */
+/*< 999 CONTINUE >*/
+#line 683 "t_conv.f"
/* L999: */
+/*< RETURN >*/
+#line 684 "t_conv.f"
return 0;
+/*< END >*/
} /* flip_case__ */
diff --git a/t_conv.f b/t_conv.f
index 827a608e..cd0e0287 100755
--- a/t_conv.f
+++ b/t_conv.f
@@ -2,6 +2,11 @@ C------------------------------------------------------------------------
C slightly edited version for inclusion into SICS
C
C Mark Koennecke, November 2000
+C
+C Found that ERRESO looks error messages up in a 2D array. Modified IER
+C values to refer to a 1D array.
+C
+C Mark Koennecke, January 2002
C-------------------------------------------------------------------------
SUBROUTINE INICURVE(MIDX, MRX1, MRX2, AIDX, ARX1, ARX2,
+ MMIN, MMAX, AMIN, AMAX)
@@ -165,12 +170,12 @@ C
IF (LDK(2*IFX-1) .OR. LDK(2*IFX)) THEN
LMOAN(IFX) = .TRUE.
IF (LDK(2*IFX-1)) THEN
- IER = 1
+ IER = 1 + 8
IF(EDEF(1) .LT. EPS1) GOTO 999
IER = 0
AKDEF(1) = SQRT(EDEF(1)/F)
ELSE
- IER = 1
+ IER = 1 + 8
IF(AKDEF(1) .LT. EPS1) GOTO 999
IER = 0
EDEF(1) = F*AKDEF(1)**2
@@ -193,12 +198,12 @@ C
IF (LDK(5-2*IFX) .OR. LDK(6-2*IFX)) THEN
LMOAN(3-IFX) = .TRUE.
IF (LDK(5-2*IFX)) THEN
- IER = 1
+ IER = 1 + 8
IF(EDEF(2) .LT. EPS4) GOTO 999
IER = 0
AKDEF(2) = SQRT(EDEF(2)/F)
ELSE
- IER = 1
+ IER = 1 + 8
IF(AKDEF(2) .LT. EPS4) GOTO 999
IER = 0
EDEF(2) = F*AKDEF(2)**2
@@ -207,7 +212,7 @@ C
ELSEIF (LQHKLE) THEN
LMOAN(3-IFX) = .TRUE.
EDEF(2) = EDEF(1)+(2*IFX-3)*EN
- IER = 1
+ IER = 1 + 8
IF(EDEF(2) .LT. EPS4) GOTO 999
IER = 0
AKDEF(2) = SQRT(EDEF(2)/F)
@@ -235,6 +240,7 @@ C
LDR_ALM = .TRUE.
endif
ELSE
+ IER = IER + 8
GOTO 999
ENDIF
ENDIF
@@ -254,6 +260,7 @@ C
LDR_RA = .TRUE.
endif
ELSE
+ IER = IER + 8
GOTO 999
ENDIF
ENDIF
@@ -276,6 +283,7 @@ C
LDRA(4) = .TRUE.
QM = DQM
ELSE
+ IER = IER + 4
GOTO 999
ENDIF
ENDIF
diff --git a/t_update.c b/t_update.c
index 11c42d3d..7aef7fc1 100644
--- a/t_update.c
+++ b/t_update.c
@@ -4,7 +4,6 @@
*/
#include "f2c.h"
-#include
/* Subroutine */ int t_update__(real *p_a__, real *p_ih__, real *c_ih__,
logical *lpa, real *dm, real *da, integer *isa, real *helm, real *f1h,
@@ -116,7 +115,7 @@
/* ----------------------------------------------------------------------- */
ieri = 0;
- ieru = 1;
+ ieru = 9;
ex_up__(&ddm, &dei, &daki, &da2, &df, &ieri);
if (ieri == 0) {
*ei = dei;
@@ -125,9 +124,10 @@
} else {
imod = 3;
erreso_(&imod, &ieri);
+ ieru = ieri + 8;
}
ieri = 0;
- ieru = 1;
+ ieru = 9;
ex_up__(&dda, &def, &dakf, &da6, &df, &ieri);
if (ieri == 0) {
*ef = def;
@@ -143,6 +143,7 @@
} else {
imod = 3;
erreso_(&imod, &ieri);
+ ieru = ieri + 8;
}
}
if (*isa == 0) {
@@ -156,7 +157,7 @@
*en = dei - def;
}
ieri = 0;
- ieru = 1;
+ ieru = 5;
sam_up__(dbqhkl, &dqm, &dqs, &dphi, &daki, &dakf, &da3, &da4, &ieri);
if (ieri == 0) {
for (id = 1; id <= 3; ++id) {
@@ -168,6 +169,7 @@
} else {
imod = 2;
erreso_(&imod, &ieri);
+ ieru = ieri + 4;
}
ieri = 0;
@@ -222,6 +224,9 @@
/* ----------------------------------------------------------------------- */
/* Error - IER=1 if DX OR AX TOO SMALL */
/* ----------------------------------------------------------------------- */
+/* !!!!!!!!!! This has to be fixed manually after conversion by f2c. */
+/* !!!!!!!!!! The reason is a different definition of the abs function. */
+/* !!!!!!!!!! MK, May 2001 */
arg = *dx * sin(*ax2 / 114.59155902616465);
if(arg < .0)
arg = -arg;
diff --git a/t_update.f b/t_update.f
index 0ff78f9d..f5c63d81 100755
--- a/t_update.f
+++ b/t_update.f
@@ -89,7 +89,7 @@ C
C-----------------------------------------------------------------------
C
IERI=0
- IERU=1
+ IERU=1 + 8
CALL EX_UP(DDM,DEI,DAKI,DA2,DF,IERI)
IF (IERI.EQ.0) THEN
EI=DEI
@@ -98,9 +98,10 @@ C
ELSE
IMOD=3
CALL ERRESO(IMOD,IERI)
+ IERU = IERI + 8
ENDIF
IERI=0
- IERU=1
+ IERU=1 + 8
CALL EX_UP(DDA,DEF,DAKF,DA6,DF,IERI)
IF (IERI.EQ.0) THEN
EF=DEF
@@ -116,6 +117,7 @@ C
ELSE
IMOD=3
CALL ERRESO(IMOD,IERI)
+ IERU = 8 + IERI
ENDIF
ENDIF
IF (ISA.EQ.0) THEN
@@ -127,7 +129,7 @@ C
ENDIF
IF (IERU.EQ.0) EN=DEI-DEF
IERI=0
- IERU=1
+ IERU=1 + 4
CALL SAM_UP(DBQHKL,DQM,DQS,DPHI,DAKI,DAKF,DA3,DA4,IERI)
IF (IERI.EQ.0) THEN
DO ID=1,3
@@ -139,6 +141,7 @@ C
ELSE
IMOD=2
CALL ERRESO(IMOD,IERI)
+ IERU = IERI + 4
ENDIF
C
IERI=0
diff --git a/tas.h b/tas.h
index 535a1266..b4c02d51 100644
--- a/tas.h
+++ b/tas.h
@@ -107,8 +107,21 @@
#define LOC 91
#define SWUNIT 92
#define SINFO 93
+#define TEI 94
+#define TKI 95
+#define TEF 96
+#define TKF 97
+#define TQH 98
+#define TQK 99
+#define TQL 100
+#define TEN 101
+#define TQM 102
+#define HX 34
+#define HY 35
+#define HZ 36
-#define MAXPAR 94
+
+#define MAXPAR 103
#define MAXADD 20
#define MAXEVAR 10
diff --git a/tascom.tcl b/tascom.tcl
index 4d858e55..df8f0219 100644
--- a/tascom.tcl
+++ b/tascom.tcl
@@ -148,7 +148,7 @@ proc scatSense {par {val -1000} } {
if { $val == -1000 } {
return [eval $par]
}
- if {$val != 1 && $val != -1 } {
+ if {$val != 1 && $val != -1 && $val != 0 } {
error "ERROR: invalid scattering sense $val"
}
switch $par {
@@ -162,14 +162,43 @@ proc scatSense {par {val -1000} } {
}
sa {
set oldzero [tasSplit [madZero $mot]]
- set newZero [expr $val*180 + $oldzero]
- madZero $mot $newZero
- a5 softlowerlim $newZero
- $par $val
+ set oldupper [tasSplit [$mot softupperlim]]
+ set oldlower [tasSplit [$mot softlowerlim]]
+ set oldsa [tasSplit [sa]]
+ if { $val == 0 && $oldsa == 1} {
+ set newzero [expr $oldzero - 90.]
+ set newlower [expr $oldlower - 90.]
+ set newupper [expr $oldupper - 90.]
+ } elseif {$val == 0 && $oldsa == -1} {
+ set newzero [expr $oldzero + 90.]
+ set newlower [expr $oldlower + 90.]
+ set newupper [expr $oldupper + 90.]
+ } elseif { $val == 1 && $oldsa == 0} {
+ set newzero [expr $oldzero + 90.]
+ set newlower [expr $oldlower + 90.]
+ set newupper [expr $oldupper + 90.]
+ } elseif { $val == -1 && $oldsa == 0} {
+ set newzero [expr $oldzero - 90.]
+ set newlower [expr $oldlower - 90.]
+ set newupper [expr $oldupper - 90.]
+ } elseif { $val == 1 && $oldsa == -1} {
+ set newzero [expr $oldzero + 180. ]
+ set newlower [expr $oldlower + 180. ]
+ set newupper [expr $oldupper + 180. ]
+ } elseif {$val == -1 && $oldsa == 1} {
+ set newzero [expr $oldzero - 180. ]
+ set newlower [expr $oldlower - 180. ]
+ set newupper [expr $oldupper - 180. ]
+ } else {
+ error "Unknown SA setting combination"
+ }
+ $par $val
+ madZero $mot $newzero
+ $mot softupperlim $newupper
+ $mot softlowerlim $newlower
}
}
}
-
#--------------------------------------------------------------------------
# The output command
@@ -190,7 +219,7 @@ proc ou args {
}
#--------------------------------------------------------------------------
-# typeATokenizer extracts tokens from acpmmand string. Tokens can be
+# typeATokenizer extracts tokens from a command string. Tokens can be
# either variable names or - indicating a series of variables.
# Returns the token value or END if the end of the string text is
# reached. Uses and updates a variable pos which indicates the current
@@ -431,7 +460,6 @@ proc varSet { command } {
} else {
clientput [format " %s = %s" $token $value]
}
-
} else {
set ret [catch {eval $token $value} msg]
if { $ret != 0 } {
@@ -443,8 +471,8 @@ proc varSet { command } {
set token [varToken $command $pos]
set value [varToken $command $pos]
}
+ catch {updateqe} msg
}
-
#--------------------------------------------------------------------------
# co for count is the funny MAD count procedure. Please note, that the
# count mode is automatically set through the last MN or TI variable.
@@ -570,7 +598,7 @@ proc se args {
#------- is it the only command line case?
if {[llength $args] > 0 } {
set line [join $args]
- return [varSet $line]
+ return [varSet $line]
} else {
#------- we are prompting
while { 1== 1} {
@@ -728,14 +756,26 @@ proc le args {
set v7 [tasSplit [ql]]
set val [format " %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f \n" \
$v1 $v2 $v3 $v4 $v5 $v6 $v7]
+ set v1 [tasSplit [tei]]
+ set v2 [tasSplit [tki]]
+ set v3 [tasSplit [tef]]
+ set v4 [tasSplit [tkf]]
+ set v5 [tasSplit [tqh]]
+ set v6 [tasSplit [tqk]]
+ set v7 [tasSplit [tql]]
+ set val2 [format " %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f \n" \
+ $v1 $v2 $v3 $v4 $v5 $v6 $v7]
append output [format "POSN: %s" $val]
- append output [format "TARG: %s" $val]
+ append output [format "TARG: %s" $val2]
append output [format " EN QM\n"]
set v1 [tasSplit [en]]
set v2 [tasSplit [qm]]
set val [format " %9.4f %9.4f\n" $v1 $v2]
+ set v1 [tasSplit [ten]]
+ set v2 [tasSplit [qm]]
+ set val2 [format " %9.4f %9.4f\n" $v1 $v2]
append output [format "POSN: %s" $val]
- append output [format "TARG: %s" $val]
+ append output [format "TARG: %s" $val2]
return $output
}
@@ -885,6 +925,7 @@ proc sz args {
#-------action
set newZero [expr $zero + ($val - $pos)]
madZero $mot $newZero
+ catch {updateqe} msg
#-------- more output
set zero [tasSplit [madZero $mot]]
set loh [tasSplit [eval $mot hardlowerlim]]
diff --git a/tasdrive.c b/tasdrive.c
index 44a54f00..cb5963ef 100644
--- a/tasdrive.c
+++ b/tasdrive.c
@@ -148,6 +148,7 @@ int TASDrive(SConnection *pCon, SicsInterp *pSics, void *pData,
tasMask[varPointer] = 1;
oldEnergy[varPointer] = self->tasPar[EMIN+varPointer]->fVal;
self->tasPar[EMIN + varPointer]->fVal = atof(pToken);
+ self->tasPar[ETARGET + varPointer]->fVal = atof(pToken);
}
else
{
diff --git a/tasinit.c b/tasinit.c
index 987b6043..5b295462 100644
--- a/tasinit.c
+++ b/tasinit.c
@@ -150,6 +150,15 @@ extern char *tasVariableOrder[] = {
"local",
"swunit",
"scaninfo",
+ "tei",
+ "tki",
+ "tef",
+ "tkf",
+ "tqh",
+ "tqk",
+ "tql",
+ "ten",
+ "tqm",
NULL};
/*---------------------------------------------------------------------
There is a special feauture in MAD where the count mode is determined
@@ -183,6 +192,25 @@ static int TimerCallback(int iEvent, void *pEvent, void *pUser)
SetCounterPreset(self->pScan->pCounterData,self->tasPar[TI]->fVal);
return 1;
}
+/*-----------------------------------------------------------------------
+ This is an interpreter wrapper function which allows to call for the
+ recalculation of the energy variables from scripts.
+--------------------------------------------------------------------------*/
+extern int TASUpdate(pTASdata self,SConnection *pCon); /* tasutil.c */
+
+static int RecalcAction(SConnection *pCon, SicsInterp *pSics, void *pData,
+ int argc, char *argv[])
+{
+ pTASdata self = NULL;
+
+ assert(pCon);
+ assert(pSics);
+ self = (pTASdata)pData;
+ assert(self);
+
+ return TASUpdate(self,pCon);
+}
+
/*-----------------------------------------------------------------------
A function for killing the TAS data structure is needed
-------------------------------------------------------------------------*/
@@ -293,10 +321,17 @@ int TASFactory(SConnection *pCon, SicsInterp *pSics, void *pData,
TASKill(pNew);
return 0;
}
- iError = AddCommand(pSics,"sf",TASScan,NULL,pNew);
+ iError = AddCommand(pSics,"fs",TASScan,NULL,pNew);
if(!iError)
{
- SCWrite(pCon,"ERROR: duplicate set command not created",eError);
+ SCWrite(pCon,"ERROR: duplicate sf command not created",eError);
+ TASKill(pNew);
+ return 0;
+ }
+ iError = AddCommand(pSics,"updateqe",RecalcAction,NULL,pNew);
+ if(!iError)
+ {
+ SCWrite(pCon,"ERROR: duplicate updateqe command not created",eError);
TASKill(pNew);
return 0;
}
diff --git a/tasscan.c b/tasscan.c
index ebc76225..597ccc27 100644
--- a/tasscan.c
+++ b/tasscan.c
@@ -552,6 +552,8 @@ static int TASScanDrive(pScanData self, int iPoint)
iTAS = 1;
pTAS->tasPar[EI+iPtr]->fVal =
pVar->fStart + iPoint * pVar->fStep;
+ pTAS->tasPar[ETARGET+iPtr]->fVal =
+ pVar->fStart + iPoint * pVar->fStep;
tasMask[iPtr] = 1;
}
else
@@ -575,14 +577,16 @@ static int TASScanDrive(pScanData self, int iPoint)
{
status = TASCalc(pTAS,self->pCon,tasMask,
tasTargets, tasTargetMask);
- if(!status)
- return 0;
- TASStart(pTAS,self->pCon,
+ if(status)
+ {
+ /*
+ Errors both in calculation or in starting motors are
+ ignored here on purpose. There is a slight chance that
+ other points in the scan fit the bill.
+ */
+ TASStart(pTAS,self->pCon,
self->pSics,tasTargets,tasTargetMask);
- /*
- again ignore errors, the scan shall continue if an error was
- found
- */
+ }
}
/*
diff --git a/tastest.tcl b/tastest.tcl
index 40e2b3d1..d161cbba 100644
--- a/tastest.tcl
+++ b/tastest.tcl
@@ -57,8 +57,8 @@ TokenInit connan
SicsUser Spy 007 1
#---------------------------------------------------------------------------
# M O T O R S
-Motor A1 SIM 0. 111. -.1 2. # Monochromator Theta
-Motor A2 SIM 33.1 120. -.1 2. # Monochromator Two-Theta
+Motor A1 SIM -87. 6.1 -.1 2. # Monochromator Theta
+Motor A2 SIM -129.1 -22. -.1 2. # Monochromator Two-Theta
Motor A3 SIM -177.3 177.3 -.1 2. # Sample theta or omega
Motor A4 SIM -135.1 123.4 -.1 2. # Sample Two-Theta
Motor A5 SIM -200 200 -.1 2. # Analyzer Theta
@@ -122,6 +122,18 @@ VarMake QH Float User
VarMake QK Float User
VarMake QL Float User
VarMake EN Float User
+
+#-------- energy Q targets
+VarMake TEI Float User
+VarMake TKI Float User
+VarMake TEF Float User
+VarMake TKF Float User
+VarMake TQH Float User
+VarMake TQK Float User
+VarMake TQL Float User
+VarMake TEN Float User
+VarMake TQM Float User
+
#---------------------------------------------------------------------------
# I N S T R U M E N T V A R I A B L E S
# DM, DA d-spacing monochromator, analyzer
@@ -142,7 +154,7 @@ instrument lock
VarMake DM Float Mugger
VarMake DA Float Mugger
VarMake SM Int User
-SM 1
+SM -1
SM lock
VarMake SS Int User
VarMake SA Int User
diff --git a/tasu.h b/tasu.h
index 3bf1eaab..cbf5cf6a 100644
--- a/tasu.h
+++ b/tasu.h
@@ -21,11 +21,11 @@ extern char *tasVariableOrder[];
/*
Note: the defines below MUST map the range between EI - HZ in the list
of variables as defined in tas.h. Otherwise quite interesting things
- can happen.
+ can happen. ETARGET is the variable order index for the energy targets.
*/
#define EMIN 25
#define EMAX 36
-
+#define ETARGET 94
int isTASMotor(char *val);
int isTASVar(char *val);
diff --git a/tasutil.c b/tasutil.c
index 190fa11f..e98aa2ec 100644
--- a/tasutil.c
+++ b/tasutil.c
@@ -131,7 +131,66 @@ void prepare2Parse(char *line)
}
return 1;
}
-
+/*--------------------------------------------------------------------
+ print calculation errors.
+*/
+static int printError(int ier, SConnection *pCon)
+{
+ /*
+ error messages taken from erreso
+ */
+ switch(ier)
+ {
+ case 1:
+ SCWrite(pCon,"ERROR: Bad lattice parameters(AS,BS,CS)",eError);
+ return 0;
+ break;
+ case 2:
+ SCWrite(pCon,"ERROR: Bad cell angles",eError);
+ return 0;
+ break;
+ case 3:
+ SCWrite(pCon,"ERROR: Bad scattering plane ",eError);
+ return 0;
+ break;
+ case 4:
+ SCWrite(pCon,"ERROR: Bad lattice or lattice plane",eError);
+ return 0;
+ break;
+ case 5:
+ SCWrite(pCon,"ERROR: Check Lattice and Scattering Plane",eError);
+ return 0;
+ break;
+ case 6:
+ SCWrite(pCon,"ERROR: Q not in scattering plane",eError);
+ return 0;
+ break;
+ break;
+ case 7:
+ SCWrite(pCon,"ERROR: Q-modulus to small",eError);
+ return 0;
+ break;
+ case 8:
+ case 12:
+ SCWrite(pCon,"ERROR: KI,KF,Q triangle cannot close",eError);
+ return 0;
+ break;
+ case 9:
+ SCWrite(pCon,"ERROR: KI or K, check d-spacings",eError);
+ return 0;
+ break;
+ case 10:
+ SCWrite(pCon,"ERROR: KI or KF cannot be obtained",eError);
+ return 0;
+ break;
+ case 11:
+ SCWrite(pCon,"ERROR: KI or KF to small",eError);
+ return 0;
+ break;
+ default:
+ break;
+ }
+}
/*----------------------------------------------------------------------
TASCalc does the triple axis spectrometer calculations. This function is
invoked whenever energy or Q needs to be driven. The parameters:
@@ -266,14 +325,14 @@ int TASCalc(pTASdata self, SConnection *pCon,
This done, we painstackingly initialize the tons of parameters required
by the t_conv
*/
- ei = (real)self->tasPar[EI]->fVal;
- aki = (real)self->tasPar[KI]->fVal;
- ef = (real)self->tasPar[EF]->fVal;
- akf = (real)self->tasPar[KF]->fVal;
- qhkl[0] = (real)self->tasPar[QH]->fVal;
- qhkl[1] = (real)self->tasPar[QK]->fVal;
- qhkl[2] = (real)self->tasPar[QL]->fVal;
- en = (real)self->tasPar[EN]->fVal;
+ ei = (real)self->tasPar[TEI]->fVal;
+ aki = (real)self->tasPar[TKI]->fVal;
+ ef = (real)self->tasPar[TEF]->fVal;
+ akf = (real)self->tasPar[TKF]->fVal;
+ qhkl[0] = (real)self->tasPar[TQH]->fVal;
+ qhkl[1] = (real)self->tasPar[TQK]->fVal;
+ qhkl[2] = (real)self->tasPar[TQL]->fVal;
+ en = (real)self->tasPar[TEN]->fVal;
hx = (real)self->tasPar[HX]->fVal;
hy = (real)self->tasPar[HY]->fVal;
hz = (real)self->tasPar[HZ]->fVal;
@@ -288,6 +347,9 @@ int TASCalc(pTASdata self, SConnection *pCon,
ldf = (logical)tasMask[9];
if(self->tasPar[LPA]->iVal > 0)
lpa = (logical)1;
+ else
+ lpa = (logical)0;
+
dm = (real)self->tasPar[DM]->fVal;
da = (real)self->tasPar[DA]->fVal;
qm = (real)self->tasPar[QM]->fVal;
@@ -323,7 +385,7 @@ int TASCalc(pTASdata self, SConnection *pCon,
{
ldra[i] = 0;
}
- l_RA = l_RM = l_ALM = 0;
+ l_RA = l_RM = l_ALM = ier = 0;
/* now we can call */
t_conv__(&ei, &aki, &ef, &akf,
@@ -334,58 +396,11 @@ int TASCalc(pTASdata self, SConnection *pCon,
angles, &tRM, &tALM, &tRA, &qm, ldra,
&l_RM, &l_ALM,&l_RA, currents,helmconv,
&ier);
- /*
- error messages taken from erreso
- */
- switch(ier)
+
+ if(ier != 0)
{
- case 1:
- SCWrite(pCon,"ERROR: Bad lattice parameters",eError);
- return 0;
- break;
- case 2:
- SCWrite(pCon,"ERROR: Bad cell angles",eError);
- return 0;
- break;
- case 3:
- SCWrite(pCon,"ERROR: Bad scattering plane ",eError);
- return 0;
- break;
- case 4:
- SCWrite(pCon,"ERROR: Bad lattice or lattice plane",eError);
- return 0;
- break;
- case 5:
- SCWrite(pCon,"ERROR: Q not in scattering plane",eError);
- return 0;
- break;
- case 6:
- SCWrite(pCon,"ERROR: Q modulus to small",eError);
- return 0;
- break;
- case 7:
- case 12:
- SCWrite(pCon,"ERROR: KI, KF, Q triangle cannot close ",eError);
- return 0;
- break;
- case 8:
- SCWrite(pCon,"ERROR: in KI, KF, check d-spacings",eError);
- return 0;
- break;
- case 9:
- SCWrite(pCon,"ERROR: KI or KF cannot be obtained",eError);
- return 0;
- break;
- case 10:
- SCWrite(pCon,"ERROR: KI or KF to small",eError);
- return 0;
- break;
- case 11:
- SCWrite(pCon,"ERROR: KI or KF cannot be obtained",eError);
- return 0;
- break;
- default:
- break;
+ printError(ier,pCon);
+ return 0;
}
/*
@@ -409,12 +424,12 @@ int TASCalc(pTASdata self, SConnection *pCon,
{
motorTargets[17+i] = helmconv[i];
}
- self->tasPar[EI]->fVal = (float)ei;
- self->tasPar[KI]->fVal = (float)aki;
- self->tasPar[EF]->fVal = (float)ef;
- self->tasPar[KF]->fVal = (float)akf;
- self->tasPar[EN]->fVal = (float)en;
-
+ self->tasPar[TEI]->fVal = (float)ei;
+ self->tasPar[TKI]->fVal = (float)aki;
+ self->tasPar[TEF]->fVal = (float)ef;
+ self->tasPar[TKF]->fVal = (float)akf;
+ self->tasPar[TEN]->fVal = (float)en;
+ self->tasPar[QM]->fVal = (float)qm;
return 1;
}
@@ -578,6 +593,9 @@ int TASUpdate(pTASdata self, SConnection *pCon)
*/
if(self->tasPar[LPA]->iVal > 0)
lpa = (logical)1;
+ else
+ lpa = (logical)0;
+
da = (real)self->tasPar[DA]->fVal;
dm = (real)self->tasPar[DM]->fVal;
isa = (integer)self->tasPar[SA]->iVal;
@@ -593,76 +611,28 @@ int TASUpdate(pTASdata self, SConnection *pCon)
/*
now call t_update to do the calculation
*/
+ ier = 0;
t_update__(amot, helmCurrent, convH, &lpa, &dm, &da, &isa, &helm,
&f1h, &f1v, &f2h, &f2v, &f, &ei, &aki, &ef, &akf,
qhkl, &en, &hx, &hy, &hz, &if1, &if2, &qm, &ier);
/*
!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- t_update's function ex_up must be dited manually after conversion
+ t_update's function ex_up must be edited manually after conversion
to c by f2c. The reason is a different definition of the abs function.
The line:
arg = (d__1 = *dx * sin(*ax2 / 114.59155902616465, abs(d__1));
- has to be changed to:
+ has to be changed to:
arg = *dx * sin(*ax2 / 114.59155902616465);
if(arg < .0)
arg = -arg;
-
+
!!!!!!!!!!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!!!!!!!!!!!!!*/
-
- /*
- error messages taken from erreso
- */
- switch(ier)
+ if(ier != 0)
{
- case 1:
- SCWrite(pCon,"ERROR: Bad lattice parameters",eError);
- return 0;
- break;
- case 2:
- SCWrite(pCon,"ERROR: Bad cell angles",eError);
- return 0;
- break;
- case 3:
- SCWrite(pCon,"ERROR: Bad scattering plane ",eError);
- return 0;
- break;
- case 4:
- SCWrite(pCon,"ERROR: Bad lattice or lattice plane",eError);
- return 0;
- break;
- case 5:
- SCWrite(pCon,"ERROR: Q not in scattering plane",eError);
- return 0;
- break;
- case 6:
- SCWrite(pCon,"ERROR: Q modulus to small",eError);
- return 0;
- break;
- case 7:
- case 12:
- SCWrite(pCon,"ERROR: KI, KF, Q triangle cannot close ",eError);
- return 0;
- break;
- case 8:
- SCWrite(pCon,"ERROR: in KI, KF, check d-spacings",eError);
- return 0;
- break;
- case 9:
- SCWrite(pCon,"ERROR: KI or KF cannot be obtained",eError);
- return 0;
- break;
- case 10:
- SCWrite(pCon,"ERROR: KI or KF to small",eError);
- return 0;
- break;
- case 11:
- SCWrite(pCon,"ERROR: KI or KF cannot be obtained",eError);
- return 0;
- break;
- default:
- break;
+ printError(ier,pCon);
+ return 0;
}
/*
@@ -711,3 +681,4 @@ int TASUpdate(pTASdata self, SConnection *pCon)
return 1;
}
+
diff --git a/trics.dic b/trics.dic
index 2d2a794b..9b29ad5c 100644
--- a/trics.dic
+++ b/trics.dic
@@ -49,7 +49,8 @@ framepreset = \
/$(framename),NXentry/TRICS,NXinstrument/count_control,NXcounter/SDS preset
framemode = /$(framename),NXentry/TRICS,NXinstrument/count_control,NXcounter/SDS countmode \
-type DFNT_UINT8 -rank 1 -dim {132}
-framemonitor = /$(framename),NXentry/TRICS,NXinstrument/count_control,NXcounter/SDS monitor \
+framemonitor = /$(framename),NXentry/TRICS,NXinstrument/count_control,NXcounter/SDS monitor -type DFNT_INT32
+sinqmonitor= /$(framename),NXentry/TRICS,NXinstrument/count_control,NXcounter/SDS beam_monitor \
-type DFNT_INT32
#------------------------ Detector
dnumber = detector1
@@ -72,7 +73,7 @@ frametilt = /$(framename),NXentry/TRICS,NXinstrument/$(dnumber),NXdetector/SDS
-attr {units,degrees}
framecounts = /$(framename),NXentry/TRICS,NXinstrument/$(dnumber),NXdetector/SDS counts \
-attr {signal,1} -attr {units,counts} -type DFNT_INT32 \
- -LZW -chunk {256,256} -rank 2 -dim {$(framedim1),$(framedim2)}
+ -LZW -rank 2 -dim {$(framedim1),$(framedim2)}
detzerox = \
/frame0000,NXentry/TRICS,NXinstrument/$(dnumber),NXdetector/SDS x_zero_point \
-attr {units,pixel}