/* -------------------------------------------------------------------

   -----  d02cbf.c  -----

  Seiwi's 2nd tiny example for using NAG Fortran-Lib with gcc 2.7.x

  Comments from NAG:

  D02CJF ODEs, IVP, Adams method, until function of solution is zero,
         intermediate output (simple driver)
  (D02CJF seems to do the same as D02CJF)

  This demo is the same as d02cbfe.f

  To compile with:
                   f90 d02cbf.c -o d02cbf -lnag -lm

  Tested with: NAG fllux19d9 and NAGWare_f90

  You should get the following output after running d02cbf:

D02CBF Example Program Results
Calculation with TOL = 0.000100
    X         Y(1)         Y(2)         Y(3)
   0.00      0.00000      0.50000      0.62832
   1.00      0.62716      0.44642      0.48440
   2.00      1.04927      0.40548      0.30662
   3.00      1.25713      0.38056      0.09771
   4.00      1.24211      0.37433     -0.12885
   5.00      0.99519      0.38729     -0.35136
   6.00      0.50618      0.41727     -0.55062
   7.00     -0.23638      0.46048     -0.71787
   8.00     -1.24534      0.51297     -0.85370
Calculation with TOL = 0.000010
    X         Y(1)         Y(2)         Y(3)
   0.00      0.00000      0.50000      0.62832
   1.00      0.62716      0.44642      0.48440
   2.00      1.04931      0.40548      0.30662
   3.00      1.25733      0.38055      0.09768
   4.00      1.24229      0.37433     -0.12890
   5.00      0.99500      0.38731     -0.35143
   6.00      0.50558      0.41730     -0.55068
   7.00     -0.23705      0.46052     -0.71791
   8.00     -1.24587      0.51299     -0.85371


------------------------------------------------------------------- */

#include <f90.h>

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


struct {
  double xend;
  double h;
  int i;
} COMMON;


extern double x01aaf_();
extern void d02cbf_();

void fcn(double *, double *, double *);
void out(double *, double *);



int main(int argc, char *argv[])
{
  double Tmp1;
  double Tmp2;
  int Tmp3;
  auto double pi_;
  auto double tol;
  auto double x;
  auto int ifail;
  auto int ir;
  auto int j;
  auto double w[86];
  auto double y[3];


  f90_init(argc,argv);

  printf("D02CBF Example Program Results\n");

  /* PI from NAG Fortran-lib  */
  Tmp2 = 0.0;
  pi_ = x01aaf_(&Tmp2);

  ir = 0;

  COMMON.xend = 8.0;

  for(j = 4; j <= 5; j++) {
    tol = pow( 10.,-(j));
    printf("Calculation with TOL = %lf\n", tol);
    printf("    X         Y(1)         Y(2)         Y(3)\n");

    x = 0.0e+00;
    y[0] = 0.0;
    y[1] = 0.5;
    y[2] = pi_/ 5.0;

    COMMON.i = 7;
    COMMON.h = (COMMON.xend - x)/(double)((COMMON.i + 1));
    ifail = 0;

    Tmp3 = 3;
    d02cbf_(&x, &COMMON.xend, &Tmp3, y, &tol, &ir, fcn, out, w, &ifail);
  }

  f90_finish(0);

  return(0);
}




void fcn(double *t, double y[3], double f[3])
{
  double Tmp1;

  f[0] = tan(y[2]);
  f[1] = -0.032*tan(y[2])/y[1] - 0.02*y[1]/cos(y[2]);

  Tmp1 = y[1];
  f[2] = -0.032/(Tmp1*Tmp1);

  return;
}



void out(double *x, double y[3])
{
  double Tmp1;
  auto int j;

  printf("%7.2lf", *x);

  for(j = 0; j < 3; j++) {
   printf("%13.5lf", y[j]);
  }
  printf("\n");

 *x = COMMON.xend - (double)((COMMON.i))*COMMON.h;

 COMMON.i--;   /* = COMMON.i - 1;*/

 return;
}