/*	vel2d.c:
 *	  Creates starting 2-D velocity profile as a 3D model. Model
 *        originates in upper left corner, and variation is in x and z.
 *        Specify x-locataions along the profile where velocity depth
 *        profiles are specified....
 *
 *	ari, oct 02
 *
 *	V 1.2
 *	  Option to build the model in y-direction rather than in x....
 *
 *	ari, feb 04
 *
 *	V 1.2.3
 *	  Small fix to allow vertical columns to be "double" (like in horizontal layering) - to allow
 *	  horizontal discontinuities as well.
 *
 *	ari, sep 11
 */

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

#define VER 1.2
#define FIX 3
#define AUTHOR "Ari Tryggvason"
#define DATE "110908"

#define MAX 250
#define MAXL 25
#define NINT(x) (int) floor((x)+0.5)

char           *progname;
int main(argc, argv)
   int             argc;
   char          **argv;
{
   char           *progname,
                   string[MAX], vfname[MAX], cfname[MAX];

   int             i, j, k, l, m, nx, ny, nz, nxy, nxyz, nc, nd, yes, index,
                   xflg, jnk_i;

   float           dx, dy, dz, x0, y0, z0, xmax, ymax, zmax, *x, this_x, this_z,
                  *vel, *zl, *vl, *zu, *vu, v1, v2, zxu, zxl, velocity, jnk_f;

   FILE           *f1_out, *f2_out;

   int             get_line();

/* command line argument(s) */
   progname = *argv;

fprintf(stderr,"%s V %.1f.%1d by %s %s\n", progname, VER, FIX, AUTHOR, DATE);

/* read form stdin */
fprintf(stdout,"Output 2-D velocity file name: ");
get_line(stdin,string, MAX);
sscanf(string,"%s", cfname);
fprintf(stdout,"Output 3-D velocity file name: ");
get_line(stdin,string, MAX);
sscanf(string,"%s", vfname);
fprintf(stdout,"Spatial sample interval in km (x y z): ");
get_line(stdin,string, MAX);
sscanf(string,"%f %f %f", &dx, &dy, &dz);
fprintf(stdout,"Origin of model in km (x y z): ");
get_line(stdin,string, MAX);
sscanf(string,"%f %f %f", &x0, &y0, &z0);
fprintf(stdout,"Opposite corner of model in km (x y z): ");
get_line(stdin,string, MAX);
sscanf(string,"%f %f %f", &xmax, &ymax, &zmax);
fprintf(stdout,"Velocities to vary in X (1) or Y (0) direction: ");
get_line(stdin,string, MAX);
sscanf(string,"%d", &xflg);
fprintf(stdout,"Number of depth nodes and columns: ");
get_line(stdin, string, MAX);
sscanf(string,"%d %d", &nd, &nc);
fprintf(stdout, "x0=%.1f y0=%.1f z0=%.1f\n", x0, y0, z0);
fprintf(stdout, "xm=%.1f ym=%.1f zm=%.1f\n", xmax, ymax, zmax);
fprintf(stdout, "dx=%.1f dy=%.1f dz=%.1f\n", dx, dy, dz);
fprintf(stdout, "nc=%d nd=%d\n", nc, nd);

 FILE *f;
 f = fopen("new2.txt", "a");

nx = NINT((xmax - x0)/dx);
ny = NINT((ymax - y0)/dy);
nz = NINT((zmax - z0)/dz);


if (!xflg)
{
  jnk_i = nx;
  nx = ny;
  ny = jnk_i;
  jnk_f = x0;
  x0 = y0;
  y0 = jnk_f;
  jnk_f = xmax;
  xmax = ymax;
  ymax = jnk_f;
  jnk_f = dx;
  dx = dy;
  dy = jnk_f;
}
nxy = nx * ny;
nxyz = nx * ny * nz;
if ((vel = (float *) malloc(nxyz*sizeof(float))) == (float *) NULL)
{
   fprintf(stderr,"%s: Couldn't allocate memory for vel.\n", progname);
   return (-1);
}

if ((vu = (float *) malloc(nc*sizeof(float))) == (float *) NULL)
{
   fprintf(stderr,"%s: Couldn't allocate memory.\n", progname);
   return (-1);
}
if ((vl = (float *) malloc(nc*sizeof(float))) == (float *) NULL)
{
   fprintf(stderr,"%s: Couldn't allocate memory.\n", progname);
   return (-1);
}
if ((zu = (float *) malloc(nc*sizeof(float))) == (float *) NULL)
{
   fprintf(stderr,"%s: Couldn't allocate memory.\n", progname);
   return (-1);
}
if ((zl = (float *) malloc(nc*sizeof(float))) == (float *) NULL)
{
   fprintf(stderr,"%s: Couldn't allocate memory.\n", progname);
   return (-1);
}
if ((x = (float *) malloc(nc*sizeof(float))) == (float *) NULL)
{
   fprintf(stderr,"%s: Couldn't allocate memory.\n", progname);
   return (-1);
}

/* open files */
if ((f1_out = fopen(vfname, "w")) == (FILE *) NULL)
{
   fprintf(stderr, "%s: Can't open file %s.\n", progname, vfname);
   return (-1);
}
if ((f2_out = fopen(cfname, "w")) == (FILE *) NULL)
{
   fprintf(stderr, "%s: Can't open file %s.\n", progname, cfname);
   return (-1);
}

for (j=0; j<nc; j++)
{
   fscanf(stdin,"%f", &x[j]);
   x[j] -= x0;
}
fprintf(stdout, "Node locations (x):\n");
for (i=0; i<nc; i++)
   fprintf(stdout,"%.1f\t", x[i]);
fprintf(stdout, "\nVelocity profiles (z v):\n");
/* read first layer before entering loop */
for (j=0; j<nc; j++)
{
   fscanf(stdin,"%f %f", &zl[j], &vl[j]);
   zl[j] -= z0;
   fprintf(stdout, "%4.1f %4.1f\t", zl[j]+z0, vl[j]);
}
fprintf(stdout,"\n");

for (i=1; i<nd; i++)
{ /* repeat over all layers */
   for (j=0; j<nc; j++)
   { /* lower boundary is now top boundary */
      zu[j] = zl[j];
      vu[j] = vl[j];
   }
   for (j=0; j<nc; j++)
   {
      fscanf(stdin,"%f %f", &zl[j], &vl[j]);
      zl[j] -= z0;
      fprintf(stdout, "%4.1f %4.1f\t", zl[j]+z0, vl[j]);
   }
   fprintf(stdout,"\n");
   for (k=0,j=1; k<nx; k++)
   {
      this_x = k*dx;
      if (this_x > x[j])
         j++;
	  /* check a second time */
      if (this_x > x[j])
         j++;
      v1 = vu[j-1] + (this_x - x[j-1]) * (vu[j] - vu[j-1]) / (x[j] - x[j-1]);
      v2 = vl[j-1] + (this_x - x[j-1]) * (vl[j] - vl[j-1]) / (x[j] - x[j-1]);
      zxu = zu[j-1] + (this_x - x[j-1]) * (zu[j] - zu[j-1]) / (x[j] - x[j-1]);
      zxl = zl[j-1] + (this_x - x[j-1]) * (zl[j] - zl[j-1]) / (x[j] - x[j-1]);
      for (m=0; m<nz; m++)
      {
         this_z = m*dz;
         if (this_z >= zxu && this_z < zxl)
         {
            velocity = v1 + (this_z - zxu) * (v2 - v1) / (zxl - zxu);
/* fprintf(stdout, "this_x=%f, x[%d]=%f\n", this_x, j-1, x[j-1]); */
            if (this_x == x[j-1])
               fprintf(f2_out, "%.1f %.1f %.1f\n", this_x, this_z, velocity);
            for (l=0; l<ny; l++)
            {
	       if (xflg){
                  index = m * nxy + l*nx + k;
                  fprintf(f,"%5d\n",index);
	      } else {
	          index = m * nxy + k*ny + l;
               vel[index] = velocity;
		      }
            }
         }
      }
   }

}  /* closer for (i) - all layers*/
yes = fwrite(vel, nxyz*sizeof(float), 1, f1_out);
if (!xflg)
{
   jnk_i = nx;
   nx = ny;
   ny = jnk_i;
}
if (yes != 1) fprintf(stdout,"%s: Hey, fwrite error, code %d.\n", progname, yes);
fprintf(stdout,"\nnx = %d\nny = %d\nnz = %d\n", nx, ny, nz);
return (0);
}				/* the end */