/* f4d.c  f(x,y,z,t)=1.4*x^2 + 0.8*y^2 + z^3 -z -(t-0.5)^2 -1.0 */
/*        data for plot4d   f4d.dat                             */

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

static double xmin = -1.0;
static double xmax = 1.0;
static int nx = 9;
static double ymin = -1.0;
static double ymax = 1.0;
static int ny = 9;
static double zmin = -1.0;
static double zmax = 1.0;
static int nz = 9;
static double tmin = -1.0;
static double tmax = 1.0;
static int nt = 9;
static char fname[] = "f4d.dat";
static FILE* fout;

static double f(double x, double y, double z, double t) 
{
  return 1.4*x*x + 0.8*y*y + z*z*z - z - (t-0.5)*(t-0.5) - 1.0;
} /* end f */


int main(int argc, char* argv[])
{
  double x, y, z, t, hx, hy, hz, ht;
  int ix, iy, iz, it;

  printf("f4c.d writing f4d.dat \n");
  printf("f(x,y,z,t)=1.4*x^2 + 0.8*y^2 + z^3 -z -(t-0.5)^2 -1.0\n");

  hx = (xmax-xmin)/(double)(nx-1);
  hy = (ymax-ymin)/(double)(ny-1);
  hz = (zmax-zmin)/(double)(nz-1);
  ht = (tmax-tmin)/(double)(nt-1);
  printf("xmin=%f, xmax=%f, hx=%f, nx=%d \n", xmin, xmax, hx, nx);
  printf("ymin=%f, ymax=%f, hy=%f, ny=%d \n", ymin, ymax, hy, ny);
  printf("zmin=%f, zmax=%f, hz=%f, nz=%d \n", zmin, zmax, hz, nz);
  printf("tmin=%f, tmax=%f, ht=%f, nt=%d \n", tmin, tmax, ht, nt);

  fout = fopen(fname, "w");
  fprintf(fout, "%d %d %d %d %d\n", 4, nx, ny, nz, nt);
  for(ix=0; ix<nx; ix++)
  {
    x = xmin+hx*(double)ix;
    for(iy=0; iy<ny; iy++)
    {
      y = ymin+hy*(double)iy;
      for(iz=0; iz<ny; iz++)
      {
        z = zmin+hz*(double)iz;
        for(it=0; it<ny; it++)
        {
          t = tmin+ht*(double)it;
          fprintf(fout, "%f %f %f %f %f\n", x, y, z, t, f(x,y,z,t));
	}
      }
    }
  }
  fclose(fout);
  printf("f4d.dat written \n");
  return 0;
} /* end main of f4d.c */