/*
	trigrid.c

	tg1.bat = cl /AL trigrid.c proj1 getfile3

	trigrid [input file name] [output file name] <options>

	where '.lbs' file format is:

		FILE_TYPE          = ELEVATION_ARRAY
		IMAGE_LINES        =    120
		LINE_SAMPLES       =    120
		UPPER_LEFT_LAT     =   36.25000
		UPPER_LEFT_LON     = -122.75000
		LOWER_RIGHT_LAT    =   36.01200 (optional)
		LOWER_RIGHT_LON    = -122.51200 (optional)
		DELTA_LAT          =   500 points/degree
		DELTA_LON          =   500 points/degree
		NUMBER_OF_LAT      =    120
		NUMBER_OF_LON      =    120
		MAXIMUM_Z          =     -2021.94263 (optional)
		MINIMUM_Z          =     -3482.91040 (optional)
		SPLINE_FLEXIBILITY = 0.20000 (optional)
		MISSING_DATA_VALUE = 0.000000 (optional)
		DATA_TYPE          = float
		DATA_SIZE          = 4 bytes
		DATA_FILE_POINTER  = 'm102.D64'
		INPUT_DATA_NAME    = 'pacific.D32' (optional)

*/
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "proj1.h"
#include "nomouse.h"

#define DEG_RAD 0.01745329251				/* converts degrees to radians */
#define E_RAD 6400
#define E_CIR 40212
#define ROW 2060
#define COL 4320

char Text[20][80];
char *Path[10];
int NumPath=0;
struct control Map;
struct grid_stat
{
	double maxlat,minlat,maxlon,minlon;
	double maxz,minz;
	double dellat,dellon;			/* array size */
	double latspace,lonspace;	/* line spacing */
	double scale;
	int numlat,numlon;
	char datafile[100];
	char sourcefile[100];
	double flex;
   double msng_val;
	char msng;
	int type,size;	/* type -- 1 = unsigned char (1 byte) */
                  /* type -- 2 = char (1 byte) */
                  /* type -- 3 = unsigned short (2 bytes) */
                  /* type -- 4 = short (2 bytes) */
                  /* type -- 5 = unsigned long (4 bytes) */
                  /* type -- 6 = long (4 bytes) */
                  /* type -- 7 = float (4 bytes) */
                  /* type -- 8 = double (8 bytes) */
}Fs;
union Data
{
	unsigned char uc[32768];
	char c[32768];
	short us[16384];
	short s[16384];
	unsigned long ul[8192];
	long l[8192];
	float f[8192];
	double d[4096];
}D;
double huge Z[3][4096];
long N[3][4096];
int Num_Data=0;
int Elev[COL],E16[COL];
unsigned char E8[COL];
char CBuff[6481];

int scan_height(double lat1,double lat2,double lon1,double lon2,
				long *max,long *min);
int wf_flat(double lat1,double lat2,double lon1,double lon2,
		double clat,double r,long dmin,long dmax,char *name,double km_in,
		double bottom_thick);
int read_lbs(char *name,struct grid_stat *fs);
int make_wvf(char *name,struct grid_stat fs);
int write_fs(FILE *fpw,struct grid_stat fs);
int read_data(FILE *fpdat,struct grid_stat fs,double *Z);
FILE *fopen_path(char *filename,char *type);

main(int argc,char *argv[])

{
	int i,j,k;
	char inname[100],outname[100];

	NumPath=argc;
	for(i=1;i<NumPath;i++)
		Path[i]=argv[i];
	if(argc>1)
	{
		strcpy(inname,argv[1]);
	}
	else if(get_file_name("*.lbs",inname)<0)
		exit(0);
	for(k=0;k<strlen(inname);k++)
		if(inname[k]=='.')
			inname[k]='\0';
	strcpy(outname,inname);
	for(k=0;k<strlen(outname);k++)
		if(outname[k]=='.')
			outname[k]='\0';
   read_lbs(inname,&Fs);
	make_wvf(outname,Fs);		
	printf("Use program 'WAVEPREP' to process data for program 'WAVE'\n");
/*
	printf("\n\n--- waveprep %s ---\n\n",outname);
	execlp("waveprep.exe","waveprep",outname,NULL);
	printf("\n\nCOULD NOT FIND PROGRAM 'WAVEPREP'\n\n\n");
	exit(0);*/
}



/**********************************************************************
**
**
**
**********************************************************************/

int make_wvf(char *name,struct grid_stat fs)

