#!/bin/perl

################################################################################
#                                                                              #
#                   PulseShaper - Guillermo Moyna - USIP, 2001                 #
#                                                                              #
# PulseShaper is a script for the generation of shaped pulses for use with     #
# selective excitation NMR experiments. The program has several pulse profile  #
# functions and the shaped pulse files it creates are for use with the EFT-90, #
# but are generic and can be easily converted to other formats suitable with   #
# other spectrometers. The user can add new functions by copying existing      #
# functions, modifying them, and adding the appropriate lines to the menu.     #
#                                                                              #
# PulseShaper also calculates the excitation profile for the shaped pulses,    #
# using both Fourier analysis as well as Bloch equation analysis. The data     #
# from these calculations, as tables of excitation versus offset, is saved as  #
# plain text, and can be read and plotted with any spreadsheet program.        #
#                                                                              #
# The 'in-line' comments are really helpful to understand how the program      #
# works, and pulseshaper.txt has other details.                                #
#                                                                              #
# Send suggestions and comments to the author at g.moyna@usip.edu              #
#                                                                              #
################################################################################

#
# This is for the Fourier transforms. A real life saver...
#

use Math::FFT;

#
# Set the version number. Depending on the OS, change certain things...
#

$version = "1.0";

if( ( $^O eq "linux" ) || ( $^O eq "irix" ) )
  {
  $clrstring = "clear";
  $rmstring = "rm";
  }
elsif( $^O eq "MSWin32" )
  {
  $clrstring = "cls";
  $rmstring = "del";
  }

#
# Define Fourier series coeficients for certain 'hardwired' pulses. E-BURP1 coeficients...
#

%EBurp1_sin = ( 1, -0.400,
                2, -1.420,
		3,  0.770,
		4,  0.060,
		5,  0.030,
		6, -0.040,
		7, -0.020,
		8,  0.010 );

%EBurp1_cos = ( 1,  0.880,
                2, -1.040,
	   	3, -0.240,
		4,  0.140,
		5,  0.030,
		6,  0.040,
		7, -0.030,
		8,  0.010 );

#
# E-BURP2 coeficients...
#

%EBurp2_sin = ( 1, -0.120,
                2, -1.790,
		3,  0.010,
		4,  0.410,
		5,  0.080,
		6,  0.070,
		7,  0.010,
		8, -0.040,
		9, -0.010,
	       10, -0.000 );

%EBurp2_cos = ( 1,  0.910,
                2,  0.450,
		3, -1.310,
		4, -0.120,
		5,  0.030,
		6,  0.010,
		7,  0.060,
		8,  0.010,
		9, -0.020,
	       10, -0.010 );

#
# SNEEZE coeficients...
#

%Sneeze_sin = ( 1, -0.197,
                2, -1.772,
		3,  0.204,
		4,  0.619,
		5,  0.076,
		6,  0.039,
		7, -0.025,
		8, -0.060,
		9,  0.005,
	       10,  0.017 );

%Sneeze_cos = ( 1,  0.943,
                2,  0.180,
		3, -1.527,
		4,  0.003,
		5,  0.143,
		6,  0.050,
		7,  0.072,
		8, -0.015,
		9, -0.040,
	       10, -0.005 );

#
# G4-CASCADE coeficients...
#

%G4_cascade_sin = ( 1,  2.185,
                    2, -0.003,
		    3,  0.440,
		    4,  0.134,
		    5, -0.339,
		    6,  0.048,
		    7,  0.065,
		    8, -0.046,
		    9,  0.029,
	           10,  0.012,
		   11,  0.009 );

%G4_cascade_cos = ( 1, -0.199,
                    2,  0.588,
		    3, -1.615,
		    4,  0.884,
		    5, -0.034,
		    6,  0.003,
		    7,  0.043,
		    8,  0.014,
		    9, -0.019,
	           10,  0.014,
		   11,  0.002 );

#
# Define the elements of the menu, and put it on the screen...
#

%PulseShapes = ( 1, "SQUARE",
                 2, "TRAPEZOIDAL",
                 3, "TRIANG-RAMP",
                 4, "GAUSSIAN",
                 5, "HALF-GAUSSIAN",
                 6, "HERMITE",
                 7, "SINC",
                 8, "SINE",
		 9, "E-BURP1",
                10, "E-BURP2",
	        11, "SNEEZE",
		12, "G4-CASCADE",
		13, "CUSTOM-FS" );

