/***********************************************************************/
/*  DSSref              7.12.2009       derived from                   */
/*	calc_wishart.dmo		12.12.2002                                     */
/***********************************************************************/
/*	Short Description :					                                       */
/*	AU program for referencing indirect dimensions 		                 */
/*	according to chemical shift ratio for DSS                          */
/*	as recommended by IUPAC for proteins:			                         */
/*	Wishart, et al., J. Biomol. NMR 6, 135-140 (1995).	               */
/***********************************************************************/
/*	Keywords :						                                             */
/***********************************************************************/
/*	Usage :	1. Calibrate 1D of sample containing DSS by	               */
/*                 setting DSS signal to 0 ppm. 		                   */
/*              2. Note the value of SR in the 1D.		                 */
/*              3. Execute this script with either                     */
/*                  'DSSref' return	and enter SR in dialog	           */
/*                     or                                              */
/*                  'DSSref <SR>'   where SR is the value from 1D      */
/*                     this 2nd format suppresses output               */
/*		4. Script sets SR & SF values of indirect dims	                 */
/*              5. After running script, reprocess the spectrum        */
/*                                                                     */
/*  NOTE: does not work on 2D projections of a 3D		                   */
/*       for projections of a 3D/ND				                             */
/*        manually copy appropriate SR values from		                 */
/*       the calibrated 3D or ND into indirect dim of		               */
/*       the projection						                                     */
/*                                                                     */
/*  This script is a modification of the script calc_wishart.dmo       */
/*  and not an official script from the Bruker library of AU scripts   */
/*                                                                     */
/*  DISCLAIMER:  This script is provided as is, without warranty       */
/*              of any kind, expressed or implied, including but not   */
/*              limited to the warranties of merchantability, fitness  */
/*              for a particular purpose and noninfringement.          */
/*              The authors shall not be held responsible for any      */
/*              claim, damages, or liability, arising in connection    */
/*              with use of this script.                               */
/***********************************************************************/
/*	Author(s) :						                                             */
/*	Name		: Ruediger Weisemann 			                                 */
/*	Organisation	: Bruker Biospin GmbH			                           */
/*                                                                     */
/*	Name		: Fred F. Damberger			                                   */
/*	Organisation	: ETH Zurich				                                 */
/***********************************************************************/
/*  Name  Date	 Modification:			                                   */
/*  rwe	  021212  created				                                       */
/*  dmo   090712  sr added to calculation			                         */
/*  ffd   091207  added ratio for 31P				                           */
/*                added ratio for 19F 				                         */
/*                added Usage description			                         */
/*        121029  fixed bug in TS3, rewrote to use SR (not SF)         */
/*        140614  added ratio for 2H and modified output info	         */
/*                correct for SI dependent error in calibration	       */
/*        160107  bug fix, now uses SR entered by user for 2D & 3D     */
/*                               and uses nuc1 to calibrate F2 in 2D   */
/*        160125  revised so it uses bf1 in direct dim for H1FREQMHZ   */
/*                only calibrates if direct dim nuc1 is 1H             */
/*        160404  revised to allow commandline entry of SR             */
/*        180308  revised so it looks for NUC=H1 to get H1FREQMHZ      */
/*                from the corresponding bf, can now cal if bf1 not H1 */
/*        201120  revised to remove 1/2 point correction               */
/*                (not needed in TS4.0.7)                              */
/*        210404  simplified and introduced ability to apply 1/2 point */
/*                correction if Topspin is below a certain version     */
/*                Topspin version detection did not work               */
/*                For old topspin versions with half pixel error:      */
/*                You must manually modify script                      */
/*                to set pixelcorrection variable:                     */
/*                  see commented out lines containing pixelcorrection */
/***********************************************************************/
/* comment in next line for full debugging output */
/* #define DEBUG */

