/******* -*-C++-*- **********************************************/
/*  au_dosy_prep    26.04.2012                                  */
/****************************************************************/
/*  Short Description :                                         */
/*  AU program for acquisition of DOSY data in automation.      */
/****************************************************************/
/*  Keywords :                                                  */
/*  DOSY, automation, acquisition                               */
/****************************************************************/
/*  Description/Usage :                                         */
/*  AU program for acquisition of DOSY data in automation.      */
/*  The AU tries to find the right diffusion parameters.        */
/*  rga is executed, then acquisition is started.               */
/****************************************************************/
/*  Author(s) :                                                 */
/*  Name            : Rainer Kerssebaum                         */
/*  Organisation    : Bruker BioSpin GmbH                       */
/*  Email           : rainer.kerssebaum@bruker.com		*/
/****************************************************************/
/*  Name    Date    Modification:                               */
/*  rke     111117  created                                     */
/*  rke     170130  limit D20 to 2 sec                          */
/****************************************************************/
/*
$Id:$
*/

AUERR = local_au(curdat);
QUIT


static double round10(double);
static int getgradcalib(double*);
static int getshapescale(const char*, double*);
static int phase1d(const char*);
static int read1ddatapoints(const char*, const char*, const int, const int, int*);

#define TLL 1024
#define MAXPEAKS 2
static const int debugflag = 0;
static const int defaultflag = 0;
static const int proconly = 0;
static const int maxpeaks = MAXPEAKS;


