/*	getv30.c

	Program to summarize velocity-profile files in Luke's 2003 format,
	placing a vertically averaged velocity value, by slowness, on the
	standard output.

	Run:	getv30 [depth] file.txt [file2.txt ...] >output

	Compile:	cc -o getv30 getv30.c

	Parameters:
		depth	averaging depth, meters, must be a float value
			>= zero, default 30.
*/

#include <strings.h>
#include <stdio.h>

struct Point {float dep; float vel;} apoint;
struct Profile
	{
	char file[360];	/* File name */
	char site[360];	/* Site name string */
	double lat;	/* Latitude */
	double lon;	/* Longitude */
	char date[360];	/* Date string */
	char perf[360];	/* Performed by string */
	char type[360];	/* Type of measurements string */
	char depth[360];	/* Depth of measurements string */
	int stepped;	/* Function stepped yes or no */
	char post[360];	/* Posted by string */
	float metdepth;	/* factor to convert depths to meters: 1, 1000, or 0.3048 */
	char vtype[10];	/* Velocity type: "Vs" or "Vp" */
	float metvel;	/* factor to convert velocities to meters: 1, 1000, or 0.3048 */
	int npoints;	/* number of points in data array */
	struct Point **vec;	/* pointer to profile data array */
	} *prof;
float depth; /* Depth of averaging */
int nfiles, nalloc;
struct Point **newvec;
float calcV();
struct Point *constPoint(), *newPoint();
struct Profile *constProfile();
	
main(ac, av)
int ac; char **av;
	{
	int ifile, ifile0, nread, ip;
	float depz, velocity;
	double lat, lon;
	char line[720];
	FILE *fp;

	/* Self document */
	if (ac < 2 || (ac == 2 && isNumber(av[1])))
		{
		fprintf(stderr, "getv30 [depth] file.txt [file2.txt ...] >output\n");
		exit(0);
		}

	/* Parse command line parameters and file names */
	if (isNumber(av[1]))
		{
		sscanf(av[1], "%f", &depth);
		if (depth <= 0.0 || depth > 7e6)
			{
			fprintf(stderr,
			"getv30: depth=%g, must be >0 and <7e6, abort\n", depth);
			exit(-1);
			}
		nfiles = ac -2;
		ifile0 = 2;
		}
	else
		{
		depth = 30.0;
		nfiles = ac - 1;
		ifile0 = 1;
		}

	/* Loop over files */
	nread = 0;
	nalloc = 0;
	prof = NULL;
	for (ifile=ifile0; ifile < ac; ifile++)
		{
		/* Open file */
		fp = fopen(av[ifile], "r");
		if (fp==NULL)
			{
			fprintf(stderr, "Can't open profile file %s, skip\n", av[ifile]);
			continue;
			}

		/* read each line from file */
		freeProfile();
		while( fgets(line, 719, fp) != NULL )
			{
			/* Remove trailing newline and return characters from line */
			rmnl(line, 720);
			/* Identify and read header lines */
			if (sscanf(line, "%f %f", &depz, &velocity) == 2
			&& prof != NULL)
				{
				/* Profile data line */
				/* Convert values to m and m/s */
				depz *= prof->metdepth;
				velocity *= prof->metvel;
				/* Add new profile point */
				addPoint(depz, velocity);
				}
			else if (prof == NULL)
				/* Header line, first in new profile */
				initHeaders(line, av[ifile]);
			else if (prof->npoints > 0)
				{
				/* Header line encountered after profile data,
				   must be a new profile in the same file */
				/* Compute depth-averaged velocity */
				/* Print file name, lat/long, average velocity to
				    standard out */
				report();
				/* Dispose of profile, init new profile */
				freeProfile();
				initHeaders(line, av[ifile]);
				}
			else
				/* Header line encountered before profile data */
				parseHeaders(line);
			}

		/* Close file */
		fclose(fp);
		nread++;
		if (prof != NULL)
			{
			/* file ended after profile data, report */
			/* Compute depth-averaged velocity */
			/* Print file name, lat/long, average velocity to
			    standard out */
			report();
			}
		}

	fprintf(stderr, "getv30: %d files could be read and processed.\n", nread);
	exit(nread);
	}

