/* nag_stat_pdf_vavilov (g01muc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.0, 2023.
 */

#include <nag.h>
#include <stdio.h>

int main(void) {

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

  exit_status = 0;

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

  printf(" nag_stat_pdf_vavilov (g01muc) Example Program Results\n\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  while (scanf("%lf%lf%lf%*[^\n] ", &x, &rkappa, &beta2) != EOF)
  {
    if ((rkappa != c1) || (beta2 != c2)) {
      /* nag_stat_init_vavilov (g01zuc).
       * Initialization function for
       * nag_stat_pdf_vavilov (g01muc) and
       * nag_stat_prob_vavilov (g01euc)
       */
      nag_stat_init_vavilov(rkappa, beta2, mode, &xl, &xu, comm_arr, &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_stat_init_vavilov (g01zuc).\n%s\n",
               fail.message);
        exit_status = 1;
        goto END;
      }
    }

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

    printf("   X     Rkappa    Beta2        Y\n\n");
    printf("  %3.1f      %3.1f      %3.1f   %13.4e\n", x, rkappa, beta2, y);
    c1 = rkappa;
    c2 = beta2;
  }
END:
  return exit_status;
}