AU shearing program xfshear for z-filter MQMAS
AU shearing program xfshear for z-filter MQMAS
/****************************************************************/
/* xfshear 08.11.1996 */
/****************************************************************/
/* Short Description: */
/* */
/* Program for shearing of 2D MQMAS spectra of odd half */
/* integer quadrupolar nuclei. Data need to be aquired in */
/* States Mode */
/* Program used for shearing of ZQ-TROSY experiments */
/* like trosyzqgpphwg */
/****************************************************************/
/* Keywords: */
/* */
/* shear, odd half spin nuclei, multi quantum, mq, MQ */
/* ZQ-TROSY */
/****************************************************************/
/* Description/Usage: */
/* */
/* Program is used */
/* */
/* a) for 2D MQ experiments on nuclei with */
/* odd half integer spin for shearing-FT of 2D spectrum */
/* including 2DFT and referencing in F1 dimension */
/* according to delta(MQ)=delta(iso)+(p-R)delta(qis) */
/* It checks for the parameter "nucleus" (AMX) or "NUC1" */
/* (Avance) and sets the spin value accordingly. */
/* If the nucleus is not in the list it asks for the spin, */
/* too. Only most common nuclei are in the list. */
/* If the pulse program name matches mpxq, where */
/* x corresponds to the multiple quantum order of */
/* the experiment, then this is assumed, otherwise */
/* it is asked for by the program. */
/* Program asks if abs2 in F2 dimension is desired. */
/* If answered with yes, then a Hilbert transform */
/* is automatically done, which requires xwinnmr1.3 */
/* or later. In older versions, this command has to */
/* be eliminated from the program code. */
/* A F1 frequency shift in ppm is asked for to compensate */
/* for malchosen o1 setting. This value is stored and */
/* suggested for any repeated processing. */
/* */
/* b) for ZQ-TROSY experiments to eliminate the */
/* contribution of proton chemical shift in the */
/* F1 direction */
/****************************************************************/
/* Author(s): */
/* */
/* Name : Christian Fernandez */
/* Organisation : Universite de Lille, CNRS-801 */
/* Email : christian.fernandez@univ.lille.fr */
/* */
/* Name : Stefan Steuernagel */
/* Organisation : Bruker Analytik */
/* Email : stefan.steuernagel@bruker.de */
/****************************************************************/
/* Name Date Modification: */
/* */
/* ste 960201 created */
/* ste 980813 check xdim's and status sizes */
/* at right time */
/* ste 980813 XHT2 included after XF1 to */
/* allow f2 phase correction after */
/* processing */
/* ste 981101 close of input data files and */
/* "mv temp-file to 2rr" etc */
/* included */
/* ste 000105 XIF2 and XFB included after */
/* shearing to allow apodisation */
/* in F2, set-up to be done prior */
/* to execution of xfshear, all */
/* handling done by this program */
/* ber/ste 010909 handling for ZQ-TROSY experiment */
/* included, allows to use strip */
/* transformation in F2 only */
/* not strip-FT for MQMAS! */
/****************************************************************/
/*
$Id: xfshear,v 1.9.2.1 2001/09/27 14:46:43 ste Exp $
*/
#include <math.h>
int parmode, si2, si1, outexpno, tempfile2rr, infile2rr,
tempfile2ir, infile2ir, in2rr[4096], in2ir[4096],
temp2rr[4096], temp2ir[4096], temp2rr2[4096], temp2ir2[4096],
temp2rrB[4096], temp2irB[4096], temp2rrB2[4096], temp2irB2[4096],
nbytes, sizeofint, nbytesread, xdim2, xdim1, loopcount4,
sig, mq, nspin, wdw2, bcmod, phmod, stsr2, stsi1, stsi2;
double sw1, sw2, ratio, ph, ph1, ph2, ph3, phc1, phc2,
in0, sfo1, sfo12, sfo11, bf11, bf12;
float Noise, spin, offset1, offset2, off2, off1, sr2, sr1, swh1;
char inname2rr[PATH_LENGTH], tempname2rr[PATH_LENGTH],
inname2ir[PATH_LENGTH], tempname2ir[PATH_LENGTH],
outname[PATH_LENGTH], nucl[80], yes[20],
pulprog[80], ti[80];
/* select current dataset
====================== */
GETCURDATA;
(void) strcpy (yes, "yes");
/* is-it a 2D spectrum ?
===================== */
FETCHPAR("PARMODE", &parmode);
if ( parmode != 1 )
{ STOPMSG(" Not a 2D spectrum!"); }
/* init some variables
=================== */
FETCHPARS("PULPROG", pulprog);
FETCHPARS("SFO1", &sfo1);
FETCHPAR("WDW", &wdw2)
STOREPAR("WDW", 0)
/* shearing ratio calculation for quadrupolar nuclei
set mq and spin to 1 for 1H
================================================= */
spin = 5./2.;
mq = 3;
FETCHPARS("NUC1", nucl)
if (!strcmp(nucl, "off")) FETCHPARS ("NUCLEUS", nucl)
if (!strcmp(nucl, "7Li")) spin = 1.5;
else if (!strcmp(nucl, "11B")) spin = 1.5;
else if (!strcmp(nucl, "17O")) spin = 2.5;
else if (!strcmp(nucl, "23Na")) spin = 1.5;
else if (!strcmp(nucl, "27Al")) spin = 2.5;
else if (!strcmp(nucl, "35Cl")) spin = 1.5;
else if (!strcmp(nucl, "39K")) spin = 1.5;
else if (!strcmp(nucl, "45Sc")) spin = 3.5;
else if (!strcmp(nucl, "45Sc")) spin = 3.5;
else if (!strcmp(nucl, "51V")) spin = 3.5;
else if (!strcmp(nucl, "55Mn")) spin = 2.5;
else if (!strcmp(nucl, "63Cu")) spin = 1.5;
else if (!strcmp(nucl, "65Cu")) spin = 1.5;
else if (!strcmp(nucl, "69Ga")) spin = 1.5;
else if (!strcmp(nucl, "71Ga")) spin = 1.5;
else if (!strcmp(nucl, "85Rb")) spin = 2.5;
else if (!strcmp(nucl, "87Rb")) spin = 1.5;
else if (!strcmp(nucl, "93Nb")) spin = 4.5;
else if (!strcmp(nucl, "121Sb")) spin = 2.5;
else if (!strcmp(nucl, "1H")) spin = 1.0;
else GETFLOAT(" Enter spin number : ", spin);
nspin = (int) ceil( (double)(2*spin) );
if ( nspin == 3 )
{ mq = -3;
ratio = 7./9.; }
else
{ if (!strncmp(pulprog, "mp3q", 4)) mq = 3;
else if(!strncmp(pulprog, "mp5q", 4)) mq = 5;
else if(!strncmp(pulprog, "mp7q", 4)) mq = 7;
else if(!strncmp(pulprog, "mp9q", 4)) mq = 9;
else if(!strcmp(pulprog, "trosyzqgpphwg")) mq = 1;
else GETINT(" Enter the pQ order : ", mq);
mq = (int) ceil( (double)mq );
if ( mq == nspin ) mq = -mq;
if ( nspin == 5 && mq == -5) ratio = 25./12.;
if ( nspin == 5 && mq == 3) ratio = 19./12.;
if ( nspin == 7 && mq == -7) ratio = 161./45.;
if ( nspin == 7 && mq == 5) ratio = 11./9.;
if ( nspin == 7 && mq == 3) ratio = 101./45.;
if ( nspin == 9 && mq == -9) ratio = 31./6.;
if ( nspin == 9 && mq == 7) ratio = 7./18.;
if ( nspin == 9 && mq == 5) ratio = 95./36.;
if ( nspin == 9 && mq == 3) ratio = 91./36.;
if ( mq == 1 ) ratio = -1.; };
/* ignore strip parameters in F1 for 1H spectra!
============================================= */
if ( mq == 1 )
{ FETCHPAR1("STSI", &stsi1);
STOREPAR1("STSR", 0);
STOREPAR1("STSI", 0);
if ( stsi1 != 0 )
Proc_err(0, "F1 Strip-FT not yet implemented, \nparameters disabled"); }
/* ignore all strip parameters for MQMAS spectra !
=============================================== */
if ( mq != 1 )
{ FETCHPAR("STSI", &stsi2);
FETCHPAR1("STSI", &stsi1);
STOREPAR("STSR", 0);
STOREPAR("STSI", 0);
STOREPAR1("STSR", 0);
STOREPAR1("STSI", 0);
if ( stsi2 != 0 || stsi1 != 0 )
Proc_err(0, "Strip-FT not yet implemented, \nparameters disabled"); }
/* 1st FFT : acquisition must be done in STATES or
echo/antiecho mode!
and automatic baseline correction ?
=============================================== */
if ( mq != 1 ) GETSTRING("Apply ABS2 ? : ", yes);
STOREPAR1("MC2", 3);
STOREPAR1S("MC2", 3);
/* additional F1 shift ?
===================== */
f1 = 0.;
FETCHPAR1S("NOISF1", &Noise);
FETCHPARS("TI", ti);
if ( !strcmp(ti,"shearing done") ) f1 = Noise;
if ( mq != 1 ) GETFLOAT(" F1 shift in ppm ? : ", f1);
/* calculate F2 transform
====================== */
XF2
if ( yes[0] == 'y' ) {ABS2 XHT2}
/* get required status parameters after F2 transform
================================================= */
FETCHPARS("SW_p", &sw2);
FETCHPARS("SI", &si2);
FETCHPARS("STSR", &stsr2);
FETCHPAR1S("SWH", &swh1);
FETCHPAR1S("SI", &si1);
FETCHPAR1S("XDIM", &xdim1);
FETCHPARS("XDIM", &xdim2);
sizeofint = sizeof(int);
nbytes = sizeofint*xdim2*2;
/* phase for shearing
================== */
ph1 = -2.*3.141529*ratio*sw2/swh1/(double)(si2);
if ( mq < 0 ) ph1 = -ph1;
d1 = fabs(ratio - (double)mq);
if ( mq == 1 ) d1 = 1.0;
sfo11 = d1*sfo1;
/* additional F1 shift
=================== */
ph2 = (double)f1;
phc2 = -2.*3.141529*ph2*sfo11/(double)swh1;
/* compensate for strip transformation of 1H spectra
================================================= */
if ( mq == 1 )
{ ph3 = -ph1*stsr2;
phc2 += ph3; }
/* Open source file RR and IR
========================== */
(void)sprintf( inname2rr, "%s/data/%s/nmr/%s/%d/pdata/%d/2rr",
disk, user, name, expno, procno);
if ( (infile2rr = open(inname2rr,0)) == -1 )
{ (void)sprintf(text, " I/O Error (Open) \n%s ", inname2rr);
STOPMSG(text); }
(void)sprintf( inname2ir, "%s/data/%s/nmr/%s/%d/pdata/%d/2ir",
disk, user, name, expno, procno);
if ( (infile2ir = open(inname2ir,0)) == -1 )
{ (void)sprintf(text, " I/O Error (Open) \n%s ", inname2ir);
STOPMSG(text); }
/* Create temporary files : 2rrtemp et 2irtemp
=========================================== */
(void)sprintf( tempname2rr, "%s/data/%s/nmr/%s/%d/2rrtemp",
disk, user, name, expno);
if ( (tempfile2rr = creat(tempname2rr,0664)) == -1 )
{ (void)sprintf(text, " I/O Error (Create) \n%s ", tempname2rr);
STOPMSG(text); }
(void)sprintf( tempname2ir, "%s/data/%s/nmr/%s/%d/2irtemp",
disk, user, name, expno);
if ( (tempfile2ir = creat(tempname2ir,0664)) == -1 )
{ (void)sprintf(text, " I/O Error (Create) \n%s ", tempname2ir);
STOPMSG(text); }
(void)sprintf(text, "shear : calculating");
Show_status(text);
/* Read data in submatrix
====================== */
TIMES(si1/xdim1)
TIMES2(si2/xdim2)
TIMES3(xdim1/2)
if ((nbytesread = read(infile2rr, in2rr, nbytes)) <= 0)
/* Read 2 rows ! */
{ (void)remove(tempname2rr);
STOPMSG(" 2rr file corrupted! "); }
if ((nbytesread = read(infile2ir, in2ir, nbytes)) <= 0)
{ (void)remove(tempname2ir);
STOPMSG(" 2ir file corrupted! "); }
/* Store Sx, I = 5/2, mq = 3 */
for (loopcount4 = 0; loopcount4 < xdim2; loopcount4++)
{ temp2rr[loopcount4] = in2rr[loopcount4]; /* Re[Sx] */
temp2ir[loopcount4] = in2ir[loopcount4]; } /* Im[Sx] */
/* Store iSy, I = 5/2, mq = 3 */
for (loopcount4 = 0; loopcount4 < xdim2; loopcount4++)
{ temp2rr2[loopcount4] = in2rr[loopcount4 + xdim2]; /* Re[Sy] */
temp2ir2[loopcount4] = in2ir[loopcount4 + xdim2]; } /* Im[Sy] */
/* combine Sx and iSy to create
============================ */
for (loopcount4 = 0; loopcount4 < xdim2; loopcount4++)
{
/* echo, I = 5/2, mq = 3 */
/* Re[echo] = 0.5*(Re[Sx] - Im[Sy]) */
temp2rrB[loopcount4] = (temp2rr[loopcount4]
- temp2ir2[loopcount4])/2;
/* Im[echo] = 0.5*(Im[Sx] + Re[Sy]) */
temp2irB[loopcount4] = (temp2ir[loopcount4]
+ temp2rr2[loopcount4])/2;
/* and antiecho, I = 5/2, mq = 3 */
/* Re[antiecho] = 0.5*(Re[Sx] + Im[Sy]) */
temp2rrB2[loopcount4] = (temp2rr[loopcount4]
+ temp2ir2[loopcount4])/2;
/* Im[antiecho] = 0.5*(Im[Sx] - Re[Sy]) */
temp2irB2[loopcount4] = (temp2ir[loopcount4]
- temp2rr2[loopcount4])/2; }
/* perform shearing transformation
=============================== */
i1 = loopcount3 + (loopcount1)*xdim1/2;
for (loopcount4 = 0; loopcount4 < xdim2; loopcount4++)
{ i2 = loopcount4 + (loopcount2)*xdim2;
phc1 = ph1*(double)((i2 - si2/2)*i1);
phc1 += phc2*(double)i1;
/* Sheared echo, I = 5/2, mq = 3 */
/* Re[sheared-echo] = 0.5*(Re[echo]*cos phc1
+ Im[echo]*sin phc1) */
in2rr[loopcount4] =
(int)(((double)temp2rrB[loopcount4])*cos(phc1)
+ ((double)temp2irB[loopcount4])*sin(phc1))/2;
/* Im[sheared-echo] = 0.5*(Im[echo]*cos phc1
- Re[echo]*sin phc1) */
in2ir[loopcount4] =
(int)(((double)temp2irB[loopcount4])*cos(phc1)
- ((double)temp2rrB[loopcount4])*sin(phc1))/2;
/* Sheared antiecho, I = 5/2, mq = 3 */
/* Re[sheared-antiecho] = 0.5*(Re[antiecho]*cos phc1
- Im[antiecho]*sin phc1) */
in2rr[loopcount4 + xdim2] =
(int)(((double)temp2rrB2[loopcount4])*cos(phc1)
- ((double)temp2irB2[loopcount4])*sin(phc1))/2;
/* Im[sheared-antiecho] = 0.5*(Im[antiecho]*cos phc1
+ Re[antiecho]*sin phc1) */
in2ir[loopcount4 + xdim2] =
(int)(((double)temp2irB2[loopcount4])*cos(phc1)
+ ((double)temp2rrB2[loopcount4])*sin(phc1))/2; }
for (loopcount4 = 0; loopcount4 < xdim2; loopcount4++)
{
/* Sheared Sx, I = 5/2, mq = 3 */
/* Re[sheared-Sx] = 0.5*(Re[sheared-echo] + Re[sheared-antiecho]) */
temp2rr[loopcount4] = 0.5*(
in2rr[loopcount4] + in2rr[loopcount4 + xdim2]);
/* Im[sheared-Sx] = 0.5*(Im[sheared-echo] + Im[sheared-antiecho]) */
temp2ir[loopcount4] = 0.5*(
in2ir[loopcount4] + in2ir[loopcount4 + xdim2]);
/* Sheared Sy, I = 5/2, mq = 3 */
/* -Im[sheared-Sy] = 0.5*(Re[sheared-echo] - Re[sheared-antiecho]) */
temp2irB[loopcount4] = 0.5*(
in2rr[loopcount4] - in2rr[loopcount4 + xdim2]);
/* Re[sheared-Sy] = 0.5*(Im[sheared-echo] - Im[sheared-antiecho]) */
temp2rrB[loopcount4] = 0.5*(
in2ir[loopcount4] - in2ir[loopcount4 + xdim2]); }
for (loopcount4 = 0; loopcount4 < xdim2; loopcount4++)
{ /** Real part for XF1 **/
in2rr[loopcount4] = temp2rr[loopcount4]; /*Re[sheared-Sx]*/
in2ir[loopcount4] = temp2ir[loopcount4]; /*Im[sheared-Sx]*/
/** Imaginary part for XF1 **/
in2rr[loopcount4 + xdim2] = temp2rrB[loopcount4];/*Re[sheared-Sy]*/
in2ir[loopcount4 + xdim2] = temp2irB2[loopcount4]; } /*Im[sheared-antiecho]*/
/* Save
==== */
write(tempfile2rr, in2rr, nbytesread);
write(tempfile2ir, in2ir, nbytesread);
END; /* loop3 */
END; /* loop2 */
END; /* loop1 */
close(tempfile2rr);
close(tempfile2ir);
close(infile2rr);
close(infile2ir);
/* Copy
==== */
(void)remove(inname2rr);
if (rename(tempname2rr, inname2rr) != 0) Proc_err(1, "error %d", errno);
(void)remove(inname2ir);
if (rename(tempname2ir, inname2ir) != 0) Proc_err(1, "error %d", errno);
(void)sprintf(text, "shear : finished");
Show_status(text);
/* F1 Scaling for MQ MAS experiments
================================= */
if ( d1 != 1.0 )
{ FETCHPARS("BF1", &bf12);
FETCHPARS("SR", &sr2);
bf11 = d1*bf12;
sr1 = (float)(d1*sr2);
STOREPAR1S("BF1", bf11);
STOREPAR1("BF1", bf11);
STOREPAR1S("SR", sr1);
STOREPAR1("SR", sr1);
STOREPAR1S("SFO1", sfo11);
STOREPAR1("SFO1", sfo11);
FETCHPAR1S("OFFSET", &offset1);
offset1 = offset1 - ph2;
STOREPAR1S("OFFSET", offset1); }
/* store TI to recall shearing and shifting
======================================== */
STOREPARS("TI", "shearing done");
STOREPAR1S("NOISF1", f1);
/* inverse F2 and subsequent F2+F1 transform
========================================= */
XHT2
XIF2
FETCHPAR("PH_mod", &phmod)
STOREPAR("PH_mod", 0)
STOREPAR("WDW", wdw2)
FETCHPAR("BC_mod", &bcmod)
STOREPAR("BC_mod", 0)
XFB
STOREPAR("PH_mod", phmod)
STOREPAR("BC_mod", bcmod)
/* Do not delete these two lines */
VIEWDATA
QUITMSG(" Shearing-FT done!")
Before executing the AU program xfshear, the first spectrum from the SER file should be phased like a positive amplitude absorption spectrum.