#define MSGSTART "DSSref (Version 2020-04-08 -FFD) updated SR values.\n\nUSAGE: \nDSSref with no argument asks user for SR value\nDSSref <SR>  with a single argument equal to SR value does not display this window.\n\nAfter running DSSref it may be necessary to reprocess spectrum in order for calibrations to update.\n\nChi values used:"

/* chemical shift ratios for DSS */
#define NRATIO 0.101329118          
#define CRATIO 0.251449530          /* Wishart, et al. (1995) see NREF below */
/*#define PRATIO 0.404807370 */     /* Edzes H (1992) */
#define PRATIO 0.404807356          /* Maurer & Kalbitzer (1996) */
/*#define PRATIO 0.404808636  */    /* Wishart, et al. (1995) see CREF below */
#define FRATIO 0.940866982          /* Maurer, T. and Kalbitzer (1996) see FREF below */
#define H2RATIO 0.153506088         /* Markley et al. (1995) see H2REF below */

/* literature references for ratios */
#define NREF "Wishart et al., J. Biomol. NMR 6, 135-140 (1995)."
#define CREF "Wishart et al., J. Biomol. NMR 6, 135-140 (1995)."
/*#define PREF "Markley et al., Pure & Appl. Chem. 70,117-142 (1995). " */
#define PREF "Maurer, T. and Kalbitzer, H. R., J. Magn. Reson. B 113, 177-178 (1996): 85% phosphoric acid in a sphere."
/*#define PREF "Edzes, H. Magn. Reson. Chem. 30, 850-854 (1992): 85% phosphoric acid in standard tube at 25C (TMS dispersed)."*/
#define FREF "Maurer, T. and Kalbitzer, H. R., J. Magn. Reson. B 113, 177-178 (1996): trifluoroacetic acid in a sphere."
#define H2REF "Markley et al., Pure & Appl. Chem. 70,117-142 (1995)."
#define REFNOTE "All values obtained from BMRB table at http://www.bmrb.wisc.edu/ref_info/cshift.shtml"
#define REFNOTE2 "Script does not work on projections of 3Ds\nCalibrate 3D before generating projections."
#define REFNOTE3 "After running script, repeat processing of spectrum to update ppm scales."

char pul[40],nuc1[8],nuc2[8],nuc3[8],nuc4[8],nuc5[8],nuc6[8],TempString[256],message[196];
double factor1,factor2,factor3,grad1,grad2,grad3,grad4;
double bf1,bf2,bf3,bf4,bf5,bf6,bfzero;
double sf1, sf2, sf3, sf4, sf5, sf6, sf;
double H1FREQMHZ;
float sr1, sr2, sr3, sr=0.0,srnew;
float hzppt1,hzppt2,hzppt3,pixelcorrection;
float swh1,swh2,swh3;
int si1,si2,si3;
/* double sr1,sr2,sr3,sr,srnew; */
int dim,quiet=0;
int TSversionMajor,TSversionMinor, TSversionPatch;
float TSversion;
float TSversionWhereHalfPointReferencingErrorWasFixed=4.0;
/* sadly, the section below had to be commented out because the commands
  GetTsVersion .. were only created in TS4, whereas in earlier TS versions they are named getxwinvers(....)
   so there is no way to make a script which checks the topspin version, how clever you are bruker!
TSversionMajor = GetTsVersionMajor();
TSversionMinor = GetTsVersionMinor();
TSversionPatch = GetTsVersionPl();

TSversion = TSversionMajor + TSversionMinor/10.0 + TSversionPatch/100.0;

#ifdef DEBUG
Proc_err(DEF_ERR_OPT,"Topspin version as floating point:  %f major: %i minor: %i patchlevel: %i \n", TSversion, TSversionMajor,TSversionMinor,TSversionPatch);
#endif

*/

GETCURDATA

FETCHPARS("PARMODE", &dim)
dim++;

