/*-------------------------------------------------------------------------- This is one implementation file for the TASMAD simulation module for SICS. The requirement is to make SICS look as much as TASMAD as possible. This includes: - TASMAD is variable driven - Sometimes variables are accessed in storage order. - A complicated calculation has to be done for getting the instruments settings right. The appropriate F77 routine from TASMAD will be reused. - The scan logic is different. - Output of ILL-formatted data files is required. This file implements a couple of utility functions mainly for parsing and second for interfacing to the F77 routines for calculating the TAS settings.. Mark Koennecke, November 2000 Polarisation support added. Mark Koennecke, April 2002 ---------------------------------------------------------------------------*/ #include #include #include #include #include #include #include #include #include "tas.h" #include "tasu.h" /*----------------------------------------------------------------------- isTASMotor finds out if the given string is a TAS motor and returns its index in the motor list if this the case. Else -1 is returned. -------------------------------------------------------------------------*/ int isTASMotor(char *val) { int iPtr; iPtr = 0; while (tasMotorOrder[iPtr] != NULL) { if (strcmp(val, tasMotorOrder[iPtr]) == 0) { return iPtr; } else { iPtr++; } } return -1; } /*------------------------------------------------------------------------ return number if val is a current variable, else -1 never used - and wrong (the number 18) M.Z. june 07 -----------------------------------------------------------------------*/ int isTASCurrent(char *val) { int test; test = isTASMotor(val); if (test < 0) { return test; } if (test > 18) { return test; } else { return -1; } } /*-------------------------------------------------------------------------- isTASVar finds out if the given string is a TAS variable and returns its index if this is the case. Else -1 will be returned. --------------------------------------------------------------------------*/ int isTASVar(char *val) { int iPtr; iPtr = 0; while (tasVariableOrder[iPtr] != NULL) { if (strcmp(val, tasVariableOrder[iPtr]) == 0) { return iPtr; } else { iPtr++; } } return -1; } /*------------------------------------------------------------------------ isTASEnergy finds out if a variable is a TAS energy variable. This would require a TAS calculation to take place. If the test is positive then the index of the Variable is returned, else -1 is returned. -------------------------------------------------------------------------*/ int isTASEnergy(char *val) { int iPtr; for (iPtr = EMIN; iPtr <= EMAX; iPtr++) { if (strcmp(val, tasVariableOrder[iPtr]) == 0) { return iPtr - EMIN; } } return -1; } /*----------------------------------------------------------------------- prepare2Parse removes all unnecessary delimiters from the parse string. These are: ' = ----------------------------------------------------------------------*/ void prepare2Parse(char *line) { int i, iLang; iLang = strlen(line); for (i = 0; i < iLang; i++) { if (line[i] == ',' || line[i] == '=') { line[i] = ' '; } } } /*---------------------------------------------------------------------*/ int tasNumeric(char *pText) { int i, ii, iGood; static char pNum[13] = { "1234567890.+-" }; for (i = 0; i < strlen(pText); i++) { for (ii = 0; ii < 13; ii++) { iGood = 0; if (pText[i] == pNum[ii]) { iGood = 1; break; } } if (!iGood) { return 0; } } 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: return 0; break; } } /*--------------------------------------------------------------------- getSRO reads the SRO motor. This motor is for the rotation of the magnet. The value has to be subtracted from the helm value during calculations -----------------------------------------------------------------------*/ int getSRO(SConnection * pCon, float *sroVal) { pMotor pMot = NULL; int status; pMot = FindMotor(pServ->pSics, "sro"); if (!pMot) { SCWrite(pCon, "ERROR: internal: motor sro not found!", eError); return 0; } status = MotorGetSoftPosition(pMot, pCon, sroVal); if (status != 1) { SCWrite(pCon, "ERROR: failed to read motor sro", eError); return 0; } return 1; } /*---------------------------------------------------------------------- TASCalc does the triple axis spectrometer calculations. This function is invoked whenever energy or Q needs to be driven. The parameters: - self a pointer to a TASdata structure holding information about all necessary parameters. - pCon a SConnection object needed for reporting errors. - tasMask a mask with 0,1 indicating which variable is driven. IT IS ASSUMED THAT THE APPROPRIATE VARIABLE HAS ALREADY BEEN SET IN THE INTERPRETER BEFORE TASCalc IS CALLED. This is because for convenience this functions collects the values from the interpreter. - motorTargets will hold the new targets for the relevant motors (A1-A6 plus some magnet related) when the call finished happily. - motorMask will hold 1's for each motor to drive, 0 for each motor not affected. This version uses the TASMAD calculation routines. These have been translated from F77 to C with f2c (f2c -A crysconv.f , tacov.f). The C is not maintainable, always edit the F77 sources in case of need and retranslate. This scheme requires the definitions going before TASCalc. The advantage of using f2c is that we get less portability problems with different calling conventions for C anf F77 on different platforms. -----------------------------------------------------------------------*/ #include "f2c.h" extern int setrlp_(real * sam, integer * ier); extern 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, real * da, real * helm, real * f1h, real * f1v, real * f2h, real * f2v, real * f, integer * ifx, integer * iss, integer * ism, integer * isa, real * t_a__, real * t_rm__, real * t_alm__, real * t_ra__, real * qm, logical * ldra, logical * ldr_rm__, logical * ldr_alm__, logical * ldr_ra__, real * p_ih__, real * c_ih__, doublereal * a4, integer * ier); extern int inicurve_(integer * midx, real * mrx1, real * mrx2, integer * aidx, real * arx1, real * arx2, real * mmin, real * mmax, real * amin, real * amax); /*-------------------------------------------------------------------*/ int TASCalc(pTASdata self, SConnection * pCon, unsigned char tasMask[MAXEVAR], float motorTargets[20], unsigned char motorMask[20]) { /* for setrlp */ real sam[13]; integer ier; /* for t_conv */ real ei, aki, ef, akf, qhkl[3], en, hx, hy, hz, dm, da, helm, f1h, f1v, f2h, f2v, f, angles[6], tRM, tRA, tALM, qm, currents[8], helmconv[4]; doublereal a4; integer if1, if2, ifx, iss, ism, isa; logical ldk[8], ldh, ldf, lpa, ldra[6], l_RM, l_RA, l_ALM; /* for inicurve */ real mrx1, mrx2, arx1, arx2, mmin, mmax, amin, amax; integer im, ia; pMotor pMot; float fVal, fSRO; /* else */ int i; assert(self); assert(pCon); /* The first thing we need to do is to set up the orientation matrix from the sample parameters. */ sam[0] = (real) self->tasPar[AS]->fVal; sam[1] = (real) self->tasPar[BS]->fVal; sam[2] = (real) self->tasPar[CS]->fVal; sam[3] = (real) self->tasPar[AA]->fVal; sam[4] = (real) self->tasPar[BB]->fVal; sam[5] = (real) self->tasPar[CC]->fVal; sam[6] = (real) self->tasPar[AX]->fVal; sam[7] = (real) self->tasPar[AY]->fVal; sam[8] = (real) self->tasPar[AZ]->fVal; sam[9] = (real) self->tasPar[BX]->fVal; sam[10] = (real) self->tasPar[BY]->fVal; sam[11] = (real) self->tasPar[BZ]->fVal; setrlp_(sam, &ier); switch (ier) { case 1: SCWrite(pCon, "ERROR: error on lattice parameters", eError); return 0; break; case 2: SCWrite(pCon, "ERROR: error on lattice angles", eError); return 0; break; case 3: SCWrite(pCon, "ERROR: error on vectors A1, A2", eError); return 0; break; default: break; } /* now we need to set the parameters for mono/analyzer curvature calculations. */ mrx1 = (real) self->tasPar[MRX1]->fVal; mrx2 = (real) self->tasPar[MRX2]->fVal; arx1 = (real) self->tasPar[ARX1]->fVal; arx2 = (real) self->tasPar[ARX2]->fVal; im = 7; ia = 9; pMot = FindMotor(pServ->pSics, "MCV"); if (!pMot) { SCWrite(pCon, "ERROR: cannot find monochromator curvature motor", eError); return 0; } MotorGetPar(pMot, "softupperlim", &fVal); mmax = (real) fVal; MotorGetPar(pMot, "softlowerlim", &fVal); mmin = (real) fVal; pMot = FindMotor(pServ->pSics, "ACH"); if (!pMot) { SCWrite(pCon, "ERROR: cannot find analyzer curvature motor", eError); return 0; } MotorGetPar(pMot, "softupperlim", &fVal); amax = (real) fVal; MotorGetPar(pMot, "softlowerlim", &fVal); amin = (real) fVal; inicurve_(&im, &mrx1, &mrx2, &ia, &arx1, &arx2, &mmin, &mmax, &amin, &amax); /* This done, we painstackingly initialize the tons of parameters required by the t_conv */ 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[THX]->fVal; hy = (real) self->tasPar[THY]->fVal; hz = (real) self->tasPar[THZ]->fVal; f = (real) 2.072; /* energy unit meV */ if1 = (integer) self->tasPar[F1]->iVal; if2 = (integer) self->tasPar[F2]->iVal; for (i = 0; i < 8; i++) { ldk[i] = (logical) tasMask[i]; } if (if1 > 0 || if2 > 0) { ldf = 1; } else { ldf = 0; } if (tasMask[9] || tasMask[10] || tasMask[11]) { ldh = 1; } else { ldh = 0; } 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; helm = (real) self->tasPar[HELM]->fVal; f1h = (real) self->tasPar[IF1H]->fVal; f1v = (real) self->tasPar[IF1V]->fVal; f2h = (real) self->tasPar[IF2H]->fVal; f2v = (real) self->tasPar[IF2V]->fVal; ifx = (integer) self->tasPar[FX]->iVal; iss = (integer) self->tasPar[SS]->iVal; ism = (integer) self->tasPar[SM]->iVal; isa = (integer) self->tasPar[SA]->iVal; /* initialize the magnet currents. */ for (i = 0; i < 4; i++) { helmconv[i] = .0; currents[i] = readDrivable(tasMotorOrder[CURMOT + i], pCon); currents[i + 4] = readDrivable(tasMotorOrder[CURMOT + i + 4], pCon); } readConversionFactors(self, helmconv); /* warning if energy fixed */ if (ifx == 0) { SCWrite(pCon, "WARNING: energy fixed", eWarning); } /* read a4 in order to make magnet calculations work */ a4 = .0; a4 = readDrivable(tasMotorOrder[A4MOT], pCon); /* initalise the motorMasks to 0. */ for (i = 0; i < 20; i++) { motorMask[i] = 0; } for (i = 0; i < 6; i++) { ldra[i] = 0; } l_RA = l_RM = l_ALM = ier = 0; /* now we can call */ t_conv__(&ei, &aki, &ef, &akf, qhkl, &en, &hx, &hy, &hz, &if1, &if2, ldk, &ldh, &ldf, &lpa, &dm, &da, &helm, &f1h, &f1v, &f2h, &f2v, &f, &ifx, &iss, &ism, &isa, angles, &tRM, &tALM, &tRA, &qm, ldra, &l_RM, &l_ALM, &l_RA, currents, helmconv, &a4, &ier); if (ier != 0) { printError(ier, pCon); return 0; } /* now we have to put those things back which interest us, including updating the variables which may have accidentally been set. */ for (i = 0; i < 6; i++) { motorTargets[i] = angles[i]; motorMask[i] = ldra[i]; } /* these additional motors are the curvature motors */ motorMask[6] = l_RM; motorTargets[6] = tRM; motorMask[8] = l_RA; motorTargets[8] = tRA; for (i = 0; i < 8; i++) { motorTargets[9 + i] = currents[i]; } if (ldh) { /* currents must be driven */ for (i = 0; i < 3; i++) { motorMask[13 + i] = 1; } } if (ldf) { /* currents must be driven */ for (i = 0; i < 4; i++) { motorMask[9 + i] = 1; } } 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; self->tasPar[TQM]->fVal = (float) qm; return 1; } /*------------------------------------------------------------------------- TASStart starts the motors calculated by TASCalc ---------------------------------------------------------------------------*/ int TASStart(pTASdata self, SConnection * pCon, SicsInterp * pSics, float motorTargets[20], unsigned char motorMask[20]) { int status, i; char pBueffel[80]; /* do the A1-A6 motors first */ for (i = 0; i < 6; i++) { if (motorMask[i] != 0) { status = StartMotor(pServ->pExecutor, pSics, pCon, tasMotorOrder[i], RUNDRIVE, motorTargets[i]); } } /* start the monochromator and analyzer curvature motors */ if (motorMask[6]) { status = StartMotor(pServ->pExecutor, pSics, pCon, "MCV", RUNDRIVE, motorTargets[6]); if (status == 0) { /* the error will have been reported but must be ignored */ SCSetInterrupt(pCon, eContinue); return 0; } } if (motorMask[8]) { status = StartMotor(pServ->pExecutor, pSics, pCon, "ACH", RUNDRIVE, motorTargets[8]); if (status == 0) { /* the error will have been reported but must be ignored */ SCSetInterrupt(pCon, eContinue); return 0; } } /* start magnet currents. */ for (i = 0; i < 8; i++) { if (motorMask[9 + i] == 1) { startCurrent(tasMotorOrder[CURMOT + i], pCon, motorTargets[9 + i]); } } return 1; } /*------------------ from t_update.f converted to t_update.c -----------*/ extern int t_update__(real * p_a__, real * p_ih__, real * c_ih__, logical * lpa, real * dm, real * da, integer * isa, real * helm, real * f1h, real * f1v, real * f2h, real * f2v, real * f, real * ei, real * aki, real * ef, real * akf, real * qhkl, real * en, real * hx, real * hy, real * hz, integer * if1, integer * if2, real * qm, integer * ier); /*--------------------------------------------------------------------------- TASUpdate calculates the Q/E variables back from the motors A1-A6 and then proceeds to update the energy variables. It also checks if the analyzer monochromator motors match up i.e tth = 2*om. ----------------------------------------------------------------------------*/ int TASUpdate(pTASdata self, SConnection * pCon) { real sam[16], amot[6], helmCurrent[8], convH[4], dm, da, helm, f1h, f1v; real f, ei, aki, ef, akf, qhkl[3], en, hx, hy, hz, qm, f2v, f2h; logical lpa; integer isa, if1, if2, ier, ifx; pMotor pMot; char pBueffel[132]; float fMot, diff, fSRO; int i, status; assert(self); assert(pCon); /* The first thing we need to do is to set up the orientation matrix from the sample parameters. */ sam[0] = (real) self->tasPar[AS]->fVal; sam[1] = (real) self->tasPar[BS]->fVal; sam[2] = (real) self->tasPar[CS]->fVal; sam[3] = (real) self->tasPar[AA]->fVal; sam[4] = (real) self->tasPar[BB]->fVal; sam[5] = (real) self->tasPar[CC]->fVal; sam[6] = (real) self->tasPar[AX]->fVal; sam[7] = (real) self->tasPar[AY]->fVal; sam[8] = (real) self->tasPar[AZ]->fVal; sam[9] = (real) self->tasPar[BX]->fVal; sam[10] = (real) self->tasPar[BY]->fVal; sam[11] = (real) self->tasPar[BZ]->fVal; setrlp_(sam, &ier); switch (ier) { case 1: SCWrite(pCon, "ERROR: error on lattice parameters", eError); return 0; break; case 2: SCWrite(pCon, "ERROR: error on lattice angles", eError); return 0; break; case 3: SCWrite(pCon, "ERROR: error on vectors A1, A2", eError); return 0; break; default: break; } /* correct helmholtz angle for SRO rotation */ if (getSRO(pCon, &fSRO)) { fSRO = self->oldSRO - fSRO; fSRO = self->tasPar[HELM]->fVal - fSRO; if (fSRO > 360) { fSRO -= 360.; } if (fSRO < -360) { fSRO += 360.; } getSRO(pCon, &self->oldSRO); self->tasPar[HELM]->fVal = fSRO; } /* get the motor angles A1-A6 */ for (i = 0; i < 6; i++) { pMot = FindMotor(pServ->pSics, tasMotorOrder[i]); if (!pMot) { sprintf(pBueffel, "ERROR: internal: motor %s NOT found", tasMotorOrder[i]); SCWrite(pCon, pBueffel, eError); return 0; } status = MotorGetSoftPosition(pMot, pCon, &fMot); if (status != 1) { sprintf(pBueffel, "ERROR: failed to read motor %s", tasMotorOrder[i]); SCWrite(pCon, pBueffel, eError); return 0; } amot[i] = fMot; } /* read magnet currents. */ for (i = 0; i < 4; i++) { convH[i] = .0; helmCurrent[i] = readDrivable(tasMotorOrder[CURMOT + i], pCon); helmCurrent[4 + i] = readDrivable(tasMotorOrder[CURMOT + i + 4], pCon); } readConversionFactors(self, convH); /* collect all the other machine parameters needed for a call to t_update */ 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; helm = (real) self->tasPar[HELM]->fVal; f1h = (real) self->tasPar[IF1H]->fVal; f1v = (real) self->tasPar[IF1V]->fVal; f2h = (real) self->tasPar[IF2H]->fVal; f2v = (real) self->tasPar[IF2V]->fVal; f = (real) 2.072; /* energy unit meV */ ifx = (integer) self->tasPar[FX]->iVal; /* now call t_update to do the calculation */ ier = 0; en = ef = ei = .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 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: arg = *dx * sin(*ax2 / 114.59155902616465); if(arg < .0) arg = -arg; !!!!!!!!!!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!!!!!!!!!!!!! */ if (ier != 0) { printError(ier, pCon); return 0; } /* update all the energy variables */ self->tasPar[EI]->fVal = ei; self->tasPar[KI]->fVal = aki; self->tasPar[EF]->fVal = ef; self->tasPar[KF]->fVal = akf; self->tasPar[QH]->fVal = qhkl[0]; self->tasPar[QK]->fVal = qhkl[1]; self->tasPar[QL]->fVal = qhkl[2]; /* self->tasPar[EN]->fVal = en; */ self->tasPar[EN]->fVal = ei - ef; self->tasPar[HX]->fVal = hx; self->tasPar[HY]->fVal = hy; self->tasPar[HZ]->fVal = hz; self->tasPar[QM]->fVal = qm; self->tasPar[HX]->fVal = hx; self->tasPar[HY]->fVal = hy; self->tasPar[HZ]->fVal = hz; SCparChange(pCon); /* now check the analyzer or monochromator angles */ if (ifx == 1) { diff = amot[1] / 2 - amot[0]; if (diff < .0) diff = -diff; if (diff > .1) { sprintf(pBueffel, "WARNING: A2 = %f is not half of A1 = %f", amot[1], amot[0]); SCWrite(pCon, pBueffel, eWarning); } } else if (ifx == 2) { diff = amot[5] / 2 - amot[4]; if (diff < .0) diff = -diff; if (diff > .1) { sprintf(pBueffel, "WARNING: A6 = %f is not half of A5 = %f", amot[5], amot[4]); SCWrite(pCon, pBueffel, eWarning); } } return 1; } /*----------------------------------------------------------------------- readDrivable tries to read the value of one of the magnet currents. All errors are ignored because most of the time currents will not be present in the system. -----------------------------------------------------------------------*/ float readDrivable(char *val, SConnection * pCon) { pIDrivable pDriv; CommandList *pCom; pDummy pDum; pMotor pMot = NULL; float fVal; /* else: read general drivable */ pCom = FindCommand(pServ->pSics, val); if (pCom != NULL) { pDriv = GetDrivableInterface(pCom->pData); if (pDriv != NULL) { return pDriv->GetValue(pCom->pData, pCon); } } return -999.99; } /*----------------------------------------------------------------------- startCurrent starts driving a current to a new value All errors are ignored because most of the time currents will not be present in the system. -----------------------------------------------------------------------*/ void startCurrent(char *val, SConnection * pCon, float fVal) { StartMotor(pServ->pExecutor, pServ->pSics, pCon, val, RUNDRIVE,fVal); } /*----------------------------------------------------------------------- readConversionFactors updates the array given as a parameter with the values of the Gauss->Ampere conversion factors -------------------------------------------------------------------------*/ void readConversionFactors(pTASdata self, float convH[4]) { convH[0] = self->tasPar[HCONV1]->fVal; convH[1] = self->tasPar[HCONV2]->fVal; convH[2] = self->tasPar[HCONV3]->fVal; convH[3] = self->tasPar[HCONV4]->fVal; }