#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define fdiv(i1,i2) ((float) i1) / ((float) 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 the Varian (Pulsetool) you can simply use
   nawk '{ if(NR>1) print $2, $1*10.23, "1.0" }' < Bruker_shape > Varian_shape
*/

#define MAX_ARRAY 16384

main()
{
	float diff, temp, scale, R, phi, theta[MAX_ARRAY];
	float pi, Tp, f1, phs[MAX_ARRAY], amp[MAX_ARRAY],cutoff;
	float phs2[MAX_ARRAY], amp2[MAX_ARRAY];
	float wc, w1, deltawc, deltaw, delta, phslock[MAX_ARRAY];
	float v1c, vc, v1, deltavc, deltav, pw90;
	float factor, amplitude, plength, offset;
	int i, nplock, np, nelem;
        char fname[255];
        FILE *fp;

/* READING INPUT FROM TERMINAL */

	printf("\n  OUTPUT FILE : ");
        	scanf("%s",fname);
        	fp=fopen(fname,"w");
        	if (fp==NULL) printf("Error opening file\n");
        printf("\n  ADIABATIC PULSE PATTERN USING tanh/TAN ");
        printf("\n  GIVE VALUE FOR R : ");
		i = scanf("%f",&R);
        printf("\n  GIVE VALUE FOR OFFSET OF ENTIRE PULSE IN HZ: ");
		i = scanf("%f",&offset);
        printf("\n  GIVE LENGTH OF ENTIRE PULSE IN MS: ");
		i = scanf("%f",&plength);
		plength = plength/1000;
        printf("\n  GIVE # OF STEPS FOR EACH RAMP (MAX = 16384) : ");
		i = scanf("%d",&np);
   	printf("\n  GIVE NUMBER OF POINTS DURING SPIN-LOCK : ");  
		i = scanf("%d",&nplock);

 
/* INITIALISATION */

	pi = 4*atan(1);
        cutoff = 50.0;
	nelem = np*2 + nplock;

    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)]; 
	}




   for (i=0;i<nelem;i++)
  	{
         phs2[i] = phs2[i] + 360*offset*plength*( fdiv(i,nelem) - 0.5 ); 
	}

 
/* OUTPUT TO FILE  */

	fprintf(fp,"RFVERSION_F\n");

    for (i=0;i<nelem;i++) 
	{	fprintf(fp,"%10.8f %11.8f\n",amp2[i],phs2[i]);
	}

        fclose(fp);
}




/* 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

*/