$go_menu = "FALSE";

while( $go_menu eq "FALSE" )

  {

  system( $clrstring );

  print "PulseShaper $version - GMB 2001\n";
  print "\n";
  print "Select the desired pulse profile from the following menu:\n";
  print "\n";
  printf "%3s%13s%19s%18s\n", "1)", "$PulseShapes{1} pulse", "2)", "$PulseShapes{2} pulse";
  printf "%3s%18s%14s%15s\n", "3)", "$PulseShapes{3} pulse", "4)", "$PulseShapes{4} pulse";
  printf "%3s%20s%12s%14s\n", "5)", "$PulseShapes{5} pulse", "6)", "$PulseShapes{6} pulse";
  printf "%3s%11s%21s%11s\n", "7)", "$PulseShapes{7} pulse", "8)", "$PulseShapes{8} pulse";
  printf "%3s%14s%18s%14s\n", "9)", "$PulseShapes{9} pulse", "10)", "$PulseShapes{10} pulse";
  printf "%3s%13s%19s%17s\n","11)", "$PulseShapes{11} pulse", "12)", "$PulseShapes{12} pulse";
  printf "%3s%22s%10s%19s\n","13)", "Custom Fourier series", "14)", "Frequency profiles";
  printf "%3s%22s\n", "15)", "Convert vendor shapes";
  print "\n";
  print "Enter pulse shape, or <q> to quit: ";

  chomp( $option = <STDIN> );

  if( ( $option eq "q" ) || ( $option eq "Q" ) )
    {
    $go_menu = "TRUE";
    }
  elsif( ( $option < 1 ) || ( $option > 15 ) )
    {
    $go_menu = "FALSE";
    }
  else
    {
    $ppname = $PulseShapes{$option};

    system( $clrstring );

    %profile = ();
    %profile_fft = ();
    %profile_re = ();
    %profile_im = ();
    %freqdata_re = ();
    %freqdata_im = ();

    $pi = 3.1415927;
    $pulseint = 0;
    $max = 0;

    unless( $option == 14 )

      {

      unless( $option == 15 )
	{
	print "PulseShaper $version - $PulseShapes{$option} pulse parameters\n";
	print "\n";
	print "Number of steps: ";
	chomp( $nsteps = <STDIN> );
	}
      else
	{
	print "PulseShaper $version - Bruker/Varian shaped pulse file conversion\n";
	print "\n";
	}

#
# Here enter additional parameters for certain pulses. Trapezoid...
#

      if( $option == 2 )
	{
	print "First point of trapezoid (0..$nsteps): ";
	chomp( $trap1 = <STDIN> );
	print "Last point of trapezoid (0..$nsteps): ";
	chomp( $trap2 = <STDIN> );
	}

#
# Triangle...
#

      if( $option == 3 )
	{
	print "Triangle maximum (0..$nsteps): ";
	chomp( $triang = <STDIN> );
	}

#
# Hermite...
#

      if( $option == 6 )
	{
	print "Power term coefficient: ";
	chomp( $hcoeff = <STDIN> );
	if( $hcoeff == 1 ) { $hcoeff = 1.0001 }
	}

#
# Gaussian, Half-gaussian, or Hermite...
#

      if( ( $option == 4 ) || ( $option == 5 ) || ( $option == 6 ) )
	{
	print "Truncation cutoff (%): ";
	chomp( $gausscutoff = <STDIN> );
	}

#
# Sinc...
#

      if( ( $option == 7 ) || ( $option == 8 ) )
	{
	print "Number of cycles or half-cycles: ";
	chomp( $sincycles = <STDIN> );
	}

#
# If the user decides to painstakingly enter a custom Fourier Series, let him do it...
#

      if( $option == 13 )
	{
	print "Number of terms in Fourier series: ";
	chomp( $fourierterms = <STDIN> );
	print "\n";
	print "Enter the terms following: A0 + SUM[ AN*cos(Nwt) + BN*sin(Nwt) ]\n";
	print "\n";
	print "A0: ";
	chomp( $cstmFS_A0 = <STDIN> );

	for( $i = 1; $i <= $fourierterms; ++$i )
	  {
	  print "A$i, B$i: ";
	  @Rows = split( /\,/, <STDIN> );
	  $cosFS[$i] = $Rows[0];
	  $sinFS[$i] = $Rows[1];
	  }

	print "\n";

	}

#
# Convert a Bruker/Varian shape file...
#

      if( $option == 15 )
	{
	print "Name of Bruker/Varian shape file: ";
	chomp( $vendorshapefile = <STDIN> );

	while( ! -e $vendorshapefile )
	  {
	  print "File $vendorshapefile not found. Enter a valid pulse profile: ";
	  chomp( $vendorshapefile = <STDIN> );
	  }

	}

#
# Finish by asking a name for the file...
#

      print "Name of PulseShaper shape file: ";
      chomp( $shapefile = <STDIN> );

#
# At this point we have everything we need to calculate the pulse and save the new shape file...
#

      unless( $option == 15 )
	{

	for( $i = 0; $i <= $nsteps - 1; ++$i )

	  {

	  if( $option == 1 )
            {
            $profile[$i] = 1;
            }

	  if( $option == 3 )
            {
            if( $i <= $triang )
	      {
	      $profile[$i] = $i / $triang;
	      }
            else
              {
	      $profile[$i] = 1 - ( $i - $triang ) / ( ( $nsteps - 1 ) - $triang );
	      }
            }

#
# Trapezoid...
#

	  if( $option == 2 )
            {
            if( $i < $trap1 )
 	      {
	      $profile[$i] = $i / $trap1;
	      }
            elsif( ( $i >= $trap1 ) && ( $i < $trap2 ) )
	      {
	      $profile[$i] = 1;
	      }
            else
	      {
	      $profile[$i] = 1 - ( $i - $trap2 ) / ( ( $nsteps - 1 )- $trap2 );
	      }
            }

#
# Gaussian...
#

	  if( $option == 4 )
            {
            $percent_to_constant = 4 * log( $gausscutoff / 100 ) / $nsteps ** 2;
            $profile[$i] = exp( $percent_to_constant * ( $i - ( $nsteps - 1 ) / 2 ) ** 2 );
            }

#
# Half-gaussian...
#

	  if( $option == 5 )
            {
            $percent_to_constant = log( $gausscutoff / 100 ) / $nsteps ** 2;
            $profile[$i] = exp( $percent_to_constant * ( $i - ( $nsteps - 1 ) ) ** 2 );
            }

#
# Hermite...
#

	  if( $option == 6 )
            {
            $percent_to_constant = 4 * log( ( $gausscutoff / 100 ) / abs( $hcoeff ** 2 - 1 ) ) / $nsteps ** 2;
            $profile[$i] = ( 1 - ( $hcoeff * ( 2 * $i / ( $nsteps - 1 ) - 1 ) ) ** 2 ) * 
                             exp( $percent_to_constant * ( $i - ( $nsteps - 1 ) / 2 ) ** 2 );
            }

#
# Sinc...
#

	  if( $option == 7 )
            {
            if( $i == $nsteps / 2 )
              {
	      $profile[$i] = sin( 2 * $sincycles * $pi * ( 1e-37 ) / $nsteps )
	                	  / ( $pi * ( 1e-37 ) / $nsteps );
	      }
            else
              { 
              $profile[$i] = sin( 2 * $sincycles * $pi * ( $i - ( $nsteps - 1 ) / 2 ) / $nsteps )
	                	  / ( $pi * ( $i - ( $nsteps - 1 ) / 2 ) / $nsteps );

              }

	   }

#
# Sine...
#

	  if( $option == 8 )
            {
            $profile[$i] = sin( $sincycles * $pi * $i / ( $nsteps - 1 ) );
            }

#
# E-BURP1
#

	  if( $option == 9 )
            {

            $omega = 2 * $pi / ( $nsteps - 1 );
            $profile[$i] = $omega * 0.230;

            for( $j = 1; $j <= 8; ++$j )
              {
              $profile[$i] = $profile[$i] + $omega * ( $EBurp1_cos{$j} * cos( $j * $omega * $i ) +
	                                               $EBurp1_sin{$j} * sin( $j * $omega * $i ) );
              }

            }

#
# E-BURP2...
#

	  if( $option == 10 )
            {

            $omega = 2 * $pi / ( $nsteps - 1 );
            $profile[$i] = $omega * 0.260;

            for( $j = 1; $j <= 10; ++$j )
              {
              $profile[$i] = $profile[$i] + $omega * ( $EBurp2_cos{$j} * cos( $j * $omega * $i ) +
	                                               $EBurp2_sin{$j} * sin( $j * $omega * $i ) );
              }

            }

#
# SNEEZE...
#

	  if( $option == 11 )
            {

            $omega = 2 * $pi / ( $nsteps - 1);
            $profile[$i] = $omega * 0.250;

            for( $j = 1; $j <= 10; ++$j )
              {
              $profile[$i] = $profile[$i] + $omega * ( $Sneeze_cos{$j} * cos( $j * $omega * $i ) +
	                                               $Sneeze_sin{$j} * sin( $j * $omega * $i ) );
              }

            }

#
# G4-CASCADE...
#

	  if( $option == 12 )
            {

            $omega = 2 * $pi / ( $nsteps - 1);
            $profile[$i] = $omega * 0.250;

            for( $j = 1; $j <= 11; ++$j )
              {
              $profile[$i] = $profile[$i] + $omega * ( $G4_cascade_cos{$j} * cos( $j * $omega * $i ) +
	                                               $G4_cascade_sin{$j} * sin( $j * $omega * $i ) );
              }

            }

#
# Custom Fourier series...
#

	  if( $option == 13 )
            {

            $omega = 2 * $pi / ( $nsteps - 1);
            $profile[$i] = $omega * $cstmFS_A0;

            for( $j = 1; $j <= $fourierterms; ++$j )
              {
              $profile[$i] = $profile[$i] + $omega * ( $cosFS[$j] * cos( $j * $omega * $i ) +
	                                               $sinFS[$j] * sin( $j * $omega * $i ) );
              }

            }

	  if( abs( $profile[$i] ) > $max )
            {
            $max = abs( $profile[$i] );
            }

	  $pulseint = $pulseint + $profile[$i];

	  }

	}

#
# Read in Bruker/Varian file. Read in the first line of the header and decide if
# it is Bruker or Varian...
#

      else

	{

	open( VENDORSHAPE, $vendorshapefile );
        @Rows = split( /\s+/, <VENDORSHAPE> );

#
# If the first character is a '#', its a Varian file...
#

        if( $Rows[0] eq "#" )
	  {
          $ppname = $Rows[1];
	  $nsteps = $Rows[4];

#
# Discard the empty row between header and pulse...
#

	  $line = <VENDORSHAPE>;

#
# In the Varian files we have two rows per pulse step. One is the reverse pulse (i.e.,
# the pulse maximum amplitude minus the the variation, line with a '0'), and the other
# one is the regular pulse profile (line with a '1023')...
#

          for( $i = 0; $i <= $nsteps - 1; ++$i )
	    {
	    $line = <VENDORSHAPE>;
  	    $line = <VENDORSHAPE>;
	    @Rows = split( /\s+/, $line );
	    $profile[$i] = $Rows[2];
	    $phase = $Rows[0];
	    if( $phase > 0 ) { $profile[$i] = -1 * $profile[$i] }
            }
	  }

#
# Not Varian, must be Bruker. Read the J-CAMP header first...
#

	else
	  {
	  for( $i = 2; $i <= 19; ++$i )
	    {
	    $line = <VENDORSHAPE>;
	    if( $i == 8 )
	      {
	      @Rows = split( /\s+/, $line );
	      $ppname = $Rows[2];
              }
	    if( $i == 18 )
	      {
	      @Rows = split( /\s+/, $line );
	      $nsteps = $Rows[1];
              }
	    }

#
# Now read in shape...
#        

	  for( $i = 0; $i <= $nsteps - 1; ++$i )
	    {
	    @Rows = split( /\,/, <VENDORSHAPE> );
	    $profile[$i] = $Rows[0];
	    $phase = $Rows[1];
	    if( $phase > 0 ) { $profile[$i] = -1 * $profile[$i] }
            }

	  }

	close( VENDORSHAPE );

#
# Now get the values of pulse integral and maximum intensity to normalize the profiles...
#

        for( $i = 0; $i <= $nsteps - 1; ++$i )
	  {
          if( abs( $profile[$i] ) > $max )
            {
            $max = abs( $profile[$i] );
            }
          $pulseint = $pulseint + $profile[$i];
          }
	  
	}     

#
# Now the profile is made or converted, normalize it and make the file. Junk it if
# it is already there...
#

      if( -e $shapefile ) { system( "$rmstring $shapefile" ) }

      open( OUTPUT, ">$shapefile" );

      print OUTPUT "$ppname\n";
      print OUTPUT "PULSES $nsteps\n";
      $pulseint = $pulseint / $max;
      print OUTPUT "SUM ";
      printf OUTPUT "%3.3f\n", $pulseint;

      for( $i = 0; $i <= $nsteps - 1; ++$i )
	{

	$profile[$i] = $profile[$i] / $max;

	if( $profile[$i] < 0 )
	  {
	  $profvalue = - 1 * $profile[$i];
	  $phase = "180";
	  }
	else
	  {
	  $profvalue = $profile[$i];
	  $phase = "0";
	  }
	printf OUTPUT "%0.3f", $profvalue;
	print OUTPUT " $phase\n";
	}

      close( OUTPUT );

      }

    else

      {

#
# If we selected to read a file, ask for the name and read the pulse profile...
#

      print "PulseShaper $version - Shaped pulse file selecion\n";
      print "\n";
      print "Enter the pulse profile file you want to read into memory. Provide the\n";
      print "entire path name to the file.\n";
      print "\n";
      print "Pulse profile file: ";

      chomp( $shapefile = <STDIN> );
      $fileok = "FALSE"; 

#
# Check if the file exists. Otherwise, loop...
#  

      while( ! -e $shapefile )
	{
	print "File $shapefile not found. Enter a valid pulse profile file: ";
	chomp( $shapefile = <STDIN> );
	}

#
# Now we have a valid file name. Read the damn thing...
#

      open( PPFILE, $shapefile );

      chomp( $ppname = <PPFILE> );

      @Rows = split( /\s+/, <PPFILE> );
      $nsteps = $Rows[1];
      @Rows = split( /\s+/, <PPFILE> );
      $pulseint = $Rows[1];

      for( $i = 0; $i <= $nsteps - 1; ++$i )
	{
	@Rows = split( /\s+/, <PPFILE> );
	$profile[$i] = $Rows[0];
	$phase = $Rows[1];
	if( $phase eq "180" ) { $profile[$i] = -1 * $profile[$i] }
	}

      close( PPFILE );

      }  

#
# After the file has been saved or read, allow the user to do some estimation of the 
# excitation bandwitdh depending on the pulse width. Calculate the FFT of the excitation
# pulse once - Do it before entering the Bloch loop, because it won't change...
#

#
# Zero-fill to the ( closest power of 2 of nsteps ) * 16 to get a nice profile. First calculate
# the ( closest power of 2 of nsteps ) * 16...
#

    $div = $nsteps / 2;
    $powers = 1;
    
    while( $div > 1 )
      {
      $div = $div / 2;
      ++$powers;
      }

    $ftsize = 2 ** ( $powers + 4 );

    for( $i = 0; $i <= $ftsize - 1; ++$i )
      {
      if( $i <= $nsteps - 1 )
	{
	$profile_fft[2*$i] = $profile[$i];
	$profile_fft[2*$i+1] = 0;
	}
      else
	{
	$profile_fft[2*$i] = 0;
	$profile_fft[2*$i+1] = 0;
	}
      }

#
# Destroy the '$data' object before using it so that we can change the
# size of the FFT without any problems...
#

    $data = shift;

    for( $i = 0; $i <= 2 * $ftsize - 1; ++$i ) { $data->[$i] = $profile_fft[$i] }

    $fft = new Math::FFT( $data );
    $coeff = $fft->cdft( );

#
# Retrieve the data and find the maximum...
#

    $fftprofmax = 0;
    $val = 0;
    
    for( $i = 0; $i <= $ftsize - 1; ++$i )
      {
      $freqdata_re[$i] = $coeff->[2*$i];
      $freqdata_im[$i] = $coeff->[2*$i+1];
      $val = sqrt( $freqdata_re[$i] ** 2 + $freqdata_im[$i] ** 2 );
      if( $val > $fftprofmax ) { $fftprofmax = $val }
      }

#
# Rotate FT array and normalize it...
#

    for( $i = 0; $i <= $ftsize / 2 - 1; ++$i )
      {
      $temp_re = $freqdata_re[$i];
      $temp_im = $freqdata_im[$i];

      $freqdata_re[$i] = $freqdata_re[$ftsize/2+$i] / $fftprofmax;
      $freqdata_im[$i] = $freqdata_im[$ftsize/2+$i] / $fftprofmax;

      $freqdata_re[$ftsize/2+$i] = $temp_re / $fftprofmax;
      $freqdata_im[$ftsize/2+$i] = $temp_im / $fftprofmax;

      }

#
# Calculate FT width at half-height. Do it by finding the first point that reaches the maximum and
# then double that around zero...
#

    $index = 0;
    $val = 0;
    $ftcounter = 0;
    
    while( $val <= 0.5 )
      {
      $val = sqrt( $freqdata_re[$index] ** 2 + $freqdata_im[$index] ** 2 );
      ++$index;
      }

    $ftcounter = 2 * ( $ftsize / 2 - $index );

    $go = "FALSE";

    while( $go eq "FALSE" )

      {

      system( $clrstring );

      print "PulseShaper $version - $ppname excitation profile calculation\n";
      print "\n";
      print "You can now calculate the frequency profile for the $ppname pulse\n";
      print "Enter the pulse width to get a frequency profile, or <q> to quit.\n";
      print "\n";
      print "Shaped pulse width (mS): ";
      chomp( $pwidth = <STDIN> );

      if( ( $pwidth eq "q" ) || ( $pwidth eq "Q" ) )
	{
	$go = "TRUE";
	}
      else
	{

        print "Shaped pulse tip angle (deg.): ";
        chomp( $tipangledeg = <STDIN> );
	$tipangle = $pi / 180 * $tipangledeg;

#
# Allow the user to start with magnetization anywere. Z-axis is the default...
#

        print "Initial magnetization [0 0 1]: ";
        chomp( $initmag = <STDIN> );

        if( $initmag eq "" )
	  {
	  $Mx_ini = 0;
  	  $My_ini = 0;
  	  $Mz_ini = 1;
	  }
	else
	  {
	  @Rows = split( /\s+/, $initmag );
	  
#
# Normalize whatever the user entered to a unit vector...
#	  
	  
	  $Mtot = sqrt( $Rows[0] ** 2 + $Rows[1] ** 2 + $Rows[2] ** 2 );
	  
	  $Mx_ini = $Rows[0] / $Mtot;
          $My_ini = $Rows[1] / $Mtot;
          $Mz_ini = $Rows[2] / $Mtot;
	  }

#
# Calculate the dwell time necessary for all offset estimations...
#        

	$dwtime = $pwidth / ( 1000 * $nsteps );

#
# Calculate bandwidth from width at half-height for the FT profile...
#

	$ftwidthathh = $ftcounter / ( $dwtime * $ftsize );

#
# Calculate gammaB1 / 2pi...
#
# gammaB1 (Hz) = gammaB1 / 2pi (rad), so:
# gammaB1 (Hz, tipangledeg pulse) = tipangle / ( 2 * pi * delta(pw) * sum{ shape(t) } ), and
# delta(pw) = pulsewidth / nsteps, so:
# gammaB1 (Hz, 90 deg. pulse) = nsteps * tipangle / ( 2 * pi * pulsewidth *  sum{ shape(t) } )...
#

	$gammaB1 = $nsteps * $tipangle / ( 2 * $pi * ( $pwidth / 1000 ) * $pulseint );

#
# Simulate Bloch equations. For the shaped pulse we have to apply the cummulative
# effects of every 'slice' of the pulse in the time domain as a short square pulse,
# and evaluate the effect of all the summed rotations (with varying gammaB1) for all
# the offsets we will have (the number determined by ftsize). Lets see if this works...
#

        %Mx = ();
	%My = ();
	%Mz = ();
	
	for( $i = 0; $i <= $ftsize - 1; ++$i )
	  
	  {
	  
	  $Mx[$i] = $Mx_ini;
  	  $My[$i] = $My_ini;
  	  $Mz[$i] = $Mz_ini;
	  	  
          $offset = 1 / $dwtime * ( $i - $ftsize / 2 ) / $ftsize;

  	  for(  $j = 0; $j <= $nsteps - 1; ++$j )
	    
	    {
	    
	    $gB1 = $gammaB1 * $profile[$j];
	    $theta = atan2 ( $gB1, $offset );
	    $beff = $tipangle * sqrt( $gB1 ** 2 + ( $offset ) ** 2 ) / ( $gammaB1 * $pulseint );

#
# With all the angles, do the rotations using analytical solution of the Bloch equations. The solutions
# are from:
#
# 1) Liu, Weisz, and James, JMR-A 1993, 105, 184-192.*
# 2) Lunati, Cofrancesco, Villa, Marzola, and Osculati, JMR 1998, 134, 223-235.
#
# * They have a small error in the rotation matrix - The last 'My' should be an 'Mz'...
#

	    $Mx_temp = ( ( sin( $theta ) ** 2 ) + ( cos( $theta ) ** 2 ) * cos( $beff ) ) * $Mx[$i] +
	               cos( $theta ) * sin( $beff ) * $My[$i] +
		       sin( $theta ) * cos( $theta ) * ( 1 - cos( $beff ) ) * $Mz[$i];
	    
	    $My_temp = -1 * cos( $theta ) * sin( $beff ) * $Mx[$i] +
	               cos( $beff ) * $My[$i] +
		       sin( $theta ) * sin( $beff ) * $Mz[$i];
		       
	    $Mz_temp = sin( $theta ) * cos( $theta ) * ( 1 - cos( $beff ) ) * $Mx[$i] -
	               sin( $theta ) * sin( $beff ) * $My[$i] +
		       ( ( cos( $theta ) ** 2 ) + ( sin( $theta ) ** 2 ) * cos( $beff ) ) * $Mz[$i];
		       
            $Mx[$i] = $Mx_temp;
	    $My[$i] = $My_temp;
	    $Mz[$i] = $Mz_temp;

	    }
	  
          }

#
# Now get the excitation bandwidth by estimating when the transverse magnetization
# (the square root of Mx^2+My^2) has droped to 70.8%...
#

        $index = 0;
        $val = 0;
        $excounter = 0;
    
        while( $val <= 0.708 )
          {
          $val = sqrt( $Mx[$index] ** 2 + $My[$index] ** 2 );
          ++$index;
          }

        $excounter = 2 * ( $ftsize / 2 - $index );
	$exwidthathh = $excounter / ( $dwtime * $ftsize );

#
# Save the frequency profile...
#

	$freqprofile = $shapefile. ".frq";
	if( -e $freqprofile ) { system( "$rmstring $freqprofile" ) }  

	open( FREQPROF, ">$freqprofile" );

	print FREQPROF "Hz\tReal\tImaginary\tPower\tMx\tMy\tMz\tMx^2+My^2\n";

	for( $i = 0; $i <= $ftsize - 1; ++$i )
	  {
	  $pow = sqrt( $freqdata_re[$i] ** 2 + $freqdata_im[$i] ** 2 );
	  $Mx2My2 = sqrt( $Mx[$i] ** 2 + $My[$i] ** 2 );
	  $freq = 1 / $dwtime * ( $i - $ftsize / 2 ) / $ftsize;
	  print FREQPROF "$freq\t$freqdata_re[$i]\t$freqdata_im[$i]\t$pow\t$Mx[$i]\t$My[$i]\t$Mz[$i]\t$Mx2My2\n";
	  }

	close( FREQPROF );

#
# Print some stuff regarding the profile...
#

	print "gammaB1 for $tipangledeg deg. pulse: ";
	printf "%0.2f", $gammaB1;
	print " Hz\n";
	print "Fourier bandwidth: ";
	printf "%0.2f", $ftwidthathh;
	print " Hz\n";
	print "Excitation bandwith: ";
	printf "%0.2f", $exwidthathh;
	print " Hz\n";
	print "Frequency profile saved as: $freqprofile\n";
	print "\n";
	print "Press <Enter> to continue...";

	$dummy = <STDIN>;

        }

      }

    }
    
  }

system( $clrstring );

print "Thanks for using PulseShaper! Send comments and suggestions to g.moyna\@usip.edu\n";
print "\n";

exit