float calcV(depz)
float depz;
	{
	float V, sum, thick;
	int i;
	V = 0.0;
	if (prof == NULL)
		{
		fprintf(stderr, "getv30:calcV: Error: no profile data to summarize.\n");
		return(V);
		}
	if (prof->npoints < 1)
		{
		/*fprintf(stderr, "getv30:calcV: Error: no profile data points to summarize.\n");*/
		return(V);
		}
	sum = 0.0;
	for (i=0; i<prof->npoints; i++)
		{
		if (i == 0)
			{
			if (((prof->vec)[i])->dep > depz)
				{
				sum += depz/((prof->vec)[i])->vel;
				break;
				}
			else
				sum += ((prof->vec)[i])->dep/((prof->vec)[i])->vel;
			}
		else
			{
			if (((prof->vec)[i])->dep > depz)
				thick = depz - ((prof->vec)[i-1])->dep;
			else
				thick = ((prof->vec)[i])->dep - ((prof->vec)[i-1])->dep;
			sum += thick/((prof->vec)[i])->vel;
			sum += thick/2.0*
			(1.0/((prof->vec)[i])->vel - 1.0/((prof->vec)[i-1])->vel);
			if (((prof->vec)[i])->dep > depz)
				break;
			}
		}
	if (prof->vec[prof->npoints - 1]->dep < depz)
		sum += (depz - ((prof->vec)[prof->npoints - 1])->dep)
			/((prof->vec)[prof->npoints - 1])->vel;
	V = depz/sum;
	return(V);
	}

/* Allocate and initialize prof structure, parse the header line */
initHeaders(str, file)
char *str, *file;
	{
	prof = constProfile();
	strncpy(prof->file, file, 719);
	(prof->file)[719] = '\0';
	parseHeaders(str);
	}

