/****************************************************************/
/* qcpmgaddjt 01.02.2003                      */
/****************************************************************/
/* Short Description : */
/* Program to store echos from qcpmgall.av experiment */
/* by adding them up */
/****************************************************************/
/* Keywords : */
/* quadrupolar echo, CPMG */
/****************************************************************/
/* Description/Usage: */
/* After data acquisition using qcpmgall.av */
/* if digmod is used please convert data to analog with convdta command */
/* WARNING : if program FAILS data will be CORRUPTED */
/* ALWAYS use this AU program on a COPY of dataset */
/****************************************************************/
/* TESTED ON avance II with DRU ONLY */
/* Program takes data points with echo information */
/* and stores as result summation of all echos. */
/* Each echo is apodised by a gaussian function using parameter LB */
/* Program stores result as time domain data into 1r */
/* and 1i files for further processing, */
/* raw data are restored after program completes ! */
/****************************************************************/
/* Author(s) : */
/* Name : Stefan Steuernagel */
/* Organisation : Bruker Analytik */
/* Email : stefan.steuernagel@bruker.de */
/* Name : Julien TREBOSC */
/* Organisation : UCCS, University of Lille 1 */
/* Email : julien.trebosc@univ-lille1.fr */
/****************************************************************/
/* Name Date Modification: */
/* ste 030201 created */
/* jt  070501 add gaussian apodisation of individual echos */
/* jt  070501 add first 1/2 echo */
/* jt  070501 add noise normalisation */
/****************************************************************/



int size, sizeout, si2, tdout, td2, ireal, sizeofint, parmod, zero=0;
char binfile[255], pulprog[255];
int wdw2,pk2,bc2,ft2,mem2;
int *in1, *in2;
int order_a_BIT, first, echopoints, echos, l22, dead, pivot, zeroth;
float del3, del6, pul3, pul4, dwell, glb;
double swh;
char s1[255];

GETCURDATA;

/* ---------------------------------------------*/
/* make sure program works on proper data set */

FETCHPARS("PARMODE",&parmod);
if (parmod != 0) STOPMSG("Program only works on 1D data");
FETCHPARS("PULPROG",pulprog);
//if (strcmp(pulprog,"qcpmgall.av")!=0) { Proc_err(0,"Wrong AU program for %s pulse program",pulprog);
//return(AUERR);}

/* ---------------------------------------------*/
/* make sure that pivot point can be entered */

(void) sprintf(s1,"%s %s","au_qcpmgadd AU program","- Note \n\
- If you don't know the pivot point yet,\n\
  press Cancel and check with\n\
  'calibrate'\n\
- otherwise press OK and continue");
pivot=0;
//FETCHPARS("FL1",&pivot);
if (pivot>80) {
 AUERR=Proc_err(ERROPT_AK_CAN|ERROPT_BEEP_DEF,"%s",s1);
 if (AUERR == ERR_CANCEL) ABORT
}

/* ---------------------------------------------*/
/* get relevant acquisition parameters */

FETCHPARS("TD",&td2);
FETCHPARS("SI",&si2);
FETCHPARS("SW_h",&swh);
FETCHPARS("D 6",&del6);
FETCHPARS("D 3",&del3);
FETCHPARS("P 3",&pul3);
FETCHPARS("P 4",&pul4);
FETCHPARS("L 22",&l22);
FETCHPARS("DW",&dwell);
FETCHPARS("BYTORDA",&order_a_BIT);

/* ---------------------------------------------*/
/* check for starting point */

if (pivot>80) pivot=(int)((10/dwell));
if (strlen(cmd) == 0) 
{
    GETINT("Enter pivot point : ",pivot);
}
else
{
    if (1 != sscanf(cmd, "%d", &pivot))
    {
STOPMSG("illegal input");
    }
}

/* ---------------------------------------------*/
/* calculate required loop parameters */

zeroth=(int)((pul3+del3*2000000.0+pul4)/dwell+0.5)+pivot;
zeroth=(zeroth/2)*2;
first=(int)((del6*1000000.0+del3*2000000.0+pul4)/dwell+0.5)+zeroth;
first=(first/2)*2;
echopoints=(int)(del6*2000000.0/dwell+0.5);
echopoints=(echopoints/2)*2;
dead=(int)((del3*2000000.0+pul4)/dwell+0.5)+2;
dead=(dead/2)*2;
echos=l22;
tdout=echopoints;
sizeout=tdout*sizeof(int);
    GETINT("Enter number of echoes to add : ",echos);
if (td2 < ((((del6*2e6+del3*2e6+pul4)*echos)+(del6+del3)*2e6)/dwell)) {
Proc_err(0,"Inconsistent FID file : too many echos for TD ");
exit(1);
}
// ask for apodization of individual echo
FETCHPAR("LB",&glb);
// GETFLOAT("Enter apodization : ",glb);