int local_au(const char* curdat)
{
int     i, numPeaks, peakindex, pmod, s_expno, s_procno, td1, ymax;
int     datapointval[MAXPEAKS + 1][2 + 1], maxi[MAXPEAKS + 1];
float   offset, ph0, ph1, DELTAbig = 0.06, deltalittle = 0.0012;
double  corr=1., corr2, deltalittlenew, DELTAbignew, gcalib;
double  maxfp, peakintensity, shapescale, sw=0., D[MAXPEAKS + 1], decayfact[MAXPEAKS + 1], maxpi[MAXPEAKS + 1];
char    nuc[10], solvent[PATH_MAX], text[TLL];
static double   decaytarget = 0.025, gamma = 4257.7 * 2. * 3.14159, gradlimmin = 600., gradlimmax = 3000.;
static double   gstart = 2., gstop = 98.;
static const char   diffshape[] = "Difftrap";


FETCHPAR("PARMODE",&pmod)
if (pmod != 1)
    STOPMSG("Not a 2D dataset !\naborting ...")

FETCHPAR("NUC1", nuc)
if (strcmp("1H", nuc) )
    STOPMSG("only 1H supported !\naborting ...")

s_expno = expno;
s_procno = procno;

if (proconly == 0)
    {
    WRA(99997)
    ERRORABORT

    DATASET(name, 99997, procno, disk, user)
    }

if (defaultflag == 1)
    {
    STOREPAR("P 30", deltalittle * 500000.0)
    STOREPAR("D 20", DELTAbig)
    STOREPAR("D 21", 0.005)
    STOREPAR("GPNAM6", diffshape)
    STOREPAR("GPNAM7", "SMSQ10.50")
    STOREPAR("GPNAM8", "SMSQ10.50")
    STOREPAR("GPZ 7", -17.13)
    STOREPAR("GPZ 8", -13.17)
    STOREPAR("DS", 4)
    STOREPAR("NS", 16)
    }
FETCHPAR("P 30", &deltalittle)
deltalittle /= 500000.0;
FETCHPAR("D 20", &DELTAbig)
if (debugflag > 0)
    Proc_err(ERROPT_AK_OK, "DEBUG\n BIG_DELTA = %.1f ms\n little_delta = %.1f us", DELTAbig * 1000, deltalittle * 1e6);

STOREPAR("PSCAL", 4)
sprintf(text, "dosy %.1f %.1f 2 l", gstart, gstop);
XCMD(text)

if (proconly == 0)
    {
    RGA
    ZG
    }

STOREPAR("CY", 1.0)
STOREPAR1("SI", 4)
XF2
ERRORABORT
RSR(1, 111)
ERRORABORT
RPROCNO(111)
phase1d(curdat);

FETCHPARS("PHC0", &ph0)
FETCHPARS("PHC1", &ph1)
RPROCNO(s_procno)
STOREPAR("PHC0", ph0)
STOREPAR("PHC1", ph1)
STOREPAR("PH_mod", 1)

XF2
ERRORABORT
/* STOREPAR("SR", 0.0) */
SREF
RSR(1, 111)
RPROCNO(111)

/* Set some parameters for peak picking */
FETCHPARS("OFFSET", &offset)
FETCHPARS("SW", &sw)
FETCHPARS("YMAX_p", &ymax)
FETCHPARS("SOLVENT", solvent)
STOREPAR("MI", 0.02)
STOREPAR("MAXI", 1.1)
STOREPAR("F1P", offset)
STOREPAR("F2P", offset - sw)
sprintf(text, "1H.%s", solvent);
STOREPAR("SREGLST", text)

PP
ERRORABORT

/* open peak picking result file */
numPeaks = readPeakList(PROCPATH(0));

if (numPeaks <= 0)
    {
    Proc_err(DEF_ERR_OPT, "No peak found !\nExecution aborted !");
    ABORT
    }

for (i = 1; i <= maxpeaks; i++)
    maxpi[i] = 0.0;

for (i = 0; i < numPeaks; i++)
    {
    peakintensity = getPeakIntensity(i);
    if (peakintensity > maxpi[1])
        {
        maxpi[1] = peakintensity;
        maxi[1] = getPeakAddress(i);
        maxfp = getPeakFreqPPM(i);
        }
    else if ((peakintensity > maxpi[2]) && (numPeaks > 1))
        {
        maxpi[2] = peakintensity;
        maxi[2] = getPeakAddress(i);
        }
    }
freePeakList();


if (debugflag > 1)
    Proc_err(ERROPT_AK_OK, "DEBUG\nnumPeaks %i\nrel.Int: %.2f %.2f\nPPM: %.2f\nIndex: %i %i\nabs.Int: %i %i",
            numPeaks,maxpi[1],maxpi[2],maxfp,maxi[1],maxi[2],(int)(ymax * maxpi[1]),(int)(ymax * maxpi[2]));

/*  extract datapoints and calculate decay ratio */
for (i = 1; i <= maxpeaks; i++)
    {
    RPROCNO(s_procno)
    RSC(maxi[i], 112)
    RPROCNO(112)
    read1ddatapoints(curdat, PROCPATH("1r"), 2, 0, datapointval[i]);

    if (datapointval[i][2] < 0) datapointval[i][2] = 0;

    decayfact[i] = (double)datapointval[i][2] / (double)datapointval[i][1];

    if (debugflag > 0)
        Proc_err(ERROPT_AK_OK, "DEBUG\n peak %i\n value 1: %i\n value 2: %i\n rel: %.4f",
            i, datapointval[i][1], datapointval[i][2], decayfact[i]);
    }

if ((decayfact[1] == 0) && (decayfact[2] == 0))
    STOPMSG("no decay detectable")

/*  use second biggest peak */
if (decayfact[1] > decayfact[2])    peakindex = 1;
else    peakindex = 2;

/* get gradient shape scaling factor */
if (getshapescale(diffshape, &shapescale) < 0)  return -1;

/*  get gradient calibration */
if (getgradcalib(&gcalib) < 0)  return -1;

/*  estimate D value */
for (i = 1; i <= maxpeaks; i++)
    D[i] = log(1 / decayfact[i]) / (pow((gstop - gstart) * shapescale * gcalib * gamma * deltalittle, 2)
            * (DELTAbig - deltalittle / 3));

if (debugflag > 1)
    {
    Proc_err(ERROPT_AK_OK, "DEBUG\n gradient strength\n value 1: %.2f\n value 2: %.2f",
        gstart * shapescale * gcalib / 100, gstop * shapescale * gcalib / 100);
    Proc_err(ERROPT_AK_OK, "DEBUG\n D\n value 1: %.5g\n value 2: %.5g\n BIG_DELTA %.1f ms\n little_delta %.1f ms",
        D[1], D[2], DELTAbig * 1000, deltalittle * 1000);
    }

/*  get correction if decay too weak or too strong */
if ( (decayfact[peakindex] < (decaytarget - 0.015) ) ||  (decayfact[peakindex] > (decaytarget + 0.015) ) )
    {
    corr = 1 - ( log(decaytarget / decayfact[peakindex]) ) / ( D[peakindex] * pow( gamma * deltalittle * gcalib * shapescale
            * gstop, 2) * (DELTAbig - deltalittle / 3) );
    corr = sqrt(corr);

    if (debugflag > 0)
        Proc_err(ERROPT_AK_OK, "DEBUG %s/%d\n peak %d\n Corr = %.2f\n Corr^2 = %.2f", name, expno, peakindex, corr, corr*corr);

    DATASET(name, s_expno, s_procno, disk, user)

    /*  correct little delta */
    DELTAbignew = DELTAbig * 1000.;
    deltalittlenew = round10(deltalittle * 5e5 * corr);
    if (debugflag > 0)
        Proc_err(ERROPT_AK_OK, "DEBUG\n P30\n old = %.1f us\n new = %.1f us\n actual = %.1f us",
                                deltalittle * 5e5, deltalittlenew, deltalittle * 5e5 * corr);

    /* if littledelta gets too small or too big correct BIGdelta */
    if ( (deltalittlenew < gradlimmin) || (deltalittlenew > gradlimmax) )
        {
        if (deltalittlenew < gradlimmin)
            {
            corr2 = deltalittlenew / gradlimmin;
            deltalittlenew = gradlimmin;
            }
        if (deltalittlenew > gradlimmax)
            {
            corr2 = deltalittlenew / gradlimmax;
            deltalittlenew = gradlimmax;
            }
        DELTAbignew = round10(DELTAbig * 1000 * corr2 * corr2);
        if (DELTAbignew > 2000.0)
        {
        DELTAbignew = 2000.0;
        Proc_err(DEF_ERR_OPT, "D20 too large, max value 2s !");
        }

        if (debugflag > 0)
            Proc_err(ERROPT_AK_OK, "DEBUG\n BIGdelta = %.1f ms\n P30 = %.1f us\n corr2 = %.2f\n corr = %.2f",
                                    DELTAbignew, deltalittlenew, corr2*corr2, corr);
        }

    if (proconly == 0)
        {
        STOREPAR("P 30", deltalittlenew)
        STOREPAR("D 20", DELTAbignew / 1000.0)
        Show_meta(SM_RAWP);
        }
    }

DATASET(name, s_expno, s_procno, disk, user)

FETCHPAR1("TD",&td1);
sprintf(text,"dosy %.1f %.1f %d l y y y", gstart, gstop, td1);
if (proconly == 0)
    {
    XCMD(text)
    }

if (debugflag > 0)
    Proc_err(ERROPT_AK_OK, "DEBUG\n %s\n P30 = %.1f us\n D20 = %.1f", text, deltalittlenew, DELTAbignew);

return 0;
}

