/*
	splnbin.c

	spb1.bat = cl /AL splnbin.c spline vimage3 graphlib

splnbin [input] [output] [lat] [lon] [dlat] [dlon] [nlat] [nlon] [flex] <optns>

where options:
                             Use          Default

	Extrapolation -------- /e OFF USE     ON NO_USE
   Isopleth Violations -- /v ON          OFF
	Spline Border -------- /b

*/	
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "nomouse.h"
#include "vimage.h"
#include "spline.h"
#include "font.h"

#define DEG_RAD 0.01745329251				/* converts degrees to radians */
/* VGA color stuff */

#define BLACK      0
#define D_GREY     1
#define M_GREY     2
#define L_GREY     3
#define WHITE      4
#define RED        5
#define ORANGE     6
#define YELLOW     7
#define GREEN      8
#define CYAN       9
#define BLUE      10
#define PURPLE    11
#define BUFF2     12
#define BROWN     13
#define D_BLUE    14
#define BUFF      15

#define MAX_PTS 5000
#define MAX_COL 15000

int  UseExt=1,Border,Square=4;
float LL[MAX_PTS],Z[MAX_PTS],Y2[MAX_PTS],U[MAX_PTS],Temp[MAX_PTS];
unsigned Order[MAX_PTS];
float huge Array[MAX_COL][2];

unsigned char Buffer[MAX_PTS];
struct Color Vcolor[16]=
{
	{   0,  0,  0},
	{  46, 46, 46},
	{ 128,128,128},
	{ 191,191,191},
	{ 255,255,255},
	{ 255,  0,  0},
	{ 255,128,  0},
	{ 191,191,  0},
	{   0,255,  0},
	{   0,191,191},
	{  40, 40,255},
	{ 128,  0,255},
	{  95, 81, 67},
	{  96, 28, 14},
	{   0, 20, 60},
	{ 105, 90, 75}
};
long DelAtr;
int ScreenXs=640,ScreenYs=480;
struct grid_stat
{
	float maxlat,minlat,maxlon,minlon;
	float maxz,minz,midz;
	float midlat,midlon,cos;
	float dellat,dellon;
	float latspace,lonspace;
	float flex,scale,zscale;
	int numlat,numlon;
}Fs;
int Zht=128;
FILE *SB_Fp,*SB_Fphdr,*SB_Fpatr;
char InName[100],OutName[100];
extern float MsngVal;
extern int NumMsng;

void splint(float xa[],float ya[],float y2a[],int n,float x,float *y);
void spline(float x[],float y[],int n,float yp1,float ypn,
		float y2[],float u[]);
int dsplint(float xa[],float ya[],float y2a[],int n,float x,float *y,
	float *dist,float dely);
long check_atr(struct grid_stat *fs);
int earth_screen(struct grid_stat fs,float lat, float lon,float z,
	float *x,float *y);
int plot_grid(struct grid_stat fs);
int reorder(float *x,unsigned *order,int num);
int compare(int *a,int *b);
int get_lon(float clat,float *glon,float *z,int num);
int get_lat(float clon,float *glat,float *z,int num);
int open_files(char *name);
int dedupe(float *LL,float *Z,int n);
int video_off(void);

main(int argc,char *argv[])

