TopSpin2.1
AU shearing program xfshear




Home and Applets > Pulse Program > MQMAS > TopSpin AU Xfshear

TopSpin2.1 AU program xfshear

/****************************************************************/
/*  xfshear              01.04.2005                             */
/****************************************************************/
/*  Short Description:                                          */
/*  Program used for shearing of 2D MQMAS and STMAS spectra     */
/*  of half integer quadrupolar nuclei. Data must be            */
/*  aquired in States, States/TPPI, or echo/antiecho mode       */
/*  Program used for shearing of ZQ-TROSY experiments           */
/*  like trosyzqgpphwg                                          */
/****************************************************************/
/*  Keywords:                                                   */
/*  shear, half integer spin nuclei, multi quantum, mq, MQ      */
/*  ZQ-TROSY                                                    */
/****************************************************************/
/*  Description/Usage:                                          */
/*  Program is used                                             */
/*  a) for 2D MQMAS and STMAS experiments of nuclei with        */
/*  half integer spin for shearing-FT of 2D spectrum            */
/*  including 2D FT and referencing in F1 dimension             */
/*  according to delta(MQ)=delta(iso)+(p-R)delta(qis)           */
/*                                                              */
/*  b) for ZQ-TROSY experiments to eliminate the                */
/*  contribution of proton chemical shift in the                */
/*  F1 direction                                                */
/*                                                              */
/*  When the program is started without any options it asks     */
/*  to apply abs2 and to enter a F1 shift in ppm.               */
/*  abs2 will be applied after the first xf2, the F1 shift      */
/*  allows to rotate the spectrum about the F2 axis             */
/*  in case of folded peaks.                                    */
/*                                                              */
/*  The program can be started with the following options:      */
/*  abs:    perform abs2 after F2 transform (before             */
/*          shearing operation)                                 */
/*  noabs:  do not perform abs2 after F2 transform              */
/*  rotate: no shearing is done and only the F1 shift is        */
/*          applied. This option can be useful for all          */
/*          types of States, States/TPPI, and echo/antiecho     */
/*          experiments in which case a number must be          */
/*          specified with this option.                         */
/*  ratio:  allows to shear at a different angle, which is      */
/*          defined by the parameter R.                         */
/*  lastf1: applies the same F1 shift as applied at             */
/*          previous processing.                                */
/****************************************************************/
/*  Author(s):                                                  */
/*  Name        : Christian Fernandez                           */
/*  Organisation: ENSICAEN LCS UMR6506                          */
/*  Email       : christian.fernandez@ensicaen.fr               */
/*                                                              */
/*  Name        : Stefan Steuernagel                            */
/*  Organisation: Bruker BioSpin GmbH                           */
/*  Email       : stefan.steuernagel@bruker-biospin.de          */
/****************************************************************/
/*  Name          Date         Modification:                    */
/*  wem           20050401     created                          */
/*                                                              */
/****************************************************************/
/*
$Id: xfshear,v 1.22.2.1 2007/12/21 16:34:01 wem Exp $
*/

  static const double rtab[5][10] =
/* mq 0        1    2     3        4    5        6     7        8    9         spin */
  { { 0.,     -1.,  0.,   0.,      0.,  0.,      0.,   0.,      0.,  0.   }, /* 1/2 */
    {-8./9.,  -1.,  0.,   7./9.,   0.,  0.,      0.,   0.,      0.,  0.   }, /* 3/2 */
    { 7./24., -1.,  0.,  19./12.,  0., 25./12.,  0.,   0.,      0.,  0.   }, /* 5/2 */
    {28./45., -1.,  0., 101./45.,  0., 11./9.,   0., 161./45.,  0.,  0.   }, /* 7/2 */
    {55./72., -1.,  0.,  91./36.,  0., 95./36.,  0.,   7./18.,  0., 31./6.}  /* 9/2 */
  };
  int *in2rr;
  int pparmod, bytordp, si2, si1, infile2rr, size[2], xdim[2],
      nbytes, nbytesread, xdim2, xdim1, loopcount4, i,
      mq, nspin = 0, wdw, bcmod, phmod, stsi1, stsi2, tdeff1, fnmode;
  int absflag = 0, f1flag = 0, ratioflag = 0, rotateflag = 0;
  double dra, ph1 = 0, phc2, sfo1, ratio = 0, f1shift = 0;
  float swh1;

  off_t pos = 0;
  FILE* fp;
  char  name2rr[PATH_MAX], audittext[512], nucl[80], pulprog[80];
  char* auditp = audittext + sprintf(audittext, "shear");