/* if cmd line has at least 2 arguments, evaluate 2nd argument of cmd line. E.g. 'DSSref 82.83' , 2nd argument = 82.83 */
    if (i_argc > 2)
    {
      sscanf(i_argv[2],"%f",&sr);
      quiet=1;
    }
    else
    {
GETDOUBLE("Enter SR value obtained by calibrating 1H 1D to DSS with a very large SI value e.g. 64k",sr);
    }

/* insert section here to ensure that sr has a value - found no way to check whether user entered text or not, no text seems to set sr=0 */

#ifdef DEBUG
Proc_err(DEF_ERR_OPT,"data dimension %d", dim);
#endif
/* Get the spectrometer 1H frequency */
FETCHPAR("NUC1", nuc1);
FETCHPAR("BF1", &bf1);
FETCHPAR("NUC2", nuc2);
FETCHPAR("BF2", &bf2);

FETCHPAR("NUC3", nuc3);
FETCHPAR("BF3", &bf3);

FETCHPAR("NUC4", nuc4);
FETCHPAR("BF4", &bf4);

FETCHPAR("NUC5", nuc5);
FETCHPAR("BF5", &bf5);

FETCHPAR("NUC6", nuc6);
FETCHPAR("BF6", &bf6);

if (strcmp(nuc1, "1H") == 0) H1FREQMHZ = bf1;
else if (strcmp(nuc2, "1H") == 0) H1FREQMHZ = bf2;
else if (strcmp(nuc3, "1H") == 0) H1FREQMHZ = bf3;
else if (strcmp(nuc4, "1H") == 0) H1FREQMHZ = bf4;
else if (strcmp(nuc5, "1H") == 0) H1FREQMHZ = bf5;
else if (strcmp(nuc6, "1H") == 0) H1FREQMHZ = bf6;
else {
  sprintf(text, "Cannot find a channel defined as 1H nucleus, aborting.");
  STOPMSG(text)
}
/* Could request H1FREQMHZ from user in this case but it might disturb the queueing capability */

/* comment out if (dim == 1) bracket, always calibrate direct dimension */
/* if (dim == 1) { */
	FETCHPAR("NUC1", nuc1);
	FETCHPAR("BF1", &bf1);
	/*FETCHPAR("SR", &sr1);*/
	FETCHPAR("SWH", &swh1);
	FETCHPAR("SI", &si1);
	hzppt1=swh1/si1;

	pixelcorrection=0; /* note: if you are in Topspin versions below 4, you can instead define it as follows so spectra with different SI overlay
	pixelcorrection=hzppt1*0.5 */
	/* here is how i tried to do it before discovering that there is no consistent name for the command to determine the topspin version 
if (TSversion < TSversionWhereHalfPointReferencingErrorWasFixed) pixelcorrection=hzppt1*0.5; */
if (strcmp(nuc1, "1H") == 0) sr1 = sr;
else if (strcmp(nuc1, "13C") == 0) sr1 = ((H1FREQMHZ+sr*1e-6)*CRATIO-bf1)*1e6;
else if (strcmp(nuc1, "15N") == 0) sr1 = ((H1FREQMHZ+sr*1e-6)*NRATIO-bf1)*1e6;
else if (strcmp(nuc1, "31P") == 0) sr1 = ((H1FREQMHZ+sr*1e-6)*PRATIO-bf1)*1e6;
else if (strcmp(nuc1, "19F") == 0) sr1 = ((H1FREQMHZ+sr*1e-6)*FRATIO-bf1)*1e6;
else if (strcmp(nuc1, "2H") == 0) sr1 = ((H1FREQMHZ+sr*1e-6)*H2RATIO-bf1)*1e6;
else {
  sprintf(text, "No reference available for dim1 nuc1 %s", nuc1);
  STOPMSG(text)
}

STOREPAR("SR", sr1-pixelcorrection)
#ifdef DEBUG
Proc_err(DEF_ERR_OPT,"new SR(F1)  %f \n old SR(F1)  %f\n", sr,sr1);
#endif
/*
            }
*/