/* Allocate and initialize prof structure, parse the header line */
parseHeaders(str)
char *str;
	{
	char word1[150], word2[150];
	if (prof == NULL)
		{
		fprintf(stderr,
		"getv30:parseHeaders: Error: no profile allocated for header line:\n%s\n",
		str);
		return(-1);
		}
	if (strncasecmp(str, "Site", 4) == 0)
		{
		strncpy(prof->site, str, 359);
		(prof->site)[359] = '\0';
		}
	if (strncasecmp(str, "Lat", 3) == 0)
		{
		if (sscanf(str, "%s %s", word1, word2) != 2)
			{
			fprintf(stderr,
			"getv30:parseHeaders: Error: could not parse Lat/long for header line in file %s :\n%s\n",
			prof->file, str);
			return(-1);
			}
		if (sscanf(word2, "%lf,%lf", &(prof->lat), &(prof->lon)) != 2)
			if (sscanf(word2, "%lf/%lf", &(prof->lat), &(prof->lon)) != 2)
				{
				fprintf(stderr,
				"getv30:parseHeaders: Error: could not parse Lat/long for header line in file %s :\n%s\n",
				prof->file, str);
				return(-1);
				}
		if ((prof->lat == 0.0 && prof->lon == 0.0)
		|| prof->lat < -90.0 || prof->lat > 90.0
		|| prof->lon < -180.0 || prof->lon > 180.0)
			{
			fprintf(stderr,
			"getv30:parseHeaders: Error: could not parse Lat/long for header line in file %s :\n%s\n",
			prof->file, str);
			return(-1);
			}
		}
	if (strncasecmp(str, "Date", 4) == 0)
		{
		strncpy(prof->date, str, 359);
		(prof->date)[359] = '\0';
		}
	if (strncasecmp(str, "Perf", 4) == 0)
		{
		strncpy(prof->perf, str, 359);
		(prof->perf)[359] = '\0';
		}
	if (strncasecmp(str, "Type", 4) == 0)
		{
		strncpy(prof->type, str, 359);
		(prof->type)[359] = '\0';
		}
	if (strncasecmp(str, "Depth of", 8) == 0)
		{
		strncpy(prof->depth, str, 359);
		(prof->depth)[359] = '\0';
		}
	if (strncasecmp(str, "Function stepped", 16) == 0)
		{
		if (strncasecmp(str+27, "y", 1) == 0)
			prof->stepped = 1;
		else if (strncasecmp(str+27, "yes", 3) == 0)
			prof->stepped = 1;
		else if (strncasecmp(str+27, "n", 1) == 0)
			prof->stepped = 0;
		else if (strncasecmp(str+27, "no", 2) == 0)
			prof->stepped = 0;
		else
			{
			fprintf(stderr,
			"getv30:parseHeaders: Error: could not parse Function stepped switch for header line in file %s :\n%s\nAssuming Function stepped.\n",
			prof->file, str);
			return(-1);
			}
		}
	if (strncasecmp(str, "Posted", 6) == 0)
		{
		strncpy(prof->post, str, 359);
		(prof->post)[359] = '\0';
		}
	if (strncasecmp(str, "Depth(", 6) == 0)
		{
		if (sscanf(str, "%s %s", word1, word2) != 2)
			{
			fprintf(stderr,
			"getv30:parseHeaders: Error: could not parse Vtype and units for header line in file %s :\n%s\nAssuming Vs and meters and m/s.\n",
			prof->file, str);
			return(-1);
			}
		if (strncasecmp(word1, "Depth(m)", 8) == 0)
			prof->metdepth = 1.0;
		else if  (strncasecmp(word1, "Depth(km)", 9) == 0)
			prof->metdepth = 1000.0;
		else if  (strncasecmp(word1, "Depth(ft)", 9) == 0)
			prof->metdepth = 0.3048;
		else
			{
			fprintf(stderr,
			"getv30:parseHeaders: Error: could not parse Depth units for header line in file %s :\n%s\nAssuming depths in meters.\n",
			prof->file, str);
			return(-1);
			}
		if (strncasecmp(word2, "Vs(m/s)", 7) == 0)
			{
			strcpy(prof->vtype, "Vs");
			prof->metvel = 1.0;
			}
		else if (strncasecmp(word2, "Vs(km/s)", 8) == 0)
			{
			strcpy(prof->vtype, "Vs");
			prof->metvel = 1000.0;
			}
		else if (strncasecmp(word2, "Vs(ft/s)", 8) == 0)
			{
			strcpy(prof->vtype, "Vs");
			prof->metvel = 0.3048;
			}
		else if (strncasecmp(word2, "Vp(m/s)", 7) == 0)
			{
			strcpy(prof->vtype, "Vp");
			prof->metvel = 1.0;
			}
		else if (strncasecmp(word2, "Vp(km/s)", 8) == 0)
			{
			strcpy(prof->vtype, "Vp");
			prof->metvel = 1000.0;
			}
		else if (strncasecmp(word2, "Vp(ft/s)", 8) == 0)
			{
			strcpy(prof->vtype, "Vp");
			prof->metvel = 0.3048;
			}
		else
			{
			fprintf(stderr,
			"getv30:parseHeaders: Error: could not parse Vtype and units for header line in file %s :\n%s\nAssuming Vs in m/s.\n",
			prof->file, str);
			return(-1);
			}
		}
	}

/* Print report on profile data */
report()
	{
	float V;
	if (prof == NULL)
		{
		fprintf(stderr, "getv30:report: Error: no profile data to summarize.\n");
		return(-1);
		}
	if (prof->npoints < 1)
		return(-1);
	V = calcV(depth);
	printf("%s\t%s\t%.6f\t%.6f\t%s%g= %.0f m/s\n",
		prof->file, prof->site, prof->lat, prof->lon, prof->vtype, depth, V);
	}

/* Free memory allocated to and referenced by the prof structure */
freeProfile()
	{
	if (prof != NULL)
		{
		int i;
		if (prof->vec != NULL)
			{
			for (i=0; i<nalloc; i++)
				{
				free((prof->vec)[i]);
				(prof->vec)[i] = NULL;
				}
			free(prof->vec);
			/*prof->vec = NULL;*/
			}
		free(prof);
		prof = NULL;
		}
	}

initPoint(point)
struct Point *point;
	{
	point->dep = -1.0;
	point->vel = -1.0;
	}
struct Point *constPoint()
	{
	struct Point *point;
	point = (struct Point *)(malloc(sizeof(struct Point)));
	allocerr(point, sizeof(struct Point), "constPoint:point");
	initPoint(point);
	return(point);
	}
setPoint(point, dep, vel)
struct Point *point;
float dep, vel;
	{
	point->dep = dep;
	point->vel = vel;
	}