{
	int i,j,k,n;
	char inname[100],outname[100];
	short numlat,numlon;
	float xscale,yscale,lat,lon,z;
	float x1,y1,x2,y2;
	
	if(argc<10)
	{
		printf("\n\n\n\n\n");
		printf(
"splnbin [in] [out] [lat] [lon] [dlat] [dlon] [nlat] [nlon] [flex] <optns>\n\n");
		printf("where:\n");
		printf("    in -- input file name (without extention)\n");
		printf("   out -- output file name (without extention)\n");
		printf("   lat -- upper left corner latitude\n");
		printf("   lon -- upper left corner longitude\n");
		printf("  dlat -- latitude spacing of lines in degrees\n");
		printf("  dlon -- longitude spacing of lines in degrees\n");
		printf("  nlat -- number of latitude lines\n");
		printf("  nlon -- number of longitude lines\n");
		printf(
"  flex -- spline flexibility (0.0 for no curve to 1.0 for max curve)\n");
		printf("options:\n");
		printf("                          To Use        Default\n\n");
		printf("  Extrapolation -------- /e ON NO_USE   OFF USE\n");
		printf("  Isopleth Violations -- /v OFF         ON\n\n");
		printf("  Spline Border -------- /b OFF         ON\n\n");
		printf(
"or use program 'mapsplin' to set up data and call this program automatically.\n\n");
		exit(0);
	}
	printf("\n\n\n\n\n");
	UseExt=ON;
	Extrapolate=ON;
	IsoV=ON;
	Border=0;
	for(i=1;i<argc;i++)
	{
		if(strcmpi(argv[i],"/b")==0)
		{
			Border=1;
			printf("Spline Border = ON\n");
		}
		if(strcmpi(argv[i],"/e")==0)
		{
			if(i<argc-1&&strcmpi(argv[i+1],"ON")==0)
			{
				Extrapolate=ON;
				printf(
"Extrapolate=ON   Exrapolated values will be used.\n");
			}
			if(i<argc-1&&strcmpi(argv[i+1],"NO_USE")==0)
			{
				UseExt=OFF;
				printf(
"     UseExt=OFF    Extrapolated values will not be used\n");
			}
			if(i+2<argc)
			{
				if(strcmpi(argv[i+2],"ON")==0)
				{
					Extrapolate=ON;
					printf(
"Extrapolate=ON   Exrapolated values will be used.\n");
				}
				if(strcmpi(argv[i+2],"NO_USE")==0)
				{
					UseExt=OFF;
					printf(
"     UseExt=OFF    Extrapolated values will not be used\n");
				}
			}
		}
		if(strcmpi(argv[i],"/v")==0)
		{
			if(strcmpi(argv[i+1],"OFF")==0)
			{
				IsoV=OFF;
				printf(
"       IsoV=OFF    Values may now exceed isopleth limits.\n");
			}
		}
	}
	strcpy(inname,argv[1]);
	strcpy(InName,argv[1]);
	open_files(inname);
	strcpy(outname,argv[2]);
	strcpy(OutName,argv[2]);
	sscanf(argv[3],"%f",&Fs.maxlat);	
	sscanf(argv[4],"%f",&Fs.minlon);	
	sscanf(argv[5],"%f",&Fs.dellat);	
	sscanf(argv[6],"%f",&Fs.dellon);	
	sscanf(argv[7],"%d",&Fs.numlat);	
	sscanf(argv[8],"%d",&Fs.numlon);	
	sscanf(argv[9],"%f",&Fs.flex);	
	Fs.minlat=Fs.maxlat-(float)(Fs.numlat-1)*Fs.dellat;
	Fs.maxlon=Fs.minlon+(float)(Fs.numlon-1)*Fs.dellon;
	printf("\n    Input File Name = '%s'\n",inname);
	printf("   Output File Name = '%s'\n",outname);
	printf("     Upper Left Lat = %8.4f\n",Fs.maxlat);
	printf("     Upper Left Lon = %9.4f\n",Fs.minlon);
	printf("    Lower Right Lat = %8.4f\n",Fs.minlat);
	printf("    Lower Right Lon = %9.4f\n",Fs.maxlon);
	printf("       Lat Spaacing = %0.8g\n",Fs.dellat);
	printf("        Lon Spacing = %0.8g\n",Fs.dellon);
	printf("Number of Lat Lines = %d\n",Fs.numlat);
	printf("Number of Lon Lines = %d\n",Fs.numlon);
	printf("        Flexibility = %0.6g\n",Fs.flex);
	UScale=Fs.flex;
	Fs.midlon=(Fs.maxlon+Fs.minlon)/2.0;
	Fs.midlat=(Fs.maxlat+Fs.minlat)/2.0;
	Fs.cos=cos((double)(Fs.midlat*DEG_RAD));
	DelAtr=check_atr(&Fs);
	Fs.midz=(Fs.maxz+Fs.minz)/2.0;
	Fs.zscale=Zht/(Fs.maxz-Fs.minz);
	Fs.scale=1.0;
	printf("     Bin File Min Z = %0.6g   Max Z = %0.6g\n",Fs.minz,Fs.maxz);
	lat=Fs.maxlat;
	lon=Fs.minlon;
	z=Fs.midz;
	earth_screen(Fs,lat,lon,z,&x1,&y1);
	lat=Fs.minlat;
	lon=Fs.maxlon;
	earth_screen(Fs,lat,lon,z,&x2,&y2);
	yscale=(480-Zht)/(y2-y1);
	lat=Fs.minlat;
	lon=Fs.minlon;
	earth_screen(Fs,lat,lon,z,&x1,&y1);
	lat=Fs.maxlat;
	lon=Fs.maxlon;
	earth_screen(Fs,lat,lon,z,&x2,&y2);
	xscale=640/(x2-x1);
	Fs.scale=xscale;
	if(yscale<xscale)
		Fs.scale=yscale;
	printf("\n                   >>> HIT ANY KEY TO CONTINUE <<<");
	getch();
	plot_grid(Fs);
	fcloseall();
	printf("\n\ntrigrid %s\n\n",OutName);
	execlp("trigrid.exe","trigrid",OutName,NULL);
	printf("\n\nCOULD NOT FIND PROGRAM 'TRIGRID'\n\n\n");
	exit(0);
}

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

int plot_grid(struct grid_stat fs)

