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