/* nag_prob_density_vavilov (g01muc) Example Program.
 *
 * Copyright 2002 Numerical Algorithms Group.
 *
 * Mark 7, 2002.
 */

#include <stdio.h>
#include <nag.h>
#include <nagx04.h>
#include <nag_stdlib.h>
#include <nagg01.h>
#include <nagx02.h>

int main(int argc, char *argv[])
{
  FILE     *fpin, *fpout;

  /* Scalars */
  double   c1, c2, x, rkappa, beta2, xl, xu, y;
  Integer  exit_status, mode;
  NagError fail;

#define WKMAX 322

  double   comm_arr[WKMAX];

  mode = 0;

  INIT_FAIL(fail);

  /* Check for command-line IO options */
  fpin = nag_example_file_io(argc, argv, "-data", NULL);
  fpout = nag_example_file_io(argc, argv, "-results", NULL);
  exit_status = 0;

  /* nag_real_largest_number (x02alc).
   * The largest positive model number
   */
  c1 = -X02ALC;
  /* nag_real_largest_number (x02alc), see above. */
  c2 = -X02ALC;

  fprintf(fpout,
          " nag_prob_density_vavilov (g01muc) Example Program Results\n\n");

  /* Skip heading in data file */
  fscanf(fpin, "%*[^\n] ");

  while (fscanf(fpin, "%lf%lf%lf%*[^\n] ", &x, &rkappa, &beta2) != EOF)
    {
      if ((rkappa != c1) || (beta2 != c2))
        {
          /* nag_init_vavilov (g01zuc).
           * Initialization function for
           * nag_prob_density_vavilov (g01muc) and
           * nag_prob_vavilov (g01euc)
           */
          nag_init_vavilov(rkappa, beta2, mode, &xl, &xu, comm_arr, &fail);
          if (fail.code != NE_NOERROR)
            {
              fprintf(fpout, "Error from nag_init_vavilov (g01zuc).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }
        }

      /* nag_prob_density_vavilov (g01muc).
       * Vavilov density function phi_V((lambda~;~kappa~)beta^2)
       */
      y = nag_prob_density_vavilov(x, comm_arr);

      fprintf(fpout, "   X     Rkappa    Beta2        Y\n\n");
      fprintf(fpout, "  %3.1f      %3.1f      %3.1f   %13.4e\n", x, rkappa,
              beta2, y);
      c1 = rkappa;
      c2 = beta2;
    }
 END:
  if (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  return exit_status;
}