{
	int i,j,k,n,m;
	long lk,lj;
	float lat0=fs.maxlat,lon0=fs.minlon,lat1,lon1;
	float lat,lon,x,y,z,xo,yo,zo,dist,coslat;
	int ix0,ix1,iy0,iy1,ix2,ix3,iy2,iy3;
	int jx0,jx1,jy0,jy1,jx2,jx3,jy2,jy3;
	int status,val;
	FILE *fptemp,*fp1,*fp2,*fp3,*fplbl,*fpdat,*fplbs,*fpbdr;
	long row=fs.numlat,col=fs.numlon,subcol=(long)MAX_COL/row;
	long numsub=(col+subcol-1)/subcol,load=0,a1=0,a2=col*2,a3=col*4;
	int prntcol=subcol*2;
	char string[100];
	float z1,z2,d1,d2,*array=Array[a3],zm=fs.minz,zr=fs.maxz-fs.minz;
	float zmin=MsngVal,zmax=MsngVal;
	long red1=0,red2=0,red12=0,yellow1=0,yellow2=0,offset;
	int msng;

	lat1=lat0-fs.dellat*(float)(fs.numlat-1);
	lon1=lon0+fs.dellon*(float)(fs.numlon-1);
	if(prntcol>(int)(col*2))
		prntcol=col*2;
	video_on();
	set_lut();
	paint_box(9,BUFF,Buffer,0,0,ScreenXs,ScreenYs);
	z=fs.maxz;
	earth_screen(fs,lat0,lon0,z,&x,&y);
	ix0=x;iy0=y;
	earth_screen(fs,lat0,lon1,z,&x,&y);
	ix1=x;iy1=y;
	earth_screen(fs,lat1,lon1,z,&x,&y);
	ix2=x;iy2=y;
	earth_screen(fs,lat1,lon0,z,&x,&y);
	ix3=x;iy3=y;
	z=fs.minz;
	earth_screen(fs,lat0,lon0,z,&x,&y);
	jx0=x;jy0=y;
	earth_screen(fs,lat0,lon1,z,&x,&y);
	jx1=x;jy1=y;
	earth_screen(fs,lat1,lon1,z,&x,&y);
	jx2=x;jy2=y;
	earth_screen(fs,lat1,lon0,z,&x,&y);
	jx3=x;jy3=y;
	plotln(0,ix0,iy0,ix1,iy1,BLACK);
	plotln(0,ix1,iy1,ix2,iy2,BLACK);
	plotln(0,ix2,iy2,ix3,iy3,BLACK);
	plotln(0,ix3,iy3,ix0,iy0,BLACK);
	plotln(0,jx0,jy0,jx1,jy1,BLACK);
	plotln(0,jx1,jy1,jx2,jy2,BLACK);
	plotln(0,jx2,jy2,jx3,jy3,BLACK);
	plotln(0,jx3,jy3,jx0,jy0,BLACK);
	plotln(0,ix0,iy0,jx0,jy0,BLACK);
	plotln(0,ix1,iy1,jx1,jy1,BLACK);
	plotln(0,ix2,iy2,jx2,jy2,BLACK);
	plotln(0,ix3,iy3,jx3,jy3,BLACK);
	box(0,2,4,YELLOW,5,5);
	plot_font_h(L_GREY,12,2,"Point Violated Isopleths",2,Font);	
	if(Border==0)
	{
		box(0,2,19,RED,5,5);
		plot_font_h(L_GREY,12,17,"Extrapolated Point",2,Font);	
	}
	fptemp=fopen("temp1.d32","wb");
	if(!fptemp)
		term_error("Could not open 'temp1.d32' to write.");
/* PROCESSING VERTICAL LINES */
	if(Border>0)	/* make a spline border */
	{
		paint_box(0,PURPLE,Buffer,0,31,168,13);
		plot_font_h(WHITE,4,32," MAKING SPLINE BORDER  ",2,Font);	
		fpbdr=fopen("border.tmp","w+b");
		if(!fpbdr)
		{
			printf(
				"Spline Border selected and 'temp file could not be opened.\n\n");
			exit(0);
		}
		n=get_lon(lat0,LL,Z,MAX_PTS);
		if(n<3)
			term_error("Not enough points at lat0");
		for(m=0;m<n;m++)
			Order[m]=m;
		qsort((char *)Order,n,sizeof(int),compare);
		reorder(Z,Order,n);
		reorder(LL,Order,n);
		n=dedupe(LL,Z,n);
		spline(LL,Z,n,1.0e30,1.0e30,Y2,U);
		lon=lon0;
		for(j=0;j<fs.numlon;j++)
		{
			status=dsplint(LL,Z,Y2,n,lon,&z,&dist,(float)DelAtr);	
			lon+=fs.dellon;
			offset=(long)j*8;
			fseek(fpbdr,offset,SEEK_SET);
			fwrite((char *)&z,sizeof(float),1,fpbdr);
		}
		n=get_lon(lat1,LL,Z,MAX_PTS);
		if(n<3)
			term_error("Not enough points at lat1");
		for(m=0;m<n;m++)
			Order[m]=m;
		qsort((char *)Order,n,sizeof(int),compare);
		reorder(Z,Order,n);
		reorder(LL,Order,n);
		n=dedupe(LL,Z,n);
		spline(LL,Z,n,1.0e30,1.0e30,Y2,U);
		lon=lon0;
		for(j=0;j<fs.numlon;j++)
		{
			status=dsplint(LL,Z,Y2,n,lon,&z,&dist,(float)DelAtr);	
			lon+=fs.dellon;
			offset=(long)j*8+4;
			fseek(fpbdr,offset,SEEK_SET);
			fwrite((char *)&z,sizeof(float),1,fpbdr);
		}
		rewind(fpbdr);
	}
	lon=lon0;
	paint_box(0,BLUE,Buffer,0,31,168,13);
	plot_font_h(WHITE,4,32,"  PROCESSING LONGITUDE ",2,Font);	
	load=0;
	lj=lk=0;
	for(j=0;j<fs.numlon;j++)
	{
		n=get_lat(lon,LL,Z,MAX_PTS);
		if(Border>0)
		{
			if(fread((char *)&Z[n],sizeof(float),2,fpbdr)<2)
				term_error("Could not read from 'border.tmp'.");
			LL[n]=lat0;
			LL[n+1]=lat1;
			n+=2;
		}
		if(n>2)
		{
			for(m=0;m<n;m++)
				Order[m]=m;
			qsort((char *)Order,n,sizeof(int),compare);
			reorder(Z,Order,n);
			reorder(LL,Order,n);
			n=dedupe(LL,Z,n);
			spline(LL,Z,n,1.0e30,1.0e30,Y2,U);
			lat=lat0;
		}
		for(i=0;i<fs.numlat;i++)
		{
			status=dsplint(LL,Z,Y2,n,lat,&z,&dist,(float) DelAtr);
			earth_screen(fs,lat,lon,z,&x,&y);
			if(i>0&&z!=MsngVal&&zo!=MsngVal)
				plotln(0,(int)x,(int)y,(int)xo,(int)yo,L_GREY);
			if(status<1&&z!=MsngVal)
			{
				box(0,(int)x-2,(int)y-2,YELLOW,5,5);
				yellow1+=1;
			}
			Array[subcol*(long)i+lj][0]=z;
			Array[subcol*(long)i+lj][1]=dist;
			lat-=fs.dellat;
			xo=x;
			yo=y;
			zo=z;
		}
		lj+=1;
		if((j+1)%subcol==0)
		{
			paint_box(0,CYAN,Buffer,0,31,168,13);
			plot_font_h(BLACK,4,32,"    SAVING LONGITUDE ",2,Font);	
			fseek(fptemp,lk*(long)8,SEEK_SET);
			offset=(long)(col-subcol)*(long)8;
			for(i=0;i<row;i++)
			{
				if(fwrite((char *)Array[(long)i*subcol],sizeof(float),
						prntcol,fptemp)<prntcol)
					term_error(
					"Could not write to 'temp1.d32'.  Disk may be full.\n\n");
				if(offset>0)
					fseek(fptemp,offset,SEEK_CUR);
			}
			lj=0;
			load+=1;
			lk=load*subcol;
			if((lk+subcol)>col)
				prntcol=(col-lk)*2;
			paint_box(0,BLUE,Buffer,0,31,168,13);
			plot_font_h(WHITE,4,32,"  PROCESSING LONGITUDE ",2,Font);	
		}
		lon+=fs.dellon;
	}
	if(load==0||prntcol<(int)subcol*2)
	{
		paint_box(0,CYAN,Buffer,0,31,168,13);
		plot_font_h(BLACK,4,32,"    SAVING LONGITUDE ",2,Font);	
		offset=(long)(col-prntcol/2)*(long)8;
		fseek(fptemp,lk*(long)8,SEEK_SET);
		for(i=0;i<row;i++)
		{
			if(fwrite((char *)Array[(long)i*subcol],sizeof(float),
					prntcol,fptemp)<prntcol)
				term_error(
					"Could not write to 'temp1.d32'.  Disk may be full.\n\n");
			if(offset>0)
				fseek(fptemp,offset,SEEK_CUR);
		}
	}
	fclose(fptemp);
	fptemp=fopen("temp2.d32","wb");
	if(!fptemp)
		term_error("Could not open 'temp2.d32' to write.");
/* PROCESSING HORIZONTAL LINES */
	if(Border>0)	/* make a spline border */
	{
		paint_box(0,PURPLE,Buffer,0,31,168,13);
		plot_font_h(WHITE,4,32," MAKING SPLINE BORDER  ",2,Font);	
		rewind(fpbdr);
		n=get_lat(lon0,LL,Z,MAX_PTS);
		for(m=0;m<n;m++)
			Order[m]=m;
		if(n>2)
		{
			qsort((char *)Order,n,sizeof(int),compare);
			reorder(Z,Order,n);
			reorder(LL,Order,n);
			n=dedupe(LL,Z,n);
			spline(LL,Z,n,1.0e30,1.0e30,Y2,U);
		}
		lat=lat0;
		for(j=0;j<fs.numlat;j++)
		{
			status=dsplint(LL,Z,Y2,n,lat,&z,&dist,(float)DelAtr);	
			lat-=fs.dellat;
			offset=(long)j*8;
			fseek(fpbdr,offset,SEEK_SET);
			if(fwrite((char *)&z,sizeof(float),1,fpbdr)<1)
				term_error("Could not write Z for lon0.\n");
		}
		n=get_lat(lon1,LL,Z,MAX_PTS);
		if(n>2)
		{
			for(m=0;m<n;m++)
				Order[m]=m;
			qsort((char *)Order,n,sizeof(int),compare);
			reorder(Z,Order,n);
			reorder(LL,Order,n);
			n=dedupe(LL,Z,n);
			spline(LL,Z,n,1.0e30,1.0e30,Y2,U);
		}
		lat=lat0;
		for(j=0;j<fs.numlat;j++)
		{
			status=dsplint(LL,Z,Y2,n,lat,&z,&dist,(float)DelAtr);	
			lat-=fs.dellat;
			offset=(long)j*8+4;
			fseek(fpbdr,offset,SEEK_SET);
			if(fwrite((char *)&z,sizeof(float),1,fpbdr)<1)
				term_error("Could not write Z for lon1.\n");
		}
		rewind(fpbdr);
	}
	paint_box(0,GREEN,Buffer,0,31,168,13);
	plot_font_h(BLACK,4,32,"  PROCESSING LATITUDE  ",2,Font);	
	lat=lat0;
	for(j=0;j<fs.numlat;j++)
	{
		coslat=cos((double)(lat*DEG_RAD));
		n=get_lon(lat,LL,Z,MAX_PTS);
		if(Border>0)
		{
			if(fread((char *)&Z[n],sizeof(float),2,fpbdr)<2)
				term_error("Could not read from 'border.tmp'.");
			LL[n]=lon0;
			LL[n+1]=lon1;
			n+=2;
		}
		if(n>2)
		{
			for(m=0;m<n;m++)
				Order[m]=m;
			qsort((char *)Order,n,sizeof(int),compare);
			reorder(Z,Order,n);
			reorder(LL,Order,n);
			n=dedupe(LL,Z,n);
			spline(LL,Z,n,1.0e30,1.0e30,Y2,U);
		}
		lon=lon0;
		for(i=0;i<fs.numlon;i++)
		{
			status=dsplint(LL,Z,Y2,n,lon,&z,&dist,(float) DelAtr);
			earth_screen(fs,lat,lon,z,&x,&y);
			if(i>0&&z!=MsngVal&&zo!=MsngVal)
				plotln(0,(int)x,(int)y,(int)xo,(int)yo,L_GREY);
			if(status<1&&z!=MsngVal)
			{
				box(0,(int)x-2,(int)y-2,YELLOW,5,5);
				yellow2+=1;
			}
			Array[i][0]=z;
			Array[i][1]=dist*coslat;
			lon+=fs.dellon;
			xo=x;
			yo=y;
			zo=z;
		}
		if(fwrite((char *)Array,sizeof(float),(int)(col*2),fptemp)<(int)(col*2))
			term_error(
				"Could not write to 'temp2.d32'.  Disk may be full.\n\n");
		lat-=fs.dellat;
	}
	fclose(fptemp);
	paint_box(0,YELLOW,Buffer,0,31,168,13);
	plot_font_h(BLACK,8,32,"    COMBINING FILES    ",2,Font);	
	fp1=fopen("temp1.d32","rb");
	if(!fp1)
		term_error("Could not open 'temp1.d32' to read.");
	fp2=fopen("temp2.d32","rb");
	if(!fp2)
		term_error("Could not open 'temp2.d32' to read.");
	sprintf(string,"%s.d32",OutName);
	fp3=fopen(string,"wb");
	if(!fp3)
	{
		sprintf(string,"Could not open '%s.d32' to write.",OutName);
		term_error(string);
	}
	fpdat=fopen("temp.dat","wb");
	if(!fpdat)
		term_error("Could not open 'temp.dat' to write.");
	for(i=0;i<row;i++)
	{
		if(fread((char *)Array[a1],sizeof(float),(int)(col*2),fp1)<(int)(col*2))
			term_error("Hit EOF on 'temp1.d32'\n");
		if(fread((char *)Array[a2],sizeof(float),(int)(col*2),fp2)<(int)(col*2))
			term_error("Hit EOF on 'temp2.d32'\n");
		for(j=0;j<col;j++)
		{
			z1=Array[a1+j][0];
			d1=Array[a1+j][1];
			z2=Array[a2+j][0];
			d2=Array[a2+j][1];
			if(Border>0)
			{
				if(d1<0.0)
					d1*=-1.0;
				if(d2<0.0)
					d2*=-1.0;
			}
			if(d1<0.0&&d2<0.0&&Border==0)
				red12+=1;
			if(d1<0.0&&d2<0.0)
				array[j]=MsngVal;
			else if(d1<0.0&&d2>=0.0)
				array[j]=z2;
			else if(d2<0.0&&d2>=0.0)
				array[j]=z1;
			else if(d1==0.0||(UseExt==OFF&&d2<0.0))
				array[j]=z1;
			else if(d2==0.0||(UseExt==OFF&&d1<0.0))
				array[j]=z2;
			else
			{
				if(Square>0)
					array[j]=(z2*d1*d1*d1*d1+z1*d2*d2*d2*d2)/
								(d1*d1*d1*d1+d2*d2*d2*d2);
				else
					array[j]=(z1*d2+z2*d1)/(d1+d2);
			}
			Buffer[j]=(array[j]-zm)/zr*255.0;
			if(zmax==MsngVal&&zmin==MsngVal&&array[j]!=MsngVal)
			{
				zmax=zmin=array[j];
			}
			else if(array[j]!=MsngVal)
			{
				if(array[j]>zmax)
					zmax=array[j];
				if(array[j]<zmin)
					zmin=array[j];
			}
		}
		if(fwrite((char *)array,sizeof(float),(int)col,fp3)<col)
		{
			sprintf(string,"Could not write to '%s.d32'",OutName);
			term_error(string);
		}
		if(fwrite((char *)Buffer,sizeof(char),(int)col,fpdat)<col)
		{
			term_error("Could not write to 'temp.dat'");
		}
	}
	paint_box(0,RED,Buffer,0,29,175,15);
	plot_font_h(WHITE,8,32,"< HIT ANY KEY TO EXIT >",2,Font);	
getch();
	video_off();
	if(red1+red2+red12+yellow1+yellow2>0)
		printf("\n\nWARNING!\n\n");
	if(yellow1>0)
		printf(
"longitude (vertical) lines contained %ld isopleth violations (YELLOW)\n",
			yellow1);
	if(yellow2>0)
		printf(
"latitude (horizontal) lines contained %ld isopleth violations (YELLOW)\n",
			yellow2);
	if(red1>0)
		printf(
"longitude (vertical) lines contained %ld extrapolations (RED)\n",red1);
	if(red2>0)
		printf(
"latitude (horizontal) lines contained %ld extrapolations (RED)\n",red2);
	if(red12>0)
		printf("   and %ld extrapolations coincided\n",red12);
	if(red1+red2+red12+yellow1+yellow2>0)
		printf("\n\n");
	fclose(fp1);
	unlink("temp1.d32");
	fclose(fp2);
	unlink("temp2.d32");
	fclose(fp3);
	fclose(fpdat);
	fplbl=fopen("temp.lbl","wt");
	if(!fplbl)
		term_error("Could not open 'temp.lbl' to write.");
	fprintf(fplbl,"FILE_TYPE          = IMAGE\n");
	fprintf(fplbl,"IMAGE_LINES        = %6ld\n",row);
	fprintf(fplbl,"LINE_SAMPLES       = %6ld\n",col);
	fprintf(fplbl,"IMAGE_POINTER      = 'TEMP.DAT'\n");
	fclose(fplbl);
	sprintf(string,"%s.lbs",OutName);
	fplbs=fopen(string,"wt");
	if(!fplbs)
	{
		sprintf(string,"Could not open '%s.lbs' to write.",OutName);
		term_error(string);
	}
	fprintf(fplbs,"FILE_TYPE          = ELEVATION_ARRAY\n");
	fprintf(fplbs,"IMAGE_LINES        = %6ld\n",row);
	fprintf(fplbs,"LINE_SAMPLES       = %6ld\n",col);
	fprintf(fplbs,"UPPER_LEFT_LAT     = %10.5f\n",fs.maxlat);
	fprintf(fplbs,"UPPER_LEFT_LON     = %10.5f\n",fs.minlon);
	fprintf(fplbs,"LOWER_RIGHT_LAT    = %10.5f (optional)\n",fs.minlat);
	fprintf(fplbs,"LOWER_RIGHT_LON    = %10.5f (optional)\n",fs.maxlon);
	fprintf(fplbs,"DELTA_LAT          = %0.8f\n",fs.dellat);
	fprintf(fplbs,"DELTA_LON          = %0.8f\n",fs.dellon);
	fprintf(fplbs,"NUMBER_OF_LAT      = %6d\n",fs.numlat);
	fprintf(fplbs,"NUMBER_OF_LON      = %6d\n",fs.numlon);
	fprintf(fplbs,"MAXIMUM_Z          = %0.5f (optional)\n",zmax);
	fprintf(fplbs,"MINIMUM_Z          = %0.5f (optional)\n",zmin);
	fprintf(fplbs,"MISSING_DATA_VALUE = %0.1f (optional)\n",MsngVal);
	fprintf(fplbs,"SPLINE_FLEXIBILITY = %0.5f (information only)\n",fs.flex);
	fprintf(fplbs,"   INPUT_DATA_NAME = '%s' (information only)\n",InName);
	fprintf(fplbs,"DATA_TYPE          = float\n");
	fprintf(fplbs,"DATA_SIZE          = 4 bytes\n");
	fprintf(fplbs,"DATA_FILE_POINTER  = '%s.D32'\n",OutName);
	fclose(fplbs);
	return(1);
}

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

