/* Kaiser Window function */

/*Usage: -fn MAC -macro KW.M -var kw value offset value
    kw is the Kaiser window constant, typically a multiple of PI/2
    offset is the value of the first sampling point, typically 0 or 0.5
*/

/*Bessel function constants*/
p1 = 1.0e1;
p2 = 3.5156229e1;
p3 = 3.0899424e1;
p4 = 1.2067492e1;
p5 = 0.2659732e1;
p6 = 0.360768e-1;
p7 = 0.45813e-2;

q1 = 0.39894228e1;
q2 = 0.1328592e-1;
q3 = 0.225319e-2;
q4 = 0.157565e-2;
q5 = 0.916281e-2;
q6 = 0.2057706e-1;
q7 = 0.2635537e-1;
q8 = 0.1647633e-1;
q9 = 0.392377e-2;

/* spectra specific parameters*/
if( yLoc == 1) {
/*    (void) printf( "Kaiser window function parameter: %.4f \n", kw );
    (void) printf( "Kaiser window offset parameter: %.4f \n", offset );
*/
    };

kw = abs(kw);

/*Calculate normalization*/
    if( kw < 3.75 )
      {
        w = (kw/3.75)^2;
        norm = p1 + w*(p2 + w*(p3 + w*(p4 + w*(p5 + w*(p6 + w*p7)))));
       }

    else
       {
        y = 3.75/kw;
        norm = (exp(kw)/sqrt(kw))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9))))))));
       };

/*Begin actual calculations*/
   for(i = 0; i < xSize; i++) {

    x = kw * sqrt(1.0 - ((i+offset)/(xSize+offset))^2);

    if( x < 3.75 )
      {
        w = (x/3.75)^2;
        b = p1 + w*(p2 + w*(p3 + w*(p4 + w*(p5 + w*(p6 + w*p7)))));
       }

    else
       {
        y = 3.75/x;
        b = (exp(x)/sqrt(x))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9))))))));
       };
    g = b / norm;
    rdata[i] *= g;
    idata[i] *= g;
    if (yLoc == 1) { 
        /* (void) printf( " %.0f  %.6f\n", i, g); */
        }; 
   };

