/*** ^^A -*-C++-*- **********************************************/
/*	zg_dfs			17.08.2000                      */
/****************************************************************/
/*	Short Description :					*/
/*	Program to calculate shape file for double frequency	*/
/*	sweep and subsequent data-acquisition			*/
/*	The first and the last 10% of the sweep are multiplied	*/
/*	by a sine or cosine function, respectively.		*/
/****************************************************************/
/*	Keywords :						*/
/*	adiabatic sweep, shaped pulse, MQMAS			*/
/****************************************************************/
/*	Description/Usage:					*/
/*	This program runs with xaua from within poptau or	*/
/*	paropt etc.						*/
/*	This means it can be used to optimise the parameters	*/
/*	for the sweep in an automatic way.			*/
/*	Do ased first to make sure that the parameters listed	*/
/*	below are set appropriately:				*/
/*								*/
/*	p2 : = 1s / (cnst31 * cnst0), duration of the sweep	*/
/*	where cnst31 = spinning frequency			*/
/*	cnst0 > 1 : fraction of a rotor period			*/
/*	cnst0 < 1 : multiple of a rotor period			*/
/*	sp1 : power level of the sweep				*/
/*	spnam1 : dfs is always used as default name		*/
/*	cnst1 : (in kHz) Startfrequency of the sweep		*/
/*	cnst2 : (in kHz) Endfrequency of the sweep		*/
/*	cnst3 : (in ns) used to define the timing resolution 	*/
/*			of the sweep				*/
/*	Make sure that a customer-made pulse program uses	*/
/*	these parameters for the same purpose!			*/
/*	This AU program is suitable to be used with the 	*/ 
/*	following library pulse programs:			*/
/*	dfs90sel :	sweep followed by 90 degree sel. pulse	*/
/*	mp3qdfs.av or mp3qdfs :					*/
/*			phase-modulated 3-pulse			*/
/*			shifted echo split-t1 3QMAS		*/
/*	mp3qdfsz.av or mp3qdfsz :				*/
/*			amplitude-modulated 4-pulse z-filtered	*/
/*			split-t1 3QMAS for spin-3/2		*/
/****************************************************************/
/*	Author(s) :						*/
/*	Name		: Stefan Steuernagel			*/
/*	Organisation	: Bruker BioSpin GmbH			*/
/*	Email		: stefan.steuernagel@bruker.de		*/
/****************************************************************/
/*	Name		Date	Modification:			*/
/*	ste		000817	created				*/
/*	ste		031128	changed to allow sweep of 	*/
/*				fractions of rotor periods	*/
/*				by setting l0 = 1, 2, or 4	*/
/*				blank character included for	*/
/*				npoints in parameter part	*/
/*	ste		060504	change l0 to cnst0 to allow	*/
/*				multiples AND fractions of a	*/
/*				rotor period			*/
/****************************************************************/
/*
$Id: zg_dfs,v 1.7 2006/12/06 13:14:02 wem Exp $
*/

AUERR = zg_dfs(curdat);
QUIT

#include <ShapeIO/ShapeIOC.h>

int zg_dfs(const char* curdat)
{
  char   outputfile[PATH_MAX];
  double sweep, Startfreq, Endfreq;
  double* wbuf;
  double* wp;
  float fval;
  int i, npoints;

  /* read parameters for sweep calculation from acqu file	*/
  /* and store filename for shaped pulse 			*/

  /* used to calculate duration of sweep as multiple or		*/
  /* fraction of a rotor period					*/ 
  FETCHPAR("CNST 31",&fval) 
  sweep = 1. / fval;
  FETCHPAR("CNST 0",&fval)
  sweep /= fval;

  /* get start- and end-frequencies */
  FETCHPAR("CNST 1",&fval)
  Startfreq = (double)fval * 1000;
  FETCHPAR("CNST 2",&fval)
  Endfreq = (double)fval * 1000;

  /* used for timing resolution of sweep */
  /* hardware dependent down to 25 ns    */
  FETCHPAR("CNST 3",&fval)

  /* calculate number points in the file */
  npoints = (int)((sweep * 1.e9) / fval);

  /* generate shaped pulse filename	 */
  if (getParfileDirForWrite("dfs", SHAPE_DIRS, outputfile) < 0)
  {
    Proc_err(DEF_ERR_OPT, "dfs: %s", outputfile);
    return -1;
  }

  STOREPAR("SPNAM1","dfs")

  /* allocate memory buffer for amplitude and phase */
  wp = wbuf = (double*)calloc(2 * npoints, sizeof(double));
  if (wbuf == 0)
  {
    Proc_err(DEF_ERR_OPT, "Not enough buffer space");
    return -1;
  }

  sweep *= (2. * M_PI) / npoints;
  Endfreq = (Startfreq - Endfreq) / (2 * (npoints - 1));

  for (i = 0; i < npoints; i++)
  {
    double a;

    if (i < npoints / 10)
      a = sin((10. * M_PI * i) / (2 * npoints));
    else if (i < (9 * npoints) / 10)
      a = 1.0;
    else
      a = cos((10. * M_PI * (i - (9 * npoints) / 10)) / (2 * npoints));

    *wp++ = a * 100. * cos(i * sweep * (Startfreq - Endfreq * i));
  }

  /* store amplitude values in shaped pulse file	*/
  if (writeShapeC(outputfile, wbuf, npoints, 1, 0, 0) == 0)
    return -1;

  free(wbuf);

  /* start acquisition */
  ZG
  ERRORABORT

  return 0;
}