int term_error(char *message)

{
	video_off();
	printf("\n%s\n\n",message);
	exit(0);
}



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

int earth_screen(struct grid_stat fs,float lat, float lon,float z,
	float *x,float *y)

{
	float ylat=fs.midlat-lat,xlon=(lon-fs.midlon)*fs.cos;

	*x=(xlon-ylat)*fs.scale+(float)(ScreenXs/2);
	if(z!=MsngVal)
		*y=(xlon+ylat)*fs.scale/1.414+(fs.midz-z)*fs.zscale+(float)(ScreenYs/2);
	else
		*y=(xlon+ylat)*fs.scale/1.414+(float)(ScreenYs/2);
}


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

int video_off()
{
	SetVideoMode(0);
}

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

int video_on()

{
	int i,row;

	i=GetVideoBoardID();
	if(i==0)
		VideoType='E';
	else if(i==1&&(VideoType=='X'||VideoType=='Y'))
		VideoType='V';
	if(VideoType=='X')
	{
		row=SetVideoMode(480,&B_Id);
		if(B_Id.color==256&&B_Id.row==480)
		{
			ScreenXs=640;
			ScreenYs=480;
		}
		else if(B_Id.color==256&&B_Id.row==400)
			VideoType='Y';
		else
			VideoType='V';
	}
	if(VideoType=='Y')
	{
		row=SetVideoMode(400,&B_Id);
		if(B_Id.row<400)
			VideoType='V';
		else
		{
			ScreenXs=640;
			ScreenYs=row;
		}
	}
	if(VideoType=='V')
	{
		row=SetVideoMode(0x12,&B_Id);
		if(B_Id.row!=480||row!=480||B_Id.color!=16||B_Id.col!=640)
			VideoType='E';
		else
		{
			ScreenXs=640;
			ScreenYs=480;
		}
	}
	if(VideoType=='S')
	{
		row=SetVideoMode(0x13,&B_Id);
		if(B_Id.row!=200)
		{
			SetVideoMode(0,&B_Id);
			printf("Could not boot color board.\n");
			exit(0);
		}
		ScreenXs=320;
		ScreenYs=200;
	}
	if(VideoType=='E')
	{
		row=SetVideoMode(0x10,&B_Id);
		if(B_Id.row!=350||row<350||B_Id.col!=640)
		{
			SetVideoMode(0,&B_Id);
			printf("Could not boot color board.\n");
			exit(0);
		}
		ScreenXs=640;
		ScreenYs=350;
	}
	if(B_Id.row<0)
	{
		SetVideoMode(0,&B_Id);
		printf("Could not boot board -- row = %d\n\n",row);
		exit(0);
	}
	set_lut();
}


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