{
	int i,j,k,l,m,n;
	int idm1,idp1,id,numt;
	long row=fs.numlat,col=fs.numlon,r,c,ij[2][2],ijk[2][2][2];
	FILE *fptxt,*fpw,*fpv,*fpdat,*fpwvf;
	char string[100],cr=0xd;
	double lat,lon,x,y,z;
	double x1,x2,mx,y1,y2,my;	
	double xscale,yscale,scale,vscale;
	double midlat=(fs.maxlat+fs.minlat)/2.0,midlon=(fs.maxlon+fs.minlon)/2.0;
	double mz=(fs.maxz+fs.minz)/2.0;
	double thick=0.0;
	long nv=0,tnv=0;

	fpdat=fopen_path(fs.datafile,"rb");
	if(!fpdat)
	{
		printf("Could not open '%s' to read.\n\n",fs.datafile);
		exit(0);
	}
	fptxt=fopen_path("proj1.txt","w+t");
	if(!fptxt)
	{
		printf("Could not open 'proj1.txt'\n");
		exit(0);
	}
	if(choose_projection(&Map,Text[0])<0)
	{
		printf("No projection chosen.\n\n");
		exit(0);
	}
	fprintf(fptxt,"%s",Text[0]);
	get_constants(&Map,Text,&numt);
	for(k=0;k<numt;k++)
		fprintf(fptxt,"%s",Text[k]);
	sprintf(string,"%s.wvf",name);
	fpwvf=fopen_path(string,"wt");
	if(!fpwvf)
	{
		printf("Could not open '%s' to write.\n\n");
		exit(0);
	}
	sprintf(string,"%s.w",name);
	fpw=fopen_path(string,"wt");
	if(!fpw)
	{
		printf("Could not open '%s' to write.\n\n");
		exit(0);
	}
	sprintf(string,"%s.v",name);
	fpv=fopen_path(string,"wt");
	if(!fpv)
	{
		printf("Could not open '%s' to write.\n\n");
		exit(0);
	}
	rewind(fptxt);
	write_fs(fpv,fs);
	fprintf(fpv,"# MAP PROJECTION INFORMATION:\n");
	while(fgets(string,100,fptxt))
		fprintf(fpv,"# %s",string);
	fclose(fptxt);
	ll2xy((double)fs.maxlat,(double)fs.minlon,&x1,&y1,Map);
	ll2xy((double)fs.minlat,(double)fs.maxlon,&x2,&y2,Map);
	mx=(x1+x2)/2.0;
	my=(y1+y2)/2.0;
	xscale=9.0/(x2-x1);
	yscale=9.0/(y1-y2);
	scale=xscale<yscale?xscale:yscale;
	ll2xy((double)fs.maxlat,(double)fs.minlon,&x,&y,Map);
	x1=x2=x;y1=y2=y;
	ll2xy((double)fs.minlat,(double)fs.maxlon,&x,&y,Map);
	if(x<x1)x1=x;if(y>y1)y1=y;if(x>x2)x2=x;if(y<y2)y2=y;
	ll2xy((double)fs.minlat,(double)fs.minlon,&x,&y,Map);
	if(x<x1)x1=x;if(y>y1)y1=y;if(x>x2)x2=x;if(y<y2)y2=y;
	ll2xy((double)fs.maxlat,(double)fs.maxlon,&x,&y,Map);
	if(x<x1)x1=x;if(y>y1)y1=y;if(x>x2)x2=x;if(y<y2)y2=y;
	ll2xy((double)fs.maxlat,midlon,&x,&y,Map);
	if(x<x1)x1=x;if(y>y1)y1=y;if(x>x2)x2=x;if(y<y2)y2=y;
	ll2xy((double)fs.minlat,midlon,&x,&y,Map);
	if(x<x1)x1=x;if(y>y1)y1=y;if(x>x2)x2=x;if(y<y2)y2=y;
	ll2xy(midlat,(double)fs.minlon,&x,&y,Map);
	if(x<x1)x1=x;if(y>y1)y1=y;if(x>x2)x2=x;if(y<y2)y2=y;
	ll2xy(midlat,(double)fs.maxlon,&x,&y,Map);
	if(x<x1)x1=x;if(y>y1)y1=y;if(x>x2)x2=x;if(y<y2)y2=y;
	xscale=9.0/(x2-x1);
	yscale=9.0/(y1-y2);
	scale=xscale<yscale?xscale:yscale;
	printf("\n\n\n\n\n");
	printf("The figure has been scaled to be about 9 inches in size and\n");
	printf("centered at the origin.  If you wish to have an absolute\n");
	printf("coordinate system so that several pieces will fit together or two\n");
	printf("figures overlay each other, you can reset the SCALE and the X0, Y0\n");
	printf("and Z0 to match the other figures.\n");
	printf("\n");	
	printf("Scale = %lf km/in  (xyz will be in inches)\n",1.0/scale);
	printf("X0 = %11.5lf km\n",mx);
	printf("Y0 = %11.5lf km\n",my);
	printf("Z0 = %11.5lf m\n",mz);
	printf("Do you want to change these?  (y or n)\n");
	if(getch()=='y')
	{
		printf("Give new scale.\n");
		scanf("%lf",&scale);
		scale=1.0/scale;
		printf("Give new X0.\n");
		scanf("%lf",&my);
		printf("Give new Y0.\n");
		scanf("%lf",&my);
		printf("Give new Z0.\n");
		scanf("%lf",&mz);
	}
	fprintf(fpv,"# Scale = %lf km/in  (xyz below in inches)\n",1.0/scale);
	fprintf(fpv,"# X0 = %11.5lf km\n",mx);
	fprintf(fpv,"# Y0 = %11.5lf km\n",my);
	fprintf(fpv,"# Z0 = %11.5lf m\n",mz);
	if(fs.maxz!=0.0||fs.minz!=0.0)
		fprintf(fpv,"# Z min = %11.5f m  Z max = %11.5f m  Z range = %11.5f m\n",
			fs.minz,fs.maxz,fs.maxz-fs.minz);
	printf("\n\nGive Vertical Exaggeration.\n\n");
	scanf("%lf",&vscale);	
	fprintf(fpv,"# Vertical Exaggeration = %lf\n",vscale);
	fprintf(fpv,"# x = (X - X0)/Scale\n");
	fprintf(fpv,"# y = (Y - Y0)/Scale\n");
	fprintf(fpv,"# z = (Z - Z0)/Scale*Vert.Exagg./1000.0\n");
	fprintf(fpv,"# \n");
	fprintf(fpv,"# Use program 'mapsplin' to start the process\n");
	fprintf(fpv,"# 'mapsplin' calls 'splnbin' which does the spline fitting\n");
	fprintf(fpv,"# 'splnbin' calls 'trigrid' which makes the wave front file\n");
	fprintf(fpv,"# 'trigrid' calls 'waveprep' which prepares it for viewing\n");
	fprintf(fpv,"# Use program 'wave' to view the data\n");
	fprintf(fpv,"#\n");
	vscale/=1000.0;	/* convert meters to km */
	printf(
"\n\nDo you want a two sided surface (volume rather than surface)?  (y or n)\n");
	if(getch()=='y')
	{
		printf("Give thickness in inches (1 in = %lf km).\n",scale);
		scanf("%lf",&thick);
		if(thick<=0.0)
		{
			thick=0.0;
			printf("invalid thickness -- figure will be single sided.\n\n");
		}
		if(thick>0.0)
			fprintf(fpv,
		"# Double sided figure thickness = %lf in (%lf km)\n",thick,thick*scale);
	}
	printf("Writing Verticies i 'V' and tiangles in 'W'\n");
	lat=fs.maxlat;
	for(i=0;i<3;i++)
	{
		for(j=0;j<fs.numlon;j++)
		{
			Z[i][j]=fs.msng_val;
			N[i][j]=-1;
		}
	}
	read_data(fpdat,fs,Z[0]);	/* read 1st line */
	for(i=0;i<fs.numlat;i++)
	{
		idm1=(i+2)%3;
		id=(i+3)%3;
		idp1=(i+4)%3;
		printf("%c%4d of %d",cr,i+1,fs.numlat);
		lon=fs.minlon;
		if(i<fs.numlat-1)	/* read next line */
			read_data(fpdat,fs,Z[idp1]);
		else			/* load line after last line with missing value */
			for(j=0;j<fs.numlon;j++)
			{
				Z[idp1][j]=fs.msng_val;
				N[idp1][j]=-1;
			}
		for(j=0;j<fs.numlon;j++)
		{
			tnv+=1;
			if(Z[id][j]!=fs.msng_val)	/* valid point */
			{
				if(							/* point has neighbors */
					(j<fs.numlon-1&&					/* top right */
					 Z[idm1][j  ]!=fs.msng_val&&
					 Z[idm1][j+1]!=fs.msng_val&&
					 Z[id  ][j+1]!=fs.msng_val)||
					(j<fs.numlon-1&&					/* bottom right */
					 Z[idp1][j  ]!=fs.msng_val&&
					 Z[idp1][j+1]!=fs.msng_val&&
					 Z[id  ][j+1]!=fs.msng_val)||
					(j>0&&								/* top left */
					 Z[idm1][j  ]!=fs.msng_val&&
					 Z[idm1][j-1]!=fs.msng_val&&
					 Z[id  ][j-1]!=fs.msng_val)||
					(j>0&&								/* bottom left */
					 Z[idp1][j  ]!=fs.msng_val&&
					 Z[idp1][j-1]!=fs.msng_val&&
					 Z[id  ][j-1]!=fs.msng_val))
				{
					ll2xy(lat,lon,&x,&y,Map);
					z=Z[id][j];
					N[id][j]=nv++;
					fprintf(fpv,"v  %lf %lf %lf\n",
						(x-mx)*scale,(y-my)*scale,(z-mz)*scale*vscale);
					if(thick>0.0)
						fprintf(fpv,"v  %lf %lf %lf\n",
						(x-mx)*scale,(y-my)*scale,(z-mz)*scale*vscale-thick);
					if(j>0&&				/* top left box present */
						Z[idm1][j  ]!=fs.msng_val&&
						Z[idm1][j-1]!=fs.msng_val&&
						Z[id  ][j-1]!=fs.msng_val)
					{
						if(thick==0.0)
						{
							fprintf(fpw,"f %ld %ld %ld\n",
								N[idm1][j  ]+1,N[idm1][j-1]+1,N[id  ][j-1]+1);
							fprintf(fpw,"f %ld %ld %ld\n",
								N[idm1][j  ]+1,N[id  ][j-1]+1,N[id  ][j  ]+1);
						}
						else
						{
							fprintf(fpw,"f %ld %ld %ld\n",
								N[idm1][j  ]*2+1,N[idm1][j-1]*2+1,N[id  ][j-1]*2+1);
							fprintf(fpw,"f %ld %ld %ld\n",
								N[idm1][j  ]*2+1,N[id  ][j-1]*2+1,N[id  ][j  ]*2+1);
							fprintf(fpw,"f %ld %ld %ld\n",
								N[idm1][j-1]*2+2,N[idm1][j  ]*2+2,N[id  ][j-1]*2+2);
							fprintf(fpw,"f %ld %ld %ld\n",
								N[id  ][j-1]*2+2,N[idm1][j  ]*2+2,N[id  ][j  ]*2+2);
						}
					}
					if(thick!=0.0)	/* do edges */
					{
						if((i>0&&j<fs.numlon-1&&	/* top right box present */
							 Z[idm1][j  ]!=fs.msng_val&&
							 Z[idm1][j+1]!=fs.msng_val&&
							 Z[id  ][j+1]!=fs.msng_val)&&
							(j==0||				/* top left box missing */
							 Z[idm1][j-1]==fs.msng_val||
							 Z[id  ][j-1]==fs.msng_val))/* edge left side */
						{
							fprintf(fpw,"f %ld %ld %ld\n",
								N[idm1][j  ]*2+1,N[idm1][j  ]*2+2,N[id  ][j  ]*2+2);
							fprintf(fpw,"f %ld %ld %ld\n",
								N[idm1][j  ]*2+1,N[id  ][j  ]*2+2,N[id][j  ]*2+1);
						}
						if((j>0&&	/* top left box present */
							 Z[idm1][j  ]!=fs.msng_val&&
							 Z[idm1][j-1]!=fs.msng_val&&
							 Z[id  ][j-1]!=fs.msng_val)&&
							(j==fs.numlon-1||				/* top right box missing */
							 Z[idm1][j+1]==fs.msng_val||
							 Z[id  ][j+1]==fs.msng_val))/* edge right side */
						{
							fprintf(fpw,"f %ld %ld %ld\n",
								N[idm1][j  ]*2+2,N[idm1][j  ]*2+1,N[id  ][j  ]*2+2);
							fprintf(fpw,"f %ld %ld %ld\n",
								N[id  ][j  ]*2+2,N[idm1][j  ]*2+1,N[id][j  ]*2+1);
						}
						if((j>0&&	/* top left box present */
							 Z[idm1][j  ]!=fs.msng_val&&
							 Z[idm1][j-1]!=fs.msng_val&&
							 Z[id  ][j-1]!=fs.msng_val)&&
							(i==fs.numlat-1||		/* bottom left box missing */
							 Z[idp1][j  ]==fs.msng_val||
							 Z[idp1][j-1]==fs.msng_val))/* edge bottom side */
						{
							fprintf(fpw,"f %ld %ld %ld\n",
								N[id  ][j  ]*2+1,N[id  ][j-1]*2+1,N[id  ][j-1]*2+2);
							fprintf(fpw,"f %ld %ld %ld\n",
								N[id  ][j  ]*2+1,N[id  ][j-1]*2+2,N[id  ][j  ]*2+2);
						}
						if((j>0&&	/* bottom left box present */
							 Z[idp1][j  ]!=fs.msng_val&&
							 Z[idp1][j-1]!=fs.msng_val&&
							 Z[id  ][j-1]!=fs.msng_val)&&
							(				/* top left box missing */
							 Z[idm1][j  ]==fs.msng_val||
							 Z[idm1][j-1]==fs.msng_val))/* edge top side */
						{
							fprintf(fpw,"f %ld %ld %ld\n",
								N[id  ][j-1]*2+1,N[id  ][j  ]*2+1,N[id  ][j-1]*2+2);
							fprintf(fpw,"f %ld %ld %ld\n",
								N[id  ][j-1]*2+2,N[id  ][j  ]*2+1,N[id  ][j  ]*2+2);
						}
					}
				}
			}
			lon+=fs.lonspace;
		}
		lat-=fs.latspace;
	}
	if(fs.msng!=0)
	{
		if(thick==0.0)
			printf("\n\nThere are %ld valid vertexes out of %ld\n",nv,tnv);
		else
			printf("\n\nThere are %ld valid vertexes out of %ld\n",nv*2,tnv*2);
	}
	else
		printf("/n/nThere are %ld vertiexes\n",nv,tnv);
	printf("Creating 'WVF' file.\n\n");
	fclose(fpv);
	fclose(fpw);
	fclose(fpdat);	
	sprintf(string,"%s.v",name);
	fpv=fopen_path(string,"rt");
	if(!fpv)
	{
		printf("Could not open '%s' to read.\n\n");
		exit(0);
	}
	while(fgets(string,100,fpv))
		fprintf(fpwvf,"%s",string);
	fclose(fpv);
	sprintf(string,"%s.w",name);
	fpw=fopen_path(string,"rt");
	if(!fpw)
	{
		printf("Could not open '%s' to read.\n\n");
		exit(0);
	}
	while(fgets(string,100,fpw))
		fprintf(fpwvf,"%s",string);
	flushall();
	fclose(fpw);
	printf("Completed '%s.wvf'\n",name);
	return(1);
}