/* is it a 2D spectrum either States, States-TPPI or EA?
   ===================== */
  FETCHPARS("PPARMOD",&pparmod)
  if (pparmod != 1)
    STOPMSG("Program is only suitable for 2D data!");
  FETCHPAR1S("FnMODE",&fnmode)
  if (fnmode == 0)
  {
    int mc2;
    FETCHPAR("MC2",&mc2)
    if (mc2 != 3 && mc2 != 4 && mc2 != 5)
      STOPMSG("MC2 must be States, States-TPPI, or Echo-Antiecho!");
    fnmode = mc2 + 1;
  }
  else if (fnmode != 4 && fnmode != 5 && fnmode != 6)
    STOPMSG("FnMode must be in States, States-TPPI, or Echo-Antiecho!");

/* check command line options
   ========================== */
  for (i = 2; i < i_argc; i++)
  {
    char* endptr;
    auditp += sprintf(auditp, " %s", i_argv[i]);
    if (strcmp(i_argv[i], "lastf1") == 0)
      f1flag = 1;
    if (strcmp(i_argv[i], "abs") == 0)
      absflag = 1;
    if (strcmp(i_argv[i], "noabs") == 0)
      absflag = -1;
    if (strcmp(i_argv[i], "rotate") == 0)
      rotateflag = 1;
    if (strcmp(i_argv[i], "ratio") == 0)
    {
      ratioflag = 1;
      if (i + 1 < i_argc)
      {
        ratio = strtod(i_argv[i + 1], &endptr);
        if (*endptr == 0)
        {
          ratioflag = 2;
          i++;
        }
      }
    }
    if (f1flag == 0)
    {
      f1shift = strtod(i_argv[i], &endptr);
      if (*endptr == 0)
      f1flag = 2;
    }
  }
  auditp += sprintf(auditp, ":");

/* init some variables 
   =================== */
  /* carrier frequency in F2 dimension: SFO1 */
  FETCHPARS("SFO1",&sfo1)
  FETCHPARS("NUC1",nucl)
  if (strcmp(nucl,"off") == 0)
    FETCHPARS ("NUCLEUS",nucl)
  sprintf(name2rr, "%s/nmr/lists/nuclei.all", PathXWinNMRExpStan());
  fp = fopen(name2rr, "r");
  if (fp == 0)
    Proc_err(DEF_ERR_OPT, "cannot open\n%s", name2rr);
  else
  {                          /* search spin quantum number of  */
    size_t nlen = strlen(nucl);   /* nucleus in nuclei table   */
    for (;;)
    {
      char* cp;
      char* ep;
      name2rr[0] = 0;
                         /* read a single line of nuclei table */
      if (fgets(name2rr, sizeof(name2rr), fp) == 0)
        break;
      for (cp = name2rr; isspace((int)*cp); cp++)
        ;                     /* skip leading blanks in line   */
      if (strncmp(nucl, cp, nlen) == 0)
      {                               /* found nucleus in line */
         for (cp += nlen; isdigit((int)*cp) == 0; cp++)
         ;                          /* skip non digits in line */
         strtod(cp, &ep);
         cp = ep;                    /* skip Receptivity entry */
         nspin = (int)strtol(cp, &ep, 10);
         if (*ep != '/')
           nspin *= 2;                /* integer spin detected */
         break;                     /* got spin quantum number */
      }
    }
    fclose(fp);
  }
  if (nspin == 0)
  {
    double spin = 1.5;
    GETDOUBLE("Enter spin quantum number:", spin);
    auditp += sprintf(auditp, " spin quantum number %.1g", spin);
    nspin = (int)(2. * spin);
  }
  if (nspin <= 0  ||  nspin > 9  ||  nspin % 2 == 0)
    STOPMSG("Invalid spin quantum number")
  FETCHPARS("PULPROG",pulprog)
  (void)strlwr(pulprog);
  /* parameter p: mq*/
  if (strstr(pulprog,"3q"))                 mq = 3;
  else if (strstr(pulprog,"5q"))            mq = 5;
  else if (strstr(pulprog,"7q"))            mq = 7;
  else if (strstr(pulprog,"9q"))            mq = 9;
  else if (strncmp(pulprog,"stmas",5) == 0) mq = 0;
  else if (nspin == 1)                      mq = 1;
  else
  {
    mq = 3;
    GETINT("Enter the pQ order for MQMAS or enter 0 for STMAS:", mq);
    auditp += sprintf(auditp, " pQ order %d", mq);
  }
  if (mq < 0  ||  mq > (rotateflag ? 9 : nspin)  ||
     (mq == 0  &&  nspin == 1)  ||  (mq  &&  mq % 2 == 0))
    STOPMSG("Invalid pQ order")
  if (f1flag < 2  &&  mq != 1)
  {
    char ti[80];
    FETCHPARS("TI",ti)
    if (strcmp(ti, "shearing done") == 0)
    {
      float noisf1;
      FETCHPAR1S("NOISF1",&noisf1)
      f1shift = noisf1;
    }
  }

