From beba0d4644631cf60530ccef7df21a97c5f1887f Mon Sep 17 00:00:00 2001 From: koennecke Date: Wed, 23 Mar 2005 08:19:47 +0000 Subject: [PATCH] - Initial commit of a UB calculation setup for four circle diffractometers --- cell.c | 150 +++++++++++++++++++++ cell.h | 46 +++++++ exeman.c | 3 +- macro.c | 4 +- make_gen | 6 +- ofac.c | 7 +- scan.c | 47 +------ stdscan.c | 65 +++++++++ stdscan.h | 4 + trigd.c | 92 +++++++++++++ trigd.h | 17 +++ ubcalc.c | 389 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ ubcalc.h | 19 +++ ubcalc.w | 28 ++++ ubfour.c | 136 +++++++++++++++++++ ubfour.h | 40 ++++++ vector.c | 145 ++++++++++++++++++++ vector.h | 93 +++++++++++++ 18 files changed, 1236 insertions(+), 55 deletions(-) create mode 100644 cell.c create mode 100644 cell.h create mode 100644 trigd.c create mode 100644 trigd.h create mode 100644 ubcalc.c create mode 100644 ubcalc.h create mode 100644 ubcalc.w create mode 100644 ubfour.c create mode 100644 ubfour.h create mode 100644 vector.c create mode 100644 vector.h diff --git a/cell.c b/cell.c new file mode 100644 index 00000000..9779d656 --- /dev/null +++ b/cell.c @@ -0,0 +1,150 @@ +/** + * this is a little library for performing crystallographic cell transformations + * for SICS. Some of the actual code was lifted from the Risoe program tascom. + * + * copyright: see file COPYRIGHT + * + * Mark Koennecke, March 2005 + */ +#include +#include +#include +#include "trigd.h" +#include "cell.h" + +/* define constants */ +#ifndef PI +#define PI (3.1415926536) /* pi */ +#endif +#define TWOPI (2*PI) /* 2*pi */ + +/******************************************************************************* +* Transform direct lattice to reciprocal lattice. +*******************************************************************************/ +int directToReciprocalLattice(lattice direct, plattice reciprocal) +{ + double alfa, beta, gamma; + double cos_alfa, cos_beta, cos_gamma; + double sin_alfa, sin_beta, sin_gamma; + double ad, bd, cd; + double arg, vol; + + alfa = direct.alpha; + beta = direct.beta; + gamma = direct.gamma; + + cos_alfa = Cosd (alfa); + cos_beta = Cosd (beta); + cos_gamma = Cosd (gamma); + + sin_alfa = Sind (alfa); + sin_beta = Sind (beta); + sin_gamma = Sind (gamma); + + reciprocal->alpha = Acosd ((cos_beta*cos_gamma - cos_alfa)/sin_beta/sin_gamma); + reciprocal->beta =Acosd ((cos_alfa*cos_gamma - cos_beta)/sin_alfa/sin_gamma); + reciprocal->gamma = Acosd ((cos_alfa*cos_beta - cos_gamma)/sin_alfa/sin_beta); + + ad = direct.a; + bd = direct.b; + cd = direct.c; + + arg = 1 + 2*cos_alfa*cos_beta*cos_gamma - cos_alfa*cos_alfa - + cos_beta*cos_beta - + cos_gamma*cos_gamma; + if (arg < 0.0) + { + return REC_NO_VOLUME; + } + + vol = ad*bd*cd*sqrt (arg); + reciprocal->a = bd*cd*sin_alfa/vol; + reciprocal->b = ad*cd*sin_beta/vol; + reciprocal->c = bd*ad*sin_gamma/vol; + + return (0); +} +/******************************************************************************* +* Transform reciprocal lattice to direct lattice. +*******************************************************************************/ +int reciprocalToDirectLattice(lattice reciprocal, plattice direct) +{ + double alfa, beta, gamma; + double cos_alfa, cos_beta, cos_gamma; + double sin_alfa, sin_beta, sin_gamma; + double ar, br, cr; + double arg, vol; + + alfa = reciprocal.alpha; + beta = reciprocal.beta; + gamma = reciprocal.gamma; + + cos_alfa = Cosd (alfa); + cos_beta = Cosd (beta); + cos_gamma = Cosd (gamma); + + sin_alfa = Sind (alfa); + sin_beta = Sind (beta); + sin_gamma = Sind (gamma); + + direct->alpha = Acosd ((cos_beta*cos_gamma - cos_alfa)/sin_beta/sin_gamma); + direct->beta = Acosd ((cos_alfa*cos_gamma - cos_beta)/sin_alfa/sin_gamma); + direct->gamma = Acosd ((cos_alfa*cos_beta - cos_gamma)/sin_alfa/sin_beta); + + ar = reciprocal.a; + br = reciprocal.b; + cr = reciprocal.c; + + arg = 1 + 2*cos_alfa*cos_beta*cos_gamma - cos_alfa*cos_alfa - + cos_beta*cos_beta - + cos_gamma*cos_gamma; + if (arg < 0.0) + { + return REC_NO_VOLUME; + } + + vol = ar*br*cr*sqrt (arg); + direct->a = br*cr*sin_alfa/vol; + direct->b = ar*cr*sin_beta/vol; + direct->c = br*ar*sin_gamma/vol; + + return (0); +} +/*************************************************************************************** + * Build a B matrix + ***************************************************************************************/ +int calculateBMatrix(lattice direct, MATRIX B) { + lattice reciprocal; + int status; + + assert(MatRow(B) == 3); + assert(MatCol(B) == 3); + + status = directToReciprocalLattice(direct,&reciprocal); + if(status < 0) { + return status; + } + + mat_fill(B,ZERO_MATRIX); + + /* + top row + */ + B[0][0] = reciprocal.a; + B[0][1] = reciprocal.b*Cosd(reciprocal.gamma); + B[0][2] = reciprocal.c*Cosd(reciprocal.beta); + + /* + middle row + */ + B[1][1] = reciprocal.b*Sind(reciprocal.gamma); + B[1][2] = -reciprocal.c*Sind(reciprocal.beta)*Cosd(direct.alpha); + + /* + bottom row + */ + B[2][2] = 1./direct.c; + + return 1; +} + diff --git a/cell.h b/cell.h new file mode 100644 index 00000000..8accfffc --- /dev/null +++ b/cell.h @@ -0,0 +1,46 @@ +/** + * this is a little library for performing crystallographic cell transformations + * for SICS. Some of the actual code was lifted from the Risoe program tascom. + * + * copyright: see file COPYRIGHT + * + * Mark Koennecke, March 2005 + */ +#ifndef SICSCELL +#define SICSCELL +#include "matrix/matrix.h" + +/** + * error codes + */ +#define REC_NO_VOLUME -100 + +/** + * lattice parameters: either reciprocal or direct + */ + typedef struct { + double a,b,c; + double alpha, beta, gamma; + }lattice, *plattice; +/** + * conversion from a direct lattice to the recipcrocal one. + * @param direct The input direct cell parameters. + * @param reciprocal The output reciprocal cell constants + * @return 0 on success, > 0 else + */ + int directToReciprocalLattice(lattice direct, plattice reciprocal); +/** + * conversion from a reciprocal lattice to the directone. + * @param reciprocal The input reciprocal cell parameters. + * @param direct The output direct cell constants + * @return 0 on success, > 0 else + */ + int reciprocalToDirectLattice(lattice reciprocal, plattice direct); +/** + * calculate a crystallographic B matrix from the cell constants + * @param direct The direct cell lattice to calculate B from + * @param B will be filled with the B matrix. MUST be 3x3 + * @return 1 on success, an negative error code else + */ +int calculateBMatrix(lattice direct, MATRIX B); +#endif diff --git a/exeman.c b/exeman.c index f26b4b91..5df8f048 100644 --- a/exeman.c +++ b/exeman.c @@ -721,8 +721,9 @@ int ExeManagerWrapper(SConnection *pCon, SicsInterp *pSics, void *pData, self = (pExeMan)pData; assert(self != NULL); - strncpy(pBufferName,argv[1],255); + if(argc > 1){ + strncpy(pBufferName,argv[1],255); strtolower(argv[1]); status = handleBatchPath(self,pCon,argc,argv); if(status >= 0){ diff --git a/macro.c b/macro.c index c79fd1bc..78cbe346 100644 --- a/macro.c +++ b/macro.c @@ -958,8 +958,8 @@ static int ProtectedExec(ClientData clientData, Tcl_Interp *interp, Arg2Text(argc-1,&argv[1],pCommand,1023); strtolower(argv[0]); if(strcmp(argv[0],"fulltransact") == 0){ - snprintf(pStart,1024,"TRANSACTIONSTART %s",pCommand); - SCWrite(pCon,pStart,eError); + snprintf(pStart,1023,"TRANSACTIONSTART %s",pCommand); + SCWrite(pCon,pStart,eError); } iRet = InterpExecute(pSics,pCon,pCommand); SCWrite(pCon,"TRANSACTIONFINISHED",eError); diff --git a/make_gen b/make_gen index 887afd3d..4fefe724 100644 --- a/make_gen +++ b/make_gen @@ -15,13 +15,13 @@ SOBJ = network.o ifile.o conman.o SCinter.o splitter.o passwd.o \ script.o o2t.o alias.o napi.o nxdata.o stringdict.o sdynar.o\ histmem.o histdriv.o histsim.o interface.o callback.o \ event.o emon.o evcontroller.o evdriver.o simev.o perfmon.o \ - danu.o nxdict.o varlog.o stptok.o nread.o \ + danu.o nxdict.o varlog.o stptok.o nread.o trigd.o cell.o\ scan.o fitcenter.o telnet.o token.o wwildcard.o hklmot.o\ tclev.o hkl.o integrate.o optimise.o dynstring.o nxutil.o \ mesure.o uubuffer.o commandlog.o udpquieck.o fourtable.o\ - rmtrail.o help.o nxupdate.o confvirtualmot.o \ + rmtrail.o help.o nxupdate.o confvirtualmot.o vector.o\ simchop.o choco.o chadapter.o trim.o scaldate.o \ - hklscan.o xytable.o exebuf.o exeman.o\ + hklscan.o xytable.o exebuf.o exeman.o ubfour.o ubcalc.o\ circular.o maximize.o sicscron.o scanvar.o \ d_sign.o d_mod.o tcldrivable.o stdscan.o diffscan.o\ synchronize.o definealias.o oscillate.o \ diff --git a/ofac.c b/ofac.c index 9375868f..ed91d189 100644 --- a/ofac.c +++ b/ofac.c @@ -109,6 +109,7 @@ #include "oscillate.h" #include "diffscan.h" #include "hklmot.h" +#include "ubcalc.h" /*----------------------- Server options creation -------------------------*/ static int IFServerOption(SConnection *pCon, SicsInterp *pSics, void *pData, int argc, char *argv[]) @@ -294,6 +295,8 @@ MakeDiffScan,NULL,NULL); AddCommand(pInter,"MakeHKLMot", HKLMotInstall,NULL,NULL); + AddCommand(pInter,"MakeUBCalc", + MakeUBCalc,NULL,NULL); /* @@ -356,6 +359,7 @@ RemoveCommand(pSics,"MakeOscillator"); RemoveCommand(pSics,"MakeDiffScan"); RemoveCommand(pSics,"MakeHKLMot"); + RemoveCommand(pSics,"MakeUBCalc"); /* remove site specific installation commands @@ -365,10 +369,7 @@ site->RemoveSiteCommands(pSics); } } - - /*--------------------------------------------------------------------------*/ - int InitObjectCommands(pServer pServ, char *file) { SConnection *pCon = NULL; diff --git a/scan.c b/scan.c index ae3148f2..07562243 100644 --- a/scan.c +++ b/scan.c @@ -42,37 +42,6 @@ #include "lld.h" #include "stdscan.h" -/*-----------------------------------------------------------------------*/ -static char *fixExtension(char *filename) -{ - if(strstr(filename,".hdf") != NULL) - { - changeExtension(filename,"dat"); - } - return filename; -} -/*------------------------------------------------------------------------*/ - char *ScanMakeFileName(SicsInterp *pSics, SConnection *pCon) - { - pSicsVariable pPath = NULL, pPref = NULL, pEnd = NULL; - char *pRes = NULL; - int iLen, iNum, iYear; - char pNumText[10]; - CommandList *pCom = NULL; - - /* - make a simulated filename if in simulation mode - */ - if(pServ->simMode) - return strdup("sim001001901.sim"); - - pRes = makeFilename(pSics,pCon); - if(pRes == NULL) - { - pRes = strdup("emergency.scn"); - } - return fixExtension(pRes); - } /*---------------------------------------------------------------------------*/ static void DeleteCountEntry(void *pData) { @@ -688,21 +657,7 @@ CountEntry CollectCounterData(pScanData self) return 0; } - /* allocate a new data file */ - pPtr = ScanMakeFileName(self->pSics,self->pCon); - if(!pPtr) - { - SCWrite(self->pCon, - "ERROR: cannot allocate new data filename, Scan aborted", - eError); - self->pCon = NULL; - self->pSics = NULL; - return 0; - } - sprintf(pBueffel,"Writing data file: %s ...",pPtr); - SCWrite(self->pCon,pBueffel,eWarning); - strcpy(self->pFile,pPtr); - free(pPtr); + iRet = self->WriteHeader(self); if(!iRet) { diff --git a/stdscan.c b/stdscan.c index 7c6354a4..80197749 100644 --- a/stdscan.c +++ b/stdscan.c @@ -30,6 +30,37 @@ #include "site.h" #include "lld.h" +/*-----------------------------------------------------------------------*/ +static char *fixExtension(char *filename) +{ + if(strstr(filename,".hdf") != NULL) + { + changeExtension(filename,"dat"); + } + return filename; +} +/*------------------------------------------------------------------------*/ + char *ScanMakeFileName(SicsInterp *pSics, SConnection *pCon) + { + pSicsVariable pPath = NULL, pPref = NULL, pEnd = NULL; + char *pRes = NULL; + int iLen, iNum, iYear; + char pNumText[10]; + CommandList *pCom = NULL; + + /* + make a simulated filename if in simulation mode + */ + if(pServ->simMode) + return strdup("sim001001901.sim"); + + pRes = makeFilename(pSics,pCon); + if(pRes == NULL) + { + pRes = strdup("emergency.scn"); + } + return fixExtension(pRes); + } /*---------------------------------------------------------------------------*/ int WriteHeader(pScanData self) { @@ -343,6 +374,7 @@ float fVal; char pBueffel[512]; char pMessage[1024]; + char *pPtr = NULL; assert(self); assert(self->iNP > 0); @@ -390,6 +422,22 @@ SetCounterMode((pCounter)self->pCounterData,self->iMode); SetCounterPreset((pCounter)self->pCounterData, self->fPreset); self->iCounts = 0; + + /* allocate a new data file */ + pPtr = ScanMakeFileName(self->pSics,self->pCon); + if(!pPtr) + { + SCWrite(self->pCon, + "ERROR: cannot allocate new data filename, Scan aborted", + eError); + self->pCon = NULL; + self->pSics = NULL; + return 0; + } + snprintf(pBueffel,511,"Writing data file: %s ...",pPtr); + SCWrite(self->pCon,pBueffel,eWarning); + strcpy(self->pFile,pPtr); + free(pPtr); return 1; } @@ -402,6 +450,7 @@ float fVal; char pBueffel[512]; char pMessage[1024]; + char *pPtr = NULL; assert(self); assert(self->iNP > 0); @@ -431,6 +480,22 @@ SetCounterMode((pCounter)self->pCounterData,self->iMode); SetCounterPreset((pCounter)self->pCounterData, self->fPreset); self->iCounts = 0; + + /* allocate a new data file */ + pPtr = ScanMakeFileName(self->pSics,self->pCon); + if(!pPtr) + { + SCWrite(self->pCon, + "ERROR: cannot allocate new data filename, Scan aborted", + eError); + self->pCon = NULL; + self->pSics = NULL; + return 0; + } + snprintf(pBueffel,511,"Writing data file: %s ...",pPtr); + SCWrite(self->pCon,pBueffel,eWarning); + strcpy(self->pFile,pPtr); + free(pPtr); return 1; } diff --git a/stdscan.h b/stdscan.h index 135ffbe7..f348dac6 100644 --- a/stdscan.h +++ b/stdscan.h @@ -11,6 +11,10 @@ #ifndef SICSSTDSCAN #define SICSSTDSCAN + /** + * make a scan file name + */ + char *ScanMakeFileName(SicsInterp *pSics, SConnection *pCon); /** * write the header of the scan file */ diff --git a/trigd.c b/trigd.c new file mode 100644 index 00000000..b3cb8953 --- /dev/null +++ b/trigd.c @@ -0,0 +1,92 @@ +/** + * This is a library of trigonmetric functions acting upon proper + * angles and not radians. Lifted from the Risoe tascom code + * + * March 2005 + */ +#include +#include "trigd.h" + +/* define constants */ +#ifndef PI +#define PI (3.1415926536) /* pi */ +#endif +#define DEGREE_RAD (PI/180.0) /* Radians per degree */ +/*******************************************************************************/ +double Sign(double d) +{ + if(d < .0){ + return -1; + } else { + return 1; + } +} +/******************************************************************************* +* Sinus of angle in degrees. +*******************************************************************************/ +extern double Sind (double x) +{ + return (sin (x*DEGREE_RAD)); +} +/******************************************************************************* +* Tangens of angle in degrees. +*******************************************************************************/ +extern double Tand(double x) +{ + return (tan(x*DEGREE_RAD)); +} +/******************************************************************************* +* Cosinus of angle in degrees. +*******************************************************************************/ +extern double Cosd (double x) +{ + return (cos (x*DEGREE_RAD)); +} +/******************************************************************************* +* Atan of angle in degrees. +*******************************************************************************/ +extern double Atand (double x) +{ + double data; + + data = (atan(x)/DEGREE_RAD); + return (data); +} +/******************************************************************************* +* Atan2 of angle in degrees. +*******************************************************************************/ +extern double Atand2 (double x) +{ + double data; + + data = (atan(x)/DEGREE_RAD); + return (data); +} +/******************************************************************************* +* Acos of angle in degrees. +*******************************************************************************/ +extern double Acosd (double x) +{ + double data; + + data = acos(x)/DEGREE_RAD; + return (data); +} +/******************************************************************************* +* Asin of angle in degrees. +*******************************************************************************/ +extern double Asind (double x) +{ + double data; + + data = x*x; + if (data == 1.0) + return (180.00 - Sign (x)*90.00); + else if (data > 1) + { + return (0); + } + else + return (Atand (x/sqrt (1 - data))); + +} diff --git a/trigd.h b/trigd.h new file mode 100644 index 00000000..bf3c6c47 --- /dev/null +++ b/trigd.h @@ -0,0 +1,17 @@ +/** + * This is a library of trigonmetric functions acting upon proper + * angles and not radians. Lifted from the Risoe tascom code + * + * March 2005 + */ +#ifndef SICSTRIGD +#define SICSTRIGD + double Sign(double d); + double Sind (double x); + double Tand(double x); + double Cosd (double x); + double Atand (double x); + double Atand2 (double x); + double Acosd (double x); + double Asind (double x); +#endif diff --git a/ubcalc.c b/ubcalc.c new file mode 100644 index 00000000..3d66027f --- /dev/null +++ b/ubcalc.c @@ -0,0 +1,389 @@ +/*---------------------------------------------------------------------- + UB calculation routines for four circle diffraction. + This is the interpreter interface to functionality implemented + in fourlib.c + + copyright: see file COPYRIGHT + + Mark Koennecke, March 2005 +-----------------------------------------------------------------------*/ +#include +#include +#include +#include "tcl.h" +#include "fortify.h" +#include "sics.h" +#include "ubfour.h" +#include "ubcalc.h" +#include "motor.h" +#include "hkl.h" +/*---------------------------------------------------------------------*/ +typedef struct { + pObjectDescriptor pDes; + lattice direct; + reflection r1, r2; + MATRIX UB; +} UBCALC, *pUBCALC; +/*----------------------------------------------------------------------*/ +static void killUBCALC(void *pData){ + pUBCALC self = (pUBCALC)pData; + if(self == NULL){ + return; + } + if(self->pDes != NULL){ + DeleteDescriptor(self->pDes); + } + if(self->UB != NULL){ + mat_free(self->UB); + } + free(self); +} +/*--------------------------------------------------------------------*/ +static int SaveUBCalc(void *data, char *name, FILE *fd){ + pUBCALC self = (pUBCALC)data; + if(self == NULL){ + return 0; + } + fprintf(fd,"%s cell %f %f %f %f %f %f\n", name, self->direct.a, + self->direct.b, self->direct.c, self->direct.alpha, + self->direct.beta, self->direct.gamma); + fprintf(fd,"%s ref1 %f %f %f %f %f %f %f\n",name, + self->r1.h, self->r1.k, self->r1.l, + self->r1.s2t, self->r1.om, self->r1.chi, self->r1.phi); + fprintf(fd,"%s ref2 %f %f %f %f %f %f %f\n",name, + self->r2.h, self->r2.k, self->r2.l, + self->r2.s2t, self->r2.om, self->r2.chi, self->r2.phi); + return 1; +} +/*---------------------------------------------------------------------*/ +static pUBCALC makeUBCALC(){ + pUBCALC pNew = NULL; + + pNew = (pUBCALC)malloc(sizeof(UBCALC)); + if(pNew == NULL){ + return NULL; + } + memset(pNew,0,sizeof(UBCALC)); + pNew->pDes = CreateDescriptor("UBcalc"); + pNew->UB = mat_creat(3,3,UNIT_MATRIX); + if(pNew->pDes == NULL || pNew->UB == NULL){ + return NULL; + } + pNew->pDes->SaveStatus = SaveUBCalc; + return pNew; +} +/*----------------------------------------------------------------------*/ +int MakeUBCalc(SConnection *pCon, SicsInterp *pSics, void *pData, + int argc, char *argv[]){ + pUBCALC pNew = NULL; + int status; + + pNew = makeUBCALC(); + if(pNew == NULL){ + SCWrite(pCon,"ERROR: out of memory creating UBCALC",eError); + return 0; + } + if(argc > 1){ + status = AddCommand(pSics,argv[1],UBCalcWrapper,killUBCALC,pNew); + } else { + status = AddCommand(pSics,"ubcalc",UBCalcWrapper,killUBCALC,pNew); + } + if(status != 1){ + SCWrite(pCon,"ERROR: failed to create duplicate UBCALC module",eError); + } + return status; +} +/*---------------------------------------------------------------------*/ +static int readCell(SConnection *pCon, pUBCALC self, int argc, char *argv[]){ + int status; + Tcl_Interp *pTcl = InterpGetTcl(pServ->pSics); + char pBueffel[256]; + + if(argc < 8){ + SCWrite(pCon,"ERROR: insufficient number of arguments to ubcalc cell", + eError); + return 0; + } + status = Tcl_GetDouble(pTcl,argv[2],&self->direct.a); + if(status != TCL_OK){ + snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[2]); + SCWrite(pCon,pBueffel,eError); + return 0; + } + status = Tcl_GetDouble(pTcl,argv[3],&self->direct.b); + if(status != TCL_OK){ + snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[3]); + SCWrite(pCon,pBueffel,eError); + return 0; + } + status = Tcl_GetDouble(pTcl,argv[4],&self->direct.c); + if(status != TCL_OK){ + snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[4]); + SCWrite(pCon,pBueffel,eError); + return 0; + } + status = Tcl_GetDouble(pTcl,argv[5],&self->direct.alpha); + if(status != TCL_OK){ + snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[5]); + SCWrite(pCon,pBueffel,eError); + return 0; + } + status = Tcl_GetDouble(pTcl,argv[6],&self->direct.beta); + if(status != TCL_OK){ + snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[6]); + SCWrite(pCon,pBueffel,eError); + return 0; + } + status = Tcl_GetDouble(pTcl,argv[7],&self->direct.gamma); + if(status != TCL_OK){ + snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[7]); + SCWrite(pCon,pBueffel,eError); + return 0; + } + SCparChange(pCon); + SCSendOK(pCon); + return 1; +} +/*--------------------------------------------------------------------*/ +static int readRefMotors(SConnection *pCon, SicsInterp *pSics, + reflection *r){ + pMotor pMot = NULL; + float val; + + pMot = FindMotor(pSics,"stt"); + if(pMot == NULL){ + SCWrite(pCon,"ERROR: cannot find stt motor",eError); + return 0; + } + MotorGetSoftPosition(pMot,pCon,&val); + r->s2t = val; + + pMot = FindMotor(pSics,"om"); + if(pMot == NULL){ + SCWrite(pCon,"ERROR: cannot find om motor",eError); + return 0; + } + MotorGetSoftPosition(pMot,pCon,&val); + r->om = val; + + pMot = FindMotor(pSics,"chi"); + if(pMot == NULL){ + SCWrite(pCon,"ERROR: cannot find chi motor",eError); + return 0; + } + MotorGetSoftPosition(pMot,pCon,&val); + r->chi = val; + + pMot = FindMotor(pSics,"phi"); + if(pMot == NULL){ + SCWrite(pCon,"ERROR: cannot find phi motor",eError); + return 0; + } + MotorGetSoftPosition(pMot,pCon,&val); + r->phi = val; + + SCparChange(pCon); + SCSendOK(pCon); + return 1; +} +/*---------------------------------------------------------------------*/ +static int readReflection(SConnection *pCon, SicsInterp *pSics, reflection *r, + int argc, char *argv[]){ + char pBueffel[255]; + int status; + Tcl_Interp *pTcl = InterpGetTcl(pSics); + + if(argc < 5){ + SCWrite(pCon,"ERROR: not enough arguments to ubcalc ref1,2",eError); + return 0; + } + status = Tcl_GetDouble(pTcl,argv[2],&r->h); + if(status != TCL_OK){ + snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[2]); + SCWrite(pCon,pBueffel,eError); + return 0; + } + status = Tcl_GetDouble(pTcl,argv[3],&r->k); + if(status != TCL_OK){ + snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[3]); + SCWrite(pCon,pBueffel,eError); + return 0; + } + status = Tcl_GetDouble(pTcl,argv[4],&r->l); + if(status != TCL_OK){ + snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[4]); + SCWrite(pCon,pBueffel,eError); + return 0; + } + + if(argc >= 9){ + status = Tcl_GetDouble(pTcl,argv[5],&r->s2t); + if(status != TCL_OK){ + snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[5]); + SCWrite(pCon,pBueffel,eError); + return 0; + } + status = Tcl_GetDouble(pTcl,argv[6],&r->om); + if(status != TCL_OK){ + snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[6]); + SCWrite(pCon,pBueffel,eError); + return 0; + } + status = Tcl_GetDouble(pTcl,argv[7],&r->chi); + if(status != TCL_OK){ + snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[7]); + SCWrite(pCon,pBueffel,eError); + return 0; + } + status = Tcl_GetDouble(pTcl,argv[8],&r->phi); + if(status != TCL_OK){ + snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[8]); + SCWrite(pCon,pBueffel,eError); + return 0; + } + SCparChange(pCon); + SCSendOK(pCon); + return 1; + } else { + return readRefMotors(pCon,pSics,r); + } +} +/*---------------------------------------------------------------------*/ +static void listUB(SConnection *pCon, MATRIX UB){ + Tcl_DString list; + char pBueffel[255]; + int i; + + + Tcl_DStringInit(&list); + if(UB == NULL){ + Tcl_DStringAppend(&list,"NO UB",-1); + } else { + Tcl_DStringAppend(&list,"UB = ", -1); + for(i = 0; i < 3; i++){ + snprintf(pBueffel,255,"%f %f %f\n", UB[i][0], + UB[i][1],UB[i][2]); + Tcl_DStringAppend(&list,pBueffel,-1); + } + } + SCWrite(pCon,Tcl_DStringValue(&list),eValue); + Tcl_DStringFree(&list); +} +/*---------------------------------------------------------------------*/ +static int calcUB(pUBCALC self, SConnection *pCon){ + MATRIX newUB = NULL; + int err = 1; + + newUB = calcUBFromCellAndReflections(self->direct, self->r1, self->r2, &err); + if(newUB == NULL){ + switch(err){ + case UBNOMEMORY: + SCWrite(pCon,"ERROR: out of memory while calculating UB",eError); + return 0; + break; + case REC_NO_VOLUME: + SCWrite(pCon,"ERROR: bad cell constants",eError); + return 0; + break; + default: + SCWrite(pCon,"ERROR: unknown error on UB matrix calculation",eError); + return 0; + } + } else { + if(self->UB != NULL){ + mat_free(self->UB); + } + self->UB = newUB; + SCSendOK(pCon); + return 1; + } +} +/*---------------------------------------------------------------------*/ +static void listCell(SConnection *pCon, char *name, lattice direct){ + char pBueffel[255]; + + snprintf(pBueffel,255,"%s.cell = %f %f %f %f %f %f", + name,direct.a, direct.b,direct.c, + direct.alpha,direct.beta,direct.gamma); + SCWrite(pCon,pBueffel,eValue); + +} +/*---------------------------------------------------------------------*/ +static void listReflection(SConnection *pCon, char *name, + char *refName, reflection r){ + char pBueffel[255]; + + snprintf(pBueffel,255,"%s.%s = %f %f %f %f %f %f %f", + name,refName, r.h, r.k, r.l, r.s2t, r.om, r.chi, r.phi); + SCWrite(pCon,pBueffel,eValue); +} +/*---------------------------------------------------------------------*/ +static int sendUBToHKL(SConnection *pCon, SicsInterp *pSics, + char *name, MATRIX UB){ + pHKL hkl = NULL; + char pBueffel[256]; + float ub[9]; + int i; + + hkl = FindCommandData(pSics,name,"4-Circle-Calculus"); + if(hkl == NULL){ + snprintf(pBueffel,255,"ERROR: %s not found or of wrong type",name); + SCWrite(pCon,pBueffel,eError); + return 0; + } + for(i = 0; i < 3; i++){ + ub[i] = UB[0][i]; + ub[i+3] = UB[1][i]; + ub[i+6] = UB[2][i]; + } + SetUB(hkl,ub); + SCSendOK(pCon); + return 1; +} +/*----------------------------------------------------------------------*/ +int UBCalcWrapper(SConnection *pCon, SicsInterp *pSics, void *pData, + int argc, char *argv[]){ + pUBCALC self = (pUBCALC)pData; + + assert(self); + if(argc < 2){ + SCWrite(pCon,"Insuffcient number of arguments to ubcalc",eError); + return 0; + } + strtolower(argv[1]); + + if(strcmp(argv[1],"cell") == 0){ + return readCell(pCon, self, argc, argv); + } else if(strcmp(argv[1],"ref1") == 0){ + return readReflection(pCon,pSics, &self->r1,argc,argv); + } else if(strcmp(argv[1],"ref2") ==0){ + return readReflection(pCon,pSics, &self->r2,argc,argv); + } else if(strcmp(argv[1],"listub") == 0){ + listUB(pCon,self->UB); + return 1; + } else if(strcmp(argv[1],"calcub") == 0){ + return calcUB(self,pCon); + } else if(strcmp(argv[1],"listcell") == 0){ + listCell(pCon,argv[0],self->direct); + return 1; + } else if(strcmp(argv[1],"listref1") == 0){ + listReflection(pCon,argv[0],"ref1",self->r1); + return 1; + } else if(strcmp(argv[1],"listref2") == 0){ + listReflection(pCon,argv[0],"ref2",self->r2); + return 1; + } else if(strcmp(argv[1],"list") == 0){ + listCell(pCon,argv[0],self->direct); + listReflection(pCon,argv[0],"ref1",self->r1); + listReflection(pCon,argv[0],"ref2",self->r2); + listUB(pCon,self->UB); + return 1; + } else if(strcmp(argv[1],"sendto") == 0){ + if(argc < 3){ + SCWrite(pCon,"ERROR: Nothing to send UB too",eError); + return 0; + } + return sendUBToHKL(pCon,pSics,argv[2],self->UB); + } + return 1; +} diff --git a/ubcalc.h b/ubcalc.h new file mode 100644 index 00000000..11f6594b --- /dev/null +++ b/ubcalc.h @@ -0,0 +1,19 @@ + +/*---------------------------------------------------------------------- + UB caclualtion routines for four circle diffraction. + This is the interpreter interface to functionality implemented + in fourlib.c + + copyright: see file COPYRIGHT + + Mark Koennecke, March 2005 +-----------------------------------------------------------------------*/ +#ifndef SICSUBCALC +#define SICSUBCALC + +int MakeUBCalc(SConnection *pCon,SicsInterp *pSics, void *pData, + int argc, char *argv[]); +int UBCalcWrapper(SConnection *pCon, SicsInterp *pSics, void *pData, + int argc, char *argv[]); + +#endif diff --git a/ubcalc.w b/ubcalc.w new file mode 100644 index 00000000..c21a082a --- /dev/null +++ b/ubcalc.w @@ -0,0 +1,28 @@ +\subsection{UB Matrix Calculation for Four Circle Diffractometers} +This module helps in the calculation of UB matrices for four + circle diffraction. This is only an interpreter interface, the + actual functionality is in the ubfour.h, .c files. These are documented + in ubfour.h in order to be able to give this away separatly as a + library. + +@o ubcalc.h @{ +/*---------------------------------------------------------------------- + UB caclualtion routines for four circle diffraction. + This is the interpreter interface to functionality implemented + in fourlib.c + + copyright: see file COPYRIGHT + + Mark Koennecke, March 2005 +-----------------------------------------------------------------------*/ +#ifndef SICSUBCALC +#define SICSUBCALC + +int MakeUBCalc(SConnection *pCon,SicsInterp *pSics, void *pData, + int argc, char *argv[]); +int UBCalcWrapper(SConnection *pCon, SicsInterp *pSics, void *pData, + int argc, char *argv[]); + +#endif +@} + \ No newline at end of file diff --git a/ubfour.c b/ubfour.c new file mode 100644 index 00000000..ac5ec6e1 --- /dev/null +++ b/ubfour.c @@ -0,0 +1,136 @@ +/** + * This is a library for calculating UB matrices for four circle diffraction. + * The algorithm and setting definitions is from: + * + * Busing & Levy, Acta Cryst. (1967), 22, 457ff + * + * Implemented: + * - UB from cell cell constants and two reflections. + * + * Mark Koennecke, March 2005 + */ +#include "ubfour.h" +#include +#include "vector.h" +#include "trigd.h" +/*--------------------------------------------------------------------------------------*/ +static MATRIX calcUVectorFromAngles(reflection r){ + MATRIX u; + + u = makeVector(); + if(u == NULL){ + return NULL; + } + vectorSet(u,0,Cosd(r.om)*Cosd(r.chi)*Cosd(r.phi) - Sind(r.om)*Sind(r.phi)); + vectorSet(u,1,Cosd(r.om)*Cosd(r.chi)*Sind(r.phi) + Sind(r.om)*Cosd(r.phi)); + vectorSet(u,2,Cosd(r.om)*Sind(r.chi)); + + return u; +} +/*--------------------------------------------------------------------------------------*/ +static MATRIX reflectionToHC(reflection r, MATRIX B){ + MATRIX h = NULL, hc = NULL; + + h = makeVector(); + if(h == NULL){ + return NULL; + } + vectorSet(h,0,r.h); + vectorSet(h,1,r.k); + vectorSet(h,2,r.l); + + hc = mat_mul(B,h); + killVector(h); + return hc; +} +/*---------------------------------------------------------------------------------------*/ +MATRIX calcUBFromCellAndReflections(lattice direct, reflection r1, + reflection r2, int *errCode){ + MATRIX B, HT, UT, U, UB, HTT ; + MATRIX u1, u2, h1, h2; + double ud[3]; + int status; + + *errCode = 1; + + /* + calculate the B matrix and the HT matrix + */ + B = mat_creat(3,3,ZERO_MATRIX); + status = calculateBMatrix(direct,B); + if(status < 0){ + *errCode = status; + return NULL; + } + h1 = reflectionToHC(r1,B); + h2 = reflectionToHC(r2,B); + if(h1 == NULL || h2 == NULL){ + *errCode = UBNOMEMORY; + return NULL; + } + HT = matFromTwoVectors(h1,h2); + if(HT == NULL){ + *errCode = UBNOMEMORY; + return NULL; + } + + /* + calculate U vectors and UT matrix + */ + u1 = calcUVectorFromAngles(r1); + u2 = calcUVectorFromAngles(r2); + if(u1 == NULL || u2 == NULL){ + *errCode = UBNOMEMORY; + return NULL; + } + UT = matFromTwoVectors(u1,u2); + if(UT == NULL){ + *errCode = UBNOMEMORY; + return NULL; + } + + /* + debugging output + printf("B-matrix\n"); + mat_dump(B); + printf("HT-matrix\n"); + mat_dump(HT); + printf("UT-matrix\n"); + mat_dump(UT); + */ + + /* + UT = U * HT + */ + HTT = mat_tran(HT); + if(HTT == NULL){ + *errCode = UBNOMEMORY; + return NULL; + } + U = mat_mul(UT,HTT); + if(U == NULL){ + *errCode = UBNOMEMORY; + return NULL; + } + UB = mat_mul(U,B); + if(UB == NULL){ + *errCode = UBNOMEMORY; + } + + /* + clean up + */ + killVector(h1); + killVector(h2); + mat_free(HT); + mat_free(HTT); + + killVector(u1); + killVector(u2); + mat_free(UT); + + mat_free(U); + mat_free(B); + + return UB; +} diff --git a/ubfour.h b/ubfour.h new file mode 100644 index 00000000..84134b68 --- /dev/null +++ b/ubfour.h @@ -0,0 +1,40 @@ +/** + * This is a library for calculating UB matrices for four circle diffraction. + * The algorithm and setting definitions is from: + * Busing & Levy, Acta Cryst. (1967), 22, 457ff + * + * This initial version only supports the calculation of the UB matrix from + * cell constants and two reflections. + * + * Mark Koennecke, march 2005 + */ +#ifndef SICSUBFOUR +#define SICSUBFOUR + +#include +#include "matrix/matrix.h" +#include "cell.h" + +/** + * error codes: also see cell.h + */ +#define UBNOMEMORY -200 + +/** + * a reflection data structure holding the indices h,k,l and + * the magic four circle angles two_theta, om, chi and phi. + */ +typedef struct{ + double h,k,l; + double s2t,om,chi,phi; +}reflection; +/** + * calculate a UB matrix from cell constants and two reflections + * @param direct The direct cell constants + * @param r1 The first reflection. + * @param r2 The second reflection. + * @param errCode an error indicator if things go wrong. + * @return The resulting UB matrix or NULL on errors + */ +MATRIX calcUBFromCellAndReflections(lattice direct, reflection r1, reflection r2, int *errCode); +#endif diff --git a/vector.c b/vector.c new file mode 100644 index 00000000..e9598f62 --- /dev/null +++ b/vector.c @@ -0,0 +1,145 @@ +/** + * This is a little collection of vector functions for use with the UB + * calculation routines. Therefore vectors are mapped onto the matrix + * package. Strictly speaking all routines only work properly on 3D + * vectors. + * + * copyright: see file COPYRIGHT + * + * Mark Koennecke, March 2005 + */ +#include +#include +#include +#include "vector.h" +/*----------------------------------------------------------------------*/ +MATRIX makeVector(){ + return mat_creat(3,1,ZERO_MATRIX); +} +/*---------------------------------------------------------------------*/ +MATRIX makeVectorInit(double val[3]){ + int i; + MATRIX result = NULL; + + result = makeVector(); + if(result == NULL){ + return result; + } + for(i = 0; i < 3; i++){ + result[i][0] = val[i]; + } + return result; +} +/*---------------------------------------------------------------------*/ +void vectorToArray(MATRIX v, double val[3]){ + int i; + + assert(MatCol(v) == 3); + for(i = 0; i < 3; i++){ + val[i] = v[i][0]; + } +} +/*----------------------------------------------------------------------*/ +void killVector(MATRIX v){ + mat_free(v); +} +/*----------------------------------------------------------------------*/ +void vectorSet(MATRIX v, int idx, double value){ + assert(idx >= 0 && idx < 3); + + v[idx][0] = value; +} +/*----------------------------------------------------------------------*/ +double vectorGet(MATRIX v, int idx){ + assert(idx >= 0 && idx < 3); + + return v[idx][0]; +} +/*----------------------------------------------------------------------*/ +double vectorLength(MATRIX v){ + assert(MatRow(v) == 3); + + return sqrt(v[0][0]*v[0][0] + v[1][0]*v[1][0] + v[2][0]*v[2][0]); +} +/*---------------------------------------------------------------------*/ +void normalizeVector(MATRIX v){ + int i; + double norm; + + norm = vectorLength(v); + if(norm > .001) { + for(i = 0; i < 3; i++){ + v[i][0] /= norm; + } + } else { + for(i = 0; i < 3; i++){ + v[i][0] = .0; + } + } +} +/*----------------------------------------------------------------------*/ +double vectorDotProduct(MATRIX v1, MATRIX v2){ + double sum; + int i; + + assert(MatRow(v1) == MatRow(v2)); + + sum = .0; + for(i = 0; i < MatRow(v1); i++){ + sum += v1[i][0]*v2[i][0]; + } + return sum; +} +/*----------------------------------------------------------------------*/ +MATRIX vectorCrossProduct(MATRIX v1, MATRIX v2){ + MATRIX result; + + assert(MatRow(v1) == 3); + assert(MatRow(v2) == 3); + + result = makeVector(); + if(result == NULL){ + return NULL; + } + vectorSet(result,0,vectorGet(v1,1)*vectorGet(v2,2) - vectorGet(v1,2)*vectorGet(v2,1)); + vectorSet(result,1,vectorGet(v1,2)*vectorGet(v2,0) - vectorGet(v1,0)*vectorGet(v2,2)); + vectorSet(result,2,vectorGet(v1,0)*vectorGet(v2,1) - vectorGet(v1,1)*vectorGet(v2,0)); + return result; +} +/*-------------------------------------------------------------------------*/ +MATRIX matFromTwoVectors(MATRIX v1, MATRIX v2){ + MATRIX a1, a2, a3, result; + int i; + + a1 = mat_copy(v1); + if(a1 == NULL){ + return NULL; + } + normalizeVector(a1); + + a3 = vectorCrossProduct(v1,v2); + if(a3 == NULL){ + return NULL; + } + normalizeVector(a3); + + a2 = vectorCrossProduct(a1,a3); + if(a2 == NULL){ + return NULL; + } + + result = mat_creat(3,3,ZERO_MATRIX); + if(result == NULL){ + return NULL; + } + + for(i = 0; i < 3; i++){ + result[i][0] = vectorGet(a1,i); + result[i][1] = vectorGet(a2,i); + result[i][2] = vectorGet(a3,i); + } + killVector(a1); + killVector(a2); + killVector(a3); + return result; +} diff --git a/vector.h b/vector.h new file mode 100644 index 00000000..79d3c888 --- /dev/null +++ b/vector.h @@ -0,0 +1,93 @@ +/** + * This is a little collection of vector functions for use with the UB + * calculation routines. Therefore vectors are mapped onto the matrix + * package. Strictly speaking all routine only work properly on 3D + * vectors. + * + * copyright: see file COPYRIGHT + * + * Mark Koennecke, March 2005 + */ +#ifndef SICSVECTOR +#define SICSVECTOR +#include +#include "matrix/matrix.h" + +/** + * make a 3D vector + * @return a 3D vector expressed as a MATRIX. The caller is resonsible + * for removing any memory associated with this vector. + */ +MATRIX makeVector(); +/** + * make a vector and initialize with values given in val + * @param val Array of values for the new vector + * @return a new vector with the values in val + */ +MATRIX makeVectorInit(double val[3]); +/** + * copy the value in v to the array val + * @param v The vector to copy + * @param val The array to which to copy the values + */ +void vectorToArray(MATRIX v, double val[3]); +/** + * delete the vector v + * @param v The vector to delete. + */ +void killVector(MATRIX v); +/** + * set a component of a vector + * @param v The vector to which to apply the new value + * @param idx The index of the vector component to set + * @param val The new value for the component. + */ +void vectorSet(MATRIX v, int idx, double value); +/** + * get the value of a vector component + * @param v The vector to query + * @param idx The index of the component + * @return The value of the vectors component. + */ +double vectorGet(MATRIX v, int idx); +/** + * calculate the length of the vector v + * @param v The vector to calculate the length for + * @return The length of the vector + */ +double vectorLength(MATRIX v); +/** + * normalize the vector by making it unit length + * @param v The vector to normalize. + */ +void normalizeVector(MATRIX v); +/** + * calculate the dot or scalar product of two vectors + * @param v1 First vector + * @param v2 Second vector + * @return The dot product v1*v2 + */ +double vectorDotProduct(MATRIX v1, MATRIX v2); +/** + * calculate the cross product of vectors v1 and v2. This is only + * correct when both vectors are expressed in terms of a + * cartesian coordinate system. + * @param v1 The first vector + * @param v2 The second vector + * @return The cross product of the vectors v1, v2 + */ +MATRIX vectorCrossProduct(MATRIX v1, MATRIX v2); +/** + * this is a special function used in UB matrix calculations. + * It first calculates a coordinate system from the two vectors + * where: + * a1 = v1/|v1| + * a2 = v1*v2/|v1*v2| + * a3 = a1*a2 + * and then constructs a matrix with the ai as columns. + * @param v1 The first vector + * @param v2 The second vector + * @return A matrix as descibed above. + */ +MATRIX matFromTwoVectors(MATRIX v1, MATRIX v2); +#endif