/**********************************************************************
**
**
**
**********************************************************************/

int read_data(FILE *fpdat,struct grid_stat fs,double *z)

{
	int j;

	if(fs.type==1)
	{
		if(fread((char *)D.uc,sizeof(char),fs.numlon,fpdat)<fs.numlon)
		{
			printf("EOF hit in '%s'\n\n",fs.datafile);
			exit(0);
		}
		for(j=0;j<fs.numlon;j++)
			z[j]=D.uc[j];
	}
	else if(fs.type==2)
	{
		if(fread((char *)D.c,sizeof(char),fs.numlon,fpdat)<fs.numlon)
		{
			printf("EOF hit in '%s'\n\n",fs.datafile);
			exit(0);
		}
		for(j=0;j<fs.numlon;j++)
			z[j]=D.c[j];
	}
	else if(fs.type==3)
	{
		if(fread((char *)D.us,sizeof(short),fs.numlon,fpdat)<fs.numlon)
		{
			printf("EOF hit in '%s'\n\n",fs.datafile);
			exit(0);
		}
		for(j=0;j<fs.numlon;j++)
			z[j]=D.us[j];
	}
	else if(fs.type==4)
	{
		if(fread((char *)D.s,sizeof(short),fs.numlon,fpdat)<fs.numlon)
		{
			printf("EOF hit in '%s'\n\n",fs.datafile);
			exit(0);
		}
		for(j=0;j<fs.numlon;j++)
			z[j]=D.s[j];
	}
	else if(fs.type==5)
	{
		if(fread((char *)D.ul,sizeof(long),fs.numlon,fpdat)<fs.numlon)
		{
			printf("EOF hit in '%s'\n\n",fs.datafile);
			exit(0);
		}
		for(j=0;j<fs.numlon;j++)
			z[j]=D.ul[j];
	}
	else if(fs.type==6)
	{
		if(fread((char *)D.l,sizeof(long),fs.numlon,fpdat)<fs.numlon)
		{
			printf("EOF hit in '%s'\n\n",fs.datafile);
			exit(0);
		}
		for(j=0;j<fs.numlon;j++)
			z[j]=D.l[j];
	}
	else if(fs.type==7)
	{
		if(fread((char *)D.f,sizeof(float),fs.numlon,fpdat)<fs.numlon)
		{
			printf("EOF hit in '%s'\n\n",fs.datafile);
			exit(0);
		}
		for(j=0;j<fs.numlon;j++)
			z[j]=D.f[j];
	}
	else if(fs.type==8)
	{
		if(fread((char *)D.d,sizeof(double),fs.numlon,fpdat)<fs.numlon)
		{
			printf("EOF hit in '%s'\n\n",fs.datafile);
			exit(0);
		}
		for(j=0;j<fs.numlon;j++)
			z[j]=D.d[j];
	}
	else
	{
		printf("\n\nInvalid data type = %d\n",fs.type);
		exit(0);
	}
}



