/* D03PCF Example Program Text 
 * C Version.  NAG Copyright 1997, 2009, 2012. 
 * NAG C Header File example
 */

#include <stdio.h>
#include <math.h>
#include "nagmk23.h"

/* Global data */
double alpha;
extern void uinit(double uu[20][2], double *, Integer );
extern void NAG_CALL pdedef(const Integer *npde, const double *t,
                            const double *x, const double *u,
                            const double dudx[], double p[], double q[],
                            double r[], Integer *ires);
extern void NAG_CALL bndary(const Integer *npde, const double *t,
                            const double u[], const double *ux,
                            const Integer *ibnd, double beta[], double gamma[],
                            Integer *ires);

/* Routine for PDE initial condition */
void uinit(double u[20][2], double *x, Integer npts)
{
  Integer i;
  for (i=0; i<npts; i++)
    {
      u[i][0] = 2.0*alpha*x[i];
      u[i][1] = 1.0;
    }
}

/* Define the system of PDEs */
void NAG_CALL pdedef(const Integer *npde, const double *t, const double *x,
		     const double *u, const double dudx[], double p[],
		     double q[], double r[], Integer *ires)
{
/* Macro to do subscripting into 2 dimensional array p, which
   is passed as double p[]. */
#define P(I,J) p[(I)*(*npde) + (J)]
  q[0] = 4.0*alpha*(u[1] + *x * dudx[1]);
  q[1] = 0.0;
  r[0] = *x * dudx[0];
  r[1] = dudx[1] - u[0]*u[1];
  P(0,0) = 0.0;
  P(0,1) = 0.0;
  P(1,0) = 0.0;
  P(1,1) = 1.0 - (*x)*(*x);
}

void NAG_CALL bndary(const Integer *npde, const double *t, const double u[],
		     const double *ux, const Integer *ibnd, double beta[],
		     double gamma[], Integer *ires)
{
  if (*ibnd==0)
    {
      beta[0] = 0.0;
      beta[1] = 1.0;
      gamma[0] = u[0];
      gamma[1] = -u[0]*u[1];
    }
  else
    {
      beta[0] = 1.0;
      beta[1] = 0.0;
      gamma[0] = -u[0];
      gamma[1] = u[1];
    }
}

int main(void)
{
  const Integer npde=2, npts=20, intpts=6, itype=1, niw=64, nw=1128;
  const double acc = 1.0e-3;
  double hx, pi, piby2, tout, ts;
  const Integer m = 1, itask = 1, itrace = 0;
  Integer i, ifail, ind, j;
  /* u[npts][npde], uout[itype][intpts][npde], w[nw], x[npts] */
  double u[20][2], uout[1][6][2], w[1128], x[20];
  /* xout[intpts] */
  double xout[] = {0.0,0.4,0.6,0.8,0.9,1.0};
  Integer iw[64];
  void NAG_CALL bndary(const Integer *npde, const double *t, const double *u,
		       const double *ux, const Integer *ibnd, double *beta,
		       double *gamma, Integer *ires);
  void NAG_CALL pdedef(const Integer *npde, const double *t, const double *x,
		       const double *u, const double *dudx, double p[],
		       double *q, double *r, Integer *ires);
  void uinit(double u[20][2], double *x, Integer npts);
  
  printf("D03PCF Example Program Results\n\n");
  alpha = 1.0;
  ind = 0;

  /* Set spatial mesh points */

  piby2 = 0.5*x01aaf_(&pi);
  hx = piby2/(npts-1);
  x[0] = 0.0;
  x[npts-1] = 1.0;
  for (i=1; i<=npts-2; i++)
    x[i] = sin(hx*i);

  /* Set initial conditions */

  ts = 0.0;
  tout = 0.1e-4;
  printf("Accuracy requirement = %lf\n",acc);
  printf("Parameter alpha = %lf\n\n",alpha);
  printf("  T  /  X  ");
  for (i=0; i<=5; i++)
    printf("%8.4lf",xout[i]);
  printf("\n");

  /* Set the initial values */

  uinit(u,x,npts);
  for (i=0; i<=4; i++)
    {
      ifail = -1;
      tout = 10.0*tout;
      /* Cast 2-D array u to double * to match prototype */
      d03pcf_(&npde,&m,&ts,&tout,pdedef,bndary,(double *)u,&npts,x,&acc,w,
              &nw,iw,&niw,&itask,&itrace,&ind,&ifail);

      /* Interpolate at required spatial points */

      /* Cast 2-D array u and 3-D array uout 
         to double * to match prototypes */
      d03pzf_(&npde,&m,(double *)u,&npts,x,xout,&intpts,&itype,
              (double *)uout,&ifail);
      printf("%6.4lf U(1)",tout);
      for (j=0; j<intpts; j++)
        printf("%8.4lf",uout[0][j][0]);
      printf("\n       U(2)");
      for (j=0; j<intpts; j++)
        printf("%8.4lf",uout[0][j][1]);
      printf("\n");
    }
  printf("\n");

  /* Print integration statistics */

  printf("Number of integration steps in time                    %"NAG_IFMT"\n",iw[0]);
  printf("Number of residual evaluations of resulting ODE system %"NAG_IFMT"\n",iw[1]);
  printf("Number of Jacobian evaluations                         %"NAG_IFMT"\n",iw[2]);
  printf("Number of iterations of nonlinear solver               %"NAG_IFMT"\n",iw[4]);
  return 0;
}