//Proc_err(INFO_OPT,"zero=%d\n first=%d\n",zeroth,first);
/* ---------------------------------------------*/
/* open file and allocate array  */

(void)sprintf(binfile,"%s/data/%s/nmr/%s/%d/fid",disk,user,name,expno);
sizeofint=sizeof(int);
size=td2*sizeofint;
in1=calloc(td2,sizeof(int));
in2=calloc(td2,sizeof(int));

if ( (ireal=open(binfile, O_RDWR)) == -1) {
Perror(DEF_ERR_OPT,binfile);
ABORT;
}
if ( read(ireal,in1,size) == -1 ) {
Perror(DEF_ERR_OPT,binfile);
ABORT;
}

/* ---------------------------------------------*/
/* correct for byteorder   */

local_swap4(in1,size,order_a_BIT);

/* ---------------------------------------------*/
/* create final output array   */

i1=first;
i2=0;

TIMES(echopoints);
in2[i2]=0;
i2++;
END;

i2=0;

// calculate lbfactor
int centerpoint=(echopoints/4)*2;
float lbfactor=1.8867*1e-6*dwell*2.0;
float apod=1.0;
apod=exp(-pow(abs(i2-centerpoint)*glb*lbfactor,2));

// printf("dwell=%f\n apod max=%f\n",dwell,apod);

// sum echo and apply apodisation

TIMES(echopoints);
TIMES2(echos);
if (i2%2) { apod=exp(-pow(abs(i2-centerpoint)*glb*lbfactor,2)); }
in2[i2]+=(int) (apod* (float) (in1[i1]));
i1+=(dead+echopoints);
END;
i2++;
i1=first+i2;
END;

int yes=1;
GETINT("add first 1/2 echo ? 0=no 1=yes",yes);
if(yes) {
i2=centerpoint;
i1=zeroth;
TIMES2(centerpoint);
if (i2%2) { apod=exp(-pow(abs(i2-centerpoint)*glb*lbfactor,2)); }
in2[i2]+=(int) (apod* (float) (in1[i1]));
i1++;
i2++;
END;
}

int nyes=1;
GETINT("normalise noise ? 0=no 1=yes",nyes);
if (nyes) {
 float nfactor;
 nfactor=1/sqrt(echos*2+yes);
 i2=0;
 TIMES(echopoints)
  in2[i2]=(int) (nfactor* (float) (in2[i2]));
  i2++;
 END
}

if ( lseek(ireal, 0, SEEK_SET) == -1 ) {
Perror(DEF_ERR_OPT,binfile);
ABORT;
}

/* ---------------------------------------------*/
/* correct for byteorder   */

local_swap4(in2,sizeout,order_a_BIT);

if ( write(ireal,in2,sizeout) == -1 ) {
Perror(DEF_ERR_OPT,binfile);
ABORT;
}

/* ---------------------------------------------*/
/* correct relevant acquisition parameters */

STOREPARS("TD",tdout);
STOREPAR("TDeff",tdout);
STOREPAR("SI",tdout*2);

/* ---------------------------------------------*/
/* prepare for processing */

FETCHPAR("WDW",&wdw2)
FETCHPAR("PH_mod",&pk2)
FETCHPAR("BC_mod",&bc2)
FETCHPAR("FT_mod",&ft2)
FETCHPAR("ME_mod",&mem2)

/* ---------------------------------------------*/
/* set required parameters to NO */

STOREPAR("WDW",0)
STOREPAR("PH_mod",0)
STOREPAR("BC_mod",0)
STOREPAR("FT_mod",0)
STOREPAR("ME_mod",0)

/* ---------------------------------------------*/
/* store new FID into 1r and 1i */

TRF;

/* ---------------------------------------------*/
/* restore original FID */

if ( lseek(ireal, 0, SEEK_SET) == -1 ) {
Perror(DEF_ERR_OPT,binfile);
ABORT;
}

/* ---------------------------------------------*/
/* correct for byteorder   */

local_swap4(in1,size,order_a_BIT);

if ( write(ireal,in1,size) == -1 ) {
Perror(DEF_ERR_OPT,binfile);
ABORT;
}

/* ---------------------------------------------*/
/* restore original parameters */

STOREPAR("WDW",wdw2)
STOREPAR("PH_mod",pk2)
STOREPAR("BC_mod",bc2)
STOREPAR("FT_mod",ft2)
STOREPAR("ME_mod",mem2)
STOREPARS("TD",td2);
STOREPAR("TD",td2);
STOREPAR("TDeff",td2);
//STOREPARS("FL1",pivot);

free(in1);
free(in2);
(void)close(ireal);
(void)sprintf(text,"QCPMG echo summation done");
Show_status(text);
QUIT