/* if lastf1 selected, don't ask for ABS2 and additional F1 shift
   ============================================================== */

  if (absflag == 0  &&  f1flag != 1  &&  mq != 1)
  {
/* automatic baseline correction for MQMAS and STMAS ? 
   =================================================== */
    char yes[80];
    (void)strcpy(yes, "yes");
    GETSTRING("Apply abs2?", yes);
    if (yes[0] == 'y')
      absflag = 1;
    else
      absflag = -1;
  }
  if (absflag == 0)
  {
    if (mq != 1)
      absflag = 1;
    else
      absflag = -1;
  }

/* additional F1 shift ?
   ===================== */
  if (f1flag == 0  &&  (mq != 1  ||  rotateflag))
    GETDOUBLE("Enter F1 shift in ppm:", f1shift);
  auditp += sprintf(auditp, " F1 shift %.7g ppm", f1shift);
  /* parameter R (position of echo): ratio */
  if (ratioflag < 2)
    ratio = rtab[nspin / 2][mq];
  if (ratioflag == 1)
    GETDOUBLE("Enter different shearing ratio:", ratio);
  if (ratioflag)
    auditp += sprintf(auditp, " shearing ratio %.15g", ratio);

/* ignore strip parameters in F1 for 1H spectra!
   ============================================= */
  FETCHPAR1("STSI",&stsi1)
  STOREPAR1("STSR",0)
  STOREPAR1("STSI",0)
  if (mq == 1)
  {
    if (stsi1 != 0)
      Proc_err(DEF_ERR_OPT,"F1 Strip-FT not yet implemented,\nparameters disabled");
  }

/* ignore all strip parameters for MQMAS spectra !
   =============================================== */
  if (mq != 1)
  {
    if (mq == nspin) mq = -mq; 
    FETCHPAR("STSI",&stsi2)
    STOREPAR("STSR",0)
    STOREPAR("STSI",0)
    if (stsi2 != 0 || stsi1 != 0)
      Proc_err(DEF_ERR_OPT,"Strip-FT not yet implemented,\nparameters disabled");
  }

/* calculate F2 transform 
   ====================== */
  FETCHPAR("WDW",&wdw)
  STOREPAR("WDW",0)
  i = CPR_exec("xf2 same raw n",WAIT_TERM);
  STOREPAR("WDW",wdw)
  if (i < 0)
    return -1;
  if (absflag > 0)
  {
    ABS2
    if (AUERR < 0)
      ABORT
  }

/* get required status parameters after F2 transform 
   ================================================= */
  FETCHPARS("BYTORDP",&bytordp)
  FETCHPARS("SI",&si2)
  FETCHPARS("XDIM",&xdim2)
  FETCHPAR1S("SWH",&swh1)
  FETCHPAR1S("SI",&si1)
  FETCHPAR1S("XDIM",&xdim1)
  FETCHPAR1S("TDeff",&tdeff1)
  nbytes = xdim2 * (2 * (int)sizeof(int));

/* calculate phase for shearing and scaling for MQMAS/STMAS F1
   =========================================================== */
  if (rotateflag == 0)
  {
    double sw2;
    FETCHPARS("SW_p",&sw2)
    ph1 = -M_PI * ratio * sw2 / swh1 / (double)si2;
  }
  /* parameter R - p: dra */
  dra = fabs(ratio - (double)mq);
/* for MQMAS */
  if (mq < 0) ph1 = -ph1;
/* for TROSY */
  if (mq == 1) dra = 1.0;
/* for STMAS of inner satellite */
  if (mq == 0) dra = fabs(ratio - 1.0);
/* for rotate */
  if (mq > nspin) dra = 1.0;
  /* carrier frequency in F1 dimension: SFO1*(R - p) */
  sfo1 *= dra;

/* additional F1 shift
   =================== */
  phc2 = -M_PI * f1shift * sfo1 / (double)swh1; 

/* compensate for strip transformation of 1H spectra
   ================================================= */
  if (mq == 1)
  {
    int stsr;
    FETCHPARS("STSR",&stsr)
    phc2 -= ph1 * stsr;
  }

/* allocate memory for I/O buffer
   ============================== */
  in2rr = (int*)malloc(nbytes);
  if (in2rr == 0)
    STOPMSG("Not enough memory");

/* Open source file RR
   =================== */
  strcpy(name2rr, PROCPATH("2rr"));
  infile2rr = open(name2rr, O_RDWR);
  if (infile2rr == -1)
  {                          /* on Windows GUI may access */
    sleep(1);                /* 2rr exclusively           */
    infile2rr = open(name2rr, O_RDWR);   /* second attempt*/
  }
  if (infile2rr == -1)
  {
    Perror(DEF_ERR_OPT, name2rr);
    return -1;
  }
  Show_status("shearing calculation");