long check_atr(struct grid_stat *fs)

{
	int i;
	float mll[4],mnln,mxln,mnlt,mxlt;
	char string[100];
	int numsec;
	long attrib,maxz,minz,atr0,delatr,mindel=1000000000;

	rewind(SB_Fpatr);
	if(fread((char *)&attrib,sizeof(long),1,SB_Fpatr)<1)
	{
		video_off();
		exit(0);
	}
	maxz=minz=atr0=attrib;
	while(fread((char *)&attrib,sizeof(long),1,SB_Fpatr)>0)
	{
		if(attrib==(long)MsngVal)
		{
			video_off();
			printf("attribute (%ld) = MISSING_DATA_VALUE\n",attrib);
			exit(0);
		}
		if(attrib>maxz)
			maxz=attrib;
		if(attrib<minz)
			minz=attrib;
		delatr=atr0-attrib;
		if(delatr<0)
			delatr*=-1;
		if(delatr>0&&delatr<mindel)
			mindel=delatr;
	}
	fs->maxz=maxz;
	fs->minz=minz;
	return(mindel);
}

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

int set_lut()

{
	int i;
	struct Color lut[256];

	for(i=16;i<256;i++)
	{
		lut[i].r=i;
		lut[i].g=i;
		lut[i].b=i;
	}
	for(i=0;i<16;i++)
	{
		lut[i].r=Vcolor[i].r;
		lut[i].g=Vcolor[i].g;
		lut[i].b=Vcolor[i].b;
	}
	WritePalette(lut); 
}

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