/**********************************************************************
**
**
**
**********************************************************************/

int write_fs(FILE *fpv,struct grid_stat fs)

{
	fprintf(fpv,
"# File made by cubic spline interpolation using program 'splnbin'\n");
	fprintf(fpv,"# and program 'trigrid'\n");
	fprintf(fpv,"# DATA ARRAY INFORMATION:\n");
	fprintf(fpv,"# Upper Left Corner  %10.5f %10.5f\n",fs.maxlat,fs.minlon);
	fprintf(fpv,"# Lower Right Corner %10.5f %10.5f\n",fs.minlat,fs.maxlon);
	fprintf(fpv,"# Latitude Range    = %10.5f\n",fs.dellat);
	fprintf(fpv,"# Longitude Range   = %10.5f\n",fs.dellon);
	fprintf(fpv,"# Latitude Spacing  = %10.8f\n",fs.latspace);
	fprintf(fpv,"# Longitude Spacing = %10.8f\n",fs.lonspace);
	fprintf(fpv,"# Number of Lat Lines (Rows)    = %5d\n",fs.numlat);
	fprintf(fpv,"# Number of Lon Lines (Columns) = %5d\n",fs.numlon);
	if(fs.msng>(char)0)
		fprintf(fpv,"# Missing Value = %lf\n",fs.msng_val);
	if(fs.flex>-1.0)
		fprintf(fpv,"# Spline Fleibility = %f\n",fs.flex);
	if(fs.type==1)fprintf(fpv,"# Data Type = 'unsigned char (1 byte)'\n");
	if(fs.type==2)fprintf(fpv,"# Data Type = 'char (1 byte)'\n");
	if(fs.type==3)fprintf(fpv,"# Data Type = 'unsigned short (2 bytes)'\n");
	if(fs.type==4)fprintf(fpv,"# Data Type = 'short (2 bytes)'\n");
	if(fs.type==5)fprintf(fpv,"# Data Type = 'unsigned long (4 bytes)'\n");
	if(fs.type==6)fprintf(fpv,"# Data Type = 'long (4 bytes)'\n");
	if(fs.type==7)fprintf(fpv,"# Data Type = 'float (4 bytes)'\n");
	if(fs.type==8)fprintf(fpv,"# Data Type = 'double (8 bytes)'\n");
	fprintf(fpv,"# Vector Source File = '%s'\n",fs.sourcefile);
	fprintf(fpv,"# Data File = '%s'\n",fs.datafile);
}