/*Proc_err(DEF_ERR_OPT,"1D data set\nNo referencing needed");*/

if ((dim < 1) || (dim > 3)) {
  Proc_err(DEF_ERR_OPT, "Dimensionality of data set is out of scripts range: allowed dim 1,2, or 3!");
}

	
if ((dim ==2) || (dim ==3)) {
	FETCHPAR1("NUC1", nuc2);
	FETCHPAR1("BF1", &bf2);
	FETCHPAR1("SWH", &swh2);
  FETCHPAR1("SI", &si2);
  hzppt2=swh2/si2;
  pixelcorrection=0; /* for TS below 4 you can instead define pixelcorrection=hzppt2*0.5; */
  /* removed next line because it is not recognized in TS3.2 */
/* if (TSversion < TSversionWhereHalfPointReferencingErrorWasFixed) pixelcorrection=hzppt2*0.5; */
/* indirect dimension 1 */
  if (strcmp(nuc2, "1H") == 0) sr2 = sr;
  else if (strcmp(nuc2, "13C") == 0) sr2 = ((H1FREQMHZ+sr*1e-6)*CRATIO-bf2)*1e6;
  else if (strcmp(nuc2, "15N") == 0) sr2 = ((H1FREQMHZ+sr*1e-6)*NRATIO-bf2)*1e6;
  else if (strcmp(nuc2, "31P") == 0) sr2 = ((H1FREQMHZ+sr*1e-6)*PRATIO-bf2)*1e6;
  else if (strcmp(nuc2, "19F") == 0) sr2 = ((H1FREQMHZ+sr*1e-6)*FRATIO-bf2)*1e6;
  else if (strcmp(nuc2, "2H") == 0) sr2 = ((H1FREQMHZ+sr*1e-6)*H2RATIO-bf2)*1e6;
  else {
    sprintf(text, "No reference available for 2D or 3D in dim2 nuc1 %s", nuc2);
    STOPMSG(text)
  }
  STOREPAR1("SR", sr2-pixelcorrection)
  #ifdef DEBUG
  Proc_err(DEF_ERR_OPT,"new SR(F2)  %f \n old SR(F2)  %f\n", sr,sr2);
  #endif
}

if (dim==3) {
	FETCHPAR3("NUC1", nuc3);
	FETCHPAR3("BF1", &bf3);
	/*FETCHPAR3("SR", &sr3);*/
	FETCHPAR3("SWH", &swh3);
  FETCHPAR3("SI", &si3);
  hzppt3=swh3/si3;
  pixelcorrection=0; /* for TS below 4 you can instead define pixelcorrection=hzppt3*0.5; */
  /* removed next line because it is not recognized in TS3.2 */
/* if (TSversion < TSversionWhereHalfPointReferencingErrorWasFixed) pixelcorrection=hzppt3*0.5; */
	if (strcmp(nuc3, "1H") == 0) sr3 = sr;
        else if (strcmp(nuc3, "13C") == 0) sr3 = ((H1FREQMHZ+sr*1e-6)*CRATIO-bf3)*1e6;
        else if (strcmp(nuc3, "15N") == 0) sr3 = ((H1FREQMHZ+sr*1e-6)*NRATIO-bf3)*1e6;
        else if (strcmp(nuc3, "31P") == 0) sr3 = ((H1FREQMHZ+sr*1e-6)*PRATIO-bf3)*1e6;
        else if (strcmp(nuc3, "19F") == 0) sr3 = ((H1FREQMHZ+sr*1e-6)*FRATIO-bf3)*1e6;
        else if (strcmp(nuc3, "2H") == 0) sr3 = ((H1FREQMHZ+sr*1e-6)*H2RATIO-bf3)*1e6;
	else {
		 sprintf(text, "No reference available for 2D or 3D in dim3 nuc3 %s", nuc3);
  		 STOPMSG(text)
	}
  #ifdef DEBUG
  Proc_err(DEF_ERR_OPT,"bf1+sr1= %f",bf1+sr1*1e-6);
  if (strcmp(nuc3, "13C") == 0) bfzero = (H1FREQMHZ+sr*1e-6)*CRATIO;
  if (strcmp(nuc3, "15N") == 0) bfzero = (H1FREQMHZ+sr*1e-6)*NRATIO;
  if (strcmp(nuc3, "31P") == 0) bfzero = (H1FREQMHZ+sr*1e-6)*PRATIO;
  if (strcmp(nuc3, "19F") == 0) bfzero = (H1FREQMHZ+sr*1e-6)*FRATIO;
  if (strcmp(nuc3, "2H") == 0) bfzero = (H1FREQMHZ+sr*1e-6)*H2RATIO;
  srnew=(bfzero-bf3)*1e6;
  Proc_err(DEF_ERR_OPT,"(bf1+sr1)*XRATIO= %f\nsr(F1)=bfzero-bf3= %f",bfzero,srnew);
  #endif
  	STOREPAR3("SR", sr3-pixelcorrection)
  #ifdef DEBUG
  Proc_err(DEF_ERR_OPT,"new SR(F1)  %f \n old SR(F1)  %f\n", sr,sr3);
  #endif

}