struct Point *newPoint(dep, vel)
float dep, vel;
	{
	struct Point *point;
	point = constPoint();
	setPoint(point, dep, vel);
	return(point);
	}

/* Add a new point to prof->vec, extending the vector if needed */
addPoint(dep, vel)
float dep, vel;
	{
	int newalloc, i;
	if (prof->npoints > 0)
	if (dep < ((prof->vec)[prof->npoints - 1])->dep)
		{
		/* Profile data not in depth order, ignore */
		fprintf(stderr,
		"getv30:addPoint: Error: new profile point dep=%g vel=%g in file %s\n\t%s\n\thas depth less than previous point dep=%g vel=%g, ignoring this point.\n",
		dep, vel, prof->file, prof->site, ((prof->vec)[prof->npoints - 1])->dep,
		((prof->vec)[prof->npoints - 1])->vel);
		return(-1);
		}
	if (vel <= 0.0)
		{
		/* Profile velocity value not >0, ignore */
		fprintf(stderr,
		"getv30:addPoint: Error: new profile point dep=%g vel=%g in file %s\n\t%s\n\thas a velocity <= zero, ignoring this point.\n",
		dep, vel, prof->file, prof->site, ((prof->vec)[prof->npoints - 1])->dep,
		((prof->vec)[prof->npoints - 1])->vel);
		return(-1);
		}
	/* if sufficient pre-allocated points, set value to new profile point */
	if (prof->npoints < nalloc)
		setPoint((prof->vec)[prof->npoints], dep, vel);
	else
		{
		/* Allocate an extended profile vector */
		newalloc = nalloc + 300;
		newvec = (struct Point **)(malloc(newalloc*sizeof(struct Point *)));
		allocerr(newvec, newalloc*sizeof(struct Point *), "addPoint:newvec");
		/* Copy profile->vec data to new vector */
		for (i=0; i<nalloc; i++)
			{
			newvec[i] = newPoint(((prof->vec)[i])->dep,
				((prof->vec)[i])->vel);
			free((prof->vec)[i]);
			}
		/* assign prof-vec to new vector */
		free(prof->vec);
		prof->vec = newvec;
		/* pre-allocate the remainder of the new vector */
		for (i=nalloc; i<newalloc; i++)
			(prof->vec)[i] = constPoint();
		nalloc = newalloc;
		}
	(prof->npoints)++;
	}

initProfile(prf)
struct Profile *prf;
	{
	int i;
	sprintf(prf->file, "");
	sprintf(prf->site, "");
	prf->lat = 0.0;
	prf->lon = 0.0;
	sprintf(prf->date, "");
	sprintf(prf->perf, "");
	sprintf(prf->type, "");
	sprintf(prf->depth, "");
	prf->stepped = 1;
	sprintf(prf->post, "");
	prf->metdepth = 1.0;
	sprintf(prf->vtype, "Vs");
	prf->metvel = 1.0;
	prf->npoints = 0;
	nalloc = 300;
	prf->vec = (struct Point **)(malloc(nalloc*sizeof(struct Point *)));
	allocerr(prf->vec, nalloc*sizeof(struct Point *), "initProfile:prf->vec");
	for (i=0; i<nalloc; i++)
		(prf->vec)[i] = constPoint();
	}
struct Profile *constProfile()
	{
	struct Profile *prf;
	prf = (struct Profile *)(malloc(sizeof(struct Profile)));
	allocerr(prf, sizeof(struct Profile), "constProfile:prf");
	initProfile(prf);
	return(prf);
	}

int isNumber(str)
char *str;
	{
	float f;
	int i;
	if (sscanf(str, "%f", &f) > 0)
		return(1);
	if (sscanf(str, "%d", &i) > 0)
		return(1);
	return(0);
	}
rmnl(str, n)
char *str;
int n;
	{
	int i;
	for (i=0; i<n; i++)
		if (str[i] == '\n' || str[i] == '\r')
			str[i] = '\0';
	}
		
allocerr(p, nbytes, s)			/* error handling for malloc */
int p, nbytes;
char s[120];
	{
	if(p==NULL)
		{
		fprintf(stderr, "allocerr: cannot allocate %d bytes of memory for %s\n",
		nbytes, s);
		if (nbytes > 2000000)
			fprintf(stderr, "Don't be so greedy!\n");
		exit(-1);
		}
	return(0);
	}