Modified AU shearing program xfshear for SPAM echo/antiecho MQMAS from Bruker BioSpin
Modified AU shearing program xfshear for SPAM echo/antiecho MQMAS from Bruker BioSpin
The three-pulse z-filter MQMAS pulse sequence generates echo and antiecho, which have the same amplitude with the same sign. In contrast, the SPAM echo/antiecho pulse program generates echo amplitude and antiecho amplitude with opposite signs. Furthermore, the first sequence involves States acquisition procedure whereas the second sequence involves the echo/antiecho procedure. As a result, the shearing AU program should be modified accordingly.
/****************************************************************/
/* xfshear2005 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! */
/* ste 031031 shearing and calibration */
/* of STMAS included */
/* "xf2 raw" command used to avoid */
/* error if partly processed data */
/* exist */
/* ste 040112 check pparmod for 2D, so it can */
/* be applied to planes from 3D */
/* echo/anti-echo treatment incl. */
/* ste 041224 auto function included to */
/* suppress abs2 and F1-shift query */
/* dynamic memory allocation incl. */
/****************************************************************/
/*
$Id: xfshear,v 1.15.2.2 2005/02/07 14:06:38 es Exp $
*/
#include <math.h>
int parmode, si2, si1, outexpno, tempfile2rr, infile2rr,
tempfile2ir, infile2ir;
int *in2rr, *in2ir, *temp2rr, *temp2ir, *temp2rr2,
*temp2ir2, *temp2rrB, *temp2irB, *temp2rrB2,
*temp2irB2;
int dummy, nbytes, sizeofint, nbytesread, xdim2, xdim1, loopcount4,
sig, mq, nspin, wdw2, bcmod, phmod, stsr2, stsi1, stsi2, fnmode, mc2;
int option1=0, option2=0;
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],
imfile[8],
outname[PATH_LENGTH], nucl[80], yes[20],
pulprog[80], ti[80], ShearOpt[PATH_MAX];
/* select current dataset
====================== */
GETCURDATA;
(void) strcpy (yes,"yes");
/* is-it a 2D spectrum and taken in States, States-TPPI or EA?
===================== */
FETCHPARS("PPARMOD",&parmode);
if ( parmode != 1 )
{ STOPMSG("Program is only suitable for 2D data!");}
FETCHPAR1S("FnMODE",&fnmode);
if ( fnmode == 0 ) {
FETCHPAR("MC2",&mc2);
if ( mc2 != 3 && mc2 != 4 && mc2 != 5 )
STOPMSG("mc2 must be States, States-TPPI, or Echo-Antiecho!");
}
else if ( fnmode != 4 && fnmode != 5 && fnmode != 6)
STOPMSG("FnMode must be in States, States-TPPI, or Echo-Antiecho!");
/* prepare for automatic execution and other ratio
=============================================== */
sprintf(ShearOpt,"%s",cmd);
if ( strstr(ShearOpt,"auto")) option1 = 1;
if ( strstr(ShearOpt,"ratio")) option2 = 1;
/* if auto selected, don't aks for ABS2 and additional F1 shift
============================================================ */
if (option1 == 1)
{
f1=0.;
FETCHPAR1S("NOISF1",&Noise);
FETCHPARS("TI",ti);
if (!strcmp(ti,"shearing done")) f1=Noise;
}
if (option1 == 0)
{
/* automatic baseline correction for MQMAS and STMAS ?
=================================================== */
if (mq != 1) GETSTRING("Apply ABS2 ? : ",yes);
/* 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);
}
/* 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,"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,"209Bi")) spin=4.5;
else if(!strcmp(nucl,"1H")) spin=1.0;
else GETFLOAT(" Enter spin number : ",spin);
nspin = (int) ceil ( (double)(2*spin) );
if ( nspin == 3 && !strncmp(pulprog,"mp3q",4) )
{ 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 if(!strncmp(pulprog,"stmas",5))
{
mq=0;
if(nspin == 3) ratio=-8./9.;
if(nspin == 5) ratio=7./24.;
if(nspin == 7) ratio=28./45.;
if(nspin == 9) ratio=55./72.;
}
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.; } ;
if ( option2 == 1) GETDOUBLE(" Enter different shearing ratio :",ratio);
/* 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");
}
/* calculate F2 transform
====================== */
SETCURDATA AUERR=CPR_exec( "xf2 raw",WAIT_TERM);
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;
/* calculate phase for shearing and scaling for MQMAS/STMAS F1
=========================================================== */
ph1= -2.*3.141529*ratio*sw2 / swh1 / (double)(si2);
/* for MQMAS */
if ( mq<0 ) ph1= -ph1;
d1= fabs(ratio - (double)mq);
/* for TROSY */
if ( mq==1 ) d1=1.0;
/* for STMAS of inner satellite */
if ( mq==0 ) d1=fabs(ratio - 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;
}
/* decide 2ii or 2ir file
====================== */
if ( fnmode == 1 ) (void)sprintf(imfile,"2ii");
else (void)sprintf(imfile,"2ir");
/* 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/%s",
disk,user,name,expno,procno,imfile);
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);
/* allocate memory for array */
in2rr = (int*)malloc(nbytes*10);
if (in2rr == 0) STOPMSG("Not enough memory");
in2ir = in2rr + nbytes/sizeofint;
temp2rr= in2ir + nbytes/sizeofint;
temp2ir= temp2rr + nbytes/sizeofint;
temp2rr2= temp2ir + nbytes/sizeofint;
temp2ir2= temp2rr2 + nbytes/sizeofint;
temp2rrB= temp2ir2 + nbytes/sizeofint;
temp2irB= temp2rrB + nbytes/sizeofint;
temp2rrB2= temp2irB + nbytes/sizeofint;
temp2irB2= temp2rrB2 + nbytes/sizeofint;
/* 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 */
for (loopcount4=0;loopcount4<xdim2;loopcount4++)
{ temp2rr[loopcount4]= in2rr[loopcount4];
temp2ir[loopcount4]= in2ir[loopcount4];}
/* Store iSy */
for (loopcount4=0;loopcount4<xdim2;loopcount4++)
{ temp2rr2[loopcount4]= in2rr[loopcount4+xdim2];
temp2ir2[loopcount4]= in2ir[loopcount4+xdim2];}
/* combine Sx and iSy to create ...
================================ */
for (loopcount4=0;loopcount4<xdim2;loopcount4++)
{
/* echo
==== */
temp2rrB[loopcount4]=(temp2rr[loopcount4]
-temp2ir2[loopcount4])/2;
temp2irB[loopcount4]=(temp2ir[loopcount4]
+temp2rr2[loopcount4])/2;
/* and antiecho
============ */
temp2rrB2[loopcount4]=(temp2rr[loopcount4]
+temp2ir2[loopcount4])/2;
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;
in2rr[loopcount4]=
(int)(((double)temp2rrB[loopcount4])*cos(phc1)
+((double)temp2irB[loopcount4])*sin(phc1))/2;
in2ir[loopcount4]=
(int)(((double)temp2irB[loopcount4])*cos(phc1)
-((double)temp2rrB[loopcount4])*sin(phc1))/2;
in2rr[loopcount4+xdim2]=
(int)(((double)temp2rrB2[loopcount4])*cos(phc1)
-((double)temp2irB2[loopcount4])*sin(phc1))/2;
in2ir[loopcount4+xdim2]=
(int)(((double)temp2irB2[loopcount4])*cos(phc1)
+((double)temp2rrB2[loopcount4])*sin(phc1))/2;
}
for (loopcount4=0;loopcount4<xdim2;loopcount4++)
{temp2rr[loopcount4]=0.5*(
in2rr[loopcount4]+in2rr[loopcount4+xdim2]);
temp2ir[loopcount4]=0.5*(
in2ir[loopcount4]+in2ir[loopcount4+xdim2]);
temp2irB[loopcount4]=0.5*(
in2rr[loopcount4]-in2rr[loopcount4+xdim2]);
temp2rrB[loopcount4]=0.5*(
in2ir[loopcount4]-in2ir[loopcount4+xdim2]);
}
for (loopcount4=0;loopcount4<xdim2;loopcount4++)
{in2rr[loopcount4]= temp2rr[loopcount4];
in2ir[loopcount4]= temp2ir[loopcount4];
in2rr[loopcount4+xdim2]= temp2rrB[loopcount4];
in2ir[loopcount4+xdim2]= temp2irB2[loopcount4];
}
/* 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 ( option2 == 0) {
if ( d1 != 1.0 ) {
FETCHPARS("BF1",&bf12);
FETCHPARS("SR",&sr2);
bf11= d1*bf12;
sr1= (float)(d1*sr2);
STOREPAR1S("BF1",bf11);
STOREPAR1("BF1",bf11);
#if 1
STOREPAR1S("SF",bf11 + sr1 * 1e-6);
STOREPAR1("SF",bf11 + sr1 * 1e-6);
#endif
#if 0
STOREPAR1S("SR",sr1 );
STOREPAR1("SR",sr1);
#endif
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)
CPR_exec("xfb nc_proc 0",WAIT_TERM);
STOREPAR("PH_mod",phmod)
STOREPAR("BC_mod",bcmod)
VIEWDATA
QUITMSG(" Shearing-FT done!")