#ifdef DEBUG
Proc_err(DEF_ERR_OPT,"NUC1(aq) %s \nNUC1(F2) %s\nNUC1(F1) %s", nuc1,nuc2,nuc3);
Proc_err(DEF_ERR_OPT,"BF1(aq) %f \nBF1(F2) %f\nBF1(F1) %f", bf1,bf2,bf3);
Proc_err(DEF_ERR_OPT,"SR(aq) %f \nSR(F2) %f\nSR(F1) %f", sr1,sr2,sr3);
Proc_err(DEF_ERR_OPT,"SR(aq) %f \nSR(F2) %f\nSR(F1) %f", sr1,sr2,sr3);
Proc_err(DEF_ERR_OPT,"SW(aq) %f \nSW(F2) %f\nSW(F1) %f", swh1,swh2,swh3);
Proc_err(DEF_ERR_OPT,"SI(aq) %i \nSI(F2) %i\nSI(F1) %i", si1,si2,si3);
Proc_err(DEF_ERR_OPT,"HZpPT(aq) %f \nHZpPT(F2) %f\nHZpPT(F1) %f", hzppt1,hzppt2,hzppt3);
#endif



#ifdef DEBUG
Proc_err(DEF_ERR_OPT,"bf1+sr1= %f",bf1+sr1*1e-6);
if (strcmp(nuc2, "13C") == 0) bfzero = (H1FREQMHZ+sr*1e-6)*CRATIO;
if (strcmp(nuc2, "15N") == 0) bfzero = (H1FREQMHZ+sr*1e-6)*NRATIO;
if (strcmp(nuc2, "31P") == 0) bfzero = (H1FREQMHZ+sr*1e-6)*PRATIO;
if (strcmp(nuc2, "19F") == 0) bfzero = (H1FREQMHZ+sr*1e-6)*FRATIO;
if (strcmp(nuc2, "2H") == 0) bfzero = (H1FREQMHZ+sr*1e-6)*H2RATIO;
srnew=(bfzero-bf2)*1e6;
Proc_err(DEF_ERR_OPT,"(bf1+sr1)*XRATIO= %f\nsr(F2)=bfzero-bf2= %f",bfzero,srnew);
#endif
      
VIEWDATA_SAMEWIN

sprintf(message,"%s\n13C  %f   %s\n15N  %f   %s\n31P  %f   %s\n19F  %f   %s\n2H  %f %s\n\n%s\n\n%s\n\n%s",MSGSTART,CRATIO,CREF,NRATIO,NREF,PRATIO,PREF,FRATIO,FREF,H2RATIO,H2REF,REFNOTE,REFNOTE2,REFNOTE3);

if (quiet==1) STOP;
if (quiet==0) QUITMSG(message);