#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define fdiv(i1,i2) ((double) i1) / ((double) i2)
#define htan(x) (exp(x)-exp(-x))/(exp(x)+exp(-x))		
	

/* Written by Robin de Graaf and Frans Mulder
   Bijvoet center for biomolecular research
   University of Utrecht, the Netherlands
   http://www.nmr.chem.uu.nl/bijvoet.html     

   Major modification: new arrays amp2/phs2 run
   really chronlogically and additional offset of
   the entire pulse can be achieved through an
   extra linear phase ramp

   For Varian (Pulsetool) .DEC pattern

   Version 2 (Frans Mulder, Feb 2004) 
    - some cleaning up of unused variables
    - allow for command line arguments instead of queries
    - divide up on-resonance spin-lock in 255 deg elements
      to allow higher definition of ramps, avoiding
      waveform overflow
*/

#define MAX_ARRAY 65536
#define MAX_RAMP 4096

main(argc, argv)
        int argc;
        char *argv[];

{
	double R, theta[MAX_ARRAY];
	double pi, phs[MAX_ARRAY], amp[MAX_ARRAY],cutoff;
	double phs2[MAX_ARRAY], amp2[MAX_ARRAY];
	double plength, offset, tu, rlength;
	int i, nplock, np, nelem, rest, nbig;
        int int_rlength, int_tu, int_plength;

/* READING INPUT FROM COMMAND LINE */

        if (argc != 5) goto usage;
        rlength = atof(argv[1]);
        plength = atof(argv[2]);
        R = atof(argv[3]);
        tu = atof(argv[4]);

        nplock = (int) (1000*plength/tu + 0.5);
        np = (int) (1000*rlength/tu + 0.5); 
 
        /* convert to integer and nanoseconds to check parameter values */
        int_rlength = (int) (rlength*1.0e6 + 0.5);
        int_plength = (int) (plength*1.0e6 + 0.5);
        int_tu = (int) (tu*1.0e3 + 0.5);

        if (int_tu%50 != 0.0) {
                fprintf(stderr,
          " Error: time unit %f us is not a muliplier of 50 ns\n",tu);
                return 1;
        }

        if (int_rlength%int_tu != 0.0) {
                fprintf(stderr,
          " Error: ramp time %4.1f ms is not a muliplier of time unit %4.1f us\n",rlength,tu);
                return 1;
        }

        if (int_plength%int_tu != 0.0) {
                fprintf(stderr,
          " Error: spin lock period %4.1f ms is not a muliplier of time unit %4.1f us\n",plength,tu);
                return 1;
        }

        if (np > MAX_RAMP) {
                fprintf(stderr,
          " Error: too many steps in ramp, must be less than %d \n",MAX_RAMP);
                return 1;
        }

        rlength = rlength/1000.0;
        plength = plength/1000.0;
        tu = tu/1000000.0;
	nelem = np*2 + nplock;

/*
                fprintf(stderr," plength-tu*nplock = %18.16e \n",plength-tu*nplock);
                fprintf(stderr," rlength-tu*np = %18.16f \n",rlength-tu*np);
                fprintf(stderr," length = %18.16f \n",rlength);
                fprintf(stderr," tu = %18.16f \n",tu);
                fprintf(stderr," np = %d \n",np);
                fprintf(stderr," nplock = %1d \n",nplock);
*/

                fprintf(stderr,"\n  ADIABATIC PULSE PATTERN USING tanh/TAN \n");
                fprintf(stderr,"  Number of ramp elements: %d \n",np);
                fprintf(stderr,"  Number of spin lock elements: %d \n",nplock);
                fprintf(stderr,"  Set DRES = 1.0 \n");
                fprintf(stderr,"  Set PW90 (=1/DMF) for pattern to: %8.4f us \n\n",90.0*tu*1.0e6);


        rest = nplock%255;
        nbig = (nplock-rest)/255;

/* INITIALISATION */

	pi = 4*atan(1);
        cutoff = 50.0;

    for (i=0;i<np;i++)
	{
   	theta[i] = fdiv(i,np-1)*atan(cutoff);
	} 

  	amp[0] = 100*tanh(10.0);
	phs[0] = (360*R/np)*(tan(theta[0])/cutoff); 


/* CALCULATION  */

   for (i=1;i<np;i++)
  	{
	amp[i] = 100*(tanh(10-(fdiv(10*i,np-1))));
        phs[i] = phs[i-1] - (360*R/np)*(tan(theta[i])/cutoff);
	}
   for (i=np;i<2*np;i++)
  	{
         amp[i+nplock+1] = amp[(2*np-1)-i];
         phs[i+nplock+1] = phs[(2*np-1)-i]; 
	}
   for (i=np;i<np+nplock+1;i++)
  	{
         amp[i] = 100;
         phs[i] = 0; 
	}

/* CHANGE THE ORDER INTO A NEW ARRAY  */

     for (i=0;i<np;i++) 
  	{
         amp2[i] = amp[(np-1)-i];
         phs2[i] = phs[(np-1)-i]; 
	}
     for (i=np;i<np+nplock+1;i++)
  	{
         amp2[i] = amp[i];
         phs2[i] = phs[i]; 
	}
     for (i=np+nplock+1;i<nelem;i++)
  	{
         amp2[i] = amp[(nelem+np+nplock-i)];
         phs2[i] = -phs[(nelem+np+nplock-i)]; 
	}


/* OUTPUT TO FILE  */

     for (i=0;i<np;i++) 
	{	printf("   1.0\t%6.1f\t%6.1f\n",phs2[i],10.23*amp2[i]);
  	}
     for (i=0;i<nbig;i++)
	{	printf(" 255.0     0.0  1023.0\n");
	}
     if (rest > 0) {
	 	printf(" %3d.0     0.0  1023.0\n",rest);
        }
     for (i=np+nplock;i<nelem;i++)
	{	printf("   1.0\t%6.1f\t%6.1f\n",phs2[i],10.23*amp2[i]);
  	}


        return 0;
usage:
        fprintf(stderr,"\n  ADIABATIC PULSE PATTERN USING tanh/TAN \n");
        fprintf(stderr,"arg 1: ramp time (ms)\n");
        fprintf(stderr,"arg 2: spin lock time (ms)\n");
        fprintf(stderr,"arg 3: R = dW*T [150]\n");
        fprintf(stderr,"arg 4: time step (us) (multiple of 0.05)\n");
        fprintf(stderr,"\n  Result will be redirected to standard output\n");
        return 1;

}




/* BACKGROUND:

offset = time derivative of phase

example A: What is the frequency sweep of an adiabatic pulse?
For an on-resonance AHP with TAN shaped frequency modulation the phase
difference between the first points equals -1878.26 - -1747.35 = - 130.91
degrees = -130.91/360 = -0.3636 cycles.
The time of the first element is given by the length of the pulse / number
of elements. In our case the ramp is 4 ms containing 275 elements, thus
each element is 14.5454 microseconds. The offset = dy/dx = -25,000 cycles/s.

example B: What is the phase ramp necessary for -1000 Hz offset?
This linear ramp needs to satisfy: phase/element = -1000 cycles/s. The number
of elements of the pulse is 3988, and its length 58 ms. The total phase
difference over the entire pulse is -1000*0.058 = -58 cycles = -20,880 degrees
This means -20,880/3988 = -5.236 degrees per element

*/