/* subroutines ****************************************************************/
static double round10(double value)
{
return floor(value / 10. + .5) * 10.;
} /* end subroutine */


static int phase1d(const char* curdat)
{
int     digmod;

FETCHPARS("DIGMOD", &digmod)
if (digmod == 3)
    APK0
else
    APK
ERRORABORT

return 0;
} /* end subroutine */


static int getshapescale(const char* shapename, double* scal)
{
FILE    *fptr;
char    inputline[TLL], path[PATH_MAX];

if (getParfileDirForRead(shapename, GP_DIRS, path) < 0)
    {
    Proc_err(DEF_ERR_OPT, "%s: %s", shapename, path);
    return -1;
    }

if ((fptr = fopen(path, "r")) == 0)
    {
    Proc_err(DEF_ERR_OPT, "Cannot open shape file for reading\n%s", path);
    return -1;
    }

while (fgets(inputline, sizeof(inputline), fptr) != 0)
    {
    if (strncmp("##$SHAPE_INTEGFAC=", inputline, 18) == 0)
        {
        sscanf(inputline, "%*s%lf", scal);
        }
    }
fclose(fptr);

if (debugflag > 1)
    Proc_err(ERROPT_AK_OK,"DEBUG_sub\nshapename = %s\nshape integfactor = %.2f", shapename, *scal);

return 0;
} /* end subroutine */


static int read1ddatapoints(const char* curdat, const char* path, const int numpoints, const int points2skip, int* value)
{
int     i, endian, *fin, si, sizeofint;
FILE    *fpin;

FETCHPARS("BYTORDP", &endian)
FETCHPARS("SI", &si)
sizeofint = sizeof(int);
fin = calloc(si, sizeofint);

if ((fpin = fopen(path, "rb")) == 0)
    {
    Proc_err(DEF_ERR_OPT, "Cannot open file for reading\n%s", path);
    ABORT
    }

if (fread(fin, sizeofint, si, fpin) != si)
    STOPMSG("Can not read 1r file")
fclose(fpin);
local_swap4(fin, sizeofint * si, endian);

for (i = 1; i <= numpoints; i++)
    value[i] = fin[i - 1 + points2skip];

free(fin);

return 0;
} /* end subroutine */


static int getgradcalib(double* scal)
{
double  gcalib=53.5;
char    unit[TLL], path[PATH_MAX];
FILE    *fptr;

sprintf(path, "%s/conf/instr/gradient_calib", PathXWinNMRInst());
if ((fptr = fopen(path, "r")) != 0 )
    {
    if (fscanf(fptr, "%lf %s", &gcalib, unit) == 2)
        {
        if (strcmp(unit, "G/cm") != 0)
            {
            if (strcmp(unit, "G/mm") == 0)
                gcalib *= 10;
            else
                gcalib = 53.5;
            }
        }
    else
        {
        gcalib=53.5;
        }
    fclose(fptr);
    }

*scal = gcalib;

if (debugflag > 1)
    Proc_err(ERROPT_AK_OK,"DEBUG_sub\n gradient calibration\n %.2f", *scal);

return 0;
} /* end subroutine */
