/* nag_1d_pade (e02rac) Example Program.

 *

 * Copyright 2001 Numerical Algorithms Group.

 *

 * Mark 7, 2001.

 * Mark 7b revised, 2004.

 */



#include <stdio.h>

#include <nag.h>

#include <nag_stdlib.h>

#include <nagc02.h>

#include <nage02.h>



int main(void)

{

  /* Scalars */

  Integer exit_status, i, l, m, ia, ib, ic;

  NagError fail;



  /* Arrays */

  double *aa = 0, *bb = 0, *cc = 0, *dd = 0;

  Complex *z = 0;



  INIT_FAIL(fail);

  exit_status = 0;

  Vprintf("nag_1d_pade (e02rac) Example Program Results\n");



  l = 4;

  m = 4;

  ia = l + 1;

  ib = m + 1;

  ic = ia + ib - 1;



  /* Allocate memory */

  if ( !(aa = NAG_ALLOC(ia, double)) ||

       !(bb = NAG_ALLOC(ib, double)) ||

       !(cc = NAG_ALLOC(ic, double)) ||

       !(dd = NAG_ALLOC(ia + ib, double)) ||

       !(z = NAG_ALLOC(l+m, Complex)) )

    {

      Vprintf("Allocation failure\n");

      exit_status = -1;

      goto END;

    }



  /* Power series coefficients in cc */

  cc[0] = 1.0;

  for (i = 1; i <= 8; ++i)

    cc[i] = cc[i-1] / (double) i;



  Vprintf("\n");



  Vprintf("The given series coefficients are\n");



  for (i = 1; i <= ic; ++i)

    {

      Vprintf("%13.4e", cc[i-1]);

      Vprintf(i%5 == 0 || i == ic ? "\n" : " ");

    }



  /* nag_1d_pade (e02rac).

   * Pade-approximants

   */

  nag_1d_pade(ia, ib, cc, aa, bb, &fail);

  if (fail.code != NE_NOERROR)

    {

      Vprintf("Error from nag_1d_pade (e02rac).\n%s\n", fail.message);

      exit_status = 1;

      goto END;

    }





  Vprintf("\n");

  Vprintf("Numerator coefficients\n");



  for (i = 1; i <= ia; ++i)

    {

      Vprintf("%13.4e", aa[i-1]);

      Vprintf(i%5 == 0 || i == ia ?"\n":" ");

    }



  Vprintf("\n");

  Vprintf("Denominator coefficients\n");



  for (i = 1; i <= ib; ++i)

    {

      Vprintf("%13.4e", bb[i-1]);

      Vprintf(i%5 == 0 || i == ib ? "\n" : " ");

    }



  /* Calculate zeros of the approximant using nag_zeros_real_poly (c02agc) */

  /* First need to reverse order of coefficients */

  for (i = 1; i <= ia; ++i)

    dd[ia-i] = aa[i-1];



  /* nag_zeros_real_poly (c02agc).

   * Zeros of a polynomial with real coefficients

   */

  nag_zeros_real_poly(l, dd, Nag_TRUE, z, &fail);

  if (fail.code != NE_NOERROR)

    {

      Vprintf("Error from c02agc, 1st call.\n%s\n", fail.message);

      exit_status = 1;

      goto END;

    }



  Vprintf("\n");

  Vprintf("Zeros of approximant are at\n");

  Vprintf("    Real part    Imag part\n");

  for (i = 1; i <= l; ++i)

    Vprintf("%13.4e%13.4e\n", z[i-1].re, z[i-1].im);



  /* Calculate poles of the approximant using nag_zeros_real_poly (c02agc) */

  /* Reverse order of coefficients */

  for (i = 1; i <= ib; ++i)

    dd[ib-i] = bb[i-1];



  /* nag_zeros_real_poly (c02agc), see above. */

  nag_zeros_real_poly(m, dd,  Nag_TRUE, z, &fail);

  if (fail.code != NE_NOERROR)

    {

      Vprintf("Error from c02agc, 2nd call.\n%s\n", fail.message);

      exit_status = 1;

      goto END;

    }



  Vprintf("\n");

  Vprintf("Poles of approximant are at\n");

  Vprintf("    Real part    Imag part\n");

  for (i = 1; i <= m; ++i)

    Vprintf("%13.4e%13.4e\n", z[i-1].re, z[i-1].im);



 END:

  if (aa) NAG_FREE(aa);

  if (bb) NAG_FREE(bb);

  if (cc) NAG_FREE(cc);

  if (dd) NAG_FREE(dd);

  if (z) NAG_FREE(z);



  return exit_status;

}