/* Read data in submatrix 
   ====================== */
  for (loopcount1 = 0; loopcount1 < tdeff1; loopcount1 += xdim1)
    for (loopcount2 = -si2 / 2; loopcount2 < si2 / 2; loopcount2 += xdim2)
    {
      for (loopcount3 = loopcount1; loopcount3 < loopcount1 + xdim1  &&
        loopcount3 < tdeff1; loopcount3 += 2)
      {
        if ((nbytesread=read(infile2rr,in2rr,nbytes)) != nbytes) /* Read 2 rows ! */
        {
          if (nbytesread < 0)
            Perror(DEF_ERR_OPT, name2rr);
          else
            Proc_err(DEF_ERR_OPT, "%s\nread off file limits", name2rr);
          return -1;
        }
        for (loopcount4 = 0; loopcount4 < xdim2; loopcount4++)
        {
          double ev = in2rr[loopcount4];
          double od = in2rr[loopcount4 + xdim2];
          double ph = (ph1 * (loopcount2 + loopcount4) + phc2) * loopcount3;
          double co = cos(ph);
          double si = sin(ph);

/* perform shearing transformation  
   =============================== */
          in2rr[loopcount4]         = (int)(ev * co + od * si);
          in2rr[loopcount4 + xdim2] = (int)(od * co - ev * si);
        }

/* Save 
   ==== */
        if (lseek(infile2rr, pos, SEEK_SET) == -1  ||
          write(infile2rr,in2rr,nbytes) < 0)
        {
          Perror(DEF_ERR_OPT, name2rr);
          return -1;
        }
        pos += nbytes;
      }
      if (loopcount1 + xdim1 > tdeff1)
        pos += nbytes * ((loopcount1 + xdim1 - (tdeff1 - 1)) / 2);
    }
  close(infile2rr);

/* Audit file entry
   ================ */
  size[0] = si2;
  size[1] = si1;
  xdim[0] = xdim2;
  xdim[1] = xdim1;
  if (CheckSumFile(name2rr, 0, auditp, 0, 0, bytordp, 0, 2, size[0], size, xdim) > 0)
    AuditAppend(PROCPATH("auditp.txt"), audittext);

/* F1 Scaling for MQ MAS experiments
   ================================= */
  if (ratioflag == 0  &&  dra != 1.0)
  {
    double bf1;
    float  sr;
    FETCHPARS("SR",&sr)
    FETCHPARS("BF1",&bf1)
    /* BF1 in F1 dimension: BF1*(R - p) */
    bf1 *= dra;  
    STOREPAR1S("BF1",bf1)
    /* Larmor frequency of reference sample in F1 dimension: SF*(R - p) */
    bf1 += (dra * sr + f1shift * sfo1) * 1e-6;
    STOREPAR1S("SF",bf1)
    STOREPAR1S("SFO1",sfo1)
  }

 /* store TI to recall shearing and shifting 
    ======================================== */
  STOREPARS("TI","shearing done")
  STOREPAR1S("NOISF1",f1shift)

 /* inverse F2 and subsequent F2+F1 transform 
    ========================================= */
  XHT2
  if (AUERR < 0)
    ABORT
  XIF2
  if (AUERR < 0)
    ABORT
  FETCHPAR("PH_mod",&phmod)
  STOREPAR("PH_mod",0)
  FETCHPAR("BC_mod",&bcmod)
  STOREPAR("BC_mod",0)
  STOREPAR1S("AXUNIT","")
  XFB
  STOREPAR("PH_mod",phmod)
  STOREPAR("BC_mod",bcmod)
  VIEWDATA
  Show_status("shearing calculation finished");
  QUIT

(top)

[Contact me] - Last updated December 16, 2012
Solid-state NMR bibliography for
Aluminum-27
Antimony-121/123
Arsenic-75
Barium-135/137
Beryllium-9
Bismuth-209
Boron-11
Bromine-79/81
Calcium-43
Cesium-133
Chlorine-35/37
Chromium-53
Cobalt-59
Copper-63/65
Deuterium-2
Gallium-69/71
Germanium-73
Gold-197
Hafnium-177/179
Indium-113/115
Iodine-127
Iridium-191/193
Krypton-83
Lanthanum-139
Lithium-7
Magnesium-25
Manganese-55
Mercury-201
Molybdenum-95/97
Neon-21
Nickel-61
Niobium-93
Nitrogen-14
Osmium-189
Oxygen-17
Palladium-105
Potassium-39/41
Rhenium-185/187
Rubidium-85/87
Ruthenium-99/101
Scandium-45
Sodium-23
Strontium-87
Sulfur-33
Tantalum-181
Titanium-47/49
Vanadium-51
Xenon-131
Zinc-67
Zirconium-91



Copyright 2002- pascal-man.com. All rights reserved.