/**********************************************************************
**
**
**
**********************************************************************/

int read_lbs(char *name,struct grid_stat *fs)

{
	int i,j,k;
	FILE *fp;
	char string[100],str1[50],str2[50],str3[50],str4[50];
	char *type;

	sprintf(string,"%s.lbs",name);
	fp=fopen_path(string,"rt");
	if(!fp)
	{
		printf("Could not open '%s' to read label information.\n\n",string);
		exit(0);
	}
	fs->maxlat=-360.0;
	fs->maxlon=-360.0;
	fs->minlat=360.0;
	fs->minlon=360.0;
	fs->dellat=0.0;
	fs->dellon=0.0;
	fs->latspace=0.0;
	fs->lonspace=0.0;
	fs->scale=0.0;
	fs->numlat=0;
	fs->numlon=0;
	fs->size=0;
	fs->type=0;
	fs->flex=-1.0;
	fs->msng_val=-999999.0;
	fs->msng=0;
	while(fgets(string,100,fp))
	{
		sscanf(string,"%s%s%s%s",str1,str2,str3,str4);
		if(strcmpi(str1,"MISSING_DATA_VALUE")==0)
		{
			if(sscanf(str3,"%lf",&fs->msng_val)<1)
			{
				printf("Invalid MISSING_DATA_VALUE '%s'\n\n",str3);
				exit(0);
			}
			fs->msng=1;
			printf("MISSING_DATA_VALUE = %lf\n",fs->msng_val);
		}
		if(strcmpi(str1,"IMAGE_LINES")==0)
		{
			if(sscanf(str3,"%d",&fs->numlat)<1)
			{
				printf("Invalid IMAGE_LINES '%s'\n\n",str3);
				exit(0);
			}
			printf("IMAGE_LINES = %d\n",fs->numlat);
		}
		if(fs->numlat==0&&strcmpi(str1,"NUMBER_OF_LAT")==0)
		{
			if(sscanf(str3,"%d",&fs->numlat)<1)
			{
				printf("Invalid NUMBER_OF_LAT'%s'\n\n",str3);
				exit(0);
			}
			printf("NUMBER_OF_LAT = %d\n",fs->numlat);
		}
		if(strcmpi(str1,"LINE_SAMPLES")==0)
		{
			if(sscanf(str3,"%d",&fs->numlon)<1)
			{
				printf("Invalid LINE_SAMPLES '%s'\n\n",str3);
				exit(0);
			}
			printf("LINE_SAMPLES = %d\n",fs->numlon);
		}
		if(fs->numlon==0&&strcmpi(str1,"NUMBER_OF_LON")==0)
		{
			if(sscanf(str3,"%d",&fs->numlon)<1)
			{
				printf("Invalid NUMBER_OF_LON'%s'\n\n",str3);
				exit(0);
			}
			printf("NUMBER_OF_LON = %d\n",fs->numlon);
		}
		if(strcmpi(str1,"UPPER_LEFT_LAT")==0)
		{
			if(sscanf(str3,"%lf",&fs->maxlat)<1)
			{
				printf("Invalid UPPER_LEFT_LAT '%s'\n\n",str3);
				exit(0);
			}
			printf("UPPER_LEFT_LAT = %lf\n",fs->maxlat);
		}
		if(strcmpi(str1,"UPPER_LEFT_LON")==0)
		{
			if(sscanf(str3,"%lf",&fs->minlon)<1)
			{
				printf("Invalid UPPER_LEFT_LON '%s'\n\n",str3);
				exit(0);
			}
			printf("UPPER_LEFT_LON = %lf\n",fs->minlon);
		}
		if(strcmpi(str1,"LOWER_RIGHT_LAT")==0)
		{
			if(sscanf(str3,"%lf",&fs->minlat)<1)
			{
				printf("Invalid LOWER_RIGHT_LAT '%s'\n\n",str3);
				exit(0);
			}
			printf("LOWER_RIGHT_LAT = %lf\n",fs->minlat);
		}
		if(strcmpi(str1,"LOWER_RIGHT_LON")==0)
		{
			if(sscanf(str3,"%lf",&fs->maxlon)<1)
			{
				printf("Invalid LOWER_RIGHT_LON '%s'\n\n",str3);
				exit(0);
			}
			printf("LOWER_RIGHT_LON = %lf\n",fs->maxlon);
		}
		if(strcmpi(str1,"DELTA_LAT")==0)
		{
			if(sscanf(str3,"%lf",&fs->latspace)<1||fs->latspace<=0.0)
			{
				printf("Invalid DELTA_LAT '%s'\n\n",str3);
				exit(0);
			}
			if(strcmpi(str4,"points/degree")==0)
			fs->latspace=1.0/fs->latspace;
			printf("DELTA_LAT = %lf\n",fs->latspace);
		}
		if(strcmpi(str1,"DELTA_LON")==0)
		{
			if(sscanf(str3,"%lf",&fs->lonspace)<1||fs->lonspace<=0.0)
			{
				printf("Invalid DELTA_LON '%s'\n\n",str3);
				exit(0);
			}
			if(strcmpi(str4,"points/degree")==0)
			fs->lonspace=1.0/fs->lonspace;
			printf("DELTA_LON = %lf\n",fs->lonspace);
		}
		if(strcmpi(str1,"MAXIMUM_Z")==0)
		{
			if(sscanf(str3,"%lf",&fs->maxz)<1)
			{
				printf("Invalid MAXIMUM_Z '%s'\n\n",str3);
				exit(0);
			}
			printf("MAXIMUM_Z = %lf\n",fs->maxz);
		}
		if(strcmpi(str1,"MINIMUM_Z")==0)
		{
			if(sscanf(str3,"%lf",&fs->minz)<1)
			{
				printf("Invalid MINIMUM_Z '%s'\n\n",str3);
				exit(0);
			}
			printf("MINIMUM_Z = %lf\n",fs->minz);
		}
		if(strcmpi(str1,"DATA_FILE_POINTER")==0)
		{
			k=0;
			for(i=0;i<strlen(str3);i++)
				if(str3[i]!=39)		/* single quote */
					fs->datafile[k++]=str3[i];
			fs->datafile[k]='\0';
			printf("DATA_FILE_POINTER = '%s'\n",fs->datafile);
		}
		if(strcmpi(str1,"DATA_TYPE")==0)
		{
			strcpy(type,str3);
			printf("DATA_TYPE = '%s'\n",type);
		}
		if(strcmpi(str1,"DATA_SIZE")==0)
		{
			if(sscanf(str3,"%d",&fs->size)<1)
			{
				printf("Invalid DATA_SIZE '%s'\n\n",str3);
				exit(0);
			}
			printf("DATA_SIZE = %d\n",fs->size);
		}
		if(strcmpi(str1,"SPLINE_FLEXIBILITY")==0)
		{
			sscanf(str3,"%lf",&fs->flex);
			printf("SPLINE_FLEXIBILITY = %lf\n",fs->flex);
		}
		if(strcmpi(str1,"INPUT_DATA_NAME")==0)
		{
			strcpy(fs->sourcefile,str3);
			printf("INPUT_DATA_NAME = %s\n",fs->sourcefile);
		}
	}
	fclose(fp);

	if(strcmpi(type,"unsigned")==0&&fs->size==1)
	{
		fs->type=1;
		printf("TYPE = 'unsigned char'\n");
	}
	if(strcmpi(type,"char")==0&&fs->size==1)
	{
		fs->type=2;
		printf("TYPE = 'char'\n");
	}
	if(strcmpi(type,"unsigned")==0&&fs->size==2)
	{
		fs->type=3;
		printf("TYPE = 'unsigned short'\n");
	}
	if(strcmpi(type,"short")==0&&fs->size==2)
	{
		fs->type=4;
		printf("TYPE = 'short'\n");
	}
	if(strcmpi(type,"unsigned long")==0&&fs->size==4)
	{
		fs->type=5;
		printf("TYPE = 'unsigned long'\n");
	}
	if(strcmpi(type,"long")==0&&fs->size==1)
	{
		fs->type=6;
		printf("TYPE = 'long'\n");
	}
	if(strcmpi(type,"float")==0&&fs->size==4)
	{
		fs->type=7;
		printf("TYPE = 'float'\n");
	}
	if(strcmpi(type,"double")==0&&fs->size==8)
	{
		fs->type=8;
		printf("TYPE = 'double'\n");
	}
	if(fs->type==0)
	{
		printf("Invalid DATA_TYPE '%s' + %d bytes\n\n",type,fs->size);
		exit(0);
	}
	if(fs->minlat==360.0)
	{
		fs->minlat=fs->maxlat-(double)(fs->numlat-1)*fs->latspace;
		printf("LOWER_RIGHT_LAT = %lf\n",fs->minlat);
	}
	if(fs->maxlon==-360.0)
	{
		fs->maxlon=fs->minlon+(double)(fs->numlon-1)*fs->lonspace;
		printf("LOWER_RIGHT_LON = %lf\n",fs->maxlon);
	}
	fs->dellat=fs->maxlat-fs->minlat;
	fs->dellon=fs->maxlon-fs->minlon;
	printf("LAT_RANGE = %lf\n",fs->dellat);
	printf("LON_RANGE = %lf\n",fs->dellon);
	if(fs->numlat==0&&fs->dellat>0.0)
	{
		fs->numlat=(fs->dellat+fs->latspace/1000.0)/fs->latspace;
		printf("NUMBER_OF_LAT = %d\n",fs->numlat);
	}
	if(fs->numlon==0&&fs->dellon>0.0)
	{
		fs->numlon=(fs->dellon+fs->lonspace/1000.0)/fs->lonspace;
		printf("NUMBER_OF_LON = %d\n",fs->numlon);
	}
	if(fs->numlon==0)
	{
		printf("Could not find or calculate NUMBER_OF_LON\n");
		exit(0);
	}
	if(fs->numlat==0)
	{
		printf("Could not find or calculate NUMBER_OF_LAT\n");
		exit(0);
	}
	if(fs->minlon>=360.0)
	{
		printf("Could not find or calculate UPPER_LEFT_LON\n");
		exit(0);
	}
	if(fs->maxlat<=-360.0)
	{
		printf("Could not find or calculate UPPER_LEFT_LAT\n");
		exit(0);
	}
	if(fs->numlon==0)
	{
		printf("Could not find or calculate DELTA_LON\n");
		exit(0);
	}
	if(fs->numlat==0)
	{
		printf("Could not find or calculate DELTA_LAT\n");
		exit(0);
	}
}

/**********************************************************************
**
**
**
**********************************************************************/

FILE *fopen_path(char *filename,char *type)

{
	int i,j,k;
	FILE *fp;
	char string[100];

	fp=fopen(filename,type);
	if(fp)
		return(fp);
	for(i=0;i<NumPath;i++)
	{
		sprintf(string,"%s%s",Path[i],filename);
		fp=fopen(string,type);
		if(fp)
			return(fp);
	}
	return(NULL);
}