/************************************************************/
/* shim1h 26/11/2010 */
/************************************************************/
/* Description:                                             */
/* AU program to Topshim on the most prominent solvent peak */
/* in a mix of protonated solvents.                         */
/* The program turns the lock off, then runs a scout 1H     */
/* spectrum in expno 99999.  After processing and peak-     */
/* picking, it finds the largest available  well-separated  */
/* peak to shim on and does a 1D 1H Topshim on it.          */
/* If it can't find a well-separated peak, it shims on the  */
/* largest peak and hopes for the best.                     */
/************************************************************/
/* Use:                                                     */
/* shim1h (or XAU shim1h)                                   */
/* can be used as a shim routine in ICON-NMR as XAU shim1h  */
/************************************************************/
/* Name   Date    Modification                              */
/* rss    101124  created                                   */
/* rss    101125  added handling for case where no well-    */
/*                separated peak found                      */
/* rss    101126  error message added if no separated peak  */
/*                found, hwcal.out replaced hwcal.txt       */
/* rss    101202  lowered peak picking threshold            */ 
/* rss    101210  corrected for hwcal failure causing       */
/*                of shim1h                                 */
/************************************************************/

/* declare variables */
int curexpno, numofpeaks, peak, peakmax, i, successflag = 0, compareflag = 0 ;
float ppmvals[128], intensities[128];
double sf, sw_p, f1p, f2p, f1porig, f2porig;
float offset, adr, hz, ppmmax, intensmax, hwmax, leftlimit, rightlimit, selwid;
float ppmtest, intenstest;
char ppfile[PATH_MAX], line[PATH_MAX], hwfile[PATH_MAX], topshimcmd[256];
FILE *f_pp, *f_hw;

GETCURDATA
/* set locknuc to off and turn lock and sweep off */
STOREPAR("LOCNUC","off");
LOCK_AND_SWEEP_OFF;

/* run scout spectrum (pulprog=zg30, ns=1, locnuc=off, */
/* expno=99997, rg=1) */
curexpno = expno;
DATASET(name,99997,procno,disk,user)
RPAR("PROTON","all")
STOREPAR("SOLVENT","no")
STOREPAR("LOCNUC","off")
STOREPAR("DS",0)
STOREPAR("NS",1)
STOREPAR("RG",1.0)

GETPROSOL

ZG

/* process */
STOREPAR("LB",4.0)
EFP
APK
ABS

/* pick peaks */
STOREPAR("MI",10.0);
STOREPAR("PSCAL",0);
STOREPAR("MAXI",1500.0);
STOREPAR("CY",1000.0);
STOREPAR("PSIGN",0);
FETCHPARS("OFFSET",&offset)
FETCHPARS("SW_p",&sw_p)
FETCHPARS("SF",&sf)
f1p = offset;
f2p = f1p - sw_p / sf;
STOREPAR("F1P",f1p)
STOREPAR("F2P",f2p)
PP

/* process peak picking list, looking for largest well separated peak */
XCMD("sendgui convertpeaklist txt")
ERRORABORT
numofpeaks = 0;

(void) strcpy(ppfile,PROCPATH("peak.txt"));

if ((f_pp = fopen(ppfile,"r")) == NULL)
{
    STOPMSG("shim1h failed\n\n peaklist does not exist ")
}
/* Read comment lines and discard them */

TIMES(4)
  fgets (line,sizeof(line),f_pp);
END

/* read all the peak information into variables */
while(( AUERR = fscanf( f_pp,"%d %f %f %f %f",&peak,&adr,&hz,&ppmvals[numofpeaks],&intensities[numofpeaks])) > 0)
{
  numofpeaks = numofpeaks + 1;
} 

if ( numofpeaks == 0 )
{
	STOPMSG("shim1h failed\n\nno peaks found")
} 
fclose(f_pp);

/* starting by setting various flags / values to zero: */
successflag = 0;
compareflag = 0 ;
ppmmax = 0.0;
intensmax = 0.0;
int problemflags[numofpeaks], problemflagsorig[numofpeaks];
for (i = 0; i < numofpeaks; i++)
{
	problemflags[i] = 0;
	problemflagsorig[i] = 0;
}

do {

  /* inner loop: look for greatest possible peak, if it's not been flagged */
	for (i = 0; i < numofpeaks; i++)
	{
		if ( ( intensities[i] > intensmax ) && ( problemflags[i] == 0 ) )
		{
			ppmmax = ppmvals[i];
			intensmax = intensities[i];
			peakmax = i;
		}
	}

	/* if no peak is suitable, set successflag to -1 */
	if ( ppmmax == 0.0 )
	{
		successflag = -1;
	}
	else
	{
		/* once a possible peak is found, measure its half width */
		f1p = ppmmax + 1;
		f2p = ppmmax - 1;
		STOREPAR("F1P",f1p)
		STOREPAR("F2P",f2p)
		XAU("hwcal","")
		(void) strcpy(hwfile,PROCPATH("hwcal.out"));

		if ((f_hw = fopen(hwfile,"r")) == NULL)
		{
			hwmax = 20.0 ; 
		}
		else
		{
			fscanf(f_hw,"%f",&hwmax);
			fclose(f_hw);
		}

		/* then check the peaks to see whether there's a peak nearby */
		leftlimit = ppmmax + 3 * hwmax / sf ;
		rightlimit = ppmmax - 3 * hwmax / sf ;
		for (i = 0; i < numofpeaks; i++)
		{
			if ( ppmvals[i] != ppmmax )
			{
				if ( ( ppmvals[i] < leftlimit ) && ( ppmvals[i] > rightlimit ) )
				{
					problemflags[peakmax] = 1 ;  // potential maximum is too close to another peak
					problemflags[i]  = 1 ;       // the other peak is too close to the current best candidate
				}
			}
		}
 
		/* compare problemflags and problemflagsorig to see whether any problems have occurred */
		for (i = 0; i < numofpeaks ; i ++ )
		{
			if ( problemflags[i] != problemflagsorig[i] )
			{
				compareflag = 1 ;
			}
		}

		/* if there are problems, reset ppmmax and intensmax to 0.0 */
		if ( compareflag == 1 )
		{
			ppmmax = 0.0 ;
			intensmax = 0.0 ;
			compareflag = 0 ;
			for ( i = 0; i < numofpeaks ; i ++ )
			{
				problemflagsorig[i] = problemflags[i] ;
			}
		}
		else /* if no problems have occurred, set successflag to 1 */
		{
			successflag = 1 ;
		}
	}
} while ( successflag == 0 );

/* create Topshim command */
if ( successflag == 1 )
{
	selwid = 2 * hwmax / sf ;
	if ( selwid < 0.1 )
	{
		selwid = 0.1 ;
	}
}
else /* no well-separated peak was found, so just shim on the biggest peak */
{
	intensmax = 0.0 ;
	for (i = 0; i < numofpeaks; i++)
	{
		if ( intensities[i] > intensmax )
		{
			ppmmax = ppmvals[i];
			intensmax = intensities[i];
		}
	}
	selwid = 0.05 ;
	Proc_err(INFO_OPT, "No well-separated peak was found, so TopShim \n"
	"will just shim on the largest peak.");
}

//sprintf(topshimcmd, "topshim 3d 1h lockoff o1p=%f selwid=%f", ppmmax, selwid);
sprintf(topshimcmd, "topshim 1h lockoff o1p=%f selwid=%f", ppmmax, selwid);
//sprintf(topshimcmd, "topshim plotall convcomp 1h lockoff o1p=%f selwid=%f", ppmmax, selwid);

/* start Topshim */
XCMD(topshimcmd)

QUIT