int dedupe(float *x,float *y,int n)

{
	int i,j,k;
	int dup=0;

	for(i=1;i<n;i++)
	{
		if(x[i]==x[i-1])
		{
			dup+=1;
		}
		if(dup>0)
		{
			x[i-dup]=x[i];
			y[i-dup]=y[i];
		}
	}
	return(n-dup);
}




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

int reorder(float *x,unsigned *order,int num)

{
	int i;

	for(i=0;i<num;i++)
		Temp[i]=x[order[i]];
	for(i=0;i<num;i++)
		x[i]=Temp[i];
}

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

int compare(int *a,int *b)

{
	if(LL[*a]>LL[*b])
		return(1);
	if(LL[*a]<LL[*b])
		return(-1);
	else
		return(0);

}


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

int get_lon(float clat,float *glon,float *z,int num)

{
	int i,j,k,m,n=0,p,oo=0;
	double x,y,xo,yo;
	int ix,iy,ixo=-1,iyo=-1;
	double dx,dy,dist;
	int secnum,numpt,numsec,ival;
	float latlon[2],mll[4],mnln,mxln,mnlt,mxlt;
	float lat,lon,lato,lono,clon;
	char string[40],pth[40],type[40];
	int doit,numplot=0;
	int first_seg=-1,last_seg=-1,atr=0;
	long offset,attrib=-1;

	rewind(SB_Fp);
	rewind(SB_Fphdr);
	rewind(SB_Fpatr);
	if(fread((char *)&numsec,sizeof(int),1,SB_Fphdr)<1)
	{
		video_off();
		exit(0);
	}
	if(fread((char *)&numsec,sizeof(int),1,SB_Fp)<1)
	{
		video_off();
		exit(0);
	}
	for(j=0;j<numsec;j++)
	{
		if(kbhit()>0)
		{
			ival=getch();
			if(ival==27)
			{
				video_off();
				unlink("temp1.d32");
				unlink("temp2.d32");
				exit(0);
			}
		}
		if(fread((char *)&offset,sizeof(long),1,SB_Fphdr)<1)
		{
			video_off();
			exit(0);
		}
		if(fread((char *)mll,sizeof(float),4,SB_Fphdr)<1)
		{
			video_off();
			exit(0);
		}
		if(fread((char *)&attrib,sizeof(long),1,SB_Fpatr)<1)
		{
			video_off();
			exit(0);
		}
		doit=1;
		mnln=mll[0];mxln=mll[1];mnlt=mll[2];mxlt=mll[3];
		if(mnlt>clat||mxlt<clat)
			doit=0;
		if(doit==1)
		{
			numplot+=1;
			fseek(SB_Fp,offset,SEEK_SET);
			if(fread((char *)&numpt,sizeof(int),1,SB_Fp)<1)
			{
				video_off();
				exit(0);
			}
			if(fread((char *)mll,sizeof(float),4,SB_Fp)<1)
			{
				video_off();
				exit(0);
			}
			m=0;
			for(i=0;i<numpt;i++)
			{
				if(fread((char *)latlon,sizeof(float),2,SB_Fp)<1)
				{
					video_off();
					exit(0);
				}
				lat=latlon[0];
				lon=latlon[1];
				if(doit==1)
				{
					if(lat==clat)
					{
						glon[n]=lon;
						z[n]=attrib;
						n+=1;
					}
					else if(lat==lato)
						clon=(lon+lono)/2.0;
					else if(i>0)
					{
	/* check for crossing line */
						if((lat>clat&&lato<clat)||(lat<clat&&lato>clat))
						{
							glon[n]=lono+(lon-lono)*(clat-lato)/(lat-lato);
							z[n]=attrib;
							n+=1;
						}
					}
					if(n>=num)
					{
						video_off();
						printf("Too many hits\n");
						exit(0);
					}
					lato=lat;
					lono=lon;
				}
			}
		}
	}
	return(n);
}

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

int get_lat(float clon,float *glat,float *z,int num)

{
	int i,j,k,m,n=0,p,oo=0;
	double x,y,xo,yo;
	int ix,iy,ixo=-1,iyo=-1;
	double dx,dy,dist;
	int secnum,numpt,numsec,ival;
	float latlon[2],mll[4],mnln,mxln,mnlt,mxlt;
	float lat,lon,lato,lono,clat;
	char string[40],pth[40],type[40];
	int doit,numplot=0;
	int first_seg=-1,last_seg=-1,atr=0;
	long offset,attrib=-1;

	rewind(SB_Fp);
	rewind(SB_Fphdr);
	rewind(SB_Fpatr);
	if(fread((char *)&numsec,sizeof(int),1,SB_Fphdr)<1)
	{
		video_off();
		exit(0);
	}
	if(fread((char *)&numsec,sizeof(int),1,SB_Fp)<1)
	{
		video_off();
		exit(0);
	}
	for(j=0;j<numsec;j++)
	{
		if(kbhit()>0)
		{
			ival=getch();
			if(ival==27)
			{
				video_off();
				unlink("temp1.d32");
				unlink("temp2.d32");
				exit(0);
			}
		}
		if(fread((char *)&offset,sizeof(long),1,SB_Fphdr)<1)
		{
			video_off();
			exit(0);
		}
		if(fread((char *)mll,sizeof(float),4,SB_Fphdr)<1)
		{
			video_off();
			exit(0);
		}
		if(fread((char *)&attrib,sizeof(long),1,SB_Fpatr)<1)
		{
			video_off();
			exit(0);
		}
		doit=1;
		mnln=mll[0];mxln=mll[1];mnlt=mll[2];mxlt=mll[3];
		if(mnln>clon||mxln<clon)
			doit=0;
		if(doit==1)
		{
			numplot+=1;
			fseek(SB_Fp,offset,SEEK_SET);
			if(fread((char *)&numpt,sizeof(int),1,SB_Fp)<1)
			{
				video_off();
				exit(0);
			}
			if(fread((char *)mll,sizeof(float),4,SB_Fp)<1)
			{
				video_off();
				exit(0);
			}
			m=0;
			for(i=0;i<numpt;i++)
			{
				if(fread((char *)latlon,sizeof(float),2,SB_Fp)<1)
				{
					video_off();
					exit(0);
				}
				lat=latlon[0];
				lon=latlon[1];
				if(doit==1)
				{
					if(lon==clon)
					{
						glat[n]=lat;
						z[n]=attrib;
						n+=1;
					}
					else if(lon==lono)
						clat=(lat+lato)/2.0;
					else if(i>0)
					{
	/* check for crossing line */
						if((lon>clon&&lono<clon)||(lon<clon&&lono>clon))
						{
							glat[n]=lato+(lat-lato)*(clon-lono)/(lon-lono);
							z[n]=attrib;
							n+=1;
						}
					}
					if(n>=num)
					{
						video_off();
						printf("Too many hits\n");
						exit(0);
					}
					lato=lat;
					lono=lon;
				}
			}
		}
	}
	return(n);
}

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

int open_files(char *name)

{
	char string[100];

	sprintf(string,"%s.bin",name);
	SB_Fp=fopen(string,"rb");
	if(!SB_Fp)
	{
		printf("could not open '%s'\n",string);
		exit(0);
	}
	sprintf(string,"%s.hdr",name);
	SB_Fphdr=fopen(string,"rb");
	if(!SB_Fphdr)
	{
		printf("could not open '%s'\n",string);
		exit(0);
	}
	sprintf(string,"%s.atr",name);
	SB_Fpatr=fopen(string,"rb");
	if(!SB_Fpatr)
	{
		printf("could not open '%s'\n",string);
		exit(0);
	}
}